porousmediumflow/mpnc/volumevariables.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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
 
#ifndef DUMUX_MPNC_VOLUME_VARIABLES_HH
#define DUMUX_MPNC_VOLUME_VARIABLES_HH
 
#include <dumux/porousmediumflow/volumevariables.hh>
#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
 
#include <dumux/material/constraintsolvers/compositionfromfugacities.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
 
#include "pressureformulation.hh"
 
namespace Dumux {
 
// forward declaration
template <class Traits, bool enableChemicalNonEquilibrium>
class MPNCVolumeVariablesImplementation;
template <class Traits>
using MPNCVolumeVariables =  MPNCVolumeVariablesImplementation<Traits, Traits::ModelTraits::enableChemicalNonEquilibrium()>;
 
template <class Traits>
class MPNCVolumeVariablesImplementation<Traits, false>
: public PorousMediumFlowVolumeVariables<Traits>
, public EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> >
{
    using ParentType = PorousMediumFlowVolumeVariables<Traits>;
    using EnergyVolVars = EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> >;
 
    using Scalar = typename Traits::PrimaryVariables::value_type;
    using PermeabilityType = typename Traits::PermeabilityType;
 
    using ModelTraits = typename Traits::ModelTraits;
    static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
 
    static constexpr bool enableThermalNonEquilibrium = ModelTraits::enableThermalNonEquilibrium();
    static constexpr bool enableChemicalNonEquilibrium = ModelTraits::enableChemicalNonEquilibrium();
    static constexpr bool enableDiffusion = ModelTraits::enableMolecularDiffusion();
 
    using ComponentVector = Dune::FieldVector<Scalar, ModelTraits::numFluidComponents()>;
    using CompositionFromFugacities = Dumux::CompositionFromFugacities<Scalar, typename Traits::FluidSystem>;
    using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
    using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
 
public:
    using Indices = typename Traits::ModelTraits::Indices;
    using FluidSystem = typename Traits::FluidSystem;
    using FluidState = typename Traits::FluidState;
    using SolidState = typename Traits::SolidState;
    using SolidSystem = typename Traits::SolidSystem;

    static constexpr int numFluidPhases() { return ModelTraits::numFluidPhases(); }
    static constexpr int numFluidComps = ParentType::numFluidComponents();

    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol& elemSol,
                const Problem& problem,
                const Element& element,
                const Scv& scv)
    {
        ParentType::update(elemSol, problem, element, scv);
 
        completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
 
        // calculate the remaining quantities
        const auto& spatialParams = problem.spatialParams();
        const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol);
        const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
        const auto relPerm = fluidMatrixInteraction.relativePermeabilities(fluidState_, wPhaseIdx);
        std::copy(relPerm.begin(), relPerm.end(), relativePermeability_.begin());
 
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState_);
 
        // porosity
        updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
 
        if constexpr (enableDiffusion)
        {
            auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
            { return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); };
 
            effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
        }
 
        EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
        EnergyVolVars::updateEffectiveThermalConductivity();
    }

 
    template<class ElemSol, class Problem, class Element, class Scv>
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
                            FluidState& fluidState,
                            SolidState& solidState)
    {

        // set the fluid phase temperatures
        EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);

        // set the phase saturations
        auto&& priVars = elemSol[scv.localDofIndex()];
        Scalar sumSat = 0;
        for (int phaseIdx = 0; phaseIdx < numFluidPhases() - 1; ++phaseIdx) {
            sumSat += priVars[Indices::s0Idx + phaseIdx];
            fluidState.setSaturation(phaseIdx, priVars[Indices::s0Idx + phaseIdx]);
        }
        fluidState.setSaturation(numFluidPhases() - 1, 1.0 - sumSat);

        // set the phase pressures
        // capillary pressure parameters
        const auto& spatialParams = problem.spatialParams();
        const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol);
        fluidState.setWettingPhase(wPhaseIdx);
        // capillary pressures
 
        const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
        const auto capPress = fluidMatrixInteraction.capillaryPressures(fluidState, wPhaseIdx);
 
        // add to the pressure of the first fluid phase
        // depending on which pressure is stored in the primary variables
        if(pressureFormulation == MpNcPressureFormulation::mostWettingFirst){
            // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector.
            // For two phases this means that there is one pressure as primary variable: pw
            const Scalar pw = priVars[Indices::p0Idx];
            for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx)
                fluidState.setPressure(phaseIdx, pw - capPress[0] + capPress[phaseIdx]);
        }
        else if(pressureFormulation == MpNcPressureFormulation::leastWettingFirst){
            // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector.
            // For two phases this means that there is one pressure as primary variable: pn
            const Scalar pn = priVars[Indices::p0Idx];
            for (int phaseIdx = numFluidPhases()-1; phaseIdx >= 0; --phaseIdx)
                fluidState.setPressure(phaseIdx, pn - capPress[numFluidPhases()-1] + capPress[phaseIdx]);
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "MPNCVolumeVariables do not support the chosen pressure formulation");

        // set the fluid compositions
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState);
 
        ComponentVector fug;
        // retrieve component fugacities
        for (int compIdx = 0; compIdx < numFluidComps; ++compIdx)
            fug[compIdx] = priVars[Indices::fug0Idx + compIdx];
 
        // calculate phase compositions
        for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) {
            // initial guess
            for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) {
                // sanity check to make sure we have non-zero fugacity coefficients
                if (FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx) == 0.0)
                    DUNE_THROW(NumericalProblem, "MPNCVolumeVariables do not support fluidsystems with fugacity coefficients of 0");
                
                // set initial guess of the component's mole fraction
                fluidState.setMoleFraction(phaseIdx, compIdx, 1.0/numFluidComps);
            }
 
            // calculate the phase composition from the component
            // fugacities
            CompositionFromFugacities::guessInitial(fluidState, paramCache, phaseIdx, fug);
            CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fug);
        }
 
        // dynamic viscosities
        for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) {
            // viscosities
            Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
            fluidState.setViscosity(phaseIdx, mu);
 
            // compute and set the enthalpy
            Scalar h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
            fluidState.setEnthalpy(phaseIdx, h);
        }
    }

    const FluidState &fluidState() const
    { return fluidState_; }

    const SolidState &solidState() const
    { return solidState_; }

    Scalar saturation(int phaseIdx) const
    { return fluidState_.saturation(phaseIdx); }

    Scalar massFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.massFraction(phaseIdx, compIdx); }

    Scalar moleFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.moleFraction(phaseIdx, compIdx); }

    Scalar molarity(const int phaseIdx, int compIdx) const
    { return fluidState_.molarity(phaseIdx, compIdx); }

    Scalar molarDensity(const int phaseIdx) const
    { return fluidState_.molarDensity(phaseIdx);}

    Scalar pressure(const int phaseIdx) const
    { return fluidState_.pressure(phaseIdx); }

    Scalar density(const int phaseIdx) const
    { return fluidState_.density(phaseIdx); }

    Scalar temperature() const
    { return fluidState_.temperature(0/* phaseIdx*/); }
 
    Scalar temperature(const int phaseIdx) const
    { return fluidState_.temperature(phaseIdx); }

    Scalar enthalpy(const int phaseIdx) const
    { return fluidState_.enthalpy(phaseIdx); }

    Scalar internalEnergy(const int phaseIdx) const
    { return fluidState_.internalEnergy(phaseIdx); }

    Scalar fluidThermalConductivity(const int phaseIdx) const
    { return FluidSystem::thermalConductivity(fluidState_, phaseIdx); }

    Scalar fugacity(const int compIdx) const
    { return fluidState_.fugacity(compIdx); }

    Scalar averageMolarMass(const int phaseIdx) const
    { return fluidState_.averageMolarMass(phaseIdx); }

    Scalar mobility(const unsigned int phaseIdx) const
    {
        return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx);
    }

    Scalar viscosity(const unsigned int phaseIdx) const
    { return fluidState_.viscosity(phaseIdx); }

    Scalar relativePermeability(const unsigned int phaseIdx) const
    { return relativePermeability_[phaseIdx]; }

    Scalar porosity() const
    { return solidState_.porosity(); }

    const PermeabilityType& permeability() const
    { return permeability_; }

    bool isPhaseActive(const unsigned int phaseIdx) const
    {
        return
            phasePresentIneq(fluidState(), phaseIdx) -
            phaseNotPresentIneq(fluidState(), phaseIdx)
            >= 0;
    }

    Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
    {
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updatePhase(fluidState_, phaseIdx);
        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
    }

    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
    { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }

    Scalar phaseNcp(const unsigned int phaseIdx) const
    {
        Scalar aEval = phaseNotPresentIneq(fluidState(), phaseIdx);
        Scalar bEval = phasePresentIneq(fluidState(), phaseIdx);
        if (aEval > bEval)
            return phasePresentIneq(fluidState(), phaseIdx);
        return phaseNotPresentIneq(fluidState(), phaseIdx);
    };

    Scalar phasePresentIneq(const FluidState &fluidState,
                            const unsigned int phaseIdx) const
    { return fluidState.saturation(phaseIdx); }

    Scalar phaseNotPresentIneq(const FluidState &fluidState,
                               const unsigned int phaseIdx) const
    {
        // difference of sum of mole fractions in the phase from 100%
        Scalar a = 1;
        for (int compIdx = 0; compIdx < numFluidComps; ++compIdx)
            a -= fluidState.moleFraction(phaseIdx, compIdx);
        return a;
    }
 
