diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index b8a70ee7a3c..4c76060ba6b 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -95,7 +95,6 @@ SET_BOOL_PROP(EclFlowProblem, EnableBrine, false); SET_TYPE_PROP(EclFlowProblem, EclWellModel, Opm::BlackoilWellModel); SET_TAG_PROP(EclFlowProblem, LinearSolverSplice, FlowIstlSolver); - END_PROPERTIES namespace Opm { @@ -145,6 +144,10 @@ namespace Opm { typedef Dune::BlockVector BVector; typedef ISTLSolverEbos ISTLSolverType; + + typedef Dune::BCRSMatrix< Dune::FieldMatrix> PressureMatrixType; + typedef Dune::BlockVector> PressureVectorType; + typedef Dune::FlexibleSolver PressureSolverType; //typedef typename SolutionVector :: value_type PrimaryVariables ; // --------- Public methods --------- @@ -488,27 +491,29 @@ namespace Opm { > makePressureSolver(Dune::BCRSMatrix< Dune::FieldMatrix>& pmatrix) { boost::property_tree::ptree prm; - if (std::filesystem::exists("pressuresolver.json") ) { + if (std::filesystem::exists(param_.pressure_solver_json_) ) { try{ - boost::property_tree::read_json("pressuresolver.json", prm); + boost::property_tree::read_json(param_.pressure_solver_json_, prm); }catch(...){ - OPM_THROW(std::logic_error,"failed in passing configuration file pressuresolver.json"); + OPM_THROW(std::logic_error,"failed in passing configuration file:" + param_.pressure_solver_json_); } } else { // using a default setup for pressure solver std::stringstream ss; + prm.put("solver","umfpack"); + prm.put("blocksize",1); + prm.put("verbosity",1); + /* ss << "{\n" " \"solver\": \"umfpack\",\n" " \"blocksize\": \"1\"\n" "}"; boost::property_tree::read_json(ss, prm); + */ } - using PressureMatrixType = Dune::BCRSMatrix< Dune::FieldMatrix>; - using PressureVectorType = Dune::BlockVector>; std::any parallelInformation; extractParallelGridInformationToISTL(ebosSimulator_.vanguard().grid(), parallelInformation); using AbstractOperatorType = Dune::AssembledLinearOperator; - using PressureSolverType = Dune::FlexibleSolver; std::unique_ptr operator_for_flexiblesolver; std::function weightsCalculator;// dummy std::unique_ptr pressureSolver; @@ -542,7 +547,18 @@ namespace Opm { return pressureSolver; } - + bool shouldCreatePressureSolver(){ + if(!pressureSolver_){ + return true; + } + if(param_.reuse_pressure_solver_ > 0){ + return false; + }else{ + return true; + } + assert(false); + + } /// Solve the Jacobian system Jx = r where J is the Jacobian and /// r is the residual. void solveJacobianSystem(BVector& x) @@ -551,28 +567,43 @@ namespace Opm { auto& ebosJac = ebosSimulator_.model().linearizer().jacobian(); auto& ebosResid = ebosSimulator_.model().linearizer().residual(); // NB when tested the linear solver an the martrix could be part of linearizer - using PressureMatrixType = Dune::BCRSMatrix< Dune::FieldMatrix>; - using PressureVectorType = Dune::BlockVector>; + int pressureVarIndex=1; // true impes case could add case with trivial weighs i equations is modifind with weights // this would make a newton based method if derivative of the weights are takein into account BVector weights = this->getPressureWeights(pressureVarIndex); - PressureMatrixType pmatrix = - PressureHelper::makePressureMatrix( - ebosJac.istlMatrix(), - pressureVarIndex, - weights); + if(!pmatrix_){ + // make matrix structure and fill with elements + pmatrix_ = + std::make_unique(PressureHelper::makePressureMatrix( + ebosJac.istlMatrix(), + pressureVarIndex, + weights)); + }else{ + //only update elements of matrix + // assume matrix structure never change + PressureHelper::makePressureMatrixEntries( + *pmatrix_, + ebosJac.istlMatrix(), + pressureVarIndex, + weights); + } PressureVectorType rhs(ebosResid.size(),0); PressureHelper::moveToPressureEqn(ebosResid, rhs, weights); - - using PressureSolverType = Dune::FlexibleSolver; - std::unique_ptr pressureSolver = makePressureSolver(pmatrix); + if(this->shouldCreatePressureSolver()){ + pressureSolver_ = makePressureSolver(*pmatrix_); + }else{ + pressureSolver_->preconditioner().update(); + } + PressureVectorType xp(x.size(),0); Dune::InverseOperatorResult res; - pressureSolver->apply(xp,rhs, res); + pressureSolver_->apply(xp,rhs, res); /* bool write_pressure_system = false; if(write_pressure_system){ @@ -1114,6 +1145,9 @@ namespace Opm { double drMaxRel() const { return param_.dr_max_rel_; } double maxResidualAllowed() const { return param_.max_residual_allowed_; } double linear_solve_setup_time_; + std::unique_ptr pressureSolver_; + std::unique_ptr pmatrix_; + public: std::vector wasSwitched_; }; diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index e182935aefb..62daae2a514 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -47,13 +47,15 @@ NEW_PROP_TAG(UpdateEquationsScaling); NEW_PROP_TAG(UseUpdateStabilization); NEW_PROP_TAG(MatrixAddWellContributions); NEW_PROP_TAG(EnableWellOperabilityCheck); -NEW_PROP_TAG(SimulationType); // for sequential NEW_PROP_TAG(MaxStrictIterSeq); NEW_PROP_TAG(ToleranceCnvSeq); NEW_PROP_TAG(ToleranceCnvRelaxedSeq); NEW_PROP_TAG(TolerancePressure); +NEW_PROP_TAG(ReusePressureSolver); +NEW_PROP_TAG(SimulationType); +NEW_PROP_TAG(PressureSolverJson); // parameters for multisegment wells NEW_PROP_TAG(TolerancePressureMsWells); @@ -92,7 +94,6 @@ SET_INT_PROP(FlowModelParameters, MaxInnerIterWells, 50); SET_INT_PROP(FlowModelParameters, StrictInnerIterMsWells, 40); SET_SCALAR_PROP(FlowModelParameters, RegularizationFactorMsw, 1); SET_BOOL_PROP(FlowModelParameters, EnableWellOperabilityCheck, true); -SET_STRING_PROP(FlowModelParameters, SimulationType, "implicit"); SET_SCALAR_PROP(FlowModelParameters, RelaxedFlowTolInnerIterMsw, 1); SET_SCALAR_PROP(FlowModelParameters, RelaxedPressureTolInnerIterMsw, 0.5e5); @@ -102,6 +103,9 @@ SET_INT_PROP(FlowModelParameters, MaxStrictIterSeq, 8); SET_SCALAR_PROP(FlowModelParameters, ToleranceCnvSeq,1e-2); SET_SCALAR_PROP(FlowModelParameters, ToleranceCnvRelaxedSeq, 1e9); SET_SCALAR_PROP(FlowModelParameters, TolerancePressure, 1e-4); +SET_INT_PROP(FlowModelParameters, ReusePressureSolver, 0); +SET_STRING_PROP(FlowModelParameters, SimulationType, "implicit"); +SET_STRING_PROP(FlowModelParameters, PressureSolverJson, "pressuresolver.json"); // if openMP is available, determine the number threads per process automatically. #if _OPENMP @@ -207,7 +211,12 @@ namespace Opm double tolerance_cnv_relaxed_seq_; /// pressure tolerance double tolerance_pressure_; + /// policy for reuse of pressure solver + int reuse_pressure_solver_; + // name of json file for pressure setup + std::string pressure_solver_json_; /// Construct from user parameters or defaults. + BlackoilModelParametersEbos() { dbhp_max_rel_= EWOMS_GET_PARAM(TypeTag, Scalar, DbhpMaxRel); @@ -244,6 +253,9 @@ namespace Opm tolerance_cnv_seq_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnvSeq); tolerance_cnv_relaxed_seq_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnvRelaxedSeq); tolerance_pressure_ = EWOMS_GET_PARAM(TypeTag, Scalar, TolerancePressure); + reuse_pressure_solver_ = EWOMS_GET_PARAM(TypeTag, int, ReusePressureSolver); + pressure_solver_json_ = EWOMS_GET_PARAM(TypeTag, std::string,PressureSolverJson); + } @@ -276,13 +288,16 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, UseUpdateStabilization, "Try to detect and correct oscillations or stagnation during the Newton method"); EWOMS_REGISTER_PARAM(TypeTag, bool, MatrixAddWellContributions, "Explicitly specify the influences of wells between cells in the Jacobian and preconditioner matrices"); EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheck, "Enable the well operability checking"); - EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationType, "Define simulation type: default implicit"); + // for sequential EWOMS_REGISTER_PARAM(TypeTag, int, MaxStrictIterSeq, "Maximum number of sequential iterations before relaxed tolerances are used for the CNV convergence criterion"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceCnvSeq, "Local convergence tolerance for sequentialtransport(Maximum of local saturation errors)"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceCnvRelaxedSeq, "Relaxed local convergence tolerance that applies for iterations after the iterations with the strict tolerance for sequential transport"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, TolerancePressure, "Relaxed local convergence tolerance that applies for iterations after the iterations with the strict tolerance for sequential transport"); + EWOMS_REGISTER_PARAM(TypeTag, int, ReusePressureSolver, "Reuse policy for pressure"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, PressureSolverJson, "Json setupfile for pressure solver"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationType, "Define simulation type: default implicit"); } };