blackoildiffusionmodule.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 OPM_BLACKOIL_DIFFUSION_MODULE_HH
#define OPM_BLACKOIL_DIFFUSION_MODULE_HH

#include <dune/common/fvector.hh>

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

#include <opm/input/eclipse/EclipseState/EclipseState.hpp>

#include <opm/material/common/Valgrind.hpp>

#include <opm/models/blackoil/blackoilbioeffectsmodules.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>

#include <array>
#include <stdexcept>

namespace Opm {

template <class TypeTag, bool enableDiffusion>
class BlackOilDiffusionModule;

template <class TypeTag, bool enableDiffusion>
class BlackOilDiffusionExtensiveQuantities;

template <class TypeTag>
class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/false>
{
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;

public:
    static void initFromState(const EclipseState&)
    {}

    static void registerParameters()
    {}

    template <class Context>
    OPM_HOST_DEVICE static void addDiffusiveFlux(RateVector&,
                                                 const Context&,
                                                 unsigned,
                                                 unsigned)
    {}
};

template <class TypeTag>
class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/true>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using BioeffectsModule = BlackOilBioeffectsModule<TypeTag>;
    using BioeffectsParams = BlackOilBioeffectsParams<TypeTag>;

    enum { numPhases = FluidSystem::numPhases };
    enum { conti0EqIdx = Indices::conti0EqIdx };

    enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
    enum { enableMICP = Indices::enableMICP };

    static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
    static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
    static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;

    using Toolbox = MathToolbox<Evaluation>;

public:
    using ExtensiveQuantities = BlackOilDiffusionExtensiveQuantities<TypeTag,true>;

    static void initFromState(const EclipseState& eclState)
    {
        // TODO: support use_move_fraction on GPUs
        #if OPM_IS_INSIDE_HOST_FUNCTION
        use_mole_fraction_ = eclState.getTableManager().diffMoleFraction();
        #else
        OPM_THROW(std::runtime_error, "Diffusion module: Calling initFromState not expected on GPU.");
        #endif
    }

    static void registerParameters()
    {}

    template <class Context>
    OPM_HOST_DEVICE static void addDiffusiveFlux(RateVector& flux, const Context& context,
                                                 unsigned spaceIdx, unsigned timeIdx)
    {
        // Only work if diffusion is enabled run-time by DIFFUSE in the deck
        if (!FluidSystem::enableDiffusion()) {
            return;
        }

        const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
        const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
        const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
        const auto& diffusivity = extQuants.diffusivity();
        const auto& effectiveDiffusionCoefficient = extQuants.effectiveDiffusionCoefficient();
        addDiffusiveFlux(flux, inIq, exIq, diffusivity, effectiveDiffusionCoefficient);
        if constexpr(enableBioeffects) {
            const auto& effectiveBioDiffCoefficient = extQuants.effectiveBioDiffCoefficient();
            addBioDiffFlux(flux, inIq, exIq, diffusivity, effectiveBioDiffCoefficient);
        }
    }

    template<class IntensiveQuantities,class EvaluationArray, class RateVectorT>
    OPM_HOST_DEVICE static void addDiffusiveFlux(RateVectorT& flux,
                                                 const IntensiveQuantities& inIq,
                                                 const IntensiveQuantities& exIq,
                                                 const Evaluation& diffusivity,
                                                 const EvaluationArray& effectiveDiffusionCoefficient)
    {
        const FluidSystem& fsys = inIq.getFluidSystem();

        const auto& inFs = inIq.fluidState();
        const auto& exFs = exIq.fluidState();
        unsigned pvtRegionIndex = inFs.pvtRegionIndex();
        Evaluation diffR = 0.0;

        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!fsys.phaseIsActive(phaseIdx)) {
                continue;
            }

            // no diffusion in water for blackoil models
            if (!fsys.enableDissolvedGasInWater() && fsys.waterPhaseIdx == phaseIdx) {
                continue;
            }

            // no diffusion in gas phase in water + gas system.
            if (fsys.gasPhaseIdx == phaseIdx && !fsys.phaseIsActive(fsys.oilPhaseIdx)) {
                continue;
            }

