OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-28 02:00:59 -06:00
Code Issues Packages Projects Releases Wiki Activity

SimulatorBase: only care about the solver in the run() method

Browse Source 
this is necessary because some older simulations only provide the
full-fledged solver class but no physical model.

(also, this allows to use something else than the standard newton
solver.)
This commit is contained in:
Andreas Lauser 2015-05-28 12:24:39 +02:00
parent a154c8394d
commit 662cd9791e
5 changed files with 51 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
opm/autodiff/NewtonSolver.hpp
View File

@ -66,7 +66,7 @@ namespace Opm {
/// \param[in] param parameters controlling nonlinear Newton process
/// \param[in, out] model physical simulation model
explicit NewtonSolver(const SolverParameters& param,
PhysicalModel& model);
std::shared_ptr<PhysicalModel> model);
/// Take a single forward step, after which the states will be modified
/// according to the physical model.
@ -94,7 +94,7 @@ namespace Opm {
private:
// --------- Data members ---------
SolverParameters param_;
PhysicalModel& model_;
std::shared_ptr<PhysicalModel> model_;
unsigned int newtonIterations_;
unsigned int linearIterations_;
unsigned int newtonIterationsLast_;

32
opm/autodiff/NewtonSolver_impl.hpp
View File

@ -29,7 +29,7 @@ namespace Opm
{
template <class PhysicalModel>
NewtonSolver<PhysicalModel>::NewtonSolver(const SolverParameters& param,
PhysicalModel& model)
std::shared_ptr<PhysicalModel> model)
: param_(param),
model_(model),
newtonIterations_(0),
@ -58,21 +58,21 @@ namespace Opm
WellState& well_state)
{
// Do model-specific once-per-step calculations.
model_.prepareStep(dt, reservoir_state, well_state);
model_->prepareStep(dt, reservoir_state, well_state);
// For each iteration we store in a vector the norms of the residual of
// the mass balance for each active phase, the well flux and the well equations.
std::vector<std::vector<double>> residual_norms_history;
// Assemble residual and Jacobian, store residual norms.
model_.assemble(reservoir_state, well_state, true);
residual_norms_history.push_back(model_.computeResidualNorms());
model_->assemble(reservoir_state, well_state, true);
residual_norms_history.push_back(model_->computeResidualNorms());
// Set up for main Newton loop.
double omega = 1.0;
int iteration = 0;
bool converged = model_.getConvergence(dt, iteration);
const int sizeNonLinear = model_.sizeNonLinear();
bool converged = model_->getConvergence(dt, iteration);
const int sizeNonLinear = model_->sizeNonLinear();
V dxOld = V::Zero(sizeNonLinear);
bool isOscillate = false;
bool isStagnate = false;
@ -82,17 +82,17 @@ namespace Opm
// ---------- Main Newton loop ----------
while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
// Compute the Newton update to the primary variables.
V dx = model_.solveJacobianSystem();
V dx = model_->solveJacobianSystem();
// Store number of linear iterations used.
linearIterations += model_.linearIterationsLastSolve();
linearIterations += model_->linearIterationsLastSolve();
// Stabilize the Newton update.
detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate);
if (isOscillate) {
omega -= relaxIncrement();
omega = std::max(omega, relaxMax());
if (model_.terminalOutputEnabled()) {
if (model_->terminalOutputEnabled()) {
std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl;
}
}
@ -100,20 +100,20 @@ namespace Opm
// Apply the update, the model may apply model-dependent
// limitations and chopping of the update.
model_.updateState(dx, reservoir_state, well_state);
model_->updateState(dx, reservoir_state, well_state);
// Assemble residual and Jacobian, store residual norms.
model_.assemble(reservoir_state, well_state, false);
residual_norms_history.push_back(model_.computeResidualNorms());
model_->assemble(reservoir_state, well_state, false);
residual_norms_history.push_back(model_->computeResidualNorms());
// increase iteration counter
++iteration;
converged = model_.getConvergence(dt, iteration);
converged = model_->getConvergence(dt, iteration);
}
if (!converged) {
if (model_.terminalOutputEnabled()) {
if (model_->terminalOutputEnabled()) {
std::cerr << "WARNING: Failed to compute converged solution in " << iteration << " iterations." << std::endl;
}
return -1; // -1 indicates that the solver has to be restarted
@ -125,7 +125,7 @@ namespace Opm
newtonIterationsLast_ = iteration;
// Do model-specific post-step actions.
model_.afterStep(dt, reservoir_state, well_state);
model_->afterStep(dt, reservoir_state, well_state);
return linearIterations;
}
@ -197,7 +197,7 @@ namespace Opm
const std::vector<double>& F0 = residual_history[it];
const std::vector<double>& F1 = residual_history[it - 1];
const std::vector<double>& F2 = residual_history[it - 2];
for (int p= 0; p < model_.numPhases(); ++p){
for (int p= 0; p < model_->numPhases(); ++p){
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);

38
opm/autodiff/SimulatorBase.hpp
View File

@ -88,7 +88,7 @@ namespace Opm
typedef typename Traits::WellState WellState;
typedef typename Traits::OutputWriter OutputWriter;
typedef typename Traits::Grid Grid;
typedef typename Traits::Model Model;
typedef typename Traits::Solver Solver;
/// Initialise from parameters and objects to observe.
/// \param[in] param parameters, this class accepts the following:
@ -150,20 +150,32 @@ namespace Opm
const Wells* wells)
{};
std::shared_ptr<Model> createModel(const typename Model::ModelParameters &modelParams,
const Wells* wells)
std::shared_ptr<Solver> createSolver(const Wells* wells)
{
return std::make_shared<Model>(modelParams,
grid_,
props_,
geo_,
rock_comp_props_,
wells,
solver_,
has_disgas_,
has_vapoil_,
terminal_output_);
typedef typename Traits::Model Model;
typedef typename Model::ModelParameters ModelParams;
ModelParams modelParams( param_ );
typedef NewtonSolver<Model> Solver;
auto model = std::make_shared<Model>(modelParams,
grid_,
props_,
geo_,
rock_comp_props_,
wells,
solver_,
has_disgas_,
has_vapoil_,
terminal_output_);
if (!threshold_pressures_by_face_.empty()) {
model->setThresholdPressures(threshold_pressures_by_face_);
}
typedef typename Solver::SolverParameters SolverParams;
SolverParams solverParams( param_ );
return std::make_shared<Solver>(solverParams, model);
}
void

20
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -83,12 +83,6 @@ namespace Opm
std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
typedef typename Model::ModelParameters ModelParams;
ModelParams modelParams( param_ );
typedef NewtonSolver<Model> Solver;
typedef typename Solver::SolverParameters SolverParams;
SolverParams solverParams( param_ );
// adaptive time stepping
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
if( param_.getDefault("timestep.adaptive", true ) )
@ -150,11 +144,7 @@ namespace Opm
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
auto model = asImp_().createModel(modelParams, wells);
if (!threshold_pressures_by_face_.empty()) {
model->setThresholdPressures(threshold_pressures_by_face_);
}
Solver solver(solverParams, *model);
auto solver = asImp_().createSolver(wells);
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are to large for the solver to converge
@ -162,19 +152,19 @@ namespace Opm
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
if( adaptiveTimeStepping ) {
adaptiveTimeStepping->step( timer, solver, state, well_state, output_writer_ );
adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ );
}
else {
// solve for complete report step
solver.step(timer.currentStepLength(), state, well_state);
solver->step(timer.currentStepLength(), state, well_state);
}
// take time that was used to solve system for this reportStep
solver_timer.stop();
// accumulate the number of Newton and Linear Iterations
totalNewtonIterations += solver.newtonIterations();
totalLinearIterations += solver.linearIterations();
totalNewtonIterations += solver->newtonIterations();
totalLinearIterations += solver->linearIterations();
// Report timing.
const double st = solver_timer.secsSinceStart();

3
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -22,6 +22,8 @@
#include "SimulatorBase.hpp"
#include "NewtonSolver.hpp"
namespace Opm {
template<class GridT>
class SimulatorFullyImplicitBlackoil;
@ -34,6 +36,7 @@ struct SimulatorTraits<SimulatorFullyImplicitBlackoil<GridT> >
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef BlackoilModel<Grid> Model;
typedef NewtonSolver<Model> Solver;
};
/// a simulator for the blackoil model