FlowProblemComp.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright 2024, 2026 SINTEF Digital
 
  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 2 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/>.
 
  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef OPM_FLOW_PROBLEM_COMP_HPP
#define OPM_FLOW_PROBLEM_COMP_HPP
 
 
#include <opm/models/ptflash/flashparameters.hh>
#include <opm/simulators/flow/ExplicitInitialSaturations.hpp>
#include <opm/simulators/flow/FlowProblem.hpp>
#include <opm/simulators/flow/OutputCompositionalModule.hpp>
#include <opm/simulators/flow/equil/InitStateEquilComp.hpp>
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
#include <opm/simulators/utils/DeferredLogger.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 
#include <opm/grid/LookUpData.hh>
 
#include <opm/material/constraintsolvers/PTFlashMethod.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
 
#include <opm/material/thermal/EclThermalLawManager.hpp>
 
#include <opm/common/OpmLog/OpmLog.hpp>
 
#include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
 
#include <dune/grid/common/rangegenerators.hh>
 
#include <fmt/format.h>
 
#include <algorithm>
#include <array>
#include <cstdint>
#include <functional>
#include <optional>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace Opm::Parameters {
 
// Treat explicit compositional saturations as given before the hydrocarbon
// flash and rescale SWAT accordingly. By default, keep the input SWAT.
struct NormalizeExplicitSwat { static constexpr bool value = false; };
 
} // namespace Opm::Parameters
 
namespace Opm {
 
template <class TypeTag>
class FlowProblemComp : public FlowProblem<TypeTag>
{
    // TODO: the naming of the Types will be adjusted
    using FlowProblemType = FlowProblem<TypeTag>;
 
    using typename FlowProblemType::Scalar;
    using typename FlowProblemType::Simulator;
    using typename FlowProblemType::GridView;
    using typename FlowProblemType::FluidSystem;
    using typename FlowProblemType::Vanguard;
 
    // might not be needed
    using FlowProblemType::dim;
    using FlowProblemType::dimWorld;
 
    using FlowProblemType::numPhases;
    using FlowProblemType::numComponents;
 
    using FlowProblemType::gasPhaseIdx;
    using FlowProblemType::oilPhaseIdx;
    using FlowProblemType::waterPhaseIdx;
 
    using typename FlowProblemType::Indices;
    using typename FlowProblemType::PrimaryVariables;
    using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
    using typename FlowProblemType::Evaluation;
    using typename FlowProblemType::IntensiveQuantities;
    using typename FlowProblemType::MaterialLaw;
    using typename FlowProblemType::RateVector;
 
    using InitialFluidState = CompositionalFluidState<Scalar, FluidSystem>;
    using EclWriterType = EclWriter<TypeTag, OutputCompositionalModule<TypeTag> >;
 
public:
    using FlowProblemType::porosity;
    using FlowProblemType::pvtRegionIndex;
 
    static void registerParameters()
    {
        FlowProblemType::registerParameters();
        Parameters::SetDefault<Parameters::InputSkipMode>("300");
 
        EclWriterType::registerParameters();
 
        // tighter tolerance is needed for compositional modeling here
        Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1e-7);
        // the CPR solves of compositional systems converge far more reliably
        // when the fine smoother also runs before the pressure correction.
        // Only ISTLSolver and its runtime proxy register the parameter;
        // flowexp_comp uses another linear solver.
        using Backend = GetPropType<TypeTag, Properties::LinearSolverBackend>;
        if constexpr (std::is_same_v<Backend, ISTLSolver<TypeTag>>
                      || std::is_same_v<Backend, ISTLSolverRuntimeOptionProxy<TypeTag>>) {
            Parameters::SetDefault<Parameters::CprPreSmooth>(1);
        }
 
        Parameters::Register<Parameters::NormalizeExplicitSwat>(
            "Treat the explicit initial saturations of compositional runs as "
            "saturations before the hydrocarbon flash: SGAS and the oil "
            "saturation 1 - SWAT - SGAS each hold their input composition at "
            "the vapour or liquid equation-of-state root, and water and the "
            "flashed phases are then rescaled to fill the pore volume. When "
            "disabled (the default), the input SWAT is retained.");
    }
 
