Merge pull request #2851 from osae/extboSimulators

Alternative solvent extension for the black oil model.

CMakeLists.txt

@@ -228,7 +228,7 @@ add_dependencies(moduleVersion opmsimulators)
 set(FLOW_TGTS)
 foreach(OBJ blackoil brine energy foam gasoil oilwater
             oilwater_brine oilwater_polymer
-            oilwater_polymer_injectivity polymer solvent)
+            oilwater_polymer_injectivity polymer solvent extbo)
   add_library(flow_lib${OBJ} OBJECT flow/flow_ebos_${OBJ}.cpp)
   list(APPEND FLOW_TGTS $<TARGET_OBJECTS:flow_lib${OBJ}>)
 endforeach()
@@ -342,7 +342,7 @@ endif()
 # the research oriented general-purpose ECL simulator ("ebos" == &ecl
 # &black-&oil &simulator)
 set(MEBOS_TARGETS "")
-foreach(OBJ blackoil solvent polymer foam brine gasoil oilwater oilwaterpolymer thermal)
+foreach(OBJ blackoil solvent extbo polymer foam brine gasoil oilwater oilwaterpolymer thermal)
   add_library(ebos_lib${OBJ} OBJECT EXCLUDE_FROM_ALL ebos/ebos_${OBJ}.cc)
   list(APPEND MEBOS_TARGETS $<TARGET_OBJECTS:ebos_lib${OBJ}>)
 endforeach()
@@ -366,7 +366,7 @@ else()
   set(EBOS_EXTENSIONS_DEFAULT_ENABLE_IF "TRUE")
 endif()
 
