Merge remote-tracking branch 'upstream/master' into PR/ewoms-can-use-dune-fem

applications/ebos/eclequilinitializer.hh

@@ -0,0 +1,208 @@
+/*
+  Copyright (C) 2014 by Andreas Lauser
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+/**
+ * \file
+ *
+ * \copydoc Ewoms::EclEquilInitializer
+ */
+#ifndef EWOMS_ECL_EQUIL_INITIALIZER_HH
+#define EWOMS_ECL_EQUIL_INITIALIZER_HH
+
+#include <opm/material/fluidstates/CompositionalFluidState.hpp>
+
+// the ordering of these includes matters. do not touch it if you're not prepared to deal
+// with some trouble!
+#include <dune/grid/cpgrid/GridHelpers.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/simulator/initStateEquil.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+
+#include <vector>
+
+namespace Ewoms {
+/*!
+ * \ingroup EclBlackOilSimulator
+ *
+ * \brief Computes the initial condition based on the EQUIL keyword from ECL.
+ *
+ * So far, it uses the "initStateEquil()" function from opm-core. Since this method is
+ * very much glued into the opm-core data structures, it should be reimplemented in the
+ * medium to long term for some significant memory savings and less significant
+ * performance improvements.
+ */
+template <class TypeTag>
+class EclEquilInitializer
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+
+    typedef Opm::CompositionalFluidState<Scalar, FluidSystem> ScalarFluidState;
+
+    enum { numPhases = FluidSystem::numPhases };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+
+    enum { numComponents = FluidSystem::numComponents };
+    enum { oilCompIdx = FluidSystem::oilCompIdx };
+    enum { gasCompIdx = FluidSystem::gasCompIdx };
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+
+    enum { dimWorld = GridView::dimensionworld };
+
+public:
+    template <class MaterialLawManager>
+    EclEquilInitializer(const Simulator& simulator,
+                        std::shared_ptr<MaterialLawManager> materialLawManager)
+        : simulator_(simulator)
+    {
+        const auto& gridManager = simulator.gridManager();
+        const auto& grid = gridManager.grid();
+
+        // create the data structures which are used by initStateEquil()
+        Opm::parameter::ParameterGroup tmpParam;
+        Opm::BlackoilPropertiesFromDeck opmBlackoilProps(
+            gridManager.deck(),
+            gridManager.eclState(),
+            materialLawManager,
+            Opm::UgGridHelpers::numCells(grid),
+            Opm::UgGridHelpers::globalCell(grid),
+            Opm::UgGridHelpers::cartDims(grid),
+            Opm::UgGridHelpers::beginCellCentroids(grid),
+            Opm::UgGridHelpers::dimensions(grid),
+            tmpParam);
+
+        // initialize the boiler plate of opm-core the state structure.
+        Opm::BlackoilState opmBlackoilState;
+        opmBlackoilState.init(grid.size(/*codim=*/0),
+                              /*numFaces=*/0, // we don't care here
+                              numPhases);
+
+        // do the actual computation.
+        Opm::initStateEquil(gridManager.grid(),
+                            opmBlackoilProps,
+                            gridManager.deck(),
+                            gridManager.eclState(),
+                            simulator.problem().gravity()[dimWorld - 1],
+                            opmBlackoilState);
+
+        const Scalar rhooRef = FluidSystem::referenceDensity(oilPhaseIdx, /*regionIdx=*/0);
+        const Scalar rhogRef = FluidSystem::referenceDensity(gasPhaseIdx, /*regionIdx=*/0);
+        const Scalar MG = FluidSystem::molarMass(gasCompIdx);
+        const Scalar MO = FluidSystem::molarMass(oilCompIdx);
+
+        // copy the result into the array of initial fluid states
+        int numElems = gridManager.gridView().size(/*codim=*/0);
+        initialFluidStates_.resize(numElems);
+        for (int elemIdx = 0; elemIdx < numElems; ++elemIdx) {
+            auto &fluidState = initialFluidStates_[elemIdx];
+
+            // get the PVT region index of the current element
+            int regionIdx = simulator_.problem().pvtRegionIndex(elemIdx);
+
+            // set the phase saturations
+            for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                Scalar S = opmBlackoilState.saturation()[elemIdx*numPhases + phaseIdx];
+                fluidState.setSaturation(phaseIdx, S);
+            }
+
+            // set the temperature
+            const auto& temperatureVector = opmBlackoilState.temperature();
+            Scalar T = FluidSystem::surfaceTemperature;
+            if (!temperatureVector.empty())
+                T = temperatureVector[elemIdx];
+            fluidState.setTemperature(T);
+
+            // set the phase pressures. the Opm::BlackoilState only provides the oil
+            // phase pressure, so we need to calculate the other phases' pressures
+            // ourselfs.
+            Scalar pC[numPhases];
+            const auto& matParams = simulator.problem().materialLawParams(elemIdx);
+            MaterialLaw::capillaryPressures(pC, matParams, fluidState);
+            Scalar po = opmBlackoilState.pressure()[elemIdx];
+            for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                fluidState.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
+            Scalar pg = fluidState.pressure(gasPhaseIdx);
+
+            // reset the phase compositions
+            for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                    fluidState.setMoleFraction(phaseIdx, compIdx, 0.0);
+
+            // the composition of the water phase is simple: it only consists of the
+            // water component.
+            fluidState.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
+
+            if (gridManager.deck()->hasKeyword("DISGAS")) {
+                // 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 = opmBlackoilState.gasoilratio()[elemIdx];
+
+                // dissolved gas surface volume to mass fraction
+                Scalar XoG = Rs/(rhooRef/rhogRef + Rs);
+                // mass fraction to mole fraction
+                Scalar xoG = XoG*MO / (MG*(1 - XoG) + XoG*MO);
+
+                Scalar xoGMax = FluidSystem::saturatedOilGasMoleFraction(T, pg, regionIdx);
+                if (fluidState.saturation(gasPhaseIdx) > 0.0 || xoG > xoGMax)
+                    xoG = xoGMax;
+
+                fluidState.setMoleFraction(oilPhaseIdx, oilCompIdx, 1 - xoG);
+                fluidState.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
+            }
+
+            // retrieve the surface volume of vaporized gas
+            if (gridManager.deck()->hasKeyword("VAPOIL")) {
+                Scalar Rv = opmBlackoilState.rv()[elemIdx];
+
+                // vaporized oil surface volume to mass fraction
+                Scalar XgO = Rv/(rhogRef/rhooRef + Rv);
+                // mass fraction to mole fraction
+                Scalar xgO = XgO*MG / (MO*(1 - XgO) + XgO*MG);
+
+                Scalar xgOMax = FluidSystem::saturatedGasOilMoleFraction(T, pg, regionIdx);
+                if (fluidState.saturation(oilPhaseIdx) > 0.0 || xgO > xgOMax)
+                    xgO = xgOMax;
+
+                fluidState.setMoleFraction(gasPhaseIdx, oilCompIdx, xgO);
+                fluidState.setMoleFraction(gasPhaseIdx, gasCompIdx, 1 - xgO);
+            }
+        }
+    }
+
+    /*!
+     * \brief Return the initial thermodynamic state which should be used as the initial
+     *        condition.
+     *
+     * This is supposed to correspond to hydrostatic conditions.
+     */
+    const ScalarFluidState& initialFluidState(int elemIdx) const
+    { return initialFluidStates_[elemIdx]; }
+
+protected:
+    const Simulator& simulator_;
+
+    std::vector<ScalarFluidState> initialFluidStates_;
+};
+} // namespace Ewoms
+
+#endif

