24#ifndef OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED
25#define OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED
27#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED
29#include <opm/simulators/flow/BlackoilModel.hpp>
32#include <dune/common/timer.hh>
34#include <opm/common/ErrorMacros.hpp>
35#include <opm/common/OpmLog/OpmLog.hpp>
37#include <opm/simulators/flow/countGlobalCells.hpp>
39#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
47#include <fmt/format.h>
51template <
class TypeTag>
57 : simulator_(simulator)
58 , grid_(simulator_.vanguard().grid())
63 , current_relaxation_(1.0)
64 , dx_old_(simulator_.model().numGridDof())
65 , conv_monitor_(param_.monitor_params_)
69 convergence_reports_.reserve(300);
73 OpmLog::info(
"Using Non-Linear Domain Decomposition solver (nldd).");
75 nlddSolver_ = std::make_unique<BlackoilModelNldd<TypeTag>>(*this);
78 OpmLog::info(
"Using Newton nonlinear solver.");
81 OPM_THROW(std::runtime_error,
"Unknown nonlinear solver option: " +
86template <
class TypeTag>
96 if (grid_.comm().size() > 1 && lastStepFailed != grid_.comm().min(lastStepFailed)) {
98 fmt::format(
"Misalignment of the parallel simulation run in prepareStep "
99 "- the previous step succeeded on rank {} but failed on the "
100 "other ranks.", grid_.comm().rank()));
102 if (lastStepFailed) {
103 simulator_.model().updateFailed();
106 simulator_.model().advanceTimeLevel();
115 simulator_.model().newtonMethod().setIterationIndex(0);
117 simulator_.problem().beginTimeStep();
119 unsigned numDof = simulator_.model().numGridDof();
120 wasSwitched_.resize(numDof);
121 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
123 if (param_.update_equations_scaling_) {
124 OpmLog::error(
"Equation scaling not supported");
129 nlddSolver_->prepareStep();
132 report.pre_post_time +=
perfTimer.stop();
136 if (FluidSystem::phaseIsActive(
phaseIdx)) {
137 const unsigned sIdx = FluidSystem::solventComponentIndex(
phaseIdx);
138 return Indices::canonicalToActiveComponentIndex(
sIdx);
143 const auto& schedule = simulator_.vanguard().schedule();
144 auto&
rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
145 rst_conv.init(simulator_.vanguard().globalNumCells(),
147 {getIdx(FluidSystem::oilPhaseIdx),
148 getIdx(FluidSystem::gasPhaseIdx),
149 getIdx(FluidSystem::waterPhaseIdx),
157template <
class TypeTag>
171 report.total_linearizations = 1;
175 report += assembleReservoir(timer,
iteration);
176 report.assemble_time +=
perfTimer.stop();
179 report.assemble_time +=
perfTimer.stop();
180 failureReport_ += report;
192 ConvergenceReport::Severity severity =
convrep.severityOfWorstFailure();
193 convergence_reports_.back().report.push_back(std::move(
convrep));
196 if (severity == ConvergenceReport::Severity::NotANumber) {
197 failureReport_ += report;
199 }
else if (severity == ConvergenceReport::Severity::TooLarge) {
200 failureReport_ += report;
202 }
else if (severity == ConvergenceReport::Severity::ConvergenceMonitorFailure) {
203 failureReport_ += report;
206 "Total penalty count exceeded cut-off-limit of {}",
207 param_.monitor_params_.cutoff_
215template <
class TypeTag>
216template <
class NonlinearSolverType>
226 residual_norms_history_.clear();
227 conv_monitor_.reset();
228 current_relaxation_ = 1.0;
231 convergence_reports_.back().report.reserve(11);
235 if ((this->param_.nonlinear_solver_ !=
"nldd") ||
248 auto&
rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
249 rst_conv.update(simulator_.model().linearizer().residual());
254template <
class TypeTag>
255template <
class NonlinearSolverType>
266 this->initialLinearization(report,
273 if (!report.converged) {
276 report.total_newton_iterations = 1;
279 unsigned nc = simulator_.model().numGridDof();
283 linear_solve_setup_time_ = 0.0;
288 wellModel().linearize(simulator().model().linearizer().jacobian(),
289 simulator().model().linearizer().residual());
292 solveJacobianSystem(x);
294 report.linear_solve_setup_time += linear_solve_setup_time_;
295 report.linear_solve_time +=
perfTimer.stop();
296 report.total_linear_iterations += linearIterationsLastSolve();
299 report.linear_solve_setup_time += linear_solve_setup_time_;
300 report.linear_solve_time +=
perfTimer.stop();
301 report.total_linear_iterations += linearIterationsLastSolve();
303 failureReport_ += report;
313 wellModel().postSolve(x);
315 if (param_.use_update_stabilization_) {
320 residual_norms_history_.size() - 1,
325 current_relaxation_ = std::max(current_relaxation_,
327 if (terminalOutputEnabled()) {
328 std::string
msg =
" Oscillating behavior detected: Relaxation set to "
329 + std::to_string(current_relaxation_);
347template <
class TypeTag>
355 simulator_.problem().endTimeStep();
356 report.pre_post_time +=
perfTimer.stop();
360template <
class TypeTag>
367 simulator_.model().newtonMethod().setIterationIndex(
iterationIdx);
368 simulator_.problem().beginIteration();
369 simulator_.model().linearizer().linearizeDomain();
370 simulator_.problem().endIteration();
371 return wellModel().lastReport();
374template <
class TypeTag>
375typename BlackoilModel<TypeTag>::Scalar
382 const auto&
elemMapper = simulator_.model().elementMapper();
383 const auto& gridView = simulator_.gridView();
393 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
394 FluidSystem::numActivePhases() > 1 &&
395 priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
401 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
402 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
403 priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
405 assert(Indices::compositionSwitchIdx >= 0);
410 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
427 if (FluidSystem::numActivePhases() > 1) {
428 if (
priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
434 if (
priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
436 assert(Indices::compositionSwitchIdx >= 0 );
442 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
461template <
class TypeTag>
466 auto& jacobian = simulator_.model().linearizer().jacobian().istlMatrix();
467 auto& residual = simulator_.model().linearizer().residual();
468 auto&
linSolver = simulator_.model().newtonMethod().linearSolver();
472 if (terminal_output_) {
473 OpmLog::debug(
"\nRunning speed test for comparing available linear solvers.");
482 for (
int solver = 0; solver <
numSolvers; ++solver) {
494 if (terminal_output_) {
495 OpmLog::debug(fmt::format(
"Solver time {}: {}", solver,
times[solver]));
513 linear_solve_setup_time_ =
perfTimer.stop();
523template <
class TypeTag>
529 auto& newtonMethod = simulator_.model().newtonMethod();
530 SolutionVector& solution = simulator_.model().solution(0);
532 newtonMethod.update_(solution,
542 simulator_.model().invalidateAndUpdateIntensiveQuantities(0);
546template <
class TypeTag>
547std::tuple<typename BlackoilModel<TypeTag>::Scalar,
548 typename BlackoilModel<TypeTag>::Scalar>
553 std::vector< Scalar >&
R_sum,
555 std::vector< Scalar >&
B_avg)
562 if (comm.size() > 1) {
606template <
class TypeTag>
607std::pair<typename BlackoilModel<TypeTag>::Scalar,
608 typename BlackoilModel<TypeTag>::Scalar>
612 std::vector<Scalar>&
B_avg,
618 const auto& model = simulator_.model();
619 const auto& problem = simulator_.problem();
621 const auto& residual = simulator_.model().linearizer().residual();
623 ElementContext elemCtx(simulator_);
624 const auto& gridView = simulator().gridView();
626 OPM_BEGIN_PARALLEL_TRY_CATCH();
627 for (
const auto&
elem :
elements(gridView, Dune::Partitions::interior)) {
628 elemCtx.updatePrimaryStencil(
elem);
629 elemCtx.updatePrimaryIntensiveQuantities(0);
631 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
632 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
633 const auto& fs = intQuants.fluidState();
647 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModel::localConvergenceData() failed: ", grid_.comm());
651 for (
int i = 0; i <
bSize; ++i) {
652 B_avg[i] /= Scalar(global_nc_);
658template <
class TypeTag>
659std::pair<std::vector<double>, std::vector<int>>
668 constexpr auto numPvGroups = std::vector<double>::size_type{3};
670 auto cnvPvSplit = std::pair<std::vector<double>, std::vector<int>> {
671 std::piecewise_construct,
679 std::inner_product(residual.begin(), residual.end(),
680 B_avg.begin(), Scalar{0},
681 [](
const Scalar
m,
const auto& x)
684 return std::max(m, abs(x));
685 }, std::multiplies<>{});
690 const auto& model = this->simulator().model();
691 const auto& problem = this->simulator().problem();
692 const auto& residual = model.linearizer().residual();
693 const auto& gridView = this->simulator().gridView();
697 ElementContext elemCtx(this->simulator());
699 OPM_BEGIN_PARALLEL_TRY_CATCH();
700 for (
const auto&
elem :
elements(gridView, Dune::Partitions::interior)) {
706 elemCtx.updatePrimaryStencil(
elem);
708 const unsigned cell_idx = elemCtx.globalSpaceIndex(0, 0);
714 const auto ix = (
maxCnv > this->param_.tolerance_cnv_)
715 + (
maxCnv > this->param_.tolerance_cnv_relaxed_);
721 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilModel::characteriseCnvPvSplit() failed: ",
730template <
class TypeTag>
735 this->param_.tolerance_mb_ =
tuning.XXXMBE;
737 if (terminal_output_) {
738 OpmLog::debug(fmt::format(
"Setting BlackoilModel mass "
739 "balance limit (XXXMBE) to {:.2e}",
744template <
class TypeTag>
746BlackoilModel<TypeTag>::
747getReservoirConvergence(
const double reportTime,
751 std::vector<Scalar>&
B_avg,
755 using Vector = std::vector<Scalar>;
757 ConvergenceReport report{reportTime};
770 this->convergenceReduction(this->grid_.comm(),
775 report.setCnvPoreVolSplit(this->characteriseCnvPvSplit(
B_avg,
dt),
786 this->param_.min_strict_mb_iter_ < 0 &&
iteration == maxIter;
789 (this->param_.min_strict_mb_iter_ >= 0 &&
790 iteration >= this->param_.min_strict_mb_iter_);
799 this->param_.min_strict_cnv_iter_ < 0 &&
iteration == maxIter;
802 this->param_.min_strict_cnv_iter_ >= 0 &&
803 iteration >= this->param_.min_strict_cnv_iter_;
811 const auto& cnvPvSplit = report.cnvPvSplit().first;
815 return static_cast<Scalar
>(cnvPvSplit[1] + cnvPvSplit[2]) <
816 this->param_.relaxed_max_pv_fraction_ *
eligible;
824 this->terminal_output_)
826 std::string message =
827 "Number of newton iterations reached its maximum "
828 "try to continue with relaxed tolerances:";
831 message += fmt::format(
" MB: {:.1e}", param_.tolerance_mb_relaxed_);
835 message += fmt::format(
" CNV: {:.1e}", param_.tolerance_cnv_relaxed_);
838 OpmLog::debug(message);
844 const auto tol_eb =
use_relaxed_mb ? param_.tolerance_energy_balance_relaxed_ : param_.tolerance_energy_balance_;
856 using CR = ConvergenceReport;
858 const Scalar
res[2] = {
862 const CR::ReservoirFailure::Type
types[2] = {
863 CR::ReservoirFailure::Type::MassBalance,
864 CR::ReservoirFailure::Type::Cnv,
868 if (has_energy_ &&
compIdx == contiEnergyEqIdx) {
873 for (
int ii : {0, 1}) {
874 if (std::isnan(
res[
ii])) {
875 report.setReservoirFailed({
types[
ii], CR::Severity::NotANumber,
compIdx});
876 if (this->terminal_output_) {
877 OpmLog::debug(
"NaN residual for " + this->compNames_.name(
compIdx) +
" equation.");
880 else if (
res[
ii] > maxResidualAllowed()) {
881 report.setReservoirFailed({
types[
ii], CR::Severity::TooLarge,
compIdx});
882 if (this->terminal_output_) {
883 OpmLog::debug(
"Too large residual for " + this->compNames_.name(
compIdx) +
" equation.");
886 else if (
res[
ii] < 0.0) {
887 report.setReservoirFailed({
types[
ii], CR::Severity::Normal,
compIdx});
888 if (this->terminal_output_) {
889 OpmLog::debug(
"Negative residual for " + this->compNames_.name(
compIdx) +
" equation.");
893 report.setReservoirFailed({
types[
ii], CR::Severity::Normal,
compIdx});
901 if (this->terminal_output_) {
904 std::string
msg =
"Iter";
920 std::ostringstream
ss;
921 const std::streamsize
oprec =
ss.precision(3);
922 const std::ios::fmtflags
oflags =
ss.setf(std::ios::scientific);
936 OpmLog::debug(
ss.str());
942template <
class TypeTag>
952 std::vector<Scalar>
B_avg(numEq, 0.0);
958 report += wellModel().getWellConvergence(
B_avg,
962 conv_monitor_.checkPenaltyCard(report,
iteration);
967template <
class TypeTag>
968std::vector<std::vector<typename BlackoilModel<TypeTag>::Scalar> >
975 std::vector<std::vector<Scalar> >
regionValues(0, std::vector<Scalar>(0,0.0));
979template <
class TypeTag>
984 return nlddSolver_ ? nlddSolver_->localAccumulatedReports()
988template <
class TypeTag>
993 if (this->nlddSolver_ !=
nullptr) {
994 this->nlddSolver_->writePartitions(
odir);
998 const auto& elementMapper = this->simulator().model().elementMapper();
999 const auto&
cartMapper = this->simulator().vanguard().cartesianIndexMapper();
1001 const auto& grid = this->simulator().vanguard().grid();
1002 const auto& comm = grid.comm();
1003 const auto nDigit = 1 +
static_cast<int>(std::floor(std::log10(comm.size())));
1005 std::ofstream
pfile {
odir / fmt::format(
"{1:0>{0}}",
nDigit, comm.rank())};
1008 pfile << comm.rank() <<
' '
1009 <<
cartMapper.cartesianIndex(elementMapper.index(cell)) <<
' '
1010 << comm.rank() <<
'\n';
1014template <
class TypeTag>
1015template<
class Flu
idState,
class Res
idual>
1017BlackoilModel<TypeTag>::
1018getMaxCoeff(
const unsigned cell_idx,
1019 const IntensiveQuantities& intQuants,
1020 const FluidState& fs,
1023 std::vector<Scalar>&
B_avg,
1024 std::vector<Scalar>&
R_sum,
1030 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
1034 const unsigned sIdx = FluidSystem::solventComponentIndex(
phaseIdx);
1035 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(
sIdx);
1048 if constexpr (has_solvent_) {
1049 B_avg[contiSolventEqIdx] +=
1050 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
1052 R_sum[contiSolventEqIdx] +=
R2;
1056 if constexpr (has_extbo_) {
1057 B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1059 R_sum[ contiZfracEqIdx ] +=
R2;
1063 if constexpr (has_polymer_) {
1064 B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1066 R_sum[contiPolymerEqIdx] +=
R2;
1070 if constexpr (has_foam_) {
1071 B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
1077 if constexpr (has_brine_) {
1078 B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1085 if constexpr (has_polymermw_) {
1086 static_assert(has_polymer_);
1088 B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1093 R_sum[contiPolymerMWEqIdx] +=
R2;
1098 if constexpr (has_energy_) {
1099 B_avg[contiEnergyEqIdx] += 1.0 / (4.182e1);
1101 R_sum[contiEnergyEqIdx] +=
R2;
1106 if constexpr (has_micp_) {
1107 B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1109 R_sum[contiMicrobialEqIdx] +=
R1;
1112 B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1114 R_sum[contiOxygenEqIdx] +=
R2;
1117 B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1122 B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1124 R_sum[contiBiofilmEqIdx] +=
R4;
1127 B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
1129 R_sum[contiCalciteEqIdx] +=
R5;
A model implementation for three-phase black oil.
Definition BlackoilModel.hpp:62
BlackoilModel(Simulator &simulator, const ModelParameters ¶m, BlackoilWellModel< TypeTag > &well_model, const bool terminal_output)
Construct the model.
Definition BlackoilModel_impl.hpp:53
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:245
SimulatorReportSingle nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver)
Called once per nonlinear iteration.
Definition BlackoilModel_impl.hpp:219
SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &, const int iterationIdx)
Assemble the residual and Jacobian of the nonlinear system.
Definition BlackoilModel_impl.hpp:363
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_impl.hpp:945
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_impl.hpp:610
SimulatorReportSingle afterStep(const SimulatorTimerInterface &)
Called once after each time step.
Definition BlackoilModel_impl.hpp:350
void updateSolution(const BVector &dx)
Apply an update to the primary variables.
Definition BlackoilModel_impl.hpp:526
long int global_nc_
The number of cells of the global grid.
Definition BlackoilModel.hpp:327
std::unique_ptr< BlackoilModelNldd< TypeTag > > nlddSolver_
Non-linear DD solver.
Definition BlackoilModel.hpp:336
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_impl.hpp:661
SimulatorReportSingle localAccumulatedReports() const
return the statistics if the nonlinearIteration() method failed
Definition BlackoilModel_impl.hpp:982
void solveJacobianSystem(BVector &x)
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition BlackoilModel_impl.hpp:464
SimulatorReportSingle prepareStep(const SimulatorTimerInterface &timer)
Called once before each time step.
Definition BlackoilModel_impl.hpp:89
Class for handling the blackoil well model.
Definition BlackoilWellModel.hpp:94
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:109
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.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
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:174
std::string nonlinear_solver_
Nonlinear solver type: newton or nldd.
Definition BlackoilModelParameters.hpp:321
Definition AquiferGridUtils.hpp:35
A struct for returning timing data from a simulator to its caller.
Definition SimulatorReport.hpp:34