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