WIP Add flow_dyn

CMakeLists.txt

@@ -22,6 +22,7 @@ option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON)
 set( USE_OPENMP_DEFAULT OFF ) # Use of OpenMP is considered experimental
 option(BUILD_FLOW "Build the production oriented flow simulator?" ON)
 option(BUILD_FLOW_BLACKOIL_ONLY "Build the production oriented flow simulator only supporting the blackoil model?" OFF)
+option(BUILD_FLOW_DYN "Build the flow simulator with dynamic indicies?" OFF)
 option(BUILD_FLOW_VARIANTS "Build the variants for flow by default?" OFF)
 option(BUILD_EBOS "Build the research oriented ebos simulator?" ON)
 option(BUILD_EBOS_EXTENSIONS "Build the variants for various extensions of ebos by default?" OFF)
@@ -293,6 +294,29 @@ opm_add_test(flow_blackoil_dunecpr
   DEPENDS opmsimulators
   LIBRARIES opmsimulators)
 
+if (NOT BUILD_FLOW_DYN)
+  set(FLOW_DYN_DEFAULT_ENABLE_IF "FALSE")
+else()
+  set(FLOW_DYN_DEFAULT_ENABLE_IF "TRUE")
+endif()
+
+# the production oriented general-purpose ECL simulator
+opm_add_test(flow_dyn
+  ONLY_COMPILE
+  ALWAYS_ENABLE
+  DEFAULT_ENABLE_IF ${FLOW_DYN_DEFAULT_ENABLE_IF}
+  DEPENDS opmsimulators
+  LIBRARIES opmsimulators
+  SOURCES
+  flow/flow.cpp
+  flow/flow_1.cpp
+  flow/flow_2.cpp
+  flow/flow_3.cpp
+  flow/flow_4.cpp
+  flow/flow_ebos_blackoil.cpp
+  $<TARGET_OBJECTS:moduleVersion>)
+target_compile_definitions(flow_dyn PRIVATE "FLOW_DYN")
+
 opm_add_test(flow_onephase
   ONLY_COMPILE
   DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}

CMakeLists_files.cmake

@@ -238,6 +238,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/MSWellHelpers.hpp
   opm/simulators/wells/BlackoilWellModel.hpp
   opm/simulators/wells/BlackoilWellModel_impl.hpp
+  opm/simulators/wells/WellIndices.hpp
   )
 
 list (APPEND EXAMPLE_SOURCE_FILES

ebos/eclequilinitializer.hh

@@ -82,7 +82,7 @@ public:
                                     FluidSystem,
                                     enableTemperature,
                                     enableEnergy,
-                                    Indices::gasEnabled,
+                                    true,
                                     enableBrine,
                                     Indices::numPhases
                                     > ScalarFluidState;
@@ -116,18 +116,16 @@ public:
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                 if (FluidSystem::phaseIsActive(phaseIdx))
                     fluidState.setSaturation(phaseIdx, initialState.saturation()[phaseIdx][elemIdx]);
-                else if (Indices::numPhases == 3)
-                    fluidState.setSaturation(phaseIdx, 0.0);
             }
 
             if (FluidSystem::enableDissolvedGas())
                 fluidState.setRs(initialState.rs()[elemIdx]);
-            else if (Indices::gasEnabled)
+            else if (Indices::gasIsActive() && Indices::oilIsActive())
                 fluidState.setRs(0.0);
 
             if (FluidSystem::enableVaporizedOil())
                 fluidState.setRv(initialState.rv()[elemIdx]);
-            else if (Indices::gasEnabled)
+            else if (Indices::gasIsActive() && Indices::oilIsActive())
                 fluidState.setRv(0.0);
 
 

ebos/eclfluxmodule.hh

@@ -106,11 +106,11 @@ class EclTransExtensiveQuantities
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
     using GridView = GetPropType<TypeTag, Properties::GridView>;
     using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
 
     enum { dimWorld = GridView::dimensionworld };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
     enum { numPhases = FluidSystem::numPhases };
-    enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
 
     typedef Opm::MathToolbox<Evaluation> Toolbox;
@@ -438,7 +438,7 @@ protected:
                 volumeFlux_[phaseIdx] =
                     pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea);
 
-                if (enableSolvent && phaseIdx == gasPhaseIdx)
+                if (Indices::solventIsActive() && phaseIdx == gasPhaseIdx)
                     asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-transModified/faceArea));
             }
             else {
@@ -456,7 +456,7 @@ protected:
                     pressureDifference_[phaseIdx]*mob*(-transModified/faceArea);
 
                 // Solvent inflow is not yet supported
-                if (enableSolvent && phaseIdx == gasPhaseIdx)
+                if (Indices::solventIsActive() && phaseIdx == gasPhaseIdx)
                     asImp_().setSolventVolumeFlux(0.0);
             }
         }

