blackoillocalresidualtpfa.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_LOCAL_TPFA_RESIDUAL_HH
#define EWOMS_BLACK_OIL_LOCAL_TPFA_RESIDUAL_HH

#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/input/eclipse/Schedule/BCProp.hpp>

#include <opm/material/common/ConditionalStorage.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>

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

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

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

#include <array>
#include <cassert>
#include <limits>
#include <stdexcept>
#include <string>

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

namespace Opm
{
template <class TypeTag>
class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLocalResidual>
{
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    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 GridView = GetPropType<TypeTag, Properties::GridView>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using FluidState = typename IntensiveQuantities::FluidState;

    enum { conti0EqIdx = Indices::conti0EqIdx };
    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
    enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };

    enum { dimWorld = GridView::dimensionworld };
    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };

    enum { gasCompIdx = FluidSystem::gasCompIdx };
    enum { oilCompIdx = FluidSystem::oilCompIdx };
    enum { waterCompIdx = FluidSystem::waterCompIdx };
    static constexpr unsigned compositionSwitchIdx = Indices::compositionSwitchIdx;

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

    static constexpr bool blackoilConserveSurfaceVolume
        = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();

    static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
    static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
    static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
    static constexpr EnergyModules energyModuleType =
        getPropValue<TypeTag, Properties::EnergyModuleType>();
    static constexpr bool enableFullyImplicitThermal =
        energyModuleType == EnergyModules::FullyImplicitThermal;
    static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
    static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
    static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
    static constexpr bool enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>();
    static constexpr bool enableConvectiveMixing
        = getPropValue<TypeTag, Properties::EnableConvectiveMixing>();
    static constexpr bool enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();
    static constexpr bool enableSaltPrecipitation
        = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
    static constexpr bool enableMICP = Indices::enableMICP;
    static constexpr bool runAssemblyOnGpu = getPropValue<TypeTag, Properties::RunAssemblyOnGpu>();

    using BioeffectsModule = BlackOilBioeffectsModule<TypeTag, enableBioeffects>;
    using BrineModule = BlackOilBrineModule<TypeTag, enableBrine>;
    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
    using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
    using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
    using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
    using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
    using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
    using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
    using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;

    using ModuleParams = BlackoilModuleParams<ConvectiveMixingModuleParam<Scalar>>;
    using Toolbox = MathToolbox<Evaluation>;

