OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity
918b6cc594
opm-simulators/opm/simulators/flow/BlackoilModel_impl.hpp
Arne Morten Kvarving dcd739fa09 BlackoilModel: add impl file
2025-01-15 11:51:39 +01:00

1139 lines
43 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_BLACKOILMODEL_IMPL_HEADER_INCLUDED
#define OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED
#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED
#include <config.h>
#include <opm/simulators/flow/BlackoilModel.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 <opm/simulators/utils/phaseUsageFromDeck.hpp>
#include <algorithm>
#include <iomanip>
#include <limits>
#include <stdexcept>
#include <sstream>
#include <fmt/format.h>
namespace Opm {
template <class TypeTag>
BlackoilModel<TypeTag>::
BlackoilModel(Simulator& simulator,
const ModelParameters& param,
BlackoilWellModel<TypeTag>& well_model,
const bool terminal_output)
: simulator_(simulator)
, grid_(simulator_.vanguard().grid())
, phaseUsage_(phaseUsageFromDeck(eclState()))
, param_( param )
, well_model_ (well_model)
, rst_conv_(simulator_.problem().eclWriter()->collectOnIORank().localIdxToGlobalIdxMapping(),
grid_.comm())
, 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<BlackoilModelNldd<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
BlackoilModel<TypeTag>::
prepareStep(const SimulatorTimerInterface& timer)
{
SimulatorReportSingle report;
Dune::Timer perfTimer;
perfTimer.start();
// update the solution variables in the model
int lastStepFailed = timer.lastStepFailed();
if (grid_.comm().size() > 1 && lastStepFailed != grid_.comm().min(lastStepFailed)) {
OPM_THROW(std::runtime_error,
fmt::format("Misalignment of the parallel simulation run in prepareStep "
"- the previous step succeeded on rank {} but failed on the "
"other ranks.", grid_.comm().rank()));
}
if (lastStepFailed) {
simulator_.model().updateFailed();
}
else {
simulator_.model().advanceTimeLevel();
}
// Set the timestep size, episode index, and non-linear iteration 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_.model().newtonMethod().setIterationIndex(0);
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 (nlddSolver_) {
nlddSolver_->prepareStep();
}
report.pre_post_time += perfTimer.stop();
auto getIdx = [](unsigned phaseIdx) -> int
{
if (FluidSystem::phaseIsActive(phaseIdx)) {
const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
return Indices::canonicalToActiveComponentIndex(sIdx);
}
return -1;
};
const auto& schedule = simulator_.vanguard().schedule();
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
BlackoilModel<TypeTag>::
initialLinearization(SimulatorReportSingle& report,
const int iteration,
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, iteration);
report.assemble_time += perfTimer.stop();
}
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, iteration, maxIter, residual_norms);
report.converged = convrep.converged() && 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
BlackoilModel<TypeTag>::
nonlinearIteration(const int iteration,
const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver)
{
if (iteration == 0) {
// For each iteration we store in a vector the norms of the residual of
// the mass balance for each active phase, the well flux and the well equations.
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") ||
(iteration < this->param_.nldd_num_initial_newton_iter_))
{
result = this->nonlinearIterationNewton(iteration,
timer,
nonlinear_solver);
}
else {
result = this->nlddSolver_->nonlinearIterationNldd(iteration,
timer,
nonlinear_solver);
}
rst_conv_.update(simulator_.model().linearizer().residual());
return result;
}
template <class TypeTag>
template <class NonlinearSolverType>
SimulatorReportSingle
BlackoilModel<TypeTag>::
nonlinearIterationNewton(const int iteration,
const SimulatorTimerInterface& timer,
NonlinearSolverType& nonlinear_solver)
{
// ----------- Set up reports and timer -----------
SimulatorReportSingle report;
Dune::Timer perfTimer;
this->initialLinearization(report,
iteration,
nonlinear_solver.minIter(),
nonlinear_solver.maxIter(),
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
BlackoilModel<TypeTag>::
afterStep(const SimulatorTimerInterface&)
{
SimulatorReportSingle report;
Dune::Timer perfTimer;
perfTimer.start();
simulator_.problem().endTimeStep();
simulator_.problem().setConvData(rst_conv_.getData());
report.pre_post_time += perfTimer.stop();
return report;
}
template <class TypeTag>
SimulatorReportSingle
BlackoilModel<TypeTag>::
assembleReservoir(const SimulatorTimerInterface& /* timer */,
const int iterationIdx)
{
// -------- Mass balance equations --------
simulator_.model().newtonMethod().setIterationIndex(iterationIdx);
simulator_.problem().beginIteration();
simulator_.model().linearizer().linearizeDomain();
simulator_.problem().endIteration();
return wellModel().lastReport();
}
template <class TypeTag>
typename BlackoilModel<TypeTag>::Scalar
BlackoilModel<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
BlackoilModel<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("Solver time {}: {}", solver, times[solver]));
}
}
int fastest_solver = std::min_element(times.begin(), times.end()) - 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
BlackoilModel<TypeTag>::
updateSolution(const BVector& dx)
{
OPM_TIMEBLOCK(updateSolution);
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);
}
}
template <class TypeTag>
std::tuple<typename BlackoilModel<TypeTag>::Scalar,
typename BlackoilModel<TypeTag>::Scalar>
BlackoilModel<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 BlackoilModel<TypeTag>::Scalar,
typename BlackoilModel<TypeTag>::Scalar>
BlackoilModel<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("BlackoilModel::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>
std::pair<std::vector<double>, std::vector<int>>
BlackoilModel<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());
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];
}
OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::characteriseCnvPvSplit() failed: ",
this->grid_.comm());
this->grid_.comm().sum(splitPV .data(), splitPV .size());
this->grid_.comm().sum(cellCntPV.data(), cellCntPV.size());
return cnvPvSplit;
}
template <class TypeTag>
void
BlackoilModel<TypeTag>::
updateTUNING(const Tuning& tuning)
{
this->param_.tolerance_mb_ = tuning.XXXMBE;
if (terminal_output_) {
OpmLog::debug(fmt::format("Setting BlackoilModel mass "
"balance limit (XXXMBE) to {:.2e}",
tuning.XXXMBE));
}
}
template <class TypeTag>
ConvergenceReport
BlackoilModel<TypeTag>::
getReservoirConvergence(const double reportTime,
const double dt,
const int iteration,
const int maxIter,
std::vector<Scalar>& B_avg,
std::vector<Scalar>& residual_norms)
{
OPM_TIMEBLOCK(getReservoirConvergence);
using Vector = std::vector<Scalar>;
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);
report.setCnvPoreVolSplit(this->characteriseCnvPvSplit(B_avg, dt),
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 && iteration == maxIter;
const bool use_relaxed_mb = relax_final_iteration_mb ||
(this->param_.min_strict_mb_iter_ >= 0 &&
iteration >= this->param_.min_strict_mb_iter_);
// 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 && iteration == maxIter;
const bool relax_iter_cnv =
this->param_.min_strict_cnv_iter_ >= 0 &&
iteration >= 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
return static_cast<Scalar>(cnvPvSplit[1] + cnvPvSplit[2]) <
this->param_.relaxed_max_pv_fraction_ * eligible;
}();
const bool use_relaxed_cnv = relax_final_iteration_cnv
|| relax_pv_fraction_cnv
|| relax_iter_cnv;
if ((relax_final_iteration_mb || relax_final_iteration_cnv) &&
this->terminal_output_)
{
std::string message =
"Number of newton iterations reached its maximum "
"try to continue with relaxed tolerances:";
if (relax_final_iteration_mb) {
message += fmt::format(" MB: {:.1e}", param_.tolerance_mb_relaxed_);
}
if (relax_final_iteration_cnv) {
message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_);
}
OpmLog::debug(message);
}
const auto tol_cnv = use_relaxed_cnv ? param_.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]);
}
}
// Output of residuals.
if (this->terminal_output_) {
// Only rank 0 does print to std::cout
if (iteration == 0) {
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 += ") ";
}
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) << 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];
}
ss.precision(oprec);
ss.flags(oflags);
OpmLog::debug(ss.str());
}
return report;
}
template <class TypeTag>
ConvergenceReport
BlackoilModel<TypeTag>::
getConvergence(const SimulatorTimerInterface& timer,
const int iteration,
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(),
iteration, maxIter, B_avg, residual_norms);
{
OPM_TIMEBLOCK(getWellConvergence);
report += wellModel().getWellConvergence(B_avg,
/*checkWellGroupControls*/report.converged());
}
conv_monitor_.checkPenaltyCard(report, iteration);
return report;
}
template <class TypeTag>
std::vector<std::vector<typename BlackoilModel<TypeTag>::Scalar> >
BlackoilModel<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>
SimulatorReportSingle
BlackoilModel<TypeTag>::
localAccumulatedReports() const
{
return nlddSolver_ ? nlddSolver_->localAccumulatedReports()
: SimulatorReportSingle{};
}
template <class TypeTag>
void
BlackoilModel<TypeTag>::
writePartitions(const std::filesystem::path& odir) const
{
if (this->nlddSolver_ != nullptr) {
this->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
BlackoilModel<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 = Indices::canonicalToActiveComponentIndex(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 / (4.182e1); // converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water
const auto R2 = modelResid[cell_idx][contiEnergyEqIdx];
R_sum[contiEnergyEqIdx] += R2;
maxCoeff[contiEnergyEqIdx] = std::max(maxCoeff[contiEnergyEqIdx],
std::abs(R2) / pvValue);
}
if constexpr (has_micp_) {
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[contiOxygenEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R2 = modelResid[cell_idx][contiOxygenEqIdx];
R_sum[contiOxygenEqIdx] += R2;
maxCoeff[contiOxygenEqIdx] = std::max(maxCoeff[contiOxygenEqIdx],
std::abs(R2) / pvValue);
B_avg[contiUreaEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R3 = modelResid[cell_idx][contiUreaEqIdx];
R_sum[contiUreaEqIdx] += R3;
maxCoeff[contiUreaEqIdx] = std::max(maxCoeff[contiUreaEqIdx],
std::abs(R3) / pvValue);
B_avg[contiBiofilmEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R4 = modelResid[cell_idx][contiBiofilmEqIdx];
R_sum[contiBiofilmEqIdx] += R4;
maxCoeff[contiBiofilmEqIdx] = std::max(maxCoeff[contiBiofilmEqIdx],
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_BLACKOILMODEL_IMPL_HEADER_INCLUDED
Reference in New Issue View Git Blame Copy Permalink