172 static const int numEq = Indices::numEq;
173 static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
174 static const int contiZfracEqIdx = Indices::contiZfracEqIdx;
175 static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
176 static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
177 static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
178 static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
179 static const int contiBrineEqIdx = Indices::contiBrineEqIdx;
180 static const int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
181 static const int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
182 static const int contiUreaEqIdx = Indices::contiUreaEqIdx;
183 static const int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
184 static const int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
185 static const int solventSaturationIdx = Indices::solventSaturationIdx;
186 static const int zFractionIdx = Indices::zFractionIdx;
187 static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
188 static const int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
189 static const int temperatureIdx = Indices::temperatureIdx;
190 static const int foamConcentrationIdx = Indices::foamConcentrationIdx;
191 static const int saltConcentrationIdx = Indices::saltConcentrationIdx;
192 static const int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
193 static const int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
194 static const int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
195 static const int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
196 static const int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
198 typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
199 typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
200 typedef typename SparseMatrixAdapter::IstlMatrix Mat;
201 typedef Dune::BlockVector<VectorBlockType> BVector;
212 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
221 if constexpr (has_solvent_) {
222 names_[solventSaturationIdx] =
"Solvent";
225 if constexpr (has_extbo_) {
226 names_[zFractionIdx] =
"ZFraction";
229 if constexpr (has_polymer_) {
230 names_[polymerConcentrationIdx] =
"Polymer";
233 if constexpr (has_polymermw_) {
235 names_[polymerMoleWeightIdx] =
"MolecularWeightP";
238 if constexpr (has_energy_) {
239 names_[temperatureIdx] =
"Energy";
242 if constexpr (has_foam_) {
243 names_[foamConcentrationIdx] =
"Foam";
246 if constexpr (has_brine_) {
247 names_[saltConcentrationIdx] =
"Brine";
250 if constexpr (has_micp_) {
251 names_[microbialConcentrationIdx] =
"Microbes";
252 names_[oxygenConcentrationIdx] =
"Oxygen";
253 names_[ureaConcentrationIdx] =
"Urea";
254 names_[biofilmConcentrationIdx] =
"Biofilm";
255 names_[calciteConcentrationIdx] =
"Calcite";
259 const std::string& name(
const int compIdx)
const
265 std::vector<std::string> names_{};
286 : ebosSimulator_(ebosSimulator)
287 , grid_(ebosSimulator_.
vanguard().grid())
292 , current_relaxation_(1.0)
293 , dx_old_(ebosSimulator_.model().
numGridDof())
297 convergence_reports_.reserve(300);
300 bool isParallel()
const
301 {
return grid_.comm().size() > 1; }
303 const EclipseState& eclState()
const
304 {
return ebosSimulator_.vanguard().eclState(); }
314 if (
timer.lastStepFailed() ) {
315 ebosSimulator_.model().updateFailed();
317 ebosSimulator_.model().advanceTimeLevel();
324 ebosSimulator_.setTime(
timer.simulationTimeElapsed());
325 ebosSimulator_.setTimeStepSize(
timer.currentStepLength());
326 ebosSimulator_.model().newtonMethod().setIterationIndex(0);
328 ebosSimulator_.problem().beginTimeStep();
330 unsigned numDof = ebosSimulator_.model().numGridDof();
331 wasSwitched_.resize(
numDof);
332 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
334 if (param_.update_equations_scaling_) {
335 std::cout <<
"equation scaling not supported yet" << std::endl;
338 report.pre_post_time +=
perfTimer.stop();
353 template <
class NonlinearSolverType>
366 residual_norms_history_.clear();
367 current_relaxation_ = 1.0;
369 convergence_reports_.push_back({
timer.reportStepNum(),
timer.currentStepNum(), {}});
370 convergence_reports_.back().report.reserve(11);
373 report.total_linearizations = 1;
377 report.assemble_time +=
perfTimer.stop();
380 report.assemble_time +=
perfTimer.stop();
381 failureReport_ += report;
393 ConvergenceReport::Severity severity =
convrep.severityOfWorstFailure();
394 convergence_reports_.back().report.push_back(std::move(
convrep));
397 if (severity == ConvergenceReport::Severity::NotANumber) {
399 }
else if (severity == ConvergenceReport::Severity::TooLarge) {
405 if (!report.converged) {
408 report.total_newton_iterations = 1;
414 unsigned nc = ebosSimulator_.model().numGridDof();
418 linear_solve_setup_time_ = 0.0;
424 ebosSimulator().model().
linearizer().residual());
427 report.linear_solve_setup_time += linear_solve_setup_time_;
428 report.linear_solve_time +=
perfTimer.stop();
432 report.linear_solve_setup_time += linear_solve_setup_time_;
433 report.linear_solve_time +=
perfTimer.stop();
436 failureReport_ += report;
448 if (param_.use_update_stabilization_) {
455 current_relaxation_ = std::max(current_relaxation_,
nonlinear_solver.relaxMax());
457 std::string
msg =
" Oscillating behavior detected: Relaxation set to "
458 + std::to_string(current_relaxation_);
475 void printIf(
int c,
double x,
double y,
double eps, std::string type) {
476 if (std::abs(x-
y) >
eps) {
477 std::cout << type <<
" " <<
c <<
": "<<x <<
" " <<
y << std::endl;
490 ebosSimulator_.problem().endTimeStep();
491 report.pre_post_time +=
perfTimer.stop();
503 ebosSimulator_.model().newtonMethod().setIterationIndex(
iterationIdx);
504 ebosSimulator_.problem().beginIteration();
505 ebosSimulator_.model().linearizer().linearizeDomain();
506 ebosSimulator_.problem().endIteration();
512 double relativeChange()
const
517 const auto&
elemMapper = ebosSimulator_.model().elementMapper();
518 const auto&
gridView = ebosSimulator_.gridView();
527 Scalar oilSaturationNew = 1.0;
528 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
529 FluidSystem::numActivePhases() > 1 &&
530 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
531 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
532 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
535 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
536 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
537 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
538 assert(Indices::compositionSwitchIdx >= 0 );
539 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
540 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
543 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
544 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
547 const auto& priVarsOld = ebosSimulator_.model().solution(1)[globalElemIdx];
550 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
552 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
553 Scalar oilSaturationOld = 1.0;
556 Scalar tmp = pressureNew - pressureOld;
557 resultDelta += tmp*tmp;
558 resultDenom += pressureNew*pressureNew;
560 if (FluidSystem::numActivePhases() > 1) {
561 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
562 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
563 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
566 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
568 assert(Indices::compositionSwitchIdx >= 0 );
569 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
570 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
573 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
574 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
576 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
577 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
578 resultDelta += tmpSat*tmpSat;
579 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
580 assert(std::isfinite(resultDelta));
581 assert(std::isfinite(resultDenom));
586 resultDelta = gridView.comm().sum(resultDelta);
587 resultDenom = gridView.comm().sum(resultDenom);
589 if (resultDenom > 0.0)
590 return resultDelta/resultDenom;
598 return ebosSimulator_.model().newtonMethod().linearSolver().iterations ();
606 auto&
ebosJac = ebosSimulator_.model().linearizer().jacobian();
607 auto&
ebosResid = ebosSimulator_.model().linearizer().residual();
612 auto&
ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
616 linear_solve_setup_time_ =
perfTimer.stop();
645 ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(0);
646 ebosSimulator_.problem().eclWriter()->mutableEclOutputModule().invalidateLocalData();
656 template <
class CollectiveCommunication>
660 std::vector< Scalar >&
R_sum,
662 std::vector< Scalar >&
B_avg)
669 if( comm.size() > 1 )
685 sumBuffer.push_back( pvSum );
686 sumBuffer.push_back( numAquiferPvSum );
689 comm.sum( sumBuffer.data(), sumBuffer.size() );
692 comm.max( maxBuffer.data(), maxBuffer.size() );
695 for(
int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
697 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
700 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
703 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
705 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
709 pvSum = sumBuffer[sumBuffer.size()-2];
710 numAquiferPvSum = sumBuffer.back();
714 return {pvSum, numAquiferPvSum};
722 std::vector<Scalar>&
B_avg)
727 const auto&
ebosModel = ebosSimulator_.model();
728 const auto&
ebosProblem = ebosSimulator_.problem();
730 const auto&
ebosResid = ebosSimulator_.model().linearizer().residual();
733 const auto&
gridView = ebosSimulator().gridView();
734 OPM_BEGIN_PARALLEL_TRY_CATCH();
737 elemCtx.updatePrimaryIntensiveQuantities(0);
752 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
756 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
765 if constexpr (has_solvent_) {
766 B_avg[ contiSolventEqIdx ] += 1.0 /
intQuants.solventInverseFormationVolumeFactor().value();
768 R_sum[ contiSolventEqIdx ] +=
R2;
771 if constexpr (has_extbo_) {
772 B_avg[ contiZfracEqIdx ] += 1.0 /
fs.invB(FluidSystem::gasPhaseIdx).value();
774 R_sum[ contiZfracEqIdx ] +=
R2;
777 if constexpr (has_polymer_) {
778 B_avg[ contiPolymerEqIdx ] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
780 R_sum[ contiPolymerEqIdx ] +=
R2;
783 if constexpr (has_foam_) {
784 B_avg[ contiFoamEqIdx ] += 1.0 /
fs.invB(FluidSystem::gasPhaseIdx).value();
789 if constexpr (has_brine_) {
790 B_avg[ contiBrineEqIdx ] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
792 R_sum[ contiBrineEqIdx ] +=
R2;
796 if constexpr (has_polymermw_) {
797 static_assert(has_polymer_);
799 B_avg[contiPolymerMWEqIdx] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
804 R_sum[contiPolymerMWEqIdx] +=
R2;
808 if constexpr (has_energy_) {
809 B_avg[ contiEnergyEqIdx ] += 1.0 / (4.182e1);
811 R_sum[ contiEnergyEqIdx ] +=
R2;
815 if constexpr (has_micp_) {
816 B_avg[ contiMicrobialEqIdx ] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
818 R_sum[ contiMicrobialEqIdx ] +=
R1;
820 B_avg[ contiOxygenEqIdx ] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
822 R_sum[ contiOxygenEqIdx ] +=
R2;
824 B_avg[ contiUreaEqIdx ] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
828 B_avg[ contiBiofilmEqIdx ] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
830 R_sum[ contiBiofilmEqIdx ] +=
R4;
832 B_avg[ contiCalciteEqIdx ] += 1.0 /
fs.invB(FluidSystem::waterPhaseIdx).value();
834 R_sum[ contiCalciteEqIdx ] +=
R5;
839 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm());
857 const auto&
ebosModel = ebosSimulator_.model();
858 const auto&
ebosProblem = ebosSimulator_.problem();
859 const auto&
ebosResid = ebosSimulator_.model().linearizer().residual();
860 const auto&
gridView = ebosSimulator().gridView();
863 OPM_BEGIN_PARALLEL_TRY_CATCH();
873 elemCtx.updatePrimaryIntensiveQuantities(0);
892 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModelEbos::ComputeCnvError() failed: ", grid_.comm());
894 return grid_.comm().sum(
errorPV);
900 std::vector<Scalar>&
B_avg,
904 typedef std::vector< Scalar > Vector;
913 convergenceReduction(grid_.comm(),
pvSumLocal,
920 const double tol_mb = param_.tolerance_mb_;
925 const double tol_cnv =
use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
938 ConvergenceReport report{reportTime};
939 using CR = ConvergenceReport;
940 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
941 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
942 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
943 CR::ReservoirFailure::Type::Cnv };
944 double tol[2] = { tol_mb, tol_cnv };
945 for (
int ii : {0, 1}) {
946 if (std::isnan(res[ii])) {
947 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
949 OpmLog::debug(
"NaN residual for " + this->compNames_.name(compIdx) +
" equation.");
951 }
else if (res[ii] > maxResidualAllowed()) {
952 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
954 OpmLog::debug(
"Too large residual for " + this->compNames_.name(compIdx) +
" equation.");
956 }
else if (res[ii] < 0.0) {
957 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
959 OpmLog::debug(
"Negative residual for " + this->compNames_.name(compIdx) +
" equation.");
961 }
else if (res[ii] > tol[ii]) {
962 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
964 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
972 if (iteration == 0) {
973 std::string msg =
"Iter";
974 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
976 msg += this->compNames_.name(compIdx)[0];
979 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
981 msg += this->compNames_.name(compIdx)[0];
986 std::ostringstream ss;
987 const std::streamsize oprec = ss.precision(3);
988 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
989 ss << std::setw(4) << iteration;
990 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
991 ss << std::setw(11) << mass_balance_residual[compIdx];
993 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
994 ss << std::setw(11) << CNV[compIdx];
998 OpmLog::debug(ss.str());
1015 std::vector<Scalar>
B_avg(numEq, 0.0);
1016 auto report = getReservoirConvergence(
timer.simulationTimeElapsed(),
1017 timer.currentStepLength(),
1021 report +=
wellModel().getWellConvergence(
B_avg, report.converged());
1030 return phaseUsage_.num_phases;
1035 std::vector<std::vector<double> >
1042 std::vector<std::vector<double> >
1048 std::vector<std::vector<double> >
regionValues(0, std::vector<double>(0,0.0));
1052 const Simulator& ebosSimulator()
const
1053 {
return ebosSimulator_; }
1055 Simulator& ebosSimulator()
1056 {
return ebosSimulator_; }
1060 {
return failureReport_; }
1062 const std::vector<StepReport>& stepReports()
const
1064 return convergence_reports_;
1067 std::vector<StepReport> getStepReportsDestructively()
const
1069 return std::move(this->convergence_reports_);
1075 Simulator& ebosSimulator_;
1077 const PhaseUsage phaseUsage_;
1078 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1079 static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1080 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1081 static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1082 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1083 static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1084 static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1085 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1087 ModelParameters param_;
1088 SimulatorReportSingle failureReport_;
1091 BlackoilWellModel<TypeTag>& well_model_;
1098 std::vector<std::vector<double>> residual_norms_history_;
1099 double current_relaxation_;
1102 std::vector<StepReport> convergence_reports_;
1111 wellModel()
const {
return well_model_; }
1113 void beginReportStep()
1115 ebosSimulator_.problem().beginEpisode();
1118 void endReportStep()
1120 ebosSimulator_.problem().endEpisode();
1125 bool isNumericalAquiferCell(
const Dune::CpGrid& grid,
const T& elem)
1127 const auto& aquiferCells = grid.sortedNumAquiferCells();
1128 if (aquiferCells.empty())
1132 auto candidate = std::lower_bound(aquiferCells.begin(), aquiferCells.end(),
1134 return candidate != aquiferCells.end() && *candidate == elem.index();
1137 template<
class G,
class T>
1138 typename std::enable_if<!std::is_same<G,Dune::CpGrid>::value,
bool>::type
1139 isNumericalAquiferCell(
const G&,
const T&)
1144 double dpMaxRel()
const {
return param_.dp_max_rel_; }
1145 double dsMax()
const {
return param_.ds_max_; }
1146 double drMaxRel()
const {
return param_.dr_max_rel_; }
1147 double maxResidualAllowed()
const {
return param_.max_residual_allowed_; }
1148 double linear_solve_setup_time_;
1151 std::vector<bool> wasSwitched_;