applications/ebos/eclproblem.hh

@@ -30,6 +30,7 @@
 
 #include "eclgridmanager.hh"
 #include "eclwellmanager.hh"
+#include "eclequilinitializer.hh"
 #include "eclwriter.hh"
 #include "eclsummarywriter.hh"
 #include "ecloutputblackoilmodule.hh"
@@ -343,9 +344,17 @@ public:
             ++ nextEpisodeIdx;
         }
 
+        Scalar episodeLength = timeMap->getTimeStepLength(nextEpisodeIdx);
+        Scalar dt = episodeLength;
+        if (nextEpisodeIdx == 0) {
+            // allow the size of the initial time step to be set via an external parameter
+            Scalar initialDt = EWOMS_GET_PARAM(TypeTag, Scalar, InitialTimeStepSize);
+            dt = std::min(dt, initialDt);
+        }
+
         if (nextEpisodeIdx < numReportSteps) {
-            simulator.startNextEpisode(timeMap->getTimeStepLength(nextEpisodeIdx));
+            simulator.startNextEpisode(episodeLength);
-            simulator.setTimeStepSize(timeMap->getTimeStepLength(nextEpisodeIdx));
+            simulator.setTimeStepSize(dt);
         }
 
         // set up the wells
@@ -555,14 +564,23 @@ public:
                                                int spaceIdx, int timeIdx) const
     {
         int globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
-        return materialLawParams_(globalSpaceIdx);
+        return materialLawParams(globalSpaceIdx);
     }
 
+    const MaterialLawParams& materialLawParams(int globalDofIdx) const
+    { return materialLawManager_->materialLawParams(globalDofIdx); }
+
     /*!
      * \brief Returns the index of the relevant region for thermodynmic properties
      */
     template <class Context>
     int pvtRegionIndex(const Context &context, int spaceIdx, int timeIdx) const
