169 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
170 using Grid = GetPropType<TypeTag, Properties::Grid>;
171 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
172 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
173 using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
174 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
175 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
176 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
177 using Indices = GetPropType<TypeTag, Properties::Indices>;
178 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
179 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
180 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
182 static constexpr int numEq = Indices::numEq;
183 static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
184 static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx;
185 static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
186 static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
187 static constexpr int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
188 static constexpr int contiFoamEqIdx = Indices::contiFoamEqIdx;
189 static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
190 static constexpr int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
191 static constexpr int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
192 static constexpr int contiUreaEqIdx = Indices::contiUreaEqIdx;
193 static constexpr int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
194 static constexpr int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
195 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
196 static constexpr int zFractionIdx = Indices::zFractionIdx;
197 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
198 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
199 static constexpr int temperatureIdx = Indices::temperatureIdx;
200 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
201 static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
202 static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
203 static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
204 static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
205 static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
206 static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
208 using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
209 using MatrixBlockType =
typename SparseMatrixAdapter::MatrixBlock;
210 using Mat =
typename SparseMatrixAdapter::IstlMatrix;
211 using BVector = Dune::BlockVector<VectorBlockType>;
228 const ModelParameters&
param,
230 const bool terminal_output)
231 : simulator_(simulator)
232 , grid_(simulator_.vanguard().grid())
235 , well_model_ (well_model)
236 ,
rst_conv_(simulator_.problem().eclWriter()->collectOnIORank().localIdxToGlobalIdxMapping(),
239 , current_relaxation_(1.0)
240 , dx_old_(simulator_.model().numGridDof())
244 convergence_reports_.reserve(300);
246 if (param_.nonlinear_solver_ ==
"nldd") {
247 if (!param_.matrix_add_well_contributions_) {
248 OPM_THROW(std::runtime_error,
"The --nonlinear-solver=nldd option can only be used with --matrix-add-well-contributions=true");
250 if (terminal_output) {
251 OpmLog::info(
"Using Non-Linear Domain Decomposition solver (nldd).");
253 nlddSolver_ = std::make_unique<BlackoilModelNldd<TypeTag>>(*this);
254 }
else if (param_.nonlinear_solver_ ==
"newton") {
255 if (terminal_output) {
256 OpmLog::info(
"Using Newton nonlinear solver.");
259 OPM_THROW(std::runtime_error,
"Unknown nonlinear solver option: " + param_.nonlinear_solver_);
264 bool isParallel()
const
265 {
return grid_.comm().size() > 1; }
268 const EclipseState& eclState()
const
269 {
return simulator_.vanguard().eclState(); }
277 Dune::Timer perfTimer;
281 simulator_.model().updateFailed();
283 simulator_.model().advanceTimeLevel();
292 simulator_.model().newtonMethod().setIterationIndex(0);
294 simulator_.problem().beginTimeStep();
296 unsigned numDof = simulator_.model().numGridDof();
297 wasSwitched_.resize(numDof);
298 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
300 if (param_.update_equations_scaling_) {
301 OpmLog::error(
"Equation scaling not supported");
309 report.pre_post_time += perfTimer.stop();
311 auto getIdx = [](
unsigned phaseIdx) ->
int
313 if (FluidSystem::phaseIsActive(phaseIdx)) {
314 const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
315 return Indices::canonicalToActiveComponentIndex(sIdx);
320 const auto& schedule = simulator_.vanguard().schedule();
323 {getIdx(FluidSystem::oilPhaseIdx),
324 getIdx(FluidSystem::gasPhaseIdx),
325 getIdx(FluidSystem::waterPhaseIdx),
342 Dune::Timer perfTimer;
345 report.total_linearizations = 1;
350 report.assemble_time += perfTimer.stop();
353 report.assemble_time += perfTimer.stop();
354 failureReport_ += report;
359 std::vector<double> residual_norms;
364 auto convrep =
getConvergence(timer, iteration, maxIter, residual_norms);
365 report.converged = convrep.converged() && iteration >= minIter;
366 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
367 convergence_reports_.back().report.push_back(std::move(convrep));
370 if (severity == ConvergenceReport::Severity::NotANumber) {
371 OPM_THROW(NumericalProblem,
"NaN residual found!");
372 }
else if (severity == ConvergenceReport::Severity::TooLarge) {
373 OPM_THROW_NOLOG(NumericalProblem,
"Too large residual found!");
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 current_relaxation_ = 1.0;
400 convergence_reports_.back().report.reserve(11);
404 if ((this->param_.nonlinear_solver_ !=
"nldd") ||
405 (iteration < this->param_.nldd_num_initial_newton_iter_))
407 result = this->nonlinearIterationNewton(iteration, timer, nonlinear_solver);
410 result = this->
nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver);
419 template <
class NonlinearSolverType>
422 NonlinearSolverType& nonlinear_solver)
427 Dune::Timer perfTimer;
429 this->initialLinearization(report, iteration, nonlinear_solver.minIter(), nonlinear_solver.maxIter(), timer);
432 if (!report.converged) {
435 report.total_newton_iterations = 1;
438 unsigned nc = simulator_.model().numGridDof();
442 linear_solve_setup_time_ = 0.0;
447 wellModel().linearize(simulator().model().linearizer().jacobian(),
448 simulator().model().linearizer().residual());
453 report.linear_solve_setup_time += linear_solve_setup_time_;
454 report.linear_solve_time += perfTimer.stop();
458 report.linear_solve_setup_time += linear_solve_setup_time_;
459 report.linear_solve_time += perfTimer.stop();
462 failureReport_ += report;
474 if (param_.use_update_stabilization_) {
476 bool isOscillate =
false;
477 bool isStagnate =
false;
478 nonlinear_solver.detectOscillations(residual_norms_history_, residual_norms_history_.size() - 1, isOscillate, isStagnate);
480 current_relaxation_ -= nonlinear_solver.relaxIncrement();
481 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
483 std::string msg =
" Oscillating behavior detected: Relaxation set to "
484 + std::to_string(current_relaxation_);
488 nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
496 report.update_time += perfTimer.stop();
509 Dune::Timer perfTimer;
511 simulator_.problem().endTimeStep();
513 report.pre_post_time += perfTimer.stop();
519 const int iterationIdx)
522 simulator_.model().newtonMethod().setIterationIndex(iterationIdx);
523 simulator_.problem().beginIteration();
524 simulator_.model().linearizer().linearizeDomain();
525 simulator_.problem().endIteration();
530 double relativeChange()
const
532 Scalar resultDelta = 0.0;
533 Scalar resultDenom = 0.0;
535 const auto& elemMapper = simulator_.model().elementMapper();
536 const auto& gridView = simulator_.gridView();
537 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
538 unsigned globalElemIdx = elemMapper.index(elem);
539 const auto& priVarsNew = simulator_.model().solution(0)[globalElemIdx];
542 pressureNew = priVarsNew[Indices::pressureSwitchIdx];
544 Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
545 Scalar oilSaturationNew = 1.0;
546 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
547 FluidSystem::numActivePhases() > 1 &&
548 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
549 saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
550 oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
553 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
554 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
555 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
556 assert(Indices::compositionSwitchIdx >= 0 );
557 saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
558 oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
561 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
562 saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
565 const auto& priVarsOld = simulator_.model().solution(1)[globalElemIdx];
568 pressureOld = priVarsOld[Indices::pressureSwitchIdx];
570 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
571 Scalar oilSaturationOld = 1.0;
574 Scalar tmp = pressureNew - pressureOld;
575 resultDelta += tmp*tmp;
576 resultDenom += pressureNew*pressureNew;
578 if (FluidSystem::numActivePhases() > 1) {
579 if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
580 saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSwitchIdx];
581 oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
584 if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
586 assert(Indices::compositionSwitchIdx >= 0 );
587 saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
588 oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
591 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
592 saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
594 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
595 Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
596 resultDelta += tmpSat*tmpSat;
597 resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
598 assert(std::isfinite(resultDelta));
599 assert(std::isfinite(resultDenom));
604 resultDelta = gridView.comm().sum(resultDelta);
605 resultDenom = gridView.comm().sum(resultDenom);
607 if (resultDenom > 0.0)
608 return resultDelta/resultDenom;
616 return simulator_.model().newtonMethod().linearSolver().iterations ();
621 double& linearSolveSetupTime()
623 return linear_solve_setup_time_;
631 auto& jacobian = simulator_.model().linearizer().jacobian().istlMatrix();
632 auto& residual = simulator_.model().linearizer().residual();
633 auto& linSolver = simulator_.model().newtonMethod().linearSolver();
635 const int numSolvers = linSolver.numAvailableSolvers();
636 if ((numSolvers > 1) && (linSolver.getSolveCount() % 100 == 0)) {
639 OpmLog::debug(
"\nRunning speed test for comparing available linear solvers.");
642 Dune::Timer perfTimer;
643 std::vector<double> times(numSolvers);
644 std::vector<double> setupTimes(numSolvers);
647 std::vector<BVector> x_trial(numSolvers, x);
648 for (
int solver = 0; solver < numSolvers; ++solver) {
650 linSolver.setActiveSolver(solver);
652 linSolver.prepare(jacobian, residual);
653 setupTimes[solver] = perfTimer.stop();
655 linSolver.setResidual(residual);
657 linSolver.solve(x_trial[solver]);
658 times[solver] = perfTimer.stop();
661 OpmLog::debug(fmt::format(
"Solver time {}: {}", solver, times[solver]));
665 int fastest_solver = std::min_element(times.begin(), times.end()) - times.begin();
667 grid_.comm().broadcast(&fastest_solver, 1, 0);
668 linear_solve_setup_time_ = setupTimes[fastest_solver];
669 x = x_trial[fastest_solver];
670 linSolver.setActiveSolver(fastest_solver);
675 Dune::Timer perfTimer;
677 linSolver.prepare(jacobian, residual);
678 linear_solve_setup_time_ = perfTimer.stop();
679 linSolver.setResidual(residual);
693 auto& newtonMethod = simulator_.model().newtonMethod();
694 SolutionVector& solution = simulator_.model().solution(0);
696 newtonMethod.update_(solution,
705 OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
706 simulator_.model().invalidateAndUpdateIntensiveQuantities(0);
707 simulator_.problem().eclWriter()->mutableOutputModule().invalidateLocalData();
717 std::tuple<double,double> convergenceReduction(Parallel::Communication comm,
718 const double pvSumLocal,
719 const double numAquiferPvSumLocal,
720 std::vector< Scalar >& R_sum,
721 std::vector< Scalar >& maxCoeff,
722 std::vector< Scalar >& B_avg)
724 OPM_TIMEBLOCK(convergenceReduction);
726 double pvSum = pvSumLocal;
727 double numAquiferPvSum = numAquiferPvSumLocal;
729 if( comm.size() > 1 )
732 std::vector< Scalar > sumBuffer;
733 std::vector< Scalar > maxBuffer;
734 const int numComp = B_avg.size();
735 sumBuffer.reserve( 2*numComp + 2 );
736 maxBuffer.reserve( numComp );
737 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
739 sumBuffer.push_back( B_avg[ compIdx ] );
740 sumBuffer.push_back( R_sum[ compIdx ] );
741 maxBuffer.push_back( maxCoeff[ compIdx ] );
745 sumBuffer.push_back( pvSum );
746 sumBuffer.push_back( numAquiferPvSum );
749 comm.sum( sumBuffer.data(), sumBuffer.size() );
752 comm.max( maxBuffer.data(), maxBuffer.size() );
755 for(
int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
757 B_avg[ compIdx ] = sumBuffer[ buffIdx ];
760 R_sum[ compIdx ] = sumBuffer[ buffIdx ];
763 for(
int compIdx = 0; compIdx < numComp; ++compIdx )
765 maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
769 pvSum = sumBuffer[sumBuffer.size()-2];
770 numAquiferPvSum = sumBuffer.back();
774 return {pvSum, numAquiferPvSum};
781 std::vector<Scalar>& maxCoeff,
782 std::vector<Scalar>& B_avg,
783 std::vector<int>& maxCoeffCell)
786 double pvSumLocal = 0.0;
787 double numAquiferPvSumLocal = 0.0;
788 const auto& model = simulator_.model();
789 const auto& problem = simulator_.problem();
791 const auto& residual = simulator_.model().linearizer().residual();
793 ElementContext elemCtx(simulator_);
794 const auto& gridView = simulator().gridView();
796 OPM_BEGIN_PARALLEL_TRY_CATCH();
797 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
798 elemCtx.updatePrimaryStencil(elem);
799 elemCtx.updatePrimaryIntensiveQuantities(0);
801 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
802 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
803 const auto& fs = intQuants.fluidState();
805 const auto pvValue = problem.referencePorosity(cell_idx, 0) *
806 model.dofTotalVolume(cell_idx);
807 pvSumLocal += pvValue;
809 if (isNumericalAquiferCell(elem))
811 numAquiferPvSumLocal += pvValue;
814 this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue,
815 B_avg, R_sum, maxCoeff, maxCoeffCell);
818 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModel::localConvergenceData() failed: ", grid_.comm());
821 const int bSize = B_avg.size();
822 for (
int i = 0; i<bSize; ++i )
827 return {pvSumLocal, numAquiferPvSumLocal};
837 const auto& model = simulator_.model();
838 const auto& problem = simulator_.problem();
839 const auto& residual = simulator_.model().linearizer().residual();
840 const auto& gridView = simulator().gridView();
841 ElementContext elemCtx(simulator_);
844 OPM_BEGIN_PARALLEL_TRY_CATCH();
845 for (
const auto& elem : elements(gridView, Dune::Partitions::interiorBorder))
848 if (isNumericalAquiferCell(elem))
852 elemCtx.updatePrimaryStencil(elem);
854 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
855 const double pvValue = problem.referencePorosity(cell_idx, 0) * model.dofTotalVolume(cell_idx);
856 const auto& cellResidual = residual[cell_idx];
857 bool cnvViolated =
false;
859 for (
unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx)
862 Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
863 cnvViolated = cnvViolated || (abs(CNV) > param_.tolerance_cnv_);
872 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModel::ComputeCnvError() failed: ", grid_.comm());
874 return grid_.comm().sum(errorPV);
878 void updateTUNING(
const Tuning& tuning) {
879 param_.tolerance_mb_ = tuning.XXXMBE;
881 OpmLog::debug(fmt::format(
"Setting BlackoilModel mass balance limit (XXXMBE) to {:.2e}", tuning.XXXMBE));
886 ConvergenceReport getReservoirConvergence(
const double reportTime,
890 std::vector<Scalar>& B_avg,
891 std::vector<Scalar>& residual_norms)
893 OPM_TIMEBLOCK(getReservoirConvergence);
894 using Vector = std::vector<Scalar>;
896 const int numComp = numEq;
897 Vector R_sum(numComp, 0.0 );
898 Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
899 std::vector<int> maxCoeffCell(numComp, -1);
900 const auto [ pvSumLocal, numAquiferPvSumLocal] =
localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);
903 const auto [ pvSum, numAquiferPvSum ] =
904 convergenceReduction(grid_.comm(), pvSumLocal,
905 numAquiferPvSumLocal,
906 R_sum, maxCoeff, B_avg);
909 cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
917 const bool relax_final_iteration_mb = (param_.min_strict_mb_iter_ < 0) && (iteration == maxIter);
918 const bool use_relaxed_mb = relax_final_iteration_mb || (param_.min_strict_mb_iter_ >= 0 && iteration >= param_.min_strict_mb_iter_);
925 const bool relax_final_iteration_cnv = (param_.min_strict_cnv_iter_ < 0) && (iteration == maxIter);
926 const bool relax_iter_cnv = param_.min_strict_mb_iter_ >= 0 && iteration >= param_.min_strict_mb_iter_;
927 const bool relax_pv_fraction_cnv = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_;
928 const bool use_relaxed_cnv = relax_final_iteration_cnv || relax_pv_fraction_cnv || relax_iter_cnv;
930 if (relax_final_iteration_mb || relax_final_iteration_cnv) {
932 std::string message =
"Number of newton iterations reached its maximum try to continue with relaxed tolerances:";
933 if (relax_final_iteration_mb)
934 message += fmt::format(
" MB: {:.1e}", param_.tolerance_mb_relaxed_);
935 if (relax_final_iteration_cnv)
936 message += fmt::format(
" CNV: {:.1e}", param_.tolerance_cnv_relaxed_);
938 OpmLog::debug(message);
941 const double tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
942 const double tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
945 std::vector<Scalar> CNV(numComp);
946 std::vector<Scalar> mass_balance_residual(numComp);
947 for (
int compIdx = 0; compIdx < numComp; ++compIdx )
949 CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
950 mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
951 residual_norms.push_back(CNV[compIdx]);
955 ConvergenceReport report{reportTime};
956 using CR = ConvergenceReport;
957 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
958 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
959 CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
960 CR::ReservoirFailure::Type::Cnv };
961 double tol[2] = { tol_mb, tol_cnv };
962 for (
int ii : {0, 1}) {
963 if (std::isnan(res[ii])) {
964 report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
966 OpmLog::debug(
"NaN residual for " + this->compNames_.name(compIdx) +
" equation.");
968 }
else if (res[ii] > maxResidualAllowed()) {
969 report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
971 OpmLog::debug(
"Too large residual for " + this->compNames_.name(compIdx) +
" equation.");
973 }
else if (res[ii] < 0.0) {
974 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
976 OpmLog::debug(
"Negative residual for " + this->compNames_.name(compIdx) +
" equation.");
978 }
else if (res[ii] > tol[ii]) {
979 report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
981 report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
989 if (iteration == 0) {
990 std::string msg =
"Iter";
991 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
993 msg += this->compNames_.name(compIdx)[0];
996 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
998 msg += this->compNames_.name(compIdx)[0];
1003 std::ostringstream ss;
1004 const std::streamsize oprec = ss.precision(3);
1005 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1006 ss << std::setw(4) << iteration;
1007 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1008 ss << std::setw(11) << mass_balance_residual[compIdx];
1010 for (
int compIdx = 0; compIdx < numComp; ++compIdx) {
1011 ss << std::setw(11) << CNV[compIdx];
1013 ss.precision(oprec);
1015 OpmLog::debug(ss.str());
1028 const int iteration,
1030 std::vector<double>& residual_norms)
1034 std::vector<Scalar> B_avg(numEq, 0.0);
1037 iteration, maxIter, B_avg, residual_norms);
1039 OPM_TIMEBLOCK(getWellConvergence);
1040 report +=
wellModel().getWellConvergence(B_avg, report.converged());
1049 return phaseUsage_.num_phases;
1054 std::vector<std::vector<double> >
1061 std::vector<std::vector<double> >
1067 std::vector<std::vector<double> > regionValues(0, std::vector<double>(0,0.0));
1068 return regionValues;
1071 const Simulator& simulator()
const
1072 {
return simulator_; }
1074 Simulator& simulator()
1075 {
return simulator_; }
1079 {
return failureReport_; }
1088 const std::vector<StepReport>& stepReports()
const
1090 return convergence_reports_;
1093 void writePartitions(
const std::filesystem::path& odir)
const
1100 const auto& elementMapper = this->simulator().model().elementMapper();
1101 const auto& cartMapper = this->simulator().vanguard().cartesianIndexMapper();
1103 const auto& grid = this->simulator().vanguard().grid();
1104 const auto& comm = grid.comm();
1105 const auto nDigit = 1 +
static_cast<int>(std::floor(std::log10(comm.size())));
1107 std::ofstream pfile { odir / fmt::format(
"{1:0>{0}}", nDigit, comm.rank()) };
1109 for (
const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
1110 pfile << comm.rank() <<
' '
1111 << cartMapper.cartesianIndex(elementMapper.index(cell)) <<
' '
1112 << comm.rank() <<
'\n';
1116 const std::vector<std::vector<int>>& getConvCells()
const
1122 Simulator& simulator_;
1124 const PhaseUsage phaseUsage_;
1125 static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
1126 static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
1127 static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
1128 static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
1129 static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnableEnergy>();
1130 static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
1131 static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
1132 static constexpr bool has_micp_ = getPropValue<TypeTag, Properties::EnableMICP>();
1134 ModelParameters param_;
1135 SimulatorReportSingle failureReport_;
1138 BlackoilWellModel<TypeTag>& well_model_;
1147 std::vector<std::vector<double>> residual_norms_history_;
1148 double current_relaxation_;
1151 std::vector<StepReport> convergence_reports_;
1162 wellModel()
const {
return well_model_; }
1164 void beginReportStep()
1166 simulator_.problem().beginEpisode();
1169 void endReportStep()
1171 simulator_.problem().endEpisode();
1174 template<
class Flu
idState,
class Res
idual>
1175 void getMaxCoeff(
const unsigned cell_idx,
1176 const IntensiveQuantities& intQuants,
1177 const FluidState& fs,
1178 const Residual& modelResid,
1179 const Scalar pvValue,
1180 std::vector<Scalar>& B_avg,
1181 std::vector<Scalar>& R_sum,
1182 std::vector<Scalar>& maxCoeff,
1183 std::vector<int>& maxCoeffCell)
1185 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1187 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1191 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1193 B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
1194 const auto R2 = modelResid[cell_idx][compIdx];
1196 R_sum[compIdx] += R2;
1197 const double Rval = std::abs(R2) / pvValue;
1198 if (Rval > maxCoeff[compIdx]) {
1199 maxCoeff[compIdx] = Rval;
1200 maxCoeffCell[compIdx] = cell_idx;
1204 if constexpr (has_solvent_) {
1205 B_avg[contiSolventEqIdx] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1206 const auto R2 = modelResid[cell_idx][contiSolventEqIdx];
1207 R_sum[contiSolventEqIdx] += R2;
1208 maxCoeff[contiSolventEqIdx] = std::max(maxCoeff[contiSolventEqIdx],
1209 std::abs(R2) / pvValue);
1211 if constexpr (has_extbo_) {
1212 B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1213 const auto R2 = modelResid[cell_idx][contiZfracEqIdx];
1214 R_sum[ contiZfracEqIdx ] += R2;
1215 maxCoeff[contiZfracEqIdx] = std::max(maxCoeff[contiZfracEqIdx],
1216 std::abs(R2) / pvValue);
1218 if constexpr (has_polymer_) {
1219 B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1220 const auto R2 = modelResid[cell_idx][contiPolymerEqIdx];
1221 R_sum[contiPolymerEqIdx] += R2;
1222 maxCoeff[contiPolymerEqIdx] = std::max(maxCoeff[contiPolymerEqIdx],
1223 std::abs(R2) / pvValue);
1225 if constexpr (has_foam_) {
1226 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1227 const auto R2 = modelResid[cell_idx][contiFoamEqIdx];
1228 R_sum[contiFoamEqIdx] += R2;
1229 maxCoeff[contiFoamEqIdx] = std::max(maxCoeff[contiFoamEqIdx],
1230 std::abs(R2) / pvValue);
1232 if constexpr (has_brine_) {
1233 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1234 const auto R2 = modelResid[cell_idx][contiBrineEqIdx];
1235 R_sum[contiBrineEqIdx] += R2;
1236 maxCoeff[contiBrineEqIdx] = std::max(maxCoeff[contiBrineEqIdx],
1237 std::abs(R2) / pvValue);
1240 if constexpr (has_polymermw_) {
1241 static_assert(has_polymer_);
1243 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1247 const auto R2 = modelResid[cell_idx][contiPolymerMWEqIdx] / 100.;
1248 R_sum[contiPolymerMWEqIdx] += R2;
1249 maxCoeff[contiPolymerMWEqIdx] = std::max(maxCoeff[contiPolymerMWEqIdx],
1250 std::abs(R2) / pvValue);
1253 if constexpr (has_energy_) {
1254 B_avg[contiEnergyEqIdx] += 1.0 / (4.182e1);
1255 const auto R2 = modelResid[cell_idx][contiEnergyEqIdx];
1256 R_sum[contiEnergyEqIdx] += R2;
1257 maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
1258 std::abs(R2) / pvValue);
1261 if constexpr (has_micp_) {
1262 B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1263 const auto R1 = modelResid[cell_idx][contiMicrobialEqIdx];
1264 R_sum[contiMicrobialEqIdx] += R1;
1265 maxCoeff[contiMicrobialEqIdx] = std::max(maxCoeff[contiMicrobialEqIdx],
1266 std::abs(R1) / pvValue);
1267 B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1268 const auto R2 = modelResid[cell_idx][contiOxygenEqIdx];
1269 R_sum[contiOxygenEqIdx] += R2;
1270 maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
1271 std::abs(R2) / pvValue);
1272 B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1273 const auto R3 = modelResid[cell_idx][contiUreaEqIdx];
1274 R_sum[contiUreaEqIdx] += R3;
1275 maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
1276 std::abs(R3) / pvValue);
1277 B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1278 const auto R4 = modelResid[cell_idx][contiBiofilmEqIdx];
1279 R_sum[contiBiofilmEqIdx] += R4;
1280 maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
1281 std::abs(R4) / pvValue);
1282 B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1283 const auto R5 = modelResid[cell_idx][contiCalciteEqIdx];
1284 R_sum[contiCalciteEqIdx] += R5;
1285 maxCoeff[contiCalciteEqIdx] = std::max(maxCoeff[contiCalciteEqIdx],
1286 std::abs(R5) / pvValue);
1303 double dpMaxRel()
const {
return param_.dp_max_rel_; }
1304 double dsMax()
const {
return param_.ds_max_; }
1305 double drMaxRel()
const {
return param_.dr_max_rel_; }
1306 double maxResidualAllowed()
const {
return param_.max_residual_allowed_; }
1307 double linear_solve_setup_time_;
1310 std::vector<bool> wasSwitched_;