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

Merge pull request #2498 from blattms/cherry-pick-hnil-flexible-clean-interface-rebased

Browse Source 
Cleaned up flexible solver improvements.
This commit is contained in:
Atgeirr Flø Rasmussen
2020-04-02 15:28:10 +02:00
committed by GitHub
parent 74368c5f98 a8184df3c2
commit 0c8ce54b71
16 changed files with 513 additions and 234 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@@ -28,7 +28,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/timestepping/SimulatorReport.cpp
opm/simulators/flow/MissingFeatures.cpp
opm/simulators/linalg/ExtractParallelGridInformationToISTL.cpp
opm/simulators/linalg/setupPropertyTree.cpp
opm/simulators/timestepping/TimeStepControl.cpp
opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp
opm/simulators/timestepping/SimulatorTimer.cpp
@@ -169,6 +168,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/findOverlapRowsAndColumns.hpp
opm/simulators/linalg/getQuasiImpesWeights.hpp
opm/simulators/linalg/setupPropertyTree.hpp
opm/simulators/linalg/setupPropertyTree_impl.hpp
opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp
opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp
opm/simulators/timestepping/ConvergenceReport.hpp

1
flow/flow_blackoil_dunecpr.cpp
View File

@@ -55,6 +55,7 @@ SET_INT_PROP(EclFlowProblemSimple, CprMaxEllIter, 1);
SET_INT_PROP(EclFlowProblemSimple, CprEllSolvetype, 3);
SET_INT_PROP(EclFlowProblemSimple, CprReuseSetup, 3);
SET_INT_PROP(EclFlowProblemSimple, CprSolverVerbose, 0);
SET_STRING_PROP(EclFlowProblemSimple, LinearSolverConfiguration, "ilu0");
SET_STRING_PROP(EclFlowProblemSimple, SystemStrategy, "quasiimpes");
END_PROPERTIES

65
opm/simulators/linalg/FlexibleSolver.hpp
View File

@@ -33,6 +33,17 @@
namespace Dune
{
template<class C>
struct IsComm : std::false_type
{};
template<>
struct IsComm<Dune::Amg::SequentialInformation> : std::true_type
{};
template<class Index>
struct IsComm<Dune::OwnerOverlapCopyCommunication<Index>> : std::true_type
{};
/// A solver class that encapsulates all needed objects for a linear solver
/// (operator, scalar product, iterative solver and preconditioner) and sets
@@ -46,16 +57,27 @@ public:
using VectorType = VectorTypeT;
/// Create a sequential solver.
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix)
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix,
const std::function<VectorTypeT()>& weightsCalculator = std::function<VectorTypeT()>())
{
init(prm, matrix, Dune::Amg::SequentialInformation());
init(prm, matrix, weightsCalculator, Dune::Amg::SequentialInformation());
}
/// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
template <class Comm>
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix, const Comm& comm)
FlexibleSolver(const boost::property_tree::ptree& prm,
const MatrixType& matrix,
const typename std::enable_if<IsComm<Comm>::value, Comm>::type& comm)
{
init(prm, matrix, comm);
init(prm, matrix, std::function<VectorTypeT()>(), comm);
}
/// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
template <class Comm>
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix,
const std::function<VectorTypeT()>& weightsCalculator, const Comm& comm)
{
init(prm, matrix, weightsCalculator, comm);
}
virtual void apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) override
@@ -89,33 +111,41 @@ private:
// Machinery for making sequential or parallel operators/preconditioners/scalar products.
template <class Comm>
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, const Comm& comm)
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
const std::function<VectorTypeT()> weightsCalculator, const Comm& comm)
{
// Parallel case.
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
using pt = const boost::property_tree::ptree;
auto linop = std::make_shared<ParOperatorType>(matrix, comm);
linearoperator_ = linop;
auto child = prm.get_child_optional("preconditioner");
preconditioner_
= Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, prm.get_child("preconditioner"), comm);
= Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, child? *child : pt(),
weightsCalculator, comm);
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, linearoperator_->category());
}
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, const Dune::Amg::SequentialInformation&)
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
const std::function<VectorTypeT()> weightsCalculator, const Dune::Amg::SequentialInformation&)
{
// Sequential case.
using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
using pt = const boost::property_tree::ptree;
auto linop = std::make_shared<SeqOperatorType>(matrix);
linearoperator_ = linop;
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, prm.get_child("preconditioner"));
auto child = prm.get_child_optional("preconditioner");
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, child? *child : pt(),
weightsCalculator);
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
}
void initSolver(const boost::property_tree::ptree& prm)
{
const double tol = prm.get<double>("tol");
const int maxiter = prm.get<int>("maxiter");
const int verbosity = prm.get<int>("verbosity");
const std::string solver_type = prm.get<std::string>("solver");
const double tol = prm.get<double>("tol", 1e-2);
const int maxiter = prm.get<int>("maxiter", 200);
const int verbosity = prm.get<int>("verbosity", 0);
const std::string solver_type = prm.get<std::string>("solver", "bicgstab");
if (solver_type == "bicgstab") {
linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*linearoperator_,
*scalarproduct_,
@@ -131,7 +161,7 @@ private:
maxiter, // maximum number of iterations
verbosity));
} else if (solver_type == "gmres") {
int restart = prm.get<int>("restart");
int restart = prm.get<int>("restart", 15);
linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_,
*scalarproduct_,
*preconditioner_,
@@ -145,9 +175,7 @@ private:
linsolver_.reset(new Dune::UMFPack<MatrixType>(linearoperator_->getmat(), verbosity, dummy));
#endif
} else {
std::string msg("Solver not known ");
msg += solver_type;
throw std::runtime_error(msg);
OPM_THROW(std::invalid_argument, "Properties: Solver " << solver_type << " not known.");
}
}
@@ -155,9 +183,10 @@ private:
// Main initialization routine.
// Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
template <class Comm>
void init(const boost::property_tree::ptree& prm, const MatrixType& matrix, const Comm& comm)
void init(const boost::property_tree::ptree& prm, const MatrixType& matrix,
const std::function<VectorTypeT()> weightsCalculator, const Comm& comm)
{
initOpPrecSp(matrix, prm, comm);
initOpPrecSp(matrix, prm, weightsCalculator, comm);
initSolver(prm);
}

