BlackoilModel_impl.hpp

/*
  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 <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::BlackoilModel<TypeTag>::DebugFlags f)
    {
        using F = typename Opm::BlackoilModel<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>
BlackoilModel<TypeTag>::
BlackoilModel(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<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 && 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
BlackoilModel<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
BlackoilModel<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
BlackoilModel<TypeTag>::
nonlinearIterationNewton(const SimulatorTimerInterface& timer,
                         NonlinearSolverType& nonlinear_solver)
{
    // -----------   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
BlackoilModel<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 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(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
BlackoilModel<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
BlackoilModel<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
BlackoilModel<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 BlackoilModel<TypeTag>::MaxSolutionUpdateData
BlackoilModel<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 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>
typename BlackoilModel<TypeTag>::CnvPvSplitData
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());

    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("BlackoilModel::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
BlackoilModel<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
BlackoilModel<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
BlackoilModel<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
BlackoilModel<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("BlackoilModel::convergencePerCell() failed: ",
                               this->grid_.comm());
    rst_conv.updateNewton(convNewt);
}

template <class TypeTag>
ConvergenceReport
BlackoilModel<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 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>
const SimulatorReport&
BlackoilModel<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>&
BlackoilModel<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
BlackoilModel<TypeTag>::
writeNonlinearIterationsPerCell(const std::filesystem::path& odir) const
{
    if (hasNlddSolver()) {
        nlddSolver_->writeNonlinearIterationsPerCell(odir);
    }
}

template <class TypeTag>
void
BlackoilModel<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
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 = 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_BLACKOILMODEL_IMPL_HEADER_INCLUDED