blackoilprimaryvariables.hh

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:
/*
  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 EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH
#define EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH

#include <dune/common/fvector.hh>

#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/utility/gpuDecorators.hpp>

#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>

#include <opm/models/blackoil/blackoilmodules.hpp>
#include <opm/models/blackoil/blackoilmeanings.hh>
#include <opm/models/blackoil/blackoilproperties.hh>

#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>

#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <limits>
#include <stdexcept>
#include <type_traits>

#include <fmt/format.h>

namespace Opm::Parameters {

template<class Scalar>
struct PressureScale
{ static constexpr Scalar value = 1.0; };

} // namespace Opm::Parameters

namespace Opm {

template <class TypeTag, template<class, int> class VectorType = Dune::FieldVector>
class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag, VectorType>
{
    using ParentType = FvBasePrimaryVariables<TypeTag, VectorType>;
    using Implementation = GetPropType<TypeTag, Properties::PrimaryVariables>;

    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
    using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;

    // number of equations
    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };

    // primary variable indices
    static constexpr unsigned waterSwitchIdx = Indices::waterSwitchIdx;
    static constexpr unsigned pressureSwitchIdx = Indices::pressureSwitchIdx;
    static constexpr unsigned compositionSwitchIdx = Indices::compositionSwitchIdx;
    static constexpr unsigned saltConcentrationIdx  = Indices::saltConcentrationIdx;
    static constexpr unsigned solventSaturationIdx  = Indices::solventSaturationIdx;

    static constexpr bool compositionSwitchEnabled =
        Indices::compositionSwitchIdx != std::numeric_limits<unsigned>::max();
    static constexpr bool waterEnabled = Indices::waterEnabled;
    static constexpr bool gasEnabled = Indices::gasEnabled;
    static constexpr bool oilEnabled = Indices::oilEnabled;

    // phase indices from the fluid system
    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };

    // component indices from the fluid system
    enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };

    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();

    static constexpr bool enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();
    static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
    static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
    static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
    static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
    static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
    static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();

    enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
    enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
    enum { enableMICP = Indices::enableMICP };
    enum { gasCompIdx = FluidSystem::gasCompIdx };
    enum { waterCompIdx = FluidSystem::waterCompIdx };
    enum { oilCompIdx = FluidSystem::oilCompIdx };

    using Toolbox = MathToolbox<Evaluation>;
    using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
    using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
    using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
    using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
    using BrineModule = BlackOilBrineModule<TypeTag, enableBrine>;
    using BioeffectsModule = BlackOilBioeffectsModule<TypeTag, enableBioeffects>;

    static_assert(numPhases == 3, "The black-oil model assumes three phases!");
    static_assert(numComponents == 3, "The black-oil model assumes three components!");

public:
    // We are instantiating this class with different TypeTags and VectorTypes,
    // but we still want to be able to copy between them. Therefore, we
    // need to define the same type aliases in all specializations.
    using WaterMeaning = ::Opm::BlackOil::WaterMeaning;
    using PressureMeaning = ::Opm::BlackOil::PressureMeaning;
    using GasMeaning = ::Opm::BlackOil::GasMeaning;
    using BrineMeaning = ::Opm::BlackOil::BrineMeaning;
    using SolventMeaning = ::Opm::BlackOil::SolventMeaning;

    // Allow conversion between different vector types
    template<class OtherTypeTag, template<class, int> class OtherVectorType>
    friend class BlackOilPrimaryVariables;

    template <class OtherTypeTag, template <class, int> class OtherVectorType>
    explicit OPM_HOST_DEVICE
    BlackOilPrimaryVariables(const BlackOilPrimaryVariables<OtherTypeTag, OtherVectorType>& other)
        : ParentType(other)
        , primaryVarsMeaningWater_(other.primaryVarsMeaningWater_)
        , primaryVarsMeaningPressure_(other.primaryVarsMeaningPressure_)
        , primaryVarsMeaningGas_(other.primaryVarsMeaningGas_)
        , primaryVarsMeaningBrine_(other.primaryVarsMeaningBrine_)
        , primaryVarsMeaningSolvent_(other.primaryVarsMeaningSolvent_)
        , pvtRegionIdx_(other.pvtRegionIdx_)
        , pcFactor_(other.pcFactor_)
    {
    }

    OPM_HOST_DEVICE BlackOilPrimaryVariables()
    {
        Valgrind::SetUndefined(*this);
        pvtRegionIdx_ = 0;
    }

    BlackOilPrimaryVariables(const BlackOilPrimaryVariables& value) = default;

    static BlackOilPrimaryVariables serializationTestObject()
    {
        BlackOilPrimaryVariables result;
        result.pvtRegionIdx_ = 1;
        result.primaryVarsMeaningBrine_ = BrineMeaning::Sp;
        result.primaryVarsMeaningGas_ = GasMeaning::Rv;
        result.primaryVarsMeaningPressure_ = PressureMeaning::Pg;
        result.primaryVarsMeaningWater_ = WaterMeaning::Rsw;
        result.primaryVarsMeaningSolvent_ = SolventMeaning::Ss;
        for (std::size_t i = 0; i < result.size(); ++i) {
            result[i] = i + 1;
        }

        return result;
    }

    static void init()
    {
        // TODO: these parameters have undocumented non-trivial dependencies
        pressureScale_ = Parameters::Get<Parameters::PressureScale<Scalar>>();

// currently for GPU we do not have a support for the pressureScale_ variable,
// so we issue a warning if the user has set it to something else than 1.0
// if we are building with GPU support. Note that we can not only test if we are compiling
// with a GPU compiler, because we might initialize the parameters from a .cpp compilation unit
// compiled by a normal C++ compiler.
#if HAVE_CUDA
        if (pressureScale_ != Scalar {1.0}) {
            OpmLog::warning(fmt::format(
                "Using a pressure scaling different from 1.0 is not supported "
                "when running with GPU support. We have detected that you are compiling with GPU support, but we can "
                "not detect whether you are running with GPU support for the assembly or property evaluation. If you "
                "are doing property evaluation or assembly on the GPU, pressure scaling will be ignored."
                "Read value of pressure scale: {}",
                pressureScale_));
        }
#endif
    }

    static void registerParameters()
    {
        Parameters::Register<Parameters::PressureScale<Scalar>>
            ("Scaling of pressure primary variable");
    }

    OPM_HOST_DEVICE Evaluation
    makeEvaluation(unsigned varIdx, unsigned timeIdx,
                   LinearizationType linearizationType = LinearizationType()) const
    {
        const Scalar scale = varIdx == pressureSwitchIdx ? this->getPressureScale() : Scalar{1.0};
        if (std::is_same_v<Evaluation, Scalar>) {
            return (*this)[varIdx] * scale; // finite differences
        }
        else {
            // automatic differentiation
            if (timeIdx == linearizationType.time) {
                return Toolbox::createVariable((*this)[varIdx], varIdx) * scale;
            }
            else {
                return Toolbox::createConstant((*this)[varIdx]) * scale;
            }
        }
    }

     OPM_HOST_DEVICE void setPvtRegionIndex(unsigned value)
    { pvtRegionIdx_ = static_cast<unsigned short>(value); }

    OPM_HOST_DEVICE unsigned pvtRegionIndex() const
    { return pvtRegionIdx_; }

    OPM_HOST_DEVICE WaterMeaning primaryVarsMeaningWater() const
    { return primaryVarsMeaningWater_; }

    OPM_HOST_DEVICE void setPrimaryVarsMeaningWater(WaterMeaning newMeaning)
    { primaryVarsMeaningWater_ = newMeaning; }

    OPM_HOST_DEVICE PressureMeaning primaryVarsMeaningPressure() const
    { return primaryVarsMeaningPressure_; }

     OPM_HOST_DEVICE void setPrimaryVarsMeaningPressure(PressureMeaning newMeaning)
    { primaryVarsMeaningPressure_ = newMeaning; }

     OPM_HOST_DEVICE GasMeaning primaryVarsMeaningGas() const
    { return primaryVarsMeaningGas_; }

    OPM_HOST_DEVICE void setPrimaryVarsMeaningGas(GasMeaning newMeaning)
    { primaryVarsMeaningGas_ = newMeaning; }

    OPM_HOST_DEVICE BrineMeaning primaryVarsMeaningBrine() const
    { return primaryVarsMeaningBrine_; }

    OPM_HOST_DEVICE void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning)
    { primaryVarsMeaningBrine_ = newMeaning; }

    OPM_HOST_DEVICE SolventMeaning primaryVarsMeaningSolvent() const
    { return primaryVarsMeaningSolvent_; }

     OPM_HOST_DEVICE void setPrimaryVarsMeaningSolvent(SolventMeaning newMeaning)
    { primaryVarsMeaningSolvent_ = newMeaning; }

    template <class FluidState>
    OPM_HOST_DEVICE void assignNaive(const FluidState& fluidState)
    {
        using ConstEvaluation = std::remove_reference_t<typename FluidState::ValueType>;
        using FsEvaluation = std::remove_const_t<ConstEvaluation>;
        using FsToolbox = MathToolbox<FsEvaluation>;

        const bool gasPresent =
            fluidState.fluidSystem().phaseIsActive(gasPhaseIdx)
                ? fluidState.saturation(gasPhaseIdx) > 0.0
                : false;
        const bool oilPresent =
            fluidState.fluidSystem().phaseIsActive(oilPhaseIdx)
                ? fluidState.saturation(oilPhaseIdx) > 0.0
                : false;
        const bool waterPresent =
            fluidState.fluidSystem().phaseIsActive(waterPhaseIdx)
                ? fluidState.saturation(waterPhaseIdx) > 0.0
                : false;
        const auto& saltSaturation =
            BlackOil::getSaltSaturation_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
        const bool precipitatedSaltPresent = enableSaltPrecipitation ? saltSaturation > 0.0 : false;
        const bool oneActivePhases = fluidState.fluidSystem().numActivePhases() == 1;
        // deal with the primary variables for the energy extension
        if constexpr (enableFullyImplicitThermal) {
            EnergyModule::assignPrimaryVars(*this, fluidState);
        }

        // Determine the meaning of the pressure primary variables
        // Depending on the phases present, this variable is either interpreted as the
        // pressure of the oil phase, gas phase (if no oil) or water phase (if only water)
        if (gasPresent && fluidState.fluidSystem().enableVaporizedOil() && !oilPresent) {
            primaryVarsMeaningPressure_ = PressureMeaning::Pg;
        }
        else if (fluidState.fluidSystem().phaseIsActive(oilPhaseIdx)) {
            primaryVarsMeaningPressure_ = PressureMeaning::Po;
        }
        else if (waterPresent && fluidState.fluidSystem().enableDissolvedGasInWater() && !gasPresent) {
            primaryVarsMeaningPressure_ = PressureMeaning::Pw;
        }
        else if (fluidState.fluidSystem().phaseIsActive(gasPhaseIdx)) {
            primaryVarsMeaningPressure_ = PressureMeaning::Pg;
        }
        else {
            assert(fluidState.fluidSystem().phaseIsActive(waterPhaseIdx));
            primaryVarsMeaningPressure_ = PressureMeaning::Pw;
        }

        // Determine the meaning of the water primary variables
        // Depending on the phases present, this variable is either interpreted as
        // water saturation or vapporized water in the gas phase
        // For two-phase gas-oil models and one-phase case the variable is disabled.
        if (waterPresent && gasPresent) {
            primaryVarsMeaningWater_ = WaterMeaning::Sw;
        }
        else if (gasPresent && fluidState.fluidSystem().enableVaporizedWater()) {
            primaryVarsMeaningWater_ = WaterMeaning::Rvw;
        }
        else if (waterPresent && fluidState.fluidSystem().enableDissolvedGasInWater()) {
            primaryVarsMeaningWater_ = WaterMeaning::Rsw;
        }
        else if (fluidState.fluidSystem().phaseIsActive(waterPhaseIdx) && !oneActivePhases) {
            primaryVarsMeaningWater_ = WaterMeaning::Sw;
        }
        else {
            primaryVarsMeaningWater_ = WaterMeaning::Disabled;
        }

        // Determine the meaning of the gas primary variables
        // Depending on the phases present, this variable is either interpreted as the
        // saturation of the gas phase, as the fraction of the gas component in the oil
        // phase (Rs) or as the  fraction of the oil component (Rv) in the gas phase.
        // For two-phase water-oil and water-gas models and one-phase case the variable is disabled.
        if (gasPresent && oilPresent) {
            primaryVarsMeaningGas_ = GasMeaning::Sg;
        }
        else if (oilPresent && fluidState.fluidSystem().enableDissolvedGas()) {
            primaryVarsMeaningGas_ = GasMeaning::Rs;
        }
        else if (gasPresent && fluidState.fluidSystem().enableVaporizedOil()) {
            primaryVarsMeaningGas_ = GasMeaning::Rv;
        }
        else if (fluidState.fluidSystem().phaseIsActive(gasPhaseIdx) && fluidState.fluidSystem().phaseIsActive(oilPhaseIdx)) {
            primaryVarsMeaningGas_ = GasMeaning::Sg;
        }
        else {
            primaryVarsMeaningGas_ = GasMeaning::Disabled;
        }

        // Determine the meaning of the brine primary variables
        if constexpr (enableSaltPrecipitation) {
            if (precipitatedSaltPresent) {
                primaryVarsMeaningBrine_ = BrineMeaning::Sp;
            }
            else {
                primaryVarsMeaningBrine_ = BrineMeaning::Cs;
            }
        }
        else {
            primaryVarsMeaningBrine_ = BrineMeaning::Disabled;
        }

        // assign the actual primary variables
        switch (primaryVarsMeaningPressure()) {
        case PressureMeaning::Po:
            this->setScaledPressure_(FsToolbox::value(fluidState.pressure(oilPhaseIdx)));
            break;
        case PressureMeaning::Pg:
            this->setScaledPressure_(FsToolbox::value(fluidState.pressure(gasPhaseIdx)));
            break;
        case PressureMeaning::Pw:
            this->setScaledPressure_(FsToolbox::value(fluidState.pressure(waterPhaseIdx)));
            break;
        default:
            OPM_THROW(std::logic_error, "No valid primary variable selected for pressure");
        }

        switch (primaryVarsMeaningWater()) {
        case WaterMeaning::Sw:
        {
            (*this)[waterSwitchIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
            break;
        }
        case WaterMeaning::Rvw:
        {
            const auto& rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
            (*this)[waterSwitchIdx] = rvw;
            break;
        }
        case WaterMeaning::Rsw:
        {
            const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
            (*this)[waterSwitchIdx] = Rsw;
            break;
        }
        case WaterMeaning::Disabled:
            break;
        default:
            OPM_THROW(std::logic_error, "No valid primary variable selected for water");
        }
        switch (primaryVarsMeaningGas()) {
        case GasMeaning::Sg:
            (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
            break;
        case GasMeaning::Rs:
        {
            const auto& rs = BlackOil::getRs_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
            (*this)[compositionSwitchIdx] = rs;
            break;
        }
        case GasMeaning::Rv:
        {
            const auto& rv = BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
            (*this)[compositionSwitchIdx] = rv;
            break;
        }
        case GasMeaning::Disabled:
            break;
        default:
            OPM_THROW(std::logic_error, "No valid primary variable selected for composision");
        }
    }

    bool adaptPrimaryVariables(const Problem& problem,
                               unsigned globalDofIdx,
                               [[maybe_unused]] Scalar swMaximum,
                               Scalar thresholdWaterFilledCell, Scalar eps = 0.0)
    {
        // this function accesses quite a few black-oil specific low-level functions
        // directly for better performance (instead of going the canonical way through
        // the IntensiveQuantities). The reason is that most intensive quantities are not
        // required to be able to decide if the primary variables needs to be switched or
        // not, so it would be a waste to compute them.

        // Both the primary variable meaning of water and gas are disabled i.e.
        // It is a one-phase case and we no variable meaning switch is needed.
        if (primaryVarsMeaningWater() == WaterMeaning::Disabled &&
            primaryVarsMeaningGas() == GasMeaning::Disabled)
        {
            return false;
        }

        // Read the current saturation from the primary variables
        Scalar sw = 0.0;
        Scalar sg = 0.0;
        Scalar saltConcentration = 0.0;
        const Scalar& T = asImp_().temperature_(problem, globalDofIdx);
        if (primaryVarsMeaningWater() == WaterMeaning::Sw) {
            sw = (*this)[waterSwitchIdx];
        }
        if (primaryVarsMeaningGas() == GasMeaning::Sg) {
            sg = (*this)[compositionSwitchIdx];
        }
        if (primaryVarsMeaningWater() == WaterMeaning::Rsw) {
            sw = 1.0;
        }

        if (primaryVarsMeaningGas() == GasMeaning::Disabled && gasEnabled) {
            sg = 1.0 - sw; // water + gas case
        }

        // if solid phase disappeares:  Sp (Solid salt saturation) -> Cs (salt concentration)
        // if solid phase appears: Cs (salt concentration) ->  Sp (Solid salt saturation)
        if constexpr (enableSaltPrecipitation) {
            const Scalar saltSolubility = BrineModule::saltSol(pvtRegionIndex());
            if (primaryVarsMeaningBrine() == BrineMeaning::Sp) {
                saltConcentration = saltSolubility;
                const Scalar saltSat = (*this)[saltConcentrationIdx];
                if (saltSat < -eps) { // precipitated salt dissappears
                    setPrimaryVarsMeaningBrine(BrineMeaning::Cs);
                    (*this)[saltConcentrationIdx] = saltSolubility; // set salt concentration to solubility limit
                }
            }
            else if (primaryVarsMeaningBrine() == BrineMeaning::Cs) {
                saltConcentration = (*this)[saltConcentrationIdx];
                if (saltConcentration > saltSolubility + eps) { // salt concentration exceeds solubility limit
                    setPrimaryVarsMeaningBrine(BrineMeaning::Sp);
                    (*this)[saltConcentrationIdx] = 0.0;
                }
            }
        }

        // if solvent saturation disappeares:  Ss (Solvent saturation) -> Rsolw (solvent dissolved in water)
        // if solvent saturation appears: Rsolw (solvent dissolved in water) ->  Ss (Solvent saturation)
        // Scalar rsolw = 0.0; // not needed at the moment since we dont allow for vapwat in combination with rsolw
        if constexpr (enableSolvent) {
            if (SolventModule::isSolubleInWater()) {
                const Scalar p = (*this)[pressureSwitchIdx]; // cap-pressure?
                const Scalar solLimit =
                    SolventModule::solubilityLimit(pvtRegionIndex(), T , p, saltConcentration);
                if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
                    const Scalar solSat = (*this)[solventSaturationIdx];
                    if (solSat < -eps) { // solvent dissappears
                        setPrimaryVarsMeaningSolvent(SolventMeaning::Rsolw);
                        (*this)[solventSaturationIdx] = solLimit; // set rsolw to solubility limit
                    }
                }
                else if (primaryVarsMeaningSolvent() == SolventMeaning::Rsolw) {
                    const Scalar rsolw = (*this)[solventSaturationIdx];
                    if (rsolw > solLimit + eps) { // solvent appears as phase
                        setPrimaryVarsMeaningSolvent(SolventMeaning::Ss);
                        (*this)[solventSaturationIdx] = 0.0;
                    }
                }
            }
        }

        // keep track if any meaning has changed
        bool changed = false;

        // Special case for cells with almost only water
        // for these cells both saturations (if the phase is enabled) is used
        // to avoid singular systems.
        // If dissolved gas in water is enabled we shouldn't enter
        // here but instead switch to Rsw as primary variable
        // as sw >= 1.0 -> gas <= 0 (i.e. gas phase disappears)
        if (sw >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {
            // make sure water saturations does not exceed sw_maximum. Default to 1.0
            if constexpr (waterEnabled) {
                (*this)[Indices::waterSwitchIdx] = std::min(swMaximum, sw);
                assert(primaryVarsMeaningWater() == WaterMeaning::Sw);
            }
            // the hydrocarbon gas saturation is set to 0.0
            if constexpr (compositionSwitchEnabled) {
                (*this)[Indices::compositionSwitchIdx] = 0.0;
            }

            changed = primaryVarsMeaningGas() != GasMeaning::Sg;
            if (changed) {
                if constexpr (compositionSwitchEnabled) {
                    setPrimaryVarsMeaningGas(GasMeaning::Sg);
                }

                // use water pressure?
            }
            return changed;
        }

        pcFactor_ = 1.0;
        if constexpr (enableBrine) {
            if (BrineModule::hasPcfactTables() && primaryVarsMeaningBrine() == BrineMeaning::Sp) {
                unsigned satnumRegionIdx = problem.satnumRegionIndex(globalDofIdx);
                Scalar Sp = saltConcentration_();
                Scalar porosityFactor  = std::min(1.0 - Sp, 1.0); //phi/phi_0
                const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
                pcFactor_ = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
            }
        }
        else if constexpr (enableBioeffects) {
            if (BioeffectsModule::hasPcfactTables() && problem.referencePorosity(globalDofIdx, 0) > 0) {
                unsigned satnumRegionIdx = problem.satnumRegionIndex(globalDofIdx);
                Scalar Sb = biofilmVolumeFraction_() /
                            problem.referencePorosity(globalDofIdx, 0);
                Scalar porosityFactor  = std::min(1.0 - Sb, 1.0); //phi/phi_0
                const auto& pcfactTable = BioeffectsModule::pcfactTable(satnumRegionIdx);
                pcFactor_ = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
            }
        }

        switch (primaryVarsMeaningWater()) {
            case WaterMeaning::Sw:
            {
                // if water phase disappeares:  Sw (water saturation) -> Rvw (fraction of water in gas phase)
                if (sw < -eps && sg > eps && FluidSystem::enableVaporizedWater()) {
                    Scalar p = this->pressure_();
                    if (primaryVarsMeaningPressure() == PressureMeaning::Po) {
                        std::array<Scalar, numPhases> pC{};
                        const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
                        const Scalar so = 1.0 - sg - solventSaturation_();
                        computeCapillaryPressures_(pC, so, sg + solventSaturation_(), /*sw=*/ 0.0, matParams);
                        p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
                    }
                    const Scalar rvwSat =
                        FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
                                                                               T,
                                                                               p,
                                                                               saltConcentration);
                    setPrimaryVarsMeaningWater(WaterMeaning::Rvw);
                    (*this)[Indices::waterSwitchIdx] = rvwSat; // primary variable becomes Rvw
                    changed = true;
                    break;
                }
                // if gas phase disappeares:  Sw (water saturation) -> Rsw (fraction of gas in water phase)
                // and Pg (gas pressure) -> Pw ( water pressure)
                if (sg < -eps && sw > eps && FluidSystem::enableDissolvedGasInWater()) {
                    const Scalar pg = this->pressure_();
                    assert(primaryVarsMeaningPressure() == PressureMeaning::Pg);
                    std::array<Scalar, numPhases> pC = { 0.0 };
                    const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
                    const Scalar so = 1.0 - sw - solventSaturation_();
                    computeCapillaryPressures_(pC, so,  /*sg=*/ 0.0, sw, matParams);
                    const Scalar pw = pg + pcFactor_ * (pC[waterPhaseIdx] - pC[gasPhaseIdx]);
                    const Scalar rswSat =
                        FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
                                                                              T,
                                                                              pw,
                                                                              saltConcentration);
                    setPrimaryVarsMeaningWater(WaterMeaning::Rsw);
                    const Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
                    (*this)[Indices::waterSwitchIdx] = std::min(rswSat, rswMax); //primary variable becomes Rsw
                    setPrimaryVarsMeaningPressure(PressureMeaning::Pw);
                    this->setScaledPressure_(pw);
                    changed = true;
                    break;
                }
                break;
            }
            case WaterMeaning::Rvw:
            {
                const Scalar& rvw = (*this)[waterSwitchIdx];
                Scalar p = this->pressure_();
                if (primaryVarsMeaningPressure() == PressureMeaning::Po) {
                    std::array<Scalar, numPhases> pC{};
                    const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
                    const Scalar so = 1.0 - sg - solventSaturation_();
                    computeCapillaryPressures_(pC, so, sg + solventSaturation_(), /*sw=*/ 0.0, matParams);
                    p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
                }
                const Scalar rvwSat =
                    FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
                                                                           T,
                                                                           p,
                                                                           saltConcentration);
                // if water phase appears: Rvw (fraction of water in gas phase) -> Sw (water saturation)
                if (rvw > rvwSat * (1.0 + eps)) {
                    setPrimaryVarsMeaningWater(WaterMeaning::Sw);
                    (*this)[Indices::waterSwitchIdx] = 0.0; // water saturation
                    changed = true;
                }
                break;
            }
            case WaterMeaning::Rsw:
            {
                // Gas phase not present. The hydrocarbon gas phase
                // appears as soon as more of the gas component is present in the water phase
                // than what saturated water can hold.
                const Scalar& pw = this->pressure_();
                assert(primaryVarsMeaningPressure() == PressureMeaning::Pw);
                const Scalar rswSat =
                    FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
                                                                          T,
                                                                          pw,
                                                                          saltConcentration);

                const Scalar rsw = (*this)[Indices::waterSwitchIdx];
                const Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
                if (rsw > std::min(rswSat, rswMax)) {
                    // the gas phase appears, i.e., switch the primary variables to WaterMeaning::Sw
                    setPrimaryVarsMeaningWater(WaterMeaning::Sw);
                    (*this)[Indices::waterSwitchIdx] = 1.0; // hydrocarbon water saturation
                    setPrimaryVarsMeaningPressure(PressureMeaning::Pg);
                    std::array<Scalar, numPhases> pC{};
                    const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
                    computeCapillaryPressures_(pC, /*so=*/ 0.0,  /*sg=*/ 0.0, /*sw=*/ 1.0, matParams);
                    const Scalar pg = pw + pcFactor_ * (pC[gasPhaseIdx] - pC[waterPhaseIdx]);
                    this->setScaledPressure_(pg);
                    changed = true;
                }
                break;
            }
            case WaterMeaning::Disabled:
                break;
            default:
                throw std::logic_error("No valid primary variable selected for water");
        }

        // if gas phase disappeares:  Sg (gas saturation) -> Rs (fraction of gas in oil phase)
        // if oil phase disappeares:  Sg (gas saturation) -> Rv (fraction of oil in gas phase)
        //                            Po (oil pressure )  -> Pg (gas pressure)

        // if gas phase appears: Rs (fraction of gas in oil phase) -> Sg (gas saturation)
        // if oil phase appears: Rv (fraction of oil in gas phase) -> Sg (gas saturation)
        //                       Pg (gas pressure )                -> Po (oil pressure)
        switch (primaryVarsMeaningGas()) {
            case GasMeaning::Sg:
            {
                const Scalar s = 1.0 - sw - solventSaturation_();
                if (sg < -eps && s > 0.0 && FluidSystem::enableDissolvedGas()) {
                    const Scalar po = this->pressure_();
                    setPrimaryVarsMeaningGas(GasMeaning::Rs);
                    const Scalar soMax = std::max(s, problem.maxOilSaturation(globalDofIdx));
                    const Scalar rsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
                    Scalar rsSat;
                    if constexpr (enableExtbo) {
                        rsSat = ExtboModule::rs(pvtRegionIndex(), po, zFraction_());
                    }
                    else {
                        rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
                                                                                    T,
                                                                                    po,
                                                                                    s,
                                                                                    soMax);
                    }
                    (*this)[Indices::compositionSwitchIdx] = std::min(rsMax, rsSat);
                    changed = true;
                }
                const Scalar so = 1.0 - sw - solventSaturation_() - sg;
                if (so < -eps && sg > 0.0 && FluidSystem::enableVaporizedOil()) {
                    // the oil phase disappeared and some hydrocarbon gas phase is still
                    // present, i.e., switch the primary variables to GasMeaning::Rv.
                    // we only have the oil pressure readily available, but we need the gas
                    // pressure, i.e. we must determine capillary pressure
                    const Scalar po = this->pressure_();
                    std::array<Scalar, numPhases> pC{};
                    const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
                    computeCapillaryPressures_(pC, /*so=*/0.0, sg + solventSaturation_(), sw, matParams);
                    const Scalar pg = po + pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);

                    // we start at the GasMeaning::Rv value that corresponds to that of oil-saturated
                    // hydrocarbon gas
                    setPrimaryVarsMeaningPressure(PressureMeaning::Pg);
                    this->setScaledPressure_(pg);
                    const Scalar soMax = problem.maxOilSaturation(globalDofIdx);
                    const Scalar rvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
                    Scalar rvSat;
                    if constexpr (enableExtbo) {
                        rvSat = ExtboModule::rv(pvtRegionIndex(), pg, zFraction_());
                    }
                    else {
                        rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
                                                                                     T,
                                                                                     pg,
                                                                                     Scalar(0),
                                                                                     soMax);
                    }
                    setPrimaryVarsMeaningGas(GasMeaning::Rv);
                    (*this)[Indices::compositionSwitchIdx] = std::min(rvMax, rvSat);
                    changed = true;
                }
                break;
            }
            case GasMeaning::Rs:
            {
                // Gas phase not present. The hydrocarbon gas phase
                // appears as soon as more of the gas component is present in the oil phase
                // than what saturated oil can hold.
                const Scalar po = this->pressure_();
                const Scalar so = 1.0 - sw - solventSaturation_();
                const Scalar soMax = std::max(so, problem.maxOilSaturation(globalDofIdx));
                const Scalar rsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
                Scalar rsSat;
                if constexpr (enableExtbo) {
                    rsSat = ExtboModule::rs(pvtRegionIndex(), po, zFraction_());
                }
                else {
                    rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
                                                                                T,
                                                                                po,
                                                                                so,
                                                                                soMax);
                }
                const Scalar rs = (*this)[Indices::compositionSwitchIdx];
                if (rs > std::min(rsMax, rsSat * (Scalar{1.0} + eps))) {
                    // the gas phase appears, i.e., switch the primary variables to GasMeaning::Sg
                    setPrimaryVarsMeaningGas(GasMeaning::Sg);
                    (*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
                    changed = true;
                }
                break;
            }
            case GasMeaning::Rv:
            {
                // The oil phase appears as
                // soon as more of the oil component is present in the hydrocarbon gas phase
                // than what saturated gas contains. Note that we use the blackoil specific
                // low-level PVT objects here for performance reasons.
                const Scalar pg = this->pressure_();
                const Scalar soMax = problem.maxOilSaturation(globalDofIdx);
                const Scalar rvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
                Scalar rvSat;
                if constexpr (enableExtbo) {
                    rvSat = ExtboModule::rv(pvtRegionIndex(), pg, zFraction_());
                }
                else {
                    rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
                                                                                 T,
                                                                                 pg,
                                                                                 /*so=*/Scalar(0.0),
                                                                                 soMax);
                }
                const Scalar rv = (*this)[Indices::compositionSwitchIdx];
                if (rv > std::min(rvMax, rvSat * (Scalar{1.0} + eps))) {
                    // switch to phase equilibrium mode because the oil phase appears. here
                    // we also need the capillary pressures to calculate the oil phase
                    // pressure using the gas phase pressure
                    const Scalar sg2 = 1.0 - sw - solventSaturation_();
                    std::array<Scalar, numPhases> pC{};
                    const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
                    computeCapillaryPressures_(pC,
                                            /*so=*/0.0,
                                            /*sg=*/sg2 + solventSaturation_(),
                                            sw,
                                            matParams);
                    const Scalar po = pg + pcFactor_ * (pC[oilPhaseIdx] - pC[gasPhaseIdx]);

                    setPrimaryVarsMeaningGas(GasMeaning::Sg);
                    setPrimaryVarsMeaningPressure(PressureMeaning::Po);
                    this->setScaledPressure_(po);
                    (*this)[Indices::compositionSwitchIdx] = sg2; // hydrocarbon gas saturation
                    changed = true;
                }
                break;
            }
            case GasMeaning::Disabled:
                break;
            default:
                throw std::logic_error("No valid primary variable selected for water");
        }
        return changed;
    }

    OPM_HOST_DEVICE bool chopAndNormalizeSaturations()
    {
        if (primaryVarsMeaningWater() == WaterMeaning::Disabled &&
            primaryVarsMeaningGas() == GasMeaning::Disabled)
        {
            return false;
        }
        Scalar sw = 0.0;
        if (primaryVarsMeaningWater() == WaterMeaning::Sw) {
            sw = (*this)[Indices::waterSwitchIdx];
        }
        Scalar sg = 0.0;
        if (primaryVarsMeaningGas() == GasMeaning::Sg) {
            sg = (*this)[Indices::compositionSwitchIdx];
        }

        Scalar ssol = 0.0;
        if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
            ssol =(*this) [Indices::solventSaturationIdx];
        }

        Scalar so = 1.0 - sw - sg - ssol;
        sw = std::min(std::max(sw, Scalar{0.0}), Scalar{1.0});
        so = std::min(std::max(so, Scalar{0.0}), Scalar{1.0});
        sg = std::min(std::max(sg, Scalar{0.0}), Scalar{1.0});
        ssol = std::min(std::max(ssol, Scalar{0.0}), Scalar{1.0});
        const Scalar st = sw + so + sg + ssol;
        sw = sw / st;
        sg = sg / st;
        ssol = ssol / st;
        assert(st > 0.5);
        if (primaryVarsMeaningWater() == WaterMeaning::Sw) {
            (*this)[Indices::waterSwitchIdx] = sw;
        }
        if (primaryVarsMeaningGas() == GasMeaning::Sg) {
            (*this)[Indices::compositionSwitchIdx] = sg;
        }
        if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
            (*this) [Indices::solventSaturationIdx] = ssol;
        }

        return (st != 1);
    }

    BlackOilPrimaryVariables& operator=(const BlackOilPrimaryVariables& other) = default;

    using ParentType::operator=; 

     OPM_HOST_DEVICE void checkDefined() const
    {
#ifndef NDEBUG
        // check the "real" primary variables
        for (unsigned i = 0; i < this->size(); ++i) {
            Valgrind::CheckDefined((*this)[i]);
        }

        // check the "pseudo" primary variables
        Valgrind::CheckDefined(primaryVarsMeaningWater_);
        Valgrind::CheckDefined(primaryVarsMeaningGas_);
        Valgrind::CheckDefined(primaryVarsMeaningPressure_);
        Valgrind::CheckDefined(primaryVarsMeaningBrine_);
        Valgrind::CheckDefined(primaryVarsMeaningSolvent_);

        Valgrind::CheckDefined(pvtRegionIdx_);
#endif // NDEBUG
    }

    template<class Serializer>
    void serializeOp(Serializer& serializer)
    {
        using FV = Dune::FieldVector<Scalar, getPropValue<TypeTag, Properties::NumEq>()>;
        serializer(static_cast<FV&>(*this));
        serializer(primaryVarsMeaningWater_);
        serializer(primaryVarsMeaningPressure_);
        serializer(primaryVarsMeaningGas_);
        serializer(primaryVarsMeaningBrine_);
        serializer(primaryVarsMeaningSolvent_);
        serializer(pvtRegionIdx_);
    }

    OPM_HOST_DEVICE bool operator==(const BlackOilPrimaryVariables& rhs) const
    {
        return
               static_cast<const FvBasePrimaryVariables<TypeTag>&>(*this) == rhs
            && this->primaryVarsMeaningWater_ == rhs.primaryVarsMeaningWater_
            && this->primaryVarsMeaningPressure_ == rhs.primaryVarsMeaningPressure_
            && this->primaryVarsMeaningGas_ == rhs.primaryVarsMeaningGas_
            && this->primaryVarsMeaningBrine_ == rhs.primaryVarsMeaningBrine_
            && this->primaryVarsMeaningSolvent_ == rhs.primaryVarsMeaningSolvent_
            && this->pvtRegionIdx_ == rhs.pvtRegionIdx_;
    }