5
opm/simulators/linalg/FlowLinearSolverParameters.hpp
View File

@@ -66,6 +66,7 @@ NEW_PROP_TAG(CprUseDrs);
NEW_PROP_TAG(CprMaxEllIter);
NEW_PROP_TAG(CprEllSolvetype);
NEW_PROP_TAG(CprReuseSetup);
NEW_PROP_TAG(LinearSolverConfiguration);
NEW_PROP_TAG(LinearSolverConfigurationJsonFile);
NEW_PROP_TAG(UseGpu);
@@ -92,6 +93,7 @@ SET_BOOL_PROP(FlowIstlSolverParams, CprUseDrs, false);
SET_INT_PROP(FlowIstlSolverParams, CprMaxEllIter, 20);
SET_INT_PROP(FlowIstlSolverParams, CprEllSolvetype, 0);
SET_INT_PROP(FlowIstlSolverParams, CprReuseSetup, 0);
SET_STRING_PROP(FlowIstlSolverParams, LinearSolverConfiguration, "ilu0");
SET_STRING_PROP(FlowIstlSolverParams, LinearSolverConfigurationJsonFile, "none");
SET_BOOL_PROP(FlowIstlSolverParams, UseGpu, false);
@@ -164,6 +166,7 @@ namespace Opm
bool use_cpr_;
std::string system_strategy_;
bool scale_linear_system_;
std::string linear_solver_configuration_;
std::string linear_solver_configuration_json_file_;
bool use_gpu_;
@@ -192,6 +195,7 @@ namespace Opm
cpr_max_ell_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxEllIter);
cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype);
cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup);
linear_solver_configuration_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolverConfiguration);
linear_solver_configuration_json_file_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolverConfigurationJsonFile);
use_gpu_ = EWOMS_GET_PARAM(TypeTag, bool, UseGpu);
}
@@ -220,6 +224,7 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxEllIter, "MaxIterations of the elliptic pressure part of the cpr solver");
EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "Solver type of elliptic pressure solve (0: bicgstab, 1: cg, 2: only amg preconditioner)");
EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse Amg Setup");
EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolverConfiguration, "Configuration of solver valid is: ilu0 (default), cpr_quasiimpes, cpr_trueimpes or file (specified in LinearSolverConfigurationJsonFile) ");
EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolverConfigurationJsonFile, "Filename of JSON configuration for flexible linear solver system.");
EWOMS_REGISTER_PARAM(TypeTag, bool, UseGpu, "Use GPU cusparseSolver as the linear solver");
}

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

@@ -902,40 +902,10 @@ protected:
Vector getStorageWeights() const
{
Vector weights(rhs_->size());
BlockVector rhs(0.0);
rhs[pressureVarIndex] = 1.0;
int index = 0;
ElementContext elemCtx(simulator_);
const auto& vanguard = simulator_.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
Dune::FieldVector<Evaluation, numEq> storage;
unsigned threadId = ThreadManager::threadId();
simulator_.model().localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0);
Scalar extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor();
Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor;
Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize();
MatrixBlockType block;
double pressure_scale = 50e5;
for (int ii = 0; ii < numEq; ++ii) {
for (int jj = 0; jj < numEq; ++jj) {
block[ii][jj] = storage[ii].derivative(jj)/storage_scale;
if (jj == pressureVarIndex) {
block[ii][jj] *= pressure_scale;
}
}
}
BlockVector bweights;
MatrixBlockType block_transpose = Opm::transposeDenseMatrix(block);
block_transpose.solve(bweights, rhs);
bweights /= 1000.0; // given normal densities this scales weights to about 1.
weights[index] = bweights;
++index;
}
Opm::Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
elemCtx, simulator_.model(),
ThreadManager::threadId());
return weights;
}

69
opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
View File

