OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-25 23:20:12 -05:00
Code Issues Packages Projects Releases Wiki Activity
Files
86bcfbd5fa6dcedbe9828d426801262269d438b3
opm-simulators/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp
T

1394 lines
54 KiB
C++
Raw Blame History

/*
Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2014, 2015 Statoil ASA.
Copyright 2015 NTNU
Copyright 2015, 2016, 2017 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_IMPL_HEADER_INCLUDED
#define OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_IMPL_HEADER_INCLUDED
#ifndef OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_HEADER_INCLUDED
#include <config.h>
#include <opm/simulators/flow/NonlinearSystemBlackOilReservoir.hpp>
#endif
#include <dune/common/timer.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/simulators/flow/countGlobalCells.hpp>
#include <algorithm>
#include <cmath>
#include <filesystem>
#include <functional>
#include <iomanip>
#include <limits>
#include <memory>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <string>
#include <string_view>
#include <tuple>
#include <utility>
#include <vector>
#include <fmt/format.h>
namespace {
template <typename TypeTag>
std::string_view
make_string(const typename Opm::NonlinearSystemBlackOilReservoir<TypeTag>::DebugFlags f)
{
using F = typename Opm::NonlinearSystemBlackOilReservoir<TypeTag>::DebugFlags;
switch (f) {
case F::STRICT: return "Strict";
case F::RELAXED: return "Relaxed";
case F::TUNINGDP: return "TuningDP";
}
return "< ??? >";
}
} // Anonymous namespace
namespace Opm {
template <class TypeTag>
NonlinearSystemBlackOilReservoir<TypeTag>::
NonlinearSystemBlackOilReservoir(Simulator& simulator,
const ModelParameters& param,
BlackoilWellModel<TypeTag>& well_model,
const bool terminal_output)
: simulator_(simulator)
, grid_(simulator_.vanguard().grid())
, param_( param )
, well_model_ (well_model)
, terminal_output_ (terminal_output)
, current_relaxation_(1.0)
, dx_old_(simulator_.model().numGridDof())
, conv_monitor_(param_.monitor_params_)
{
// compute global sum of number of cells
global_nc_ = detail::countGlobalCells(grid_);
convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
// TODO: remember to fix!
if (param_.nonlinear_solver_ == "nldd") {
if (terminal_output) {
OpmLog::info("Using Non-Linear Domain Decomposition solver (nldd).");
}
nlddSolver_ = std::make_unique<NonlinearSystemNldd<TypeTag>>(*this);
} else if (param_.nonlinear_solver_ == "newton") {
if (terminal_output) {
OpmLog::info("Using Newton nonlinear solver.");
}
} else {
OPM_THROW(std::runtime_error, "Unknown nonlinear solver option: " +
param_.nonlinear_solver_);
}
}
template <class TypeTag>
SimulatorReportSingle
NonlinearSystemBlackOilReservoir<TypeTag>::
prepareStep(const SimulatorTimerInterface& timer)
{
OPM_TIMEFUNCTION();
SimulatorReportSingle report;
Dune::Timer perfTimer;
perfTimer.start();
// update the solution variables in the model
int lastStepFailed = timer.lastStepFailed();
if (grid_.comm().size() > 1 && grid_.comm().max(lastStepFailed) != grid_.comm().min(lastStepFailed)) {
OPM_THROW(std::runtime_error, "Misalignment of the parallel simulation run in prepareStep " +
"- the previous step succeeded on some ranks but failed on others.");
}
if (lastStepFailed) {
simulator_.model().updateFailed();
}
else {
simulator_.model().advanceTimeLevel();
}
// Set the timestep size and episode index for the model explicitly.
// The model needs to know the report step/episode index because of
// timing dependent data despite the fact that flow uses its own time stepper.
// (The length of the episode does not matter, though.)
simulator_.setTime(timer.simulationTimeElapsed());
simulator_.setTimeStepSize(timer.currentStepLength());
simulator_.problem().resetIterationForNewTimestep();
simulator_.problem().beginTimeStep();
unsigned numDof = simulator_.model().numGridDof();
wasSwitched_.resize(numDof);
std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
if (param_.update_equations_scaling_) {
OpmLog::error("Equation scaling not supported");
//updateEquationsScaling();
}
if (hasNlddSolver()) {
nlddSolver_->prepareStep();
}
report.pre_post_time += perfTimer.stop();
auto getIdx = [](unsigned phaseIdx) -> int
{
if (FluidSystem::phaseIsActive(phaseIdx)) {
const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
return FluidSystem::canonicalToActiveCompIdx(sIdx);
}
return -1;
};
const auto& schedule = simulator_.vanguard().schedule();
auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
rst_conv.init(simulator_.vanguard().globalNumCells(),
schedule[timer.reportStepNum()].rst_config(),
{getIdx(FluidSystem::oilPhaseIdx),
getIdx(FluidSystem::gasPhaseIdx),
getIdx(FluidSystem::waterPhaseIdx),
contiPolymerEqIdx,
contiBrineEqIdx,
contiSolventEqIdx});
return report;
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
initialLinearization(SimulatorReportSingle& report,
const int minIter,
const int maxIter,
const SimulatorTimerInterface& timer)
{
// ----------- Set up reports and timer -----------
failureReport_ = SimulatorReportSingle();
Dune::Timer perfTimer;
perfTimer.start();
report.total_linearizations = 1;
// ----------- Assemble -----------
try {
report += assembleReservoir(timer);
report.assemble_time += perfTimer.stop();
// Mark timestep initialized after assembleReservoir(), because the well model's
// assemble() also uses needsTimestepInit() to trigger prepareTimeStep().
simulator_.problem().markTimestepInitialized();
}
catch (...) {
report.assemble_time += perfTimer.stop();
failureReport_ += report;
throw; // continue throwing the stick
}
// ----------- Check if converged -----------
std::vector<Scalar> residual_norms;
perfTimer.reset();
perfTimer.start();
// the step is not considered converged until at least minIter iterations is done
{
auto convrep = getConvergence(timer, maxIter, residual_norms);
report.converged = convrep.converged() &&
simulator_.problem().iterationContext().iteration() >= minIter;
ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
convergence_reports_.back().report.push_back(std::move(convrep));
// Throw if any NaN or too large residual found.
if (severity == ConvergenceReport::Severity::NotANumber) {
failureReport_ += report;
OPM_THROW_PROBLEM(NumericalProblem, "NaN residual found!");
} else if (severity == ConvergenceReport::Severity::TooLarge) {
failureReport_ += report;
OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!");
} else if (severity == ConvergenceReport::Severity::ConvergenceMonitorFailure) {
failureReport_ += report;
OPM_THROW_PROBLEM(ConvergenceMonitorFailure,
fmt::format(
"Total penalty count exceeded cut-off-limit of {}",
param_.monitor_params_.cutoff_
));
}
}
report.update_time += perfTimer.stop();
residual_norms_history_.push_back(residual_norms);
}
template <class TypeTag>
template <class NonlinearSolverType>
SimulatorReportSingle
NonlinearSystemBlackOilReservoir<TypeTag>::
nonlinearIteration(const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver)
{
// Model-level timestep initialization (once per timestep).
// markTimestepInitialized() is called later in initialLinearization(),
// after assembleReservoir() has triggered the well model's prepareTimeStep().
if (simulator_.problem().iterationContext().needsTimestepInit()) {
residual_norms_history_.clear();
conv_monitor_.reset();
current_relaxation_ = 1.0;
dx_old_ = 0.0;
convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
convergence_reports_.back().report.reserve(11);
}
SimulatorReportSingle result;
if (this->param_.nonlinear_solver_ != "nldd")
{
result = this->nonlinearIterationNewton(timer, nonlinear_solver);
}
else {
result = this->nlddSolver_->nonlinearIterationNldd(timer, nonlinear_solver);
}
auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
rst_conv.update(simulator_.model().linearizer().residual());
simulator_.problem().advanceIteration();
return result;
}
template <class TypeTag>
template <class NonlinearSolverType>
SimulatorReportSingle
NonlinearSystemBlackOilReservoir<TypeTag>::
nonlinearIterationNewton(const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver)
{
OPM_TIMEFUNCTION();
// ----------- Set up reports and timer -----------
SimulatorReportSingle report;
Dune::Timer perfTimer;
this->initialLinearization(report,
this->param_.newton_min_iter_,
this->param_.newton_max_iter_,
timer);
// ----------- If not converged, solve linear system and do Newton update -----------
if (!report.converged) {
perfTimer.reset();
perfTimer.start();
report.total_newton_iterations = 1;
// Compute the nonlinear update.
unsigned nc = simulator_.model().numGridDof();
BVector x(nc);
// Solve the linear system.
linear_solve_setup_time_ = 0.0;
try {
// Apply the Schur complement of the well model to
// the reservoir linearized equations.
// Note that linearize may throw for MSwells.
wellModel().linearize(simulator().model().linearizer().jacobian(),
simulator().model().linearizer().residual());
// ---- Solve linear system ----
solveJacobianSystem(x);
report.linear_solve_setup_time += linear_solve_setup_time_;
report.linear_solve_time += perfTimer.stop();
report.total_linear_iterations += linearIterationsLastSolve();
}
catch (...) {
report.linear_solve_setup_time += linear_solve_setup_time_;
report.linear_solve_time += perfTimer.stop();
report.total_linear_iterations += linearIterationsLastSolve();
failureReport_ += report;
throw; // re-throw up
}
perfTimer.reset();
perfTimer.start();
// handling well state update before oscillation treatment is a decision based
// on observation to avoid some big performance degeneration under some circumstances.
// there is no theorectical explanation which way is better for sure.
wellModel().postSolve(x);
if (param_.use_update_stabilization_) {
// Stabilize the nonlinear update.
bool isOscillate = false;
bool isStagnate = false;
nonlinear_solver.detectOscillations(residual_norms_history_,
residual_norms_history_.size() - 1,
isOscillate,
isStagnate);
if (isOscillate) {
current_relaxation_ -= nonlinear_solver.relaxIncrement();
current_relaxation_ = std::max(current_relaxation_,
nonlinear_solver.relaxMax());
if (terminalOutputEnabled()) {
std::string msg = " Oscillating behavior detected: Relaxation set to "
+ std::to_string(current_relaxation_);
OpmLog::info(msg);
}
}
nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
}
// ---- Newton update ----
// Apply the update, with considering model-dependent limitations and
// chopping of the update.
updateSolution(x);
report.update_time += perfTimer.stop();
}
return report;
}
template <class TypeTag>
SimulatorReportSingle
NonlinearSystemBlackOilReservoir<TypeTag>::
assembleReservoir(const SimulatorTimerInterface& /* timer */)
{
// -------- Mass balance equations --------
simulator_.problem().beginIteration();
simulator_.model().linearizer().linearizeDomain();
simulator_.problem().endIteration();
return wellModel().lastReport();
}
template <class TypeTag>
typename NonlinearSystemBlackOilReservoir<TypeTag>::Scalar
NonlinearSystemBlackOilReservoir<TypeTag>::
relativeChange() const
{
Scalar resultDelta = 0.0;
Scalar resultDenom = 0.0;
const auto& elemMapper = simulator_.model().elementMapper();
const auto& gridView = simulator_.gridView();
for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
unsigned globalElemIdx = elemMapper.index(elem);
const auto& priVarsNew = simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
Scalar pressureNew;
pressureNew = priVarsNew[Indices::pressureSwitchIdx];
Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
Scalar oilSaturationNew = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) &&
FluidSystem::numActivePhases() > 1 &&
priVarsNew.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
{
saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSwitchIdx];
oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
priVarsNew.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
{
assert(Indices::compositionSwitchIdx >= 0);
saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
}
const auto& priVarsOld = simulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
Scalar pressureOld;
pressureOld = priVarsOld[Indices::pressureSwitchIdx];
Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
Scalar oilSaturationOld = 1.0;
// NB fix me! adding pressures changes to saturation changes does not make sense
Scalar tmp = pressureNew - pressureOld;
resultDelta += tmp*tmp;
resultDenom += pressureNew*pressureNew;
if (FluidSystem::numActivePhases() > 1) {
if (priVarsOld.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
saturationsOld[FluidSystem::waterPhaseIdx] =
priVarsOld[Indices::waterSwitchIdx];
oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
}
if (priVarsOld.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
{
assert(Indices::compositionSwitchIdx >= 0 );
saturationsOld[FluidSystem::gasPhaseIdx] =
priVarsOld[Indices::compositionSwitchIdx];
oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
}
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
resultDelta += tmpSat*tmpSat;
resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
assert(std::isfinite(resultDelta));
assert(std::isfinite(resultDenom));
}
}
}
resultDelta = gridView.comm().sum(resultDelta);
resultDenom = gridView.comm().sum(resultDenom);
return resultDenom > 0.0 ? resultDelta / resultDenom : 0.0;
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
solveJacobianSystem(BVector& x)
{
auto& jacobian = simulator_.model().linearizer().jacobian().istlMatrix();
auto& residual = simulator_.model().linearizer().residual();
auto& linSolver = simulator_.model().newtonMethod().linearSolver();
const int numSolvers = linSolver.numAvailableSolvers();
if (numSolvers > 1 && (linSolver.getSolveCount() % 100 == 0)) {
if (terminal_output_) {
OpmLog::debug("\nRunning speed test for comparing available linear solvers.");
}
Dune::Timer perfTimer;
std::vector<double> times(numSolvers);
std::vector<double> setupTimes(numSolvers);
x = 0.0;
std::vector<BVector> x_trial(numSolvers, x);
for (int solver = 0; solver < numSolvers; ++solver) {
BVector x0(x);
linSolver.setActiveSolver(solver);
perfTimer.start();
linSolver.prepare(jacobian, residual);
setupTimes[solver] = perfTimer.stop();
perfTimer.reset();
linSolver.setResidual(residual);
perfTimer.start();
linSolver.solve(x_trial[solver]);
times[solver] = perfTimer.stop();
perfTimer.reset();
if (terminal_output_) {
OpmLog::debug(fmt::format(fmt::runtime("Solver time {}: {}"), solver, times[solver]));
}
}
int fastest_solver = std::ranges::min_element(times) - times.begin();
// Use timing on rank 0 to determine fastest, must be consistent across ranks.
grid_.comm().broadcast(&fastest_solver, 1, 0);
linear_solve_setup_time_ = setupTimes[fastest_solver];
x = x_trial[fastest_solver];
linSolver.setActiveSolver(fastest_solver);
}
else {
// set initial guess
x = 0.0;
Dune::Timer perfTimer;
perfTimer.start();
linSolver.prepare(jacobian, residual);
linear_solve_setup_time_ = perfTimer.stop();
linSolver.setResidual(residual);
// actually, the error needs to be calculated after setResidual in order to
// account for parallelization properly. since the residual of ECFV
// discretizations does not need to be synchronized across processes to be
// consistent, this is not relevant for OPM-flow...
linSolver.solve(x);
}
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
updateSolution(const BVector& dx)
{
OPM_TIMEBLOCK(updateSolution);
// Prepare to store solution update for convergence check
if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0
|| this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) {
prepareStoringSolutionUpdate();
}
auto& newtonMethod = simulator_.model().newtonMethod();
SolutionVector& solution = simulator_.model().solution(/*timeIdx=*/0);
newtonMethod.update_(/*nextSolution=*/solution,
/*curSolution=*/solution,
/*update=*/dx,
/*resid=*/dx); // the update routines of the black
// oil model do not care about the
// residual
// if the solution is updated, the intensive quantities need to be recalculated
{
OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
}
// Store solution update
if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0
|| this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) {
storeSolutionUpdate(dx);
}
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
prepareStoringSolutionUpdate()
{
// Init. solution update vector
unsigned nc = simulator_.model().numGridDof();
solUpd_.resize(nc);
const auto& elemMapper = simulator_.model().elementMapper();
const auto& gridView = simulator_.gridView();
for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
// Copy solution vector to transfer primary variables meaning
unsigned globalElemIdx = elemMapper.index(elem);
solUpd_[globalElemIdx] = simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
// Ensure each element is zero
std::ranges::fill(solUpd_[globalElemIdx], 0.0);
}
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
storeSolutionUpdate(const BVector& dx)
{
const auto& elemMapper = simulator_.model().elementMapper();
const auto& gridView = simulator_.gridView();
for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
// Get cell vectors
unsigned globalElemIdx = elemMapper.index(elem);
auto& value = solUpd_[globalElemIdx];
const auto& update = dx[globalElemIdx];
assert(value.size() == update.size());
// Transfer update from dx to solution update container (SolutionVector type)
std::ranges::copy(update, value.begin());
}
}
template <class TypeTag>
typename NonlinearSystemBlackOilReservoir<TypeTag>::MaxSolutionUpdateData
NonlinearSystemBlackOilReservoir<TypeTag>::
getMaxSolutionUpdate(const std::vector<unsigned>& ixCells)
{
static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
// Init output
Scalar dPMax = 0.0;
Scalar dSMax = 0.0;
Scalar dRsMax = 0.0;
Scalar dRvMax = 0.0;
// Loop over solution update, get the correct variables and calculate max.
for (const auto& ix : ixCells) {
const auto& value = solUpd_[ix];
for (int pvIdx = 0; pvIdx < static_cast<int>(value.size()); ++pvIdx) {
if (pvIdx == Indices::pressureSwitchIdx) {
dPMax = std::max(dPMax, std::abs(value[pvIdx]));
}
else if ( (pvIdx == Indices::waterSwitchIdx
&& value.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
|| (pvIdx == Indices::compositionSwitchIdx
&& value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
|| (enableSolvent && pvIdx == Indices::solventSaturationIdx
&& value.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss)
|| (enableBrine && enableSaltPrecipitation && pvIdx == Indices::saltConcentrationIdx
&& value.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) ) {
dSMax = std::max(dSMax, std::abs(value[pvIdx]));
}
else if (pvIdx == Indices::compositionSwitchIdx
&& value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
dRsMax = std::max(dRsMax, std::abs(value[pvIdx]));
}
else if (pvIdx == Indices::compositionSwitchIdx
&& value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
dRvMax = std::max(dRvMax, std::abs(value[pvIdx]));
}
}
}
// Communicate max values
dPMax = this->grid_.comm().max(dPMax);
dSMax = this->grid_.comm().max(dSMax);
dRsMax = this->grid_.comm().max(dRsMax);
dRvMax = this->grid_.comm().max(dRvMax);
return { dPMax, dSMax, dRsMax, dRvMax };
}
template <class TypeTag>
std::tuple<typename NonlinearSystemBlackOilReservoir<TypeTag>::Scalar,
typename NonlinearSystemBlackOilReservoir<TypeTag>::Scalar>
NonlinearSystemBlackOilReservoir<TypeTag>::
convergenceReduction(Parallel::Communication comm,
const Scalar pvSumLocal,
const Scalar numAquiferPvSumLocal,
std::vector< Scalar >& R_sum,
std::vector< Scalar >& maxCoeff,
std::vector< Scalar >& B_avg)
{
OPM_TIMEBLOCK(convergenceReduction);
// Compute total pore volume (use only owned entries)
Scalar pvSum = pvSumLocal;
Scalar numAquiferPvSum = numAquiferPvSumLocal;
if (comm.size() > 1) {
// global reduction
std::vector< Scalar > sumBuffer;
std::vector< Scalar > maxBuffer;
const int numComp = B_avg.size();
sumBuffer.reserve( 2*numComp + 2 ); // +2 for (numAquifer)pvSum
maxBuffer.reserve( numComp );
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
sumBuffer.push_back(B_avg[compIdx]);
sumBuffer.push_back(R_sum[compIdx]);
maxBuffer.push_back(maxCoeff[ compIdx]);
}
// Compute total pore volume
sumBuffer.push_back( pvSum );
sumBuffer.push_back( numAquiferPvSum );
// compute global sum
comm.sum( sumBuffer.data(), sumBuffer.size() );
// compute global max
comm.max( maxBuffer.data(), maxBuffer.size() );
// restore values to local variables
for (int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx) {
B_avg[compIdx] = sumBuffer[buffIdx];
++buffIdx;
R_sum[compIdx] = sumBuffer[buffIdx];
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
maxCoeff[compIdx] = maxBuffer[compIdx];
}
// restore global pore volume
pvSum = sumBuffer[sumBuffer.size()-2];
numAquiferPvSum = sumBuffer.back();
}
// return global pore volume
return {pvSum, numAquiferPvSum};
}
template <class TypeTag>
std::pair<typename NonlinearSystemBlackOilReservoir<TypeTag>::Scalar,
typename NonlinearSystemBlackOilReservoir<TypeTag>::Scalar>
NonlinearSystemBlackOilReservoir<TypeTag>::
localConvergenceData(std::vector<Scalar>& R_sum,
std::vector<Scalar>& maxCoeff,
std::vector<Scalar>& B_avg,
std::vector<int>& maxCoeffCell)
{
OPM_TIMEBLOCK(localConvergenceData);
Scalar pvSumLocal = 0.0;
Scalar numAquiferPvSumLocal = 0.0;
const auto& model = simulator_.model();
const auto& problem = simulator_.problem();
const auto& residual = simulator_.model().linearizer().residual();
ElementContext elemCtx(simulator_);
const auto& gridView = simulator().gridView();
IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) *
model.dofTotalVolume(cell_idx);
pvSumLocal += pvValue;
if (isNumericalAquiferCell(elem)) {
numAquiferPvSumLocal += pvValue;
}
this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue,
B_avg, R_sum, maxCoeff, maxCoeffCell);
}
OPM_END_PARALLEL_TRY_CATCH("NonlinearSystemBlackOilReservoir::localConvergenceData() failed: ", grid_.comm());
// compute local average in terms of global number of elements
const int bSize = B_avg.size();
for (int i = 0; i < bSize; ++i) {
B_avg[i] /= Scalar(global_nc_);
}
return {pvSumLocal, numAquiferPvSumLocal};
}
template <class TypeTag>
typename NonlinearSystemBlackOilReservoir<TypeTag>::CnvPvSplitData
NonlinearSystemBlackOilReservoir<TypeTag>::
characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt)
{
OPM_TIMEBLOCK(computeCnvErrorPv);
// 0: cnv <= tolerance_cnv
// 1: tolerance_cnv < cnv <= tolerance_cnv_relaxed
// 2: tolerance_cnv_relaxed < cnv
constexpr auto numPvGroups = std::vector<double>::size_type{3};
auto cnvPvSplit = std::pair<std::vector<double>, std::vector<int>> {
std::piecewise_construct,
std::forward_as_tuple(numPvGroups),
std::forward_as_tuple(numPvGroups)
};
auto maxCNV = [&B_avg, dt](const auto& residual, const double pvol)
{
return (dt / pvol) *
std::inner_product(residual.begin(), residual.end(),
B_avg.begin(), Scalar{0},
[](const Scalar m, const auto& x)
{
using std::abs;
return std::max(m, abs(x));
}, std::multiplies<>{});
};
auto& [splitPV, cellCntPV] = cnvPvSplit;
const auto& model = this->simulator().model();
const auto& problem = this->simulator().problem();
const auto& residual = model.linearizer().residual();
const auto& gridView = this->simulator().gridView();
const IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
ElementContext elemCtx(this->simulator());
std::vector<unsigned> ixCells;
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
// Skip cells of numerical Aquifer
if (isNumericalAquiferCell(elem)) {
continue;
}
elemCtx.updatePrimaryStencil(elem);
const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0)
* model.dofTotalVolume(cell_idx);
const auto maxCnv = maxCNV(residual[cell_idx], pvValue);
const auto ix = (maxCnv > this->param_.tolerance_cnv_)
+ (maxCnv > this->param_.tolerance_cnv_relaxed_);
splitPV[ix] += static_cast<double>(pvValue);
++cellCntPV[ix];
// For dP and dS check, we need cell indices of [1] violations
if ( ix > 0 &&
(this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0
|| this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0 ) ) {
ixCells.push_back(cell_idx);
}
}
OPM_END_PARALLEL_TRY_CATCH("NonlinearSystemBlackOilReservoir::characteriseCnvPvSplit() failed: ",
this->grid_.comm());
this->grid_.comm().sum(splitPV .data(), splitPV .size());
this->grid_.comm().sum(cellCntPV.data(), cellCntPV.size());
return { cnvPvSplit, ixCells };
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
updateTUNING(const Tuning& tuning)
{
this->param_.tolerance_cnv_ = tuning.TRGCNV;
this->param_.tolerance_cnv_relaxed_ = tuning.XXXCNV;
this->param_.tolerance_mb_ = tuning.TRGMBE;
this->param_.tolerance_mb_relaxed_ = tuning.XXXMBE;
this->param_.newton_max_iter_ = tuning.NEWTMX;
this->param_.newton_min_iter_ = tuning.NEWTMN;
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
updateTUNINGDP(const TuningDp& tuning_dp)
{
// NOTE: If TUNINGDP item is _not_ set it should be 0.0
this->param_.tolerance_max_dp_ = tuning_dp.TRGDDP;
this->param_.tolerance_max_ds_ = tuning_dp.TRGDDS;
this->param_.tolerance_max_drs_ = tuning_dp.TRGDDRS;
this->param_.tolerance_max_drv_ = tuning_dp.TRGDDRV;
}
template <class TypeTag>
ConvergenceReport
NonlinearSystemBlackOilReservoir<TypeTag>::
getReservoirConvergence(const double reportTime,
const double dt,
const int maxIter,
std::vector<Scalar>& B_avg,
std::vector<Scalar>& residual_norms)
{
OPM_TIMEBLOCK(getReservoirConvergence);
using Vector = std::vector<Scalar>;
const auto& iterCtx = simulator_.problem().iterationContext();
ConvergenceReport report{reportTime};
const int numComp = numEq;
Vector R_sum(numComp, Scalar{0});
Vector maxCoeff(numComp, std::numeric_limits<Scalar>::lowest());
std::vector<int> maxCoeffCell(numComp, -1);
const auto [pvSumLocal, numAquiferPvSumLocal] =
this->localConvergenceData(R_sum, maxCoeff, B_avg, maxCoeffCell);
// compute global sum and max of quantities
const auto& [pvSum, numAquiferPvSum] =
this->convergenceReduction(this->grid_.comm(),
pvSumLocal,
numAquiferPvSumLocal,
R_sum, maxCoeff, B_avg);
auto cnvSplitData = this->characteriseCnvPvSplit(B_avg, dt);
report.setCnvPoreVolSplit(cnvSplitData.cnvPvSplit,
pvSum - numAquiferPvSum);
// For each iteration, we need to determine whether to use the
// relaxed tolerances. To disable the usage of relaxed
// tolerances, you can set the relaxed tolerances as the strict
// tolerances. If min_strict_mb_iter = -1 (default) we use a
// relaxed tolerance for the mass balance for the last
// iterations. For positive values we use the relaxed tolerance
// after the given number of iterations
const bool relax_final_iteration_mb =
this->param_.min_strict_mb_iter_ < 0 && iterCtx.iteration() == maxIter;
const bool relax_iter_mb = this->param_.min_strict_mb_iter_ >= 0 &&
iterCtx.shouldRelax(this->param_.min_strict_mb_iter_);
const bool use_relaxed_mb = relax_final_iteration_mb
|| relax_iter_mb;
// If min_strict_cnv_iter = -1 we use a relaxed tolerance for
// the cnv for the last iterations. For positive values we use
// the relaxed tolerance after the given number of iterations.
// We also use relaxed tolerances for cells with total
// pore-volume less than relaxed_max_pv_fraction_. Default
// value of relaxed_max_pv_fraction_ is 0.03
const bool relax_final_iteration_cnv =
this->param_.min_strict_cnv_iter_ < 0 && iterCtx.iteration() == maxIter;
const bool relax_iter_cnv = this->param_.min_strict_cnv_iter_ >= 0 &&
iterCtx.shouldRelax(this->param_.min_strict_cnv_iter_);
// Note trailing parentheses here, just before the final
// semicolon. This is an immediately invoked function
// expression which calculates a single boolean value.
const auto relax_pv_fraction_cnv =
[&report, this, eligible = pvSum - numAquiferPvSum]()
{
const auto& cnvPvSplit = report.cnvPvSplit().first;
// [1]: tol < cnv <= relaxed
// [2]: relaxed < cnv
Scalar cnvPvSum = static_cast<Scalar>(cnvPvSplit[1] + cnvPvSplit[2]);
return cnvPvSum < this->param_.relaxed_max_pv_fraction_ * eligible &&
cnvPvSum > 0.0;
}();
// If tolerances for solution changes are met, we use the
// relaxed cnv tolerance. Note that all tolerances > 0.0
// must be met to enable use of relaxed cnv tolerance.
MaxSolutionUpdateData maxSolUpd;
const bool use_dp_tol = this->param_.tolerance_max_dp_ > 0.0;
const bool use_ds_tol = this->param_.tolerance_max_ds_ > 0.0;
const bool use_drs_tol = this->param_.tolerance_max_drs_ > 0.0;
const bool use_drv_tol = this->param_.tolerance_max_drv_ > 0.0;
const bool use_dsol_tol = use_dp_tol || use_ds_tol || use_drs_tol || use_drv_tol;
bool relax_dsol_cnv = false;
if (!iterCtx.isFirstGlobalIteration() && use_dsol_tol) {
maxSolUpd = getMaxSolutionUpdate(cnvSplitData.ixCells);
relax_dsol_cnv =
(!use_dp_tol || (maxSolUpd.dPMax > 0.0 && maxSolUpd.dPMax < this->param_.tolerance_max_dp_)) &&
(!use_ds_tol || (maxSolUpd.dSMax > 0.0 && maxSolUpd.dSMax < this->param_.tolerance_max_ds_)) &&
(!use_drs_tol || (maxSolUpd.dRsMax > 0.0 && maxSolUpd.dRsMax < this->param_.tolerance_max_drs_)) &&
(!use_drv_tol || (maxSolUpd.dRvMax > 0.0 && maxSolUpd.dRvMax < this->param_.tolerance_max_drv_));
}
// Determine if relaxed CNV tolerances should be used
const bool use_relaxed_cnv = relax_final_iteration_cnv
|| relax_pv_fraction_cnv
|| relax_iter_cnv
|| relax_dsol_cnv;
// Ensure that CNV convergence criteria is met when max.
// solution change tolerances have been fulfilled
Scalar tolerance_cnv_relaxed = relax_dsol_cnv ? 1e20 : param_.tolerance_cnv_relaxed_;
const auto tol_cnv = use_relaxed_cnv ? tolerance_cnv_relaxed : param_.tolerance_cnv_;
const auto tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
const auto tol_cnv_energy = use_relaxed_cnv ? param_.tolerance_cnv_energy_relaxed_ : param_.tolerance_cnv_energy_;
const auto tol_eb = use_relaxed_mb ? param_.tolerance_energy_balance_relaxed_ : param_.tolerance_energy_balance_;
// Finish computation
std::vector<Scalar> CNV(numComp);
std::vector<Scalar> mass_balance_residual(numComp);
for (int compIdx = 0; compIdx < numComp; ++compIdx)
{
CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
residual_norms.push_back(CNV[compIdx]);
}
using CR = ConvergenceReport;
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
const Scalar res[2] = {
mass_balance_residual[compIdx], CNV[compIdx],
};
const CR::ReservoirFailure::Type types[2] = {
CR::ReservoirFailure::Type::MassBalance,
CR::ReservoirFailure::Type::Cnv,
};
Scalar tol[2] = { tol_mb, tol_cnv, };
if (has_energy_ && compIdx == contiEnergyEqIdx) {
tol[0] = tol_eb;
tol[1] = tol_cnv_energy;
}
for (int ii : {0, 1}) {
if (std::isnan(res[ii])) {
report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
if (this->terminal_output_) {
OpmLog::debug("NaN residual for " + this->compNames_.name(compIdx) + " equation.");
}
}
else if (res[ii] > maxResidualAllowed()) {
report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
if (this->terminal_output_) {
OpmLog::debug("Too large residual for " + this->compNames_.name(compIdx) + " equation.");
}
}
else if (res[ii] < 0.0) {
report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
if (this->terminal_output_) {
OpmLog::debug("Negative residual for " + this->compNames_.name(compIdx) + " equation.");
}
}
else if (res[ii] > tol[ii]) {
report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
}
report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
}
}
// Compute the Newton convergence per cell.
this->convergencePerCell(B_avg, dt, tol_cnv, tol_cnv_energy);
// Output of residuals.
if (this->terminal_output_) {
// Only rank 0 does print to std::cout
if (iterCtx.isFirstGlobalIteration()) {
std::string msg = "Iter";
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " MB(";
msg += this->compNames_.name(compIdx)[0];
msg += ") ";
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " CNV(";
msg += this->compNames_.name(compIdx)[0];
msg += ") ";
}
if (use_dsol_tol) {
msg += use_dp_tol ? " DP " : "";
msg += use_ds_tol ? " DS " : "";
msg += use_drs_tol ? " DRS " : "";
msg += use_drv_tol ? " DRV " : "";
}
msg += " MBFLAG";
msg += " CNVFLAG";
OpmLog::debug(msg);
}
std::ostringstream ss;
const std::streamsize oprec = ss.precision(3);
const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
ss << std::setw(4) << iterCtx.iteration();
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
ss << std::setw(11) << mass_balance_residual[compIdx];
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
ss << std::setw(11) << CNV[compIdx];
}
if (use_dsol_tol) {
auto print_dsol =
[&] (bool use_tol, Scalar dsol) {
if (!use_tol) {
return;
}
if (iterCtx.isFirstGlobalIteration() || dsol <= 0.0) {
ss << std::string(5, ' ') << "-" << std::string(5, ' ');
}
else {
ss << std::setw(11) << dsol;
}
};
print_dsol(use_dp_tol, maxSolUpd.dPMax);
print_dsol(use_ds_tol, maxSolUpd.dSMax);
print_dsol(use_drs_tol, maxSolUpd.dRsMax);
print_dsol(use_drv_tol, maxSolUpd.dRvMax);
}
const auto mb_flag = use_relaxed_mb
? DebugFlags::RELAXED
: DebugFlags::STRICT;
const auto cnv_flag = relax_dsol_cnv ?
DebugFlags::TUNINGDP
: (use_relaxed_cnv
? DebugFlags::RELAXED
: DebugFlags::STRICT);
ss << std::setw(9) << make_string<TypeTag>(mb_flag)
<< std::setw(9) << make_string<TypeTag>(cnv_flag);
ss.precision(oprec);
ss.flags(oflags);
OpmLog::debug(ss.str());
}
return report;
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
convergencePerCell(const std::vector<Scalar>& B_avg,
const double dt,
const double tol_cnv,
const double tol_cnv_energy)
{
auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
if (!rst_conv.hasConv()) {
return;
}
if (simulator_.problem().iterationContext().isFirstGlobalIteration()) {
rst_conv.prepareConv();
}
const auto& residual = simulator_.model().linearizer().residual();
const auto& gridView = this->simulator().gridView();
const IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
ElementContext elemCtx(this->simulator());
std::vector<int> convNewt(residual.size(), 0);
OPM_BEGIN_PARALLEL_TRY_CATCH();
unsigned idx = 0;
const int numComp = B_avg.size();
for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
elemCtx.updatePrimaryStencil(elem);
const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto pvValue = simulator_.problem().referencePorosity(cell_idx, /*timeIdx=*/0) *
simulator_.model().dofTotalVolume(cell_idx);
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
const auto tol = (has_energy_ && compIdx == contiEnergyEqIdx) ? tol_cnv_energy : tol_cnv;
const Scalar cnv = std::abs(B_avg[compIdx] * residual[cell_idx][compIdx]) * dt / pvValue;
if (std::isnan(cnv) || cnv > maxResidualAllowed() || cnv < 0.0 || cnv > tol) {
convNewt[idx] = 1;
break;
}
}
++idx;
}
OPM_END_PARALLEL_TRY_CATCH("NonlinearSystemBlackOilReservoir::convergencePerCell() failed: ",
this->grid_.comm());
rst_conv.updateNewton(convNewt);
}
template <class TypeTag>
ConvergenceReport
NonlinearSystemBlackOilReservoir<TypeTag>::
getConvergence(const SimulatorTimerInterface& timer,
const int maxIter,
std::vector<Scalar>& residual_norms)
{
OPM_TIMEBLOCK(getConvergence);
// Get convergence reports for reservoir and wells.
std::vector<Scalar> B_avg(numEq, 0.0);
auto report = getReservoirConvergence(timer.simulationTimeElapsed(),
timer.currentStepLength(),
maxIter, B_avg, residual_norms);
{
OPM_TIMEBLOCK(getWellConvergence);
report += wellModel().getWellConvergence(B_avg,
/*checkWellGroupControlsAndNetwork*/report.converged());
}
conv_monitor_.checkPenaltyCard(report, simulator_.problem().iterationContext().iteration());
return report;
}
template <class TypeTag>
std::vector<std::vector<typename NonlinearSystemBlackOilReservoir<TypeTag>::Scalar> >
NonlinearSystemBlackOilReservoir<TypeTag>::
computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
{
OPM_TIMEBLOCK(computeFluidInPlace);
// assert(true)
// return an empty vector
std::vector<std::vector<Scalar> > regionValues(0, std::vector<Scalar>(0,0.0));
return regionValues;
}
template <class TypeTag>
const SimulatorReport&
NonlinearSystemBlackOilReservoir<TypeTag>::
localAccumulatedReports() const
{
if (!hasNlddSolver()) {
OPM_THROW(std::runtime_error, "Cannot get local reports from a model without NLDD solver");
}
return nlddSolver_->localAccumulatedReports();
}
template <class TypeTag>
const std::vector<SimulatorReport>&
NonlinearSystemBlackOilReservoir<TypeTag>::
domainAccumulatedReports() const
{
if (!nlddSolver_)
OPM_THROW(std::runtime_error, "Cannot get domain reports from a model without NLDD solver");
return nlddSolver_->domainAccumulatedReports();
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
writeNonlinearIterationsPerCell(const std::filesystem::path& odir) const
{
if (hasNlddSolver()) {
nlddSolver_->writeNonlinearIterationsPerCell(odir);
}
}
template <class TypeTag>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
writePartitions(const std::filesystem::path& odir) const
{
if (hasNlddSolver()) {
nlddSolver_->writePartitions(odir);
return;
}
const auto& elementMapper = this->simulator().model().elementMapper();
const auto& cartMapper = this->simulator().vanguard().cartesianIndexMapper();
const auto& grid = this->simulator().vanguard().grid();
const auto& comm = grid.comm();
const auto nDigit = 1 + static_cast<int>(std::floor(std::log10(comm.size())));
std::ofstream pfile {odir / fmt::format("{1:0>{0}}", nDigit, comm.rank())};
for (const auto& cell : elements(grid.leafGridView(), Dune::Partitions::interior)) {
pfile << comm.rank() << ' '
<< cartMapper.cartesianIndex(elementMapper.index(cell)) << ' '
<< comm.rank() << '\n';
}
}
template <class TypeTag>
template<class FluidState, class Residual>
void
NonlinearSystemBlackOilReservoir<TypeTag>::
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)
{
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
{
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
const unsigned compIdx = FluidSystem::canonicalToActiveCompIdx(sIdx);
B_avg[compIdx] += 1.0 / fs.invB(phaseIdx).value();
const auto R2 = modelResid[cell_idx][compIdx];
R_sum[compIdx] += R2;
const Scalar Rval = std::abs(R2) / pvValue;
if (Rval > maxCoeff[compIdx]) {
maxCoeff[compIdx] = Rval;
maxCoeffCell[compIdx] = cell_idx;
}
}
if constexpr (has_solvent_) {
B_avg[contiSolventEqIdx] +=
1.0 / intQuants.solventInverseFormationVolumeFactor().value();
const auto R2 = modelResid[cell_idx][contiSolventEqIdx];
R_sum[contiSolventEqIdx] += R2;
maxCoeff[contiSolventEqIdx] = std::max(maxCoeff[contiSolventEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_extbo_) {
B_avg[contiZfracEqIdx] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
const auto R2 = modelResid[cell_idx][contiZfracEqIdx];
R_sum[ contiZfracEqIdx ] += R2;
maxCoeff[contiZfracEqIdx] = std::max(maxCoeff[contiZfracEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_polymer_) {
B_avg[contiPolymerEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R2 = modelResid[cell_idx][contiPolymerEqIdx];
R_sum[contiPolymerEqIdx] += R2;
maxCoeff[contiPolymerEqIdx] = std::max(maxCoeff[contiPolymerEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_foam_) {
B_avg[ contiFoamEqIdx ] += 1.0 / fs.invB(FluidSystem::gasPhaseIdx).value();
const auto R2 = modelResid[cell_idx][contiFoamEqIdx];
R_sum[contiFoamEqIdx] += R2;
maxCoeff[contiFoamEqIdx] = std::max(maxCoeff[contiFoamEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_brine_) {
B_avg[ contiBrineEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R2 = modelResid[cell_idx][contiBrineEqIdx];
R_sum[contiBrineEqIdx] += R2;
maxCoeff[contiBrineEqIdx] = std::max(maxCoeff[contiBrineEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_polymermw_) {
static_assert(has_polymer_);
B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
// the residual of the polymer molecular equation is scaled down by a 100, since molecular weight
// can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations
// TODO: there should be a more general way to determine the scaling-down coefficient
const auto R2 = modelResid[cell_idx][contiPolymerMWEqIdx] / 100.;
R_sum[contiPolymerMWEqIdx] += R2;
maxCoeff[contiPolymerMWEqIdx] = std::max(maxCoeff[contiPolymerMWEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_energy_) {
B_avg[contiEnergyEqIdx] += 1.0;
const auto R2 = modelResid[cell_idx][contiEnergyEqIdx];
R_sum[contiEnergyEqIdx] += R2;
maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_bioeffects_) {
B_avg[contiMicrobialEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R1 = modelResid[cell_idx][contiMicrobialEqIdx];
R_sum[contiMicrobialEqIdx] += R1;
maxCoeff[contiMicrobialEqIdx] = std::max(maxCoeff[contiMicrobialEqIdx],
std::abs(R1) / pvValue);
B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R2 = modelResid[cell_idx][contiBiofilmEqIdx];
R_sum[contiBiofilmEqIdx] += R2;
maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
std::abs(R2) / pvValue);
if constexpr (has_micp_) {
B_avg[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R3 = modelResid[cell_idx][contiOxygenEqIdx];
R_sum[contiOxygenEqIdx] += R3;
maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
std::abs(R3) / pvValue);
B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R4 = modelResid[cell_idx][contiUreaEqIdx];
R_sum[contiUreaEqIdx] += R4;
maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
std::abs(R4) / pvValue);
B_avg[contiCalciteEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R5 = modelResid[cell_idx][contiCalciteEqIdx];
R_sum[contiCalciteEqIdx] += R5;
maxCoeff[contiCalciteEqIdx] = std::max(maxCoeff[contiCalciteEqIdx],
std::abs(R5) / pvValue);
}
}
}
} // namespace Opm
#endif // OPM_NONLINEAR_SYSTEM_BLACK_OIL_RESERVOIR_IMPL_HEADER_INCLUDED
Reference in new issue View Git Blame Copy Permalink