Merge pull request #4119 from atgeirr/refactor-boundary-flux

Refactor boundary flux

ebos/eclfluxmodule.hh

@@ -43,6 +43,8 @@
 #include <dune/common/fmatrix.hh>
 #include <fmt/format.h>
 
+#include <array>
+
 namespace Opm {
 
 template <class TypeTag>
@@ -212,8 +214,8 @@ protected:
 
 public:
 
-    static void volumeAndPhasePressureDifferences(short (&upIdx)[numPhases],
+    static void volumeAndPhasePressureDifferences(std::array<short, numPhases>& upIdx,
-                                                  short (&dnIdx)[numPhases],
+                                                  std::array<short, numPhases>& dnIdx,
                                                   Evaluation (&volumeFlux)[numPhases],
                                                   Evaluation (&pressureDifferences)[numPhases],
                                                   const ElementContext& elemCtx,
@@ -418,34 +420,75 @@ protected:
                                      unsigned timeIdx,
                                      const FluidState& exFluidState)
     {
+        const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
+        const Scalar faceArea = scvf.area();
+        const Scalar zEx = scvf.integrationPos()[dimWorld - 1];
         const auto& problem = elemCtx.problem();
+        const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(0, timeIdx);
+        const auto& intQuantsIn = elemCtx.intensiveQuantities(0, timeIdx);
+
+        calculateBoundaryGradients_(problem,
+                                    globalSpaceIdx,
+                                    intQuantsIn,
+                                    scvfIdx,
+                                    faceArea,
+                                    zEx,
+                                    exFluidState,
+                                    upIdx_,
+                                    dnIdx_,
+                                    volumeFlux_,
+                                    pressureDifference_);
+
+        // Treating solvent here and not in the static method, since that would require more
+        // extensive refactoring. It means that the TpfaLinearizer will not support bcs for solvent until this is
+        // addressed.
+        if constexpr (enableSolvent) {
+            if (upIdx_[gasPhaseIdx] == 0) {
+                const Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, scvfIdx);
+                const Scalar transModified = trans * Toolbox::value(intQuantsIn.rockCompTransMultiplier());
+                const auto solventFlux = pressureDifference_[gasPhaseIdx] * intQuantsIn.mobility(gasPhaseIdx) * (-transModified/faceArea);
+                asImp_().setSolventVolumeFlux(solventFlux);
+            } else {
+                asImp_().setSolventVolumeFlux(0.0);
+            }
+        }
+    }
+
+public:
+    /*!
+     * \brief Update the required gradients for boundary faces
+     */
+    template <class Problem, class FluidState, class EvaluationContainer>
+    static void calculateBoundaryGradients_(const Problem& problem,
+                                            const unsigned globalSpaceIdx,
+                                            const IntensiveQuantities& intQuantsIn,
+                                            const unsigned bfIdx,
+                                            const double faceArea,
+                                            const double zEx,
+                                            const FluidState& exFluidState,
+                                            std::array<short, numPhases>& upIdx,
+                                            std::array<short, numPhases>& dnIdx,
+                                            EvaluationContainer& volumeFlux,
+                                            EvaluationContainer& pressureDifference)
+    {
 
         bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
         if (!enableBoundaryMassFlux)
             return;
 
-        const auto& stencil = elemCtx.stencil(timeIdx);
+        Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx);
-        const auto& scvf = stencil.boundaryFace(scvfIdx);
-
-        unsigned interiorDofIdx = scvf.interiorIndex();
-
-        Scalar trans = problem.transmissibilityBoundary(elemCtx, scvfIdx);
-        Scalar faceArea = scvf.area();
 
         // 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.)
-        Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
+        Scalar g = problem.gravity()[dimWorld - 1];
-
-        const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, 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...
-        Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
+        Scalar zIn = problem.dofCenterDepth(globalSpaceIdx);
-        Scalar zEx = scvf.integrationPos()[dimWorld - 1];
 
         // the distances from the DOF's depths. (i.e., the additional depth of the
         // exterior DOF)
@@ -465,62 +508,54 @@ protected:
             Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
             pressureExterior += rhoAvg*(distZ*g);
 
-            pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;
+            pressureDifference[phaseIdx] = pressureExterior - pressureInterior;
 
             // decide the upstream index for the phase. for this we make sure that the
             // degree of freedom which is regarded upstream if both pressures are equal
             // is always the same: if the pressure is equal, the DOF with the lower
             // global index is regarded to be the upstream one.
-            if (pressureDifference_[phaseIdx] > 0.0) {
+            const unsigned interiorDofIdx = 0; // Valid only for cell-centered FV.
-                upIdx_[phaseIdx] = -1;
+            if (pressureDifference[phaseIdx] > 0.0) {
-                dnIdx_[phaseIdx] = interiorDofIdx;
+                upIdx[phaseIdx] = -1;
+                dnIdx[phaseIdx] = interiorDofIdx;
             }
             else {
-                upIdx_[phaseIdx] = interiorDofIdx;
+                upIdx[phaseIdx] = interiorDofIdx;
-                dnIdx_[phaseIdx] = -1;
+                dnIdx[phaseIdx] = -1;
             }
 
             Evaluation transModified = trans;
 
-            unsigned upstreamIdx = upstreamIndex_(phaseIdx);
+            if (upIdx[phaseIdx] == interiorDofIdx) {
-            if (upstreamIdx == interiorDofIdx) {
 
                 // this is slightly hacky because in the automatic differentiation case, it
                 // only works for the element centered finite volume method. for ebos this
                 // does not matter, though.
-                const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
+                const auto& up = intQuantsIn;
 
                 // deal with water induced rock compaction
-                const double transMult = Toolbox::value(up.rockCompTransMultiplier());
+                const Scalar transMult = Toolbox::value(up.rockCompTransMultiplier());
                 transModified *= transMult;
 
-                volumeFlux_[phaseIdx] =
+                volumeFlux[phaseIdx] =
-                    pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea);
+                    pressureDifference[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea);
-
-                if (enableSolvent && phaseIdx == gasPhaseIdx)
-                    asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-transModified/faceArea));
             }
             else {
                 // compute the phase mobility using the material law parameters of the
                 // interior element. TODO: this could probably be done more efficiently
-                const auto& matParams =
+                const auto& matParams = problem.materialLawParams(globalSpaceIdx);
-                    elemCtx.problem().materialLawParams(elemCtx,
-                                                        interiorDofIdx,
-                                                        /*timeIdx=*/0);
                 std::array<typename FluidState::Scalar,numPhases> kr;
                 MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
 
                 const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
-                volumeFlux_[phaseIdx] =
+                volumeFlux[phaseIdx] =
-                    pressureDifference_[phaseIdx]*mob*(-transModified/faceArea);
+                    pressureDifference[phaseIdx]*mob*(-transModified/faceArea);
-
-                // Solvent inflow is not yet supported
-                if (enableSolvent && phaseIdx == gasPhaseIdx)
-                    asImp_().setSolventVolumeFlux(0.0);
             }
         }
     }
 