    Opm::CompositionalConfig::EOSType getEosType() const
    {
        auto& simulator = this->simulator();
        const auto& eclState = simulator.vanguard().eclState();
        return eclState.compositionalConfig().eosType(0);
    }
 
    explicit FlowProblemComp(Simulator& simulator)
        : FlowProblemType(simulator)
    {
        eclWriter_ = std::make_unique<EclWriterType>(simulator);
        enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
    }
 
    void finishInit()
    {
        FlowProblemType::finishInit();
 
        auto& simulator = this->simulator();
        const bool transmissibilitiesFinished =
            this->prepareTransmissibilityOutput_(*eclWriter_, enableEclOutput_);
        if (!transmissibilitiesFinished) {
            this->finishTransmissibilities_();
        }
 
        const auto& eclState = simulator.vanguard().eclState();
        const auto& schedule = simulator.vanguard().schedule();
        this->initializeSimulatorTime_();
 
        this->initFluidSystem_();
        if (eclState.getSimulationConfig().rock_config().active()) {
            throw std::runtime_error("Rock compaction (ROCKCOMP) is not supported "
                                     "by compositional modeling yet");
        }
        this->initializeModelProperties_();
 
        // write the static output files (EGRID, INIT)
        if (enableEclOutput_) {
            eclWriter_->writeInit();
        }
 
        const auto& initconfig = eclState.getInitConfig();
        if (initconfig.restartRequested())
            readEclRestartSolution_();
        else
            this->readInitialCondition_();
 
        FlowProblemType::updatePffDofData_();
 
        if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
            const auto& vanguard = this->simulator().vanguard();
            const auto& gridView = vanguard.gridView();
            int numElements = gridView.size(/*codim=*/0);
            this->polymer_.maxAdsorption.resize(numElements, 0.0);
        }
 
        /* readBoundaryConditions_();
 
        // compute and set eq weights based on initial b values
        computeAndSetEqWeights_();
 
        if (enableDriftCompensation_) {
            drift_.resize(this->model().numGridDof());
            drift_ = 0.0;
        } */
 
        // TODO: check wether the following can work with compostional
        if (this->enableVtkOutput_() && eclState.getIOConfig().initOnly()) {
            simulator.setTimeStepSize(0.0);
            FlowProblemType::writeOutput(true);
        }
 