protected:
    Scalar porosity_; 
    std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_; 
    PermeabilityType permeability_;

    FluidState fluidState_;
    SolidState solidState_;
 
    // Effective diffusion coefficients for the phases
    DiffusionCoefficients effectiveDiffCoeff_;
};
 
template <class Traits>
class MPNCVolumeVariablesImplementation<Traits, true>
    : public PorousMediumFlowVolumeVariables<Traits>
    , public EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> >
{
    using ParentType = PorousMediumFlowVolumeVariables< Traits>;
    using EnergyVolVars = EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> >;
    using Scalar = typename Traits::PrimaryVariables::value_type;
    using PermeabilityType = typename Traits::PermeabilityType;
 
    using ModelTraits = typename Traits::ModelTraits;
    static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
 
    static constexpr bool enableThermalNonEquilibrium = ModelTraits::enableThermalNonEquilibrium();
    static constexpr bool enableChemicalNonEquilibrium = ModelTraits::enableChemicalNonEquilibrium();
    static constexpr bool enableDiffusion = ModelTraits::enableMolecularDiffusion();
 
    using Indices = typename ModelTraits::Indices;
    using ComponentVector = Dune::FieldVector<Scalar,  ModelTraits::numFluidComponents()>;
    using CompositionFromFugacities = Dumux::CompositionFromFugacities<Scalar, typename Traits::FluidSystem>;
    using ParameterCache = typename Traits::FluidSystem::ParameterCache;
    using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
    using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
 
public:
    using FluidSystem = typename Traits::FluidSystem;
    using FluidState = typename Traits::FluidState;
    using SolidState = typename Traits::SolidState;
    using SolidSystem = typename Traits::SolidSystem;

    static constexpr int numFluidPhases() { return ModelTraits::numFluidPhases(); }
    static constexpr int numFluidComps = ParentType::numFluidComponents();
    using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>;

    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol& elemSol,
                const Problem& problem,
                const Element& element,
                const Scv& scv)
    {
        ParentType::update(elemSol, problem, element, scv);
 
        completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
 
        // calculate the remaining quantities
        const auto& spatialParams = problem.spatialParams();
        const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol);
        const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
        const auto relPerm = fluidMatrixInteraction.relativePermeabilities(fluidState_, wPhaseIdx);
        std::copy(relPerm.begin(), relPerm.end(), relativePermeability_.begin());
 
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState_);
 
        updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
 
        if constexpr (enableDiffusion)
        {
            auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
            { return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); };
 
            effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
        }
 
        EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
        permeability_ = spatialParams.permeability(element, scv, elemSol);
        EnergyVolVars::updateEffectiveThermalConductivity();
    }

    template<class ElemSol, class Problem, class Element, class Scv>
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
                            FluidState& fluidState,
                            SolidState& solidState)
    {
        // set the fluid phase temperatures
        EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
        // set the phase saturations
        auto&& priVars = elemSol[scv.localDofIndex()];
        Scalar sumSat = 0;
        for (int phaseIdx = 0; phaseIdx < numFluidPhases() - 1; ++phaseIdx) {
            sumSat += priVars[Indices::s0Idx + phaseIdx];
            fluidState.setSaturation(phaseIdx, priVars[Indices::s0Idx + phaseIdx]);
        }
        fluidState.setSaturation(numFluidPhases() - 1, 1.0 - sumSat);

        // set the phase pressures
        // capillary pressure parameters
        const auto& spatialParams = problem.spatialParams();
        const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol);
        fluidState.setWettingPhase(wPhaseIdx);
        const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
        const auto capPress = fluidMatrixInteraction.capillaryPressures(fluidState, wPhaseIdx);
 
        // add to the pressure of the first fluid phase
        // depending on which pressure is stored in the primary variables
        if(pressureFormulation == MpNcPressureFormulation::mostWettingFirst){
            // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector.
            // For two phases this means that there is one pressure as primary variable: pw
            const Scalar pw = priVars[Indices::p0Idx];
            for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx)
                fluidState.setPressure(phaseIdx, pw - capPress[0] + capPress[phaseIdx]);
        }
        else if(pressureFormulation == MpNcPressureFormulation::leastWettingFirst){
            // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector.
            // For two phases this means that there is one pressure as primary variable: pn
            const Scalar pn = priVars[Indices::p0Idx];
            for (int phaseIdx = numFluidPhases()-1; phaseIdx >= 0; --phaseIdx)
                fluidState.setPressure(phaseIdx, pn - capPress[numFluidPhases()-1] + capPress[phaseIdx]);
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "MPNCVolumeVariables do not support the chosen pressure formulation");
 

        // set the fluid compositions
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState);
 
        updateMoleFraction(fluidState, paramCache, priVars);
 
        // dynamic viscosities
        for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) {
            // viscosities
            Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
            fluidState.setViscosity(phaseIdx, mu);
 
            // compute and set the enthalpy
            Scalar h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
            fluidState.setEnthalpy(phaseIdx, h);
        }
   }

    void updateMoleFraction(FluidState & actualFluidState,
                            ParameterCache & paramCache,
                            const typename Traits::PrimaryVariables& priVars)
    {
        // setting the mole fractions of the fluid state
        for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx)
        {
                // set the component mole fractions
                for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) {
                    actualFluidState.setMoleFraction(phaseIdx,
                           compIdx,
                           priVars[Indices::moleFrac00Idx +
                                   phaseIdx*numFluidComps +
                                   compIdx]);
                }
        }
 
        // TODO: Can this be removed, are fugacity coefficients used in the constraint solver?
        for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx)
            for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) {
                const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState,
                                                                        paramCache,
                                                                        phaseIdx,
                                                                        compIdx);
                actualFluidState.setFugacityCoefficient(phaseIdx,
                                                      compIdx,
                                                      phi);
            }
 
        FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium
        equilFluidState.assign(actualFluidState) ;
        ConstraintSolver::solve(equilFluidState,
                                    paramCache) ;
 
        // Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction)
        for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx){
            for (int compIdx=0; compIdx< numFluidComps; ++ compIdx){
                xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx);
            }
        }
 
        // compute densities of all phases
        for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx){
            const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx);
            actualFluidState.setDensity(phaseIdx, rho);
            const Scalar rhoMolar = FluidSystem::molarDensity(actualFluidState, paramCache, phaseIdx);
            actualFluidState.setMolarDensity(phaseIdx, rhoMolar);
        }
 
    }

    const Scalar xEquil(const unsigned int phaseIdx, const unsigned int compIdx) const
    {
        return xEquil_[phaseIdx][compIdx] ;
    }

    const FluidState &fluidState() const
    { return fluidState_; }

    const SolidState &solidState() const
    { return solidState_; }

    Scalar saturation(int phaseIdx) const
    { return fluidState_.saturation(phaseIdx); }

    Scalar massFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.massFraction(phaseIdx, compIdx); }

    Scalar moleFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.moleFraction(phaseIdx, compIdx); }

    Scalar molarity(const int phaseIdx, int compIdx) const
    { return fluidState_.molarity(phaseIdx, compIdx); }

    Scalar molarDensity(const int phaseIdx) const
    { return fluidState_.molarDensity(phaseIdx);}

    Scalar pressure(const int phaseIdx) const
    { return fluidState_.pressure(phaseIdx); }

    Scalar density(const int phaseIdx) const
    { return fluidState_.density(phaseIdx); }

    Scalar temperature() const
    { return fluidState_.temperature(0/* phaseIdx*/); }

    Scalar enthalpy(const int phaseIdx) const
    { return fluidState_.enthalpy(phaseIdx); }

    Scalar internalEnergy(const int phaseIdx) const
    { return fluidState_.internalEnergy(phaseIdx); }

    Scalar fluidThermalConductivity(const int phaseIdx) const
    { return FluidSystem::thermalConductivity(fluidState_, phaseIdx); }

    Scalar fugacity(const int compIdx) const
    { return fluidState_.fugacity(compIdx); }

    Scalar averageMolarMass(const int phaseIdx) const
    { return fluidState_.averageMolarMass(phaseIdx); }

    Scalar mobility(const unsigned int phaseIdx) const
    {
        return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx);
    }

    Scalar viscosity(const unsigned int phaseIdx) const
    { return fluidState_.viscosity(phaseIdx); }

    Scalar relativePermeability(const unsigned int phaseIdx) const
    { return relativePermeability_[phaseIdx]; }

    Scalar porosity() const
    { return solidState_.porosity(); }

    const PermeabilityType& permeability() const
    { return permeability_; }

    bool isPhaseActive(const unsigned int phaseIdx) const
    {
        return
            phasePresentIneq(fluidState(), phaseIdx) -
            phaseNotPresentIneq(fluidState(), phaseIdx)
            >= 0;
    }

    Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
    {
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updatePhase(fluidState_, phaseIdx);
        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
    }

    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
    { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }

    Scalar phaseNcp(const unsigned int phaseIdx) const
    {
        Scalar aEval = phaseNotPresentIneq(fluidState(), phaseIdx);
        Scalar bEval = phasePresentIneq(fluidState(), phaseIdx);
        if (aEval > bEval)
            return phasePresentIneq(fluidState(), phaseIdx);
        return phaseNotPresentIneq(fluidState(), phaseIdx);
    };

    Scalar phasePresentIneq(const FluidState &fluidState,
                            const unsigned int phaseIdx) const
    { return fluidState.saturation(phaseIdx); }

    Scalar phaseNotPresentIneq(const FluidState &fluidState,
                               const unsigned int phaseIdx) const
    {
        // difference of sum of mole fractions in the phase from 100%
        Scalar a = 1;
        for (int compIdx = 0; compIdx < numFluidComps; ++compIdx)
            a -= fluidState.moleFraction(phaseIdx, compIdx);
        return a;
    }
 
protected:
    std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_; 
    PermeabilityType permeability_;
    std::array<std::array<Scalar, numFluidComps>, numFluidPhases()> xEquil_;

    FluidState fluidState_;
    SolidState solidState_;
 
    // Effective diffusion coefficients for the phases
    DiffusionCoefficients effectiveDiffCoeff_;
};
 
} // end namespace
 
#endif