private:
    OPM_HOST_DEVICE Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }

    OPM_HOST_DEVICE const Implementation& asImp_() const
    { return *static_cast<const Implementation*>(this); }

    OPM_HOST_DEVICE Scalar solventSaturation_() const
    {
        if constexpr (enableSolvent) {
            if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
                return (*this)[Indices::solventSaturationIdx];
            }
        }
        return 0.0;
    }

    OPM_HOST_DEVICE Scalar zFraction_() const
    {
        if constexpr (enableExtbo) {
            return (*this)[Indices::zFractionIdx];
        }
        else {
            return 0.0;
        }
    }

    OPM_HOST_DEVICE Scalar saltConcentration_() const
    {
        if constexpr (enableBrine) {
            return (*this)[Indices::saltConcentrationIdx];
        }
        else {
            return 0.0;
        }
    }

    Scalar biofilmVolumeFraction_() const
    {
        if constexpr (enableBioeffects)
            return (*this)[Indices::biofilmVolumeFractionIdx];
        else
            return 0.0;
    }

    OPM_HOST_DEVICE Scalar temperature_(const Problem& problem, [[maybe_unused]] unsigned globalDofIdx) const
    {
        if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
            return (*this)[Indices::temperatureIdx];
        }
        else if (energyModuleType == EnergyModules::NoTemperature) {
            return FluidSystem::reservoirTemperature();
        } else {
            return problem.temperature(globalDofIdx, /*timeIdx*/ 0);
        }
    }

    template <class Container>
    OPM_HOST_DEVICE void computeCapillaryPressures_(Container& result,
                                                    Scalar so,
                                                    Scalar sg,
                                                    Scalar sw,
                                                    const MaterialLawParams& matParams) const
    {
        using SatOnlyFluidState = SimpleModularFluidState<Scalar,
                                                          numPhases,
                                                          numComponents,
                                                          FluidSystem,
                                                          /*storePressure=*/false,
                                                          /*storeTemperature=*/false,
                                                          /*storeComposition=*/false,
                                                          /*storeFugacity=*/false,
                                                          /*storeSaturation=*/true,
                                                          /*storeDensity=*/false,
                                                          /*storeViscosity=*/false,
                                                          /*storeEnthalpy=*/false>;
        SatOnlyFluidState fluidState;
        fluidState.setSaturation(waterPhaseIdx, sw);
        fluidState.setSaturation(oilPhaseIdx, so);
        fluidState.setSaturation(gasPhaseIdx, sg);

        MaterialLaw::capillaryPressures(result, matParams, fluidState);
    }

    OPM_HOST_DEVICE Scalar pressure_() const
    {
        return (*this)[Indices::pressureSwitchIdx] * this->getPressureScale();
    }

    OPM_HOST_DEVICE constexpr Scalar getPressureScale() const
    {
        // In device code, we do not have access to static members, so we return 1.0
        // In most cases, the pressure scale is set to 1.0 anyway, so this is acceptable.
        #if OPM_IS_INSIDE_DEVICE_FUNCTION
        return Scalar(1.0);
        #else
        return this->pressureScale_;
        #endif
    }

    OPM_HOST_DEVICE void setScaledPressure_(Scalar pressure)
    { (*this)[Indices::pressureSwitchIdx] = pressure / (this->getPressureScale()); }

    // NOTE: When adding new member variables, be sure to update the template copy constructor,
    WaterMeaning primaryVarsMeaningWater_{WaterMeaning::Disabled};
    PressureMeaning primaryVarsMeaningPressure_{PressureMeaning::Po};
    GasMeaning primaryVarsMeaningGas_{GasMeaning::Disabled};
    BrineMeaning primaryVarsMeaningBrine_{BrineMeaning::Disabled};
    SolventMeaning primaryVarsMeaningSolvent_{SolventMeaning::Disabled};
    unsigned short pvtRegionIdx_{};
    Scalar pcFactor_{};
    inline static Scalar pressureScale_ = 1.0;
};

} // namespace Opm

#endif