@@ -25,6 +25,10 @@
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/simulators/linalg/setupPropertyTree.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <boost/property_tree/json_parser.hpp>
#include <memory>
#include <utility>
@@ -65,6 +69,19 @@ class ISTLSolverEbosFlexible
#endif
using SolverType = Dune::FlexibleSolver<MatrixType, VectorType>;
// for quasiImpesWeights
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
//typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) WellModel;
//typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
typedef typename Vector::block_type BlockVector;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
public:
static void registerParameters()
@@ -76,7 +93,7 @@ public:
: simulator_(simulator)
{
parameters_.template init<TypeTag>();
prm_ = setupPropertyTree(parameters_);
prm_ = setupPropertyTree<TypeTag>(parameters_);
extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
// For some reason simulator_.model().elementMapper() is not initialized at this stage
// Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
@@ -136,13 +153,49 @@ public:
// Never recreate solver.
}
std::function<VectorType()> 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 =
[&mat, this, transpose, pressureIndex](){
return Opm::Amg::getQuasiImpesWeights<MatrixType,
VectorType>(
mat.istlMatrix(),
pressureIndex,
transpose);
};
}else if(weightsType == "trueimpes" ){
weightsCalculator =
[this, &b, pressureIndex](){
return this->getTrueImpesWeights(b, pressureIndex);
};
}else{
OPM_THROW(std::invalid_argument, "Weights type " << weightsType << "not implemented for cpr."
<< " Please use quasiimpes or trueimpes.");
}
}
if (recreate_solver || !solver_) {
if (isParallel()) {
#if HAVE_MPI
solver_.reset(new SolverType(prm_, mat.istlMatrix(), *comm_));
solver_.reset(new SolverType(prm_, mat.istlMatrix(), weightsCalculator, *comm_));
#endif
} else {
solver_.reset(new SolverType(prm_, mat.istlMatrix()));
solver_.reset(new SolverType(prm_, mat.istlMatrix(), weightsCalculator));
}
rhs_ = b;
} else {
@@ -182,6 +235,7 @@ public:
}
protected:
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1.0).
void makeOverlapRowsInvalid(MatrixType& matrix) const
@@ -204,6 +258,15 @@ protected:
}
}
VectorType getTrueImpesWeights(const VectorType& b,const int pressureVarIndex)
{
VectorType weights(b.size());
ElementContext elemCtx(simulator_);
Opm::Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
elemCtx, simulator_.model(),
ThreadManager::threadId());
return weights;
}
const Simulator& simulator_;

77
opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
View File

@@ -26,6 +26,8 @@
#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
#include <opm/simulators/linalg/twolevelmethodcpr.hh>
#include <opm/common/ErrorMacros.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/paamg/amg.hh>
@@ -53,6 +55,21 @@ namespace Dune
template <class MatrixTypeT, class VectorTypeT>
class FlexibleSolver;
template <typename T, typename A, int i>
std::ostream& operator<<(std::ostream& out,
const BlockVector< FieldVector< T, i >, A >& vector)
{
Dune::writeMatrixMarket(vector, out);
return out;
}
template <typename T, typename A, int i>
std::istream& operator>>(std::istream& input,
BlockVector< FieldVector< T, i >, A >& vector)
{
Dune::readMatrixMarket(vector, input);
return input;
}
/// A version of the two-level preconditioner that is:
/// - Self-contained, because it owns its policy components.
@@ -69,51 +86,59 @@ public:
using MatrixType = typename OperatorType::matrix_type;
using PrecFactory = Opm::PreconditionerFactory<OperatorType, Communication>;
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm)
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm,
const std::function<VectorType()> weightsCalculator)
: linear_operator_(linearoperator)
, finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother")))
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt()))
, comm_(nullptr)
, weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(dummy_comm_, weights_, prm.get<int>("pressure_var_index"))
, coarseSolverPolicy_(prm.get_child("coarsesolver"))
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
coarseSolverPolicy_,
transpose ? 1 : 0,
transpose ? 0 : 1)
prm.get<int>("pre_smooth", transpose? 1 : 0),
prm.get<int>("post_smooth", transpose? 0 : 1))
, prm_(prm)
{
if (prm.get<int>("verbosity") > 10) {
std::ofstream outfile(prm.get<std::string>("weights_filename"));
if (prm.get<int>("verbosity", 0) > 10) {
std::string filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
std::ofstream outfile(filename);
if (!outfile) {
throw std::runtime_error("Could not write weights");
OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
}
Dune::writeMatrixMarket(weights_, outfile);
}
}
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm, const Communication& comm)
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm,
const std::function<VectorType()> weightsCalculator, const Communication& comm)
: linear_operator_(linearoperator)
, finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"), comm))
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt(), comm))
, comm_(&comm)
, weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
, levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index"))
, coarseSolverPolicy_(prm.get_child("coarsesolver"))
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index", 1))
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
coarseSolverPolicy_,
transpose ? 1 : 0,
transpose ? 0 : 1)
prm.get<int>("pre_smooth", transpose? 1 : 0),
prm.get<int>("post_smooth", transpose? 0 : 1))
, prm_(prm)
{
if (prm.get<int>("verbosity") > 10) {
std::ofstream outfile(prm.get<std::string>("weights_filename"));
if (prm.get<int>("verbosity", 0) > 10 && comm.communicator().rank() == 0) {
auto filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
std::ofstream outfile(filename);
if (!outfile) {
throw std::runtime_error("Could not write weights");
OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
}
Dune::writeMatrixMarket(weights_, outfile);
}
@@ -136,8 +161,7 @@ public:
virtual void update() override
{
Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
linear_operator_.getmat(), prm_.get<int>("pressure_var_index"), transpose, weights_);
weights_ = weightsCalculator_();
updateImpl(comm_);
}
@@ -167,20 +191,23 @@ private:
void updateImpl(const Comm*)
{
// Parallel case.
finesmoother_ = PrecFactory::create(linear_operator_, prm_.get_child("finesmoother"), *comm_);
auto child = prm_.get_child_optional("finesmoother");
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt(), *comm_);
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
}
void updateImpl(const Dune::Amg::SequentialInformation*)
{
// Serial case.
finesmoother_ = PrecFactory::create(linear_operator_, prm_.get_child("finesmoother"));
auto child = prm_.get_child_optional("finesmoother");
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt());
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
}
const OperatorType& linear_operator_;
std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> finesmoother_;
const Communication* comm_;
std::function<VectorType()> weightsCalculator_;
VectorType weights_;
LevelTransferPolicy levelTransferPolicy_;
CoarseSolverPolicy coarseSolverPolicy_;