            // arithmetic mean of the phase's b factor weighed by saturation
            Evaluation bSAvg = inFs.saturation(phaseIdx) * inFs.invB(phaseIdx);
            bSAvg += Toolbox::value(exFs.saturation(phaseIdx)) * Toolbox::value(exFs.invB(phaseIdx));
            bSAvg /= 2;

            // phase not present, skip
            if (bSAvg < 1.0e-6) {
                continue;
            }
            Evaluation convFactor = 1.0;
            if (fsys.enableDissolvedGas() &&
                fsys.phaseIsActive(fsys.gasPhaseIdx) &&
                phaseIdx == fsys.oilPhaseIdx)
            {
                const Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2;
                convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex, fsys) + rsAvg);
                diffR = inFs.Rs() - Toolbox::value(exFs.Rs());
            }
            if (fsys.enableVaporizedOil() &&
                fsys.phaseIsActive(fsys.oilPhaseIdx) &&
                phaseIdx == fsys.gasPhaseIdx)
            {
                const Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2;
                convFactor = toFractionGasOil(pvtRegionIndex, fsys) /
                             (1.0 + rvAvg * toFractionGasOil(pvtRegionIndex, fsys));
                diffR = inFs.Rv() - Toolbox::value(exFs.Rv());
            }
            if (fsys.enableDissolvedGasInWater() && phaseIdx == fsys.waterPhaseIdx) {
                const Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2;
                convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex, fsys) + rsAvg);
                diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw());
            }
            if (fsys.enableVaporizedWater() && phaseIdx == fsys.gasPhaseIdx) {
                const Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2;
                convFactor = toFractionGasWater(pvtRegionIndex, fsys) /
                             (1.0 + rvAvg * toFractionGasWater(pvtRegionIndex, fsys));
                diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw());
            }

            // mass flux of solvent component (oil in oil or gas in gas)
            const unsigned solventCompIdx = fsys.solventComponentIndex(phaseIdx);
            const unsigned activeSolventCompIdx = fsys.canonicalToActiveCompIdx(solventCompIdx);
            flux[conti0EqIdx + activeSolventCompIdx] +=
                    -bSAvg *
                    convFactor *
                    diffR *
                    diffusivity *
                    effectiveDiffusionCoefficient[phaseIdx][solventCompIdx];

            // mass flux of solute component (gas in oil or oil in gas)
            const unsigned soluteCompIdx = fsys.soluteComponentIndex(phaseIdx);
            const unsigned activeSoluteCompIdx = fsys.canonicalToActiveCompIdx(soluteCompIdx);
            flux[conti0EqIdx + activeSoluteCompIdx] +=
                    bSAvg *
                    diffR *
                    convFactor *
                    diffusivity *
                    effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx];
        }
    }

    template<class IntensiveQuantities,class EvaluationArray>
    OPM_HOST_DEVICE static void addBioDiffFlux(RateVector& flux,
                                               const IntensiveQuantities& inIq,
                                               const IntensiveQuantities& exIq,
                                               const Evaluation& diffusivity,
                                               const EvaluationArray& effectiveBioDiffCoefficient)
    {
        const auto& inFs = inIq.fluidState();
        const auto& exFs = exIq.fluidState();
        Evaluation diffR = 0.0;

        // The diffusion coefficients are given for mass concentrations
        Evaluation bAvg = (inFs.saturation(waterPhaseIdx) * inFs.invB(waterPhaseIdx) +
            Toolbox::value(exFs.saturation(waterPhaseIdx)) * Toolbox::value(exFs.invB(waterPhaseIdx))) / 2;
        diffR = inIq.microbialConcentration() - Toolbox::value(exIq.microbialConcentration());
        flux[contiMicrobialEqIdx] +=
            bAvg *
            diffR *
            diffusivity *
            effectiveBioDiffCoefficient[BioeffectsParams::micrDiffIdx];
        if constexpr(enableMICP) {
            diffR = inIq.oxygenConcentration() - Toolbox::value(exIq.oxygenConcentration());
            flux[contiOxygenEqIdx] +=
                bAvg *
                diffR *
                diffusivity *
                effectiveBioDiffCoefficient[BioeffectsParams::oxygDiffIdx];
            diffR = inIq.ureaConcentration() - Toolbox::value(exIq.ureaConcentration());
            flux[contiUreaEqIdx] +=
                bAvg *
                diffR *
                diffusivity *
                effectiveBioDiffCoefficient[BioeffectsParams::ureaDiffIdx];
        }
    }

