Prepare flow for ebos initialization

Apply swatInit in the initialization
Stop using the equilGrid in the initialization code
Keep The initialFluidState until end of first time step to make it
possible for flow to output it.

ebos/eclequilinitializer.hh

@@ -48,7 +48,6 @@ NEW_PROP_TAG(FluidSystem);
 NEW_PROP_TAG(GridView);
 NEW_PROP_TAG(Scalar);
 NEW_PROP_TAG(MaterialLaw);
-NEW_PROP_TAG(EnableSwatinit);
 }
 
 /*!
@@ -87,54 +86,19 @@ class EclEquilInitializer
 public:
     template <class EclMaterialLawManager>
     EclEquilInitializer(const Simulator& simulator,
-                        EclMaterialLawManager& materialLawManager,
+                        EclMaterialLawManager& materialLawManager)
-                        bool enableSwatinit)
         : simulator_(simulator)
     {
         const auto& gridManager = simulator.gridManager();
-        const auto& deck = gridManager.deck();
-        const auto& eclState = gridManager.eclState();
-        const auto& equilGrid = gridManager.equilGrid();
 
         unsigned numElems = gridManager.grid().size(0);
-        unsigned numEquilElems = gridManager.equilGrid().size(0);
         unsigned numCartesianElems = gridManager.cartesianSize();
         typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-        typedef Opm::ThreePhaseMaterialTraits<double,
-                                              /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
-                                              /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
-                                              /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> EquilTraits;
 
-        // create a separate instance of the material law manager just because opm-core
+        Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> initialState(materialLawManager,
-        // only supports double as the type for scalars (but ebos may use float or quad)
-        std::vector<int> compressedToCartesianEquilElemIdx(numEquilElems);
-        std::vector<int> equilCartesianToCompressed( gridManager.equilCartesianSize(), -1 );
-
-        for (unsigned equilElemIdx = 0; equilElemIdx < numEquilElems; ++equilElemIdx)
-        {
-            unsigned int equilCartesianIdx = gridManager.equilCartesianIndex(equilElemIdx);
-            compressedToCartesianEquilElemIdx[equilElemIdx] = equilCartesianIdx;
-            equilCartesianToCompressed[ equilCartesianIdx ] = equilElemIdx;
-        }
-
-        Opm::EclMaterialLawManager<EquilTraits> equilMaterialLawManager =
-            Opm::EclMaterialLawManager<EquilTraits>();
-        equilMaterialLawManager.initFromDeck(deck, eclState, compressedToCartesianEquilElemIdx);
-
-        Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> initialState(equilMaterialLawManager,
                                                                                   gridManager.eclState(),
-                                                                                  equilGrid,
+                                                                                  gridManager.grid(),
-                                                                                  simulator.problem().gravity()[dimWorld - 1],
+                                                                                  simulator.problem().gravity()[dimWorld - 1]);
-                                                                                  enableSwatinit);
-
-
-        std::vector<int> localToEquilIndex( numElems, -1 );
-        for( unsigned int elemIdx = 0; elemIdx < numElems; ++elemIdx )
-        {
-            const int cartesianIndex = gridManager.cartesianIndex( elemIdx );
-            assert( equilCartesianToCompressed[ cartesianIndex ] >= 0 );
-            localToEquilIndex[ elemIdx ] = equilCartesianToCompressed[ cartesianIndex ];
-        }
 
         // copy the result into the array of initial fluid states
         initialFluidStates_.resize(numCartesianElems);
@@ -142,8 +106,6 @@ public:
             unsigned cartesianElemIdx = gridManager.cartesianIndex(elemIdx);
             auto& fluidState = initialFluidStates_[cartesianElemIdx];
 
-            const unsigned int equilElemIdx = localToEquilIndex[ elemIdx ];
-
             // get the PVT region index of the current element
             unsigned regionIdx = simulator_.problem().pvtRegionIndex(elemIdx);
 
@@ -153,7 +115,7 @@ public:
                 if (!FluidSystem::phaseIsActive(phaseIdx))
                     S = 0.0;
                 else {
-                    S = initialState.saturation()[phaseIdx][equilElemIdx];
+                    S = initialState.saturation()[phaseIdx][elemIdx];
                 }
                 fluidState.setSaturation(phaseIdx, S);
             }
@@ -164,7 +126,7 @@ public:
 
             // set the phase pressures.
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                fluidState.setPressure(phaseIdx, initialState.press()[phaseIdx][equilElemIdx]);
+                fluidState.setPressure(phaseIdx, initialState.press()[phaseIdx][elemIdx]);
 
             // reset the phase compositions
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@@ -178,7 +140,7 @@ public:
             if (FluidSystem::enableDissolvedGas()) {
                 // for gas and oil we have to translate surface volumes to mole fractions
                 // before we can set the composition in the fluid state
-                Scalar Rs = initialState.rs()[equilElemIdx];
+                Scalar Rs = initialState.rs()[elemIdx];
                 Scalar RsSat = FluidSystem::saturatedDissolutionFactor(fluidState, oilPhaseIdx, regionIdx);
 
                 if (Rs > RsSat)
@@ -194,7 +156,7 @@ public:
 
             // retrieve the surface volume of vaporized gas
             if (FluidSystem::enableVaporizedOil()) {
-                Scalar Rv = initialState.rv()[equilElemIdx];
+                Scalar Rv = initialState.rv()[elemIdx];
                 Scalar RvSat = FluidSystem::saturatedDissolutionFactor(fluidState, gasPhaseIdx, regionIdx);
 
                 if (Rv > RvSat)
@@ -207,16 +169,6 @@ public:
                 fluidState.setMoleFraction(gasPhaseIdx, oilCompIdx, xgO);
                 fluidState.setMoleFraction(gasPhaseIdx, gasCompIdx, 1 - xgO);
             }
-
-            // deal with the changed pressure scaling due to SWATINIT if SWATINIT is
-            // requested to be applied. this is quite hacky but hey it works!
-            if (enableSwatinit) {
-                const auto& equilScalingPoints =
-                    equilMaterialLawManager.oilWaterScaledEpsPointsDrainage(equilElemIdx);
-                auto& scalingPoints =
-                    materialLawManager.oilWaterScaledEpsPointsDrainage(elemIdx);
-                scalingPoints.setMaxPcnw(equilScalingPoints.maxPcnw());
-            }
         }
     }
 

ebos/eclproblem.hh

@@ -119,9 +119,6 @@ NEW_PROP_TAG(DisableWells);
 // print statements in debug mode.
 NEW_PROP_TAG(EnableDebuggingChecks);
 
-// If this property is set to false, the SWATINIT keyword will not be handled by ebos.
-NEW_PROP_TAG(EnableSwatinit);
-
 // Set the problem property
 SET_TYPE_PROP(EclBaseProblem, Problem, Ewoms::EclProblem<TypeTag>);
 
@@ -238,8 +235,6 @@ SET_BOOL_PROP(EclBaseProblem, DisableWells, false);
 // By default, we enable the debugging checks if we're compiled in debug mode
 SET_BOOL_PROP(EclBaseProblem, EnableDebuggingChecks, true);
 
-// ebos handles the SWATINIT keyword by default
-SET_BOOL_PROP(EclBaseProblem, EnableSwatinit, true);
 } // namespace Properties
 
 /*!
@@ -556,6 +551,10 @@ public:
             // write the summary information after each time step
             summaryWriter_.write(wellManager_);
         }
+
+        // we no longer need the initial soluiton
+        if (this->simulator().episodeIndex() == 0)
+            initialFluidStates_.clear();
     }
 
     /*!
@@ -972,26 +971,8 @@ public:
         // the initial solution.
         thresholdPressures_.finishInit();
 
-        // apply SWATINIT if requested by the programmer. this is only necessary if
-        // SWATINIT has not yet been considered at the first time the EquilInitializer
-        // was used, i.e., only if threshold pressures are enabled in addition to
-        // SWATINIT.
-        const auto& deck = this->simulator().gridManager().deck();
-        const auto& eclState = this->simulator().gridManager().eclState();
-        int numEquilRegions = eclState.getTableManager().getEqldims().getNumEquilRegions();
-        bool useThpres = deck.hasKeyword("THPRES") && numEquilRegions > 1;
-        bool useSwatinit =
-            GET_PROP_VALUE(TypeTag, EnableSwatinit) &&
-            eclState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT");
-
-        if (useThpres && useSwatinit)
-            applySwatinit();
-
         // release the memory of the EQUIL grid since it's no longer needed after this point
         this->simulator().gridManager().releaseEquilGrid();
-
-        // the fluid states specifying the initial condition are also no longer required.
-        initialFluidStates_.clear();
     }
 
     /*!
@@ -1026,48 +1007,9 @@ public:
     const EclWellManager<TypeTag>& wellManager() const
     { return wellManager_; }
 
-    /*!
+    // temporary solution to facilitate output of initial state from flow
-     * \brief Apply the necessary measures mandated by the SWATINIT keyword the the
+    const ScalarFluidState& initialFluidState(unsigned globalDofIdx ) const {
-     *        material law parameters.
+        return initialFluidStates_[globalDofIdx];
-     *
-     * If SWATINIT is not used or if the method has already been called, this method is a
-     * no-op.
-     */
-    void applySwatinit()
-    {
-        const auto& deck = this->simulator().gridManager().deck();
-        const auto& eclState = this->simulator().gridManager().eclState();
-        const auto& deckProps = eclState.get3DProperties();
-
-        if (!deckProps.hasDeckDoubleGridProperty("SWATINIT"))
-            return; // SWATINIT is not in the deck
-        else if (!deck.hasKeyword("EQUIL"))
-            // SWATINIT only applies if the initial solution is specified using the EQUIL
-            // keyword.
-            return;
-
-        // SWATINIT applies. We have to do a complete re-initialization here. this is a
-        // kludge, but to calculate the threshold pressures the initial solution without
-        // SWATINIT is required!
-
-        typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
-        EquilInitializer equilInitializer(this->simulator(),
-                                          *materialLawManager_,
-                                          /*enableSwatinit=*/true);
-        auto& model = this->model();
-        size_t numElems = model.numGridDof();
-        for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
-            const auto& fs = equilInitializer.initialFluidState(elemIdx);
-            auto& priVars0 = model.solution(/*timeIdx=*/0)[elemIdx];
-            auto& priVars1 = model.solution(/*timeIdx=*/1)[elemIdx];
-
-            priVars1.assignNaive(fs);
-            priVars0 = priVars1;
-        }
-
-        // the solution (most likely) has changed, so we need to invalidate the cache for
-        // intensive quantities
-        this->simulator().model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
     }
 
 private:
@@ -1275,20 +1217,9 @@ private:
 
     void readEquilInitialCondition_()
     {
-        const auto& deck = this->simulator().gridManager().deck();
+        // initial condition corresponds to hydrostatic conditions.
-        const auto& eclState = this->simulator().gridManager().eclState();
-        int numEquilRegions = eclState.getTableManager().getEqldims().getNumEquilRegions();
-        bool useThpres = deck.hasKeyword("THPRES") && numEquilRegions > 1;
-        bool useSwatinit =
-            GET_PROP_VALUE(TypeTag, EnableSwatinit) &&
-            eclState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT");
-
-        // initial condition corresponds to hydrostatic conditions. The SWATINIT keyword
-        // can be considered here directly if threshold pressures are disabled.
         typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
-        EquilInitializer equilInitializer(this->simulator(),
+        EquilInitializer equilInitializer(this->simulator(), *materialLawManager_);
-                                          *materialLawManager_,
-                                          /*enableSwatinit=*/useThpres && useSwatinit);
 
         // since the EquilInitializer provides fluid states that are consistent with the
         // black-oil model, we can use naive instead of mass conservative determination