+    { return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
+
+    /*!
+     * \brief Returns the index the relevant PVT region given a cell index
+     */
+    int pvtRegionIndex(int elemIdx) const
     {
         Opm::DeckConstPtr deck = this->simulator().gridManager().deck();
 
@@ -571,10 +589,7 @@ public:
 
         const auto& gridManager = this->simulator().gridManager();
 
-        // this is quite specific to the ECFV discretization. But so is everything in an
+        int cartesianDofIdx = gridManager.cartesianCellId(elemIdx);
-        // ECL deck, i.e., we don't need to care here...
-        int compressedDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
-        int cartesianDofIdx = gridManager.cartesianCellId(compressedDofIdx);
 
         return deck->getKeyword("PVTNUM")->getIntData()[cartesianDofIdx] - 1;
     }
@@ -622,8 +637,12 @@ public:
 
         values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
 
-        const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
+        if (useMassConservativeInitialCondition_) {
-        values.assignMassConservative(initialFluidStates_[globalDofIdx], matParams);
+            const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
+            values.assignMassConservative(initialFluidStates_[globalDofIdx], matParams);
+        }
+        else
+            values.assignNaive(initialFluidStates_[globalDofIdx]);
     }
 
     void initialSolutionApplied()
@@ -794,7 +813,8 @@ private:
         for (unsigned elemIdx = 0; elemIdx < numDof; ++elemIdx)
             compressedToCartesianElemIdx[elemIdx] = gridManager.cartesianCellId(elemIdx);
 
-        materialLawManager_.initFromDeck(deck, eclState, compressedToCartesianElemIdx);
+        materialLawManager_ = std::make_shared<EclMaterialLawManager>();
+        materialLawManager_->initFromDeck(deck, eclState, compressedToCartesianElemIdx);
     }
 
     void initFluidSystem_()
@@ -836,6 +856,7 @@ private:
     Opm::OilPvtInterface<Scalar, Evaluation>* createOilPvt_(Opm::DeckConstPtr deck,
                                                             Opm::EclipseStateConstPtr eclState)
     {
+        const auto tableManager = eclState->getTableManager();
         Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
         int numPvtRegions = densityKeyword->size();
 
@@ -844,7 +865,7 @@ private:
             oilPvt->setNumRegions(numPvtRegions);
 
             for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx)
-                oilPvt->setPvtoTable(regionIdx, eclState->getPvtoTables()[regionIdx]);
+                oilPvt->setPvtoTable(regionIdx, tableManager->getPvtoTables()[regionIdx]);
 
             oilPvt->initEnd();
             return oilPvt;
@@ -854,7 +875,7 @@ private:
             oilPvt->setNumRegions(numPvtRegions);
 
             for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx)
-                oilPvt->setPvdoTable(regionIdx, eclState->getPvdoTables()[regionIdx]);
+                oilPvt->setPvdoTable(regionIdx, tableManager->getPvdoTables()[regionIdx]);
 
             oilPvt->initEnd();
             return oilPvt;
@@ -880,6 +901,7 @@ private:
                                                             Opm::EclipseStateConstPtr eclState)
     {
         Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
+        const auto tableManager = eclState->getTableManager();
         int numPvtRegions = densityKeyword->size();
 
         if (deck->hasKeyword("PVTG")) {
@@ -887,7 +909,7 @@ private:
             gasPvt->setNumRegions(numPvtRegions);
 
             for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx)
-                gasPvt->setPvtgTable(regionIdx, eclState->getPvtgTables()[regionIdx]);
+                gasPvt->setPvtgTable(regionIdx, tableManager->getPvtgTables()[regionIdx]);
 
             gasPvt->initEnd();
             return gasPvt;
@@ -897,7 +919,7 @@ private:
             gasPvt->setNumRegions(numPvtRegions);
 
             for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx)
-                gasPvt->setPvdgTable(regionIdx, eclState->getPvdgTables()[regionIdx]);
+                gasPvt->setPvdgTable(regionIdx, tableManager->getPvdgTables()[regionIdx]);
 
             gasPvt->initEnd();
             return gasPvt;
@@ -927,30 +949,67 @@ private:
     }
 
     void readInitialCondition_()