ebos/ecloutputblackoilmodule.hh

@@ -94,6 +94,7 @@ class EclOutputBlackOilModule
     using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementIterator = typename GridView::template Codim<0>::Iterator;
 
@@ -388,7 +389,7 @@ public:
             rstKeywords["RV"] = 0;
         }
 
-        if (getPropValue<TypeTag, Properties::EnableSolvent>())
+        if (Indices::solventIsActive())
             sSol_.resize(bufferSize, 0.0);
         if (getPropValue<TypeTag, Properties::EnablePolymer>())
             cPolymer_.resize(bufferSize, 0.0);

ebos/eclproblem.hh

@@ -876,6 +876,8 @@ public:
 
         initFluidSystem_();
 
+        initBlackoilIndices_();
+
         // deal with DRSDT
         unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables();
         size_t numDof = this->model().numGridDof();
@@ -1775,7 +1777,7 @@ public:
         values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
         values.assignNaive(initialFluidStates_[globalDofIdx]);
 
-        if (enableSolvent)
+        if (Indices::solventIsActive())
             values[Indices::solventSaturationIdx] = solventSaturation_[globalDofIdx];
 
         if (enablePolymer)
@@ -2192,11 +2194,11 @@ private:
                                          std::to_string(Indices::numPhases)+".");
 
             // make sure that the correct phases are active
-            if (FluidSystem::phaseIsActive(oilPhaseIdx) && !Indices::oilEnabled)
+            if (FluidSystem::phaseIsActive(oilPhaseIdx) && !Indices::oilIsActive())
                 throw std::runtime_error("The deck enables oil, but this simulator cannot handle it.");
-            if (FluidSystem::phaseIsActive(gasPhaseIdx) && !Indices::gasEnabled)
+            if (FluidSystem::phaseIsActive(gasPhaseIdx) && !Indices::gasIsActive())
                 throw std::runtime_error("The deck enables gas, but this simulator cannot handle it.");
-            if (FluidSystem::phaseIsActive(waterPhaseIdx) && !Indices::waterEnabled)
+            if (FluidSystem::phaseIsActive(waterPhaseIdx) && !Indices::waterIsActive())
                 throw std::runtime_error("The deck enables water, but this simulator cannot handle it.");
             // the opposite cases should be fine (albeit a bit slower than what's possible)
         }
@@ -2630,6 +2632,13 @@ private:
         FluidSystem::initFromState(eclState, schedule);
    }
 
+   void initBlackoilIndices_()
+   {
+        const auto& simulator = this->simulator();
+        const auto& eclState = simulator.vanguard().eclState();
+        Indices::init(eclState.runspec().phases());
+   }
+
     void readInitialCondition_()
     {
         const auto& simulator = this->simulator();
@@ -2698,7 +2707,7 @@ private:
 
         size_t numElems = this->model().numGridDof();
         initialFluidStates_.resize(numElems);
-        if (enableSolvent)
+        if (Indices::solventIsActive())
             solventSaturation_.resize(numElems, 0.0);
 
         if (enablePolymer)
@@ -2725,13 +2734,13 @@ private:
             // the function and discard the unchanged result.  Do not index
             // into 'solventSaturation_' unless solvent is enabled.
             {
-                auto ssol = enableSolvent
+                auto ssol = Indices::solventIsActive()
                     ? eclWriter_->eclOutputModule().getSolventSaturation(elemIdx)
                     : Scalar(0);
 
                 processRestartSaturations_(elemFluidState, ssol);
 
-                if (enableSolvent)
+                if (Indices::solventIsActive())
                     solventSaturation_[elemIdx] = ssol;
             }
 
@@ -2803,7 +2812,7 @@ private:
             }
 
         }
-        if (enableSolvent) {
+        if (Indices::solventIsActive()) {
             if (solventSaturation < smallSaturationTolerance)
                 solventSaturation = 0.0;
 
@@ -2818,7 +2827,7 @@ private:
                 elemFluidState.setSaturation(phaseIdx, saturation);
             }
         }
-        if (enableSolvent) {
+        if (Indices::solventIsActive()) {
             solventSaturation = solventSaturation / sumSaturation;
         }
     }
@@ -2934,12 +2943,12 @@ private:
 
             if (FluidSystem::enableDissolvedGas())
                 dofFluidState.setRs(rsData[dofIdx]);
-            else if (Indices::gasEnabled && Indices::oilEnabled)
+            else if (Indices::gasIsActive() && Indices::oilIsActive())
                 dofFluidState.setRs(0.0);
 
             if (FluidSystem::enableVaporizedOil())
                 dofFluidState.setRv(rvData[dofIdx]);