        // after finishing the initialization and writing the initial solution, we move
        // to the first "real" episode/report step
        // for restart the episode index and start is already set
        if (!initconfig.restartRequested()) {
            simulator.startNextEpisode(schedule.seconds(1));
            simulator.setEpisodeIndex(0);
            simulator.setTimeStepIndex(0);
        }
    }
 
    void initialSolutionApplied() override
    {
        FlowProblemType::initialSolutionApplied();
 
        // Defaulted threshold pressures are taken from the initial solution.
        this->thresholdPressures_.finishInit();
 
        reportTwoPhaseEquilCells_();
    }
 
    void endTimeStep() override
    {
        FlowProblemType::endTimeStep();
 
        // after the solution is updated, the values in output module also needs to be updated
        this->eclWriter_->mutableOutputModule().invalidateLocalData();
 
        // For CpGrid with LGRs, ecl/vtk output is not supported yet.
        const auto& grid = this->simulator().vanguard().gridView().grid();
 
        using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
        constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
        if (!isCpGrid || (grid.maxLevel() == 0)) {
            this->eclWriter_->evalSummaryState(! this->episodeWillBeOver());
        }
    }
 
    void writeReports(const SimulatorTimer& timer) {
        if (enableEclOutput_){
            eclWriter_->writeReports(timer);
        }
    }
 
    void writeOutput(bool verbose) override
    {
        FlowProblemType::writeOutput(verbose);
 
        if (! this->enableEclOutput_) {
            return;
        }
 
        const auto isSubStep = !this->episodeWillBeOver();
 
        if (!isSubStep || Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
            auto localCellData = data::Solution {};
 
            this->eclWriter_->writeOutput(std::move(localCellData), isSubStep,
                                          this->simulator().vanguard().schedule()
                                         .exitStatus().has_value());
        }
    }
 
    template <class Context>
    void boundary(BoundaryRateVector& values,
                  const Context& context,
                  unsigned spaceIdx,
                  unsigned /* timeIdx */) const
    {
        OPM_TIMEBLOCK_LOCAL(eclProblemBoundary, Subsystem::Assembly);
        if (!context.intersection(spaceIdx).boundary())
            return;
 
        values.setNoFlow();
 
        if (this->nonTrivialBoundaryConditions()) {
            throw std::logic_error("boundary condition is not supported by compostional modeling yet");
        }
    }
 
    template <class Context>
    void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
    {
        const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
        const auto& initial_fs = initialFluidStates_[globalDofIdx];
        Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
        for (unsigned p = 0; p < numPhases; ++p) { // TODO: assuming the phaseidx continuous
            // pressure
            fs.setPressure(p, initial_fs.pressure(p));
 
            // saturation
            fs.setSaturation(p, initial_fs.saturation(p));
 
            // temperature
            fs.setTemperature(initial_fs.temperature(p));
        }
 
 
        if (!zmf_initialization_) {
            for (unsigned p = 0; p < numPhases; ++p) {
                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                    fs.setMoleFraction(p, compIdx, initial_fs.moleFraction(p, compIdx));
                }
            }
 
            {
                const auto& eos_type = getEosType();
                typename FluidSystem::template ParameterCache<Scalar> paramCache(eos_type);
                paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
                paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
                fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
                fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
            }
            // determine the component fractions
            Dune::FieldVector<Scalar, numComponents> z(0.0);
            Scalar sumMoles = 0.0;
            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                if (Indices::waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)){
                    continue;
                }
                const auto saturation = fs.saturation(phaseIdx);
                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                    Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation;
                    tmp = max(tmp, 1e-8);
                    z[compIdx] += tmp;
                    sumMoles += tmp;
                }
            }
            z /= sumMoles;
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                fs.setMoleFraction(compIdx, z[compIdx]);
            }
        } else {
            // TODO: should we normalize the input?
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                fs.setMoleFraction(compIdx, initial_fs.moleFraction(compIdx));
            }
        }
 
        // Set initial K and L
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
            const auto& Ktmp = fs.wilsonK_(compIdx);
            fs.setKvalue(compIdx, Ktmp);
        }
 
        const Scalar& Ltmp = -1.0;
        fs.setLvalue(Ltmp);
 
        values.assignNaive(fs);
    }
 
    void addToSourceDense(RateVector&, unsigned, unsigned) const override
    {
        // we do nothing for now
    }
 
    const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
    { return initialFluidStates_[globalDofIdx]; }
 
    std::vector<InitialFluidState>& initialFluidStates()
    { return initialFluidStates_; }
 
    const std::vector<InitialFluidState>& initialFluidStates() const
    { return initialFluidStates_; }
 
    const EclWriterType& eclWriter() const
    { return *eclWriter_; }
 
    EclWriterType& eclWriter()
    { return *eclWriter_; }
 
    void setSubStepReport(const SimulatorReportSingle& report)
    { return eclWriter_->setSubStepReport(report); }
 
    void setSimulationReport(const SimulatorReport& report)
    { return eclWriter_->setSimulationReport(report); }
 
    void finalizeOutput()
    {
        OPM_TIMEBLOCK(finalizeOutput);
        eclWriter_.reset();
    }
 
    // TODO: do we need this one?
    template<class Serializer>
    void serializeOp(Serializer& serializer)
    {
        serializer(static_cast<FlowProblemType&>(*this));
        serializer(*eclWriter_);
    }