public:
    struct ResidualNBInfo {
        double trans;
        double faceArea;
        double thpresInToEx;
        double thpresExToIn;
        double dZg;
        FaceDir::DirEnum faceDir;
        double Vin;
        double Vex;
        ConditionalStorage<enableFullyImplicitThermal, double> inAlpha;
        ConditionalStorage<enableFullyImplicitThermal, double> outAlpha;
        ConditionalStorage<enableDiffusion, double> diffusivity;
        ConditionalStorage<enableDispersion, double> dispersivity;
    };

    template <class LhsEval>
    void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
                        const ElementContext& elemCtx,
                        unsigned dofIdx,
                        unsigned timeIdx) const
    {
        const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
        computeStorage<LhsEval>(storage, intQuants);
    }

    template <class LhsEval, class StorageType, class IntensiveQuantitiesType = IntensiveQuantities>
    OPM_HOST_DEVICE static void computeStorage(StorageType& storage,
                                               const IntensiveQuantitiesType& intQuants)
    {
        OPM_TIMEBLOCK_LOCAL(computeStorage, Subsystem::Assembly);
        // retrieve the intensive quantities for the SCV at the specified point in time
        const auto& fs = intQuants.fluidState();
        storage = 0.0;

        const FluidSystem& fsys = intQuants.getFluidSystem();

        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!fsys.phaseIsActive(phaseIdx)) {
                continue;
            }
            unsigned activeCompIdx
                = fsys.canonicalToActiveCompIdx(fsys.solventComponentIndex(phaseIdx));
            LhsEval surfaceVolume = Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
                * Toolbox::template decay<LhsEval>(fs.invB(phaseIdx))
                * Toolbox::template decay<LhsEval>(intQuants.porosity());

            storage[conti0EqIdx + activeCompIdx] += surfaceVolume;

            // account for dissolved gas
            if (phaseIdx == oilPhaseIdx && fsys.enableDissolvedGas()) {
                unsigned activeGasCompIdx = fsys.canonicalToActiveCompIdx(gasCompIdx);
                storage[conti0EqIdx + activeGasCompIdx]
                    += Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs())
                    * surfaceVolume;
            }

            // account for dissolved gas in water
            if (phaseIdx == waterPhaseIdx && fsys.enableDissolvedGasInWater()) {
                unsigned activeGasCompIdx = fsys.canonicalToActiveCompIdx(gasCompIdx);
                storage[conti0EqIdx + activeGasCompIdx]
                    += Toolbox::template decay<LhsEval>(intQuants.fluidState().Rsw())
                    * surfaceVolume;
            }

            // account for vaporized oil
            if (phaseIdx == gasPhaseIdx && fsys.enableVaporizedOil()) {
                unsigned activeOilCompIdx = fsys.canonicalToActiveCompIdx(oilCompIdx);
                storage[conti0EqIdx + activeOilCompIdx]
                    += Toolbox::template decay<LhsEval>(intQuants.fluidState().Rv())
                    * surfaceVolume;
            }

            // account for vaporized water
            if (phaseIdx == gasPhaseIdx && fsys.enableVaporizedWater()) {
                unsigned activeWaterCompIdx = fsys.canonicalToActiveCompIdx(waterCompIdx);
                storage[conti0EqIdx + activeWaterCompIdx]
                    += Toolbox::template decay<LhsEval>(intQuants.fluidState().Rvw())
                    * surfaceVolume;
            }
        }

        adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex(), fsys);

        // deal with solvents (if present)
        if constexpr (enableSolvent) {
            SolventModule::addStorage(storage, intQuants);
        }

        // deal with zFracton (if present)
        if constexpr (enableExtbo) {
            ExtboModule::addStorage(storage, intQuants);
        }

        // deal with polymer (if present)
        if constexpr (enablePolymer) {
            PolymerModule::addStorage(storage, intQuants);
        }

        // deal with energy (if present)
        if constexpr (enableFullyImplicitThermal) {
            EnergyModule::addStorage(storage, intQuants);
        }

        // deal with foam (if present)
        if constexpr (enableFoam) {
            FoamModule::addStorage(storage, intQuants);
        }

        // deal with salt (if present)
        if constexpr (enableBrine) {
            BrineModule::addStorage(storage, intQuants);
        }

        // deal with bioeffects (if present)
        if constexpr (enableBioeffects) {
            BioeffectsModule::addStorage(storage, intQuants);
        }
    }

    template <class ModuleParamsT,
              class RateVectorT,
              class IntensiveQuantitiesT,
              class ResidualNBInfoT>
    OPM_HOST_DEVICE static void computeFlux(RateVectorT& flux,
                                            RateVectorT& darcy,
                                            const unsigned globalIndexIn,
                                            const unsigned globalIndexEx,
                                            const IntensiveQuantitiesT& intQuantsIn,
                                            const IntensiveQuantitiesT& intQuantsEx,
                                            const ResidualNBInfoT& nbInfo,
                                            const ModuleParamsT& moduleParams)
    {
        OPM_TIMEBLOCK_LOCAL(computeFlux, Subsystem::Assembly);
        flux = 0.0;
        darcy = 0.0;

        calculateFluxes_(flux,
                         darcy,
                         intQuantsIn,
                         intQuantsEx,
                         globalIndexIn,
                         globalIndexEx,
                         nbInfo,
                         moduleParams);
    }

    // This function demonstrates compatibility with the ElementContext-based interface.
    // Actually using it will lead to double work since the element context already contains
    // fluxes through its stored ExtensiveQuantities.
    static void
    computeFlux(RateVector& flux, const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
    {
        OPM_TIMEBLOCK_LOCAL(computeFlux, Subsystem::Assembly);
        assert(timeIdx == 0);

        flux = 0.0;
        RateVector darcy = 0.0;
        // need for dary flux calculation
        const auto& problem = elemCtx.problem();
        const auto& stencil = elemCtx.stencil(timeIdx);
        const auto& scvf = stencil.interiorFace(scvfIdx);

        unsigned interiorDofIdx = scvf.interiorIndex();
        unsigned exteriorDofIdx = scvf.exteriorIndex();
        assert(interiorDofIdx != exteriorDofIdx);

        // unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
        // unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
        Scalar Vin = elemCtx.dofVolume(interiorDofIdx, /*timeIdx=*/0);
        Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, /*timeIdx=*/0);
        const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx);
        const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx);
        Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
        Scalar faceArea = scvf.area();
        const auto faceDir = faceDirFromDirId(scvf.dirId());
        Scalar thpresInToEx = problem.thresholdPressure(globalIndexIn, globalIndexEx);
        Scalar thpresExToIn = problem.thresholdPressure(globalIndexEx, globalIndexIn);

        // estimate the gravity correction: for performance reasons we use a simplified
        // approach for this flux module that assumes that gravity is constant and always
        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
        const Scalar g = problem.gravity()[dimWorld - 1];
        const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
        const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);

        // this is quite hacky because the dune grid interface does not provide a
        // cellCenterDepth() method (so we ask the problem to provide it). The "good"
        // solution would be to take the Z coordinate of the element centroids, but since
        // ECL seems to like to be inconsistent on that front, it needs to be done like
        // here...
        const Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
        const Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
        // the distances from the DOF's depths. (i.e., the additional depth of the
        // exterior DOF)
        const Scalar distZ = zIn - zEx;
        // for thermal harmonic mean of half trans
        const Scalar inAlpha = problem.thermalHalfTransmissibility(globalIndexIn, globalIndexEx);
        const Scalar outAlpha = problem.thermalHalfTransmissibility(globalIndexEx, globalIndexIn);
        const Scalar diffusivity = problem.diffusivity(globalIndexEx, globalIndexIn);
        const Scalar dispersivity = problem.dispersivity(globalIndexEx, globalIndexIn);

        const ResidualNBInfo res_nbinfo {trans,
                                         faceArea,
                                         thpresInToEx,
                                         thpresExToIn,
                                         distZ * g,
                                         faceDir,
                                         Vin,
                                         Vex,
                                         inAlpha,
                                         outAlpha,
                                         diffusivity,
                                         dispersivity};

        calculateFluxes_(flux,
                         darcy,
                         intQuantsIn,
                         intQuantsEx,
                         globalIndexIn,
                         globalIndexEx,
                         res_nbinfo,
                         problem.moduleParams());
    }

    template <class RateVectorT,
              class IntensiveQuantitiesT,
              class ResidualNBInfoT,
              class ModuleParamsT>
    OPM_HOST_DEVICE static void calculateFluxes_(RateVectorT& flux,
                                                 RateVectorT& darcy,
                                                 const IntensiveQuantitiesT& intQuantsIn,
                                                 const IntensiveQuantitiesT& intQuantsEx,
                                                 const unsigned& globalIndexIn,
                                                 const unsigned& globalIndexEx,
                                                 const ResidualNBInfoT& nbInfo,
                                                 const ModuleParamsT& moduleParams)
    {
        OPM_TIMEBLOCK_LOCAL(calculateFluxes, Subsystem::Assembly);
        const Scalar Vin = nbInfo.Vin;
        const Scalar Vex = nbInfo.Vex;
        const Scalar distZg = nbInfo.dZg;
        const Scalar thpresInToEx = nbInfo.thpresInToEx;
        const Scalar thpresExToIn = nbInfo.thpresExToIn;
        const Scalar trans = nbInfo.trans;
        const Scalar faceArea = nbInfo.faceArea;
        FaceDir::DirEnum facedir = nbInfo.faceDir;

        const FluidSystem& fsys = intQuantsIn.getFluidSystem();

        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!fsys.phaseIsActive(phaseIdx)) {
                continue;
            }
            // darcy flux calculation
            short dnIdx;
            //
            short upIdx;
            // fake intices should only be used to get upwind anc compatibility with old functions
            short interiorDofIdx = 0; // NB
            short exteriorDofIdx = 1; // NB
            Evaluation pressureDifference;
            ExtensiveQuantities::calculatePhasePressureDiff_(upIdx,
                                                             dnIdx,
                                                             pressureDifference,
                                                             intQuantsIn,
                                                             intQuantsEx,
                                                             phaseIdx, // input
                                                             interiorDofIdx, // input
                                                             exteriorDofIdx, // input
                                                             Vin,
                                                             Vex,
                                                             globalIndexIn,
                                                             globalIndexEx,
                                                             distZg,
                                                             thpresInToEx,
                                                             thpresExToIn,
                                                             moduleParams);

            const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
            unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
            // Use arithmetic average (more accurate with harmonic, but that requires recomputing
            // the transmissbility)
            Evaluation transMult = (intQuantsIn.rockCompTransMultiplier()
                                    + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))
                / 2;
            if constexpr (enableBioeffects || enableSaltPrecipitation) {
                transMult
                    *= (intQuantsIn.permFactor() + Toolbox::value(intQuantsEx.permFactor())) / 2;
            }

            Evaluation darcyFlux;
            if (globalUpIndex == globalIndexIn) {
                darcyFlux = pressureDifference * up.mobility(phaseIdx, facedir) * transMult
                    * (-trans / faceArea);
            } else {
                darcyFlux = pressureDifference
                    * (Toolbox::value(up.mobility(phaseIdx, facedir)) * transMult
                       * (-trans / faceArea));
            }

            unsigned activeCompIdx
                = fsys.canonicalToActiveCompIdx(fsys.solventComponentIndex(phaseIdx));
            // NB! For the FLORES fluxes without derivatives
            darcy[conti0EqIdx + activeCompIdx] = darcyFlux.value() * faceArea;

            unsigned pvtRegionIdx = up.pvtRegionIndex();
            // if (upIdx == globalFocusDofIdx){
            if (globalUpIndex == globalIndexIn) {
                const auto& invB = getInvB_<FluidSystem, FluidState, Evaluation>(
                    up.fluidState(), phaseIdx, pvtRegionIdx, fsys);
                const auto& surfaceVolumeFlux = invB * darcyFlux;

                evalPhaseFluxes_<Evaluation>(
                    flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
                if constexpr (enableFullyImplicitThermal) {
                    EnergyModule::template addPhaseEnthalpyFluxes_<Evaluation>(
                        flux, phaseIdx, darcyFlux, up.fluidState());
                }
                if constexpr (enableBioeffects) {
                    BioeffectsModule::template addBioeffectsFluxes_<Evaluation>(
                        flux, phaseIdx, darcyFlux, up);
                }
                if constexpr (enableBrine) {
                    BrineModule::template addBrineFluxes_<Evaluation, FluidState>(
                        flux, phaseIdx, darcyFlux, up.fluidState());
                }
            } else {
                const auto& invB = getInvB_<FluidSystem, FluidState, Scalar>(
                    up.fluidState(), phaseIdx, pvtRegionIdx, fsys);
                const auto& surfaceVolumeFlux = invB * darcyFlux;
                evalPhaseFluxes_<Scalar>(
                    flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
                if constexpr (enableFullyImplicitThermal) {
                    EnergyModule::template addPhaseEnthalpyFluxes_<Scalar>(
                        flux, phaseIdx, darcyFlux, up.fluidState());
                }
                if constexpr (enableBioeffects) {
                    BioeffectsModule::template addBioeffectsFluxes_<Scalar>(
                        flux, phaseIdx, darcyFlux, up);
                }
                if constexpr (enableBrine) {
                    BrineModule::template addBrineFluxes_<Scalar, FluidState>(
                        flux, phaseIdx, darcyFlux, up.fluidState());
                }
            }
        }

        // deal with solvents (if present)
        static_assert(
            !enableSolvent,
            "Relevant computeFlux() method must be implemented for this module before enabling.");
        // SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);

        // deal with zFracton (if present)
        static_assert(
            !enableExtbo,
            "Relevant computeFlux() method must be implemented for this module before enabling.");
        // ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);

        // deal with polymer (if present)
        static_assert(
            !enablePolymer,
            "Relevant computeFlux() method must be implemented for this module before enabling.");
        // PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);

        // deal with convective mixing
        if constexpr (enableConvectiveMixing) {
            ConvectiveMixingModule::addConvectiveMixingFlux(
                flux,
                intQuantsIn,
                intQuantsEx,
                globalIndexIn,
                globalIndexEx,
                nbInfo.dZg,
                nbInfo.trans,
                nbInfo.faceArea,
                moduleParams.convectiveMixingModuleParam);
        }

        // deal with energy (if present)
        if constexpr (enableFullyImplicitThermal) {
            const Scalar inAlpha = nbInfo.inAlpha;
            const Scalar outAlpha = nbInfo.outAlpha;
            Evaluation heatFlux;

            short interiorDofIdx = 0; // NB
            short exteriorDofIdx = 1; // NB

            EnergyModule::ExtensiveQuantities::updateEnergy(heatFlux,
                                                            interiorDofIdx, // focusDofIndex,
                                                            interiorDofIdx,
                                                            exteriorDofIdx,
                                                            intQuantsIn,
                                                            intQuantsEx,
                                                            intQuantsIn.fluidState(),
                                                            intQuantsEx.fluidState(),
                                                            inAlpha,
                                                            outAlpha,
                                                            faceArea);
            EnergyModule::addHeatFlux(flux, heatFlux);
        }
        // NB need to be tha last energy call since it does scaling
        // EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);

        // deal with foam (if present)
        static_assert(
            !enableFoam,
            "Relevant computeFlux() method must be implemented for this module before enabling.");
        // FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);

        // deal with diffusion (if present). opm-models expects per area flux (added in the
        // tmpdiffusivity).
        if constexpr (enableDiffusion) {
            typename DiffusionModule::ExtensiveQuantities::EvaluationArray
                effectiveDiffusionCoefficient;
            DiffusionModule::ExtensiveQuantities::update(
                effectiveDiffusionCoefficient, intQuantsIn, intQuantsEx);
            const Scalar diffusivity = nbInfo.diffusivity;
            const Scalar tmpdiffusivity = diffusivity / faceArea;
            DiffusionModule::addDiffusiveFlux(
                flux, intQuantsIn, intQuantsEx, tmpdiffusivity, effectiveDiffusionCoefficient);
        }

        // deal with dispersion (if present). opm-models expects per area flux (added in the
        // tmpdispersivity).
        if constexpr (enableDispersion) {
            typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg;
            DispersionModule::ExtensiveQuantities::update(
                normVelocityAvg, intQuantsIn, intQuantsEx);
            const Scalar dispersivity = nbInfo.dispersivity;
            const Scalar tmpdispersivity = dispersivity / faceArea;
            DispersionModule::addDispersiveFlux(
                flux, intQuantsIn, intQuantsEx, tmpdispersivity, normVelocityAvg);
        }

        // apply the scaling for the urea equation in MICP
        if constexpr (enableMICP) {
            BioeffectsModule::applyScaling(flux);
        }
    }

    template <class BoundaryConditionData>
    static void computeBoundaryFluxRate(RateVector& bdyFlux, const BoundaryConditionData& bdyInfo)
    {
        bdyFlux.setMassRate(bdyInfo.massRate, bdyInfo.pvtRegionIdx);
    }

    template <class BoundaryConditionData>
    static void computeBoundaryFluxFree(const Problem& problem,
                                        RateVector& bdyFlux,
                                        const BoundaryConditionData& bdyInfo,
                                        const IntensiveQuantities& insideIntQuants,
                                        unsigned globalSpaceIdx)
    {
        OPM_TIMEBLOCK_LOCAL(computeBoundaryFluxFree, Subsystem::Assembly);
        std::array<short, numPhases> upIdx;
        std::array<short, numPhases> dnIdx;
        std::array<Evaluation, numPhases> volumeFlux;
        std::array<Evaluation, numPhases> pressureDifference;

        ExtensiveQuantities::calculateBoundaryGradients_(problem,
                                                         globalSpaceIdx,
                                                         insideIntQuants,
                                                         bdyInfo.boundaryFaceIndex,
                                                         bdyInfo.faceArea,
                                                         bdyInfo.faceZCoord,
                                                         bdyInfo.exFluidState,
                                                         upIdx,
                                                         dnIdx,
                                                         volumeFlux,
                                                         pressureDifference);

        // advective fluxes of all components in all phases
        bdyFlux = 0.0;
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!FluidSystem::phaseIsActive(phaseIdx)) {
                continue;
            }
            const auto& pBoundary = bdyInfo.exFluidState.pressure(phaseIdx);
            const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx);
            const unsigned pvtRegionIdx = insideIntQuants.pvtRegionIndex();

            RateVector tmp(0.0);
            const auto& darcyFlux = volumeFlux[phaseIdx];
            // mass conservation
            if (pBoundary < pInside) {
                // outflux
                const auto& invB = getInvB_<FluidSystem, FluidState, Evaluation>(
                    insideIntQuants.fluidState(), phaseIdx, pvtRegionIdx);
                Evaluation surfaceVolumeFlux = invB * darcyFlux;
                evalPhaseFluxes_<Evaluation>(tmp,
                                             phaseIdx,
                                             insideIntQuants.pvtRegionIndex(),
                                             surfaceVolumeFlux,
                                             insideIntQuants.fluidState());
                if constexpr (enableFullyImplicitThermal) {
                    EnergyModule::template addPhaseEnthalpyFluxes_<Evaluation>(
                        tmp, phaseIdx, darcyFlux, insideIntQuants.fluidState());
                }
            } else if (pBoundary > pInside) {
                // influx
                using ScalarFluidState = decltype(bdyInfo.exFluidState);
                const auto& invB = getInvB_<FluidSystem, ScalarFluidState, Scalar>(
                    bdyInfo.exFluidState, phaseIdx, pvtRegionIdx);
                Evaluation surfaceVolumeFlux = invB * darcyFlux;
                evalPhaseFluxes_<Scalar>(tmp,
                                         phaseIdx,
                                         insideIntQuants.pvtRegionIndex(),
                                         surfaceVolumeFlux,
                                         bdyInfo.exFluidState);
                if constexpr (enableFullyImplicitThermal) {
                    EnergyModule::template addPhaseEnthalpyFluxes_<Scalar>(
                        tmp, phaseIdx, darcyFlux, bdyInfo.exFluidState);
                }
            }

            for (unsigned i = 0; i < tmp.size(); ++i) {
                bdyFlux[i] += tmp[i];
            }
        }

        // conductive heat flux from boundary
        if constexpr (enableFullyImplicitThermal) {
            Evaluation heatFlux;
            // avoid overload of functions with same number of elements in eclproblem
            Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(
                globalSpaceIdx, bdyInfo.boundaryFaceIndex);
            unsigned inIdx = 0; // dummy
            // always calculated with derivatives of this cell
            EnergyModule::ExtensiveQuantities::updateEnergyBoundary(heatFlux,
                                                                    insideIntQuants,
                                                                    /*focusDofIndex*/ inIdx,
                                                                    inIdx,
                                                                    alpha,
                                                                    bdyInfo.exFluidState);
            EnergyModule::addHeatFlux(bdyFlux, heatFlux);
        }

        static_assert(
            !enableSolvent,
            "Relevant treatment of boundary conditions must be implemented before enabling.");
        static_assert(
            !enablePolymer,
            "Relevant treatment of boundary conditions must be implemented before enabling.");

        const FluidSystem& fsys = insideIntQuants.getFluidSystem();

        // make sure that the right mass conservation quantities are used
        adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex(), fsys);

