24#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED
25#define OPM_BLACKOILMODEL_HEADER_INCLUDED
27#include <fmt/format.h>
29#include <opm/common/ErrorMacros.hpp>
30#include <opm/common/Exceptions.hpp>
31#include <opm/common/OpmLog/OpmLog.hpp>
33#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
34#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
36#include <opm/simulators/aquifers/AquiferGridUtils.hpp>
37#include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>
38#include <opm/simulators/flow/BlackoilModelNldd.hpp>
39#include <opm/simulators/flow/BlackoilModelParameters.hpp>
40#include <opm/simulators/flow/countGlobalCells.hpp>
42#include <opm/simulators/flow/NonlinearSolver.hpp>
43#include <opm/simulators/flow/RSTConv.hpp>
44#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
45#include <opm/simulators/timestepping/ConvergenceReport.hpp>
46#include <opm/simulators/timestepping/SimulatorReport.hpp>
47#include <opm/simulators/timestepping/SimulatorTimer.hpp>
49#include <opm/simulators/wells/BlackoilWellModel.hpp>
51#include <opm/simulators/utils/ComponentName.hpp>
52#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
53#include <opm/simulators/utils/ParallelCommunication.hpp>
54#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
56#include <dune/common/timer.hh>
58#include <fmt/format.h>
76namespace Opm::Properties {
80struct FlowProblem {
using InheritsFrom = std::tuple<FlowBaseProblemBlackoil, BlackOilModel>; };
85template<
class TypeTag>
87{
static constexpr bool value =
true; };
89template<
class TypeTag>
91{
static constexpr bool value =
false; };
93template<
class TypeTag>
99template<
class TypeTag>
101{
static constexpr bool value =
false; };
103template<
class TypeTag>
105{
static constexpr bool value =
false; };
107template<
class TypeTag>
109{
static constexpr bool value =
true; };
111template<
class TypeTag>
113{
static constexpr bool value =
false; };
115template<
class TypeTag>
117{
static constexpr bool value =
false; };
119template<
class TypeTag>
121{
static constexpr bool value =
false; };
123template<
class TypeTag>
125{
static constexpr bool value =
false; };
127template<
class TypeTag>
129{
static constexpr bool value =
false; };
131template<
class TypeTag>
133{
static constexpr bool value =
false; };
135template<
class TypeTag>
137{
static constexpr bool value =
true; };
139template<
class TypeTag>
143template<
class TypeTag>
147template<
class TypeTag>
149{
static constexpr bool value =
false; };
161 template <
class TypeTag>
180 static constexpr int numEq = Indices::numEq;
181 static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
182 static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx;
183 static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
184 static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
185 static constexpr int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
186 static constexpr int contiFoamEqIdx = Indices::contiFoamEqIdx;
187 static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
188 static constexpr int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
189 static constexpr int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
190 static constexpr int contiUreaEqIdx = Indices::contiUreaEqIdx;
191 static constexpr int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
192 static constexpr int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
193 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
194 static constexpr int zFractionIdx = Indices::zFractionIdx;
195 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
196 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
197 static constexpr int temperatureIdx = Indices::temperatureIdx;
198 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
199 static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
200 static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
201 static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
202 static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
203 static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
204 static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
206 using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
207 using MatrixBlockType =
typename SparseMatrixAdapter::MatrixBlock;
208 using Mat =
typename SparseMatrixAdapter::IstlMatrix;
209 using BVector = Dune::BlockVector<VectorBlockType>;
228 const bool terminal_output)
229 : simulator_(simulator)
230 , grid_(simulator_.vanguard().grid())
233 , well_model_ (well_model)
234 ,
rst_conv_(simulator_.problem().eclWriter()->collectOnIORank().localIdxToGlobalIdxMapping(),
237 , current_relaxation_(1.0)
238 , dx_old_(simulator_.model().numGridDof())
242 convergence_reports_.reserve(300);
245 if (terminal_output) {
246 OpmLog::info(
"Using Non-Linear Domain Decomposition solver (nldd).");
248 nlddSolver_ = std::make_unique<BlackoilModelNldd<TypeTag>>(*this);
250 if (terminal_output) {
251 OpmLog::info(
"Using Newton nonlinear solver.");
254 OPM_THROW(std::runtime_error,
"Unknown nonlinear solver option: " + param_.
nonlinear_solver_);
259 bool isParallel()
const
260 {
return grid_.comm().size() > 1; }
263 const EclipseState& eclState()
const
264 {
return simulator_.vanguard().eclState(); }
272 Dune::Timer perfTimer;
276 simulator_.model().updateFailed();
278 simulator_.model().advanceTimeLevel();
287 simulator_.model().newtonMethod().setIterationIndex(0);
289 simulator_.problem().beginTimeStep();
291 unsigned numDof = simulator_.model().numGridDof();
292 wasSwitched_.resize(numDof);
293 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
296 OpmLog::error(
"Equation scaling not supported");
304 report.pre_post_time += perfTimer.stop();
306 auto getIdx = [](
unsigned phaseIdx) ->
int
308 if (FluidSystem::phaseIsActive(phaseIdx)) {
309 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
310 return Indices::canonicalToActiveComponentIndex(sIdx);
315 const auto& schedule = simulator_.vanguard().schedule();
318 {getIdx(FluidSystem::oilPhaseIdx),
319 getIdx(FluidSystem::gasPhaseIdx),
320 getIdx(FluidSystem::waterPhaseIdx),
337 Dune::Timer perfTimer;
340 report.total_linearizations = 1;
345 report.assemble_time += perfTimer.stop();
348 report.assemble_time += perfTimer.stop();
349 failureReport_ += report;
354 std::vector<Scalar> residual_norms;
359 auto convrep =
getConvergence(timer, iteration, maxIter, residual_norms);
360 report.converged = convrep.converged() && iteration >= minIter;
361 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
362 convergence_reports_.back().report.push_back(std::move(convrep));
365 if (severity == ConvergenceReport::Severity::NotANumber) {
366 failureReport_ += report;
367 OPM_THROW_PROBLEM(NumericalProblem,
"NaN residual found!");
368 }
else if (severity == ConvergenceReport::Severity::TooLarge) {
369 failureReport_ += report;
370 OPM_THROW_NOLOG(NumericalProblem,
"Too large residual found!");
371 }
else if (severity == ConvergenceReport::Severity::ConvergenceMonitorFailure) {
372 failureReport_ += report;
373 OPM_THROW_PROBLEM(ConvergenceMonitorFailure,
"Total penalty count exceeded cut-off-limit of " + std::to_string(param_.
convergence_monitoring_cutoff_));
376 report.update_time += perfTimer.stop();
377 residual_norms_history_.push_back(residual_norms);
388 template <
class NonlinearSolverType>
391 NonlinearSolverType& nonlinear_solver)
393 if (iteration == 0) {
396 residual_norms_history_.clear();
397 total_penaltyCard_.reset();
398 prev_above_tolerance_ = 0;
399 prev_distance_ = std::numeric_limits<double>::infinity();
400 current_relaxation_ = 1.0;
403 convergence_reports_.back().report.reserve(11);
408 (iteration < this->param_.nldd_num_initial_newton_iter_))
410 result = this->nonlinearIterationNewton(iteration, timer, nonlinear_solver);
413 result = this->
nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver);
422 template <
class NonlinearSolverType>
425 NonlinearSolverType& nonlinear_solver)
430 Dune::Timer perfTimer;
432 this->initialLinearization(report, iteration, nonlinear_solver.minIter(), nonlinear_solver.maxIter(), timer);
435 if (!report.converged) {
438 report.total_newton_iterations = 1;
441 unsigned nc = simulator_.model().numGridDof();
445 linear_solve_setup_time_ = 0.0;
450 wellModel().linearize(simulator().model().linearizer().jacobian(),
451 simulator().model().linearizer().residual());
456 report.linear_solve_setup_time += linear_solve_setup_time_;
457 report.linear_solve_time += perfTimer.stop();
461 report.linear_solve_setup_time += linear_solve_setup_time_;
462 report.linear_solve_time += perfTimer.stop();
465 failureReport_ += report;
479 bool isOscillate =
false;
480 bool isStagnate =
false;
481 nonlinear_solver.detectOscillations(residual_norms_history_, residual_norms_history_.size() - 1, isOscillate, isStagnate);
483 current_relaxation_ -= nonlinear_solver.relaxIncrement();
484 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
486 std::string msg =
" Oscillating behavior detected: Relaxation set to "
487 + std::to_string(current_relaxation_);
491 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
499 report.update_time += perfTimer.stop();
512 Dune::Timer perfTimer;
514 simulator_.problem().endTimeStep();
516 report.pre_post_time += perfTimer.stop();
522 const int iterationIdx)
525 simulator_.model().newtonMethod().setIterationIndex(iterationIdx);
526 simulator_.problem().beginIteration();
527 simulator_.model().linearizer().linearizeDomain();
528 simulator_.problem().endIteration();
533 Scalar relativeChange()
const
535 Scalar resultDelta = 0.0;
536 Scalar resultDenom = 0.0;
538 const auto& elemMapper = simulator_.model().elementMapper();
539 const auto& gridView = simulator_.gridView();
540 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
541 unsigned globalElemIdx = elemMapper.index(elem);
542 const auto& priVarsNew = simulator_.model().solution(0)[globalElemIdx];
545 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
547 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
548 Scalar oilSaturationNew = 1.0;
549 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
550 FluidSystem::numActivePhases() > 1 &&
551 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
552 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
553 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
556 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
557 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
558 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
559 assert(Indices::compositionSwitchIdx >= 0 );
560 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
561 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
564 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
565 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
568 const auto& priVarsOld = simulator_.model().solution(1)[globalElemIdx];
571 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
573 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
574 Scalar oilSaturationOld = 1.0;
577 Scalar tmp = pressureNew - pressureOld;
578 resultDelta += tmp*tmp;
579 resultDenom += pressureNew*pressureNew;
581 if (FluidSystem::numActivePhases() > 1) {
582 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
583 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
584 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
587 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
589 assert(Indices::compositionSwitchIdx >= 0 );
590 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
591 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
594 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
595 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
597 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
598 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
599 resultDelta += tmpSat*tmpSat;
600 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
601 assert(std::isfinite(resultDelta));
602 assert(std::isfinite(resultDenom));
607 resultDelta = gridView.comm().sum(resultDelta);
608 resultDenom = gridView.comm().sum(resultDenom);
610 if (resultDenom > 0.0)
611 return resultDelta/resultDenom;
619 return simulator_.model().newtonMethod().linearSolver().iterations ();
624 double& linearSolveSetupTime()
626 return linear_solve_setup_time_;
634 auto& jacobian = simulator_.model().linearizer().jacobian().istlMatrix();
635 auto& residual = simulator_.model().linearizer().residual();
636 auto& linSolver = simulator_.model().newtonMethod().linearSolver();
638 const int numSolvers = linSolver.numAvailableSolvers();
639 if ((numSolvers > 1) && (linSolver.getSolveCount() % 100 == 0)) {
642 OpmLog::debug(
"\nRunning speed test for comparing available linear solvers.");
645 Dune::Timer perfTimer;
646 std::vector<double> times(numSolvers);
647 std::vector<double> setupTimes(numSolvers);
650 std::vector<BVector> x_trial(numSolvers, x);
651 for (
int solver = 0; solver < numSolvers; ++solver) {
653 linSolver.setActiveSolver(solver);
655 linSolver.prepare(jacobian, residual);
656 setupTimes[solver] = perfTimer.stop();
658 linSolver.setResidual(residual);
660 linSolver.solve(x_trial[solver]);
661 times[solver] = perfTimer.stop();
664 OpmLog::debug(fmt::format(
"Solver time {}: {}", solver, times[solver]));
668 int fastest_solver = std::min_element(times.begin(), times.end()) - times.begin();
670 grid_.comm().broadcast(&fastest_solver, 1, 0);
671 linear_solve_setup_time_ = setupTimes[fastest_solver];
672 x = x_trial[fastest_solver];
673 linSolver.setActiveSolver(fastest_solver);
678 Dune::Timer perfTimer;
680 linSolver.prepare(jacobian, residual);
681 linear_solve_setup_time_ = perfTimer.stop();
682 linSolver.setResidual(residual);
696 auto& newtonMethod = simulator_.model().newtonMethod();
697 SolutionVector& solution = simulator_.model().solution(0);
699 newtonMethod.update_(solution,
708 OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
709 simulator_.model().invalidateAndUpdateIntensiveQuantities(0);
719 std::tuple<Scalar,Scalar> convergenceReduction(Parallel::Communication comm,
720 const Scalar pvSumLocal,
721 const Scalar numAquiferPvSumLocal,
722 std::vector< Scalar >& R_sum,
723 std::vector< Scalar >& maxCoeff,
724 std::vector< Scalar >& B_avg)
726 OPM_TIMEBLOCK(convergenceReduction);
728 Scalar pvSum = pvSumLocal;
729 Scalar numAquiferPvSum = numAquiferPvSumLocal;
731 if( comm.size() > 1 )
734 std::vector< Scalar > sumBuffer;
735 std::vector< Scalar > maxBuffer;
736 const int numComp = B_avg.size();
737 sumBuffer.reserve( 2*numComp + 2 );
738 maxBuffer.reserve( numComp );
739 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
741 sumBuffer.push_back( B_avg[ compIdx ] );
742 sumBuffer.push_back( R_sum[ compIdx ] );
743 maxBuffer.push_back( maxCoeff[ compIdx ] );
747 sumBuffer.push_back( pvSum );
748 sumBuffer.push_back( numAquiferPvSum );
751 comm.sum( sumBuffer.data(), sumBuffer.size() );
754 comm.max( maxBuffer.data(), maxBuffer.size() );
757 for(
int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
759 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
762 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
765 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
767 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
771 pvSum = sumBuffer[sumBuffer.size()-2];
772 numAquiferPvSum = sumBuffer.back();
776 return {pvSum, numAquiferPvSum};
783 std::vector<Scalar>& maxCoeff,
784 std::vector<Scalar>& B_avg,
785 std::vector<int>& maxCoeffCell)
788 Scalar pvSumLocal = 0.0;
789 Scalar numAquiferPvSumLocal = 0.0;
790 const auto& model = simulator_.model();
791 const auto& problem = simulator_.problem();
793 const auto& residual = simulator_.model().linearizer().residual();
795 ElementContext elemCtx(simulator_);
796 const auto& gridView = simulator().gridView();
798 OPM_BEGIN_PARALLEL_TRY_CATCH();
799 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
800 elemCtx.updatePrimaryStencil(elem);
801 elemCtx.updatePrimaryIntensiveQuantities(0);
803 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
804 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
805 const auto& fs = intQuants.fluidState();
807 const auto pvValue = problem.referencePorosity(cell_idx, 0) *
808 model.dofTotalVolume(cell_idx);
809 pvSumLocal += pvValue;
811 if (isNumericalAquiferCell(elem))
813 numAquiferPvSumLocal += pvValue;
816 this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue,
817 B_avg, R_sum, maxCoeff, maxCoeffCell);
820 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModel::localConvergenceData() failed: ", grid_.comm());
823 const int bSize = B_avg.size();
824 for (
int i = 0; i<bSize; ++i )
829 return {pvSumLocal, numAquiferPvSumLocal};
836 std::pair<std::vector<double>, std::vector<int>>
839 OPM_TIMEBLOCK(computeCnvErrorPv);
844 constexpr auto numPvGroups = std::vector<double>::size_type{3};
846 auto cnvPvSplit = std::pair<std::vector<double>, std::vector<int>> {
847 std::piecewise_construct,
848 std::forward_as_tuple(numPvGroups),
849 std::forward_as_tuple(numPvGroups)
852 auto maxCNV = [&B_avg, dt](
const auto& residual,
const double pvol)
855 std::inner_product(residual.begin(), residual.end(),
856 B_avg.begin(), Scalar{0},
857 [](
const Scalar m,
const auto& x)
860 return std::max(m, abs(x));
861 }, std::multiplies<>{});
864 auto& [splitPV, cellCntPV] = cnvPvSplit;
866 const auto& model = this->simulator().model();
867 const auto& problem = this->simulator().problem();
868 const auto& residual = model.linearizer().residual();
869 const auto& gridView = this->simulator().gridView();
873 ElementContext elemCtx(this->simulator());
875 OPM_BEGIN_PARALLEL_TRY_CATCH();
876 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
878 if (isNumericalAquiferCell(elem)) {
882 elemCtx.updatePrimaryStencil(elem);
884 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
885 const auto pvValue = problem.referencePorosity(cell_idx, 0)
886 * model.dofTotalVolume(cell_idx);
888 const auto maxCnv = maxCNV(residual[cell_idx], pvValue);
893 splitPV[ix] +=
static_cast<double>(pvValue);
897 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModel::characteriseCnvPvSplit() failed: ",
900 this->grid_.comm().sum(splitPV .data(), splitPV .size());
901 this->grid_.comm().sum(cellCntPV.data(), cellCntPV.size());
907 void updateTUNING(
const Tuning& tuning)
912 OpmLog::debug(fmt::format(
"Setting BlackoilModel mass "
913 "balance limit (XXXMBE) to {:.2e}",
919 ConvergenceReport getReservoirConvergence(
const double reportTime,
923 std::vector<Scalar>& B_avg,
924 std::vector<Scalar>& residual_norms)
926 OPM_TIMEBLOCK(getReservoirConvergence);
927 using Vector = std::vector<Scalar>;
929 ConvergenceReport report{reportTime};
931 const int numComp = numEq;
933 Vector R_sum(numComp, Scalar{0});
934 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest());
935 std::vector<int> maxCoeffCell(numComp, -1);
937 const auto [pvSumLocal, numAquiferPvSumLocal] =
941 const auto& [pvSum, numAquiferPvSum] =
942 this->convergenceReduction(this->grid_.comm(),
944 numAquiferPvSumLocal,
945 R_sum, maxCoeff, B_avg);
948 pvSum - numAquiferPvSum);
957 const bool relax_final_iteration_mb =
959 && (iteration == maxIter);
961 const bool use_relaxed_mb = relax_final_iteration_mb ||
971 const bool relax_final_iteration_cnv =
973 && (iteration == maxIter);
981 const auto relax_pv_fraction_cnv =
982 [&report,
this, eligible = pvSum - numAquiferPvSum]()
984 const auto& cnvPvSplit = report.cnvPvSplit().first;
988 return static_cast<Scalar
>(cnvPvSplit[1] + cnvPvSplit[2]) <
989 this->param_.relaxed_max_pv_fraction_ * eligible;
992 const bool use_relaxed_cnv = relax_final_iteration_cnv
993 || relax_pv_fraction_cnv
996 if ((relax_final_iteration_mb || relax_final_iteration_cnv) &&
999 std::string message =
1000 "Number of newton iterations reached its maximum "
1001 "try to continue with relaxed tolerances:";
1003 if (relax_final_iteration_mb) {
1007 if (relax_final_iteration_cnv) {
1011 OpmLog::debug(message);
1020 std::vector<Scalar> CNV(numComp);
1021 std::vector<Scalar> mass_balance_residual(numComp);
1022 for (
int compIdx = 0; compIdx < numComp; ++compIdx)
1024 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
1025 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
1026 residual_norms.push_back(CNV[compIdx]);
1029 using CR = ConvergenceReport;
1030 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1031 const Scalar res[2] = {
1032 mass_balance_residual[compIdx], CNV[compIdx],
1035 const CR::ReservoirFailure::Type types[2] = {
1036 CR::ReservoirFailure::Type::MassBalance,
1037 CR::ReservoirFailure::Type::Cnv,
1040 Scalar tol[2] = { tol_mb, tol_cnv, };
1041 if (has_energy_ && compIdx == contiEnergyEqIdx) {
1043 tol[1] = tol_cnv_energy;
1046 for (
int ii : {0, 1}) {
1047 if (std::isnan(res[ii])) {
1048 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
1050 OpmLog::debug(
"NaN residual for " + this->compNames_.name(compIdx) +
" equation.");
1053 else if (res[ii] > maxResidualAllowed()) {
1054 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
1056 OpmLog::debug(
"Too large residual for " + this->compNames_.name(compIdx) +
" equation.");
1059 else if (res[ii] < 0.0) {
1060 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
1062 OpmLog::debug(
"Negative residual for " + this->compNames_.name(compIdx) +
" equation.");
1065 else if (res[ii] > tol[ii]) {
1066 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
1069 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
1076 if (iteration == 0) {
1077 std::string msg =
"Iter";
1078 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1080 msg += this->compNames_.name(compIdx)[0];
1084 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1086 msg += this->compNames_.name(compIdx)[0];
1093 std::ostringstream ss;
1094 const std::streamsize oprec = ss.precision(3);
1095 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1097 ss << std::setw(4) << iteration;
1098 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1099 ss << std::setw(11) << mass_balance_residual[compIdx];
1102 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1103 ss << std::setw(11) << CNV[compIdx];
1106 ss.precision(oprec);
1109 OpmLog::debug(ss.str());
1115 void checkCardPenalty(ConvergenceReport& report,
int iteration)
1118 const auto& current_metrics = report.reservoirConvergence();
1119 auto distances = std::vector<double>(current_metrics.size(), 0.0);
1120 int current_above_tolerance = 0;
1122 for (
size_t i = 0; i < current_metrics.size(); ++i) {
1123 distances[i] = std::max(std::log10(current_metrics[i].value()/current_metrics[i].tolerance()), 0.0);
1125 if (current_metrics[i].value() > current_metrics[i].tolerance()) {
1126 current_above_tolerance++;
1131 double current_distance = std::accumulate(distances.begin(), distances.end(), 0.0);
1133 if (iteration > 0) {
1135 if (current_above_tolerance > prev_above_tolerance_) {
1136 report.addNonConvergedPenalty();
1140 report.addDistanceDecayPenalty();
1144 prev_distance_ = current_distance;
1145 prev_above_tolerance_ = current_above_tolerance;
1147 if (report.wellFailures().size() > 0) {
1148 report.addLargeWellResidualsPenalty();
1151 total_penaltyCard_ += report.getPenaltyCard();
1154 report.setReservoirFailed({ConvergenceReport::ReservoirFailure::Type::ConvergenceMonitorFailure,
1155 ConvergenceReport::Severity::ConvergenceMonitorFailure,
1167 const int iteration,
1169 std::vector<Scalar>& residual_norms)
1173 std::vector<Scalar> B_avg(numEq, 0.0);
1176 iteration, maxIter, B_avg, residual_norms);
1178 OPM_TIMEBLOCK(getWellConvergence);
1179 report +=
wellModel().getWellConvergence(B_avg, report.converged());
1182 checkCardPenalty(report, iteration);
1191 return phaseUsage_.num_phases;
1196 std::vector<std::vector<Scalar> >
1203 std::vector<std::vector<Scalar> >
1209 std::vector<std::vector<Scalar> > regionValues(0, std::vector<Scalar>(0,0.0));
1210 return regionValues;
1214 {
return simulator_; }
1216 Simulator& simulator()
1217 {
return simulator_; }
1221 {
return failureReport_; }
1230 const std::vector<StepReport>& stepReports()
const
1232 return convergence_reports_;
1235 void writePartitions(
const std::filesystem::path& odir)
const
1242 const auto& elementMapper = this->simulator().model().elementMapper();
1243 const auto& cartMapper = this->simulator().vanguard().cartesianIndexMapper();
1245 const auto& grid = this->simulator().vanguard().grid();
1246 const auto& comm = grid.comm();
1247 const auto nDigit = 1 +
static_cast<int>(std::floor(std::log10(comm.size())));
1249 std::ofstream pfile { odir / fmt::format(
"{1:0>{0}}", nDigit, comm.rank()) };
1251 for (
const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
1252 pfile << comm.rank() <<
' '
1253 << cartMapper.cartesianIndex(elementMapper.index(cell)) <<
' '
1254 << comm.rank() <<
'\n';
1258 const std::vector<std::vector<int>>& getConvCells()
const
1264 Simulator& simulator_;
1266 const PhaseUsage phaseUsage_;
1276 ModelParameters param_;
1277 SimulatorReportSingle failureReport_;
1280 BlackoilWellModel<TypeTag>& well_model_;
1289 std::vector<std::vector<Scalar>> residual_norms_history_;
1290 Scalar current_relaxation_;
1293 std::vector<StepReport> convergence_reports_;
1304 wellModel()
const {
return well_model_; }
1306 void beginReportStep()
1308 simulator_.problem().beginEpisode();
1311 void endReportStep()
1313 simulator_.problem().endEpisode();
1316 template<
class Flu
idState,
class Res
idual>
1317 void getMaxCoeff(
const unsigned cell_idx,
1318 const IntensiveQuantities& intQuants,
1319 const FluidState& fs,
1320 const Residual& modelResid,
1321 const Scalar pvValue,
1322 std::vector<Scalar>& B_avg,
1323 std::vector<Scalar>& R_sum,
1324 std::vector<Scalar>& maxCoeff,
1325 std::vector<int>& maxCoeffCell)
1327 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1329 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1333 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1335 B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
1336 const auto R2 = modelResid[cell_idx][compIdx];
1338 R_sum[compIdx] += R2;
1339 const Scalar Rval = std::abs(R2) / pvValue;
1340 if (Rval > maxCoeff[compIdx]) {
1341 maxCoeff[compIdx] = Rval;
1342 maxCoeffCell[compIdx] = cell_idx;
1346 if constexpr (has_solvent_) {
1347 B_avg[contiSolventEqIdx] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1348 const auto R2 = modelResid[cell_idx][contiSolventEqIdx];
1349 R_sum[contiSolventEqIdx] += R2;
1350 maxCoeff[contiSolventEqIdx] = std::max(maxCoeff[contiSolventEqIdx],
1351 std::abs(R2) / pvValue);
1353 if constexpr (has_extbo_) {
1354 B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1355 const auto R2 = modelResid[cell_idx][contiZfracEqIdx];
1356 R_sum[ contiZfracEqIdx ] += R2;
1357 maxCoeff[contiZfracEqIdx] = std::max(maxCoeff[contiZfracEqIdx],
1358 std::abs(R2) / pvValue);
1360 if constexpr (has_polymer_) {
1361 B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1362 const auto R2 = modelResid[cell_idx][contiPolymerEqIdx];
1363 R_sum[contiPolymerEqIdx] += R2;
1364 maxCoeff[contiPolymerEqIdx] = std::max(maxCoeff[contiPolymerEqIdx],
1365 std::abs(R2) / pvValue);
1367 if constexpr (has_foam_) {
1368 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1369 const auto R2 = modelResid[cell_idx][contiFoamEqIdx];
1370 R_sum[contiFoamEqIdx] += R2;
1371 maxCoeff[contiFoamEqIdx] = std::max(maxCoeff[contiFoamEqIdx],
1372 std::abs(R2) / pvValue);
1374 if constexpr (has_brine_) {
1375 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1376 const auto R2 = modelResid[cell_idx][contiBrineEqIdx];
1377 R_sum[contiBrineEqIdx] += R2;
1378 maxCoeff[contiBrineEqIdx] = std::max(maxCoeff[contiBrineEqIdx],
1379 std::abs(R2) / pvValue);
1382 if constexpr (has_polymermw_) {
1383 static_assert(has_polymer_);
1385 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1389 const auto R2 = modelResid[cell_idx][contiPolymerMWEqIdx] / 100.;
1390 R_sum[contiPolymerMWEqIdx] += R2;
1391 maxCoeff[contiPolymerMWEqIdx] = std::max(maxCoeff[contiPolymerMWEqIdx],
1392 std::abs(R2) / pvValue);
1395 if constexpr (has_energy_) {
1396 B_avg[contiEnergyEqIdx] += 1.0 / (4.182e1);
1397 const auto R2 = modelResid[cell_idx][contiEnergyEqIdx];
1398 R_sum[contiEnergyEqIdx] += R2;
1399 maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
1400 std::abs(R2) / pvValue);
1403 if constexpr (has_micp_) {
1404 B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1405 const auto R1 = modelResid[cell_idx][contiMicrobialEqIdx];
1406 R_sum[contiMicrobialEqIdx] += R1;
1407 maxCoeff[contiMicrobialEqIdx] = std::max(maxCoeff[contiMicrobialEqIdx],
1408 std::abs(R1) / pvValue);
1409 B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1410 const auto R2 = modelResid[cell_idx][contiOxygenEqIdx];
1411 R_sum[contiOxygenEqIdx] += R2;
1412 maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
1413 std::abs(R2) / pvValue);
1414 B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1415 const auto R3 = modelResid[cell_idx][contiUreaEqIdx];
1416 R_sum[contiUreaEqIdx] += R3;
1417 maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
1418 std::abs(R3) / pvValue);
1419 B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1420 const auto R4 = modelResid[cell_idx][contiBiofilmEqIdx];
1421 R_sum[contiBiofilmEqIdx] += R4;
1422 maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
1423 std::abs(R4) / pvValue);
1424 B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1425 const auto R5 = modelResid[cell_idx][contiCalciteEqIdx];
1426 R_sum[contiCalciteEqIdx] += R5;
1427 maxCoeff[contiCalciteEqIdx] = std::max(maxCoeff[contiCalciteEqIdx],
1428 std::abs(R5) / pvValue);
1445 Scalar dpMaxRel()
const {
return param_.dp_max_rel_; }
1446 Scalar dsMax()
const {
return param_.ds_max_; }
1447 Scalar drMaxRel()
const {
return param_.dr_max_rel_; }
1448 Scalar maxResidualAllowed()
const {
return param_.max_residual_allowed_; }
1449 double linear_solve_setup_time_;
1450 ConvergenceReport::PenaltyCard total_penaltyCard_;
1451 double prev_distance_ = std::numeric_limits<double>::infinity();
1452 int prev_above_tolerance_ = 0;
1454 std::vector<bool> wasSwitched_;
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Class for handling the blackoil aquifer model.
Definition BlackoilAquiferModel.hpp:50
A model implementation for three-phase black oil.
Definition BlackoilModelNldd.hpp:75
BlackoilModel(Simulator &simulator, const ModelParameters ¶m, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
Construct the model.
Definition BlackoilModel.hpp:225
std::vector< std::vector< Scalar > > computeFluidInPlace(const T &, const std::vector< int > &fipnum) const
Wrapper required due to not following generic API.
Definition BlackoilModel.hpp:1197
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition BlackoilModel.hpp:617
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Called once per nonlinear iteration.
Definition BlackoilModel.hpp:389
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition BlackoilModel.hpp:521
ConvergenceReport getConvergence(const SimulatorTimerInterface &timer, const int iteration, const int maxIter, std::vector< Scalar > &residual_norms)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv).
Definition BlackoilModel.hpp:1166
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition BlackoilModel.hpp:714
int numPhases() const
The number of active fluid phases in the model.
Definition BlackoilModel.hpp:1189
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Called once after each time step.
Definition BlackoilModel.hpp:509
const ComponentName & compNames() const
Returns const reference to component names.
Definition BlackoilModel.hpp:1439
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition BlackoilModel.hpp:693
long int global_nc_
The number of cells of the global grid.
Definition BlackoilModel.hpp:1287
const ModelParameters & param() const
Returns const reference to model parameters.
Definition BlackoilModel.hpp:1433
std::unique_ptr< BlackoilModelNldd< TypeTag > > nlddSolver_
Non-linear DD solver.
Definition BlackoilModel.hpp:1296
std::pair< std::vector< double >, std::vector< int > > characteriseCnvPvSplit(const std::vector< Scalar > &B_avg, const double dt)
Compute pore-volume/cell count split among "converged", "relaxed converged", "unconverged" cells base...
Definition BlackoilModel.hpp:837
const SimulatorReportSingle & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModel.hpp:1220
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition BlackoilModel.hpp:1301
std::vector< std::vector< Scalar > > computeFluidInPlace(const std::vector< int > &) const
Should not be called.
Definition BlackoilModel.hpp:1204
SimulatorReportSingle localAccumulatedReports() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModel.hpp:1224
std::pair< Scalar, Scalar > localConvergenceData(std::vector< Scalar > &R_sum, std::vector< Scalar > &maxCoeff, std::vector< Scalar > &B_avg, std::vector< int > &maxCoeffCell)
Get reservoir quantities on this process needed for convergence calculations.
Definition BlackoilModel.hpp:782
void solveJacobianSystem(BVector &x)
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition BlackoilModel.hpp:632
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Called once before each time step.
Definition BlackoilModel.hpp:269
RSTConv rst_conv_
Helper class for RPTRST CONV.
Definition BlackoilModel.hpp:1282
bool terminal_output_
Whether we print something to std::cout.
Definition BlackoilModel.hpp:1285
Class for handling the blackoil well model.
Definition WellInterface.hpp:31
Definition ComponentName.hpp:34
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition ConvergenceReport.hpp:38
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:92
Class computing RPTRST CONV output.
Definition RSTConv.hpp:34
const std::vector< std::vector< int > > & getData() const
Obtain a const-ref to the accumulated data.
Definition RSTConv.hpp:58
void update(const ResidualVector &resid)
Adds the CONV output for given residual vector.
Definition RSTConv.cpp:58
void init(const std::size_t numCells, const RSTConfig &rst_config, const std::array< int, 6 > &compIdx)
Init state at beginning of step.
Definition RSTConv.cpp:35
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.
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition phaseUsageFromDeck.cpp:137
Solver parameters for the BlackoilModel.
Definition BlackoilModelParameters.hpp:162
Scalar tolerance_mb_relaxed_
Relaxed mass balance tolerance (can be used when iter >= min_strict_mb_iter_).
Definition BlackoilModelParameters.hpp:182
Scalar tolerance_energy_balance_
Relative energy balance tolerance (total energy balance error).
Definition BlackoilModelParameters.hpp:184
bool update_equations_scaling_
Update scaling factors for mass balance equations.
Definition BlackoilModelParameters.hpp:249
Scalar tolerance_energy_balance_relaxed_
Relaxed energy balance tolerance (can be used when iter >= min_strict_mb_iter_).
Definition BlackoilModelParameters.hpp:186
Scalar convergence_monitoring_decay_factor_
Decay factor used in convergence monitoring.
Definition BlackoilModelParameters.hpp:316
int min_strict_mb_iter_
Minimum number of Newton iterations before we can use relaxed MB convergence criterion.
Definition BlackoilModelParameters.hpp:243
int min_strict_cnv_iter_
Minimum number of Newton iterations before we can use relaxed CNV convergence criterion.
Definition BlackoilModelParameters.hpp:240
bool use_update_stabilization_
Try to detect oscillation or stagnation.
Definition BlackoilModelParameters.hpp:252
Scalar tolerance_cnv_energy_
Local energy convergence tolerance (max of local energy errors).
Definition BlackoilModelParameters.hpp:192
Scalar tolerance_cnv_energy_relaxed_
Relaxed local energy convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViola...
Definition BlackoilModelParameters.hpp:194
Scalar tolerance_mb_
Relative mass balance tolerance (total mass balance error).
Definition BlackoilModelParameters.hpp:180
bool convergence_monitoring_
Whether to enable convergence monitoring.
Definition BlackoilModelParameters.hpp:312
Scalar tolerance_cnv_
Local convergence tolerance (max of local saturation errors).
Definition BlackoilModelParameters.hpp:188
Scalar tolerance_cnv_relaxed_
Relaxed local convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViolatedPV <...
Definition BlackoilModelParameters.hpp:190
int convergence_monitoring_cutoff_
Cut-off limit for convergence monitoring.
Definition BlackoilModelParameters.hpp:314
std::string nonlinear_solver_
Nonlinear solver type: newton or nldd.
Definition BlackoilModelParameters.hpp:294
Definition AquiferGridUtils.hpp:35
Definition FlowBaseProblemProperties.hpp:62
Enable surface volume scaling.
Definition blackoilproperties.hh:54
Enable the ECL-blackoil extension for salt.
Definition blackoilproperties.hh:60
Enable convective mixing?
Definition multiphasebaseproperties.hh:85
Definition FlowBaseProblemProperties.hpp:73
Enable dispersive fluxes?
Definition multiphasebaseproperties.hh:82
Specify whether energy should be considered as a conservation quantity or not.
Definition multiphasebaseproperties.hh:76
Enable the ECL-blackoil extension for foam.
Definition blackoilproperties.hh:57
Enable the ECL-blackoil extension for MICP.
Definition blackoilproperties.hh:72
Enable the ECL-blackoil extension for polymer.
Definition blackoilproperties.hh:48
Enable the ECL-blackoil extension for salt precipitation.
Definition blackoilproperties.hh:63
Enable the ECL-blackoil extension for solvents. ("Second gas")
Definition blackoilproperties.hh:42
Allow the spatial and temporal domains to exhibit non-constant temperature in the black-oil model.
Definition blackoilproperties.hh:78
Definition fvbaseproperties.hh:53
Definition ISTLSolver.hpp:63
Definition BlackoilModel.hpp:80
Specify whether to use volumetric residuals or not.
Definition fvbaseproperties.hh:237
Definition ISTLSolver.hpp:69
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34