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

extend and clean up the SimulatorReport

Browse Source
This commit is contained in:
Andreas Lauser
2016-11-23 21:44:33 +01:00
parent 0429c756ba
commit 8c5f92dbc4
18 changed files with 244 additions and 198 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
opm/autodiff/BlackoilModelBase.hpp
View File

@@ -37,6 +37,7 @@
#include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <array>
@@ -175,7 +176,7 @@ namespace Opm {
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
template <class NonlinearSolverType>
IterationReport nonlinearIteration(const int iteration,
SimulatorReport nonlinearIteration(const int iteration,
const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver,
ReservoirState& reservoir_state,
@@ -194,8 +195,7 @@ namespace Opm {
/// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
/// \return well iterations.
IterationReport
SimulatorReport
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly);
@@ -396,7 +396,7 @@ namespace Opm {
assembleMassBalanceEq(const SolutionState& state);
IterationReport
SimulatorReport
solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
const ReservoirState& reservoir_state,

75
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@@ -50,8 +50,10 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <dune/common/timer.hh>
#include <cassert>
#include <cmath>
#include <iostream>
@@ -233,7 +235,7 @@ typedef Eigen::Array<double,
template <class Grid, class WellModel, class Implementation>
template <class NonlinearSolverType>
IterationReport
SimulatorReport
BlackoilModelBase<Grid, WellModel, Implementation>::
nonlinearIteration(const int iteration,
const SimulatorTimerInterface& timer,
@@ -241,6 +243,10 @@ typedef Eigen::Array<double,
ReservoirState& reservoir_state,
WellState& well_state)
{
SimulatorReport report;
Dune::Timer perfTimer;
perfTimer.start();
const double dt = timer.currentStepLength();
if (iteration == 0) {
@@ -250,16 +256,46 @@ typedef Eigen::Array<double,
current_relaxation_ = 1.0;
dx_old_ = V::Zero(sizeNonLinear());
}
IterationReport iter_report = asImpl().assemble(reservoir_state, well_state, iteration == 0);
try {
report += asImpl().assemble(reservoir_state, well_state, iteration == 0);
report.assemble_time += perfTimer.stop();
}
catch (...) {
report.assemble_time += perfTimer.stop();
throw;
}
report.total_linearizations = 1;
perfTimer.reset();
perfTimer.start();
report.converged = asImpl().getConvergence(timer, iteration);
residual_norms_history_.push_back(asImpl().computeResidualNorms());
const bool converged = asImpl().getConvergence(timer, iteration);
const bool must_solve = (iteration < nonlinear_solver.minIter()) || (!converged);
report.update_time += perfTimer.stop();
const bool must_solve = (iteration < nonlinear_solver.minIter()) || (!report.converged);
if (must_solve) {
perfTimer.reset();
perfTimer.start();
report.total_newton_iterations = 1;
// enable single precision for solvers when dt is smaller then maximal time step for single precision
residual_.singlePrecision = ( dt < param_.maxSinglePrecisionTimeStep_ );
// Compute the nonlinear update.
V dx = asImpl().solveJacobianSystem();
V dx;
try {
dx = asImpl().solveJacobianSystem();
report.linear_solve_time += perfTimer.stop();
report.total_linear_iterations += linearIterationsLastSolve();
}
catch (...) {
report.linear_solve_time += perfTimer.stop();
report.total_linear_iterations += linearIterationsLastSolve();
throw;
}
perfTimer.reset();
perfTimer.start();
if (param_.use_update_stabilization_) {
// Stabilize the nonlinear update.
@@ -281,10 +317,10 @@ typedef Eigen::Array<double,
// Apply the update, applying model-dependent
// limitations and chopping of the update.
asImpl().updateState(dx, reservoir_state, well_state);
report.update_time += perfTimer.stop();
}
const bool failed = false; // Not needed in this model.
const int linear_iters = must_solve ? asImpl().linearIterationsLastSolve() : 0;
return IterationReport{ failed, converged, linear_iters , iter_report.well_iterations};
return report;
}
@@ -702,7 +738,7 @@ typedef Eigen::Array<double,
template <class Grid, class WellModel, class Implementation>
IterationReport
SimulatorReport
BlackoilModelBase<Grid, WellModel, Implementation>::
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
@@ -710,6 +746,8 @@ typedef Eigen::Array<double,
{
using namespace Opm::AutoDiffGrid;
SimulatorReport report;
// If we have VFP tables, we need the well connection
// pressures for the "simple" hydrostatic correction
// between well depth and vfp table depth.
@@ -757,9 +795,8 @@ typedef Eigen::Array<double,
asImpl().assembleMassBalanceEq(state);
// -------- Well equations ----------
IterationReport iter_report = {false, false, 0, 0};
if ( ! wellsActive() ) {
return iter_report;
return report;
}
std::vector<ADB> mob_perfcells;
@@ -767,7 +804,7 @@ typedef Eigen::Array<double,
asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
report += asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
}
V aliveWells;
std::vector<ADB> cq_s;
@@ -783,7 +820,7 @@ typedef Eigen::Array<double,
asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
}
return iter_report;
return report;
}
@@ -962,7 +999,7 @@ typedef Eigen::Array<double,
template <class Grid, class WellModel, class Implementation>
IterationReport
SimulatorReport
BlackoilModelBase<Grid, WellModel, Implementation>::
solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
@@ -1072,9 +1109,11 @@ typedef Eigen::Array<double,
if (!converged) {
well_state = well_state0;
}
const bool failed = false; // Not needed in this method.
const int linear_iters = 0; // Not needed in this method
return IterationReport{failed, converged, linear_iters, it};
SimulatorReport report;
report.total_well_iterations = it;
report.converged = converged;
return report;
}

77
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -45,6 +45,7 @@
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/core/grid.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
@@ -64,6 +65,7 @@
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/common/parallel/collectivecommunication.hh>
#include <dune/common/timer.hh>
#include <dune/common/unused.hh>
#include <cassert>
@@ -247,12 +249,16 @@ namespace Opm {
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
template <class NonlinearSolverType>
IterationReport nonlinearIteration(const int iteration,
SimulatorReport nonlinearIteration(const int iteration,
const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver,
ReservoirState& reservoir_state,
WellState& well_state)
{
SimulatorReport report;
Dune::Timer perfTimer;
perfTimer.start();
if (iteration == 0) {
// 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.
@@ -264,16 +270,30 @@ namespace Opm {
// reset intensive quantities cache useless other options are set
// further down
invalidateIntensiveQuantitiesCache_ = true;
report.total_linearizations = 1;
try {
report += assemble(timer, iteration, reservoir_state, well_state);
report.assemble_time += perfTimer.stop();
}
catch (...) {
report.assemble_time += perfTimer.stop();
// todo (?): make the report an attribute of the class
throw; // continue throwing the stick
}
IterationReport iter_report = assemble(timer, iteration, reservoir_state, well_state);
std::vector<double> residual_norms;
const bool converged = getConvergence(timer, iteration,residual_norms);
perfTimer.reset();
perfTimer.start();
report.converged = getConvergence(timer, iteration,residual_norms);
report.update_time += perfTimer.stop();
residual_norms_history_.push_back(residual_norms);
bool must_solve = (iteration < nonlinear_solver.minIter()) || (!converged);
// don't solve if we have reached the maximum number of iteration.
must_solve = must_solve && (iteration < nonlinear_solver.maxIter());
bool must_solve = iteration < nonlinear_solver.minIter() || !report.converged;
if (must_solve) {
perfTimer.reset();
perfTimer.start();
report.total_newton_iterations = 1;
// enable single precision for solvers when dt is smaller then 20 days
//residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
@@ -284,7 +304,20 @@ namespace Opm {
BVector x(nc);
BVector xw(nw);
solveJacobianSystem(x, xw);
try {
solveJacobianSystem(x, xw);
report.linear_solve_time += perfTimer.stop();
report.total_linear_iterations += linearIterationsLastSolve();
}
catch (...) {
report.linear_solve_time += perfTimer.stop();
report.total_linear_iterations += linearIterationsLastSolve();
// todo (?): make the report an attribute of the class
throw; // re-throw up
}
perfTimer.reset();
perfTimer.start();
// Stabilize the nonlinear update.
bool isOscillate = false;
@@ -301,24 +334,20 @@ namespace Opm {
}
nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
// Apply the update, applying model-dependent
// limitations and chopping of the update.
// Apply the update, with considering model-dependent limitations and
// chopping of the update.
updateState(x,reservoir_state);
wellModel().updateWellState(xw, well_state);
// since the solution was changed, the cache for the intensive quantities are invalid
// ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
report.update_time += perfTimer.stop();
}
if( converged && (iteration >= nonlinear_solver.minIter()) )
{
// in case of convergence we do not need to reset intensive quantities
else {
// if the solution is not updated, we do not need to recalculate the
// intensive quantities in the next iteration.
assert(report.converged);
invalidateIntensiveQuantitiesCache_ = false ;
}
const bool failed = false; // Not needed in this model.
const int linear_iters = must_solve ? linearIterationsLastSolve() : 0;
return IterationReport{ failed, converged, linear_iters, iter_report.well_iterations };
return report;
}
void printIf(int c, double x, double y, double eps, std::string type) {
@@ -346,30 +375,31 @@ namespace Opm {
/// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
IterationReport assemble(const SimulatorTimerInterface& timer,
SimulatorReport assemble(const SimulatorTimerInterface& timer,
const int iterationIdx,
const ReservoirState& reservoir_state,
WellState& well_state)
{
using namespace Opm::AutoDiffGrid;
SimulatorReport report;
// -------- Mass balance equations --------
assembleMassBalanceEq(timer, iterationIdx, reservoir_state);
// -------- Well equations ----------
double dt = timer.currentStepLength();
IterationReport iter_report;
try
{
iter_report = wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state);
report = wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state);
}
catch ( const Dune::FMatrixError& e )
{
OPM_THROW(Opm::NumericalProblem,"Well equation did not converge");
}
return iter_report;
return report;
}
@@ -439,7 +469,6 @@ namespace Opm {
wellModel().applyScaleAdd(alpha, x, y);
}
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
void solveJacobianSystem(BVector& x, BVector& xw) const

4
opm/autodiff/BlackoilMultiSegmentModel.hpp
View File

@@ -106,7 +106,7 @@ namespace Opm {
/// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
IterationReport
SimulatorReport
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly);
@@ -166,7 +166,7 @@ namespace Opm {
const MultisegmentWells::MultisegmentWellOps& msWellOps() const { return well_model_.wellOps(); }
IterationReport
SimulatorReport
solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
const ReservoirState& reservoir_state,

19
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
View File

@@ -127,7 +127,7 @@ namespace Opm {
template <class Grid>
IterationReport
SimulatorReport
BlackoilMultiSegmentModel<Grid>::
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
@@ -188,9 +188,9 @@ namespace Opm {
asImpl().assembleMassBalanceEq(state);
// -------- Well equations ----------
IterationReport iter_report = {false, false, 0, 0};
if ( ! wellsActive() ) {
return iter_report;
SimulatorReport report;
return report;
}
wellModel().computeSegmentFluidProperties(state);
@@ -200,10 +200,11 @@ namespace Opm {
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
SimulatorReport report;
wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
}
// the perforation flux here are different
@@ -215,7 +216,7 @@ namespace Opm {
wellModel().addWellFluxEq(cq_s, state, residual_);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
return iter_report;
return report;
}
@@ -223,16 +224,16 @@ namespace Opm {
template <class Grid>
IterationReport
SimulatorReport
BlackoilMultiSegmentModel<Grid>::solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
const ReservoirState& reservoir_state,
SolutionState& state,
WellState& well_state)
{
IterationReport iter_report = Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
SimulatorReport report = Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
if (iter_report.converged) {
if (report.converged) {
// We must now update the state.segp and state.segqs members,
// that the base version does not know about.
const int np = numPhases();
@@ -261,7 +262,7 @@ namespace Opm {
asImpl().computeWellConnectionPressures(state, well_state);
}
return iter_report;
return report;
}

8
opm/autodiff/BlackoilPressureModel.hpp
View File

@@ -149,13 +149,14 @@ namespace Opm {
IterationReport
SimulatorReport
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly)
{
IterationReport iter_report = Base::assemble(reservoir_state, well_state, initial_assembly);
SimulatorReport report;
report += Base::assemble(reservoir_state, well_state, initial_assembly);
if (initial_assembly) {
}
@@ -186,7 +187,8 @@ namespace Opm {
assert(int(well_state.perfRates().size()) == wflux.size());
std::copy_n(wflux.data(), wflux.size(), well_state.perfRates().begin());
}
return iter_report;
return report;
}

26
opm/autodiff/BlackoilSequentialModel.hpp
View File

@@ -123,7 +123,7 @@ namespace Opm {
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
template <class NonlinearSolverType>
IterationReport nonlinearIteration(const int iteration,
SimulatorReport nonlinearIteration(const int iteration,
const SimulatorTimerInterface& timer,
NonlinearSolverType& /* nonlinear_solver */,
ReservoirState& reservoir_state,
@@ -140,7 +140,8 @@ namespace Opm {
OpmLog::info("Solving the pressure equation.");
}
ReservoirState initial_state = reservoir_state;
const int pressure_liniter = pressure_solver_.step(timer, reservoir_state, well_state);
const SimulatorReport pressure_report = pressure_solver_.step(timer, reservoir_state, well_state);
const int pressure_liniter = pressure_report.total_linear_iterations;
if (pressure_liniter == -1) {
OPM_THROW(std::runtime_error, "Pressure solver failed to converge.");
}
@@ -149,13 +150,17 @@ namespace Opm {
if (terminalOutputEnabled()) {
OpmLog::info("Solving the transport equations.");
}
const int transport_liniter = transport_solver_.step(timer, initial_state, well_state, reservoir_state, well_state);
const SimulatorReport transport_report = transport_solver_.step(timer, initial_state, well_state, reservoir_state, well_state);
const int transport_liniter = transport_report.total_linear_iterations;
if (transport_liniter == -1) {
OPM_THROW(std::runtime_error, "Transport solver failed to converge.");
}
// Report and return.
return IterationReport { false, true, pressure_liniter + transport_liniter, 0 };
SimulatorReport report;
report.converged = true;
report.total_linear_iterations = pressure_liniter + transport_liniter;
return report;
} else {
// Iterate to fully implicit solution.
// This call is just for a single iteration (one pressure and one transport solve),
@@ -168,7 +173,8 @@ namespace Opm {
if (terminalOutputEnabled()) {
OpmLog::info("Solving the pressure equation.");
}
const int pressure_liniter = pressure_solver_.step(timer, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
const SimulatorReport& pressure_report = pressure_solver_.step(timer, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
const int pressure_liniter = pressure_report.total_linear_iterations;
if (pressure_liniter == -1) {
OPM_THROW(std::runtime_error, "Pressure solver failed to converge.");
}
@@ -177,14 +183,16 @@ namespace Opm {
if (terminalOutputEnabled()) {
OpmLog::info("Solving the transport equations.");
}
const int transport_liniter = transport_solver_.step(timer, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
const SimulatorReport& transport_report = transport_solver_.step(timer, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
const int transport_liniter = transport_report.total_linear_iterations;
if (transport_liniter == -1) {
OPM_THROW(std::runtime_error, "Transport solver failed to converge.");
}
// Report and return.
const bool converged = iteration >= 3; // TODO: replace this with a proper convergence check
return IterationReport { false, converged, pressure_liniter + transport_liniter, 0 };
SimulatorReport report;
report.converged = iteration >= 3; // TODO: replace this with a proper convergence check
report.total_linear_iterations = pressure_liniter + transport_liniter;
return report;
}
}

13
opm/autodiff/BlackoilTransportModel.hpp
View File

@@ -83,14 +83,15 @@ namespace Opm {
asImpl().makeConstantState(state0_);
}
IterationReport
SimulatorReport
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly)
{
using namespace Opm::AutoDiffGrid;
SimulatorReport report;
// If we have VFP tables, we need the well connection
// pressures for the "simple" hydrostatic correction
// between well depth and vfp table depth.
@@ -125,9 +126,8 @@ namespace Opm {
asImpl().assembleMassBalanceEq(state);
// -------- Well equations ----------
IterationReport iter_report = {false, false, 0, 0};
if ( ! wellsActive() ) {
return iter_report;
return report;
}
std::vector<ADB> mob_perfcells;
@@ -135,7 +135,7 @@ namespace Opm {
asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
report += asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
}
V aliveWells;
std::vector<ADB> cq_s;
@@ -154,7 +154,8 @@ namespace Opm {
asImpl().makeConstantState(state0);
asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
}
return iter_report;
return report;
}

8
opm/autodiff/NonlinearSolver.hpp
View File

@@ -22,6 +22,7 @@
#define OPM_NONLINEARSOLVER_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/autodiff/DuneMatrix.hpp>
@@ -87,8 +88,7 @@ namespace Opm {
/// \param[in] timer simulation timer
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \return number of linear iterations used
int
SimulatorReport
step(const SimulatorTimerInterface& timer,
ReservoirState& reservoir_state,
WellState& well_state);
@@ -103,7 +103,7 @@ namespace Opm {
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \return number of linear iterations used
int
SimulatorReport
step(const SimulatorTimerInterface& timer,
const ReservoirState& initial_reservoir_state,
const WellState& initial_well_state,
@@ -178,7 +178,7 @@ namespace Opm {
double relaxRelTol() const { return param_.relax_rel_tol_; }
/// The maximum number of nonlinear iterations allowed.
double maxIter() const { return param_.max_iter_; }
int maxIter() const { return param_.max_iter_; }
/// The minimum number of nonlinear iterations allowed.
double minIter() const { return param_.min_iter_; }

40
opm/autodiff/NonlinearSolver_impl.hpp
View File

@@ -102,7 +102,7 @@ namespace Opm
}
template <class PhysicalModel>
int
SimulatorReport
NonlinearSolver<PhysicalModel>::
step(const SimulatorTimerInterface& timer,
ReservoirState& reservoir_state,
@@ -114,7 +114,7 @@ namespace Opm
template <class PhysicalModel>
int
SimulatorReport
NonlinearSolver<PhysicalModel>::
step(const SimulatorTimerInterface& timer,
const ReservoirState& initial_reservoir_state,
@@ -122,6 +122,9 @@ namespace Opm
ReservoirState& reservoir_state,
WellState& well_state)
{
SimulatorReport iterReport;
SimulatorReport report;
// Do model-specific once-per-step calculations.
model_->prepareStep(timer, initial_reservoir_state, initial_well_state);
@@ -130,44 +133,31 @@ namespace Opm
// Let the model do one nonlinear iteration.
// Set up for main solver loop.
int linIters = 0;
bool converged = false;
int wellIters = 0;
// ---------- Main nonlinear solver loop ----------
do {
// 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.
IterationReport report = model_->nonlinearIteration(iteration, timer, *this, reservoir_state, well_state);
if (report.failed) {
OPM_THROW(Opm::NumericalProblem, "Failed to complete a nonlinear iteration.");
}
iterReport = model_->nonlinearIteration(iteration, timer, *this, reservoir_state, well_state);
report += iterReport;
report.converged = iterReport.converged;
converged = report.converged;
linIters += report.linear_iterations;
wellIters += report.well_iterations;
++iteration;
} while ( (!converged && (iteration <= maxIter())) || (iteration <= minIter()));
iteration += 1;
} while ( (!converged && (iteration <= maxIter())) || (iteration < minIter()));
if (!converged) {
if (model_->terminalOutputEnabled()) {
std::cerr << "WARNING: Failed to compute converged solution in " << iteration - 1 << " iterations." << std::endl;
}
return -1; // -1 indicates that the solver has to be restarted
OPM_THROW(Opm::NumericalProblem, "Failed to complete a time step within "+std::to_string(maxIter())+" iterations.");
}
linearIterations_ += linIters;
nonlinearIterations_ += iteration - 1; // Since the last one will always be trivial.
linearizations_ += iteration;
wellIterations_ += wellIters;
linearIterationsLast_ = linIters;
nonlinearIterationsLast_ = iteration;
wellIterationsLast_ = wellIters;
// Do model-specific post-step actions.
model_->afterStep(timer, reservoir_state, well_state);
report.converged = true;
return linIters;
return report;
}

51
opm/autodiff/SimulatorBase_impl.hpp
View File

@@ -98,7 +98,6 @@ namespace Opm
// Create timers and file for writing timing info.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
@@ -133,12 +132,11 @@ namespace Opm
desiredRestoreStep );
}
unsigned int totalLinearizations = 0;
unsigned int totalNonlinearIterations = 0;
unsigned int totalLinearIterations = 0;
bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false );
std::vector<double> well_potentials;
DynamicListEconLimited dynamic_list_econ_limited;
SimulatorReport report;
SimulatorReport stepReport;
bool ooip_computed = false;
std::vector<int> fipnum_global = eclipse_state_->get3DProperties().getIntGridProperty("FIPNUM").getData();
@@ -186,9 +184,14 @@ namespace Opm
// write the inital state at the report stage
if (timer.initialStep()) {
Dune::Timer perfTimer;
perfTimer.start();
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state );
report.output_write_time += perfTimer.stop();
}
// Max oil saturation (for VPPARS), hysteresis update.
@@ -220,8 +223,7 @@ namespace Opm
step_msg << "\nTime step " << std::setw(4) <<timer.currentStepNum()
<< " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
<< "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
<< ", date = " << timer.currentDateTime()
<< "\n";
<< ", date = " << timer.currentDateTime();
OpmLog::info(step_msg.str());
}
@@ -231,11 +233,12 @@ 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_ );
report += adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ );
}
else {
// solve for complete report step
solver->step(timer, state, well_state);
stepReport = solver->step(timer, state, well_state);
report += stepReport;
if( terminal_output_ )
{
@@ -268,13 +271,9 @@ namespace Opm
// take time that was used to solve system for this reportStep
solver_timer.stop();
// accumulate the number of nonlinear and linear Iterations
totalLinearizations += solver->linearizations();
totalNonlinearIterations += solver->nonlinearIterations();
totalLinearIterations += solver->linearIterations();
// update timing.
report.solver_time += solver_timer.secsSinceStart();
// Report timing.
const double st = solver_timer.secsSinceStart();
// Compute current FIP.
std::vector<V> COIP;
COIP = solver->computeFluidInPlace(state, fipnum);
@@ -290,27 +289,17 @@ namespace Opm
}
}
// accumulate total time
stime += st;
if ( terminal_output_ )
{
std::string msg;
msg = "Fully implicit solver took: " + std::to_string(st) + " seconds. Total solver time taken: " + std::to_string(stime) + " seconds.";
msg = "Fully implicit solver took: " + std::to_string(stepReport.solver_time) + " seconds. Total solver time taken: " + std::to_string(report.solver_time) + " seconds.";
OpmLog::note(msg);
}
if ( output_writer_.output() ) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.total_newton_iterations = solver->nonlinearIterations();
step_report.total_linear_iterations = solver->linearIterations();
step_report.total_linearizations = solver->linearizations();
if ( output_writer_.isIORank() )
{
step_report.reportParam(tstep_os);
stepReport.reportParam(tstep_os);
}
}
@@ -318,8 +307,11 @@ namespace Opm
++timer;
// write simulation state at the report stage
Dune::Timer perfTimer;
perfTimer.start();
const auto& physicalModel = solver->model();
output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
report.output_write_time += perfTimer.stop();
prev_well_state = well_state;
// The well potentials are only computed if they are needed
@@ -334,13 +326,8 @@ namespace Opm
// Stop timer and create timing report
total_timer.stop();
SimulatorReport report;
report.pressure_time = stime;
report.transport_time = 0.0;
report.total_time = total_timer.secsSinceStart();
report.total_linearizations = totalLinearizations;
report.total_newton_iterations = totalNonlinearIterations;
report.total_linear_iterations = totalLinearIterations;
report.converged = true;
return report;
}

64
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@@ -161,7 +161,6 @@ public:
// Create timers and file for writing timing info.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
@@ -188,12 +187,11 @@ public:
// desiredRestoreStep );
}
unsigned int totalLinearizations = 0;
unsigned int totalNonlinearIterations = 0;
unsigned int totalLinearIterations = 0;
bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false );
std::vector<double> well_potentials;
DynamicListEconLimited dynamic_list_econ_limited;
SimulatorReport report;
SimulatorReport stepReport;
bool ooip_computed = false;
std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
@@ -253,22 +251,25 @@ public:
// write the inital state at the report stage
if (timer.initialStep()) {
Dune::Timer perfTimer;
perfTimer.start();
// calculate Intensive Quantities
// make sure that the Intensive Quantities cache is up to date
const auto& gridManager = ebosSimulator_.gridManager();
const auto& gridView = gridManager.gridView();
auto elemIt = gridView.template begin<0>();
auto elemEndIt = gridView.template end<0>();
ElementContext elemCtx(ebosSimulator_);
for (; elemIt != elemEndIt; ++ elemIt) {
// this is convenient, but slightly inefficient: one only needs to update
// the primary intensive quantities
elemCtx.updateAll(*elemIt);
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
}
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state );
report.output_write_time += perfTimer.stop();
}
// Compute orignal FIP;
@@ -283,11 +284,12 @@ public:
std::ostringstream step_msg;
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 << "\nTime step " << std::setw(4) <<timer.currentStepNum()
step_msg << "\n"
<< "Time step " << std::setw(4) <<timer.currentStepNum()
<< " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
<< "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
<< ", date = " << timer.currentDateTime()
<< "\n";
<< ", size = " << (double)unit::convert::to(timer.currentStepLength(), unit::day) << " days";
OpmLog::info(step_msg.str());
}
@@ -299,14 +301,16 @@ public:
// \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_ );
report += adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ );
}
else {
// solve for complete report step
solver->step(timer, state, well_state);
stepReport = solver->step(timer, state, well_state);
report += stepReport;
if( terminal_output_ )
{
//stepReport.briefReport();
std::ostringstream iter_msg;
iter_msg << "Stepsize " << (double)unit::convert::to(timer.currentStepLength(), unit::day);
if (solver->wellIterations() != 0) {
@@ -324,13 +328,8 @@ public:
// take time that was used to solve system for this reportStep
solver_timer.stop();
// accumulate the number of nonlinear and linear Iterations
totalLinearizations += solver->linearizations();
totalNonlinearIterations += solver->nonlinearIterations();
totalLinearIterations += solver->linearIterations();
// Report timing.
const double st = solver_timer.secsSinceStart();
// update timing.
report.solver_time += solver_timer.secsSinceStart();
// Compute current FIP.
std::vector<std::vector<double>> COIP;
@@ -345,30 +344,22 @@ public:
for (size_t reg = 0; reg < OOIP.size(); ++reg) {
outputFluidInPlace(OOIP[reg], COIP[reg], eclState().getUnits(), reg+1);
}
}
// accumulate total time
stime += st;
if ( terminal_output_ )
{
std::string msg;
msg = "Fully implicit solver took: " + std::to_string(st) + " seconds. Total solver time taken: " + std::to_string(stime) + " seconds.";
msg =
"Time step took " + std::to_string(stepReport.solver_time) + " seconds; "
"total solver time " + std::to_string(report.solver_time) + " seconds.";
OpmLog::note(msg);
}
if ( output_writer_.output() ) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.reportParam(tstep_os);
}
// Increment timer, remember well state.
++timer;
// write simulation state at the report stage
Dune::Timer perfTimer;
perfTimer.start();
output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
report.output_write_time += perfTimer.stop();
prev_well_state = well_state;
// The well potentials are only computed if they are needed
@@ -383,13 +374,8 @@ public:
// Stop timer and create timing report
total_timer.stop();
SimulatorReport report;
report.pressure_time = stime;
report.transport_time = 0.0;
report.total_time = total_timer.secsSinceStart();
report.total_linearizations = totalLinearizations;
report.total_newton_iterations = totalNonlinearIterations;
report.total_linear_iterations = totalLinearIterations;
report.converged = true;
return report;
}

5
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
View File

@@ -179,7 +179,7 @@ namespace Opm
stime += st;
if ( output_writer_.output() ) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.solver_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.reportParam(tstep_os);
}
@@ -198,9 +198,8 @@ namespace Opm
// Stop timer and create timing report
total_timer.stop();
SimulatorReport report;
report.pressure_time = stime;
report.transport_time = 0.0;
report.total_time = total_timer.secsSinceStart();
report.solver_time = stime;
report.total_newton_iterations = totalNonlinearIterations;
report.total_linear_iterations = totalLinearIterations;
return report;

27
opm/autodiff/StandardWellsDense.hpp
View File

@@ -160,14 +160,14 @@ namespace Opm {
template <typename Simulator>
IterationReport assemble(Simulator& ebosSimulator,
SimulatorReport assemble(Simulator& ebosSimulator,
const int iterationIdx,
const double dt,
WellState& well_state) {
IterationReport iter_report = {false, false, 0, 0};
SimulatorReport report;
if ( ! localWellsActive() ) {
return iter_report;
return report;
}
resetWellControlFromState(well_state);
@@ -182,7 +182,7 @@ namespace Opm {
if (param_.solve_welleq_initially_ && iterationIdx == 0) {
// solve the well equations as a pre-processing step
iter_report = solveWellEq(ebosSimulator, dt, well_state);
report = solveWellEq(ebosSimulator, dt, well_state);
}
assembleWellEq(ebosSimulator, dt, well_state, false);
@@ -190,7 +190,8 @@ namespace Opm {
if (param_.compute_well_potentials_) {
//wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
}
return iter_report;
report.converged = true;
return report;
}
template <typename Simulator>
@@ -637,9 +638,8 @@ namespace Opm {
}
}
template <typename Simulator>
IterationReport solveWellEq(Simulator& ebosSimulator,
SimulatorReport solveWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state)
{
@@ -672,9 +672,10 @@ namespace Opm {
well_state = well_state0;
}
const bool failed = false; // Not needed in this method.
const int linear_iters = 0; // Not needed in this method
return IterationReport{failed, converged, linear_iters, it};
SimulatorReport report;
report.converged = converged;
report.total_well_iterations = it;
return report;
}
void printIf(int c, double x, double y, double eps, std::string type) {
@@ -1114,9 +1115,9 @@ namespace Opm {
// We disregard terminal_ouput here as with it only messages
// for wells on one process will be printed.
std::ostringstream ss;
ss << "Switching control mode for well " << wells().name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
ss << " Switching control mode for well " << wells().name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)];
OpmLog::info(ss.str());
xw.currentControls()[w] = ctrl_index;
current = xw.currentControls()[w];

2
opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
View File

@@ -121,7 +121,7 @@ namespace Opm {
/// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
IterationReport
SimulatorReport
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly);

11
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
View File

@@ -494,13 +494,15 @@ namespace Opm {
template <class Grid>
IterationReport
SimulatorReport
BlackoilPolymerModel<Grid>::assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly)
{
using namespace Opm::AutoDiffGrid;
SimulatorReport report;
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
// updateWellControls(well_state);
@@ -531,10 +533,9 @@ namespace Opm {
// -------- Mass balance equations --------
assembleMassBalanceEq(state);
IterationReport iter_report = {false, false, 0, 0};
// -------- Well equations ----------
if ( ! wellsActive() ) {
return iter_report;
return report;
}
std::vector<ADB> mob_perfcells;
@@ -572,7 +573,9 @@ namespace Opm {
wellModel().addWellFluxEq(cq_s, state, residual_);
addWellContributionToMassBalanceEq(cq_s, state, well_state);
wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
return iter_report;
report.converged = true;
return report;
}

2
opm/simulators/WellSwitchingLogger.cpp
View File

@@ -137,7 +137,7 @@ void WellSwitchingLogger::logSwitch(const char* name, std::array<char,2> fromto,
int rank)
{
std::ostringstream ss;
ss << "Switching control mode for well " << name
ss << " Switching control mode for well " << name
<< " from " << modestring[WellControlType(fromto[0])]
<< " to " << modestring[WellControlType(fromto[1])]
<< " on rank " << rank << std::endl;

2
opm/simulators/WellSwitchingLogger.hpp
View File

@@ -74,7 +74,7 @@ public:
else
{
std::ostringstream ss;
ss << "Switching control mode for well " << name
ss << " Switching control mode for well " << name
<< " from " << modestring[from]
<< " to " << modestring[to] << std::endl;
OpmLog::info(ss.str());