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>
56#include <dune/common/timer.hh>
58#include <fmt/format.h>
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;
208 using Mat =
typename SparseMatrixAdapter::IstlMatrix;
209 using BVector = Dune::BlockVector<VectorBlockType>;
228 const bool terminal_output)
234 ,
rst_conv_(
simulator_.problem().eclWriter()->collectOnIORank().localIdxToGlobalIdxMapping(),
246 OPM_THROW(std::runtime_error,
"The --nonlinear-solver=nldd option can only be used with --matrix-add-well-contributions=true");
248 if (terminal_output) {
249 OpmLog::info(
"Using Non-Linear Domain Decomposition solver (nldd).");
251 nlddSolver_ = std::make_unique<BlackoilModelNldd<TypeTag>>(*this);
253 if (terminal_output) {
254 OpmLog::info(
"Using Newton nonlinear solver.");
263 {
return grid_.comm().size() > 1; }
275 Dune::Timer perfTimer;
290 simulator_.model().newtonMethod().setIterationIndex(0);
294 unsigned numDof =
simulator_.model().numGridDof();
299 OpmLog::error(
"Equation scaling not supported");
309 auto getIdx = [](
unsigned phaseIdx) ->
int
311 if (FluidSystem::phaseIsActive(phaseIdx)) {
312 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
313 return Indices::canonicalToActiveComponentIndex(sIdx);
318 const auto& schedule =
simulator_.vanguard().schedule();
321 {getIdx(FluidSystem::oilPhaseIdx),
322 getIdx(FluidSystem::gasPhaseIdx),
323 getIdx(FluidSystem::waterPhaseIdx),
340 Dune::Timer perfTimer;
357 std::vector<Scalar> residual_norms;
362 auto convrep =
getConvergence(timer, iteration, maxIter, residual_norms);
363 report.
converged = convrep.converged() && iteration >= minIter;
369 OPM_THROW_PROBLEM(NumericalProblem,
"NaN residual found!");
371 OPM_THROW_NOLOG(NumericalProblem,
"Too large residual found!");
386 template <
class NonlinearSolverType>
389 NonlinearSolverType& nonlinear_solver)
391 if (iteration == 0) {
408 result = this->
nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver);
417 template <
class NonlinearSolverType>
420 NonlinearSolverType& nonlinear_solver)
425 Dune::Timer perfTimer;
427 this->
initialLinearization(report, iteration, nonlinear_solver.minIter(), nonlinear_solver.maxIter(), timer);
436 unsigned nc =
simulator_.model().numGridDof();
440 linear_solve_setup_time_ = 0.0;
446 simulator().model().linearizer().residual());
474 bool isOscillate =
false;
475 bool isStagnate =
false;
481 std::string msg =
" Oscillating behavior detected: Relaxation set to "
507 Dune::Timer perfTimer;
517 const int iterationIdx)
520 simulator_.model().newtonMethod().setIterationIndex(iterationIdx);
522 simulator_.model().linearizer().linearizeDomain();
533 const auto& elemMapper =
simulator_.model().elementMapper();
535 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
536 unsigned globalElemIdx = elemMapper.index(elem);
537 const auto& priVarsNew =
simulator_.model().solution(0)[globalElemIdx];
540 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
542 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
543 Scalar oilSaturationNew = 1.0;
544 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
545 FluidSystem::numActivePhases() > 1 &&
546 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
547 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
548 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
551 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
552 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
553 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
554 assert(Indices::compositionSwitchIdx >= 0 );
555 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
556 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
559 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
560 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
563 const auto& priVarsOld =
simulator_.model().solution(1)[globalElemIdx];
566 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
568 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
569 Scalar oilSaturationOld = 1.0;
572 Scalar tmp = pressureNew - pressureOld;
573 resultDelta += tmp*tmp;
574 resultDenom += pressureNew*pressureNew;
576 if (FluidSystem::numActivePhases() > 1) {
577 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
578 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
579 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
582 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
584 assert(Indices::compositionSwitchIdx >= 0 );
585 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
586 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
589 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
590 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
592 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
593 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
594 resultDelta += tmpSat*tmpSat;
595 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
596 assert(std::isfinite(resultDelta));
597 assert(std::isfinite(resultDenom));
602 resultDelta = gridView.comm().sum(resultDelta);
603 resultDenom = gridView.comm().sum(resultDenom);
605 if (resultDenom > 0.0)
606 return resultDelta/resultDenom;
614 return simulator_.model().newtonMethod().linearSolver().iterations ();
621 return linear_solve_setup_time_;
629 auto& jacobian =
simulator_.model().linearizer().jacobian().istlMatrix();
630 auto& residual =
simulator_.model().linearizer().residual();
631 auto& linSolver =
simulator_.model().newtonMethod().linearSolver();
633 const int numSolvers = linSolver.numAvailableSolvers();
634 if ((numSolvers > 1) && (linSolver.getSolveCount() % 100 == 0)) {
637 OpmLog::debug(
"\nRunning speed test for comparing available linear solvers.");
640 Dune::Timer perfTimer;
641 std::vector<double> times(numSolvers);
642 std::vector<double> setupTimes(numSolvers);
645 std::vector<BVector> x_trial(numSolvers, x);
646 for (
int solver = 0; solver < numSolvers; ++solver) {
648 linSolver.setActiveSolver(solver);
650 linSolver.prepare(jacobian, residual);
651 setupTimes[solver] = perfTimer.stop();
653 linSolver.setResidual(residual);
655 linSolver.solve(x_trial[solver]);
656 times[solver] = perfTimer.stop();
659 OpmLog::debug(fmt::format(
"Solver time {}: {}", solver, times[solver]));
663 int fastest_solver = std::min_element(times.begin(), times.end()) - times.begin();
665 grid_.comm().broadcast(&fastest_solver, 1, 0);
666 linear_solve_setup_time_ = setupTimes[fastest_solver];
667 x = x_trial[fastest_solver];
668 linSolver.setActiveSolver(fastest_solver);
673 Dune::Timer perfTimer;
675 linSolver.prepare(jacobian, residual);
676 linear_solve_setup_time_ = perfTimer.stop();
677 linSolver.setResidual(residual);
691 auto& newtonMethod =
simulator_.model().newtonMethod();
694 newtonMethod.update_(solution,
703 OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
704 simulator_.model().invalidateAndUpdateIntensiveQuantities(0);
705 simulator_.problem().eclWriter()->mutableOutputModule().invalidateLocalData();
717 const Scalar numAquiferPvSumLocal,
718 std::vector< Scalar >& R_sum,
719 std::vector< Scalar >& maxCoeff,
720 std::vector< Scalar >& B_avg)
724 Scalar pvSum = pvSumLocal;
725 Scalar numAquiferPvSum = numAquiferPvSumLocal;
727 if( comm.size() > 1 )
730 std::vector< Scalar > sumBuffer;
731 std::vector< Scalar > maxBuffer;
732 const int numComp = B_avg.size();
733 sumBuffer.reserve( 2*numComp + 2 );
734 maxBuffer.reserve( numComp );
735 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
737 sumBuffer.push_back( B_avg[ compIdx ] );
738 sumBuffer.push_back( R_sum[ compIdx ] );
739 maxBuffer.push_back( maxCoeff[ compIdx ] );
743 sumBuffer.push_back( pvSum );
744 sumBuffer.push_back( numAquiferPvSum );
747 comm.sum( sumBuffer.data(), sumBuffer.size() );
750 comm.max( maxBuffer.data(), maxBuffer.size() );
753 for(
int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
755 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
758 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
761 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
763 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
767 pvSum = sumBuffer[sumBuffer.size()-2];
768 numAquiferPvSum = sumBuffer.back();
772 return {pvSum, numAquiferPvSum};
779 std::vector<Scalar>& maxCoeff,
780 std::vector<Scalar>& B_avg,
781 std::vector<int>& maxCoeffCell)
785 Scalar numAquiferPvSumLocal = 0.0;
789 const auto& residual =
simulator_.model().linearizer().residual();
792 const auto& gridView =
simulator().gridView();
795 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
796 elemCtx.updatePrimaryStencil(elem);
797 elemCtx.updatePrimaryIntensiveQuantities(0);
799 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
800 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
801 const auto& fs = intQuants.fluidState();
803 const auto pvValue = problem.referencePorosity(cell_idx, 0) *
804 model.dofTotalVolume(cell_idx);
805 pvSumLocal += pvValue;
807 if (isNumericalAquiferCell(elem))
809 numAquiferPvSumLocal += pvValue;
812 this->
getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue,
813 B_avg, R_sum, maxCoeff, maxCoeffCell);
819 const int bSize = B_avg.size();
820 for (
int i = 0; i<bSize; ++i )
825 return {pvSumLocal, numAquiferPvSumLocal};
832 std::pair<std::vector<double>, std::vector<int>>
835 OPM_TIMEBLOCK(computeCnvErrorPv);
840 constexpr auto numPvGroups = std::vector<double>::size_type{3};
842 auto cnvPvSplit = std::pair<std::vector<double>, std::vector<int>> {
843 std::piecewise_construct,
844 std::forward_as_tuple(numPvGroups),
845 std::forward_as_tuple(numPvGroups)
848 auto maxCNV = [&B_avg, dt](
const auto& residual,
const double pvol)
851 std::inner_product(residual.begin(), residual.end(),
853 [](
const Scalar m,
const auto& x)
856 return std::max(m, abs(x));
857 }, std::multiplies<>{});
860 auto& [splitPV, cellCntPV] = cnvPvSplit;
862 const auto& model = this->
simulator().model();
863 const auto& problem = this->
simulator().problem();
864 const auto& residual = model.linearizer().residual();
865 const auto& gridView = this->
simulator().gridView();
872 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
874 if (isNumericalAquiferCell(elem)) {
878 elemCtx.updatePrimaryStencil(elem);
880 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
881 const auto pvValue = problem.referencePorosity(cell_idx, 0)
882 * model.dofTotalVolume(cell_idx);
884 const auto maxCnv = maxCNV(residual[cell_idx], pvValue);
889 splitPV[ix] +=
static_cast<double>(pvValue);
896 this->
grid_.comm().sum(splitPV .data(), splitPV .size());
897 this->
grid_.comm().sum(cellCntPV.data(), cellCntPV.size());
908 OpmLog::debug(fmt::format(
"Setting BlackoilModel mass "
909 "balance limit (XXXMBE) to {:.2e}",
919 std::vector<Scalar>& B_avg,
920 std::vector<Scalar>& residual_norms)
923 using Vector = std::vector<Scalar>;
927 const int numComp =
numEq;
929 Vector R_sum(numComp,
Scalar{0});
930 Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest());
931 std::vector<int> maxCoeffCell(numComp, -1);
933 const auto [pvSumLocal, numAquiferPvSumLocal] =
937 const auto& [pvSum, numAquiferPvSum] =
940 numAquiferPvSumLocal,
941 R_sum, maxCoeff, B_avg);
944 pvSum - numAquiferPvSum);
953 const bool relax_final_iteration_mb =
955 && (iteration == maxIter);
957 const bool use_relaxed_mb = relax_final_iteration_mb ||
967 const bool relax_final_iteration_cnv =
969 && (iteration == maxIter);
977 const auto relax_pv_fraction_cnv =
978 [&report,
this, eligible = pvSum - numAquiferPvSum]()
980 const auto& cnvPvSplit = report.cnvPvSplit().first;
984 return static_cast<Scalar>(cnvPvSplit[1] + cnvPvSplit[2]) <
988 const bool use_relaxed_cnv = relax_final_iteration_cnv
989 || relax_pv_fraction_cnv
992 if ((relax_final_iteration_mb || relax_final_iteration_cnv) &&
995 std::string message =
996 "Number of newton iterations reached its maximum "
997 "try to continue with relaxed tolerances:";
999 if (relax_final_iteration_mb) {
1003 if (relax_final_iteration_cnv) {
1007 OpmLog::debug(message);
1016 std::vector<Scalar> CNV(numComp);
1017 std::vector<Scalar> mass_balance_residual(numComp);
1018 for (
int compIdx = 0; compIdx < numComp; ++compIdx)
1020 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
1021 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
1022 residual_norms.push_back(CNV[compIdx]);
1026 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1028 mass_balance_residual[compIdx], CNV[compIdx],
1031 const CR::ReservoirFailure::Type types[2] = {
1032 CR::ReservoirFailure::Type::MassBalance,
1033 CR::ReservoirFailure::Type::Cnv,
1036 Scalar tol[2] = { tol_mb, tol_cnv, };
1039 tol[1] = tol_cnv_energy;
1042 for (
int ii : {0, 1}) {
1043 if (std::isnan(res[ii])) {
1046 OpmLog::debug(
"NaN residual for " + this->
compNames_.
name(compIdx) +
" equation.");
1049 else if (res[ii] > maxResidualAllowed()) {
1050 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
1052 OpmLog::debug(
"Too large residual for " + this->
compNames_.
name(compIdx) +
" equation.");
1055 else if (res[ii] < 0.0) {
1056 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
1058 OpmLog::debug(
"Negative residual for " + this->
compNames_.
name(compIdx) +
" equation.");
1061 else if (res[ii] > tol[ii]) {
1062 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
1065 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
1072 if (iteration == 0) {
1073 std::string msg =
"Iter";
1074 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1080 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1089 std::ostringstream ss;
1090 const std::streamsize oprec = ss.precision(3);
1091 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1093 ss << std::setw(4) << iteration;
1094 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1095 ss << std::setw(11) << mass_balance_residual[compIdx];
1098 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1099 ss << std::setw(11) << CNV[compIdx];
1102 ss.precision(oprec);
1105 OpmLog::debug(ss.str());
1118 const int iteration,
1120 std::vector<Scalar>& residual_norms)
1124 std::vector<Scalar> B_avg(
numEq, 0.0);
1127 iteration, maxIter, B_avg, residual_norms);
1129 OPM_TIMEBLOCK(getWellConvergence);
1130 report +=
wellModel().getWellConvergence(B_avg, report.converged());
1144 std::vector<std::vector<Scalar> >
1151 std::vector<std::vector<Scalar> >
1157 std::vector<std::vector<Scalar> > regionValues(0, std::vector<Scalar>(0,0.0));
1158 return regionValues;
1190 const auto& elementMapper = this->
simulator().model().elementMapper();
1191 const auto& cartMapper = this->
simulator().vanguard().cartesianIndexMapper();
1193 const auto& grid = this->
simulator().vanguard().grid();
1194 const auto& comm = grid.comm();
1195 const auto nDigit = 1 +
static_cast<int>(std::floor(std::log10(comm.size())));
1197 std::ofstream pfile { odir / fmt::format(
"{1:0>{0}}", nDigit, comm.rank()) };
1199 for (
const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
1200 pfile << comm.rank() <<
' '
1201 << cartMapper.cartesianIndex(elementMapper.index(cell)) <<
' '
1202 << comm.rank() <<
'\n';
1215 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1216 static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1217 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1218 static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1219 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1220 static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1221 static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1222 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1264 template<
class Flu
idState,
class Res
idual>
1267 const FluidState& fs,
1268 const Residual& modelResid,
1270 std::vector<Scalar>& B_avg,
1271 std::vector<Scalar>& R_sum,
1272 std::vector<Scalar>& maxCoeff,
1273 std::vector<int>& maxCoeffCell)
1275 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1277 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1281 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1283 B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
1284 const auto R2 = modelResid[cell_idx][compIdx];
1286 R_sum[compIdx] += R2;
1287 const Scalar Rval = std::abs(R2) / pvValue;
1288 if (Rval > maxCoeff[compIdx]) {
1289 maxCoeff[compIdx] = Rval;
1290 maxCoeffCell[compIdx] = cell_idx;
1295 B_avg[
contiSolventEqIdx] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1299 std::abs(R2) / pvValue);
1302 B_avg[
contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1306 std::abs(R2) / pvValue);
1313 std::abs(R2) / pvValue);
1316 B_avg[
contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1320 std::abs(R2) / pvValue);
1323 B_avg[
contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1327 std::abs(R2) / pvValue);
1340 std::abs(R2) / pvValue);
1348 std::abs(R2) / pvValue);
1356 std::abs(R1) / pvValue);
1357 B_avg[
contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1361 std::abs(R2) / pvValue);
1362 B_avg[
contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1366 std::abs(R3) / pvValue);
1371 std::abs(R4) / pvValue);
1376 std::abs(R5) / pvValue);
1397 double linear_solve_setup_time_;
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:192
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
Class for handling the blackoil aquifer model.
Definition: BlackoilAquiferModel.hpp:50
Definition: BlackoilModel.hpp:163
static constexpr int contiEnergyEqIdx
Definition: BlackoilModel.hpp:184
SimulatorReportSingle failureReport_
Definition: BlackoilModel.hpp:1225
BlackoilModel(Simulator &simulator, const ModelParameters ¶m, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
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:1145
static constexpr bool has_energy_
Definition: BlackoilModel.hpp:1219
ModelParameters param_
Definition: BlackoilModel.hpp:1224
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition: BlackoilModel.hpp:612
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Definition: BlackoilModel.hpp:387
static constexpr bool has_extbo_
Definition: BlackoilModel.hpp:1216
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilModel.hpp:516
GetPropType< TypeTag, Properties::SparseMatrixAdapter > SparseMatrixAdapter
Definition: BlackoilModel.hpp:170
ConvergenceReport getConvergence(const SimulatorTimerInterface &timer, const int iteration, const int maxIter, std::vector< Scalar > &residual_norms)
Definition: BlackoilModel.hpp:1117
ConvergenceReport getReservoirConvergence(const double reportTime, const double dt, const int iteration, const int maxIter, std::vector< Scalar > &B_avg, std::vector< Scalar > &residual_norms)
Definition: BlackoilModel.hpp:915
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: BlackoilModel.hpp:169
static constexpr bool has_foam_
Definition: BlackoilModel.hpp:1220
static constexpr int foamConcentrationIdx
Definition: BlackoilModel.hpp:198
Scalar current_relaxation_
Definition: BlackoilModel.hpp:1238
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: BlackoilModel.hpp:173
static constexpr int microbialConcentrationIdx
Definition: BlackoilModel.hpp:200
static constexpr int numEq
Definition: BlackoilModel.hpp:180
static constexpr bool has_polymermw_
Definition: BlackoilModel.hpp:1218
static constexpr int biofilmConcentrationIdx
Definition: BlackoilModel.hpp:203
static constexpr int contiFoamEqIdx
Definition: BlackoilModel.hpp:186
static constexpr int polymerMoleWeightIdx
Definition: BlackoilModel.hpp:196
void endReportStep()
Definition: BlackoilModel.hpp:1259
static constexpr int contiCalciteEqIdx
Definition: BlackoilModel.hpp:192
const std::vector< std::vector< int > > & getConvCells() const
Definition: BlackoilModel.hpp:1206
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: BlackoilModel.hpp:168
bool isParallel() const
Definition: BlackoilModel.hpp:262
const std::vector< StepReport > & stepReports() const
Definition: BlackoilModel.hpp:1178
Scalar relativeChange() const
Definition: BlackoilModel.hpp:528
std::vector< StepReport > convergence_reports_
Definition: BlackoilModel.hpp:1241
static constexpr int contiMicrobialEqIdx
Definition: BlackoilModel.hpp:188
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModel.hpp:710
double & linearSolveSetupTime()
Definition: BlackoilModel.hpp:619
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModel.hpp:1137
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Definition: BlackoilModel.hpp:504
static constexpr int contiBiofilmEqIdx
Definition: BlackoilModel.hpp:191
const Grid & grid_
Definition: BlackoilModel.hpp:1213
static constexpr int contiUreaEqIdx
Definition: BlackoilModel.hpp:190
static constexpr int calciteConcentrationIdx
Definition: BlackoilModel.hpp:204
static constexpr int temperatureIdx
Definition: BlackoilModel.hpp:197
GetPropType< TypeTag, Properties::SolutionVector > SolutionVector
Definition: BlackoilModel.hpp:171
const ComponentName & compNames() const
Returns const reference to component names.
Definition: BlackoilModel.hpp:1387
GetPropType< TypeTag, Properties::PrimaryVariables > PrimaryVariables
Definition: BlackoilModel.hpp:172
void initialLinearization(SimulatorReportSingle &report, const int iteration, const int minIter, const int maxIter, const SimulatorTimerInterface &timer)
Definition: BlackoilModel.hpp:332
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition: BlackoilModel.hpp:688
long int global_nc_
The number of cells of the global grid.
Definition: BlackoilModel.hpp:1235
GetPropType< TypeTag, Properties::Indices > Indices
Definition: BlackoilModel.hpp:174
void updateTUNING(const Tuning &tuning)
Definition: BlackoilModel.hpp:903
typename SparseMatrixAdapter::MatrixBlock MatrixBlockType
Definition: BlackoilModel.hpp:207
const ModelParameters & param() const
Returns const reference to model parameters.
Definition: BlackoilModel.hpp:1381
static constexpr int contiOxygenEqIdx
Definition: BlackoilModel.hpp:189
void getMaxCoeff(const unsigned cell_idx, const IntensiveQuantities &intQuants, const FluidState &fs, const Residual &modelResid, const Scalar pvValue, std::vector< Scalar > &B_avg, std::vector< Scalar > &R_sum, std::vector< Scalar > &maxCoeff, std::vector< int > &maxCoeffCell)
Definition: BlackoilModel.hpp:1265
std::unique_ptr< BlackoilModelNldd< TypeTag > > nlddSolver_
Non-linear DD solver.
Definition: BlackoilModel.hpp:1244
SimulatorReportSingle nonlinearIterationNewton(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Definition: BlackoilModel.hpp:418
void writePartitions(const std::filesystem::path &odir) const
Definition: BlackoilModel.hpp:1183
static constexpr int solventSaturationIdx
Definition: BlackoilModel.hpp:193
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:833
void beginReportStep()
Definition: BlackoilModel.hpp:1254
static constexpr int polymerConcentrationIdx
Definition: BlackoilModel.hpp:195
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: BlackoilModel.hpp:166
const SimulatorReportSingle & failureReport() const
return the statistics if the nonlinearIteration() method failed
Definition: BlackoilModel.hpp:1168
Simulator & simulator()
Definition: BlackoilModel.hpp:1164
GetPropType< TypeTag, Properties::MaterialLaw > MaterialLaw
Definition: BlackoilModel.hpp:175
static constexpr int contiPolymerMWEqIdx
Definition: BlackoilModel.hpp:185
static constexpr bool has_micp_
Definition: BlackoilModel.hpp:1222
const Simulator & simulator() const
Definition: BlackoilModel.hpp:1161
BlackoilWellModel< TypeTag > & well_model_
Definition: BlackoilModel.hpp:1228
BlackoilWellModel< TypeTag > & wellModel()
return the StandardWells object
Definition: BlackoilModel.hpp:1249
std::vector< std::vector< Scalar > > residual_norms_history_
Definition: BlackoilModel.hpp:1237
typename SparseMatrixAdapter::IstlMatrix Mat
Definition: BlackoilModel.hpp:208
Dune::FieldVector< Scalar, numEq > VectorBlockType
Definition: BlackoilModel.hpp:206
static constexpr bool has_solvent_
Definition: BlackoilModel.hpp:1215
const BlackoilWellModel< TypeTag > & wellModel() const
Definition: BlackoilModel.hpp:1252
BVector dx_old_
Definition: BlackoilModel.hpp:1239
Simulator & simulator_
Definition: BlackoilModel.hpp:1212
std::vector< std::vector< Scalar > > computeFluidInPlace(const std::vector< int > &) const
Should not be called.
Definition: BlackoilModel.hpp:1152
SimulatorReportSingle localAccumulatedReports() const
return the statistics if the nonlinearIteration() method failed
Definition: BlackoilModel.hpp:1172
static constexpr bool has_polymer_
Definition: BlackoilModel.hpp:1217
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:778
const PhaseUsage phaseUsage_
Definition: BlackoilModel.hpp:1214
static constexpr bool has_brine_
Definition: BlackoilModel.hpp:1221
void solveJacobianSystem(BVector &x)
Definition: BlackoilModel.hpp:627
static constexpr int oxygenConcentrationIdx
Definition: BlackoilModel.hpp:201
ComponentName compNames_
Definition: BlackoilModel.hpp:1242
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Definition: BlackoilModel.hpp:272
RSTConv rst_conv_
Helper class for RPTRST CONV.
Definition: BlackoilModel.hpp:1230
static constexpr int ureaConcentrationIdx
Definition: BlackoilModel.hpp:202
std::tuple< Scalar, Scalar > convergenceReduction(Parallel::Communication comm, const Scalar pvSumLocal, const Scalar numAquiferPvSumLocal, std::vector< Scalar > &R_sum, std::vector< Scalar > &maxCoeff, std::vector< Scalar > &B_avg)
Definition: BlackoilModel.hpp:715
GetPropType< TypeTag, Properties::Grid > Grid
Definition: BlackoilModel.hpp:167
Dune::BlockVector< VectorBlockType > BVector
Definition: BlackoilModel.hpp:209
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: BlackoilModel.hpp:177
static constexpr int saltConcentrationIdx
Definition: BlackoilModel.hpp:199
static constexpr int zFractionIdx
Definition: BlackoilModel.hpp:194
static constexpr int contiBrineEqIdx
Definition: BlackoilModel.hpp:187
static constexpr int contiZfracEqIdx
Definition: BlackoilModel.hpp:182
static constexpr int contiSolventEqIdx
Definition: BlackoilModel.hpp:181
GetPropType< TypeTag, Properties::MaterialLawParams > MaterialLawParams
Definition: BlackoilModel.hpp:176
static constexpr int contiPolymerEqIdx
Definition: BlackoilModel.hpp:183
std::vector< bool > wasSwitched_
Definition: BlackoilModel.hpp:1400
const EclipseState & eclState() const
Definition: BlackoilModel.hpp:266
bool terminal_output_
Whether we print something to std::cout.
Definition: BlackoilModel.hpp:1233
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:98
Definition: ComponentName.hpp:34
const std::string & name(const int compIdx) const
Definition: ComponentName.hpp:38
Definition: ConvergenceReport.hpp:38
Severity
Definition: ConvergenceReport.hpp:49
void setReservoirFailed(const ReservoirFailure &rf)
Definition: ConvergenceReport.hpp:191
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowProblem.hpp:91
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.
void init(const std::size_t numCells, const RSTConfig &rst_config, const std::array< int, 6 > &compIdx)
Init state at beginning of step.
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
virtual double simulationTimeElapsed() const =0
virtual int currentStepNum() const =0
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilmodel.hh:72
std::size_t countGlobalCells(const Grid &grid)
Get the number of cells of a global grid.
Definition: countGlobalCells.hpp:82
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
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:147
Scalar tolerance_mb_relaxed_
Relaxed mass balance tolerance (can be used when iter >= min_strict_mb_iter_).
Definition: BlackoilModelParameters.hpp:161
Scalar tolerance_energy_balance_
Relative energy balance tolerance (total energy balance error).
Definition: BlackoilModelParameters.hpp:163
bool matrix_add_well_contributions_
Whether to add influences of wells between cells to the matrix and preconditioner matrix.
Definition: BlackoilModelParameters.hpp:244
bool update_equations_scaling_
Update scaling factors for mass balance equations.
Definition: BlackoilModelParameters.hpp:228
Scalar tolerance_energy_balance_relaxed_
Relaxed energy balance tolerance (can be used when iter >= min_strict_mb_iter_).
Definition: BlackoilModelParameters.hpp:165
int min_strict_mb_iter_
Minimum number of Newton iterations before we can use relaxed MB convergence criterion.
Definition: BlackoilModelParameters.hpp:222
int min_strict_cnv_iter_
Minimum number of Newton iterations before we can use relaxed CNV convergence criterion.
Definition: BlackoilModelParameters.hpp:219
bool use_update_stabilization_
Try to detect oscillation or stagnation.
Definition: BlackoilModelParameters.hpp:231
int nldd_num_initial_newton_iter_
Definition: BlackoilModelParameters.hpp:279
Scalar tolerance_cnv_energy_
Local energy convergence tolerance (max of local energy errors).
Definition: BlackoilModelParameters.hpp:171
Scalar max_residual_allowed_
Absolute max limit for residuals.
Definition: BlackoilModelParameters.hpp:154
Scalar tolerance_cnv_energy_relaxed_
Relaxed local energy convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViola...
Definition: BlackoilModelParameters.hpp:173
Scalar tolerance_mb_
Relative mass balance tolerance (total mass balance error).
Definition: BlackoilModelParameters.hpp:159
Scalar tolerance_cnv_
Local convergence tolerance (max of local saturation errors).
Definition: BlackoilModelParameters.hpp:167
Scalar tolerance_cnv_relaxed_
Relaxed local convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViolatedPV <...
Definition: BlackoilModelParameters.hpp:169
Scalar relaxed_max_pv_fraction_
Definition: BlackoilModelParameters.hpp:157
std::string nonlinear_solver_
Nonlinear solver type: newton or nldd.
Definition: BlackoilModelParameters.hpp:270
Definition: AquiferGridUtils.hpp:35
Definition: BlackoilPhases.hpp:46
int num_phases
Definition: BlackoilPhases.hpp:54
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
Definition: blackoilproperties.hh:78
Definition: fvbaseproperties.hh:53
Definition: ISTLSolver.hpp:63
Definition: BlackoilModel.hpp:80
std::tuple< FlowBaseProblemBlackoil, BlackOilModel > InheritsFrom
Definition: BlackoilModel.hpp:80
Specify whether to use volumetric residuals or not.
Definition: fvbaseproperties.hh:237
Definition: FlowBaseProblemProperties.hpp:82
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:34
double linear_solve_time
Definition: SimulatorReport.hpp:43
double assemble_time
Definition: SimulatorReport.hpp:39
bool converged
Definition: SimulatorReport.hpp:54
double pre_post_time
Definition: SimulatorReport.hpp:40
double linear_solve_setup_time
Definition: SimulatorReport.hpp:42
unsigned int total_newton_iterations
Definition: SimulatorReport.hpp:49
double update_time
Definition: SimulatorReport.hpp:44
unsigned int total_linearizations
Definition: SimulatorReport.hpp:48
unsigned int total_linear_iterations
Definition: SimulatorReport.hpp:50