private:
    template<class FluidSystemT>
    OPM_HOST_DEVICE static Scalar toFractionGasOil (unsigned regionIdx, const FluidSystemT& fsys)
    {
        const Scalar mMOil = use_mole_fraction_ ? fsys.molarMass(fsys.oilCompIdx, regionIdx) : 1;
        const Scalar rhoO = fsys.referenceDensity(fsys.oilPhaseIdx, regionIdx);
        const Scalar mMGas = use_mole_fraction_ ? fsys.molarMass(fsys.gasCompIdx, regionIdx) : 1;
        const Scalar rhoG = fsys.referenceDensity(fsys.gasPhaseIdx, regionIdx);
        return rhoO * mMGas / (rhoG * mMOil);
    }

    template<class FluidSystemT>
    OPM_HOST_DEVICE static Scalar toFractionGasWater (unsigned regionIdx, const FluidSystemT& fsys)
    {
        const Scalar mMWater = use_mole_fraction_ ? fsys.molarMass(fsys.waterCompIdx, regionIdx) : 1;
        const Scalar rhoW = fsys.referenceDensity(fsys.waterPhaseIdx, regionIdx);
        const Scalar mMGas = use_mole_fraction_ ? fsys.molarMass(fsys.gasCompIdx, regionIdx) : 1;
        const Scalar rhoG = fsys.referenceDensity(fsys.gasPhaseIdx, regionIdx);
        return rhoW * mMGas / (rhoG * mMWater);
    }

    #if OPM_IS_INSIDE_DEVICE_FUNCTION
    constexpr static bool use_mole_fraction_ = false; // currently only false i supported on GPU.
    #else
    static bool use_mole_fraction_;
    #endif
};

#if OPM_IS_INSIDE_HOST_FUNCTION
template <typename TypeTag>
bool
BlackOilDiffusionModule<TypeTag, true>::use_mole_fraction_;
#endif

template <class TypeTag, bool enableDiffusion>
class BlackOilDiffusionIntensiveQuantities;

template <class TypeTag>
class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;

public:
    OPM_HOST_DEVICE Scalar tortuosity(unsigned) const
    {
        throw std::logic_error("Method tortuosity() does not make sense "
                               "if diffusion is disabled");
    }

    OPM_HOST_DEVICE Scalar diffusionCoefficient(unsigned, unsigned) const
    {
        throw std::logic_error("Method diffusionCoefficient() does not "
                               "make sense if diffusion is disabled");
    }

    OPM_HOST_DEVICE Scalar effectiveDiffusionCoefficient(unsigned, unsigned) const
    {
        throw std::logic_error("Method effectiveDiffusionCoefficient() "
                               "does not make sense if diffusion is disabled");
    }

protected:
    template <class FluidState>
    void update_(FluidState&,
                 const unsigned,
                 const ElementContext&,
                 unsigned,
                 unsigned)
    {}
};

template <class TypeTag>
class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using BioeffectsModule = BlackOilBioeffectsModule<TypeTag>;

    enum { numPhases = FluidSystem::numPhases };
    enum { numComponents = FluidSystem::numComponents };
    enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
    enum { enableMICP = Indices::enableMICP };
    enum { numBioInWat = Indices::numBioInWat };

    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;