-foreach(OBJ solvent polymer foam brine gasoil oilwater oilwaterpolymer thermal)
+foreach(OBJ solvent extbo polymer foam brine gasoil oilwater oilwaterpolymer thermal)
   opm_add_test(ebos_${OBJ}
     ONLY_COMPILE
     DEFAULT_ENABLE_IF ${EBOS_EXTENSIONS_DEFAULT_ENABLE_IF}
@@ -406,7 +406,7 @@ opm_add_test(ebos_plain
   LIBRARIES opmsimulators)
 
 if (BUILD_EBOS_EXTENSIONS)
-  foreach(TGT ebos_solvent ebos_polymer ebos_foam ebos_brine ebos_gasoil ebos_oilwater ebos_oilwaterpolymer ebos_thermal mebos)
+  foreach(TGT ebos_solvent ebos_extbo ebos_polymer ebos_foam ebos_brine ebos_gasoil ebos_oilwater ebos_oilwaterpolymer ebos_thermal mebos)
     install(TARGETS ${TGT} DESTINATION bin)
     opm_add_bash_completion(${TGT})
   endforeach()

ebos/ebos_extbo.cc

@@ -0,0 +1,71 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief A general-purpose simulator for ECL decks using the black-oil model.
+ */
+#include "config.h"
+
+#include "ebos.hh"
+#include "startEbos.hh"
+
+namespace Opm::Properties {
+
+namespace TTag {
+struct EbosExtboTypeTag {
+    using InheritsFrom = std::tuple<EbosTypeTag>;
+};
+}
+
+// enable the polymer extension of the black oil model
+template<class TypeTag>
+struct EnableExtbo<TypeTag, TTag::EbosExtboTypeTag> {
+    static constexpr bool value = true;
+};
+
+} // namespace Opm::Properties
+
+namespace Opm {
+
+void ebosExtboSetDeck(std::unique_ptr<Opm::Deck> deck,
+                        std::unique_ptr<Opm::ParseContext> parseContext,
+                        std::unique_ptr<Opm::ErrorGuard> errorGuard,
+                        double externalSetupTime)
+{
+    using ProblemTypeTag = Properties::TTag::EbosExtboTypeTag;
+    using Vanguard = GetPropType<ProblemTypeTag, Properties::Vanguard>;
+
+    Vanguard::setExternalSetupTime(externalSetupTime);
+    Vanguard::setExternalParseContext(std::move(parseContext));
+    Vanguard::setExternalErrorGuard(std::move(errorGuard));
+    Vanguard::setExternalDeck(std::move(deck));
+}
+
+int ebosExtboMain(int argc, char **argv)
+{
+    using ProblemTypeTag = Properties::TTag::EbosExtboTypeTag;
+    return Opm::startEbos<ProblemTypeTag>(argc, argv);
+}
+
+}

ebos/ebos_extbo.hh

@@ -0,0 +1,44 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief The function prototypes required to start the solvent variant of ebos
+ */
+#ifndef EBOS_EXTBO_HH
+#define EBOS_EXTBO_HH
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/Parser/ErrorGuard.hpp>
+
+namespace Opm {
+void ebosExtboSetDeck(Opm::Deck* deck,
+                        Opm::ParseContext* parseContext,
+                        Opm::ErrorGuard* errorGuard,
+                        double externalSetupTime);
+
+int ebosExtboMain(int argc, char** argv);
+}
+
+#endif

ebos/ebos_extbo_main.cc

@@ -0,0 +1,37 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief The main function for the stand alone solvent variant of ebos.
+ *
+ * This only calls the ebosExtboMain() function.
+ */
+#include "config.h"
+
+#include "ebos_extbo.hh"
+
+int main(int argc, char** argv)
+{
+    return Opm::ebosExtboMain(argc, argv);
+}

ebos/ebos_gasoil.cc

@@ -50,6 +50,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

ebos/ebos_oilwater.cc

@@ -50,6 +50,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

ebos/ebos_oilwaterpolymer.cc

@@ -55,6 +55,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

ebos/eclfluxmodule.hh

@@ -111,6 +111,7 @@ class EclTransExtensiveQuantities
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
     enum { numPhases = FluidSystem::numPhases };
     enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+    enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
 
     typedef Opm::MathToolbox<Evaluation> Toolbox;
@@ -271,6 +272,10 @@ protected:
 
             const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
             Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
+            if (enableExtbo) // added stability; particulary useful for solvent migrating in pure water
+                             // where the solvent fraction displays a 0/1 behaviour ...
+                pressureExterior += Toolbox::value(rhoAvg)*(distZ*g);
+            else
                 pressureExterior += rhoAvg*(distZ*g);
 
             pressureDifference_[phaseIdx] = pressureExterior - pressureInterior;

ebos/ecloutputblackoilmodule.hh

@@ -396,6 +396,14 @@ public:
             cFoam_.resize(bufferSize, 0.0);
         if (getPropValue<TypeTag, Properties::EnableBrine>())
             cSalt_.resize(bufferSize, 0.0);
+        if (getPropValue<TypeTag, Properties::EnableExtbo>()) {
+            extboX_.resize(bufferSize, 0.0);
+            extboY_.resize(bufferSize, 0.0);
+            extboZ_.resize(bufferSize, 0.0);
+            mFracOil_.resize(bufferSize, 0.0);
+            mFracGas_.resize(bufferSize, 0.0);
+            mFracCo2_.resize(bufferSize, 0.0);
+        }
 
         if (simulator_.problem().vapparsActive())
             soMax_.resize(bufferSize, 0.0);
@@ -629,6 +637,11 @@ public:
                 if (viscosity_[phaseIdx].size() == 0)
                     continue;
 
+                if (extboX_.size() > 0 && phaseIdx==oilPhaseIdx)
+                    viscosity_[phaseIdx][globalDofIdx] = Opm::getValue(intQuants.oilViscosity());
+                else if (extboX_.size() > 0 && phaseIdx==gasPhaseIdx)
+                    viscosity_[phaseIdx][globalDofIdx] = Opm::getValue(intQuants.gasViscosity());
+                else
                     viscosity_[phaseIdx][globalDofIdx] = Opm::getValue(fs.viscosity(phaseIdx));
                 Opm::Valgrind::CheckDefined(viscosity_[phaseIdx][globalDofIdx]);
             }
@@ -657,6 +670,34 @@ public:
                 cSalt_[globalDofIdx] = fs.saltConcentration().value();
             }
 