-            else if (Indices::gasEnabled && Indices::oilEnabled)
+            else if (Indices::gasIsActive() && Indices::oilIsActive())
                 dofFluidState.setRv(0.0);
 
             //////
@@ -2966,7 +2975,7 @@ private:
         const auto& eclState = vanguard.eclState();
         size_t numDof = this->model().numGridDof();
 
-        if (enableSolvent) {
+        if (Indices::solventIsActive()) {
             if (eclState.fieldProps().has_double("SSOL"))
                 solventSaturation_ = eclState.fieldProps().get_double("SSOL");
             else
@@ -3149,7 +3158,7 @@ private:
                         compIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
                         break;
                     case BCComponent::SOLVENT:
-                        if (!enableSolvent)
+                        if (!Indices::solventIsActive())
                             throw std::logic_error("solvent is disabled and you're trying to add solvent to BC");
 
                         compIdx = Indices::solventSaturationIdx;

ebos/eclwriter.hh

@@ -166,6 +166,7 @@ class EclWriter
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementIterator = typename GridView::template Codim<0>::Iterator;
 
@@ -174,7 +175,6 @@ class EclWriter
     typedef std::vector<Scalar> ScalarBuffer;
 
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
-    enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
 
 
 public:
@@ -396,7 +396,7 @@ public:
             {"SWAT", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))},
             {"SGAS", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))},
             {"TEMP" , Opm::UnitSystem::measure::temperature, enableEnergy},
-            {"SSOLVENT" , Opm::UnitSystem::measure::identity, enableSolvent},
+            {"SSOLVENT" , Opm::UnitSystem::measure::identity, Indices::solventIsActive()},
             {"RS", Opm::UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
             {"RV", Opm::UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
             {"SOMAX", Opm::UnitSystem::measure::identity, simulator_.problem().vapparsActive()},

flow/flow_1.cpp

@@ -0,0 +1,92 @@
+/*
+  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 3 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/>.
+*/
+#include "config.h"
+
+#include <flow/flow_1.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/models/blackoil/blackoildynindices.hh>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/simulators/flow/FlowMainEbos.hpp>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/misc/mpimanager.hh>
+#else
+#include <dune/common/parallel/mpihelper.hh>
+#endif
+
+namespace Opm {
+namespace Properties {
+namespace TTag {
+struct EclFlow1Problem {
+    using InheritsFrom = std::tuple<EclFlowProblem>;
+};
+}
+
+//! The indices required by the model
+template<class TypeTag>
+struct Indices<TypeTag, TTag::EclFlow1Problem>
+{
+private:
+    // it is unfortunately not possible to simply use 'TypeTag' here because this leads
+    // to cyclic definitions of some properties. if this happens the compiler error
+    // messages unfortunately are *really* confusing and not really helpful.
+    using BaseTypeTag = TTag::EclFlowProblem;
+    using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
+
+public:
+    typedef Opm::BlackOilDynIndices<1,1> type;
+};
+}}
+
+namespace Opm {
+void flow1SetDeck(double setupTime, std::unique_ptr<Deck> deck,
+                             std::unique_ptr<EclipseState> eclState,
+                             std::unique_ptr<Schedule> schedule,
+                             std::unique_ptr<SummaryConfig> summaryConfig)
+{
+    using TypeTag = Properties::TTag::EclFlow1Problem;
+    using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
+
+    Vanguard::setExternalSetupTime(setupTime);
+    Vanguard::setExternalDeck(std::move(deck));
+    Vanguard::setExternalEclState(std::move(eclState));
+    Vanguard::setExternalSchedule(std::move(schedule));
+    Vanguard::setExternalSummaryConfig(std::move(summaryConfig));
+}
+
+// ----------------- Main program -----------------
+int flow1Main(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    Opm::resetLocale();
+
+#if HAVE_DUNE_FEM
+    Dune::Fem::MPIManager::initialize(argc, argv);
+#else
+    Dune::MPIHelper::instance(argc, argv);
+#endif
+
+    Opm::FlowMainEbos<Properties::TTag::EclFlow1Problem>
+        mainfunc {argc, argv, outputCout, outputFiles};
+    return mainfunc.execute();
+}
+
+}

flow/flow_1.hpp

@@ -0,0 +1,33 @@
+/*
+  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 3 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/>.
+*/
+#ifndef FLOW_1_HPP
+#define FLOW_1_HPP
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
+
+namespace Opm {
+void flow1SetDeck(double setupTime, std::unique_ptr<Deck> deck,
+                             std::unique_ptr<EclipseState> eclState,
+                             std::unique_ptr<Schedule> schedule,
+                             std::unique_ptr<SummaryConfig> summaryConfig);
+int flow1Main(int argc, char** argv, bool outputCout, bool outputFiles);
+}
+
+#endif // FLOW_1_HPP