public:
    BlackOilDiffusionIntensiveQuantities() = default;
    BlackOilDiffusionIntensiveQuantities(BlackOilDiffusionIntensiveQuantities&&) noexcept = default;
    BlackOilDiffusionIntensiveQuantities(const BlackOilDiffusionIntensiveQuantities&) = default;

    BlackOilDiffusionIntensiveQuantities& operator=(BlackOilDiffusionIntensiveQuantities&&) noexcept = default;

    BlackOilDiffusionIntensiveQuantities(const std::array<Evaluation, numPhases>& tortuosity,
                                         const std::array<std::array<Evaluation, numComponents>, numPhases>& diffusionCoefficient)
        : tortuosity_(tortuosity)
        , diffusionCoefficient_(diffusionCoefficient)
    {}

    BlackOilDiffusionIntensiveQuantities&
    operator=(const BlackOilDiffusionIntensiveQuantities& rhs)
    {
        if (this == &rhs) {
            return *this;
        }

        if (FluidSystem::enableDiffusion()) {
            tortuosity_ = rhs.tortuosity_;
            diffusionCoefficient_ = rhs.diffusionCoefficient_;
        }
        return *this;
    }

    OPM_HOST_DEVICE Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
    { return diffusionCoefficient_[phaseIdx][compIdx]; }

    OPM_HOST_DEVICE Evaluation tortuosity(unsigned phaseIdx) const
    { return tortuosity_[phaseIdx]; }

    OPM_HOST_DEVICE Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
    {
        // For the blackoil model tortuosity is disabled.
        // TODO add a run-time parameter to enable tortuosity
        static constexpr bool enableTortuosity = false;
        if constexpr (enableTortuosity) {
            return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx];
        } else {
            return diffusionCoefficient_[phaseIdx][compIdx];
        }
    }

    OPM_HOST_DEVICE Evaluation bioDiffCoefficient(unsigned compIdx) const
    { return bioDiffCoefficient_[compIdx]; }

    OPM_HOST_DEVICE Evaluation effectiveBioDiffCoefficient(unsigned compIdx) const
    {
        // TODO add a run-time parameter to enable tortuosity
        static bool enableTortuosity = false;
        if (enableTortuosity)
            return tortuosity_[waterPhaseIdx] * bioDiffCoefficient_[compIdx];

        return bioDiffCoefficient_[compIdx];
    }

protected:
    template <class FluidState>
    void update_(FluidState& fluidState,
                 const unsigned regionIdx,
                 const ElementContext& elemCtx,
                 unsigned dofIdx,
                 unsigned timeIdx)
    {
        // Only work if diffusion is enabled run-time by DIFFUSE in the deck
        if (!FluidSystem::enableDiffusion()) {
            return;
        }

        const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
        update_(fluidState, regionIdx, intQuants);
    }

    template<class FluidState>
    void update_(FluidState& fluidState,
                 const unsigned regionIdx,
                 const IntensiveQuantities& intQuants)
    {
        using Toolbox = MathToolbox<Evaluation>;

        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!FluidSystem::phaseIsActive(phaseIdx)) {
                continue;
            }

            // no diffusion in water for blackoil models
            if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
                continue;
            }

            // Based on Millington, R. J., & Quirk, J. P. (1961).
            // \Note: it is possible to use NumericalConstants later
            //  constexpr auto& numconst = GetPropValue<TypeTag, Properties::NumericalConstants>;
            constexpr double myeps = 0.0001; //numconst.blackoildiffusionmoduleeps;
            const Evaluation& base =
                Toolbox::max(myeps, //0.0001,
                             intQuants.porosity() *
                             intQuants.fluidState().saturation(phaseIdx));
            tortuosity_[phaseIdx] =
                1.0 / (intQuants.porosity() * intQuants.porosity()) *
                Toolbox::pow(base, 10.0 / 3.0);

            using PCache = typename FluidSystem::template ParameterCache<Scalar>;
            PCache pcache(regionIdx);
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                diffusionCoefficient_[phaseIdx][compIdx] =
                    FluidSystem::diffusionCoefficient(fluidState,
                                                      pcache,
                                                      phaseIdx,
                                                      compIdx);
            }
        }

        if constexpr(enableBioeffects) {
            unsigned pvtRegionIndex = intQuants.pvtRegionIndex();
            for (unsigned compIdx = 0; compIdx < numBioInWat; ++compIdx) {
                bioDiffCoefficient_[compIdx] =
                    BioeffectsModule::bioDiffCoefficient(pvtRegionIndex, compIdx);
            }
        }
    }

private:
    std::array<Evaluation, numPhases> tortuosity_{};
    std::array<std::array<Evaluation, numComponents>, numPhases> diffusionCoefficient_{};
    std::array<Evaluation, numBioInWat> bioDiffCoefficient_{};
};

template <class TypeTag, bool enableDiffusion>
class BlackOilDiffusionExtensiveQuantities;

template <class TypeTag>
class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/false>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;

protected:
    void update_(const ElementContext&,
                 unsigned,
                 unsigned)
    {}

    template <class Context, class FluidState>
    void updateBoundary_(const Context&,
                         unsigned,
                         unsigned,
                         const FluidState&)
    {}

