176 static constexpr int numEq = Indices::numEq;
177 static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
178 static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx;
179 static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
180 static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
181 static constexpr int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
182 static constexpr int contiFoamEqIdx = Indices::contiFoamEqIdx;
183 static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
184 static constexpr int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
185 static constexpr int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
186 static constexpr int contiUreaEqIdx = Indices::contiUreaEqIdx;
187 static constexpr int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
188 static constexpr int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
189 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
190 static constexpr int zFractionIdx = Indices::zFractionIdx;
191 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
192 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
193 static constexpr int temperatureIdx = Indices::temperatureIdx;
194 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
195 static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
196 static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
197 static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
198 static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
199 static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
200 static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
202 using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
203 using MatrixBlockType =
typename SparseMatrixAdapter::MatrixBlock;
204 using Mat =
typename SparseMatrixAdapter::IstlMatrix;
205 using BVector = Dune::BlockVector<VectorBlockType>;
225 : ebosSimulator_(ebosSimulator)
226 , grid_(ebosSimulator_.
vanguard().grid())
231 , current_relaxation_(1.0)
232 , dx_old_(ebosSimulator_.model().
numGridDof())
236 convergence_reports_.reserve(300);
238 if (param_.nonlinear_solver_ ==
"nldd") {
240 OpmLog::info(
"Using Non-Linear Domain Decomposition solver (nldd).");
242 nlddSolver_ = std::make_unique<BlackoilModelEbosNldd<TypeTag>>(*this);
243 }
else if (param_.nonlinear_solver_ ==
"newton") {
245 OpmLog::info(
"Using Newton nonlinear solver.");
248 OPM_THROW(std::runtime_error,
"Unknown nonlinear solver option: " + param_.nonlinear_solver_);
253 bool isParallel()
const
254 {
return grid_.comm().size() > 1; }
257 const EclipseState& eclState()
const
258 {
return ebosSimulator_.vanguard().eclState(); }
269 if (
timer.lastStepFailed() ) {
270 ebosSimulator_.model().updateFailed();
272 ebosSimulator_.model().advanceTimeLevel();
279 ebosSimulator_.setTime(
timer.simulationTimeElapsed());
280 ebosSimulator_.setTimeStepSize(
timer.currentStepLength());
281 ebosSimulator_.model().newtonMethod().setIterationIndex(0);
283 ebosSimulator_.problem().beginTimeStep();
285 unsigned numDof = ebosSimulator_.model().numGridDof();
286 wasSwitched_.resize(
numDof);
287 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
289 if (param_.update_equations_scaling_) {
290 OpmLog::error(
"Equation scaling not supported");
298 report.pre_post_time +=
perfTimer.stop();
314 report.total_linearizations = 1;
319 report.assemble_time +=
perfTimer.stop();
322 report.assemble_time += perfTimer.stop();
323 failureReport_ += report;
328 std::vector<double> residual_norms;
334 report.converged = convrep.converged() && iteration > minIter;
335 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
336 convergence_reports_.back().report.push_back(std::move(convrep));
339 if (severity == ConvergenceReport::Severity::NotANumber) {
340 OPM_THROW(NumericalProblem,
"NaN residual found!");
341 }
else if (severity == ConvergenceReport::Severity::TooLarge) {
342 OPM_THROW_NOLOG(NumericalProblem,
"Too large residual found!");
345 report.update_time += perfTimer.stop();
346 residual_norms_history_.push_back(residual_norms);
357 template <
class NonlinearSolverType>
365 residual_norms_history_.clear();
366 current_relaxation_ = 1.0;
368 convergence_reports_.push_back({
timer.reportStepNum(),
timer.currentStepNum(), {}});
369 convergence_reports_.back().report.reserve(11);
372 if ((this->param_.nonlinear_solver_ !=
"nldd") ||
383 template <
class NonlinearSolverType>
396 if (!report.converged) {
399 report.total_newton_iterations = 1;
402 unsigned nc = ebosSimulator_.model().numGridDof();
406 linear_solve_setup_time_ = 0.0;
412 ebosSimulator().model().
linearizer().residual());
417 report.linear_solve_setup_time += linear_solve_setup_time_;
418 report.linear_solve_time +=
perfTimer.stop();
422 report.linear_solve_setup_time += linear_solve_setup_time_;
423 report.linear_solve_time += perfTimer.stop();
426 failureReport_ += report;
438 if (param_.use_update_stabilization_) {
440 bool isOscillate =
false;
441 bool isStagnate =
false;
442 nonlinear_solver.detectOscillations(residual_norms_history_, residual_norms_history_.size() - 1, isOscillate, isStagnate);
444 current_relaxation_ -= nonlinear_solver.relaxIncrement();
445 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
447 std::string msg =
" Oscillating behavior detected: Relaxation set to "
448 + std::to_string(current_relaxation_);
452 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
460 report.update_time += perfTimer.stop();
475 ebosSimulator_.problem().endTimeStep();
476 report.pre_post_time +=
perfTimer.stop();
485 ebosSimulator_.model().newtonMethod().setIterationIndex(
iterationIdx);
486 ebosSimulator_.problem().beginIteration();
487 ebosSimulator_.model().linearizer().linearizeDomain();
488 ebosSimulator_.problem().endIteration();
493 double relativeChange()
const
498 const auto&
elemMapper = ebosSimulator_.model().elementMapper();
499 const auto&
gridView = ebosSimulator_.gridView();
508 Scalar oilSaturationNew = 1.0;
509 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
510 FluidSystem::numActivePhases() > 1 &&
511 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
512 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
513 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
516 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
517 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
518 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
519 assert(Indices::compositionSwitchIdx >= 0 );
520 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
521 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
524 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
525 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
528 const auto& priVarsOld = ebosSimulator_.model().solution(1)[globalElemIdx];
531 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
533 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
534 Scalar oilSaturationOld = 1.0;
537 Scalar tmp = pressureNew - pressureOld;
538 resultDelta += tmp*tmp;
539 resultDenom += pressureNew*pressureNew;
541 if (FluidSystem::numActivePhases() > 1) {
542 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
543 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
544 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
547 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
549 assert(Indices::compositionSwitchIdx >= 0 );
550 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
551 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
554 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
555 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
557 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
558 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
559 resultDelta += tmpSat*tmpSat;
560 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
561 assert(std::isfinite(resultDelta));
562 assert(std::isfinite(resultDenom));
567 resultDelta = gridView.comm().sum(resultDelta);
568 resultDenom = gridView.comm().sum(resultDenom);
570 if (resultDenom > 0.0)
571 return resultDelta/resultDenom;
579 return ebosSimulator_.model().newtonMethod().linearSolver().iterations ();
584 double& linearSolveSetupTime()
586 return linear_solve_setup_time_;
595 auto&
ebosJac = ebosSimulator_.model().linearizer().jacobian().istlMatrix();
596 auto&
ebosResid = ebosSimulator_.model().linearizer().residual();
597 auto&
ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
603 OpmLog::debug(
"\nRunning speed test for comparing available linear solvers.");
644 linear_solve_setup_time_ =
perfTimer.stop();
672 ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(0);
673 ebosSimulator_.problem().eclWriter()->mutableEclOutputModule().invalidateLocalData();
683 std::tuple<double,double> convergenceReduction(Parallel::Communication comm,
686 std::vector< Scalar >&
R_sum,
688 std::vector< Scalar >&
B_avg)
695 if( comm.size() > 1 )
711 sumBuffer.push_back( pvSum );
712 sumBuffer.push_back( numAquiferPvSum );
715 comm.sum( sumBuffer.data(), sumBuffer.size() );
718 comm.max( maxBuffer.data(), maxBuffer.size() );
721 for(
int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
723 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
726 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
729 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
731 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
735 pvSum = sumBuffer[sumBuffer.size()-2];
736 numAquiferPvSum = sumBuffer.back();
740 return {pvSum, numAquiferPvSum};
748 std::vector<Scalar>&
B_avg,
754 const auto&
ebosModel = ebosSimulator_.model();
755 const auto&
ebosProblem = ebosSimulator_.problem();
757 const auto&
ebosResid = ebosSimulator_.model().linearizer().residual();
760 const auto&
gridView = ebosSimulator().gridView();
762 OPM_BEGIN_PARALLEL_TRY_CATCH();
765 elemCtx.updatePrimaryIntensiveQuantities(0);
784 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm());
803 const auto&
ebosModel = ebosSimulator_.model();
804 const auto&
ebosProblem = ebosSimulator_.problem();
805 const auto&
ebosResid = ebosSimulator_.model().linearizer().residual();
806 const auto&
gridView = ebosSimulator().gridView();
810 OPM_BEGIN_PARALLEL_TRY_CATCH();
838 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModelEbos::ComputeCnvError() failed: ", grid_.comm());
840 return grid_.comm().sum(
errorPV);
845 param_.tolerance_mb_ =
tuning.XXXMBE;
847 OpmLog::debug(fmt::format(
"Setting BlackoilModelEbos mass balance limit (XXXMBE) to {:.2e}",
tuning.XXXMBE));
852 ConvergenceReport getReservoirConvergence(
const double reportTime,
855 std::vector<Scalar>& B_avg,
856 std::vector<Scalar>& residual_norms)
858 OPM_TIMEBLOCK(getReservoirConvergence);
859 using Vector = std::vector<Scalar>;
861 const int numComp = numEq;
862 Vector R_sum(numComp, 0.0 );
863 Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
864 std::vector<int> maxCoeffCell(numComp, -1);
865 const auto [ pvSumLocal, numAquiferPvSumLocal] =
localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);
868 const auto [ pvSum, numAquiferPvSum ] =
869 convergenceReduction(grid_.comm(), pvSumLocal,
870 numAquiferPvSumLocal,
871 R_sum, maxCoeff, B_avg);
874 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
876 const double tol_mb = param_.tolerance_mb_;
880 const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.min_strict_cnv_iter_;
881 const double tol_cnv = use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
884 std::vector<Scalar> CNV(numComp);
885 std::vector<Scalar> mass_balance_residual(numComp);
886 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
888 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
889 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
890 residual_norms.push_back(CNV[compIdx]);
894 ConvergenceReport report{reportTime};
895 using CR = ConvergenceReport;
896 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
897 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
898 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
899 CR::ReservoirFailure::Type::Cnv };
900 double tol[2] = { tol_mb, tol_cnv };
901 for (
int ii : {0, 1}) {
902 if (std::isnan(res[ii])) {
903 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
905 OpmLog::debug(
"NaN residual for " + this->compNames_.name(compIdx) +
" equation.");
907 }
else if (res[ii] > maxResidualAllowed()) {
908 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
910 OpmLog::debug(
"Too large residual for " + this->compNames_.name(compIdx) +
" equation.");
912 }
else if (res[ii] < 0.0) {
913 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
915 OpmLog::debug(
"Negative residual for " + this->compNames_.name(compIdx) +
" equation.");
917 }
else if (res[ii] > tol[ii]) {
918 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
920 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
928 if (iteration == 0) {
929 std::string msg =
"Iter";
930 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
932 msg += this->compNames_.name(compIdx)[0];
935 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
937 msg += this->compNames_.name(compIdx)[0];
942 std::ostringstream ss;
943 const std::streamsize oprec = ss.precision(3);
944 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
945 ss << std::setw(4) << iteration;
946 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
947 ss << std::setw(11) << mass_balance_residual[compIdx];
949 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
950 ss << std::setw(11) << CNV[compIdx];
954 OpmLog::debug(ss.str());
971 std::vector<Scalar>
B_avg(numEq, 0.0);
972 auto report = getReservoirConvergence(
timer.simulationTimeElapsed(),
973 timer.currentStepLength(),
977 report +=
wellModel().getWellConvergence(
B_avg, report.converged());
986 return phaseUsage_.num_phases;
991 std::vector<std::vector<double> >
998 std::vector<std::vector<double> >
1004 std::vector<std::vector<double> >
regionValues(0, std::vector<double>(0,0.0));
1008 const Simulator& ebosSimulator()
const
1009 {
return ebosSimulator_; }
1011 Simulator& ebosSimulator()
1012 {
return ebosSimulator_; }
1016 {
return failureReport_; }
1025 const std::vector<StepReport>& stepReports()
const
1027 return convergence_reports_;
1033 Simulator& ebosSimulator_;
1035 const PhaseUsage phaseUsage_;
1036 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1037 static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1038 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1039 static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1040 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1041 static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1042 static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1043 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1045 ModelParameters param_;
1046 SimulatorReportSingle failureReport_;
1049 BlackoilWellModel<TypeTag>& well_model_;
1056 std::vector<std::vector<double>> residual_norms_history_;
1057 double current_relaxation_;
1060 std::vector<StepReport> convergence_reports_;
1071 wellModel()
const {
return well_model_; }
1073 void beginReportStep()
1075 ebosSimulator_.problem().beginEpisode();
1078 void endReportStep()
1080 ebosSimulator_.problem().endEpisode();
1083 template<
class Flu
idState,
class Res
idual>
1084 void getMaxCoeff(
const unsigned cell_idx,
1085 const IntensiveQuantities& intQuants,
1086 const FluidState& fs,
1087 const Residual& ebosResid,
1088 const Scalar pvValue,
1089 std::vector<Scalar>& B_avg,
1090 std::vector<Scalar>& R_sum,
1091 std::vector<Scalar>& maxCoeff,
1092 std::vector<int>& maxCoeffCell)
1094 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1096 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1100 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1102 B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
1103 const auto R2 = ebosResid[cell_idx][compIdx];
1105 R_sum[compIdx] += R2;
1106 const double Rval = std::abs(R2) / pvValue;
1107 if (Rval > maxCoeff[compIdx]) {
1108 maxCoeff[compIdx] = Rval;
1109 maxCoeffCell[compIdx] = cell_idx;
1113 if constexpr (has_solvent_) {
1114 B_avg[contiSolventEqIdx] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1115 const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
1116 R_sum[contiSolventEqIdx] += R2;
1117 maxCoeff[contiSolventEqIdx] = std::max(maxCoeff[contiSolventEqIdx],
1118 std::abs(R2) / pvValue);
1120 if constexpr (has_extbo_) {
1121 B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1122 const auto R2 = ebosResid[cell_idx][contiZfracEqIdx];
1123 R_sum[ contiZfracEqIdx ] += R2;
1124 maxCoeff[contiZfracEqIdx] = std::max(maxCoeff[contiZfracEqIdx],
1125 std::abs(R2) / pvValue);
1127 if constexpr (has_polymer_) {
1128 B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1129 const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
1130 R_sum[contiPolymerEqIdx] += R2;
1131 maxCoeff[contiPolymerEqIdx] = std::max(maxCoeff[contiPolymerEqIdx],
1132 std::abs(R2) / pvValue);
1134 if constexpr (has_foam_) {
1135 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1136 const auto R2 = ebosResid[cell_idx][contiFoamEqIdx];
1137 R_sum[contiFoamEqIdx] += R2;
1138 maxCoeff[contiFoamEqIdx] = std::max(maxCoeff[contiFoamEqIdx],
1139 std::abs(R2) / pvValue);
1141 if constexpr (has_brine_) {
1142 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1143 const auto R2 = ebosResid[cell_idx][contiBrineEqIdx];
1144 R_sum[contiBrineEqIdx] += R2;
1145 maxCoeff[contiBrineEqIdx] = std::max(maxCoeff[contiBrineEqIdx],
1146 std::abs(R2) / pvValue);
1149 if constexpr (has_polymermw_) {
1150 static_assert(has_polymer_);
1152 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1156 const auto R2 = ebosResid[cell_idx][contiPolymerMWEqIdx] / 100.;
1157 R_sum[contiPolymerMWEqIdx] += R2;
1158 maxCoeff[contiPolymerMWEqIdx] = std::max(maxCoeff[contiPolymerMWEqIdx],
1159 std::abs(R2) / pvValue);
1162 if constexpr (has_energy_) {
1163 B_avg[contiEnergyEqIdx] += 1.0 / (4.182e1);
1164 const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx];
1165 R_sum[contiEnergyEqIdx] += R2;
1166 maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
1167 std::abs(R2) / pvValue);
1170 if constexpr (has_micp_) {
1171 B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1172 const auto R1 = ebosResid[cell_idx][contiMicrobialEqIdx];
1173 R_sum[contiMicrobialEqIdx] += R1;
1174 maxCoeff[contiMicrobialEqIdx] = std::max(maxCoeff[contiMicrobialEqIdx],
1175 std::abs(R1) / pvValue);
1176 B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1177 const auto R2 = ebosResid[cell_idx][contiOxygenEqIdx];
1178 R_sum[contiOxygenEqIdx] += R2;
1179 maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
1180 std::abs(R2) / pvValue);
1181 B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1182 const auto R3 = ebosResid[cell_idx][contiUreaEqIdx];
1183 R_sum[contiUreaEqIdx] += R3;
1184 maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
1185 std::abs(R3) / pvValue);
1186 B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1187 const auto R4 = ebosResid[cell_idx][contiBiofilmEqIdx];
1188 R_sum[contiBiofilmEqIdx] += R4;
1189 maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
1190 std::abs(R4) / pvValue);
1191 B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1192 const auto R5 = ebosResid[cell_idx][contiCalciteEqIdx];
1193 R_sum[contiCalciteEqIdx] += R5;
1194 maxCoeff[contiCalciteEqIdx] = std::max(maxCoeff[contiCalciteEqIdx],
1195 std::abs(R5) / pvValue);
1212 double dpMaxRel()
const {
return param_.dp_max_rel_; }
1213 double dsMax()
const {
return param_.ds_max_; }
1214 double drMaxRel()
const {
return param_.dr_max_rel_; }
1215 double maxResidualAllowed()
const {
return param_.max_residual_allowed_; }
1216 double linear_solve_setup_time_;
1219 std::vector<bool> wasSwitched_;