My Project
BlackoilModelEbos.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2014, 2015 Statoil ASA.
5  Copyright 2015 NTNU
6  Copyright 2015, 2016, 2017 IRIS AS
7 
8  This file is part of the Open Porous Media project (OPM).
9 
10  OPM is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  OPM is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with OPM. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 #ifndef OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
25 #define OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
26 
27 #include <ebos/eclproblem.hh>
28 #include <opm/models/utils/start.hh>
29 
30 #include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
31 
32 #include <opm/simulators/flow/NonlinearSolverEbos.hpp>
33 #include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
34 #include <opm/simulators/wells/BlackoilWellModel.hpp>
35 #include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>
36 #include <opm/simulators/wells/WellConnectionAuxiliaryModule.hpp>
37 #include <opm/simulators/flow/countGlobalCells.hpp>
38 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
39 
40 #include <opm/grid/UnstructuredGrid.h>
41 #include <opm/simulators/timestepping/SimulatorReport.hpp>
42 #include <opm/simulators/linalg/ParallelIstlInformation.hpp>
43 #include <opm/core/props/phaseUsageFromDeck.hpp>
44 #include <opm/common/ErrorMacros.hpp>
45 #include <opm/common/Exceptions.hpp>
46 #include <opm/common/OpmLog/OpmLog.hpp>
47 #include <opm/parser/eclipse/Units/Units.hpp>
48 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
49 #include <opm/common/utility/parameters/ParameterGroup.hpp>
50 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
51 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
52 
53 #include <opm/simulators/linalg/ISTLSolverEbos.hpp>
54 
55 #include <dune/istl/owneroverlapcopy.hh>
56 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7)
57 #include <dune/common/parallel/communication.hh>
58 #else
59 #include <dune/common/parallel/collectivecommunication.hh>
60 #endif
61 #include <dune/common/timer.hh>
62 #include <dune/common/unused.hh>
63 
64 #include <cassert>
65 #include <cmath>
66 #include <iostream>
67 #include <iomanip>
68 #include <limits>
69 #include <vector>
70 #include <algorithm>
71 
72 namespace Opm::Properties {
73 
74 namespace TTag {
76  using InheritsFrom = std::tuple<FlowTimeSteppingParameters, FlowModelParameters,
77  FlowNonLinearSolver, EclBaseProblem, BlackOilModel>;
78 };
79 }
80 template<class TypeTag>
81 struct OutputDir<TypeTag, TTag::EclFlowProblem> {
82  static constexpr auto value = "";
83 };
84 template<class TypeTag>
85 struct EnableDebuggingChecks<TypeTag, TTag::EclFlowProblem> {
86  static constexpr bool value = false;
87 };
88 // default in flow is to formulate the equations in surface volumes
89 template<class TypeTag>
90 struct BlackoilConserveSurfaceVolume<TypeTag, TTag::EclFlowProblem> {
91  static constexpr bool value = true;
92 };
93 template<class TypeTag>
94 struct UseVolumetricResidual<TypeTag, TTag::EclFlowProblem> {
95  static constexpr bool value = false;
96 };
97 
98 template<class TypeTag>
99 struct EclAquiferModel<TypeTag, TTag::EclFlowProblem> {
101 };
102 
103 // disable all extensions supported by black oil model. this should not really be
104 // necessary but it makes things a bit more explicit
105 template<class TypeTag>
106 struct EnablePolymer<TypeTag, TTag::EclFlowProblem> {
107  static constexpr bool value = false;
108 };
109 template<class TypeTag>
110 struct EnableSolvent<TypeTag, TTag::EclFlowProblem> {
111  static constexpr bool value = false;
112 };
113 template<class TypeTag>
114 struct EnableTemperature<TypeTag, TTag::EclFlowProblem> {
115  static constexpr bool value = true;
116 };
117 template<class TypeTag>
118 struct EnableEnergy<TypeTag, TTag::EclFlowProblem> {
119  static constexpr bool value = false;
120 };
121 template<class TypeTag>
122 struct EnableFoam<TypeTag, TTag::EclFlowProblem> {
123  static constexpr bool value = false;
124 };
125 template<class TypeTag>
126 struct EnableBrine<TypeTag, TTag::EclFlowProblem> {
127  static constexpr bool value = false;
128 };
129 template<class TypeTag>
130 struct EnableMICP<TypeTag, TTag::EclFlowProblem> {
131  static constexpr bool value = false;
132 };
133 
134 template<class TypeTag>
135 struct EclWellModel<TypeTag, TTag::EclFlowProblem> {
137 };
138 template<class TypeTag>
139 struct LinearSolverSplice<TypeTag, TTag::EclFlowProblem> {
140  using type = TTag::FlowIstlSolver;
141 };
142 
143 } // namespace Opm::Properties
144 
145 namespace Opm {
152  template <class TypeTag>
154  {
155  public:
156  // --------- Types and enums ---------
158 
159  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
160  using Grid = GetPropType<TypeTag, Properties::Grid>;
161  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
162  using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
163  using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
164  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
165  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
166  using Indices = GetPropType<TypeTag, Properties::Indices>;
167  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
168  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
169 
170  typedef double Scalar;
171  static const int numEq = Indices::numEq;
172  static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
173  static const int contiZfracEqIdx = Indices::contiZfracEqIdx;
174  static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
175  static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
176  static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
177  static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
178  static const int contiBrineEqIdx = Indices::contiBrineEqIdx;
179  static const int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
180  static const int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
181  static const int contiUreaEqIdx = Indices::contiUreaEqIdx;
182  static const int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
183  static const int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
184  static const int solventSaturationIdx = Indices::solventSaturationIdx;
185  static const int zFractionIdx = Indices::zFractionIdx;
186  static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
187  static const int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
188  static const int temperatureIdx = Indices::temperatureIdx;
189  static const int foamConcentrationIdx = Indices::foamConcentrationIdx;
190  static const int saltConcentrationIdx = Indices::saltConcentrationIdx;
191  static const int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
192  static const int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
193  static const int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
194  static const int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
195  static const int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
196 
197  typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
198  typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
199  typedef typename SparseMatrixAdapter::IstlMatrix Mat;
200  typedef Dune::BlockVector<VectorBlockType> BVector;
201 
203  //typedef typename SolutionVector :: value_type PrimaryVariables ;
204 
205  // --------- Public methods ---------
206 
217  BlackoilModelEbos(Simulator& ebosSimulator,
218  const ModelParameters& param,
219  BlackoilWellModel<TypeTag>& well_model,
220  const bool terminal_output)
221  : ebosSimulator_(ebosSimulator)
222  , grid_(ebosSimulator_.vanguard().grid())
223  , phaseUsage_(phaseUsageFromDeck(eclState()))
224  , param_( param )
225  , well_model_ (well_model)
226  , terminal_output_ (terminal_output)
227  , current_relaxation_(1.0)
228  , dx_old_(UgGridHelpers::numCells(grid_))
229  {
230  // compute global sum of number of cells
231  global_nc_ = detail::countGlobalCells(grid_);
232  convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
233  }
234 
235  bool isParallel() const
236  { return grid_.comm().size() > 1; }
237 
238  const EclipseState& eclState() const
239  { return ebosSimulator_.vanguard().eclState(); }
240 
244  {
245  SimulatorReportSingle report;
246  Dune::Timer perfTimer;
247  perfTimer.start();
248  // update the solution variables in ebos
249  if ( timer.lastStepFailed() ) {
250  ebosSimulator_.model().updateFailed();
251  } else {
252  ebosSimulator_.model().advanceTimeLevel();
253  }
254 
255  // Set the timestep size, episode index, and non-linear iteration index
256  // for ebos explicitly. ebos needs to know the report step/episode index
257  // because of timing dependent data despite the fact that flow uses its
258  // own time stepper. (The length of the episode does not matter, though.)
259  ebosSimulator_.setTime(timer.simulationTimeElapsed());
260  ebosSimulator_.setTimeStepSize(timer.currentStepLength());
261  ebosSimulator_.model().newtonMethod().setIterationIndex(0);
262 
263  ebosSimulator_.problem().beginTimeStep();
264 
265  unsigned numDof = ebosSimulator_.model().numGridDof();
266  wasSwitched_.resize(numDof);
267  std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
268 
269  if (param_.update_equations_scaling_) {
270  std::cout << "equation scaling not supported yet" << std::endl;
271  //updateEquationsScaling();
272  }
273  report.pre_post_time += perfTimer.stop();
274 
275  return report;
276  }
277 
278 
288  template <class NonlinearSolverType>
290  const SimulatorTimerInterface& timer,
291  NonlinearSolverType& nonlinear_solver)
292  {
293  SimulatorReportSingle report;
294  failureReport_ = SimulatorReportSingle();
295  Dune::Timer perfTimer;
296 
297  perfTimer.start();
298  if (iteration == 0) {
299  // For each iteration we store in a vector the norms of the residual of
300  // the mass balance for each active phase, the well flux and the well equations.
301  residual_norms_history_.clear();
302  current_relaxation_ = 1.0;
303  dx_old_ = 0.0;
304  convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
305  convergence_reports_.back().report.reserve(11);
306  }
307 
308  report.total_linearizations = 1;
309 
310  try {
311  report += assembleReservoir(timer, iteration);
312  report.assemble_time += perfTimer.stop();
313  }
314  catch (...) {
315  report.assemble_time += perfTimer.stop();
316  failureReport_ += report;
317  // todo (?): make the report an attribute of the class
318  throw; // continue throwing the stick
319  }
320 
321  std::vector<double> residual_norms;
322  perfTimer.reset();
323  perfTimer.start();
324  // the step is not considered converged until at least minIter iterations is done
325  {
326  auto convrep = getConvergence(timer, iteration,residual_norms);
327  report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();;
328  ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
329  convergence_reports_.back().report.push_back(std::move(convrep));
330 
331  // Throw if any NaN or too large residual found.
332  if (severity == ConvergenceReport::Severity::NotANumber) {
333  OPM_THROW(NumericalIssue, "NaN residual found!");
334  } else if (severity == ConvergenceReport::Severity::TooLarge) {
335  OPM_THROW(NumericalIssue, "Too large residual found!");
336  }
337  }
338  report.update_time += perfTimer.stop();
339  residual_norms_history_.push_back(residual_norms);
340  if (!report.converged) {
341  perfTimer.reset();
342  perfTimer.start();
343  report.total_newton_iterations = 1;
344 
345  // enable single precision for solvers when dt is smaller then 20 days
346  //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
347 
348  // Compute the nonlinear update.
349  const int nc = UgGridHelpers::numCells(grid_);
350  BVector x(nc);
351 
352  // Solve the linear system.
353  linear_solve_setup_time_ = 0.0;
354  try {
355  // apply the Schur compliment of the well model to the reservoir linearized
356  // equations
357  // Note that linearize may throw for MSwells.
358  wellModel().linearize(ebosSimulator().model().linearizer().jacobian(),
359  ebosSimulator().model().linearizer().residual());
360 
362  report.linear_solve_setup_time += linear_solve_setup_time_;
363  report.linear_solve_time += perfTimer.stop();
364  report.total_linear_iterations += linearIterationsLastSolve();
365  }
366  catch (...) {
367  report.linear_solve_setup_time += linear_solve_setup_time_;
368  report.linear_solve_time += perfTimer.stop();
369  report.total_linear_iterations += linearIterationsLastSolve();
370 
371  failureReport_ += report;
372  throw; // re-throw up
373  }
374 
375  perfTimer.reset();
376  perfTimer.start();
377 
378  // handling well state update before oscillation treatment is a decision based
379  // on observation to avoid some big performance degeneration under some circumstances.
380  // there is no theorectical explanation which way is better for sure.
381  wellModel().postSolve(x);
382 
383  if (param_.use_update_stabilization_) {
384  // Stabilize the nonlinear update.
385  bool isOscillate = false;
386  bool isStagnate = false;
387  nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
388  if (isOscillate) {
389  current_relaxation_ -= nonlinear_solver.relaxIncrement();
390  current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
391  if (terminalOutputEnabled()) {
392  std::string msg = " Oscillating behavior detected: Relaxation set to "
393  + std::to_string(current_relaxation_);
394  OpmLog::info(msg);
395  }
396  }
397  nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
398  }
399 
400  // Apply the update, with considering model-dependent limitations and
401  // chopping of the update.
402  updateSolution(x);
403 
404  report.update_time += perfTimer.stop();
405  }
406 
407  return report;
408  }
409 
410  void printIf(int c, double x, double y, double eps, std::string type) {
411  if (std::abs(x-y) > eps) {
412  std::cout << type << " " <<c << ": "<<x << " " << y << std::endl;
413  }
414  }
415 
416 
421  {
422  SimulatorReportSingle report;
423  Dune::Timer perfTimer;
424  perfTimer.start();
425  ebosSimulator_.problem().endTimeStep();
426  report.pre_post_time += perfTimer.stop();
427  return report;
428  }
429 
435  const int iterationIdx)
436  {
437  // -------- Mass balance equations --------
438  ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
439  ebosSimulator_.problem().beginIteration();
440  ebosSimulator_.model().linearizer().linearizeDomain();
441  ebosSimulator_.problem().endIteration();
442 
443  return wellModel().lastReport();
444  }
445 
446  // compute the "relative" change of the solution between time steps
447  double relativeChange() const
448  {
449  Scalar resultDelta = 0.0;
450  Scalar resultDenom = 0.0;
451 
452  const auto& elemMapper = ebosSimulator_.model().elementMapper();
453  const auto& gridView = ebosSimulator_.gridView();
454  auto elemIt = gridView.template begin</*codim=*/0>();
455  const auto& elemEndIt = gridView.template end</*codim=*/0>();
456  for (; elemIt != elemEndIt; ++elemIt) {
457  const auto& elem = *elemIt;
458  if (elem.partitionType() != Dune::InteriorEntity)
459  continue;
460 
461  unsigned globalElemIdx = elemMapper.index(elem);
462  const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
463 
464  Scalar pressureNew;
465  pressureNew = priVarsNew[Indices::pressureSwitchIdx];
466 
467  Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
468  Scalar oilSaturationNew = 1.0;
469  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::numActivePhases() > 1) {
470  saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
471  oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
472  }
473 
474  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
475  FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
476  priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
477  assert(Indices::compositionSwitchIdx >= 0 );
478  saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
479  oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
480  }
481 
482  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
483  saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
484  }
485 
486  const auto& priVarsOld = ebosSimulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
487 
488  Scalar pressureOld;
489  pressureOld = priVarsOld[Indices::pressureSwitchIdx];
490 
491  Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
492  Scalar oilSaturationOld = 1.0;
493 
494  // NB fix me! adding pressures changes to satutation changes does not make sense
495  Scalar tmp = pressureNew - pressureOld;
496  resultDelta += tmp*tmp;
497  resultDenom += pressureNew*pressureNew;
498 
499  if (FluidSystem::numActivePhases() > 1) {
500  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
501  saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
502  oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
503  }
504 
505  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
506  FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
507  priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
508  {
509  assert(Indices::compositionSwitchIdx >= 0 );
510  saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
511  oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
512  }
513 
514  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
515  saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
516  }
517  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
518  Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
519  resultDelta += tmpSat*tmpSat;
520  resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
521  assert(std::isfinite(resultDelta));
522  assert(std::isfinite(resultDenom));
523  }
524  }
525  }
526 
527  resultDelta = gridView.comm().sum(resultDelta);
528  resultDenom = gridView.comm().sum(resultDenom);
529 
530  if (resultDenom > 0.0)
531  return resultDelta/resultDenom;
532  return 0.0;
533  }
534 
535 
538  {
539  return ebosSimulator_.model().newtonMethod().linearSolver().iterations ();
540  }
541 
544  void solveJacobianSystem(BVector& x)
545  {
546 
547  auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
548  auto& ebosResid = ebosSimulator_.model().linearizer().residual();
549 
550  // set initial guess
551  x = 0.0;
552 
553  auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
554  Dune::Timer perfTimer;
555  perfTimer.start();
556  ebosSolver.prepare(ebosJac, ebosResid);
557  linear_solve_setup_time_ = perfTimer.stop();
558  ebosSolver.setResidual(ebosResid);
559  // actually, the error needs to be calculated after setResidual in order to
560  // account for parallelization properly. since the residual of ECFV
561  // discretizations does not need to be synchronized across processes to be
562  // consistent, this is not relevant for OPM-flow...
563  ebosSolver.setMatrix(ebosJac);
564  ebosSolver.solve(x);
565  }
566 
567 
568 
570  void updateSolution(const BVector& dx)
571  {
572  auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod();
573  SolutionVector& solution = ebosSimulator_.model().solution(/*timeIdx=*/0);
574 
575  ebosNewtonMethod.update_(/*nextSolution=*/solution,
576  /*curSolution=*/solution,
577  /*update=*/dx,
578  /*resid=*/dx); // the update routines of the black
579  // oil model do not care about the
580  // residual
581 
582  // if the solution is updated, the intensive quantities need to be recalculated
583  ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
584  }
585 
588  {
589  return terminal_output_;
590  }
591 
592  template <class CollectiveCommunication>
593  double convergenceReduction(const CollectiveCommunication& comm,
594  const double pvSumLocal,
595  std::vector< Scalar >& R_sum,
596  std::vector< Scalar >& maxCoeff,
597  std::vector< Scalar >& B_avg)
598  {
599  // Compute total pore volume (use only owned entries)
600  double pvSum = pvSumLocal;
601 
602  if( comm.size() > 1 )
603  {
604  // global reduction
605  std::vector< Scalar > sumBuffer;
606  std::vector< Scalar > maxBuffer;
607  const int numComp = B_avg.size();
608  sumBuffer.reserve( 2*numComp + 1 ); // +1 for pvSum
609  maxBuffer.reserve( numComp );
610  for( int compIdx = 0; compIdx < numComp; ++compIdx )
611  {
612  sumBuffer.push_back( B_avg[ compIdx ] );
613  sumBuffer.push_back( R_sum[ compIdx ] );
614  maxBuffer.push_back( maxCoeff[ compIdx ] );
615  }
616 
617  // Compute total pore volume
618  sumBuffer.push_back( pvSum );
619 
620  // compute global sum
621  comm.sum( sumBuffer.data(), sumBuffer.size() );
622 
623  // compute global max
624  comm.max( maxBuffer.data(), maxBuffer.size() );
625 
626  // restore values to local variables
627  for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
628  {
629  B_avg[ compIdx ] = sumBuffer[ buffIdx ];
630  ++buffIdx;
631 
632  R_sum[ compIdx ] = sumBuffer[ buffIdx ];
633  }
634 
635  for( int compIdx = 0; compIdx < numComp; ++compIdx )
636  {
637  maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
638  }
639 
640  // restore global pore volume
641  pvSum = sumBuffer.back();
642  }
643 
644  // return global pore volume
645  return pvSum;
646  }
647 
648  // Get reservoir quantities on this process needed for convergence calculations.
649  double localConvergenceData(std::vector<Scalar>& R_sum,
650  std::vector<Scalar>& maxCoeff,
651  std::vector<Scalar>& B_avg)
652  {
653  double pvSumLocal = 0.0;
654  const auto& ebosModel = ebosSimulator_.model();
655  const auto& ebosProblem = ebosSimulator_.problem();
656 
657  const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
658 
659  ElementContext elemCtx(ebosSimulator_);
660  const auto& gridView = ebosSimulator().gridView();
661  const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
662  OPM_BEGIN_PARALLEL_TRY_CATCH();
663 
664  for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
665  elemIt != elemEndIt;
666  ++elemIt)
667  {
668  const auto& elem = *elemIt;
669  elemCtx.updatePrimaryStencil(elem);
670  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
671  const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
672  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
673  const auto& fs = intQuants.fluidState();
674 
675  const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx );
676  pvSumLocal += pvValue;
677 
678  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
679  {
680  if (!FluidSystem::phaseIsActive(phaseIdx)) {
681  continue;
682  }
683 
684  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
685 
686  B_avg[ compIdx ] += 1.0 / fs.invB(phaseIdx).value();
687  const auto R2 = ebosResid[cell_idx][compIdx];
688 
689  R_sum[ compIdx ] += R2;
690  maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue );
691  }
692 
693  if constexpr (has_solvent_) {
694  B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
695  const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
696  R_sum[ contiSolventEqIdx ] += R2;
697  maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue );
698  }
699  if constexpr (has_extbo_) {
700  B_avg[ contiZfracEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
701  const auto R2 = ebosResid[cell_idx][contiZfracEqIdx];
702  R_sum[ contiZfracEqIdx ] += R2;
703  maxCoeff[ contiZfracEqIdx ] = std::max( maxCoeff[ contiZfracEqIdx ], std::abs( R2 ) / pvValue );
704  }
705  if constexpr (has_polymer_) {
706  B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
707  const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
708  R_sum[ contiPolymerEqIdx ] += R2;
709  maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
710  }
711  if constexpr (has_foam_) {
712  B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
713  const auto R2 = ebosResid[cell_idx][contiFoamEqIdx];
714  R_sum[ contiFoamEqIdx ] += R2;
715  maxCoeff[ contiFoamEqIdx ] = std::max( maxCoeff[ contiFoamEqIdx ], std::abs( R2 ) / pvValue );
716  }
717  if constexpr (has_brine_) {
718  B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
719  const auto R2 = ebosResid[cell_idx][contiBrineEqIdx];
720  R_sum[ contiBrineEqIdx ] += R2;
721  maxCoeff[ contiBrineEqIdx ] = std::max( maxCoeff[ contiBrineEqIdx ], std::abs( R2 ) / pvValue );
722  }
723 
724  if constexpr (has_polymermw_) {
725  static_assert(has_polymer_);
726 
727  B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
728  // the residual of the polymer molecular equation is scaled down by a 100, since molecular weight
729  // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations
730  // TODO: there should be a more general way to determine the scaling-down coefficient
731  const auto R2 = ebosResid[cell_idx][contiPolymerMWEqIdx] / 100.;
732  R_sum[contiPolymerMWEqIdx] += R2;
733  maxCoeff[contiPolymerMWEqIdx] = std::max( maxCoeff[contiPolymerMWEqIdx], std::abs( R2 ) / pvValue );
734  }
735 
736  if constexpr (has_energy_) {
737  B_avg[ contiEnergyEqIdx ] += 1.0;
738  const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx];
739  R_sum[ contiEnergyEqIdx ] += R2;
740  maxCoeff[ contiEnergyEqIdx ] = std::max( maxCoeff[ contiEnergyEqIdx ], std::abs( R2 ) / pvValue );
741  }
742 
743  if constexpr (has_micp_) {
744  B_avg[ contiMicrobialEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
745  const auto R1 = ebosResid[cell_idx][contiMicrobialEqIdx];
746  R_sum[ contiMicrobialEqIdx ] += R1;
747  maxCoeff[ contiMicrobialEqIdx ] = std::max( maxCoeff[ contiMicrobialEqIdx ], std::abs( R1 ) / pvValue );
748  B_avg[ contiOxygenEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
749  const auto R2 = ebosResid[cell_idx][contiOxygenEqIdx];
750  R_sum[ contiOxygenEqIdx ] += R2;
751  maxCoeff[ contiOxygenEqIdx ] = std::max( maxCoeff[ contiOxygenEqIdx ], std::abs( R2 ) / pvValue );
752  B_avg[ contiUreaEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
753  const auto R3 = ebosResid[cell_idx][contiUreaEqIdx];
754  R_sum[ contiUreaEqIdx ] += R3;
755  maxCoeff[ contiUreaEqIdx ] = std::max( maxCoeff[ contiUreaEqIdx ], std::abs( R3 ) / pvValue );
756  B_avg[ contiBiofilmEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
757  const auto R4 = ebosResid[cell_idx][contiBiofilmEqIdx];
758  R_sum[ contiBiofilmEqIdx ] += R4;
759  maxCoeff[ contiBiofilmEqIdx ] = std::max( maxCoeff[ contiBiofilmEqIdx ], std::abs( R4 ) / pvValue );
760  B_avg[ contiCalciteEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
761  const auto R5 = ebosResid[cell_idx][contiCalciteEqIdx];
762  R_sum[ contiCalciteEqIdx ] += R5;
763  maxCoeff[ contiCalciteEqIdx ] = std::max( maxCoeff[ contiCalciteEqIdx ], std::abs( R5 ) / pvValue );
764  }
765  }
766 
767  OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm());
768 
769  // compute local average in terms of global number of elements
770  const int bSize = B_avg.size();
771  for ( int i = 0; i<bSize; ++i )
772  {
773  B_avg[ i ] /= Scalar( global_nc_ );
774  }
775 
776  return pvSumLocal;
777  }
778 
779  double computeCnvErrorPv(const std::vector<Scalar>& B_avg, double dt)
780  {
781  double errorPV{};
782  const auto& ebosModel = ebosSimulator_.model();
783  const auto& ebosProblem = ebosSimulator_.problem();
784  const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
785  const auto& gridView = ebosSimulator().gridView();
786  ElementContext elemCtx(ebosSimulator_);
787 
788  OPM_BEGIN_PARALLEL_TRY_CATCH();
789 
790  for (const auto& elem: elements(gridView, Dune::Partitions::interiorBorder))
791  {
792  elemCtx.updatePrimaryStencil(elem);
793  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
794  const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
795  const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx );
796  const auto& cellResidual = ebosResid[cell_idx];
797  bool cnvViolated = false;
798 
799  for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx)
800  {
801  using std::abs;
802  Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
803  cnvViolated = cnvViolated || (abs(CNV) > param_.tolerance_cnv_);
804  }
805 
806  if (cnvViolated)
807  {
808  errorPV += pvValue;
809  }
810  }
811 
812  OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::ComputeCnvError() failed: ", grid_.comm());
813 
814  return grid_.comm().sum(errorPV);
815  }
816 
817  ConvergenceReport getReservoirConvergence(const double dt,
818  const int iteration,
819  std::vector<Scalar>& B_avg,
820  std::vector<Scalar>& residual_norms)
821  {
822  typedef std::vector< Scalar > Vector;
823 
824  const int numComp = numEq;
825  Vector R_sum(numComp, 0.0 );
826  Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
827  const double pvSumLocal = localConvergenceData(R_sum, maxCoeff, B_avg);
828 
829  // compute global sum and max of quantities
830  const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
831  R_sum, maxCoeff, B_avg);
832 
833  auto cnvErrorPvFraction = computeCnvErrorPv(B_avg, dt);
834  cnvErrorPvFraction /= pvSum;
835 
836  const double tol_mb = param_.tolerance_mb_;
837  // Default value of relaxed_max_pv_fraction_ is 1 and
838  // max_strict_iter_ is 8. Hence only iteration chooses
839  // whether to use relaxed or not.
840  // To activate only fraction use fraction below 1 and iter 0.
841  const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.max_strict_iter_;
842  const double tol_cnv = use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
843 
844  // Finish computation
845  std::vector<Scalar> CNV(numComp);
846  std::vector<Scalar> mass_balance_residual(numComp);
847  for ( int compIdx = 0; compIdx < numComp; ++compIdx )
848  {
849  CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
850  mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
851  residual_norms.push_back(CNV[compIdx]);
852  }
853 
854  // Setup component names, only the first time the function is run.
855  static std::vector<std::string> compNames;
856  if (compNames.empty()) {
857  compNames.resize(numComp);
858  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
859  if (!FluidSystem::phaseIsActive(phaseIdx)) {
860  continue;
861  }
862  const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
863  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
864  compNames[compIdx] = FluidSystem::componentName(canonicalCompIdx);
865  }
866  if constexpr (has_solvent_) {
867  compNames[solventSaturationIdx] = "Solvent";
868  }
869  if constexpr (has_extbo_) {
870  compNames[zFractionIdx] = "ZFraction";
871  }
872  if constexpr (has_polymer_) {
873  compNames[polymerConcentrationIdx] = "Polymer";
874  }
875  if constexpr (has_polymermw_) {
876  assert(has_polymer_);
877  compNames[polymerMoleWeightIdx] = "MolecularWeightP";
878  }
879  if constexpr (has_energy_) {
880  compNames[temperatureIdx] = "Energy";
881  }
882  if constexpr (has_foam_) {
883  compNames[foamConcentrationIdx] = "Foam";
884  }
885  if constexpr (has_brine_) {
886  compNames[saltConcentrationIdx] = "Brine";
887  }
888  if constexpr (has_micp_) {
889  compNames[microbialConcentrationIdx] = "Microbes";
890  compNames[oxygenConcentrationIdx] = "Oxygen";
891  compNames[ureaConcentrationIdx] = "Urea";
892  compNames[biofilmConcentrationIdx] = "Biofilm";
893  compNames[calciteConcentrationIdx] = "Calcite";
894  }
895  }
896 
897  // Create convergence report.
898  ConvergenceReport report;
899  using CR = ConvergenceReport;
900  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
901  double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
902  CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
903  CR::ReservoirFailure::Type::Cnv };
904  double tol[2] = { tol_mb, tol_cnv };
905  for (int ii : {0, 1}) {
906  if (std::isnan(res[ii])) {
907  report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
908  if ( terminal_output_ ) {
909  OpmLog::debug("NaN residual for " + compNames[compIdx] + " equation.");
910  }
911  } else if (res[ii] > maxResidualAllowed()) {
912  report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
913  if ( terminal_output_ ) {
914  OpmLog::debug("Too large residual for " + compNames[compIdx] + " equation.");
915  }
916  } else if (res[ii] < 0.0) {
917  report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
918  if ( terminal_output_ ) {
919  OpmLog::debug("Negative residual for " + compNames[compIdx] + " equation.");
920  }
921  } else if (res[ii] > tol[ii]) {
922  report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
923  }
924  }
925  }
926 
927  // Output of residuals.
928  if ( terminal_output_ )
929  {
930  // Only rank 0 does print to std::cout
931  if (iteration == 0) {
932  std::string msg = "Iter";
933  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
934  msg += " MB(";
935  msg += compNames[compIdx][0];
936  msg += ") ";
937  }
938  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
939  msg += " CNV(";
940  msg += compNames[compIdx][0];
941  msg += ") ";
942  }
943  OpmLog::debug(msg);
944  }
945  std::ostringstream ss;
946  const std::streamsize oprec = ss.precision(3);
947  const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
948  ss << std::setw(4) << iteration;
949  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
950  ss << std::setw(11) << mass_balance_residual[compIdx];
951  }
952  for (int compIdx = 0; compIdx < numComp; ++compIdx) {
953  ss << std::setw(11) << CNV[compIdx];
954  }
955  ss.precision(oprec);
956  ss.flags(oflags);
957  OpmLog::debug(ss.str());
958  }
959 
960  return report;
961  }
962 
969  const int iteration,
970  std::vector<double>& residual_norms)
971  {
972  // Get convergence reports for reservoir and wells.
973  std::vector<Scalar> B_avg(numEq, 0.0);
974  auto report = getReservoirConvergence(timer.currentStepLength(), iteration, B_avg, residual_norms);
975  report += wellModel().getWellConvergence(B_avg);
976 
977  return report;
978  }
979 
980 
982  int numPhases() const
983  {
984  return phaseUsage_.num_phases;
985  }
986 
988  template<class T>
989  std::vector<std::vector<double> >
990  computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
991  {
992  return computeFluidInPlace(fipnum);
993  }
994 
996  std::vector<std::vector<double> >
997  computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
998  {
999  //assert(true)
1000  //return an empty vector
1001  std::vector<std::vector<double> > regionValues(0, std::vector<double>(0,0.0));
1002  return regionValues;
1003  }
1004 
1005  const Simulator& ebosSimulator() const
1006  { return ebosSimulator_; }
1007 
1008  Simulator& ebosSimulator()
1009  { return ebosSimulator_; }
1010 
1013  { return failureReport_; }
1014 
1015  struct StepReport
1016  {
1017  int report_step;
1018  int current_step;
1019  std::vector<ConvergenceReport> report;
1020  };
1021 
1022  const std::vector<StepReport>& stepReports() const
1023  {
1024  return convergence_reports_;
1025  }
1026 
1027  protected:
1028  // --------- Data members ---------
1029 
1030  Simulator& ebosSimulator_;
1031  const Grid& grid_;
1032  const PhaseUsage phaseUsage_;
1033  static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1034  static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1035  static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1036  static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1037  static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1038  static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1039  static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1040  static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1041 
1042  ModelParameters param_;
1043  SimulatorReportSingle failureReport_;
1044 
1045  // Well Model
1046  BlackoilWellModel<TypeTag>& well_model_;
1047 
1051  long int global_nc_;
1052 
1053  std::vector<std::vector<double>> residual_norms_history_;
1054  double current_relaxation_;
1055  BVector dx_old_;
1056 
1057  std::vector<StepReport> convergence_reports_;
1058  public:
1061  wellModel() { return well_model_; }
1062 
1064  wellModel() const { return well_model_; }
1065 
1066  void beginReportStep()
1067  {
1068  ebosSimulator_.problem().beginEpisode();
1069  }
1070 
1071  void endReportStep()
1072  {
1073  ebosSimulator_.problem().endEpisode();
1074  }
1075 
1076  private:
1077 
1078  double dpMaxRel() const { return param_.dp_max_rel_; }
1079  double dsMax() const { return param_.ds_max_; }
1080  double drMaxRel() const { return param_.dr_max_rel_; }
1081  double maxResidualAllowed() const { return param_.max_residual_allowed_; }
1082  double linear_solve_setup_time_;
1083  public:
1084  std::vector<bool> wasSwitched_;
1085  };
1086 } // namespace Opm
1087 
1088 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
Class for handling the blackoil well model.
Definition: BlackoilAquiferModel.hpp:81
A model implementation for three-phase black oil.
Definition: BlackoilModelEbos.hpp:154
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition: BlackoilModelEbos.hpp:537
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModelEbos.hpp:587
std::vector< std::vector< double > > computeFluidInPlace(const std::vector< int > &) const
Should not be called.
Definition: BlackoilModelEbos.hpp:997
ConvergenceReport getConvergence(const SimulatorTimerInterface &timer, const int iteration, std::vector< double > &residual_norms)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv).
Definition: BlackoilModelEbos.hpp:968
std::vector< std::vector< double > > computeFluidInPlace(const T &, const std::vector< int > &fipnum) const
Wrapper required due to not following generic API.
Definition: BlackoilModelEbos.hpp:990
BlackoilModelEbos(Simulator &ebosSimulator, const ModelParameters &param, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
Construct the model.
Definition: BlackoilModelEbos.hpp:217
bool terminal_output_
Whether we print something to std::cout.
Definition: BlackoilModelEbos.hpp:1049
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Called once after each time step.
Definition: BlackoilModelEbos.hpp:420
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Called once per nonlinear iteration.
Definition: BlackoilModelEbos.hpp:289
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilModelEbos.hpp:434
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Called once before each time step.
Definition: BlackoilModelEbos.hpp:243
long int global_nc_
The number of cells of the global grid.
Definition: BlackoilModelEbos.hpp:1051
const SimulatorReportSingle & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition: BlackoilModelEbos.hpp:1012
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelEbos.hpp:982
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition: BlackoilModelEbos.hpp:1061
void solveJacobianSystem(BVector &x)
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilModelEbos.hpp:544
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition: BlackoilModelEbos.hpp:570
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:94
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition: ConvergenceReport.hpp:36
This class solves the fully implicit black-oil system by solving the reduced system (after eliminatin...
Definition: ISTLSolverEbos.hpp:78
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:34
virtual int reportStepNum() const
Current report step number. This might differ from currentStepNum in case of sub stepping.
Definition: SimulatorTimerInterface.hpp:50
virtual bool lastStepFailed() const =0
Return true if last time step failed.
virtual double currentStepLength() const =0
Current step length.
virtual double simulationTimeElapsed() const =0
Time elapsed since the start of the simulation until the beginning of the current time step [s].
virtual int currentStepNum() const =0
Current step number.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition: phaseUsageFromDeck.cpp:141
Definition: BlackoilModelEbos.hpp:1016
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParametersEbos.hpp:327
double tolerance_cnv_relaxed_
Relaxed local convergence tolerance (used when iter >= max_strict_iter_).
Definition: BlackoilModelParametersEbos.hpp:346
int max_strict_iter_
Maximum number of Newton iterations before we give up on the CNV convergence criterion.
Definition: BlackoilModelParametersEbos.hpp:392
bool use_update_stabilization_
Try to detect oscillation or stagnation.
Definition: BlackoilModelParametersEbos.hpp:401
double tolerance_cnv_
Local convergence tolerance (max of local saturation errors).
Definition: BlackoilModelParametersEbos.hpp:344
bool update_equations_scaling_
Update scaling factors for mass balance equations.
Definition: BlackoilModelParametersEbos.hpp:398
double tolerance_mb_
Relative mass balance tolerance (total mass balance error).
Definition: BlackoilModelParametersEbos.hpp:342
Definition: BlackoilPhases.hpp:45
Definition: ISTLSolverEbos.hpp:51
Definition: BlackoilModelEbos.hpp:75
Definition: ISTLSolverEbos.hpp:45
Definition: BlackoilModelParametersEbos.hpp:31
Definition: NonlinearSolverEbos.hpp:41
Definition: AdaptiveTimeSteppingEbos.hpp:24
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:31