Finalize boundary treatment for TpfaLinearizer.

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -474,32 +474,35 @@ public:
             }
             const auto& pBoundary = bdyInfo.exFluidState.pressure(phaseIdx);
             const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx);
-
+            const unsigned pvtRegionIdx = insideIntQuants.pvtRegionIndex();
-            Evaluation surfaceVolumeFlux;// = invB * volumeFlux[SOME_INDEX];
 
             RateVector tmp;
 
             // mass conservation
-            if (pBoundary < pInside)
+            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) {
+            } else if (pBoundary > pInside) {
-                using RhsEval = typename std::conditional<std::is_same<typename FluidState::Scalar, Evaluation>::value,
-                                                          Evaluation, Scalar>::type;
                 // influx
-                evalPhaseFluxes_<RhsEval>(tmp,
+                using ScalarFluidState = decltype(bdyInfo.exFluidState);
-                                          phaseIdx,
+                const auto& invB = getInvB_<FluidSystem, ScalarFluidState, Scalar>(bdyInfo.exFluidState, phaseIdx, pvtRegionIdx);
-                                          insideIntQuants.pvtRegionIndex(),
+                Evaluation surfaceVolumeFlux = invB * volumeFlux[phaseIdx];
-                                          surfaceVolumeFlux,
+                evalPhaseFluxes_<Scalar>(tmp,
-                                          bdyInfo.exFluidState);
+                                         phaseIdx,
+                                         insideIntQuants.pvtRegionIndex(),
+                                         surfaceVolumeFlux,
+                                         bdyInfo.exFluidState);
             }
 
-            for (unsigned i = 0; i < tmp.size(); ++i)
+            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.

opm/models/discretization/common/tpfalinearizer.hh

@@ -83,7 +83,6 @@ 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;
@@ -92,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>;
@@ -344,6 +344,28 @@ private:
                     }
                 }
                 neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
+                for (size_t 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 (int 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});
+                    }
+                }
             }
         }
 
@@ -472,6 +494,7 @@ private:
                 throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI));
             }
             LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, *insideIntQuants, globI, 0);
+            adres *= bdyInfo.bcdata.faceArea;
             setResAndJacobi(res, bMat, adres);
             residual_[globI] += res;
             jacobian_->addToBlock(globI, globI, bMat);
@@ -497,6 +520,7 @@ private:
     };
     SparseTable<NeighborInfo> neighborInfo_;
 
+    using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;
     struct BoundaryConditionData
     {
         BCType type;
@@ -505,7 +529,7 @@ private:
         unsigned boundaryFaceIndex;
         double faceArea;
         double faceZCoord;
-        FluidState exFluidState;
+        ScalarFluidState exFluidState;
     };
     struct BoundaryInfo
     {