+protected:
+
     /*!
      * \brief Update the volumetric fluxes for all fluid phases on the interior faces of the context
      */
@@ -545,9 +580,9 @@ private:
     Evaluation pressureDifference_[numPhases];
 
     // the local indices of the interior and exterior degrees of freedom
-    short upIdx_[numPhases];
+    std::array<short, numPhases> upIdx_;
-    short dnIdx_[numPhases];
+    std::array<short, numPhases> dnIdx_;
-};
+ };
 
 } // namespace Opm
 

ebos/eclproblem.hh

@@ -1340,6 +1340,15 @@ public:
         return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
     }
 
+    /*!
+     * \brief Direct access to a boundary transmissibility.
+     */
+    Scalar transmissibilityBoundary(const unsigned globalSpaceIdx,
+                                    const unsigned boundaryFaceIdx) const
+    {
+        return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
+    }
+
     /*!
      * \copydoc EclTransmissiblity::thermalHalfTransmissibility
      */
@@ -2048,6 +2057,29 @@ public:
         return this->rockCompTransMultWc_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true);
     }
 
+    std::pair<bool, RateVector> boundaryCondition(const unsigned int globalSpaceIdx, const int directionId)
+    {
+        if (!nonTrivialBoundaryConditions_) {
+            return { false, RateVector(0.0) };
+        }
+        switch (directionId) {
+        case 0:
+            return { freebcXMinus_[globalSpaceIdx], massratebcXMinus_[globalSpaceIdx] };
+        case 1:
+            return { freebcX_[globalSpaceIdx], massratebcX_[globalSpaceIdx] };
+        case 2:
+            return { freebcYMinus_[globalSpaceIdx], massratebcYMinus_[globalSpaceIdx] };
+        case 3:
+            return { freebcY_[globalSpaceIdx], massratebcY_[globalSpaceIdx] };
+        case 4:
+            return { freebcZMinus_[globalSpaceIdx], massratebcZMinus_[globalSpaceIdx] };
+        case 5:
+            return { freebcZ_[globalSpaceIdx], massratebcZ_[globalSpaceIdx] };
+        default:
+            return { false, RateVector(0.0) };
+        }
+    }
+
 private:
     // update the parameters needed for DRSDT and DRVDT
     void updateCompositionChangeLimits_()