public:
    OPM_HOST_DEVICE Scalar diffusivity() const
    {
        throw std::logic_error("The method diffusivity() does not "
                               "make sense if diffusion is disabled.");
    }

    OPM_HOST_DEVICE const Evaluation& effectiveDiffusionCoefficient(unsigned,
                                                                    unsigned) const
    {
        throw std::logic_error("The method effectiveDiffusionCoefficient() "
                               "does not make sense if diffusion is disabled.");
    }
};

template <class TypeTag>
class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;

    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
    enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
    enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
    enum { enableMICP = Indices::enableMICP };
    enum { numBioInWat = Indices::numBioInWat };

public:
    using EvaluationArray = std::array<std::array<Evaluation, numComponents>, numPhases>;

protected:
    void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
    {
        // Only work if diffusion is enabled run-time by DIFFUSE in the deck
        if (!FluidSystem::enableDiffusion()) {
            return;
        }

        const auto& stencil = elemCtx.stencil(timeIdx);
        const auto& face = stencil.interiorFace(faceIdx);
        const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
        const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
        const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);

        const Scalar diff = elemCtx.problem().diffusivity(elemCtx, face.interiorIndex(), face.exteriorIndex());
        const Scalar faceArea = face.area();
        diffusivity_ = diff / faceArea;
        update(effectiveDiffusionCoefficient_, intQuantsInside, intQuantsOutside);
        Valgrind::CheckDefined(diffusivity_);
        if constexpr(enableBioeffects) {
            updateBio(effectiveBioDiffCoefficient_, intQuantsInside, intQuantsOutside);
        }
    }

public:
    static void update(EvaluationArray& effectiveDiffusionCoefficient,
                       const IntensiveQuantities& intQuantsInside,
                       const IntensiveQuantities& intQuantsOutside)
    {
        const FluidSystem& fsys = intQuantsInside.getFluidSystem();

        // opm-models expects per area flux
        // use the arithmetic average for the effective
        // diffusion coefficients.
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!fsys.phaseIsActive(phaseIdx)) {
                continue;
            }
            // no diffusion in water for blackoil models
            if (!fsys.enableDissolvedGasInWater() && fsys.waterPhaseIdx == phaseIdx) {
                continue;
            }
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                effectiveDiffusionCoefficient[phaseIdx][compIdx] =
                    0.5 * (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) +
                           intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx));
                Valgrind::CheckDefined(effectiveDiffusionCoefficient[phaseIdx][compIdx]);
            }
        }
    }

    static void updateBio(std::array<Evaluation, numBioInWat>& effectiveBioDiffCoefficient,
                          const IntensiveQuantities& intQuantsInside,
                          const IntensiveQuantities& intQuantsOutside) {

        for (unsigned compIdx = 0; compIdx < numBioInWat; ++compIdx) {
            effectiveBioDiffCoefficient[compIdx] =
                0.5 * (intQuantsInside.effectiveBioDiffCoefficient(compIdx) +
                       intQuantsOutside.effectiveBioDiffCoefficient(compIdx));
            Valgrind::CheckDefined(effectiveBioDiffCoefficient[compIdx]);
        }
    }

protected:
    template <class Context, class FluidState>
    void updateBoundary_(const Context&,
                         unsigned,
                         unsigned,
                         const FluidState&)
    {
        throw std::runtime_error("Not implemented: Diffusion across boundary not implemented for blackoil");
    }

public:
    OPM_HOST_DEVICE Scalar diffusivity() const
    { return diffusivity_; }

    OPM_HOST_DEVICE const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
    { return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }

    OPM_HOST_DEVICE const auto& effectiveDiffusionCoefficient() const
    { return effectiveDiffusionCoefficient_; }

    OPM_HOST_DEVICE const Evaluation& effectiveBioDiffCoefficient(unsigned compIdx) const
    { return effectiveBioDiffCoefficient_[compIdx]; }

    OPM_HOST_DEVICE const auto& effectiveBioDiffCoefficient() const{
        return effectiveBioDiffCoefficient_;
    }

private:
    Scalar diffusivity_;
    EvaluationArray effectiveDiffusionCoefficient_;
    std::array<Evaluation, numBioInWat> effectiveBioDiffCoefficient_;
};

} // namespace Opm

#endif