+    {
+        const auto &gridManager = this->simulator().gridManager();
+        const auto deck = gridManager.deck();
+
+        if (!deck->hasKeyword("EQUIL"))
+            readExplicitInitialCondition_();
+        else
+            readEquilInitialCondition_();
+    }
+
+    void readEquilInitialCondition_()
+    {
+        // The EQUIL initializer also modifies the material law manager according to
+        // SWATINIT (although it does not belong there strictly speaking)
+        typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
+        EquilInitializer equilInitializer(this->simulator(), materialLawManager_);
+
+        // since the EquilInitializer provides fluid states that are consistent with the
+        // black-oil model, we can use naive instead of mass conservative determination
+        // of the primary variables.
+        useMassConservativeInitialCondition_ = false;
+
+        size_t numElems = this->model().numGridDof();
+        initialFluidStates_.resize(numElems);
+        for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
+            auto &elemFluidState = initialFluidStates_[elemIdx];
+            elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
+        }
+    }
+
+    void readExplicitInitialCondition_()
     {
         const auto &gridManager = this->simulator().gridManager();
         const auto deck = gridManager.deck();
         const auto eclState = gridManager.eclState();
 
+        // since the values specified in the deck do not need to be consistent, we use an
+        // initial condition that conserves the total mass specified by these values.
+        useMassConservativeInitialCondition_ = true;
+
         bool enableDisgas = deck->hasKeyword("DISGAS");
         bool enableVapoil = deck->hasKeyword("VAPOIL");
 
         // make sure all required quantities are enables
-        if (!deck->hasKeyword("SWAT") ||
+         if (!deck->hasKeyword("SWAT") ||
-            !deck->hasKeyword("SGAS"))
+             !deck->hasKeyword("SGAS"))
-            OPM_THROW(std::runtime_error,
+             OPM_THROW(std::runtime_error,
-                      "So far, the ECL input file requires the presence of the SWAT "
+                      "The ECL input file requires the presence of the SWAT "
-                      "and SGAS keywords");
+                      "and SGAS keywords if the model is initialized explicitly");
-        if (!deck->hasKeyword("PRESSURE"))
+         if (!deck->hasKeyword("PRESSURE"))
-            OPM_THROW(std::runtime_error,
+             OPM_THROW(std::runtime_error,
-                      "So far, the ECL input file requires the presence of the PRESSURE "
+                      "The ECL input file requires the presence of the PRESSURE "
-                      "keyword");
+                      "keyword if the model is initialized explicitly");
-        if (enableDisgas && !deck->hasKeyword("RS"))
+         if (enableDisgas && !deck->hasKeyword("RS"))
-            OPM_THROW(std::runtime_error,
+             OPM_THROW(std::runtime_error,
-                      "The ECL input file requires the RS keyword to be present if dissolved gas is enabled");
+                      "The ECL input file requires the RS keyword to be present if"
-        if (enableVapoil && !deck->hasKeyword("RV"))
+                      " dissolved gas is enabled");
-            OPM_THROW(std::runtime_error,
+         if (enableVapoil && !deck->hasKeyword("RV"))
-                      "The ECL input file requires the RV keyword to be present if vaporized oil is enabled");
+             OPM_THROW(std::runtime_error,
+                      "The ECL input file requires the RV keyword to be present if"
+                      " vaporized oil is enabled");
 
         size_t numDof = this->model().numGridDof();
 
@@ -1020,9 +1079,8 @@ private:
 
             // this assumes that capillary pressures only depend on the phase saturations
             // and possibly on temperature. (this is always the case for ECL problems.)
-            Scalar pc[numPhases] = { 0 };
+            Scalar pc[numPhases];
-
+            const auto& matParams = materialLawParams(dofIdx);
-            const auto& matParams = materialLawParams_(dofIdx);
             MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
             Valgrind::CheckDefined(oilPressure);
             Valgrind::CheckDefined(pc);
@@ -1130,15 +1188,10 @@ private:
         }
     }
 
-    const MaterialLawParams& materialLawParams_(int globalDofIdx) const
-    {
-        return materialLawManager_.materialLawParams(globalDofIdx);
-    }
-
     // update the hysteresis parameters of the material laws for the whole grid
     void updateHysteresis_()
     {
-        if (!materialLawManager_.enableHysteresis())
+        if (!materialLawManager_->enableHysteresis())
             return;
 
         ElementContext elemCtx(this->simulator());
@@ -1154,9 +1207,8 @@ private:
             elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
 
             int compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
-            int cartesianDofIdx = gridManager.cartesianCellId(compressedDofIdx);
             const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
-            materialLawManager_.updateHysteresis(intQuants.fluidState(), cartesianDofIdx);
+            materialLawManager_->updateHysteresis(intQuants.fluidState(), compressedDofIdx);
         }
     }
 
@@ -1164,11 +1216,12 @@ private:
     std::vector<DimMatrix> intrinsicPermeability_;
     EclTransmissibility<TypeTag> transmissibilities_;
 
-    EclMaterialLawManager materialLawManager_;
+    std::shared_ptr<EclMaterialLawManager> materialLawManager_;
 
     std::vector<unsigned short> rockTableIdx_;
     std::vector<RockParams> rockParams_;
 
+    bool useMassConservativeInitialCondition_;
     std::vector<ScalarFluidState> initialFluidStates_;
 
     EclWellManager<TypeTag> wellManager_;