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
70namespace Opm::Properties {
71
72namespace TTag {
74 using InheritsFrom = std::tuple<FlowTimeSteppingParameters, FlowModelParameters,
75 FlowNonLinearSolver, EclBaseProblem, BlackOilModel>;
76};
77}
78template<class TypeTag>
79struct OutputDir<TypeTag, TTag::EclFlowProblem> {
80 static constexpr auto value = "";
81};
82template<class TypeTag>
83struct EnableDebuggingChecks<TypeTag, TTag::EclFlowProblem> {
84 static constexpr bool value = false;
85};
86// default in flow is to formulate the equations in surface volumes
87template<class TypeTag>
88struct BlackoilConserveSurfaceVolume<TypeTag, TTag::EclFlowProblem> {
89 static constexpr bool value = true;
90};
91template<class TypeTag>
92struct UseVolumetricResidual<TypeTag, TTag::EclFlowProblem> {
93 static constexpr bool value = false;
94};
95
96template<class TypeTag>
97struct 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
103template<class TypeTag>
104struct EnablePolymer<TypeTag, TTag::EclFlowProblem> {
105 static constexpr bool value = false;
106};
107template<class TypeTag>
108struct EnableSolvent<TypeTag, TTag::EclFlowProblem> {
109 static constexpr bool value = false;
110};
111template<class TypeTag>
112struct EnableTemperature<TypeTag, TTag::EclFlowProblem> {
113 static constexpr bool value = true;
114};
115template<class TypeTag>
116struct EnableEnergy<TypeTag, TTag::EclFlowProblem> {
117 static constexpr bool value = false;
118};
119template<class TypeTag>
120struct EnableFoam<TypeTag, TTag::EclFlowProblem> {
121 static constexpr bool value = false;
122};
123template<class TypeTag>
124struct EnableBrine<TypeTag, TTag::EclFlowProblem> {
125 static constexpr bool value = false;
126};
127template<class TypeTag>
128struct EnableSaltPrecipitation<TypeTag, TTag::EclFlowProblem> {
129 static constexpr bool value = false;
130};
131template<class TypeTag>
132struct EnableMICP<TypeTag, TTag::EclFlowProblem> {
133 static constexpr bool value = false;
134};
135
136template<class TypeTag>
137struct EclWellModel<TypeTag, TTag::EclFlowProblem> {
139};
140template<class TypeTag>
141struct LinearSolverSplice<TypeTag, TTag::EclFlowProblem> {
143};
144
145} // namespace Opm::Properties
146
147namespace 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 {
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 {
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.
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 {
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
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
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition: BlackoilModelEbos.hpp:1071
std::vector< std::vector< double > > computeFluidInPlace(const std::vector< int > &) const
Should not be called.
Definition: BlackoilModelEbos.hpp:1007
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
const SimulatorReportSingle & failureReport() const
return the statistics if the nonlinearIteration() method failed
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:978
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
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
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
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelEbos.hpp:992
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
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:40
Definition: AdaptiveTimeSteppingEbos.hpp:24
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:31