231
opm/simulators/linalg/PreconditionerFactory.hpp
View File

@@ -29,6 +29,7 @@
#include <opm/simulators/linalg/amgcpr.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
#include <dune/istl/paamg/fastamg.hh>
#include <dune/istl/preconditioners.hh>
@@ -57,16 +58,30 @@ public:
using PrecPtr = std::shared_ptr<Dune::PreconditionerWithUpdate<Vector, Vector>>;
/// The type of creator functions passed to addCreator().
using Creator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&)>;
using ParCreator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&, const Comm&)>;
using Creator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&, const std::function<Vector()>&)>;
using ParCreator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&, const std::function<Vector()>&, const Comm&)>;
/// Create a new serial preconditioner and return a pointer to it.
/// \param op operator to be preconditioned.
/// \param prm parameters for the preconditioner, in particular its type.
/// \param weightsCalculator Calculator for weights used in CPR.
/// \return (smart) pointer to the created preconditioner.
static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm)
static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm,
const std::function<Vector()>& weightsCalculator = std::function<Vector()>())
{
return instance().doCreate(op, prm);
return instance().doCreate(op, prm, weightsCalculator);
}
/// Create a new parallel preconditioner and return a pointer to it.
/// \param op operator to be preconditioned.
/// \param prm parameters for the preconditioner, in particular its type.
/// \param comm communication object (typically OwnerOverlapCopyCommunication).
/// \param weightsCalculator Calculator for weights used in CPR.
/// \return (smart) pointer to the created preconditioner.
static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm,
const std::function<Vector()>& weightsCalculator, const Comm& comm)
{
return instance().doCreate(op, prm, weightsCalculator, comm);
}
/// Create a new parallel preconditioner and return a pointer to it.
@@ -76,9 +91,8 @@ public:
/// \return (smart) pointer to the created preconditioner.
static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm, const Comm& comm)
{
return instance().doCreate(op, prm, comm);
return instance().doCreate(op, prm, std::function<Vector()>(), comm);
}
/// Add a creator for a serial preconditioner to the PreconditionerFactory.
/// After the call, the user may obtain a preconditioner by
/// calling create() with the given type string as a parameter
@@ -111,13 +125,15 @@ private:
// Helpers for creation of AMG preconditioner.
static Criterion amgCriterion(const boost::property_tree::ptree& prm)
{
Criterion criterion(15, prm.get<int>("coarsenTarget"));
Criterion criterion(15, prm.get<int>("coarsenTarget", 1200));
criterion.setDefaultValuesIsotropic(2);
criterion.setAlpha(prm.get<double>("alpha"));
criterion.setBeta(prm.get<double>("beta"));
criterion.setMaxLevel(prm.get<int>("maxlevel"));
criterion.setSkipIsolated(false);
criterion.setDebugLevel(prm.get<int>("verbosity"));
criterion.setAlpha(prm.get<double>("alpha", 0.33));
criterion.setBeta(prm.get<double>("beta", 1e-5));
criterion.setMaxLevel(prm.get<int>("maxlevel", 15));
criterion.setSkipIsolated(prm.get<bool>("skip_isolated", false));
criterion.setNoPreSmoothSteps(prm.get<int>("pre_smooth", 1));
criterion.setNoPostSmoothSteps(prm.get<int>("post_smooth", 1));
criterion.setDebugLevel(prm.get<int>("verbosity", 0));
return criterion;
}
@@ -126,20 +142,30 @@ private:
{
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
SmootherArgs smootherArgs;
smootherArgs.iterations = prm.get<int>("iterations");
smootherArgs.iterations = prm.get<int>("iterations", 1);
// smootherArgs.overlap=SmootherArgs::vertex;
// smootherArgs.overlap=SmootherArgs::none;
// smootherArgs.overlap=SmootherArgs::aggregate;
smootherArgs.relaxationFactor = prm.get<double>("relaxation");
smootherArgs.relaxationFactor = prm.get<double>("relaxation", 0.9);
return smootherArgs;
}
template <class Smoother>
static PrecPtr makeAmgPreconditioner(const Operator& op, const boost::property_tree::ptree& prm)
static PrecPtr makeAmgPreconditioner(const Operator& op, const boost::property_tree::ptree& prm, bool useKamg = false)
{
auto crit = amgCriterion(prm);
auto sargs = amgSmootherArgs<Smoother>(prm);
return std::make_shared<Dune::Amg::AMGCPR<Operator, Vector, Smoother>>(op, crit, sargs);
if(useKamg){
return std::make_shared<
Dune::DummyUpdatePreconditioner<
Dune::Amg::KAMG< Operator, Vector, Smoother>
>
>(op, crit, sargs,
prm.get<size_t>("max_krylov", 1),
prm.get<double>("min_reduction", 1e-1) );
}else{
return std::make_shared<Dune::Amg::AMGCPR<Operator, Vector, Smoother>>(op, crit, sargs);
}
}
// Add a useful default set of preconditioners to the factory.
@@ -154,61 +180,63 @@ private:
using V = Vector;
using P = boost::property_tree::ptree;
using C = Comm;
doAddCreator("ILU0", [](const O& op, const P& prm, const C& comm) {
const double w = prm.get<double>("relaxation");
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const double w = prm.get<double>("relaxation", 1.0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) {
const double w = prm.get<double>("relaxation");
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const double w = prm.get<double>("relaxation", 1.0);
const int n = prm.get<int>("ilulevel", 0);
// Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner.
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ILUn", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation");
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const int n = prm.get<int>("ilulevel", 0);
const double w = prm.get<double>("relaxation", 1.0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("Jac", [](const O& op, const P& prm, const std::function<Vector()>&,
const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqJac<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("GS", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqGS<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("SOR", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSOR<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("SSOR", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("amg", [](const O& op, const P& prm, const C& comm) {
const std::string smoother = prm.get<std::string>("smoother");
if (smoother == "ILU0") {
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
auto crit = amgCriterion(prm);
auto sargs = amgSmootherArgs<Smoother>(prm);
return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
} else {
std::string msg("No such smoother: ");
msg += smoother;
throw std::runtime_error(msg);
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
}
});
doAddCreator("cpr", [](const O& op, const P& prm, const C& comm) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, comm);
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
assert(weightsCalculator);
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, comm);
});
doAddCreator("cprt", [](const O& op, const P& prm, const C& comm) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, comm);
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
assert(weightsCalculator);
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, weightsCalculator, comm);
});
}
@@ -221,45 +249,46 @@ private:
using M = Matrix;
using V = Vector;
using P = boost::property_tree::ptree;
doAddCreator("ILU0", [](const O& op, const P& prm) {
const double w = prm.get<double>("relaxation");
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&) {
const double w = prm.get<double>("relaxation", 1.0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ParOverILU0", [](const O& op, const P& prm) {
const double w = prm.get<double>("relaxation");
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ILUn", [](const O& op, const P& prm) {
const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation");
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&) {
const double w = prm.get<double>("relaxation", 1.0);
const int n = prm.get<int>("ilulevel", 0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("Jac", [](const O& op, const P& prm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&) {
const int n = prm.get<int>("ilulevel", 0);
const double w = prm.get<double>("relaxation", 1.0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("Jac", [](const O& op, const P& prm, const std::function<Vector()>&) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqJac<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("GS", [](const O& op, const P& prm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqGS<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("SOR", [](const O& op, const P& prm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqSOR<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("SSOR", [](const O& op, const P& prm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqSSOR<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("amg", [](const O& op, const P& prm) {
const std::string smoother = prm.get<std::string>("smoother");
if (smoother == "ILU0") {
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
using Smoother = SeqILU<M, V, V>;
#else
@@ -283,23 +312,45 @@ private:
#endif
return makeAmgPreconditioner<Smoother>(op, prm);
} else {
std::string msg("No such smoother: ");
msg += smoother;
throw std::runtime_error(msg);
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
}
});
doAddCreator("famg", [](const O& op, const P& prm) {
doAddCreator("kamg", [](const O& op, const P& prm, const std::function<Vector()>&) {
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
using Smoother = SeqILU0<M, V, V>;
return makeAmgPreconditioner<Smoother>(op, prm, true);
} else if (smoother == "Jac") {
using Smoother = SeqJac<M, V, V>;
return makeAmgPreconditioner<Smoother>(op, prm, true);
} else if (smoother == "SOR") {
using Smoother = SeqSOR<M, V, V>;
return makeAmgPreconditioner<Smoother>(op, prm, true);
// } else if (smoother == "GS") {
// using Smoother = SeqGS<M, V, V>;
// return makeAmgPreconditioner<Smoother>(op, prm, true);
} else if (smoother == "SSOR") {
using Smoother = SeqSSOR<M, V, V>;
return makeAmgPreconditioner<Smoother>(op, prm, true);
} else if (smoother == "ILUn") {
using Smoother = SeqILUn<M, V, V>;
return makeAmgPreconditioner<Smoother>(op, prm, true);
} else {
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
}
});
doAddCreator("famg", [](const O& op, const P& prm, const std::function<Vector()>&) {
auto crit = amgCriterion(prm);
Dune::Amg::Parameters parms;
parms.setNoPreSmoothSteps(1);
parms.setNoPostSmoothSteps(1);
return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
});
doAddCreator("cpr", [](const O& op, const P& prm) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm);
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator);
});
doAddCreator("cprt", [](const O& op, const P& prm) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(op, prm);
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(op, prm, weightsCalculator);
});
}
@@ -320,9 +371,10 @@ private:
}
// Actually creates the product object.
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm)
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm,
const std::function<Vector()> weightsCalculator)
{
const std::string& type = prm.get<std::string>("type");
const std::string& type = prm.get<std::string>("type", "ParOverILU0");
auto it = creators_.find(type);
if (it == creators_.end()) {
std::ostringstream msg;
@@ -331,14 +383,15 @@ private:
msg << prec.first << ' ';
}
msg << std::endl;
throw std::runtime_error(msg.str());
OPM_THROW(std::invalid_argument, msg.str());
}
return it->second(op, prm);
return it->second(op, prm, weightsCalculator);
}
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm, const Comm& comm)
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm,
const std::function<Vector()> weightsCalculator, const Comm& comm)
{
const std::string& type = prm.get<std::string>("type");
const std::string& type = prm.get<std::string>("type", "ParOverILU0");
auto it = parallel_creators_.find(type);
if (it == parallel_creators_.end()) {
std::ostringstream msg;
@@ -348,9 +401,9 @@ private:
msg << prec.first << ' ';
}
msg << std::endl;
throw std::runtime_error(msg.str());
OPM_THROW(std::invalid_argument, msg.str());
}
return it->second(op, prm, comm);
return it->second(op, prm, weightsCalculator, comm);
}
// Actually adds the creator.

2
opm/simulators/linalg/PreconditionerWithUpdate.hpp
View File

@@ -22,7 +22,7 @@
#include <dune/istl/preconditioner.hh>
#include <memory>
#include <boost/property_tree/ptree.hpp>
namespace Dune
{

4
opm/simulators/linalg/PressureSolverPolicy.hpp
View File

@@ -61,9 +61,9 @@ namespace Amg
: linsolver_()
{
if (op.category() == Dune::SolverCategory::overlapping) {
linsolver_.reset(new Solver(prm, op.getmat(), comm));
linsolver_.reset(new Solver(prm, op.getmat(), std::function<X()>(), comm));
} else {
linsolver_.reset(new Solver(prm, op.getmat()));
linsolver_.reset(new Solver(prm, op.getmat(), std::function<X()>()));
}
}

41
opm/simulators/linalg/getQuasiImpesWeights.hpp
View File

@@ -83,6 +83,47 @@ namespace Amg
return weights;
}
template<class Vector, class GridView, class ElementContext, class Model>
void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const GridView& gridView,
ElementContext& elemCtx, const Model& model, std::size_t threadId)
{
using VectorBlockType = typename Vector::block_type;
using Matrix = typename std::decay_t<decltype(model.linearizer().jacobian())>;
using MatrixBlockType = typename Matrix::MatrixBlock;
constexpr int numEq = VectorBlockType::size();
using Evaluation = typename std::decay_t<decltype(model.localLinearizer(threadId).localResidual().residual(0))>
::block_type;
VectorBlockType rhs(0.0);
rhs[pressureVarIndex] = 1.0;
int index = 0;
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
Dune::FieldVector<Evaluation, numEq> storage;
model.localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0);
auto extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor();
auto scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor;
auto storage_scale = scvVolume / elemCtx.simulator().timeStepSize();
MatrixBlockType block;
double pressure_scale = 50e5;
for (int ii = 0; ii < numEq; ++ii) {
for (int jj = 0; jj < numEq; ++jj) {
block[ii][jj] = storage[ii].derivative(jj)/storage_scale;
if (jj == pressureVarIndex) {
block[ii][jj] *= pressure_scale;
}
}
}
VectorBlockType bweights;
MatrixBlockType block_transpose = Details::transposeDenseMatrix(block);
block_transpose.solve(bweights, rhs);
bweights /= 1000.0; // given normal densities this scales weights to about 1.
weights[index] = bweights;
++index;
}
}
} // namespace Amg
} // namespace Opm

