OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-25 10:40:21 -06:00
Code Issues Packages Projects Releases Wiki Activity

flow: avoid wrangling around with dummy state objects

Browse Source 
these objects are only used by flow_legacy, so not having to deal with
them anymore lets non-legacy flow avoid to jump through a lot of hoops
for the sake of having a common API.

this required a fork of the NonlinearSolver and AdaptiveTimeStepping
classes. this is not a problem because the original classes would get
pruned to look like the new ones once flow_legacy gets moved out of
the opm-simulators module.
This commit is contained in:
Andreas Lauser 2018-06-06 10:59:41 +02:00
parent 2672f38375
commit dfbc24b35f
13 changed files with 1091 additions and 202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -180,11 +180,7 @@ namespace Opm {
/// Called once before each time step.
/// \param[in] timer simulation timer
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void prepareStep(const SimulatorTimerInterface& timer,
const ReservoirState& /*reservoir_state*/,
const WellState& /* well_state */)
void prepareStep(const SimulatorTimerInterface& timer)
{
// update the solution variables in ebos
@ -228,9 +224,7 @@ namespace Opm {
template <class NonlinearSolverType>
SimulatorReport nonlinearIteration(const int iteration,
const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver,
ReservoirState& /*reservoir_state*/,
WellState& /*well_state*/)
NonlinearSolverType& nonlinear_solver)
{
SimulatorReport report;
failureReport_ = SimulatorReport();
@ -341,16 +335,8 @@ namespace Opm {
/// Called once after each time step.
/// In this class, this function does nothing.
/// \param[in] timer simulation timer
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void afterStep(const SimulatorTimerInterface& timer,
const ReservoirState& reservoir_state,
WellState& well_state)
void afterStep(const SimulatorTimerInterface& OPM_UNUSED timer)
{
DUNE_UNUSED_PARAMETER(timer);
DUNE_UNUSED_PARAMETER(reservoir_state);
DUNE_UNUSED_PARAMETER(well_state);
wellModel().timeStepSucceeded();
aquiferModel().timeStepSucceeded(timer);
ebosSimulator_.problem().endTimeStep();
@ -411,8 +397,7 @@ namespace Opm {
}
// compute the "relative" change of the solution between time steps
template <class Dummy>
double relativeChange(const Dummy&, const Dummy&) const
double relativeChange() const
{
Scalar resultDelta = 0.0;
Scalar resultDenom = 0.0;
@ -484,7 +469,7 @@ namespace Opm {
resultDenom = gridView.comm().sum(resultDenom);
if (resultDenom > 0.0)
return resultDelta/resultDenom;
return resultDelta/resultDenom;
return 0.0;
}

62
opm/autodiff/BlackoilOutputEbos.hpp
View File

@ -96,13 +96,11 @@ namespace Opm
* requested output cell properties specified by the RPTRST keyword
* and write these to file.
*/
template<class SimulationDataContainer, class Model>
template<class Model>
void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulationDataContainer& /*reservoirStateDummy*/,
const Opm::WellStateFullyImplicitBlackoil& /*wellStateDummy*/,
const Model& physicalModel,
const bool substep = false,
const double nextstep = -1.0,
const bool isSubstep = false,
const double nextStepSize = -1.0,
const SimulatorReport& simulatorReport = SimulatorReport())
{
if( output_ )
@ -114,59 +112,7 @@ namespace Opm
// The writeOutput expects a local data::solution vector and a local data::well vector.
auto localWellData = localWellState.report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid()) );
ebosSimulator_.problem().writeOutput(localWellData, timer.simulationTimeElapsed(), substep, totalSolverTime, nextstep);
}
}
template <class SimulationDataContainer, class WellState>
void initFromRestartFile(const PhaseUsage& /*phaseUsage*/,
const Grid& /*grid */,
SimulationDataContainer& simulatorstate,
WellState& wellstate,
ExtraData& extra) {
std::vector<RestartKey> extra_keys = {
{"OPMEXTRA" , Opm::UnitSystem::measure::identity, false}
};
// gives a dummy dynamic_list_econ_limited
DynamicListEconLimited dummy_list_econ_limited;
const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
WellsManager wellsmanager(eclState(),
schedule(),
// The restart step value is used to identify wells present at the given
// time step. Wells that are added at the same time step as RESTART is initiated
// will not be present in a restart file. Use the previous time step to retrieve
// wells that have information written to the restart file.
std::max(eclState().getInitConfig().getRestartStep() - 1, 0),
Opm::UgGridHelpers::numCells(grid()),
Opm::UgGridHelpers::globalCell(grid()),
Opm::UgGridHelpers::cartDims(grid()),
Opm::UgGridHelpers::dimensions(grid()),
Opm::UgGridHelpers::cell2Faces(grid()),
Opm::UgGridHelpers::beginFaceCentroids(grid()),
dummy_list_econ_limited,
grid().comm().size() > 1,
defunct_well_names);
const Wells* wells = wellsmanager.c_wells();
std::vector<RestartKey> solution_keys = {};
auto restart_values = ebosSimulator_.problem().eclIO().loadRestart(solution_keys, extra_keys);
const int nw = wells->number_of_wells;
if (nw > 0) {
wellstate.resize(wells, simulatorstate, phaseUsage_ ); //Resize for restart step
wellsToState( restart_values.wells, phaseUsage_, wellstate );
}
if (restart_values.hasExtra("OPMEXTRA")) {
std::vector<double> opmextra = restart_values.getExtra("OPMEXTRA");
assert(opmextra.size() == 1);
extra.suggested_step = opmextra[0];
} else {
OpmLog::warning("Restart data is missing OPMEXTRA field, restart run may deviate from original run.");
extra.suggested_step = -1.0;
ebosSimulator_.problem().writeOutput(localWellData, timer.simulationTimeElapsed(), isSubstep, totalSolverTime, nextStepSize);
}
}

66
opm/autodiff/BlackoilWellModel.hpp
View File

@ -103,6 +103,40 @@ namespace Opm {
const ModelParameters& param,
const bool terminal_output);
void initFromRestartFile(const RestartValue& restartValues)
{
// gives a dummy dynamic_list_econ_limited
DynamicListEconLimited dummyListEconLimited;
const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames();
WellsManager wellsmanager(eclState(),
schedule(),
// The restart step value is used to identify wells present at the given
// time step. Wells that are added at the same time step as RESTART is initiated
// will not be present in a restart file. Use the previous time step to retrieve
// wells that have information written to the restart file.
std::max(eclState().getInitConfig().getRestartStep() - 1, 0),
Opm::UgGridHelpers::numCells(grid()),
Opm::UgGridHelpers::globalCell(grid()),
Opm::UgGridHelpers::cartDims(grid()),
Opm::UgGridHelpers::dimensions(grid()),
Opm::UgGridHelpers::cell2Faces(grid()),
Opm::UgGridHelpers::beginFaceCentroids(grid()),
dummyListEconLimited,
grid().comm().size() > 1,
defunctWellNames);
const Wells* wells = wellsmanager.c_wells();
const int nw = wells->number_of_wells;
if (nw > 0) {
auto phaseUsage = phaseUsageFromDeck(eclState());
size_t numCells = Opm::UgGridHelpers::numCells(grid());
well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step
wellsToState(restartValues.wells, phaseUsage, well_state_);
previous_well_state_ = well_state_;
}
}
// compute the well fluxes and assemble them in to the reservoir equations as source terms
// and in the well equations.
void assemble(const int iterationIdx,
@ -161,6 +195,32 @@ namespace Opm {
}
protected:
void extractLegacyPressure_(std::vector<double>& cellPressure) const
{
size_t nc = number_of_cells_;
std::vector<double> cellPressures(nc, 0.0);
ElementContext elemCtx(ebosSimulator_);
const auto& gridView = ebosSimulator_.vanguard().gridView();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
cellPressures[cellIdx] = p;
}
}
Simulator& ebosSimulator_;
std::unique_ptr<WellsManager> wells_manager_;
@ -204,6 +264,12 @@ namespace Opm {
const Wells* wells() const { return wells_manager_->c_wells(); }
const Grid& grid() const
{ return ebosSimulator_.vanguard().grid(); }
const EclipseState& eclState() const
{ return ebosSimulator_.vanguard().eclState(); }
const Schedule& schedule() const
{ return ebosSimulator_.vanguard().schedule(); }

2
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@ -88,7 +88,7 @@ namespace Opm {
const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
cellPressures[cellIdx] = p;
}
well_state_.init(wells(), cellPressures, previous_well_state_, phase_usage_);
well_state_.init(wells(), cellPressures, &previous_well_state_, phase_usage_);
// handling MS well related
if (param_.use_multisegment_well_) { // if we use MultisegmentWell model

348
opm/autodiff/NonlinearSolverEbos.hpp Normal file
View File

@ -0,0 +1,348 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
Copyright 2015 Statoil ASA.
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/>.
*/
#ifndef OPM_NON_LINEAR_SOLVER_EBOS_HPP
#define OPM_NON_LINEAR_SOLVER_EBOS_HPP
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <memory>
namespace Opm {
/// A nonlinear solver class suitable for general fully-implicit models,
/// as well as pressure, transport and sequential models.
template <class PhysicalModel>
class NonlinearSolverEbos
{
public:
// Available relaxation scheme types.
enum RelaxType {
Dampen,
SOR
};
// Solver parameters controlling nonlinear process.
struct SolverParametersEbos
{
RelaxType relaxType_;
double relaxMax_;
double relaxIncrement_;
double relaxRelTol_;
int maxIter_; // max nonlinear iterations
int minIter_; // min nonlinear iterations
explicit SolverParametersEbos(const ParameterGroup& param)
{
// set default values
reset();
// overload with given parameters
relaxMax_ = param.getDefault("relax_max", relaxMax_);
maxIter_ = param.getDefault("max_iter", maxIter_);
minIter_ = param.getDefault("min_iter", minIter_);
std::string relaxationType = param.getDefault("relax_type", std::string("dampen"));
if (relaxationType == "dampen") {
relaxType_ = Dampen;
} else if (relaxationType == "sor") {
relaxType_ = SOR;
} else {
OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxationType);
}
}
SolverParametersEbos()
{ reset(); }
void reset()
{
// default values for the solver parameters
relaxType_ = Dampen;
relaxMax_ = 0.5;
relaxIncrement_ = 0.1;
relaxRelTol_ = 0.2;
maxIter_ = 10;
minIter_ = 1;
}
};
// Forwarding types from PhysicalModel.
typedef typename PhysicalModel::ReservoirState ReservoirState;
typedef typename PhysicalModel::WellState WellState;
// --------- Public methods ---------
/// Construct solver for a given model.
///
/// The model is a std::unique_ptr because the object to which model points to is
/// not allowed to be deleted as long as the NonlinearSolver object exists.
///
/// \param[in] param parameters controlling nonlinear process
/// \param[in, out] model physical simulation model.
NonlinearSolverEbos(const SolverParametersEbos& param,
std::unique_ptr<PhysicalModel> model)
: param_(param)
, model_(std::move(model))
, linearizations_(0)
, nonlinearIterations_(0)
, linearIterations_(0)
, wellIterations_(0)
, nonlinearIterationsLast_(0)
, linearIterationsLast_(0)
, wellIterationsLast_(0)
{
if (!model_) {
OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver.");
}
}
SimulatorReport step(const SimulatorTimerInterface& timer)
{
SimulatorReport iterReport;
SimulatorReport report;
failureReport_ = SimulatorReport();
// Do model-specific once-per-step calculations.
model_->prepareStep(timer);
int iteration = 0;
// Let the model do one nonlinear iteration.
// Set up for main solver loop.
bool converged = false;
// ---------- Main nonlinear solver loop ----------
do {
try {
// Do the nonlinear step. If we are in a converged state, the
// model will usually do an early return without an expensive
// solve, unless the minIter() count has not been reached yet.
iterReport = model_->nonlinearIteration(iteration, timer, *this);
report += iterReport;
report.converged = iterReport.converged;
converged = report.converged;
iteration += 1;
}
catch (...) {
// if an iteration fails during a time step, all previous iterations
// count as a failure as well
failureReport_ += report;
failureReport_ += model_->failureReport();
throw;
}
}
while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
if (!converged) {
failureReport_ += report;
std::string msg = "Solver convergence failure - Failed to complete a time step within " + std::to_string(maxIter()) + " iterations.";
OPM_THROW_NOLOG(Opm::TooManyIterations, msg);
}
// Do model-specific post-step actions.
model_->afterStep(timer);
report.converged = true;
return report;
}
/// return the statistics if the step() method failed
const SimulatorReport& failureReport() const
{ return failureReport_; }
/// Number of linearizations used in all calls to step().
int linearizations() const
{ return linearizations_; }
/// Number of full nonlinear solver iterations used in all calls to step().
int nonlinearIterations() const
{ return nonlinearIterations_; }
/// Number of linear solver iterations used in all calls to step().
int linearIterations() const
{ return linearIterations_; }
/// Number of well iterations used in all calls to step().
int wellIterations() const
{ return wellIterations_; }
/// Number of nonlinear solver iterations used in the last call to step().
int nonlinearIterationsLastStep() const
{ return nonlinearIterationsLast_; }
/// Number of linear solver iterations used in the last call to step().
int linearIterationsLastStep() const
{ return linearIterationsLast_; }
/// Number of well iterations used in all calls to step().
int wellIterationsLastStep() const
{ return wellIterationsLast_; }
/// Compute fluid in place.
/// \param[in] ReservoirState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<std::vector<double> >
computeFluidInPlace(const ReservoirState& x, const std::vector<int>& fipnum) const
{ return model_->computeFluidInPlace(x, fipnum); }
std::vector<std::vector<double> >
computeFluidInPlace(const std::vector<int>& fipnum) const
{ return model_->computeFluidInPlace(fipnum); }
/// Reference to physical model.
const PhysicalModel& model() const
{ return *model_; }
/// Mutable reference to physical model.
PhysicalModel& model()
{ return *model_; }
/// Detect oscillation or stagnation in a given residual history.
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const int it, bool& oscillate, bool& stagnate) const
{
// The detection of oscillation in two primary variable results in the report of the detection
// of oscillation for the solver.
// Only the saturations are used for oscillation detection for the black oil model.
// Stagnate is not used for any treatment here.
if ( it < 2 ) {
oscillate = false;
stagnate = false;
return;
}
stagnate = true;
int oscillatePhase = 0;
const std::vector<double>& F0 = residualHistory[it];
const std::vector<double>& F1 = residualHistory[it - 1];
const std::vector<double>& F2 = residualHistory[it - 2];
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]);
oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
// Process is 'stagnate' unless at least one phase
// exhibits significant residual change.
stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
}
oscillate = (oscillatePhase > 1);
}
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
/// Implemention for Dune block vectors.
template <class BVector>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
{
// The dxOld is updated with dx.
// If omega is equal to 1., no relaxtion will be appiled.
BVector tempDxOld = dxOld;
dxOld = dx;
switch (relaxType()) {
case Dampen: {
if (omega == 1.) {
return;
}
auto i = dx.size();
for (i = 0; i < dx.size(); ++i) {
dx[i] *= omega;
}
return;
}
case SOR: {
if (omega == 1.) {
return;
}
auto i = dx.size();
for (i = 0; i < dx.size(); ++i) {
dx[i] *= omega;
tempDxOld[i] *= (1.-omega);
dx[i] += tempDxOld[i];
}
return;
}
default:
OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
}
return;
}
/// The greatest relaxation factor (i.e. smallest factor) allowed.
double relaxMax() const
{ return param_.relaxMax_; }
/// The step-change size for the relaxation factor.
double relaxIncrement() const
{ return param_.relaxIncrement_; }
/// The relaxation type (Dampen or SOR).
enum RelaxType relaxType() const
{ return param_.relaxType_; }
/// The relaxation relative tolerance.
double relaxRelTol() const
{ return param_.relaxRelTol_; }
/// The maximum number of nonlinear iterations allowed.
int maxIter() const
{ return param_.maxIter_; }
/// The minimum number of nonlinear iterations allowed.
int minIter() const
{ return param_.minIter_; }
/// Set parameters to override those given at construction time.
void setParameters(const SolverParametersEbos& param)
{ param_ = param; }
private:
// --------- Data members ---------
SimulatorReport failureReport_;
SolverParametersEbos param_;
std::unique_ptr<PhysicalModel> model_;
int linearizations_;
int nonlinearIterations_;
int linearIterations_;
int wellIterations_;
int nonlinearIterationsLast_;
int linearIterationsLast_;
int wellIterationsLast_;
};
} // namespace Opm
#endif // OPM_NON_LINEAR_SOLVER_EBOS_HPP

