ebos: refactor the way SWATINIT is dealt with

we now store the maximum oil-water capillary pressure and apply it to
the material parameters later. this simplifies things within the
EquilInitializer somewhat and also allows external code (i.e.,
flow_ebos) to choose when SWATINIT gets applied.

ebos/eclequilinitializer.hh

@@ -85,9 +85,9 @@ class EclEquilInitializer
     enum { dimWorld = GridView::dimensionworld };
 
 public:
-    template <class MaterialLawManager>
+    template <class MaxPcnwVector>
     EclEquilInitializer(const Simulator& simulator,
-                        std::shared_ptr<MaterialLawManager> materialLawManager)
+                        MaxPcnwVector& maxPcnw)
         : simulator_(simulator)
     {
         const auto& gridManager = simulator.gridManager();
@@ -153,6 +153,10 @@ public:
             localToEquilIndex[ elemIdx ] = equilCartesianToCompressed[ cartesianIndex ];
         }
 
+        // resize the array of maximum water-oil capillary pressures (needed to apply the
+        // SWATINIT keyword).
+        maxPcnw.resize(numElems);
+
         // copy the result into the array of initial fluid states
         initialFluidStates_.resize(numCartesianElems);
         for (unsigned int elemIdx = 0; elemIdx < numElems; ++elemIdx) {
@@ -248,35 +252,18 @@ public:
                 fluidState.setMoleFraction(gasPhaseIdx, oilCompIdx, xgO);
                 fluidState.setMoleFraction(gasPhaseIdx, gasCompIdx, 1 - xgO);
             }
+
+            // store the maximum oil-water capillary pressure for later (required by the
+            // SWATINIT keyword)
+            const auto& equilScalingPoints =
+                equilMaterialLawManager->oilWaterScaledEpsPointsDrainage(equilElemIdx);
+            maxPcnw[elemIdx] = equilScalingPoints.maxPcnw();
         }
 
-        // deal with the capillary pressure modification due to SWATINIT. this is
+        // immediately forget about the SWATINIT stuff if the SWATINIT keyword is not
-        // only necessary because, the fine equilibration code from opm-core requires
+        // present in the deck
-        // its own grid and its own material law manager...
+        if (!deck.hasKeyword("SWATINIT"))
-        if (GET_PROP_VALUE(TypeTag, EnableSwatinit)) {
+            maxPcnw.clear();
-            std::vector<int> cartesianToCompressedElemIdx(numCartesianElems, -1);
-            for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
-                int cartElemIdx = gridManager.cartesianIndex(elemIdx);
-                cartesianToCompressedElemIdx[cartElemIdx] = elemIdx;
-            }
-
-            for (unsigned equilElemIdx = 0; equilElemIdx < numEquilElems; ++equilElemIdx) {
-                int cartElemIdx = gridManager.equilCartesianIndex(equilElemIdx);
-                assert(cartElemIdx >= 0);
-                int elemIdx = cartesianToCompressedElemIdx[cartElemIdx];
-                if (elemIdx < 0)
-                    // the element is present in the grid for used for equilibration but
-                    // it isn't present in the one used for the simulation. the most
-                    // probable reason for this is that the simulation grid was load
-                    // balanced.
-                    continue;
-
-                auto& scalingPoints = materialLawManager->oilWaterScaledEpsPointsDrainage(elemIdx);
-                const auto& equilScalingPoints = equilMaterialLawManager->oilWaterScaledEpsPointsDrainage(equilElemIdx);
-
-                scalingPoints.setMaxPcnw(equilScalingPoints.maxPcnw());
-            }
-        }
     }
 
     /*!

ebos/eclproblem.hh

@@ -851,6 +851,27 @@ public:
     const EclWellManager<TypeTag>& wellManager() const
     { return wellManager_; }
 
+    /*!
+     * \brief Apply the necessary measures mandated by the SWATINIT keyword the the
+     *        material law parameters.
+     *
+     * If SWATINIT is not used or if the method has already been called, this method is a
+     * no-op.
+     */
+    void applySwatinit()
+    {
+        if (maxPcnw_.empty())
+            return; // SWATINIT not applicable
+
+        unsigned numElems = this->simulator().gridView().size(/*codim=*/0);
+        for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
+            auto& scalingPoints = materialLawManager_->oilWaterScaledEpsPointsDrainage(elemIdx);
+            scalingPoints.setMaxPcnw(maxPcnw_[elemIdx]);
+        }
+
+        maxPcnw_.clear();
+    }
+
 private:
     static bool enableEclOutput_()
     { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
@@ -1032,14 +1053,17 @@ private:
 
         // release the memory of the EQUIL grid since it's no longer needed after this point
         this->simulator().gridManager().releaseEquilGrid();
+
+        // apply SWATINIT if requested by the problem and if it is enabled by the deck
+        if (GET_PROP_VALUE(TypeTag, EnableSwatinit))
+            applySwatinit();
     }
 
     void readEquilInitialCondition_()
     {
-        // The EQUIL initializer also modifies the material law manager according to
+        // initial condition corresponds to hydrostatic conditions
-        // SWATINIT (although it does not belong there strictly speaking)
         typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
-        EquilInitializer equilInitializer(this->simulator(), materialLawManager_);
+        EquilInitializer equilInitializer(this->simulator(), maxPcnw_);
 
         // since the EquilInitializer provides fluid states that are consistent with the
         // black-oil model, we can use naive instead of mass conservative determination
@@ -1341,6 +1365,8 @@ private:
     EclSummaryWriter summaryWriter_;
 
     PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
+
+    std::vector<Scalar> maxPcnw_;
 };
 } // namespace Ewoms