Skip to content

Commit

Permalink
Introduced more settings for pressure solver and reuse of pressure ma…
Browse files Browse the repository at this point in the history 
…trix
  • Loading branch information
@hnil @GitPaean
hnil authored and GitPaean committed Aug 19, 2020
1 parent 0dc1f00 commit 0bea6ab
Show file tree
Hide file tree
Showing 2 changed files with 74 additions and 25 deletions.
78 changes: 56 additions & 22 deletions opm/simulators/flow/BlackoilModelEbos.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,6 @@ SET_BOOL_PROP(EclFlowProblem, EnableBrine, false);
SET_TYPE_PROP(EclFlowProblem, EclWellModel, Opm::BlackoilWellModel<TypeTag>);
SET_TAG_PROP(EclFlowProblem, LinearSolverSplice, FlowIstlSolver);


END_PROPERTIES

namespace Opm {
Expand Down Expand Up @@ -145,6 +144,10 @@ namespace Opm {
typedef Dune::BlockVector<VectorBlockType> BVector;

typedef ISTLSolverEbos<TypeTag> ISTLSolverType;

typedef Dune::BCRSMatrix< Dune::FieldMatrix<double,1,1>> PressureMatrixType;
typedef Dune::BlockVector<Dune::FieldVector<double,1>> PressureVectorType;
typedef Dune::FlexibleSolver<PressureMatrixType, PressureVectorType> PressureSolverType;
//typedef typename SolutionVector :: value_type PrimaryVariables ;

// --------- Public methods ---------
Expand Down Expand Up @@ -488,27 +491,29 @@ namespace Opm {
>
makePressureSolver(Dune::BCRSMatrix< Dune::FieldMatrix<double,1,1>>& 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<double,1,1>>;
using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double,1>>;
std::any parallelInformation;
extractParallelGridInformationToISTL(ebosSimulator_.vanguard().grid(), parallelInformation);
using AbstractOperatorType = Dune::AssembledLinearOperator<PressureMatrixType, PressureVectorType, PressureVectorType>;
using PressureSolverType = Dune::FlexibleSolver<PressureMatrixType, PressureVectorType>;
std::unique_ptr<AbstractOperatorType> operator_for_flexiblesolver;
std::function<PressureVectorType()> weightsCalculator;// dummy
std::unique_ptr<PressureSolverType> pressureSolver;
Expand Down Expand Up @@ -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)
Expand All @@ -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<double,1,1>>;
using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double,1>>;

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<Mat,
BVector,
PressureMatrixType,
PressureVectorType>(
ebosJac.istlMatrix(),
pressureVarIndex,
weights);
if(!pmatrix_){
// make matrix structure and fill with elements
pmatrix_ =
std::make_unique<PressureMatrixType>(PressureHelper::makePressureMatrix<Mat,
BVector,
PressureMatrixType,
PressureVectorType>(
ebosJac.istlMatrix(),
pressureVarIndex,
weights));
}else{
//only update elements of matrix
// assume matrix structure never change
PressureHelper::makePressureMatrixEntries<Mat,
PressureMatrixType,
BVector>(
*pmatrix_,
ebosJac.istlMatrix(),
pressureVarIndex,
weights);
}
PressureVectorType rhs(ebosResid.size(),0);
PressureHelper::moveToPressureEqn(ebosResid, rhs, weights);

using PressureSolverType = Dune::FlexibleSolver<PressureMatrixType, PressureVectorType>;
std::unique_ptr<PressureSolverType> 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){
Expand Down Expand Up @@ -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<PressureSolverType> pressureSolver_;
std::unique_ptr<PressureMatrixType> pmatrix_;

public:
std::vector<bool> wasSwitched_;
};
Expand Down
21 changes: 18 additions & 3 deletions opm/simulators/flow/BlackoilModelParametersEbos.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);


}

Expand Down Expand Up @@ -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");

}
};
Expand Down

0 comments on commit 0bea6ab

Please sign in to comment.