factor out flux output

opm/models/blackoil/blackoilintensivequantities.hh

@@ -515,10 +515,13 @@ public:
         using Dir = FaceDir::DirEnum;
         if (dirMob_) {
             switch(facedir) {
+                case Dir::XMinus:
                 case Dir::XPlus:
                     return dirMob_->mobilityX_[phaseIdx];
+                case Dir::YMinus:
                 case Dir::YPlus:
                     return dirMob_->mobilityY_[phaseIdx];
+                case Dir::ZMinus:
                 case Dir::ZPlus:
                     return dirMob_->mobilityZ_[phaseIdx];
                 default:

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -272,8 +272,6 @@ public:
         const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx);
         Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
         Scalar faceArea = scvf.area();
-        const auto& materialLawManager = problem.materialLawManager();
-
         const auto dirid = scvf.dirId();
         Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx);
 
@@ -789,25 +787,11 @@ public:
 
     static FaceDir::DirEnum faceDirFromDirId(const int dirId)
     {
-        using Dir = FaceDir::DirEnum;
         // NNC does not have a direction
         if (dirId < 0 ) {
-            return Dir::Unknown;
+            return FaceDir::DirEnum::Unknown;
-        }
-        switch(dirId) {
-            case 0:
-            case 1:
-                return Dir::XPlus;
-            case 2:
-            case 3:
-                return Dir::YPlus;
-            case 4:
-            case 5:
-                return Dir::ZPlus;
-            default:
-                OPM_THROW(std::runtime_error,
-                            "Unexpected face id" + std::to_string(dirId));
         }
+        return FaceDir::FromIntersectionIndex(dirId);
     }
 };
 

opm/models/discretization/common/fvbaselinearizer.hh

@@ -311,6 +311,10 @@ public:
         // This linearizer stores no such data.
     }
 
+    void updateFlowsInfo() {
+        // This linearizer stores no such data.
+    }
+
     /*!
      * \brief Returns the map of constraint degrees of freedom.
      *

opm/models/discretization/common/tpfalinearizer.hh

@@ -434,8 +434,6 @@ private:
         unsigned numCells = model.numTotalDof();
         neighborInfo_.reserve(numCells, 6 * numCells);
         std::vector<NeighborInfo> loc_nbinfo;
-        const auto& materialLawManager = problem_().materialLawManager();
-        using FaceDirection = FaceDir::DirEnum;
         for (const auto& elem : elements(gridView_())) {
             stencil.update(elem);
 
@@ -643,6 +641,68 @@ public:
         }
     }
 
+    void updateFlowsInfo() {
+        OPM_TIMEBLOCK(updateFlows);
+        const bool& enableFlows = simulator_().problem().eclWriter()->eclOutputModule().hasFlows() ||
+                                    simulator_().problem().eclWriter()->eclOutputModule().hasBlockFlows();
+        const bool& enableFlores = simulator_().problem().eclWriter()->eclOutputModule().hasFlores();
+        if (!enableFlows && !enableFlores) {
+            return;
+        }
+        const unsigned int numCells = model_().numTotalDof();
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+        for (unsigned globI = 0; globI < numCells; ++globI) {
+            OPM_TIMEBLOCK_LOCAL(linearizationForEachCell);
+            const auto& nbInfos = neighborInfo_[globI];
+            ADVectorBlock adres(0.0);
+            ADVectorBlock darcyFlux(0.0);
+            const IntensiveQuantities& intQuantsIn = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
+            // Flux term.
+            {
+            OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell);
+            for (const auto& nbInfo : nbInfos) {
+                OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachFace);
+                unsigned globJ = nbInfo.neighbor;
+                assert(globJ != globI);
+                adres = 0.0;
+                darcyFlux = 0.0;
+                const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
+                LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo);
+                adres *= nbInfo.res_nbinfo.faceArea;
+                if (enableFlows) {
+                    for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
+                        flowsInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = adres[phaseIdx].value();
+                    }
+                }
+                if (enableFlores) {
+                    for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
+                        floresInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = darcyFlux[phaseIdx].value();
+                    }
+                }
+            }
+            }
+        }
+
+        // Boundary terms. Only looping over cells with nontrivial bcs.
+        for (const auto& bdyInfo : boundaryInfo_) {
+            ADVectorBlock adres(0.0);
+            const unsigned globI = bdyInfo.cell;
+            const IntensiveQuantities& insideIntQuants = model_().intensiveQuantities(globI, /*timeIdx*/ 0);
+            LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, insideIntQuants, globI);
+            adres *= bdyInfo.bcdata.faceArea;
+
+            if (enableFlows) {
+                for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
+                    flowsInfo_[globI][bdyInfo.dir].flow[phaseIdx] = adres[phaseIdx].value();
+                }
+            }
+            // TODO also store Flores?
+        }
+    }
+
 private:
     template <class SubDomainType>
     void linearize_(const SubDomainType& domain)
@@ -662,9 +722,6 @@ private:
         // the full system to zero, not just our part.
         // Instead, that must be called before starting the linearization.
         const bool& enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
-        const bool& enableFlows = simulator_().problem().eclWriter()->eclOutputModule().hasFlows() ||
-                                    simulator_().problem().eclWriter()->eclOutputModule().hasBlockFlows();
-        const bool& enableFlores = simulator_().problem().eclWriter()->eclOutputModule().hasFlores();
         const unsigned int numCells = domain.cells.size();
         const bool on_full_domain = (numCells == model_().numTotalDof());
 
@@ -701,16 +758,6 @@ private:
                         velocityInfo_[globI][loc].velocity[phaseIdx] = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
                     }
                 }
-                if (enableFlows) {
-                    for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
-                        flowsInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = adres[phaseIdx].value();
-                    }
-                }
-                if (enableFlores) {
-                    for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
-                        floresInfo_[globI][nbInfo.res_nbinfo.dirId].flow[phaseIdx] = darcyFlux[phaseIdx].value();
-                    }
-                }
                 setResAndJacobi(res, bMat, adres);
                 residual_[globI] += res;
                 //SparseAdapter syntax:  jacobian_->addToBlock(globI, globI, bMat);
@@ -809,23 +856,6 @@ private:
             residual_[globI] += res;
             ////SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
             *diagMatAddress_[globI] += bMat;
-
-            if (enableFlows) {
-                for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
-                    if (adres[phaseIdx].value() == 0)
-                        continue;
-
-                    flowsInfo_[globI][bdyInfo.dir].flow[phaseIdx] = adres[phaseIdx].value();
-                }
-            }
-            // if (enableFlores) {
-            //     for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
-            //         if (adres[phaseIdx].value() == 0)
-            //             continue;
-
-            //         floresInfo_[globI][bdyInfo.dir].flow[phaseIdx] = adres[phaseIdx].value();
-            //     }
-            // }
         }
     }