#ifndef NDEBUG
        for (unsigned i = 0; i < numEq; ++i) {
            Valgrind::CheckDefined(bdyFlux[i]);
        }
        Valgrind::CheckDefined(bdyFlux);
#endif
    }

    template <class ProblemLocal, class BoundaryConditionData, class RateVectorLocal>
    OPM_HOST_DEVICE static void computeBoundaryThermal(const ProblemLocal& problem,
                                                       RateVectorLocal& bdyFlux,
                                                       const BoundaryConditionData& bdyInfo,
                                                       const IntensiveQuantities& insideIntQuants,
                                                       [[maybe_unused]] unsigned globalSpaceIdx)
    {
        OPM_TIMEBLOCK_LOCAL(computeBoundaryThermal, Subsystem::Assembly);
        // only heat is allowed to flow through this boundary
        bdyFlux = 0.0;

        // conductive heat flux from boundary
        if constexpr (enableFullyImplicitThermal) {
            Evaluation heatFlux;
            // avoid overload of functions with same numeber of elements in eclproblem
            Scalar alpha;
            if constexpr (!std::is_empty_v<GetPropType<TypeTag, Properties::FluidSystem>>) {
                alpha = problem.getAlpha(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
            } else {
                alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
            }
            unsigned inIdx = 0; // dummy
            // always calculated with derivatives of this cell
            EnergyModule::ExtensiveQuantities::updateEnergyBoundary(heatFlux,
                                                                    insideIntQuants,
                                                                    /*focusDofIndex*/ inIdx,
                                                                    inIdx,
                                                                    alpha,
                                                                    bdyInfo.exFluidState);
            EnergyModule::addHeatFlux(bdyFlux, heatFlux);
        }

#ifndef NDEBUG
        for (unsigned i = 0; i < numEq; ++i) {
            Valgrind::CheckDefined(bdyFlux[i]);
        }
        Valgrind::CheckDefined(bdyFlux);
#endif
    }

    template <class BoundaryConditionData, class RateVectorLocal, class LocalProblem>
    OPM_HOST_DEVICE static void computeBoundaryFlux(RateVectorLocal& bdyFlux,
                                                    const LocalProblem& problem,
                                                    const BoundaryConditionData& bdyInfo,
                                                    const IntensiveQuantities& insideIntQuants,
                                                    unsigned globalSpaceIdx)
    {
        if constexpr (std::is_empty_v<GetPropType<TypeTag, Properties::FluidSystem>>) {
            switch (bdyInfo.type) {
            case BCType::NONE:
                bdyFlux = 0.0;
                break;
            case BCType::RATE:
                computeBoundaryFluxRate(bdyFlux, bdyInfo);
                break;
            case BCType::FREE:
            case BCType::DIRICHLET:
                computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
                break;
            case BCType::THERMAL:
                computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
                break;
            default:
                OPM_THROW(std::logic_error, "Unknown boundary condition type "
                              + std::to_string(static_cast<int>(bdyInfo.type))
                              + " in computeBoundaryFlux().");
            }
        } else { // Non-static fluid system used in GPU assembly
            switch (bdyInfo.type) {
            case BCType::NONE:
                bdyFlux = 0.0;
                break;
            case BCType::THERMAL:
                computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
                break;
            default:
                OPM_THROW(std::logic_error,
                          "Boundary condition type " + std::to_string(static_cast<int>(bdyInfo.type))
                              + " is not supported for GPU fluid systems in computeBoundaryFlux().");
            }
        }
    }

    static void computeSource(RateVector& source,
                              const Problem& problem,
                              const IntensiveQuantities& insideIntQuants,
                              unsigned globalSpaceIdex,
                              unsigned timeIdx)
    {
        OPM_TIMEBLOCK_LOCAL(computeSource, Subsystem::Assembly);
        // retrieve the source term intrinsic to the problem
        problem.source(source, globalSpaceIdex, timeIdx);

        // deal with bioeffects (if present)
        if constexpr (enableBioeffects) {
            BioeffectsModule::addSource(source, problem, insideIntQuants, globalSpaceIdex);
        }

        // scale the source term of the energy equation
        if constexpr (enableFullyImplicitThermal) {
            source[Indices::contiEnergyEqIdx]
                *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
        }
    }

    static void computeSourceDense(RateVector& source,
                                   const Problem& problem,
                                   const IntensiveQuantities& insideIntQuants,
                                   unsigned globalSpaceIdex,
                                   unsigned timeIdx)
    {
        source = 0.0;
        problem.addToSourceDense(source, globalSpaceIdex, timeIdx);

        // deal with bioeffects (if present)
        if constexpr (enableBioeffects) {
            BioeffectsModule::addSource(source, problem, insideIntQuants, globalSpaceIdex);
        }

        // scale the source term of the energy equation
        if constexpr (enableFullyImplicitThermal) {
            source[Indices::contiEnergyEqIdx]
                *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
        }
    }

    void computeSource(RateVector& source,
                       const ElementContext& elemCtx,
                       unsigned dofIdx,
                       unsigned timeIdx) const
    {
        OPM_TIMEBLOCK_LOCAL(computeSource, Subsystem::Assembly);
        // retrieve the source term intrinsic to the problem
        elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);

        // deal with bioeffects (if present)
        if constexpr (enableBioeffects) {
            BioeffectsModule::addSource(source, elemCtx, dofIdx, timeIdx);
        }

        // scale the source term of the energy equation
        if constexpr (enableFullyImplicitThermal) {
            source[Indices::contiEnergyEqIdx]
                *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
        }
    }

    template <class UpEval, class FluidState>
    static void evalPhaseFluxes_(RateVector& flux,
                                 unsigned phaseIdx,
                                 unsigned pvtRegionIdx,
                                 const ExtensiveQuantities& extQuants,
                                 const FluidState& upFs)
    {
        const auto& invB = getInvB_<FluidSystem, FluidState, UpEval>(upFs, phaseIdx, pvtRegionIdx);
        const auto& surfaceVolumeFlux = invB * extQuants.volumeFlux(phaseIdx);
        evalPhaseFluxes_<UpEval>(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, upFs);
    }

    template <class UpEval, class Eval, class FluidState, class RateVectorT = RateVector>
    OPM_HOST_DEVICE static void evalPhaseFluxes_(RateVectorT& flux,
                                                 unsigned phaseIdx,
                                                 unsigned pvtRegionIdx,
                                                 const Eval& surfaceVolumeFlux,
                                                 const FluidState& upFs)
    {
        const FluidSystem& fsys = upFs.fluidSystem();

        unsigned activeCompIdx
            = fsys.canonicalToActiveCompIdx(fsys.solventComponentIndex(phaseIdx));

        if constexpr (blackoilConserveSurfaceVolume) {
            flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux;
        } else {
            flux[conti0EqIdx + activeCompIdx]
                += surfaceVolumeFlux * fsys.referenceDensity(phaseIdx, pvtRegionIdx);
        }

        if (phaseIdx == oilPhaseIdx) {
            // dissolved gas (in the oil phase).
            if (fsys.enableDissolvedGas()) {
                const auto& Rs
                    = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeGasCompIdx = fsys.canonicalToActiveCompIdx(gasCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeGasCompIdx] += Rs * surfaceVolumeFlux;
                } else {
                    flux[conti0EqIdx + activeGasCompIdx] += Rs * surfaceVolumeFlux
                        * fsys.referenceDensity(gasPhaseIdx, pvtRegionIdx);
                }
            }
        } else if (phaseIdx == waterPhaseIdx) {
            // dissolved gas (in the water phase).
            if (fsys.enableDissolvedGasInWater()) {
                const auto& Rsw
                    = BlackOil::getRsw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeGasCompIdx = fsys.canonicalToActiveCompIdx(gasCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeGasCompIdx] += Rsw * surfaceVolumeFlux;
                } else {
                    flux[conti0EqIdx + activeGasCompIdx] += Rsw * surfaceVolumeFlux
                        * fsys.referenceDensity(gasPhaseIdx, pvtRegionIdx);
                }
            }
        } else if (phaseIdx == gasPhaseIdx) {
            // vaporized oil (in the gas phase).
            if (fsys.enableVaporizedOil()) {
                const auto& Rv
                    = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeOilCompIdx = fsys.canonicalToActiveCompIdx(oilCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeOilCompIdx] += Rv * surfaceVolumeFlux;
                } else {
                    flux[conti0EqIdx + activeOilCompIdx] += Rv * surfaceVolumeFlux
                        * fsys.referenceDensity(oilPhaseIdx, pvtRegionIdx);
                }
            }
            // vaporized water (in the gas phase).
            if (fsys.enableVaporizedWater()) {
                const auto& Rvw
                    = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeWaterCompIdx = fsys.canonicalToActiveCompIdx(waterCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeWaterCompIdx] += Rvw * surfaceVolumeFlux;
                } else {
                    flux[conti0EqIdx + activeWaterCompIdx] += Rvw * surfaceVolumeFlux
                        * fsys.referenceDensity(waterPhaseIdx, pvtRegionIdx);
                }
            }
        }
    }

    template <class ScalarVector>
    OPM_HOST_DEVICE static void adaptMassConservationQuantities_(ScalarVector& container,
                                                                 unsigned pvtRegionIdx)
    {
        // Delegate to the generic overload using a default-constructed static FluidSystem
        // instance. Valid because the static FluidSystem is stateless (std::is_empty_v).
        adaptMassConservationQuantities_(container, pvtRegionIdx, FluidSystem{});
    }

    template <class ScalarVector, class FsysType>
    OPM_HOST_DEVICE static void adaptMassConservationQuantities_(ScalarVector& container,
                                                                 unsigned pvtRegionIdx,
                                                                 const FsysType& fsys)
    {
        if constexpr (!blackoilConserveSurfaceVolume) {
            // convert "surface volume" to mass. this is complicated a bit by the fact that
            // not all phases are necessarily enabled. (we here assume that if a fluid phase
            // is disabled, its respective "main" component is not considered as well.)

            if constexpr (waterEnabled) {
                const unsigned activeWaterCompIdx = fsys.canonicalToActiveCompIdx(waterCompIdx);
                container[conti0EqIdx + activeWaterCompIdx]
                    *= fsys.referenceDensity(waterPhaseIdx, pvtRegionIdx);
            }

            if constexpr (gasEnabled) {
                const unsigned activeGasCompIdx = fsys.canonicalToActiveCompIdx(gasCompIdx);
                container[conti0EqIdx + activeGasCompIdx]
                    *= fsys.referenceDensity(gasPhaseIdx, pvtRegionIdx);
            }

            if constexpr (oilEnabled) {
                const unsigned activeOilCompIdx = fsys.canonicalToActiveCompIdx(oilCompIdx);
                container[conti0EqIdx + activeOilCompIdx]
                    *= fsys.referenceDensity(oilPhaseIdx, pvtRegionIdx);
            }
        }
    }

    // NNC does not have a direction
    static FaceDir::DirEnum faceDirFromDirId(const int dirId)
    {
        return dirId < 0 ? FaceDir::DirEnum::Unknown : FaceDir::FromIntersectionIndex(dirId);
    }
};

} // namespace Opm

#endif