24 #ifndef OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
25 #define OPM_BLACKOILMODELEBOS_HEADER_INCLUDED
27 #include <ebos/eclproblem.hh>
28 #include <opm/models/utils/start.hh>
30 #include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
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>
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>
53 #include <opm/simulators/linalg/ISTLSolverEbos.hpp>
55 #include <dune/istl/owneroverlapcopy.hh>
56 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7)
57 #include <dune/common/parallel/communication.hh>
59 #include <dune/common/parallel/collectivecommunication.hh>
61 #include <dune/common/timer.hh>
62 #include <dune/common/unused.hh>
72 namespace Opm::Properties {
80 template<
class TypeTag>
81 struct OutputDir<TypeTag, TTag::EclFlowProblem> {
82 static constexpr
auto value =
"";
84 template<
class TypeTag>
85 struct EnableDebuggingChecks<TypeTag, TTag::EclFlowProblem> {
86 static constexpr
bool value =
false;
89 template<
class TypeTag>
90 struct BlackoilConserveSurfaceVolume<TypeTag, TTag::EclFlowProblem> {
91 static constexpr
bool value =
true;
93 template<
class TypeTag>
94 struct UseVolumetricResidual<TypeTag, TTag::EclFlowProblem> {
95 static constexpr
bool value =
false;
98 template<
class TypeTag>
99 struct EclAquiferModel<TypeTag, TTag::EclFlowProblem> {
105 template<
class TypeTag>
106 struct EnablePolymer<TypeTag, TTag::EclFlowProblem> {
107 static constexpr
bool value =
false;
109 template<
class TypeTag>
110 struct EnableSolvent<TypeTag, TTag::EclFlowProblem> {
111 static constexpr
bool value =
false;
113 template<
class TypeTag>
114 struct EnableTemperature<TypeTag, TTag::EclFlowProblem> {
115 static constexpr
bool value =
true;
117 template<
class TypeTag>
118 struct EnableEnergy<TypeTag, TTag::EclFlowProblem> {
119 static constexpr
bool value =
false;
121 template<
class TypeTag>
122 struct EnableFoam<TypeTag, TTag::EclFlowProblem> {
123 static constexpr
bool value =
false;
125 template<
class TypeTag>
126 struct EnableBrine<TypeTag, TTag::EclFlowProblem> {
127 static constexpr
bool value =
false;
129 template<
class TypeTag>
130 struct EnableSaltPrecipitation<TypeTag, TTag::EclFlowProblem> {
131 static constexpr
bool value =
false;
133 template<
class TypeTag>
134 struct EnableMICP<TypeTag, TTag::EclFlowProblem> {
135 static constexpr
bool value =
false;
138 template<
class TypeTag>
142 template<
class TypeTag>
143 struct LinearSolverSplice<TypeTag, TTag::EclFlowProblem> {
156 template <
class TypeTag>
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>;
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;
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;
224 const bool terminal_output)
225 : ebosSimulator_(ebosSimulator)
226 , grid_(ebosSimulator_.vanguard().grid())
229 , well_model_ (well_model)
231 , current_relaxation_(1.0)
232 , dx_old_(UgGridHelpers::numCells(grid_))
236 convergence_reports_.reserve(300);
239 bool isParallel()
const
240 {
return grid_.comm().size() > 1; }
242 const EclipseState& eclState()
const
243 {
return ebosSimulator_.vanguard().eclState(); }
250 Dune::Timer perfTimer;
254 ebosSimulator_.model().updateFailed();
256 ebosSimulator_.model().advanceTimeLevel();
265 ebosSimulator_.model().newtonMethod().setIterationIndex(0);
267 ebosSimulator_.problem().beginTimeStep();
269 unsigned numDof = ebosSimulator_.model().numGridDof();
270 wasSwitched_.resize(numDof);
271 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
274 std::cout <<
"equation scaling not supported yet" << std::endl;
277 report.pre_post_time += perfTimer.stop();
292 template <
class NonlinearSolverType>
295 NonlinearSolverType& nonlinear_solver)
299 Dune::Timer perfTimer;
302 if (iteration == 0) {
305 residual_norms_history_.clear();
306 current_relaxation_ = 1.0;
309 convergence_reports_.back().report.reserve(11);
312 report.total_linearizations = 1;
316 report.assemble_time += perfTimer.stop();
319 report.assemble_time += perfTimer.stop();
320 failureReport_ += report;
325 std::vector<double> 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));
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!");
342 report.update_time += perfTimer.stop();
343 residual_norms_history_.push_back(residual_norms);
344 if (!report.converged) {
347 report.total_newton_iterations = 1;
353 const int nc = UgGridHelpers::numCells(grid_);
357 linear_solve_setup_time_ = 0.0;
362 wellModel().linearize(ebosSimulator().model().linearizer().jacobian(),
363 ebosSimulator().model().linearizer().residual());
366 report.linear_solve_setup_time += linear_solve_setup_time_;
367 report.linear_solve_time += perfTimer.stop();
371 report.linear_solve_setup_time += linear_solve_setup_time_;
372 report.linear_solve_time += perfTimer.stop();
375 failureReport_ += report;
389 bool isOscillate =
false;
390 bool isStagnate =
false;
391 nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
393 current_relaxation_ -= nonlinear_solver.relaxIncrement();
394 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
396 std::string msg =
" Oscillating behavior detected: Relaxation set to "
397 + std::to_string(current_relaxation_);
401 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
408 report.update_time += perfTimer.stop();
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;
427 Dune::Timer perfTimer;
429 ebosSimulator_.problem().endTimeStep();
430 report.pre_post_time += perfTimer.stop();
439 const int iterationIdx)
442 ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
443 ebosSimulator_.problem().beginIteration();
444 ebosSimulator_.model().linearizer().linearizeDomain();
445 ebosSimulator_.problem().endIteration();
451 double relativeChange()
const
453 Scalar resultDelta = 0.0;
454 Scalar resultDenom = 0.0;
456 const auto& elemMapper = ebosSimulator_.model().elementMapper();
457 const auto& gridView = ebosSimulator_.gridView();
458 auto elemIt = gridView.template begin<0>();
459 const auto& elemEndIt = gridView.template end<0>();
460 for (; elemIt != elemEndIt; ++elemIt) {
461 const auto& elem = *elemIt;
462 if (elem.partitionType() != Dune::InteriorEntity)
465 unsigned globalElemIdx = elemMapper.index(elem);
466 const auto& priVarsNew = ebosSimulator_.model().solution(0)[globalElemIdx];
469 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
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];
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];
486 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
487 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
490 const auto& priVarsOld = ebosSimulator_.model().solution(1)[globalElemIdx];
493 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
495 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
496 Scalar oilSaturationOld = 1.0;
499 Scalar tmp = pressureNew - pressureOld;
500 resultDelta += tmp*tmp;
501 resultDenom += pressureNew*pressureNew;
503 if (FluidSystem::numActivePhases() > 1) {
504 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
505 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
506 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
509 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
510 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
511 priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
513 assert(Indices::compositionSwitchIdx >= 0 );
514 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
515 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
518 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
519 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
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));
531 resultDelta = gridView.comm().sum(resultDelta);
532 resultDenom = gridView.comm().sum(resultDenom);
534 if (resultDenom > 0.0)
535 return resultDelta/resultDenom;
543 return ebosSimulator_.model().newtonMethod().linearSolver().iterations ();
551 auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
552 auto& ebosResid = ebosSimulator_.model().linearizer().residual();
557 auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
558 Dune::Timer perfTimer;
560 ebosSolver.prepare(ebosJac, ebosResid);
561 linear_solve_setup_time_ = perfTimer.stop();
562 ebosSolver.setResidual(ebosResid);
567 ebosSolver.setMatrix(ebosJac);
576 auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod();
577 SolutionVector& solution = ebosSimulator_.model().solution(0);
579 ebosNewtonMethod.update_(solution,
587 ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(0);
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)
605 double pvSum = pvSumLocal;
606 double numAquiferPvSum = numAquiferPvSumLocal;
608 if( comm.size() > 1 )
611 std::vector< Scalar > sumBuffer;
612 std::vector< Scalar > maxBuffer;
613 const int numComp = B_avg.size();
614 sumBuffer.reserve( 2*numComp + 2 );
615 maxBuffer.reserve( numComp );
616 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
618 sumBuffer.push_back( B_avg[ compIdx ] );
619 sumBuffer.push_back( R_sum[ compIdx ] );
620 maxBuffer.push_back( maxCoeff[ compIdx ] );
624 sumBuffer.push_back( pvSum );
625 sumBuffer.push_back( numAquiferPvSum );
628 comm.sum( sumBuffer.data(), sumBuffer.size() );
631 comm.max( maxBuffer.data(), maxBuffer.size() );
634 for(
int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
636 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
639 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
642 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
644 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
648 pvSum = sumBuffer[sumBuffer.size()-2];
649 numAquiferPvSum = sumBuffer.back();
653 return {pvSum, numAquiferPvSum};
660 std::vector<Scalar>& maxCoeff,
661 std::vector<Scalar>& B_avg)
663 double pvSumLocal = 0.0;
664 double numAquiferPvSumLocal = 0.0;
665 const auto& ebosModel = ebosSimulator_.model();
666 const auto& ebosProblem = ebosSimulator_.problem();
668 const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
670 ElementContext elemCtx(ebosSimulator_);
671 const auto& gridView = ebosSimulator().gridView();
672 const auto& elemEndIt = gridView.template end<0, Dune::Interior_Partition>();
673 OPM_BEGIN_PARALLEL_TRY_CATCH();
675 for (
auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
679 const auto& elem = *elemIt;
680 elemCtx.updatePrimaryStencil(elem);
681 elemCtx.updatePrimaryIntensiveQuantities(0);
682 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
683 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
684 const auto& fs = intQuants.fluidState();
686 const double pvValue = ebosProblem.referencePorosity(cell_idx, 0) * ebosModel.dofTotalVolume( cell_idx );
687 pvSumLocal += pvValue;
689 if (isNumericalAquiferCell(gridView.grid(), elem))
691 numAquiferPvSumLocal += pvValue;
694 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
696 if (!FluidSystem::phaseIsActive(phaseIdx)) {
700 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
702 B_avg[ compIdx ] += 1.0 / fs.invB(phaseIdx).value();
703 const auto R2 = ebosResid[cell_idx][compIdx];
705 R_sum[ compIdx ] += R2;
706 maxCoeff[ compIdx ] = std::max( maxCoeff[ compIdx ], std::abs( R2 ) / pvValue );
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 );
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 );
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 );
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 );
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 );
740 if constexpr (has_polymermw_) {
741 static_assert(has_polymer_);
743 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
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 );
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 );
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 );
783 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm());
786 const int bSize = B_avg.size();
787 for (
int i = 0; i<bSize; ++i )
792 return {pvSumLocal, numAquiferPvSumLocal};
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_);
806 OPM_BEGIN_PARALLEL_TRY_CATCH();
808 for (
const auto& elem: elements(gridView, Dune::Partitions::interiorBorder))
811 if (isNumericalAquiferCell(gridView.grid(), elem))
815 elemCtx.updatePrimaryStencil(elem);
816 elemCtx.updatePrimaryIntensiveQuantities(0);
817 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
818 const double pvValue = ebosProblem.referencePorosity(cell_idx, 0) * ebosModel.dofTotalVolume( cell_idx );
819 const auto& cellResidual = ebosResid[cell_idx];
820 bool cnvViolated =
false;
822 for (
unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx)
825 Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
835 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModelEbos::ComputeCnvError() failed: ", grid_.comm());
837 return grid_.comm().sum(errorPV);
842 std::vector<Scalar>& B_avg,
843 std::vector<Scalar>& residual_norms)
845 typedef std::vector< Scalar > Vector;
847 const int numComp = numEq;
848 Vector R_sum(numComp, 0.0 );
849 Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
853 const auto [ pvSum, numAquiferPvSum ] =
854 convergenceReduction(grid_.comm(), pvSumLocal,
855 numAquiferPvSumLocal,
856 R_sum, maxCoeff, B_avg);
859 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
866 const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.
max_strict_iter_;
870 std::vector<Scalar> CNV(numComp);
871 std::vector<Scalar> mass_balance_residual(numComp);
872 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
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]);
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)) {
887 const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
888 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
889 compNames[compIdx] = FluidSystem::componentName(canonicalCompIdx);
891 if constexpr (has_solvent_) {
892 compNames[solventSaturationIdx] =
"Solvent";
894 if constexpr (has_extbo_) {
895 compNames[zFractionIdx] =
"ZFraction";
897 if constexpr (has_polymer_) {
898 compNames[polymerConcentrationIdx] =
"Polymer";
900 if constexpr (has_polymermw_) {
901 assert(has_polymer_);
902 compNames[polymerMoleWeightIdx] =
"MolecularWeightP";
904 if constexpr (has_energy_) {
905 compNames[temperatureIdx] =
"Energy";
907 if constexpr (has_foam_) {
908 compNames[foamConcentrationIdx] =
"Foam";
910 if constexpr (has_brine_) {
911 compNames[saltConcentrationIdx] =
"Brine";
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";
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});
934 OpmLog::debug(
"NaN residual for " + compNames[compIdx] +
" equation.");
936 }
else if (res[ii] > maxResidualAllowed()) {
937 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
939 OpmLog::debug(
"Too large residual for " + compNames[compIdx] +
" equation.");
941 }
else if (res[ii] < 0.0) {
942 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
944 OpmLog::debug(
"Negative residual for " + compNames[compIdx] +
" equation.");
946 }
else if (res[ii] > tol[ii]) {
947 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
956 if (iteration == 0) {
957 std::string msg =
"Iter";
958 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
960 msg += compNames[compIdx][0];
963 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
965 msg += compNames[compIdx][0];
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];
977 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
978 ss << std::setw(11) << CNV[compIdx];
982 OpmLog::debug(ss.str());
995 std::vector<double>& residual_norms)
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);
1009 return phaseUsage_.num_phases;
1014 std::vector<std::vector<double> >
1021 std::vector<std::vector<double> >
1026 std::vector<std::vector<double> > regionValues(0, std::vector<double>(0,0.0));
1027 return regionValues;
1030 const Simulator& ebosSimulator()
const
1031 {
return ebosSimulator_; }
1033 Simulator& ebosSimulator()
1034 {
return ebosSimulator_; }
1038 {
return failureReport_; }
1044 std::vector<ConvergenceReport> report;
1047 const std::vector<StepReport>& stepReports()
const
1049 return convergence_reports_;
1055 Simulator& ebosSimulator_;
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>();
1067 ModelParameters param_;
1078 std::vector<std::vector<double>> residual_norms_history_;
1079 double current_relaxation_;
1082 std::vector<StepReport> convergence_reports_;
1089 wellModel()
const {
return well_model_; }
1091 void beginReportStep()
1093 ebosSimulator_.problem().beginEpisode();
1096 void endReportStep()
1098 ebosSimulator_.problem().endEpisode();
1103 bool isNumericalAquiferCell(
const Dune::CpGrid& grid,
const T& elem)
1105 const auto& aquiferCells = grid.sortedNumAquiferCells();
1106 if (aquiferCells.empty())
1110 auto candidate = std::lower_bound(aquiferCells.begin(), aquiferCells.end(),
1112 return candidate != aquiferCells.end() && *candidate == elem.index();
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)
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_;
1128 std::vector<bool> wasSwitched_;
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 ¶m, 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