Merge pull request #718 from atgeirr/tpfa-bcs

Boundary condition support for TpfaLinearizer

opm/models/blackoil/blackoilintensivequantities.hh

@@ -131,6 +131,15 @@ public:
                                           enableBrine,
                                           enableSaltPrecipitation,
                                           Indices::numPhases>;
+    using ScalarFluidState = BlackOilFluidState<Scalar,
+                                                FluidSystem,
+                                                enableTemperature,
+                                                enableEnergy,
+                                                compositionSwitchEnabled,
+                                                enableEvaporation,
+                                                enableBrine,
+                                                enableSaltPrecipitation,
+                                                Indices::numPhases>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     BlackOilIntensiveQuantities()

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)
+    {
+        if (bdyInfo.type == BCType::RATE) {
+            computeBoundaryFluxRate(bdyFlux, bdyInfo);
+        } else if (bdyInfo.type == BCType::FREE) {
+            computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
+        } else {
+            throw std::logic_error("Unknown boundary condition type " + std::to_string(static_cast<int>(bdyInfo.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)
+    {
+        std::array<short, numPhases> upIdx;
+        std::array<short, numPhases> dnIdx;
+        RateVector volumeFlux;
+        RateVector 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;
+
+            // mass conservation
+            if (pBoundary < pInside) {
+                // outflux
+                const auto& invB = getInvB_<FluidSystem, FluidState, Evaluation>(insideIntQuants.fluidState(), phaseIdx, pvtRegionIdx);
+                Evaluation surfaceVolumeFlux = invB * volumeFlux[phaseIdx];
+                evalPhaseFluxes_<Evaluation>(tmp,
+                                             phaseIdx,
+                                             insideIntQuants.pvtRegionIndex(),
+                                             surfaceVolumeFlux,
+                                             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 * volumeFlux[phaseIdx];
+                evalPhaseFluxes_<Scalar>(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

@@ -91,6 +91,7 @@ class TpfaLinearizer
 
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>() };
+    enum { dimWorld = GridView::dimensionworld };
 
     using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock;
     using VectorBlock = Dune::FieldVector<Scalar, numEq>;
@@ -343,6 +344,28 @@ private:
                     }
                 }
                 neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
+                for (unsigned bfIndex = 0; bfIndex < stencil.numBoundaryFaces(); ++bfIndex) {
+                    const auto& bf = stencil.boundaryFace(bfIndex);
+                    const int dir_id = bf.dirId();
+                    const auto [free, massrateAD] = problem_().boundaryCondition(myIdx, dir_id);
+                    // Strip the unnecessary (and zero anyway) derivatives off massrate.
+                    VectorBlock massrate(0.0);
+                    for (size_t ii = 0; ii < massrate.size(); ++ii) {
+                        massrate[ii] = massrateAD[ii].value();
+                    }
+                    const bool nonzero_massrate = massrate != VectorBlock(0.0);
+                    if (free || nonzero_massrate) {
+                        const auto& exFluidState = problem_().initialFluidState(myIdx);
+                        BoundaryConditionData bcdata{free ? BCType::FREE : BCType::RATE,
+                                                     massrate,
+                                                     exFluidState.pvtRegionIndex(),
+                                                     bfIndex,
+                                                     bf.area(),
+                                                     bf.integrationPos()[dimWorld - 1],
+                                                     exFluidState};
+                        boundaryInfo_.push_back({myIdx, bcdata});
+                    }
+                }
             }
         }
 
@@ -458,6 +481,23 @@ 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);
+            adres *= bdyInfo.bcdata.faceArea;
+            setResAndJacobi(res, bMat, adres);
+            residual_[globI] += res;
+            jacobian_->addToBlock(globI, globI, bMat);
         }
     }
 
@@ -479,6 +519,24 @@ private:
         FaceDir::DirEnum faceDirection;
     };
     SparseTable<NeighborInfo> neighborInfo_;
+
+    using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;
+    struct BoundaryConditionData
+    {
+        BCType type;
+        VectorBlock massRate;
+        unsigned pvtRegionIdx;
+        unsigned boundaryFaceIndex;
+        double faceArea;
+        double faceZCoord;
+        ScalarFluidState exFluidState;
+    };
+    struct BoundaryInfo
+    {
+        unsigned int cell;
+        BoundaryConditionData bcdata;
+    };
+    std::vector<BoundaryInfo> boundaryInfo_;
 };
 
 } // namespace Opm

opm/models/discretization/ecfv/ecfvstencil.hh

@@ -212,7 +212,10 @@ public:
         { return area_; }
 
         /*!
-         * \brief Returns the direction of the face
+         * \brief Returns the direction id of the face w.r.t the cell.
+         *
+         * For corner point grids, this is 0-5 for I-, I+, J-, J+, K- and K+ faces,
+         * and -1 for NNC faces.
          */
         int dirId() const
         {