OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Make flexible solvers and CPR accessible for flow binary.

Browse Source 
Previously, one had to call a seperate binary called
flow_blackoil_dunecpr. Unforntunately, that was only built if users
requested that the tests (`-DBUILD_TESTING=ON`) and the flow
variants (`-DBUILD_FLOW_VARIANTS=ON`) should be built. In addition
it would use a slightly different nonlinear solver implementation.

With this commit flow can be asked to either use
CPR (`--use-cpr=true --matrix-add-well-contributions=true`) or to use
the flexible solvers bx setting the `--linear-solver-configuration`
option. In all other cases the usual solver implementations are still
used.

Note that the flexible solvers still need
`--matrix-add-well-contributions=true` and hence cannot cope with
multi-segment wells.
This commit is contained in:
Markus Blatt 2020-05-11 09:50:28 +02:00
parent 71b234ec54
commit 47acb6ecb3
1 changed files with 180 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

254
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -29,6 +29,8 @@
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
#include <opm/simulators/linalg/setupPropertyTree.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
@ -126,25 +128,9 @@ public:
WellModelMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod,
const std::any& parallelInformation OPM_UNUSED_NOMPI = std::any() )
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_()
{
#if HAVE_MPI
if( parallelInformation.type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
std::any_cast<const ParallelISTLInformation&>( parallelInformation);
comm_.reset( new communication_type( info.communicator() ) );
}
#endif
}
WellModelMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod,
std::shared_ptr<communication_type> comm )
const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >())
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_(comm)
{
}
{}
virtual void apply( const X& x, Y& y ) const override
@ -176,11 +162,6 @@ public:
virtual const matrix_type& getmat() const override { return A_for_precond_; }
std::shared_ptr<communication_type> comm()
{
return comm_;
}
protected:
const matrix_type& A_ ;
const matrix_type& A_for_precond_ ;
@ -222,19 +203,9 @@ public:
WellModelGhostLastMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod,
const size_t interiorSize,
const std::any& parallelInformation OPM_UNUSED_NOMPI = std::any() )
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), interiorSize_(interiorSize), comm_()
{
#if HAVE_MPI
if( parallelInformation.type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
std::any_cast<const ParallelISTLInformation&>( parallelInformation);
comm_.reset( new communication_type( info.communicator() ) );
}
#endif
}
const size_t interiorSize )
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), interiorSize_(interiorSize)
{}
virtual void apply( const X& x, Y& y ) const override
{
@ -269,11 +240,6 @@ public:
virtual const matrix_type& getmat() const override { return A_for_precond_; }
communication_type* comm()
{
return comm_.operator->();
}
protected:
void ghostLastProject(Y& y) const
{
@ -286,8 +252,6 @@ protected:
const matrix_type& A_for_precond_ ;
const WellModel& wellMod_;
size_t interiorSize_;
std::unique_ptr< communication_type > comm_;
};
/// This class solves the fully implicit black-oil system by
@ -313,6 +277,7 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
using FlexibleSolverType = Dune::FlexibleSolver<Matrix, Vector>;
// Due to miscibility oil <-> gas the water eqn is the one we can replace with a pressure equation.
static const bool waterEnabled = Indices::waterEnabled;
static const int pindex = (waterEnabled) ? BlackOilDefaultIndexTraits::waterCompIdx : BlackOilDefaultIndexTraits::oilCompIdx;
@ -324,6 +289,12 @@ protected:
std::unique_ptr<BdaBridge> bdaBridge;
#endif
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
#else
typedef Dune::CollectiveCommunication< int > communication_type;
#endif
public:
typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
@ -340,7 +311,16 @@ protected:
iterations_( 0 ),
converged_(false)
{
#if HAVE_MPI
comm_.reset( new communication_type( simulator_.vanguard().grid().comm() ) );
#endif
parameters_.template init<TypeTag>();
useFlexible_ = parameters_.use_cpr_ || EWOMS_PARAM_IS_SET(TypeTag, std::string, LinearSolverConfiguration);
if (useFlexible_)
{
prm_ = setupPropertyTree<TypeTag>(parameters_);
}
const auto& gridForConn = simulator_.vanguard().grid();
#if HAVE_CUDA
bool use_gpu = EWOMS_GET_PARAM(TypeTag, bool, UseGpu);
@ -359,13 +339,19 @@ protected:
}
#endif
extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
const bool useWellConn = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
if (!useWellConn_ && useFlexible_)
{
OPM_THROW(std::logic_error, "Flexible solvers and CPR need the well contribution in the matrix. Please run with"
" --matrix-add-well-contributions=true");
}
ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_);
if (!ownersFirst_ || parameters_.linear_solver_use_amg_ || parameters_.use_cpr_ ) {
detail::setWellConnections(gridForConn, simulator_.vanguard().schedule().getWellsatEnd(), useWellConn, wellConnectionsGraph_);
if (!ownersFirst_ || parameters_.linear_solver_use_amg_ || useFlexible_ ) {
detail::setWellConnections(gridForConn, simulator_.vanguard().schedule().getWellsatEnd(), useWellConn_, wellConnectionsGraph_);
// For some reason simulator_.model().elementMapper() is not initialized at this stage
// Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
// Set it up manually
@ -381,6 +367,17 @@ protected:
if (ownersFirst_)
OpmLog::warning("OwnerCellsFirst option is true, but ignored.");
}
if (useFlexible_)
{
// Print parameters to PRT/DBG logs.
if (simulator.gridView().comm().rank() == 0) {
std::ostringstream os;
os << "Property tree for linear solver:\n";
boost::property_tree::write_json(os, prm_, true);
OpmLog::note(os.str());
}
}
}
// nothing to clean here
@ -390,16 +387,51 @@ protected:
void prepare(const SparseMatrixAdapter& M, Vector& b)
{
matrix_.reset(new Matrix(M.istlMatrix()));
static bool firstcall = true;
#if HAVE_MPI
if (firstcall && parallelInformation_.type() == typeid(ParallelISTLInformation)) {
// Parallel case.
const ParallelISTLInformation* parinfo = std::any_cast<ParallelISTLInformation>(&parallelInformation_);
assert(parinfo);
const size_t size = M.istlMatrix().N();
parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
}
#endif
// update matrix entries for solvers.
if (noGhostMat_)
{
copyJacToNoGhost(M.istlMatrix(), *noGhostMat_);
}
else
{
if (firstcall)
{
//Only do this once as flexible solvers wil store a reference.
matrix_.reset(new Matrix(M.istlMatrix()));
}
else
{
// Pointers stays the same for flexible solver
*matrix_ = M.istlMatrix();
}
}
rhs_ = &b;
this->scaleSystem();
if (useFlexible_)
{
prepareFlexibleSolver();
}
else
{
this->scaleSystem();
}
firstcall = false;
}
void scaleSystem()
{
const bool matrix_cont_added = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
if (matrix_cont_added) {
if (useWellConn_) {
bool form_cpr = true;
if (parameters_.system_strategy_ == "quasiimpes") {
weights_ = getQuasiImpesWeights();
@ -456,17 +488,24 @@ protected:
bool solve(Vector& x) {
// Solve system.
if (useFlexible_)
{
Dune::InverseOperatorResult res;
assert(flexibleSolver_);
flexibleSolver_->apply(x, *rhs_, res);
iterations_ = res.iterations;
return converged_ = res.converged;
}
const WellModel& wellModel = simulator_.problem().wellModel();
if( isParallel() )
{
if ( ownersFirst_ && (!parameters_.linear_solver_use_amg_ || !parameters_.use_cpr_) ) {
typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
Operator opA(*matrix_, *matrix_, wellModel, interiorCellNum_,
parallelInformation_ );
Operator opA(*matrix_, *matrix_, wellModel, interiorCellNum_);
assert( opA.comm() );
solve( opA, x, *rhs_, *(opA.comm()) );
solve( opA, x, *rhs_, *comm_ );
}
else {
@ -474,10 +513,9 @@ protected:
copyJacToNoGhost(*matrix_, *noGhostMat_);
Operator opA(*noGhostMat_, *noGhostMat_, wellModel,
parallelInformation_ );
comm_ );
assert( opA.comm() );
solve( opA, x, *rhs_, *(opA.comm()) );
solve( opA, x, *rhs_, *comm_ );
}
}
@ -575,8 +613,7 @@ protected:
bool use_gpu = bdaBridge->getUseGpu();
if (use_gpu) {
WellContributions wellContribs;
const bool addWellContribs = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
if (!addWellContribs) {
if (!useWellConn_) {
simulator_.problem().wellModel().getWellContributions(wellContribs);
}
bdaBridge->solve_system(matrix_.get(), istlb, wellContribs, result);
@ -709,14 +746,10 @@ protected:
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
const ParallelISTLInformation& info =
std::any_cast<const ParallelISTLInformation&>( parallelInformation_);
Comm istlComm(info.communicator());
// Construct operator, scalar product and vectors needed.
typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,Comm> Operator;
Operator opA(A, istlComm);
solve( opA, x, b, istlComm );
Operator opA(A, *comm_);
solve( opA, x, b, *comm_ );
}
else
#endif
@ -744,14 +777,6 @@ protected:
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
const size_t size = opA.getmat().N();
const ParallelISTLInformation& info =
std::any_cast<const ParallelISTLInformation&>( parallelInformation_);
// As we use a dune-istl with block size np the number of components
// per parallel is only one.
info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
size, 1);
// Construct operator, scalar product and vectors needed.
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
}
@ -802,6 +827,81 @@ protected:
#endif
}
void prepareFlexibleSolver()
{
// Decide if we should recreate the solver or just do
// a minimal preconditioner update.
const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
bool recreate_solver = false;
if (this->parameters_.cpr_reuse_setup_ == 0) {
// Always recreate solver.
recreate_solver = true;
} else if (this->parameters_.cpr_reuse_setup_ == 1) {
// Recreate solver on the first iteration of every timestep.
if (newton_iteration == 0) {
recreate_solver = true;
}
} else if (this->parameters_.cpr_reuse_setup_ == 2) {
// Recreate solver if the last solve used more than 10 iterations.
if (this->iterations() > 10) {
recreate_solver = true;
}
} else {
assert(this->parameters_.cpr_reuse_setup_ == 3);
assert(recreate_solver == false);
// Never recreate solver.
}
std::function<Vector()> weightsCalculator;
auto preconditionerType = prm_.get("preconditioner.type", "cpr");
if( preconditionerType == "cpr" ||
preconditionerType == "cprt"
)
{
bool transpose = false;
if(preconditionerType == "cprt"){
transpose = true;
}
auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes");
auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1);
if(weightsType == "quasiimpes") {
// weighs will be created as default in the solver
weightsCalculator =
[this, transpose, pressureIndex](){
return Opm::Amg::getQuasiImpesWeights<Matrix,
Vector>(*this->matrix_,
pressureIndex,
transpose);
};
}else if(weightsType == "trueimpes" ){
weightsCalculator =
[this](){
return this->getStorageWeights();
};
}else{
OPM_THROW(std::invalid_argument, "Weights type " << weightsType << "not implemented for cpr."
<< " Please use quasiimpes or trueimpes.");
}
}
if (recreate_solver || !flexibleSolver_) {
if (isParallel()) {
#if HAVE_MPI
flexibleSolver_.reset(new FlexibleSolverType(prm_, *matrix_, weightsCalculator, *comm_));
#endif
} else {
flexibleSolver_.reset(new FlexibleSolverType(prm_, *matrix_, weightsCalculator));
}
}
else
{
flexibleSolver_->preconditioner().update();
}
}
/// Create sparsity pattern of matrix without off-diagonal ghost entries.
void noGhostAdjacency()
{
@ -1041,16 +1141,22 @@ protected:
Vector *rhs_;
std::unique_ptr<Matrix> matrix_for_preconditioner_;
std::unique_ptr<FlexibleSolverType> flexibleSolver_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;
std::vector<std::set<int>> wellConnectionsGraph_;
bool ownersFirst_;
bool useWellConn_;
bool useFlexible_;
size_t interiorCellNum_;
FlowLinearSolverParameters parameters_;
boost::property_tree::ptree prm_;
Vector weights_;
bool scale_variables_;
std::shared_ptr< communication_type > comm_;
}; // end ISTLSolver
} // namespace Opm