55
opm/simulators/linalg/setupPropertyTree.cpp
View File

@@ -1,55 +0,0 @@
/*
Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/linalg/setupPropertyTree.hpp>
#include <boost/version.hpp>
#include <boost/property_tree/json_parser.hpp>
namespace Opm
{
/// Set up a property tree intended for FlexibleSolver by either reading
/// the tree from a JSON file or creating a tree giving the default solver
/// and preconditioner. If the latter, the parameters --linear-solver-reduction,
/// --linear-solver-maxiter and --linear-solver-verbosity are used, but if reading
/// from file the data in the JSON file will override any other options.
boost::property_tree::ptree
setupPropertyTree(const FlowLinearSolverParameters& p)
{
boost::property_tree::ptree prm;
#if BOOST_VERSION / 100 % 1000 > 48
if (p.linear_solver_configuration_json_file_ != "none") {
boost::property_tree::read_json(p.linear_solver_configuration_json_file_, prm);
} else
#endif
{
prm.put("tol", p.linear_solver_reduction_);
prm.put("maxiter", p.linear_solver_maxiter_);
prm.put("verbosity", p.linear_solver_verbosity_);
prm.put("solver", "bicgstab");
prm.put("preconditioner.type", "ParOverILU0");
prm.put("preconditioner.relaxation", 1.0);
}
return prm;
}
} // namespace Opm

3
opm/simulators/linalg/setupPropertyTree.hpp
View File

@@ -27,8 +27,11 @@
namespace Opm
{
template<class TypeTag>
boost::property_tree::ptree setupPropertyTree(const FlowLinearSolverParameters& p);
} // namespace Opm
#include "setupPropertyTree_impl.hpp"
#endif // OPM_SETUPPROPERTYTREE_HEADER_INCLUDED

118
opm/simulators/linalg/setupPropertyTree_impl.hpp Normal file
View File

@@ -0,0 +1,118 @@
/*
Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/linalg/setupPropertyTree.hpp>
#include <boost/version.hpp>
#include <boost/property_tree/json_parser.hpp>
namespace Opm
{
/// Set up a property tree intended for FlexibleSolver by either reading
/// the tree from a JSON file or creating a tree giving the default solver
/// and preconditioner. If the latter, the parameters --linear-solver-reduction,
/// --linear-solver-maxiter and --linear-solver-verbosity are used, but if reading
/// from file the data in the JSON file will override any other options.
template<class TypeTag>
boost::property_tree::ptree
setupPropertyTree(const FlowLinearSolverParameters& p)
{
boost::property_tree::ptree prm;
if (p.linear_solver_configuration_ == "file") {
#if BOOST_VERSION / 100 % 1000 > 48
if (p.linear_solver_configuration_json_file_ == "none"){
OPM_THROW(std::invalid_argument,
"--linear-solver-configuration=file requires that a filename "
<< "is passed with "
<< "--linear-solver-configuration-json-file=filename.");
}else{
boost::property_tree::read_json(p.linear_solver_configuration_json_file_, prm);
}
#else
OPM_THROW(std::invalid_argument,
"--linear-solver-configuration=file not supported with "
<< "boost version. Needs versoin > 1.48.");
#endif
}
else
{
std::string conf = p.linear_solver_configuration_;
// Support old UseCpr if not configuration was set
if (!EWOMS_PARAM_IS_SET(TypeTag, std::string, LinearSolverConfiguration) && p.use_cpr_)
{
conf = "cpr_quasiimpes";
}
if((conf == "cpr_trueimpes") || (conf == "cpr_quasiimpes")){
prm.put("tol", p.linear_solver_reduction_);
if (EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter))
prm.put("maxiter", p.linear_solver_maxiter_);// Trust that the user knows what he does
else
prm.put("maxiter", 20); // Use our own default.
prm.put("verbosity", p.linear_solver_verbosity_);
prm.put("solver", "bicgstab");
prm.put("preconditioner.type", "cpr");
prm.put("preconditioner.weight_filename", "cpr_weights.txt");
prm.put("preconditioner.weight_type","quasiimpes");
prm.put("preconditioner.finesmoother.type", "ParOverILU0");
prm.put("preconditioner.finesmoother.relaxation", 1.0);
prm.put("preconditioner.pressure_var_index",1);
prm.put("preconditioner.verbosity",0);
prm.put("preconditioner.coarsesolver.maxiter",1);
prm.put("preconditioner.coarsesolver.tol",1e-1);
prm.put("preconditioner.coarsesolver.solver","loopsolver");
prm.put("preconditioner.coarsesolver.verbosity",0);
prm.put("preconditioner.coarsesolver.preconditioner.type","amg");
prm.put("preconditioner.coarsesolver.preconditioner.alpha",0.333333333333);
prm.put("preconditioner.coarsesolver.preconditioner.relaxation",1.0);
if (EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter))
prm.put("preconditioner.coarsesolver.preconditioner.iterations", p.cpr_max_ell_iter_);
prm.put("preconditioner.coarsesolver.preconditioner.iterations",1);
prm.put("preconditioner.coarsesolver.preconditioner.coarsenTarget",1200);
prm.put("preconditioner.coarsesolver.preconditioner.pre_smooth",1);
prm.put("preconditioner.coarsesolver.preconditioner.post_smooth",1);
prm.put("preconditioner.coarsesolver.preconditioner.beta",1e-5);
prm.put("preconditioner.coarsesolver.preconditioner.smoother","ILU0");
prm.put("preconditioner.coarsesolver.preconditioner.verbosity",0);
prm.put("preconditioner.coarsesolver.preconditioner.maxlevel",15);
prm.put("preconditioner.coarsesolver.preconditioner.skip_isolated",0);
if(p.linear_solver_configuration_ == "cpr_trueimpes"){
prm.put("preconditioner.weight_type","trueimpes");
}
} else {
if(conf != "ilu0"){
OPM_THROW(std::invalid_argument, conf << "is not a valid setting for --linear-solver-configuration."
<< " Please use ilu0, cpr_trueimpes, or cpr_quasiimpes");
}
prm.put("tol", p.linear_solver_reduction_);
prm.put("maxiter", p.linear_solver_maxiter_);
prm.put("verbosity", p.linear_solver_verbosity_);
prm.put("solver", "bicgstab");
prm.put("preconditioner.type", "ParOverILU0");
prm.put("preconditioner.relaxation", p.ilu_relaxation_);
prm.put("preconditioner.ilulevel", p.ilu_fillin_level_);
}
}
return prm;
}
} // namespace Opm

14
tests/test_flexiblesolver.cpp
View File

@@ -58,7 +58,19 @@ testSolver(const boost::property_tree::ptree& prm, const std::string& matrix_fil
}
readMatrixMarket(rhs, rhsfile);
}
Dune::FlexibleSolver<Matrix, Vector> solver(prm, matrix);
bool transpose = false;
if(prm.get<std::string>("preconditioner.type") == "cprt"){
transpose = true;
}
auto wc = [&matrix, &prm, transpose]()
{
return Opm::Amg::getQuasiImpesWeights<Matrix,
Vector>(matrix,
prm.get<int>("preconditioner.pressure_var_index"),
transpose);
};
Dune::FlexibleSolver<Matrix, Vector> solver(prm, matrix, wc);
Vector x(rhs.size());
Dune::InverseOperatorResult res;
solver.apply(x, rhs, res);

24
tests/test_preconditionerfactory.cpp
View File

@@ -93,7 +93,19 @@ testPrec(const boost::property_tree::ptree& prm, const std::string& matrix_filen
using Operator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
Operator op(matrix);
using PrecFactory = Opm::PreconditionerFactory<Operator>;
auto prec = PrecFactory::create(op, prm.get_child("preconditioner"));
bool transpose = false;
if(prm.get<std::string>("preconditioner.type") == "cprt"){
transpose = true;
}
auto wc = [&matrix, &prm, transpose]()
{
return Opm::Amg::getQuasiImpesWeights<Matrix,
Vector>(matrix,
prm.get<int>("preconditioner.pressure_var_index"),
transpose);
};
auto prec = PrecFactory::create(op, prm.get_child("preconditioner"), wc);
Dune::BiCGSTABSolver<Vector> solver(op, *prec, prm.get<double>("tol"), prm.get<int>("maxiter"), prm.get<int>("verbosity"));
Vector x(rhs.size());
Dune::InverseOperatorResult res;
@@ -183,18 +195,18 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner)
// Test with 1x1 block solvers.
{
const int bz = 1;
BOOST_CHECK_THROW(testPrec<bz>(prm, "matr33.txt", "rhs3.txt"), std::runtime_error);
BOOST_CHECK_THROW(testPrec<bz>(prm, "matr33.txt", "rhs3.txt"), std::invalid_argument);
}
// Test with 3x3 block solvers.
{
const int bz = 3;
BOOST_CHECK_THROW(testPrec<bz>(prm, "matr33.txt", "rhs3.txt"), std::runtime_error);
BOOST_CHECK_THROW(testPrec<bz>(prm, "matr33.txt", "rhs3.txt"), std::invalid_argument);
}
// Add preconditioner to factory for block size 1.
PF<1>::addCreator("nothing", [](const O<1>&, const pt::ptree&) {
PF<1>::addCreator("nothing", [](const O<1>&, const pt::ptree&, const std::function<V<1>()>&) {
return Dune::wrapPreconditioner<NothingPreconditioner<V<1>>>();
});
@@ -205,11 +217,11 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner)
// Test with 3x3 block solvers.
{
const int bz = 3;
BOOST_CHECK_THROW(testPrec<bz>(prm, "matr33.txt", "rhs3.txt"), std::runtime_error);
BOOST_CHECK_THROW(testPrec<bz>(prm, "matr33.txt", "rhs3.txt"), std::invalid_argument);
}
// Add preconditioner to factory for block size 3.
PF<3>::addCreator("nothing", [](const O<3>&, const pt::ptree&) {
PF<3>::addCreator("nothing", [](const O<3>&, const pt::ptree&, const std::function<V<3>()>&) {
return Dune::wrapPreconditioner<NothingPreconditioner<V<3>>>();
});