protected:
 
    void updateExplicitQuantities_(int /* episodeIdx*/, int /* timeStepSize */, bool /* first_step_after_restart */) override
    {
        // we do nothing here for now
    }
 
    void readEquilInitialCondition_() override
    {
        const auto& simulator = this->simulator();
        const auto& vanguard = simulator.vanguard();
        const auto& eclState = vanguard.eclState();
 
        // Zero-based equilibration region of every cell (EQLNUM, or region 0).
        // EQLNUM is given on the unrefined input grid while the equilibration
        // works on the leaf grid, so a refined cell takes its parent's region;
        // without local refinement the mapping is the identity.
        std::vector<int> eqlnum(this->model().numGridDof(), 0);
        if (eclState.fieldProps().has_int("EQLNUM")) {
            const LookUpData<typename GridView::Grid, GridView>
                lookUpData(vanguard.gridView());
            eqlnum = lookUpData.template assignFieldPropsIntOnLeaf<int>(
                eclState.fieldProps(), "EQLNUM", /*needsTranslation=*/true);
        }
 
        // The water saturation endpoints are selected by SATNUM and scaled per
        // cell, so they are read here rather than from the equilibration region.
        std::vector<Scalar> connateWater, maxWater;
        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
            const auto numDof = this->model().numGridDof();
            connateWater.resize(numDof);
            maxWater.resize(numDof);
            for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
                const auto& eps = this->materialLawManager()
                    .oilWaterScaledEpsInfoDrainage(static_cast<int>(dofIdx));
                connateWater[dofIdx] = eps.Swl;
                maxWater[dofIdx] = eps.Swu;
            }
        }
 
        EQUIL::Comp::InitialStateComputer<FluidSystem> initialState(
            eclState,
            getEosType(),
            vanguard.cellCenterDepths(),
            eqlnum,
            this->pvtnum_,
            vanguard.gridView().comm(),
            this->gravity()[dimWorld - 1],
            this->numPressurePointsEquil(),
            connateWater,
            maxWater);
 
        initialFluidStates_ = std::move(initialState.fluidStates());
        twoPhaseCheckRegions_ = std::move(initialState.twoPhaseCheckRegions());
        // The current flash has one pressure for all phases. Adapt the
        // independent equilibrium columns here, using the water column below
        // the contact. Once phase-pressure differences are supported by the
        // flash, the equilibrated phase pressures can be used directly.
        for (std::size_t cell = 0; cell < initialFluidStates_.size(); ++cell) {
            for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
                initialFluidStates_[cell].setPressure(phaseIdx,
                                                     initialState.referencePressures()[cell]);
            }
        }
        // The primary variables are formed from the total composition; see initial().
        zmf_initialization_ = true;
    }
 
    void readEclRestartSolution_()
    {
        throw std::logic_error("Restarting is not supported by compositional modeling yet");
    }
 
    void readExplicitInitialCondition_() override
    {
        readExplicitInitialConditionCompositional_();
    }
 
    void readExplicitInitialConditionCompositional_()
    {
        const auto& simulator = this->simulator();
        const auto& vanguard = simulator.vanguard();
        const auto& eclState = vanguard.eclState();
        const auto& fp = eclState.fieldProps();
        const bool has_pressure = fp.has_double("PRESSURE");
        if (!has_pressure)
            throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
                                     "keyword if the model is initialized explicitly");
 
        const bool has_xmf = fp.has_double("XMF");
        const bool has_ymf = fp.has_double("YMF");
        const bool has_zmf = fp.has_double("ZMF");
        if ( !has_zmf && !(has_xmf && has_ymf) ) {
            throw std::runtime_error("The ECL input file requires the presence of ZMF or XMF and YMF "
                                     "keyword if the model is initialized explicitly");
        }
 
        if (has_zmf && (has_xmf || has_ymf)) {
            throw std::runtime_error("The ECL input file can not handle explicit initialization "
                                     "with both ZMF and XMF or YMF");
        }
 
        if (has_xmf != has_ymf) {
            throw std::runtime_error("The ECL input file needs XMF and YMF combined to do the explicit "
                                     "initializtion when using XMF or YMF");
        }
 
        std::size_t numDof = this->model().numGridDof();
 
        initialFluidStates_.resize(numDof);
 
        const bool water_active = FluidSystem::phaseIsActive(waterPhaseIdx);
        const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
        const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
        const std::array activePhases {water_active, oil_active, gas_active};
        constexpr std::array saturationKeywords {"SWAT", "SOIL", "SGAS"};
        std::array<const std::vector<double>*, 3> saturationData {};
        bool allActiveSaturationsSupplied = true;
        for (std::size_t phase = 0; phase < activePhases.size(); ++phase) {
            if (!activePhases[phase]) {
                continue;
            }
            if (fp.has_double(saturationKeywords[phase])) {
                saturationData[phase] = &fp.get_double(saturationKeywords[phase]);
            } else {
                allActiveSaturationsSupplied = false;
            }
        }
 
        const auto& pressureData = fp.get_double("PRESSURE");
 
        // Without TEMPI in the deck, this is the temperature from RTEMPVD or RTEMP.
        const auto& tempiData = fp.get_double("TEMPI");
 
        const auto cellName = [&vanguard](const std::size_t dofIdx) {
            std::array<int, 3> ijk;
            vanguard.cartesianCoordinate(dofIdx, ijk);
            return fmt::format("cell ({}, {}, {})", ijk[0] + 1, ijk[1] + 1, ijk[2] + 1);
        };
 
        // Invalid input on one rank must not leave the others entering the
        // collective operations in normalizeExplicitSwat_().
        const auto& comm = vanguard.grid().comm();
        DeferredLogger deferredLogger;
        std::vector<bool> ownedCells;
        if (allActiveSaturationsSupplied) {
            ownedCells.resize(numDof, false);
            for (const auto& element : elements(vanguard.gridView(), Dune::Partitions::interior)) {
                ownedCells[this->elementMapper().index(element)] = true;
            }
        }
        std::uint64_t normalizedCellCount = 0;
        double maxNormalizationDeviation = 0.0;
        OPM_BEGIN_PARALLEL_TRY_CATCH();
        for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
            auto& dofFluidState = initialFluidStates_[dofIdx];
            // dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
 
            std::array<std::optional<double>, 3> supplied;
            for (std::size_t phase = 0; phase < saturationData.size(); ++phase) {
                if (saturationData[phase]) {
                    supplied[phase] = (*saturationData[phase])[dofIdx];
                }
            }
            ExplicitInitialSaturations completed;
            try {
                completed = completeExplicitSaturations(supplied, activePhases);
            } catch (const std::invalid_argument& error) {
                throw std::invalid_argument("Invalid explicit saturations in "
                                            + cellName(dofIdx) + ": " + error.what());
            }
            const auto& saturations = completed.saturations;
            if (allActiveSaturationsSupplied && ownedCells[dofIdx]
                && completed.normalizationDeviation > 0.0) {
                ++normalizedCellCount;
                maxNormalizationDeviation
                    = std::max(maxNormalizationDeviation, completed.normalizationDeviation);
            }
 
            const Scalar temperatureLoc = tempiData[dofIdx];
            if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0) {
                throw std::invalid_argument(fmt::format("Invalid temperature in {}: "
                                                       "expected a finite, positive value, got {} K",
                                                       cellName(dofIdx), temperatureLoc));
            }
            dofFluidState.setTemperature(temperatureLoc);
 
            if (gas_active) {
                dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
                                            saturations[2]);
            }
            if (oil_active) {
                dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
                                            saturations[1]);
            }
            if (water_active) {
                dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
                                            saturations[0]);
            }
 
            // set phase pressures
            const Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
 
            // TODO: zero capillary pressure for now
            const std::array<Scalar, numPhases> pc = {0};
            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                if (!FluidSystem::phaseIsActive(phaseIdx))
                    continue;
 
                if (Indices::oilEnabled)
                    dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
                else if (Indices::gasEnabled)
                    dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
                else if (Indices::waterEnabled)
                    // single (water) phase
                    dofFluidState.setPressure(phaseIdx, pressure);
            }
 
            if (has_xmf && has_ymf) {
                const auto& xmfData = fp.get_double("XMF");
                const auto& ymfData = fp.get_double("YMF");
                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                    const std::size_t data_idx = compIdx * numDof + dofIdx;
                    const Scalar xmf = xmfData[data_idx];
                    const Scalar ymf = ymfData[data_idx];
 
                    dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
                    dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
                }
            }
 
            if (has_zmf) {
                zmf_initialization_ = true;
                const auto& zmfData = fp.get_double("ZMF");
                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                    const std::size_t data_idx = compIdx * numDof + dofIdx;
                    const Scalar zmf = zmfData[data_idx];
                    dofFluidState.setMoleFraction(compIdx, zmf);
 
                    if (gas_active) {
                        const auto ymf = (dofFluidState.saturation(FluidSystem::gasPhaseIdx) > 0.) ? zmf : Scalar{0};
                        dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
                    }
                    if (oil_active) {
                        const auto xmf = (dofFluidState.saturation(FluidSystem::oilPhaseIdx) > 0.) ? zmf : Scalar{0};
                        dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
                    }
                }
            }
        }
        OPM_END_PARALLEL_TRY_CATCH_LOG(deferredLogger,
                                     "FlowProblemComp: invalid explicit initial state: ",
                                     comm.rank() == 0, comm);
 
        const auto globalNormalizedCellCount = comm.sum(normalizedCellCount);
        const auto globalMaxDeviation = comm.max(maxNormalizationDeviation);
        if (globalNormalizedCellCount > 0 && comm.rank() == 0) {
            OpmLog::warning(
                "EXPLICIT_SATURATION_NORMALIZATION",
                fmt::format("Normalized explicit saturations in {} cells; "
                            "largest absolute deviation of their input sum from 1: {:.6g}",
                            globalNormalizedCellCount,
                            globalMaxDeviation));
        }
 
        if (water_active && Parameters::Get<Parameters::NormalizeExplicitSwat>()) {
            normalizeExplicitSwat_();
        }
    }
 
    // The explicit saturations describe the cell before the hydrocarbon
    // flash. SGAS and the oil saturation So = 1 - SWAT - SGAS each hold
    // their input composition at the vapour or liquid EOS root, whatever the
    // stable state, and with volume shift applied. After flashing the total,
    // water and the flashed phases are rescaled to fill the pore volume:
    //
    //   Sw = SWAT / (SWAT + n * Vm_flash),  n = So / Vm_liquid + SGAS / Vm_vapour
    void normalizeExplicitSwat_()
    {
        using FlashSolver = GetPropType<TypeTag, Properties::FlashSolver>;
        const auto eos_type = getEosType();
        // Flash as the model does, so that this converges wherever the model's flash does.
        const Scalar flash_tolerance = Parameters::Get<Parameters::FlashTolerance<Scalar>>();
        const int flash_verbosity = Parameters::Get<Parameters::FlashVerbosity>();
        const auto flash_method
            = ptFlashMethodFromString(Parameters::Get<Parameters::FlashTwoPhaseMethod>());
        bool zmf_oil_and_gas = false;
 
        // A failing flash on one rank must not leave the others in the reduction below.
        const auto& comm = this->simulator().vanguard().grid().comm();
        OPM_BEGIN_PARALLEL_TRY_CATCH();
        for (auto& fs : initialFluidStates_) {
            const Scalar sw_in = fs.saturation(FluidSystem::waterPhaseIdx);
            if (sw_in <= 0.0 || sw_in >= 1.0) {
                continue;
            }
            if (!zmf_oil_and_gas && zmf_initialization_
                && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
                && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
                && fs.saturation(FluidSystem::oilPhaseIdx) > 0.0
                && fs.saturation(FluidSystem::gasPhaseIdx) > 0.0) {
                zmf_oil_and_gas = true;
            }
 
            typename FluidSystem::template ParameterCache<Scalar> param_cache(eos_type);
 
            // Hydrocarbon moles per unit pore volume as given. Each given phase holds its
            // input composition, ZMF or XMF/YMF, before the flash, so this is the molar
            // volume of that composition at the phase's own root.
            Dune::FieldVector<Scalar, numComponents> moles(0.0);
            Scalar total_moles = 0.0;
            for (const unsigned phaseIdx : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) {
                const Scalar saturation = fs.saturation(phaseIdx);
                if (!FluidSystem::phaseIsActive(phaseIdx) || saturation <= 0.0) {
                    continue;
                }
                param_cache.updatePhase(fs, phaseIdx);
                const Scalar phase_moles = saturation / param_cache.correctedMolarVolume(phaseIdx);
                for (unsigned c = 0; c < numComponents; ++c) {
                    moles[c] += phase_moles * fs.moleFraction(phaseIdx, c);
                }
                total_moles += phase_moles;
            }
            if (total_moles <= 0.0) {
                continue;
            }
 
            // Overall composition of the given phases, which is ZMF itself with ZMF input.
            // Floor and renormalize it as the model does before its flash.
            Dune::FieldVector<Scalar, numComponents> z(0.0);
            Scalar sum_z = 0.0;
            for (unsigned c = 0; c < numComponents; ++c) {
                z[c] = std::max(moles[c] / total_moles, IntensiveQuantities::compositionFloor);
                sum_z += z[c];
            }
 
            InitialFluidState flash_fs;
            flash_fs.setTemperature(fs.temperature(0));
            flash_fs.setPressure(FluidSystem::oilPhaseIdx, fs.pressure(FluidSystem::oilPhaseIdx));
            flash_fs.setPressure(FluidSystem::gasPhaseIdx, fs.pressure(FluidSystem::gasPhaseIdx));
            for (unsigned c = 0; c < numComponents; ++c) {
                flash_fs.setMoleFraction(c, z[c] / sum_z);
                flash_fs.setKvalue(c, flash_fs.wilsonK_(c));
            }
            flash_fs.setLvalue(-1.0);
            FlashSolver::flash_solve_scalar(
                flash_fs, flash_method, flash_tolerance, eos_type, flash_verbosity);
 
            // A single-phase result has L = 0 or 1, so only one phase counts.
            const Scalar L = flash_fs.L();
            Scalar vm_flash = 0.0;
            if (L > 0.0) {
                param_cache.updatePhase(flash_fs, FluidSystem::oilPhaseIdx);
                vm_flash += L * param_cache.correctedMolarVolume(FluidSystem::oilPhaseIdx);
            }
            if (L < 1.0) {
                param_cache.updatePhase(flash_fs, FluidSystem::gasPhaseIdx);
                vm_flash += (1.0 - L) * param_cache.correctedMolarVolume(FluidSystem::gasPhaseIdx);
            }
 
            const Scalar sw = sw_in / (sw_in + total_moles * vm_flash);
            const Scalar hc_scale = (1.0 - sw) / (1.0 - sw_in);
 
            fs.setSaturation(FluidSystem::waterPhaseIdx, sw);
            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                fs.setSaturation(FluidSystem::gasPhaseIdx,
                                 fs.saturation(FluidSystem::gasPhaseIdx) * hc_scale);
            }
            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                fs.setSaturation(FluidSystem::oilPhaseIdx,
                                 fs.saturation(FluidSystem::oilPhaseIdx) * hc_scale);
            }
        }
        OPM_END_PARALLEL_TRY_CATCH("FlowProblemComp::normalizeExplicitSwat_() failed: ", comm);
 
        if (comm.max(zmf_oil_and_gas ? 1 : 0) > 0 && comm.rank() == 0) {
            OpmLog::warning("ZMF is given with both oil and gas saturations; both phases take "
                            "the ZMF composition when SWAT is normalized.");
        }
    }
 