5
opm/autodiff/ParallelDebugOutput.hpp
View File

@ -596,8 +596,9 @@ namespace Opm
};
// gather solution to rank 0 for EclipseWriter
template <class WellState>
bool collectToIORank( const SimulationDataContainer& /*localReservoirState*/,
const WellStateFullyImplicitBlackoil& localWellState,
const WellState& localWellState,
const data::Solution& localCellData,
const int wellStateStepNumber )
{
@ -627,7 +628,7 @@ namespace Opm
std::unordered_set<std::string>());
const Wells* wells = wells_manager.c_wells();
globalWellState_.init(wells, *globalReservoirState_, globalWellState_, phaseUsage_ );
globalWellState_.initLegacy(wells, *globalReservoirState_, globalWellState_, phaseUsage_ );
globalCellData_->clear();
}

2
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -174,7 +174,7 @@ namespace Opm
defunct_well_names_);
const Wells* wells = wells_manager.c_wells();
WellState well_state;
well_state.init(wells, state, prev_well_state, props_.phaseUsage());
well_state.initLegacy(wells, state, prev_well_state, props_.phaseUsage());
// give the polymer and surfactant simulators the chance to do their stuff
asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);

195
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -24,14 +24,14 @@
#include <opm/autodiff/BlackoilOutputEbos.hpp>
#include <opm/autodiff/IterationReport.hpp>
#include <opm/autodiff/NonlinearSolver.hpp>
#include <opm/autodiff/NonlinearSolverEbos.hpp>
#include <opm/autodiff/BlackoilModelEbos.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/BlackoilAquiferModel.hpp>
#include <opm/autodiff/moduleVersion.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/common/Exceptions.hpp>
@ -64,7 +64,7 @@ public:
typedef BlackoilOutputEbos<TypeTag> OutputWriter;
typedef BlackoilModelEbos<TypeTag> Model;
typedef BlackoilModelParameters ModelParameters;
typedef NonlinearSolver<Model> Solver;
typedef NonlinearSolverEbos<Model> Solver;
typedef BlackoilWellModel<TypeTag> WellModel;
typedef BlackoilAquiferModel<TypeTag> AquiferModel;
@ -95,29 +95,26 @@ public:
SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
const ParameterGroup& param,
NewtonIterationBlackoilInterface& linsolver,
const bool has_disgas,
const bool has_vapoil,
OutputWriter& output_writer)
: ebosSimulator_(ebosSimulator),
param_(param),
model_param_(param),
solver_param_(param),
solver_(linsolver),
phaseUsage_(phaseUsageFromDeck(eclState())),
has_disgas_(has_disgas),
has_vapoil_(has_vapoil),
terminal_output_(param.getDefault("output_terminal", true)),
output_writer_(output_writer),
is_parallel_run_( false )
const bool hasDisgas,
const bool hasVapoil,
OutputWriter& outputWriter)
: ebosSimulator_(ebosSimulator)
, param_(param)
, modelParam_(param)
, solverParam_(param)
, solver_(linsolver)
, phaseUsage_(phaseUsageFromDeck(eclState()))
, hasDisgas_(hasDisgas)
, hasVapoil_(hasVapoil)
, terminalOutput_(param.getDefault("output_terminal", true))
, outputWriter_(outputWriter)
{
#if HAVE_MPI
if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
{
if (solver_.parallelInformation().type() == typeid(ParallelISTLInformation)) {
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
// Only rank 0 does print to std::cout
terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 );
is_parallel_run_ = ( info.communicator().size() > 1 );
terminalOutput_ = terminalOutput_ && (info.communicator().rank() == 0);
}
#endif
}
@ -130,43 +127,51 @@ public:
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer)
{
ReservoirState dummy_state(0,0,0);
WellState prev_well_state;
ExtraData extra;
failureReport_ = SimulatorReport();
if (output_writer_.isRestart()) {
// This is a restart, populate WellState
ReservoirState stateInit(Opm::UgGridHelpers::numCells(grid()),
Opm::UgGridHelpers::numFaces(grid()),
phaseUsage_.num_phases);
output_writer_.initFromRestartFile(phaseUsage_, grid(), stateInit, prev_well_state, extra);
// handle restarts
std::unique_ptr<RestartValue> restartValues;
if (outputWriter_.isRestart()) {
std::vector<RestartKey> extraKeys = {
{"OPMEXTRA" , Opm::UnitSystem::measure::identity, false}
};
std::vector<RestartKey> solutionKeys = {};
restartValues.reset(new RestartValue(ebosSimulator_.problem().eclIO().loadRestart(solutionKeys, extraKeys)));
}
// Create timers and file for writing timing info.
Opm::time::StopWatch solver_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
Opm::time::StopWatch solverTimer;
Opm::time::StopWatch totalTimer;
totalTimer.start();
// adaptive time stepping
const auto& events = schedule().getEvents();
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
std::unique_ptr< AdaptiveTimeSteppingEbos > adaptiveTimeStepping;
const bool useTUNING = param_.getDefault("use_TUNING", false);
if( param_.getDefault("timestep.adaptive", true ) )
{
if (param_.getDefault("timestep.adaptive", true)) {
if (useTUNING) {
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule().getTuning(), timer.currentStepNum(), param_, terminal_output_ ) );
} else {
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
adaptiveTimeStepping.reset(new AdaptiveTimeSteppingEbos(schedule().getTuning(), timer.currentStepNum(), param_, terminalOutput_));
}
else {
adaptiveTimeStepping.reset(new AdaptiveTimeSteppingEbos(param_, terminalOutput_));
}
if (output_writer_.isRestart()) {
if (extra.suggested_step > 0.0) {
adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
double suggestedStepSize = -1.0;
if (outputWriter_.isRestart()) {
// This is a restart, determine the time step size from the restart data
if (restartValues->hasExtra("OPMEXTRA")) {
std::vector<double> opmextra = restartValues->getExtra("OPMEXTRA");
assert(opmextra.size() == 1);
suggestedStepSize = opmextra[0];
}
else {
OpmLog::warning("Restart data is missing OPMEXTRA field, restart run may deviate from original run.");
suggestedStepSize = -1.0;
}
if (suggestedStepSize > 0.0) {
adaptiveTimeStepping->setSuggestedNextStep(suggestedStepSize);
}
}
}
@ -174,15 +179,13 @@ public:
SimulatorReport report;
SimulatorReport stepReport;
WellModel well_model(ebosSimulator_, model_param_, terminal_output_);
if (output_writer_.isRestart()) {
well_model.setRestartWellState(prev_well_state); // Neccessary for perfect restarts
WellModel wellModel(ebosSimulator_, modelParam_, terminalOutput_);
if (outputWriter_.isRestart()) {
wellModel.initFromRestartFile(*restartValues);
}
WellState wellStateDummy; //not used. Only passed to make the old interfaces happy
if ( model_param_.matrix_add_well_contributions_ ||
model_param_.preconditioner_add_well_contributions_ )
if (modelParam_.matrix_add_well_contributions_ ||
modelParam_.preconditioner_add_well_contributions_)
{
ebosSimulator_.model().clearAuxiliaryModules();
auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(schedule(), grid());
@ -194,19 +197,18 @@ public:
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
if ( terminal_output_ )
{
if (terminalOutput_) {
std::ostringstream ss;
timer.report(ss);
OpmLog::debug(ss.str());
}
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
solverTimer.start();
well_model.beginReportStep(timer.currentStepNum());
wellModel.beginReportStep(timer.currentStepNum());
auto solver = createSolver(well_model, aquifer_model);
auto solver = createSolver(wellModel, aquifer_model);
// write the inital state at the report stage
if (timer.initialStep()) {
@ -215,22 +217,21 @@ public:
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model() );
outputWriter_.writeTimeStep(timer, solver->model());
report.output_write_time += perfTimer.stop();
}
if( terminal_output_ )
{
std::ostringstream step_msg;
if (terminalOutput_) {
std::ostringstream stepMsg;
boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
step_msg.imbue(std::locale(std::locale::classic(), facet));
step_msg << "\nReport step " << std::setw(2) <<timer.currentStepNum()
stepMsg.imbue(std::locale(std::locale::classic(), facet));
stepMsg << "\nReport step " << std::setw(2) <<timer.currentStepNum()
<< "/" << timer.numSteps()
<< " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
<< "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
<< ", date = " << timer.currentDateTime();
OpmLog::info(step_msg.str());
OpmLog::info(stepMsg.str());
}
solver->model().beginReportStep();
@ -240,9 +241,9 @@ public:
//
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
if( adaptiveTimeStepping ) {
if (adaptiveTimeStepping) {
if (useTUNING) {
if(events.hasEvent(ScheduleEvents::TUNING_CHANGE,timer.currentStepNum())) {
if (events.hasEvent(ScheduleEvents::TUNING_CHANGE,timer.currentStepNum())) {
adaptiveTimeStepping->updateTUNING(schedule().getTuning(), timer.currentStepNum());
}
}
@ -251,18 +252,17 @@ public:
events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
stepReport = adaptiveTimeStepping->step( timer, *solver, dummy_state, wellStateDummy, event, output_writer_, nullptr );
stepReport = adaptiveTimeStepping->step(timer, *solver, event, outputWriter_, nullptr);
report += stepReport;
failureReport_ += adaptiveTimeStepping->failureReport();
}
else {
// solve for complete report step
stepReport = solver->step(timer, dummy_state, wellStateDummy);
stepReport = solver->step(timer);
report += stepReport;
failureReport_ += solver->failureReport();
if( terminal_output_ )
{
if (terminalOutput_) {
std::ostringstream ss;
stepReport.reportStep(ss);
OpmLog::info(ss.str());
@ -270,20 +270,19 @@ public:
}
solver->model().endReportStep();
well_model.endReportStep();
wellModel.endReportStep();
// take time that was used to solve system for this reportStep
solver_timer.stop();
solverTimer.stop();
// update timing.
report.solver_time += solver_timer.secsSinceStart();
report.solver_time += solverTimer.secsSinceStart();
// Increment timer, remember well state.
++timer;
if (terminal_output_ )
{
if (terminalOutput_) {
if (!timer.initialStep()) {
const std::string version = moduleVersionName();
outputTimestampFIP(timer, version);
@ -295,13 +294,12 @@ public:
perfTimer.start();
const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
output_writer_.writeTimeStep( timer, dummy_state, well_model.wellState(), solver->model(), false, nextstep, report);
outputWriter_.writeTimeStep(timer, solver->model(), false, nextstep, report);
report.output_write_time += perfTimer.stop();
if (terminal_output_ )
{
if (terminalOutput_) {
std::string msg =
"Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; "
"Time step took " + std::to_string(solverTimer.secsSinceStart()) + " seconds; "
"total solver time " + std::to_string(report.solver_time) + " seconds.";
OpmLog::debug(msg);
}
@ -309,8 +307,8 @@ public:
}
// Stop timer and create timing report
total_timer.stop();
report.total_time = total_timer.secsSinceStart();
totalTimer.stop();
report.total_time = totalTimer.secsSinceStart();
report.converged = true;
return report;
@ -318,23 +316,24 @@ public:
/** \brief Returns the simulator report for the failed substeps of the simulation.
*/
const SimulatorReport& failureReport() const { return failureReport_; };
const SimulatorReport& failureReport() const
{ return failureReport_; };
const Grid& grid() const
{ return ebosSimulator_.vanguard().grid(); }
protected:
std::unique_ptr<Solver> createSolver(WellModel& well_model, AquiferModel& aquifer_model)
std::unique_ptr<Solver> createSolver(WellModel& wellModel, AquiferModel& aquifer_model)
{
auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
model_param_,
well_model,
modelParam_,
wellModel,
aquifer_model,
solver_,
terminal_output_));
terminalOutput_));
return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
return std::unique_ptr<Solver>(new Solver(solverParam_, std::move(model)));
}
void outputTimestampFIP(const SimulatorTimer& timer, const std::string version)
@ -362,27 +361,23 @@ protected:
// Data.
Simulator& ebosSimulator_;
typedef typename Solver::SolverParameters SolverParameters;
typedef typename Solver::SolverParametersEbos SolverParametersEbos;
SimulatorReport failureReport_;
const ParameterGroup param_;
ModelParameters model_param_;
SolverParameters solver_param_;
ModelParameters modelParam_;
SolverParametersEbos solverParam_;
// Observed objects.
NewtonIterationBlackoilInterface& solver_;
PhaseUsage phaseUsage_;
// Misc. data
const bool has_disgas_;
const bool has_vapoil_;
bool terminal_output_;
// output_writer
OutputWriter& output_writer_;
// Whether this a parallel simulation or not
bool is_parallel_run_;
const bool hasDisgas_;
const bool hasVapoil_;
bool terminalOutput_;
// outputWriter
OutputWriter& outputWriter_;
};
} // namespace Opm

3
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@ -459,7 +459,8 @@ namespace Opm
std::unordered_set<std::string>());
const Wells* wells = wellsmanager.c_wells();
wellstate.resize(wells, simulatorstate, phaseUsage ); //Resize for restart step
size_t numCells = Opm::UgGridHelpers::numCells(grid);
wellstate.resize(wells, numCells, phaseUsage ); //Resize for restart step
auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
solutionToSim( restart_values, phaseUsage, simulatorstate );

6
opm/autodiff/StandardWell_impl.hpp
View File

@ -658,9 +658,8 @@ namespace Opm
{
const double target_rate = well_controls_get_current_target(well_controls_); // surface rate target
if (well_type_ == INJECTOR) {
const WellInjectionProperties& injection = well_ecl_->getInjectionProperties(current_step_);
// only handles single phase injection now
assert(injection.injectorType != WellInjector::MULTI);
assert(well_ecl_->getInjectionProperties(current_step_).injectorType != WellInjector::MULTI);
control_eq = getWQTotal() - target_rate;
} else if (well_type_ == PRODUCER) {
if (target_rate != 0.) {
@ -704,9 +703,8 @@ namespace Opm
{
const double target_rate = well_controls_get_current_target(well_controls_); // reservoir rate target
if (well_type_ == INJECTOR) {
const WellInjectionProperties& injection = well_ecl_->getInjectionProperties(current_step_);
// only handles single phase injection now
assert(injection.injectorType != WellInjector::MULTI);
assert(well_ecl_->getInjectionProperties(current_step_).injectorType != WellInjector::MULTI);
const double* distr = well_controls_get_current_distr(well_controls_);
for (int phase = 0; phase < number_of_phases_; ++phase) {
if (distr[phase] > 0.0) {

203
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@ -59,17 +59,10 @@ namespace Opm
using BaseType :: numWells;
using BaseType :: numPhases;
template <class State, class PrevWellState>
void init(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
{
init(wells, state.pressure(), prevState, pu);
}
/// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls
template <class PrevWellState>
void init(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
void init(const Wells* wells, const std::vector<double>& cellPressures, const WellStateFullyImplicitBlackoil* prevState, const PhaseUsage& pu)
{
// call init on base class
BaseType :: init(wells, cellPressures);
@ -92,10 +85,10 @@ namespace Opm
well_vaporized_oil_rates_.resize(nw, 0.0);
is_new_well_.resize(nw, true);
if ( !prevState.wellMap().empty() ) {
const auto& end = prevState.wellMap().end();
if (prevState && !prevState->wellMap().empty()) {
const auto& end = prevState->wellMap().end();
for (int w = 0; w < nw; ++w) {
const auto& it = prevState.wellMap().find( wells->name[w]);
const auto& it = prevState->wellMap().find( wells->name[w]);
if (it != end) {
is_new_well_[w] = false;
}
@ -135,6 +128,175 @@ namespace Opm
perfRateSolvent_.clear();
perfRateSolvent_.resize(nperf, 0.0);
// intialize wells that have been there before
// order may change so the mapping is based on the well name
if(prevState && !prevState->wellMap().empty()) {
typedef typename WellMapType :: const_iterator const_iterator;
const_iterator end = prevState->wellMap().end();
for (int w = 0; w < nw; ++w) {
std::string name( wells->name[ w ] );
const_iterator it = prevState->wellMap().find( name );
if( it != end )
{
const int oldIndex = (*it).second[ 0 ];
const int newIndex = w;
// bhp
bhp()[ newIndex ] = prevState->bhp()[ oldIndex ];
// thp
thp()[ newIndex ] = prevState->thp()[ oldIndex ];
// wellrates
for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
{
wellRates()[ idx ] = prevState->wellRates()[ oldidx ];
}
// perfPhaseRates
const int oldPerf_idx_beg = (*it).second[ 1 ];
const int num_perf_old_well = (*it).second[ 2 ];
const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
// copy perforation rates when the number of perforations is equal,
// otherwise initialize perfphaserates to well rates divided by the number of perforations.
if( num_perf_old_well == num_perf_this_well )
{
int old_perf_phase_idx = oldPerf_idx_beg *np;
for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np;
perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx )
{
perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ];
}
} else {
for (int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) {
for (int p = 0; p < np; ++p) {
perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well);
}
}
}
// perfPressures
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{
perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ];
}
}
// perfSolventRates
if (pu.has_solvent) {
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{
perfRateSolvent()[ perf ] = prevState->perfRateSolvent()[ oldPerf_idx ];
}
}
}
}
// If in the new step, there is no THP related target/limit anymore, its thp value should be
// set to zero.
const WellControls* ctrl = wells->ctrls[w];
const int nwc = well_controls_get_num(ctrl);
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
break;
}
}
// not finding any thp related control/limits
if (ctrl_index == nwc) {
thp()[w] = 0.;
}
}
}
{
// we need to create a trival segment related values to avoid there will be some
// multi-segment wells added later.
top_segment_index_.reserve(nw);
for (int w = 0; w < nw; ++w) {
top_segment_index_.push_back(w);
}
segpress_ = bhp();
segrates_ = wellRates();
}
}
void resize(const Wells* wells, size_t numCells, const PhaseUsage& pu)
{
std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
init(wells, tmp, nullptr, pu);
}
/// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls.
///
/// this method is only for flow_legacy!
template <class PrevWellState>
void initLegacy(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
{
// call init on base class
BaseType :: init(wells, cellPressures);
// if there are no well, do nothing in init
if (wells == 0) {
return;
}
const int nw = wells->number_of_wells;
if( nw == 0 ) return ;
// Initialize perfphaserates_, which must be done here.
const int np = wells->number_of_phases;
const int nperf = wells->well_connpos[nw];
well_reservoir_rates_.resize(nw * np, 0.0);
well_dissolved_gas_rates_.resize(nw, 0.0);
well_vaporized_oil_rates_.resize(nw, 0.0);
// Ensure that we start out with zero rates by default.
perfphaserates_.clear();
perfphaserates_.resize(nperf * np, 0.0);
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
const WellControls* ctrl = wells->ctrls[w];
if (well_controls_well_is_stopped(ctrl)) {
// Shut well: perfphaserates_ are all zero.
} else {
const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
// Open well: Initialize perfphaserates_ to well
// rates divided by the number of perforations.
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
for (int p = 0; p < np; ++p) {
perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
}
perfPress()[perf] = cellPressures[wells->well_cells[perf]];
}
}
}
// Initialize current_controls_.
// The controls set in the Wells object are treated as defaults,
// and also used for initial values.
current_controls_.resize(nw);
for (int w = 0; w < nw; ++w) {
current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
}
is_new_well_.resize(nw, true);
perfRateSolvent_.clear();
perfRateSolvent_.resize(nperf, 0.0);
// intialize wells that have been there before
// order may change so the mapping is based on the well name
if( ! prevState.wellMap().empty() )
@ -146,6 +308,9 @@ namespace Opm
const_iterator it = prevState.wellMap().find( name );
if( it != end )
{
// this is not a new added well
is_new_well_[w] = false;
const int oldIndex = (*it).second[ 0 ];
const int newIndex = w;
@ -236,10 +401,19 @@ namespace Opm
}
}
// this method is only for flow_legacy!
template <class State, class PrevWellState>
void initLegacy(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
{
initLegacy(wells, state.pressure(), prevState, pu);
}
// this method is only for flow_legacy!
template <class State>
void resize(const Wells* wells, const State& state, const PhaseUsage& pu) {
void resizeLegacy(const Wells* wells, const State& state, const PhaseUsage& pu)
{
const WellStateFullyImplicitBlackoil dummy_state{}; // Init with an empty previous state only resizes
init(wells, state, dummy_state, pu) ;
initLegacy(wells, state, dummy_state, pu) ;
}
/// One rate per phase and well connection.
@ -250,7 +424,8 @@ namespace Opm
std::vector<int>& currentControls() { return current_controls_; }
const std::vector<int>& currentControls() const { return current_controls_; }
data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override {
data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override
{
data::Wells res = WellState::report(pu, globalCellIdxMap);
const int nw = this->numWells();

2
opm/simulators/WellSwitchingLogger.hpp
View File

@ -69,7 +69,7 @@ public:
if( cc_.size() > 1 )
{
using Pair = typename SwitchMap::value_type;
switchMap_.insert(Pair(name, {char(from), char(to)}));
switchMap_.insert(Pair(name, {{char(from), char(to)}}));
}
else
{

374
opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp Normal file
View File

@ -0,0 +1,374 @@
/*
*/
#ifndef OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
#define OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP
#include <iostream>
#include <utility>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
#include <opm/simulators/timestepping/TimeStepControl.hpp>
namespace Opm {
// AdaptiveTimeStepping
//---------------------
class AdaptiveTimeSteppingEbos
{
template <class Solver>
class SolutionTimeErrorSolverWrapperEbos : public RelativeChangeInterface
{
const Solver& solver_;
public:
SolutionTimeErrorSolverWrapperEbos(const Solver& solver)
: solver_(solver)
{}
/// return || u^n+1 - u^n || / || u^n+1 ||
double relativeChange() const
{ return solver_.model().relativeChange(); }
};
template<class E>
void logException_(const E& exception, bool verbose)
{
if (verbose) {
std::string message;
message = "Caught Exception: ";
message += exception.what();
OpmLog::debug(message);
}
}
public:
//! \brief contructor taking parameter object
AdaptiveTimeSteppingEbos(const ParameterGroup& param,
const bool terminalOutput = true)
: timeStepControl_()
, restartFactor_( param.getDefault("solver.restartfactor", double(0.33)))
, growthFactor_( param.getDefault("solver.growthfactor", double(2)))
, maxGrowth_( param.getDefault("timestep.control.maxgrowth", double(3.0)))
// default is 1 year, convert to seconds
, maxTimeStep_( unit::convert::from(param.getDefault("timestep.max_timestep_in_days", 365.0), unit::day))
, solverRestartMax_( param.getDefault("solver.restart", int(10)))
, solverVerbose_( param.getDefault("solver.verbose", bool(true)) && terminalOutput)
, timestepVerbose_( param.getDefault("timestep.verbose", bool(true)) && terminalOutput)
, suggestedNextTimestep_( unit::convert::from(param.getDefault("timestep.initial_timestep_in_days", 1.0), unit::day))
, fullTimestepInitially_( param.getDefault("full_timestep_initially", bool(false)))
, timestepAfterEvent_( unit::convert::from(param.getDefault("timestep.timestep_in_days_after_event", -1.0), unit::day))
, useNewtonIteration_(false)
{
init_(param);
}
//! \brief contructor taking parameter object
//! \param tuning Pointer to ecl TUNING keyword
//! \param timeStep current report step
//! \param param The parameter object
AdaptiveTimeSteppingEbos(const Tuning& tuning,
size_t timeStep,
const ParameterGroup& param,
const bool terminalOutput = true)
: timeStepControl_()
, restartFactor_(tuning.getTSFCNV(timeStep))
, growthFactor_(tuning.getTFDIFF(timeStep))
, maxGrowth_(tuning.getTSFMAX(timeStep))
// default is 1 year, convert to seconds
, maxTimeStep_(tuning.getTSMAXZ(timeStep))
, solverRestartMax_(param.getDefault("solver.restart", int(10)))
, solverVerbose_(param.getDefault("solver.verbose", bool(true)) && terminalOutput)
, timestepVerbose_(param.getDefault("timestep.verbose", bool(true)) && terminalOutput)
, suggestedNextTimestep_(tuning.getTSINIT(timeStep))
, fullTimestepInitially_(param.getDefault("full_timestep_initially", bool(false)))
, timestepAfterEvent_(tuning.getTMAXWC(timeStep))
, useNewtonIteration_(false)
{
init_(param);
}
/** \brief step method that acts like the solver::step method
in a sub cycle of time steps
*/
template <class Solver, class Output>
SimulatorReport step(const SimulatorTimer& simulatorTimer,
Solver& solver,
const bool isEvent,
Output& outputWriter,
const std::vector<int>* fipnum = nullptr)
{
SimulatorReport report;
const double timestep = simulatorTimer.currentStepLength();
// init last time step as a fraction of the given time step
if (suggestedNextTimestep_ < 0) {
suggestedNextTimestep_ = restartFactor_ * timestep;
}
if (fullTimestepInitially_) {
suggestedNextTimestep_ = timestep;
}
// use seperate time step after event
if (isEvent && timestepAfterEvent_ > 0) {
suggestedNextTimestep_ = timestepAfterEvent_;
}
// create adaptive step timer with previously used sub step size
AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
// reset the statistics for the failed substeps
failureReport_ = SimulatorReport();
// counter for solver restarts
int restarts = 0;
// sub step time loop
while (!substepTimer.done()) {
// get current delta t
const double dt = substepTimer.currentStepLength() ;
if (timestepVerbose_) {
std::ostringstream ss;
ss <<"\nTime step " << substepTimer.currentStepNum() << ", stepsize "
<< unit::convert::to(substepTimer.currentStepLength(), unit::day) << " days.";
OpmLog::info(ss.str());
}
SimulatorReport substepReport;
std::string causeOfFailure = "";
try {
substepReport = solver.step(substepTimer);
report += substepReport;
if (solverVerbose_) {
// report number of linear iterations
OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
}
}
catch (const Opm::TooManyIterations& e) {
substepReport += solver.failureReport();
causeOfFailure = "Solver convergence failure - Iteration limit reached";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const Opm::LinearSolverProblem& e) {
substepReport += solver.failureReport();
causeOfFailure = "Linear solver convergence failure";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const Opm::NumericalIssue& e) {
substepReport += solver.failureReport();
causeOfFailure = "Solver convergence failure - Numerical problem encountered";
logException_(e, solverVerbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const std::runtime_error& e) {
substepReport += solver.failureReport();
logException_(e, solverVerbose_);
// also catch linear solver not converged
}
catch (const Dune::ISTLError& e) {
substepReport += solver.failureReport();
logException_(e, solverVerbose_);
// also catch errors in ISTL AMG that occur when time step is too large
}
catch (const Dune::MatrixBlockError& e) {
substepReport += solver.failureReport();
logException_(e, solverVerbose_);
// this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
}
if (substepReport.converged) {
// advance by current dt
++substepTimer;
// create object to compute the time error, simply forwards the call to the model
SolutionTimeErrorSolverWrapperEbos<Solver> relativeChange(solver);
// compute new time step estimate
const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
: substepReport.total_linear_iterations;
double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
substepTimer.simulationTimeElapsed());
// limit the growth of the timestep size by the growth factor
dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
// further restrict time step size growth after convergence problems
if (restarts > 0) {
dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
// solver converged, reset restarts counter
restarts = 0;
}
if (timestepVerbose_) {
std::ostringstream ss;
substepReport.reportStep(ss);
OpmLog::info(ss.str());
}
// write data if outputWriter was provided
// if the time step is done we do not need
// to write it as this will be done by the simulator
// anyway.
if (!substepTimer.done()) {
if (fipnum) {
solver.computeFluidInPlace(*fipnum);
}
Opm::time::StopWatch perfTimer;
perfTimer.start();
bool substep = true;
const auto& physicalModel = solver.model();
outputWriter.writeTimeStep(substepTimer, physicalModel, substep, -1.0, substepReport);
report.output_write_time += perfTimer.secsSinceStart();
}
// set new time step length
substepTimer.provideTimeStepEstimate(dtEstimate);
report.converged = substepTimer.done();
substepTimer.setLastStepFailed(false);
}
else { // in case of no convergence (linearIterations < 0)
substepTimer.setLastStepFailed(true);
failureReport_ += substepReport;
// increase restart counter
if (restarts >= solverRestartMax_) {
const auto msg = std::string("Solver failed to converge after cutting timestep ")
+ std::to_string(restarts) + " times.";
if (solverVerbose_) {
OpmLog::error(msg);
}
OPM_THROW_NOLOG(Opm::NumericalIssue, msg);
}
const double newTimeStep = restartFactor_ * dt;
// we need to revise this
substepTimer.provideTimeStepEstimate(newTimeStep);
if (solverVerbose_) {
std::string msg;
msg = causeOfFailure + "\nTimestep chopped to "
+ std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day)) + " days\n";
OpmLog::problem(msg);
}
++restarts;
}
}
// store estimated time step for next reportStep
suggestedNextTimestep_ = substepTimer.currentStepLength();
if (timestepVerbose_) {
std::ostringstream ss;
substepTimer.report(ss);
ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
OpmLog::debug(ss.str());
}
if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
suggestedNextTimestep_ = timestep;
}
return report;
}
/** \brief Returns the simulator report for the failed substeps of the last
* report step.
*/
const SimulatorReport& failureReport() const
{ return failureReport_; };
double suggestedNextStep() const
{ return suggestedNextTimestep_; }
void setSuggestedNextStep(const double x)
{ suggestedNextTimestep_ = x; }
void updateTUNING(const Tuning& tuning, size_t timeStep)
{
restartFactor_ = tuning.getTSFCNV(timeStep);
growthFactor_ = tuning.getTFDIFF(timeStep);
maxGrowth_ = tuning.getTSFMAX(timeStep);
maxTimeStep_ = tuning.getTSMAXZ(timeStep);
suggestedNextTimestep_ = tuning.getTSINIT(timeStep);
timestepAfterEvent_ = tuning.getTMAXWC(timeStep);
}
protected:
void init_(const ParameterGroup& param)
{
// valid are "pid" and "pid+iteration"
std::string control = param.getDefault("timestep.control", std::string("pid"));
// iterations is the accumulation of all linear iterations over all newton steops per time step
const int defaultTargetIterations = 30;
const int defaultTargetNewtonIterations = 8;
const double tol = param.getDefault("timestep.control.tol", double(1e-1));
if (control == "pid") {
timeStepControl_ = TimeStepControlType(new PIDTimeStepControl(tol));
}
else if (control == "pid+iteration") {
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations);
timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
}
else if (control == "pid+newtoniteration") {
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetNewtonIterations);
timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
useNewtonIteration_ = true;
}
else if (control == "iterationcount") {
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations);
const double decayrate = param.getDefault("timestep.control.decayrate", double(0.75));
const double growthrate = param.getDefault("timestep.control.growthrate", double(1.25));
timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
}
else if (control == "hardcoded") {
const std::string filename = param.getDefault("timestep.control.filename", std::string("timesteps"));
timeStepControl_ = TimeStepControlType(new HardcodedTimeStepControl(filename));
}
else
OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control);
// make sure growth factor is something reasonable
assert(growthFactor_ >= 1.0);
}
typedef std::unique_ptr<TimeStepControlInterface> TimeStepControlType;
SimulatorReport failureReport_; //!< statistics for the failed substeps of the last timestep
TimeStepControlType timeStepControl_; //!< time step control object
double restartFactor_; //!< factor to multiply time step with when solver fails to converge
double growthFactor_; //!< factor to multiply time step when solver recovered from failed convergence
double maxGrowth_; //!< factor that limits the maximum growth of a time step
double maxTimeStep_; //!< maximal allowed time step size
const int solverRestartMax_; //!< how many restart of solver are allowed
const bool solverVerbose_; //!< solver verbosity
const bool timestepVerbose_; //!< timestep verbosity
double suggestedNextTimestep_; //!< suggested size of next timestep
bool fullTimestepInitially_; //!< beginning with the size of the time step from data file
double timestepAfterEvent_; //!< suggested size of timestep after an event
bool useNewtonIteration_; //!< use newton iteration count for adaptive time step control
};
}
#endif // OPM_ADAPTIVE_TIME_STEPPING_EBOS_HPP