Add choice of solution strategies to simulators.

The following is changed in this commit:
 - The constructor for NewtonIterationBlackoilSimple now takes
   a parameter object instead of a linear solver.
 - The fully implicit black-oil simulators can now use the CPR
   preconditioning strategy (by passing use_cpr=true) or the
   simple strategy (the default).
Note that as of this commit, the CPR preconditioning still has
not been implemented properly, and behaves just like the simple
strategy.
This commit is contained in:
Atgeirr Flø Rasmussen 2014-05-12 11:10:59 +02:00
parent 46e4ea7736
commit 54c07f3d0c
5 changed files with 31 additions and 18 deletions

View File

@ -38,6 +38,7 @@
#include <opm/core/linalg/LinearSolverFactory.hpp> #include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp> #include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@ -141,9 +142,13 @@ try
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0; const double *grav = use_gravity ? &gravity[0] : 0;
// Linear solver. // Solver for Newton iterations.
LinearSolverFactory linsolver(param); std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
NewtonIterationBlackoilSimple fis_solver(linsolver); if (param.getDefault("use_cpr", false)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param));
} else {
fis_solver.reset(new NewtonIterationBlackoilSimple(param));
}
// Write parameters used for later reference. // Write parameters used for later reference.
bool output = param.getDefault("output", true); bool output = param.getDefault("output", true);
@ -212,7 +217,7 @@ try
*new_props, *new_props,
rock_comp->isActive() ? rock_comp.get() : 0, rock_comp->isActive() ? rock_comp.get() : 0,
wells, wells,
fis_solver, *fis_solver,
grav); grav);
SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);

View File

@ -60,6 +60,7 @@
#include <opm/core/linalg/LinearSolverFactory.hpp> #include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp> #include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@ -194,9 +195,13 @@ try
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0; const double *grav = use_gravity ? &gravity[0] : 0;
// Linear solver. // Solver for Newton iterations.
LinearSolverFactory linsolver(param); std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
NewtonIterationBlackoilSimple fis_solver(linsolver); if (param.getDefault("use_cpr", false)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param));
} else {
fis_solver.reset(new NewtonIterationBlackoilSimple(param));
}
// Write parameters used for later reference. // Write parameters used for later reference.
bool output = param.getDefault("output", true); bool output = param.getDefault("output", true);
@ -271,7 +276,7 @@ try
*new_props, *new_props,
rock_comp->isActive() ? rock_comp.get() : 0, rock_comp->isActive() ? rock_comp.get() : 0,
wells, wells,
fis_solver, *fis_solver,
grav); grav);
SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);

View File

@ -102,8 +102,7 @@ try
double grav[] = { 0.0, 0.0 }; double grav[] = { 0.0, 0.0 };
Opm::DerivedGeology geo(*g, props, grav); Opm::DerivedGeology geo(*g, props, grav);
Opm::LinearSolverFactory linsolver(param); Opm::NewtonIterationBlackoilSimple fis_solver(param);
Opm::NewtonIterationBlackoilSimple fis_solver(linsolver);
Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(*g, props, geo, 0, *wells, fis_solver); Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(*g, props, geo, 0, *wells, fis_solver);

View File

@ -22,15 +22,16 @@
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp> #include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp> #include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
namespace Opm namespace Opm
{ {
/// Construct a system solver. /// Construct a system solver.
/// \param[in] linsolver linear solver to use /// \param[in] linsolver linear solver to use
NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver) NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param)
: linsolver_(linsolver)
{ {
linsolver_.reset(new LinearSolverFactory(param));
} }
/// Solve the linear system Ax = b, with A being the /// Solve the linear system Ax = b, with A being the
@ -54,9 +55,9 @@ namespace Opm
SolutionVector dx(SolutionVector::Zero(total_residual.size())); SolutionVector dx(SolutionVector::Zero(total_residual.size()));
Opm::LinearSolverInterface::LinearSolverReport rep Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(matr.rows(), matr.nonZeros(), = linsolver_->solve(matr.rows(), matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
total_residual.value().data(), dx.data()); total_residual.value().data(), dx.data());
if (!rep.converged) { if (!rep.converged) {
OPM_THROW(std::runtime_error, OPM_THROW(std::runtime_error,
"FullyImplicitBlackoilSolver::solveJacobianSystem(): " "FullyImplicitBlackoilSolver::solveJacobianSystem(): "

View File

@ -22,7 +22,9 @@
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp> #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/linalg/LinearSolverInterface.hpp> #include <opm/core/linalg/LinearSolverInterface.hpp>
#include <memory>
namespace Opm namespace Opm
{ {
@ -35,8 +37,9 @@ namespace Opm
{ {
public: public:
/// Construct a system solver. /// Construct a system solver.
/// \param[in] linsolver linear solver to use /// \param[in] param parameters controlling the behaviour and
NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver); /// choice of linear solver.
NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param);
/// Solve the system of linear equations Ax = b, with A being the /// Solve the system of linear equations Ax = b, with A being the
/// combined derivative matrix of the residual and b /// combined derivative matrix of the residual and b
@ -46,7 +49,7 @@ namespace Opm
virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const; virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
private: private:
const LinearSolverInterface& linsolver_; std::unique_ptr<LinearSolverInterface> linsolver_;
}; };
} // namespace Opm } // namespace Opm