+            if (extboX_.size() > 0) {
+                extboX_[globalDofIdx] = intQuants.xVolume().value();
+            }
+
+            if (extboY_.size() > 0) {
+                extboY_[globalDofIdx] = intQuants.yVolume().value();
+            }
+
+            if (extboZ_.size() > 0) {
+                extboZ_[globalDofIdx] = intQuants.zFraction().value();
+            }
+
+            if (mFracCo2_.size() > 0) {
+                const Scalar stdVolOil = Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))
+                                       + Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*Opm::getValue(fs.Rv());
+                const Scalar stdVolGas = Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*(1.0-intQuants.yVolume().value())
+                                       + Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))*Opm::getValue(fs.Rs())*(1.0-intQuants.xVolume().value());
+                const Scalar stdVolCo2 = Opm::getValue(fs.saturation(gasPhaseIdx))*Opm::getValue(fs.invB(gasPhaseIdx))*intQuants.yVolume().value()
+                                       + Opm::getValue(fs.saturation(oilPhaseIdx))*Opm::getValue(fs.invB(oilPhaseIdx))*Opm::getValue(fs.Rs())*intQuants.xVolume().value();
+                const Scalar rhoO= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+                const Scalar rhoG= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+                const Scalar rhoCO2= intQuants.zRefDensity();
+                const Scalar stdMassTotal= 1.0e-10 + stdVolOil*rhoO + stdVolGas*rhoG + stdVolCo2*rhoCO2;
+                mFracOil_[globalDofIdx] = stdVolOil*rhoO/stdMassTotal;
+                mFracGas_[globalDofIdx] = stdVolGas*rhoG/stdMassTotal;
+                mFracCo2_[globalDofIdx] = stdVolCo2*rhoCO2/stdMassTotal;
+            }
+
             if (bubblePointPressure_.size() > 0) {
                 try {
                     bubblePointPressure_[globalDofIdx] = Opm::getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
@@ -1045,6 +1086,24 @@ public:
         if (sSol_.size() > 0)
             sol.insert ("SSOLVENT", Opm::UnitSystem::measure::identity, std::move(sSol_), Opm::data::TargetType::RESTART_SOLUTION);
 
+        if (extboX_.size() > 0)
+            sol.insert ("SS_X", Opm::UnitSystem::measure::identity, std::move(extboX_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (extboY_.size() > 0)
+            sol.insert ("SS_Y", Opm::UnitSystem::measure::identity, std::move(extboY_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (extboZ_.size() > 0)
+            sol.insert ("SS_Z", Opm::UnitSystem::measure::identity, std::move(extboZ_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (mFracOil_.size() > 0)
+            sol.insert ("STD_OIL", Opm::UnitSystem::measure::identity, std::move(mFracOil_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (mFracGas_.size() > 0)
+            sol.insert ("STD_GAS", Opm::UnitSystem::measure::identity, std::move(mFracGas_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (mFracCo2_.size() > 0)
+            sol.insert ("STD_CO2", Opm::UnitSystem::measure::identity, std::move(mFracCo2_), Opm::data::TargetType::RESTART_SOLUTION);
+
         if (cPolymer_.size() > 0)
             sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, std::move(cPolymer_), Opm::data::TargetType::RESTART_SOLUTION);
 
@@ -2290,6 +2349,12 @@ private:
     ScalarBuffer viscosity_[numPhases];
     ScalarBuffer relativePermeability_[numPhases];
     ScalarBuffer sSol_;
+    ScalarBuffer extboX_;
+    ScalarBuffer extboY_;
+    ScalarBuffer extboZ_;
+    ScalarBuffer mFracOil_;
+    ScalarBuffer mFracGas_;
+    ScalarBuffer mFracCo2_;
     ScalarBuffer cPolymer_;
     ScalarBuffer cFoam_;
     ScalarBuffer cSalt_;

ebos/eclproblem.hh

@@ -534,6 +534,10 @@ template<class TypeTag>
 struct EnableFoam<TypeTag, TTag::EclBaseProblem> {
     static constexpr bool value = false;
 };
+template<class TypeTag>
+struct EnableExtbo<TypeTag, TTag::EclBaseProblem> {
+    static constexpr bool value = false;
+};
 
 // disable thermal flux boundaries by default
 template<class TypeTag>
@@ -613,6 +617,7 @@ class EclProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
     enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
+    enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
@@ -647,6 +652,7 @@ class EclProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     typedef BlackOilPolymerModule<TypeTag> PolymerModule;
     typedef BlackOilFoamModule<TypeTag> FoamModule;
     typedef BlackOilBrineModule<TypeTag> BrineModule;
+    typedef BlackOilExtboModule<TypeTag> ExtboModule;
 
     typedef typename EclEquilInitializer<TypeTag>::ScalarFluidState InitialFluidState;
 
@@ -812,6 +818,7 @@ public:
         PolymerModule::initFromState(vanguard.eclState());
         FoamModule::initFromState(vanguard.eclState());
         BrineModule::initFromState(vanguard.eclState());
+        ExtboModule::initFromState(vanguard.eclState());
 
         // create the ECL writer
         eclWriter_.reset(new EclWriterType(simulator));
@@ -2176,6 +2183,11 @@ private:
             else if (!enablePolymer && deck.hasKeyword("POLYMER"))
                 throw std::runtime_error("The deck enables the polymer option, but the simulator is compiled without it.");
 
+            if (enableExtbo && !deck.hasKeyword("PVTSOL"))
+                throw std::runtime_error("The simulator requires the extendedBO option to be enabled, but the deck does not.");
+            else if (!enableExtbo && deck.hasKeyword("PVTSOL"))
+                throw std::runtime_error("The deck enables the extendedBO option, but the simulator is compiled without it.");
+
             if (deck.hasKeyword("TEMP") && deck.hasKeyword("THERMAL"))
                 throw std::runtime_error("The deck enables both, the TEMP and the THERMAL options, but they are mutually exclusive.");
 

flow/flow_ebos_extbo.cpp

@@ -0,0 +1,80 @@
+/*
+  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_ebos_extbo.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#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 EclFlowExtboProblem {
+    using InheritsFrom = std::tuple<EclFlowProblem>;
+};
+}
+template<class TypeTag>
+struct EnableExtbo<TypeTag, TTag::EclFlowExtboProblem> {
+    static constexpr bool value = true;
+};
+}}
+
+namespace Opm {
+void flowEbosExtboSetDeck(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::EclFlowExtboProblem;
+    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 flowEbosExtboMain(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();
+
+    // initialize MPI, finalize is done automatically on exit
+#if HAVE_DUNE_FEM
+    Dune::Fem::MPIManager::initialize(argc, argv);
+#else
+    Dune::MPIHelper::instance(argc, argv).rank();
+#endif
+
+    Opm::FlowMainEbos<Properties::TTag::EclFlowExtboProblem>
+        mainfunc {argc, argv, outputCout, outputFiles};
+    return mainfunc.execute();
+}
+
+}

flow/flow_ebos_extbo.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_EBOS_EXTBO_HPP
+#define FLOW_EBOS_EXTBO_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 flowEbosExtboSetDeck(double setupTime, std::unique_ptr<Deck> deck,
+                            std::unique_ptr<EclipseState> eclState,
+                            std::unique_ptr<Schedule> schedule,
+                            std::unique_ptr<SummaryConfig> summaryConfig);
+int flowEbosExtboMain(int argc, char** argv, bool outputCout, bool outputFiles);
+}
+
+#endif // FLOW_EBOS_EXTBO_HPP

flow/flow_ebos_gasoil.cpp

@@ -52,6 +52,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

flow/flow_ebos_oilwater.cpp

@@ -52,6 +52,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

flow/flow_ebos_oilwater_brine.cpp

@@ -55,6 +55,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

flow/flow_ebos_oilwater_polymer.cpp

@@ -55,6 +55,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

flow/flow_ebos_oilwater_polymer_injectivity.cpp

@@ -61,6 +61,7 @@ private:
 
 public:
     typedef Opm::BlackOilTwoPhaseIndices<0,
+                                         0,
                                          2,
                                          0,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

flow/flow_onephase.cpp

@@ -43,6 +43,7 @@ private:
 
 public:
     typedef Opm::BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

flow/flow_onephase_energy.cpp

@@ -46,6 +46,7 @@ private:
 
 public:
     typedef Opm::BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          getPropValue<TypeTag, Properties::EnableEnergy>(),
                                          getPropValue<TypeTag, Properties::EnableFoam>(),

opm/core/props/BlackoilPhases.hpp

@@ -34,10 +34,10 @@ namespace Opm
         // to and from active ones. That said, they are not considered by num_phases or
         // MaxNumPhases. The crypto phases which are currently implemented are solvent,
         // polymer, energy, polymer molecular weight, foam and brine.
-        static const int NumCryptoPhases = 6;
+        static const int NumCryptoPhases = 7;
 
         // enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 };
-        enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6, Foam = 7, Brine = 8 };
+        enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6, Foam = 7, Brine = 8, ZFraction = 9 };
     };
 
     struct PhaseUsage : public BlackoilPhases
@@ -52,6 +52,8 @@ namespace Opm
         bool has_polymermw;
         bool has_foam;
         bool has_brine;
+        bool has_zFraction;
+
     };
 
     /// Check or assign presence of a formed, free phase.  Limited to

opm/core/props/phaseUsageFromDeck.hpp

@@ -132,6 +132,16 @@ namespace Opm
         else
             pu.phase_pos[BlackoilPhases::Brine] = -1;
 
+        // Add zFraction info
+        pu.has_zFraction = phases.active(Phase::ZFRACTION);
+        if (pu.has_zFraction) {
+            pu.phase_pos[BlackoilPhases::ZFraction] = numActivePhases;
+            ++ numActivePhases;
+        }
+        else
+            pu.phase_pos[BlackoilPhases::ZFraction] = -1;
+
+
         return pu;
     }
 

opm/simulators/flow/Main.hpp

@@ -29,6 +29,7 @@
 #  include <flow/flow_ebos_oilwater.hpp>
 #  include <flow/flow_ebos_solvent.hpp>
 #  include <flow/flow_ebos_polymer.hpp>
+#  include <flow/flow_ebos_extbo.hpp>
 #  include <flow/flow_ebos_foam.hpp>
 #  include <flow/flow_ebos_brine.hpp>
 #  include <flow/flow_ebos_oilwater_brine.hpp>
@@ -292,6 +293,14 @@ namespace Opm
                                             std::move(summaryConfig_));
                 return Opm::flowEbosSolventMain(argc_, argv_, outputCout_, outputFiles_);
             }
+            // Extended BO case
+            else if ( phases.active( Opm::Phase::ZFRACTION ) ) {
+                Opm::flowEbosExtboSetDeck(setupTime_, std::move(deck_),
+                                            std::move(eclipseState_),
+                                            std::move(schedule_),
+                                            std::move(summaryConfig_));
+                return Opm::flowEbosExtboMain(argc_, argv_, outputCout_, outputFiles_);
+            }
             // Energy case
             else if (eclipseState_->getSimulationConfig().isThermal()) {
                 Opm::flowEbosEnergySetDeck(setupTime_, std::move(deck_),

opm/simulators/wells/BlackoilWellModel.hpp

@@ -296,6 +296,7 @@ namespace Opm {
             const ModelParameters param_;
             bool terminal_output_;
             bool has_solvent_;
+            bool has_zFraction_;
             bool has_polymer_;
             std::vector<int> pvt_region_idx_;
             PhaseUsage phase_usage_;

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -33,6 +33,7 @@ namespace Opm {
     BlackoilWellModel(Simulator& ebosSimulator)
         : ebosSimulator_(ebosSimulator)
         , has_solvent_(getPropValue<TypeTag, Properties::EnableSolvent>())
+        , has_zFraction_(getPropValue<TypeTag, Properties::EnableExtbo>())
         , has_polymer_(getPropValue<TypeTag, Properties::EnablePolymer>())
     {
         terminal_output_ = false;

opm/simulators/wells/StandardWell.hpp

@@ -33,6 +33,7 @@
 
 #include <opm/models/blackoil/blackoilpolymermodules.hh>
 #include <opm/models/blackoil/blackoilsolventmodules.hh>
+#include <opm/models/blackoil/blackoilextbomodules.hh>
 #include <opm/models/blackoil/blackoilfoammodules.hh>
 #include <opm/models/blackoil/blackoilbrinemodules.hh>
 
@@ -74,6 +75,7 @@ namespace Opm
         using Base::numEq;
 
         using Base::has_solvent;
+        using Base::has_zFraction;
         using Base::has_polymer;
         using Base::has_foam;
         using Base::has_brine;
@@ -89,9 +91,10 @@ namespace Opm
         static const int numEnergyEq = Indices::numEnergy;
         static const int numFoamEq = Indices::numFoam;
         static const int numBrineEq = Indices::numBrine;
+        static const int numExtbos = Indices::numExtbos;
 
         // number of the conservation equations
-        static const int numWellConservationEq = numEq - numPolymerEq - numEnergyEq - numFoamEq - numBrineEq;
+        static const int numWellConservationEq = numEq - numPolymerEq - numEnergyEq - numFoamEq - numBrineEq - numExtbos;
         // number of the well control equations
         static const int numWellControlEq = 1;
         // number of the well equations that will always be used
@@ -147,6 +150,7 @@ namespace Opm
         typedef DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + numEq + 1> EvalWell;
 
         using Base::contiSolventEqIdx;
+        using Base::contiZfracEqIdx;
         using Base::contiPolymerEqIdx;
         using Base::contiFoamEqIdx;
         using Base::contiBrineEqIdx;
@@ -523,6 +527,7 @@ namespace Opm
                                  std::vector<RateVector>& connectionRates,
                                  std::vector<EvalWell>& cq_s,
                                  EvalWell& water_flux_s,
+                                 EvalWell& cq_s_zfrac_effective,
                                  Opm::DeferredLogger& deferred_logger) const;
 
         // check whether the well is operable under the current reservoir condition

opm/simulators/wells/StandardWell_impl.hpp

@@ -388,6 +388,12 @@ namespace Opm
             b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
         }
 
+        if (has_zFraction && this->isInjector()) {
+            const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+            b_perfcells_dense[gasCompIdx] *= (1.0 - wsolvent());
+            b_perfcells_dense[gasCompIdx] += wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
+        }
+
         // Pressure drawdown (also used to determine direction of flow)
         const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
         EvalWell drawdown = pressure - well_pressure;
@@ -605,7 +611,8 @@ namespace Opm
             // Calculate perforation quantities.
             std::vector<EvalWell> cq_s(num_components_, {numWellEq_ + numEq, 0.0});
             EvalWell water_flux_s{numWellEq_ + numEq, 0.0};
-            calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, deferred_logger);
+            EvalWell cq_s_zfrac_effective{numWellEq_ + numEq, 0.0};
+            calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
 
             // Equation assembly for this perforation.
             if (has_polymer && this->has_polymermw && this->isInjector()) {
@@ -639,6 +646,12 @@ namespace Opm
                     well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
                 }
             }
+
+            if (has_zFraction) {
+                for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) {
+                    duneC_[0][cell_idx][pvIdx][contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+numEq);
+                }
+            }
         }
         // Update the connection
         connectionRates_ = connectionRates;
@@ -686,6 +699,7 @@ namespace Opm
                         std::vector<RateVector>& connectionRates,
                         std::vector<EvalWell>& cq_s,
                         EvalWell& water_flux_s,
+                        EvalWell& cq_s_zfrac_effective,
                         Opm::DeferredLogger& deferred_logger) const
     {
         const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
@@ -812,6 +826,22 @@ namespace Opm
             connectionRates[perf][contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
         }
 
+        if (has_zFraction) {
+            // TODO: the application of well efficiency factor has not been tested with an example yet
+            const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+            cq_s_zfrac_effective = cq_s[gasCompIdx];
+            if (this->isInjector()) {
+                cq_s_zfrac_effective *= wsolvent();
+            } else if (cq_s_zfrac_effective.value() != 0.0) {
+                const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value();
+                cq_s_zfrac_effective *= extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
+            }
+            well_state.perfRateSolvent()[first_perf_ + perf] = cq_s_zfrac_effective.value();
+
+            cq_s_zfrac_effective *= well_efficiency_factor_;
+            connectionRates[perf][contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective);
+        }
+
         if (has_brine) {
             // TODO: the application of well efficiency factor has not been tested with an example yet
             const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
@@ -1978,7 +2008,7 @@ namespace Opm
                     const double oilrate = std::abs(well_state.wellRates()[oilpos_well]); //in order to handle negative rates in producers
                     rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
                     if (oilrate > 0) {
-                        const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w);
+                        const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
                         double rv = 0.0;
                         if (gasrate > 0) {
                             rv = oilrate / gasrate;
@@ -2003,7 +2033,7 @@ namespace Opm
                 if (gasPresent) {
                     rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
                     const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases;
-                    const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w);
+                    const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
                     if (gasrate > 0) {
                         const double oilrate = std::abs(well_state.wellRates()[oilpos_well]);
                         double rs = 0.0;
@@ -2224,7 +2254,7 @@ namespace Opm
     {
         // the following implementation assume that the polymer is always after the w-o-g phases
         // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells
-        assert((int(B_avg.size()) == num_components_) || has_polymer || has_energy || has_foam || has_brine);
+        assert((int(B_avg.size()) == num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction);
 
         const double tol_wells = param_.tolerance_wells_;
         const double maxResidualAllowed = param_.max_residual_allowed_;
@@ -2900,7 +2930,8 @@ namespace Opm
                 primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
             }
             if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ;
+                primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]]
+                                                 - (has_solvent ? well_state.solventWellRate(well_index) : 0.0) ) / total_well_rate ;
             }
             if (has_solvent) {
                 primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ;
@@ -2919,8 +2950,9 @@ namespace Opm
 
                 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                     if (phase == InjectorType::GAS) {
-                        primary_variables_[GFrac] = 1.0 - wsolvent();
+                        primary_variables_[GFrac] = 1.0;
                         if (has_solvent) {
+                            primary_variables_[GFrac] = 1.0 - wsolvent();
                             primary_variables_[SFrac] = wsolvent();
                         }
                     } else {

opm/simulators/wells/WellInterface.hpp

@@ -92,6 +92,7 @@ namespace Opm
         typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
 
         static const bool has_solvent = getPropValue<TypeTag, Properties::EnableSolvent>();
+        static const bool has_zFraction = getPropValue<TypeTag, Properties::EnableExtbo>();
         static const bool has_polymer = getPropValue<TypeTag, Properties::EnablePolymer>();
         static const bool has_energy = getPropValue<TypeTag, Properties::EnableEnergy>();
         static const bool has_temperature = getPropValue<TypeTag, Properties::EnableTemperature>();
@@ -100,6 +101,7 @@ namespace Opm
         static const bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>();
         static const bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>();
         static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
+        static const int contiZfracEqIdx = Indices::contiZfracEqIdx;
         static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
         // index for the polymer molecular weight continuity equation
         static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;

opm/simulators/wells/WellInterface_impl.hpp

@@ -94,7 +94,7 @@ namespace Opm
 
         wsolvent_ = 0.0;
 
-        if (has_solvent && well.isInjector()) {
+        if ((has_solvent || has_zFraction) && well.isInjector()) {
             auto injectorType = well_ecl_.injectorType();
             if (injectorType == InjectorType::GAS) {
                 wsolvent_ = well_ecl_.getSolventFraction();

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -163,7 +163,6 @@ namespace Opm
 
             perfRateBrine_.clear();
             perfRateBrine_.resize(nperf, 0.0);
-
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if (prevState && !prevState->wellMap().empty()) {
@@ -576,7 +575,7 @@ namespace Opm
                     well.rates.set( rt::well_potential_gas, this->well_potentials_[well_rate_index + pu.phase_pos[Gas]] );
                 }
 
-                if ( pu.has_solvent ) {
+                if ( pu.has_solvent || pu.has_zFraction) {
                     well.rates.set( rt::solvent, solventWellRate(w) );
                 }
 
@@ -621,7 +620,7 @@ namespace Opm
                     if ( pu.has_brine ) {
                         comp.rates.set( rt::brine, this->perfRateBrine()[wt.second[1] + local_comp_index]);
                     }
-                    if ( pu.has_solvent ) {
+                    if ( pu.has_solvent || pu.has_zFraction) {
                         comp.rates.set( rt::solvent, this->perfRateSolvent()[wt.second[1] + local_comp_index]);
                     }