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