Add boundary condition treatment to TpfaLinearizer.

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -61,6 +61,7 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     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>() };
@@ -414,6 +415,114 @@ public:
 
     }
 
+    template <class BoundaryConditionData>
+    static void computeBoundaryFlux(RateVector& bdyFlux,
+                                    const Problem& problem,
+                                    const BoundaryConditionData& bdyInfo,
+                                    const IntensiveQuantities& insideIntQuants,
+                                    unsigned globalSpaceIdx,
+                                    unsigned timeIdx)
+    {
+        if (bdyInfo.type == BCType::RATE) {
+            computeBoundaryFluxRate(bdyFlux, bdyInfo);
+        } else if (bdyInfo.type == BCType::FREE) {
+            computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx, timeIdx);
+        } else {
+            throw std::logic_error("Unknown boundary condition type in computeBoundaryFlux().");
+        }
+    }
+
+    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,
+                                        unsigned timeIdx)
+    {
+        short upIdx[numPhases];
+        short dnIdx[numPhases];
+        RateVector volumeFlux;
+        RateVector pressureDifference;
+        ExtensiveQuantities::calculateBoundaryGradients_(problem,
+                                                         globalSpaceIdx,
+                                                         insideIntQuants,
+                                                         bdyInfo.boundaryFaceIndex,
+                                                         timeIdx,
+                                                         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);
+
+            Evaluation surfaceVolumeFlux;// = invB * volumeFlux[SOME_INDEX];
+
+            RateVector tmp;
+
+            // mass conservation
+            if (pBoundary < pInside)
+                // outflux
+                evalPhaseFluxes_<Evaluation>(tmp,
+                                             phaseIdx,
+                                             insideIntQuants.pvtRegionIndex(),
+                                             surfaceVolumeFlux,
+                                             insideIntQuants.fluidState());
+            else if (pBoundary > pInside) {
+                using RhsEval = typename std::conditional<std::is_same<typename FluidState::Scalar, Evaluation>::value,
+                                                          Evaluation, Scalar>::type;
+                // influx
+                evalPhaseFluxes_<RhsEval>(tmp,
+                                          phaseIdx,
+                                          insideIntQuants.pvtRegionIndex(),
+                                          surfaceVolumeFlux,
+                                          bdyInfo.exFluidState);
+            }
+
+            for (unsigned i = 0; i < tmp.size(); ++i)
+                bdyFlux[i] += tmp[i];
+
+            static_assert(!enableEnergy, "Relevant treatment of boundary conditions must be implemented before enabling.");
+            // Add energy flux treatment per phase here.
+        }
+
+        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.");
+        static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling.");
+
+        // make sure that the right mass conservation quantities are used
+        adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex());
+
+        // heat conduction
+        static_assert(!enableEnergy, "Relevant treatment of boundary conditions must be implemented before enabling.");
+
+#ifndef NDEBUG
+        for (unsigned i = 0; i < numEq; ++i) {
+            Valgrind::CheckDefined(bdyFlux[i]);
+        }
+        Valgrind::CheckDefined(bdyFlux);
+#endif
+    }
+
     static void computeSource(RateVector& source,
                               const Problem& problem,
                               unsigned globalSpaceIdex,

opm/models/discretization/common/tpfalinearizer.hh

@@ -83,6 +83,7 @@ class TpfaLinearizer
     using Stencil = GetPropType<TypeTag, Properties::Stencil>;
     using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using FluidState = typename IntensiveQuantities::FluidState;
 
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementIterator = typename GridView::template Codim<0>::Iterator;
@@ -458,6 +459,22 @@ private:
                 residual_[globI] += res;
                 jacobian_->addToBlock(globI, globI, bMat);
             }
+        } // end of loop for cell globI.
+
+        // Boundary terms. Only looping over cells with nontrivial bcs.
+        for (const auto& bdyInfo : boundaryInfo_) {
+            VectorBlock res(0.0);
+            MatrixBlock bMat(0.0);
+            ADVectorBlock adres(0.0);
+            const unsigned globI = bdyInfo.cell;
+            const IntensiveQuantities* insideIntQuants = model_().cachedIntensiveQuantities(globI, /*timeIdx*/ 0);
+            if (insideIntQuants == nullptr) {
+                throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI));
+            }
+            LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, *insideIntQuants, globI, 0);
+            setResAndJacobi(res, bMat, adres);
+            residual_[globI] += res;
+            jacobian_->addToBlock(globI, globI, bMat);
         }
     }
 
@@ -479,6 +496,23 @@ private:
         FaceDir::DirEnum faceDirection;
     };
     SparseTable<NeighborInfo> neighborInfo_;
+
+    struct BoundaryConditionData
+    {
+        BCType type;
+        VectorBlock massRate;
+        unsigned pvtRegionIdx;
+        unsigned boundaryFaceIndex;
+        double faceArea;
+        double faceZCoord;
+        FluidState exFluidState;
+    };
+    struct BoundaryInfo
+    {
+        unsigned int cell;
+        BoundaryConditionData bcdata;
+    };
+    std::vector<BoundaryInfo> boundaryInfo_;
 };
 
 } // namespace Opm