private:
 
    void reportTwoPhaseEquilCells_()
    {
        if (!twoPhaseCheckRegions_.has_value()) {
            return;
        }
 
        const auto liquidFraction = [this](const std::size_t cell) {
            const auto& intQuants
                = this->model().intensiveQuantities(static_cast<unsigned>(cell), /*timeIdx=*/0);
            return getValue(intQuants.fluidState().L());
        };
        const auto& vanguard = this->simulator().vanguard();
        EQUIL::Comp::reportTwoPhaseCells(*twoPhaseCheckRegions_,
                                         vanguard.cellCenterDepths(),
                                         liquidFraction,
                                         vanguard.eclState().getInitConfig().getEquil().size(),
                                         vanguard.gridView().comm());
        twoPhaseCheckRegions_.reset();
    }
 
    void handleSolventBC(const BCState::BCFace& /* bc */, RateVector& /* rate */) const override
    {
        throw std::logic_error("solvent is disabled for compositional modeling and you're trying to add solvent to BC");
    }
 
    void handlePolymerBC(const BCState::BCFace& /* bc */, RateVector& /* rate */) const override
    {
        throw std::logic_error("polymer is disabled for compositional modeling and you're trying to add polymer to BC");
    }
 
    void handleMicrBC(const BCState::BCFace& /* bc */, RateVector& /* rate */) const override
    {
        throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add microbes to BC");
    }
 
    void handleOxygBC(const BCState::BCFace& /* bc */, RateVector& /* rate */) const override
    {
        throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add oxygen to BC");
    }
 
    void handleUreaBC(const BCState::BCFace& /* bc */, RateVector& /* rate */) const override
    {
        throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add urea to BC");
    }
 
    std::vector<InitialFluidState> initialFluidStates_;
    std::optional<std::vector<int>> twoPhaseCheckRegions_;
 
    bool zmf_initialization_ {false};
 
    bool enableEclOutput_{false};
    std::unique_ptr<EclWriterType> eclWriter_;
};
 
} // namespace Opm
 
#endif // OPM_FLOW_PROBLEM_COMP_HPP