Merge pull request #5851 from svenn-t/dummy_water_compositional

Compositional simulator expanded to three phases

CMakeLists.txt

@@ -634,6 +634,7 @@ opm_add_test(flowexp_blackoil
 set(FLOWEXP_COMPONENTS_SOURCES)
 foreach(component IN LISTS OPM_COMPILE_COMPONENTS)
   list(APPEND FLOWEXP_COMPONENTS_SOURCES flowexperimental/comp/flowexp_comp${component}.cpp)
+  list(APPEND FLOWEXP_COMPONENTS_SOURCES flowexperimental/comp/flowexp_comp${component}_2p.cpp)
 endforeach()
 
 opm_add_test(flowexp_comp

examples/problems/co2ptflashproblem.hh

@@ -38,7 +38,8 @@
 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
 #include <opm/material/fluidmatrixinteractions/BrooksCorey.hpp>
 #include <opm/material/constraintsolvers/PTFlash.hpp> 
-#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
+#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp>
 #include <opm/material/common/Valgrind.hpp>
 #include <opm/models/immiscible/immisciblemodel.hh>
 #include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
@@ -76,6 +77,14 @@ struct NumComp<TypeTag, TTag::CO2PTBaseProblem> {
     static constexpr int value = 3;
 };
 
+template <class TypeTag, class MyTypeTag>
+struct EnableDummyWater { using type = UndefinedProperty; };
+
+template <class TypeTag>
+struct EnableDummyWater<TypeTag, TTag::CO2PTBaseProblem> {
+    static constexpr bool value = true;
+};
+
 // Set the grid type: --->2D
 template <class TypeTag>
 struct Grid<TypeTag, TTag::CO2PTBaseProblem> { using type = Dune::YaspGrid</*dim=*/2>; };
@@ -104,9 +113,10 @@ struct FluidSystem<TypeTag, TTag::CO2PTBaseProblem>
 private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
+    static constexpr bool enable_water = getPropValue<TypeTag, Properties::EnableDummyWater>();
 
 public:
-    using type = Opm::GenericOilGasFluidSystem<Scalar, num_comp>;
+    using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp, enable_water>;
 };
 
 // Set the material Law
@@ -118,8 +128,8 @@ private:
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
 
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
+    using Traits = Opm::ThreePhaseMaterialTraits<Scalar,
-                                               //  /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx, // TODO
+                                               /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
                                                /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
                                                /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
 
@@ -195,6 +205,8 @@ class CO2PTProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
     enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { pressure0Idx = Indices::pressure0Idx };
+    enum { z0Idx = Indices::z0Idx };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
@@ -237,6 +249,21 @@ public:
         porosity_ = 0.1;
     }
 
+    void initWaterPVT()
+    {
+        using WaterPvt = typename FluidSystem::WaterPvt;
+        std::shared_ptr<WaterPvt> waterPvt = std::make_shared<WaterPvt>();
+        waterPvt->setApproach(WaterPvtApproach::ConstantCompressibilityWater);
+        auto& ccWaterPvt = waterPvt->template getRealPvt<WaterPvtApproach::ConstantCompressibilityWater>();
+        ccWaterPvt.setNumRegions(/*numPvtRegions=*/1);
+        Scalar rhoRefW = 1037.0; // [kg]
+        ccWaterPvt.setReferenceDensities(/*regionIdx=*/0, /*rhoRefO=*/Scalar(0.0), /*rhoRefG=*/Scalar(0.0), rhoRefW);
+        ccWaterPvt.setViscosity(/*regionIdx=*/0, 9.6e-4);
+        ccWaterPvt.setCompressibility(/*regionIdx=*/0, 1.450377e-10);
+        waterPvt->initEnd();
+        FluidSystem::setWaterPvt(waterPvt);
+    }
+
     template <class Context>
     const DimVector&
     gravity([[maybe_unused]]const Context& context,
@@ -262,8 +289,12 @@ public:
     void finishInit()
     {
         ParentType::finishInit();
+
         // initialize fixed parameters; temperature, permeability, porosity
         initPetrophysics();
+
+        // Initialize water pvt
+        initWaterPVT();
     }
 
     /*!
@@ -358,8 +389,8 @@ public:
 
         // Calculate storage terms
          PrimaryVariables storageL, storageG;
-         this->model().globalPhaseStorage(storageL, /*phaseIdx=*/0);
+         this->model().globalPhaseStorage(storageL, /*phaseIdx=*/oilPhaseIdx);
-         this->model().globalPhaseStorage(storageG, /*phaseIdx=*/1);
+         this->model().globalPhaseStorage(storageG, /*phaseIdx=*/gasPhaseIdx);
 
          // Write mass balance information for rank 0
         //  if (this->gridView().comm().rank() == 0) {
@@ -456,12 +487,12 @@ private:
         int prod = Parameters::Get<Parameters::CellsX>() - 1;
         int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
         ComponentVector comp;
-        comp[0] = Evaluation::createVariable(0.5, 1);
+        comp[0] = Evaluation::createVariable(0.5, z0Idx);
-        comp[1] = Evaluation::createVariable(0.3, 2);
+        comp[1] = Evaluation::createVariable(0.3, z0Idx + 1);
         comp[2] = 1. - comp[0] - comp[1];
         if (spatialIdx == inj) {
-            comp[0] = Evaluation::createVariable(0.99, 1);
+            comp[0] = Evaluation::createVariable(0.99, z0Idx);
-            comp[1] = Evaluation::createVariable(0.01 - 1e-3, 2);
+            comp[1] = Evaluation::createVariable(0.01 - 1e-3, z0Idx + 1);
             comp[2] = 1. - comp[0] - comp[1];
         }
 
@@ -475,10 +506,11 @@ private:
         if (spatialIdx == prod) {
             p0 *= 0.5;
         }
-        Evaluation p_init = Evaluation::createVariable(p0, 0);
+        Evaluation p_init = Evaluation::createVariable(p0, pressure0Idx);
 
         fs.setPressure(FluidSystem::oilPhaseIdx, p_init);
         fs.setPressure(FluidSystem::gasPhaseIdx, p_init);
+        fs.setPressure(FluidSystem::waterPhaseIdx, p_init);
 
         fs.setTemperature(temperature_);
 

flowexperimental/comp/flowexp_comp.cpp

@@ -36,10 +36,13 @@
 
 template <int compileTimeComponent>
 std::tuple<bool, int>
-runComponent(int runtimeComponent, int argc, char** argv)
+runComponent(int runtimeComponent, bool water, int argc, char** argv)
 {
     if (runtimeComponent == compileTimeComponent) {
-        return std::make_tuple(true, Opm::dispatchFlowExpComp<compileTimeComponent>(argc, argv));
+        if (water)
+            return std::make_tuple(true, Opm::dispatchFlowExpComp<compileTimeComponent, true>(argc, argv));
+        else
+            return std::make_tuple(true, Opm::dispatchFlowExpComp<compileTimeComponent, false>(argc, argv));
     }
     return std::make_tuple(false, EXIT_FAILURE);
 }
@@ -63,18 +66,21 @@ runComponent(int runtimeComponent, int argc, char** argv)
  */
 template <int currentCompileTimeComponent, int nextComponent, int... components>
 std::tuple<bool, int>
-runComponent(int runtimecomponent, int argc, char** argv)
+runComponent(int runtimecomponent, bool water, int argc, char** argv)
 {
     if (currentCompileTimeComponent == runtimecomponent) {
-        return std::make_tuple(true, Opm::dispatchFlowExpComp<currentCompileTimeComponent>(argc, argv));
+        if (water)
+            return std::make_tuple(true, Opm::dispatchFlowExpComp<currentCompileTimeComponent, true>(argc, argv));
+        else
+            return std::make_tuple(true, Opm::dispatchFlowExpComp<currentCompileTimeComponent, false>(argc, argv));
     }
-    return runComponent<nextComponent, components...>(runtimecomponent, argc, argv);
+    return runComponent<nextComponent, components...>(runtimecomponent, water, argc, argv);
 }
 
 int
 main(int argc, char** argv)
 {
-    using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<0>;
+    using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<0, true>;
     Opm::registerEclTimeSteppingParameters<double>();
 
     // At the moment, this is probably as optimal as can be.
@@ -90,9 +96,11 @@ main(int argc, char** argv)
         = Opm::Parser {}.parseFile(inputFilename, Opm::ParseContext {}, std::vector {Opm::Ecl::SectionType::RUNSPEC});
     const auto runspec = Opm::Runspec(deck);
     const auto numComps = runspec.numComps();
+    const auto& phases = runspec.phases();
+    const auto wat = phases.active(Opm::Phase::WATER);
 
     auto [componentSupported, executionStatus]
-        = runComponent<OPM_COMPILE_COMPONENTS_TEMPLATE_LIST>(numComps, argc, argv);
+        = runComponent<OPM_COMPILE_COMPONENTS_TEMPLATE_LIST>(numComps, wat, argc, argv);
 
     if (!componentSupported) {
         fmt::print("Deck has {} components, not supported. In this build of the simulator, we support the "

flowexperimental/comp/flowexp_comp.hpp

@@ -20,7 +20,7 @@
 #define FLOWEXP_COMP_HPP
 
 #include <opm/material/constraintsolvers/PTFlash.hpp>
-#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
 
 #include <opm/models/discretization/common/baseauxiliarymodule.hh>
 #include <opm/models/ptflash/flashmodel.hh>
@@ -39,7 +39,7 @@
 // suggestTimeStep is taken from newton solver in problem.limitTimestep
 namespace Opm {
 
-template<int numComp>
+template<int numComp, bool EnableWater>
 int dispatchFlowExpComp(int argc, char** argv);
 
 }
@@ -47,15 +47,15 @@ int dispatchFlowExpComp(int argc, char** argv);
 namespace Opm::Properties {
 namespace TTag {
 
-template<int NumComp>
+template<int NumComp, bool EnableWater>
 struct FlowExpCompProblem {
    using InheritsFrom = std::tuple<FlowBaseProblemComp, FlashModel>;
 };
 
 }
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct SparseMatrixAdapter<TypeTag, TTag::FlowExpCompProblem<NumComp>>
+struct SparseMatrixAdapter<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>>
 {
 private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -109,30 +109,30 @@ struct LocalLinearizerSplice<TypeTag, TTag::FlowExpCompProblem>
 #endif
 
 // Set the problem property
-template <class TypeTag, int NumComp>
+template <class TypeTag, int NumComp, bool EnableWater>
-struct Problem<TypeTag, TTag::FlowExpCompProblem<NumComp>>
+struct Problem<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>>
 {
     using type = FlowProblemComp<TypeTag>;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct AquiferModel<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct AquiferModel<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     using type = EmptyModel<TypeTag>;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct WellModel<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct WellModel<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     using type = EmptyModel<TypeTag>;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct TracerModel<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct TracerModel<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     using type = EmptyModel<TypeTag>;
 };
 
 
-template <class TypeTag, int NumComp>
+template <class TypeTag, int NumComp, bool EnableWater>
-struct FlashSolver<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct FlashSolver<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
 private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
@@ -147,10 +147,18 @@ template <class TypeTag, class MyTypeTag>
 struct NumComp { using type = UndefinedProperty; };
 
 // TODO: this is unfortunate, have to check why we need to hard-code it
-template <class TypeTag, int NumComp_>
+template <class TypeTag, int NumComp_, bool EnableWater_>
-struct NumComp<TypeTag, TTag::FlowExpCompProblem<NumComp_>> {
+struct NumComp<TypeTag, TTag::FlowExpCompProblem<NumComp_, EnableWater_>> {
     static constexpr int value = NumComp_;
 };
+
+template <class TypeTag, class MyTypeTag>
+struct EnableDummyWater { using type = UndefinedProperty; };
+
+template <class TypeTag, int NumComp_, bool EnableWater_>
+struct EnableDummyWater<TypeTag, TTag::FlowExpCompProblem<NumComp_, EnableWater_>> {
+    static constexpr bool value = EnableWater_;
+};
 #if 0
 struct Temperature { using type = UndefinedProperty; };
 
@@ -161,26 +169,29 @@ struct Temperature { using type = UndefinedProperty; };
  };
 #endif
 
-template <class TypeTag, int NumComp_>
+template <class TypeTag, int NumComp_, bool EnableWater_>
-struct FluidSystem<TypeTag, TTag::FlowExpCompProblem<NumComp_>>
+struct FluidSystem<TypeTag, TTag::FlowExpCompProblem<NumComp_, EnableWater_>>
 {
 private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
+    static constexpr bool enable_water = getPropValue<TypeTag, Properties::EnableDummyWater>();
 
 public:
-    using type = Opm::GenericOilGasFluidSystem<Scalar, num_comp>;
+    using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp, enable_water>;
 };
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableMech<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableMech<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableDisgasInWater<TypeTag, TTag::FlowExpCompProblem<NumComp>> { static constexpr bool value = false; };
+struct EnableDisgasInWater<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> { 
+    static constexpr bool value = false; 
+};
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct Stencil<TypeTag, TTag::FlowExpCompProblem<NumComp>>
+struct Stencil<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>>
 {
 private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -190,62 +201,62 @@ public:
     using type = EcfvStencil<Scalar, GridView>;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableApiTracking<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableApiTracking<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableTemperature<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableTemperature<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableSaltPrecipitation<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableSaltPrecipitation<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnablePolymerMW<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnablePolymerMW<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnablePolymer<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnablePolymer<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableDispersion<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableDispersion<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableBrine<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableBrine<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableVapwat<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableVapwat<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableSolvent<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableSolvent<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableEnergy<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableEnergy<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableFoam<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableFoam<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableExtbo<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableExtbo<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
-template<class TypeTag, int NumComp>
+template<class TypeTag, int NumComp, bool EnableWater>
-struct EnableMICP<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
+struct EnableMICP<TypeTag, TTag::FlowExpCompProblem<NumComp, EnableWater>> {
     static constexpr bool value = false;
 };
 

flowexperimental/comp/flowexp_comp2.cpp

@@ -28,9 +28,9 @@
 namespace Opm {
 
 template<>
-int dispatchFlowExpComp<2>(int argc, char** argv)
+int dispatchFlowExpComp<2, true>(int argc, char** argv)
 {
-    return start<Properties::TTag::FlowExpCompProblem<2>>(argc, argv, false);
+    return start<Properties::TTag::FlowExpCompProblem<2, true>>(argc, argv, false);
 }
 
 }

flowexperimental/comp/flowexp_comp2_2p.cpp

@@ -0,0 +1,36 @@
+/*
+  Copyright 2024, SINTEF Digital
+
+  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 <opm/models/utils/start.hh>
+
+#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
+
+#include "flowexp_comp.hpp"
+
+namespace Opm {
+
+template<>
+int dispatchFlowExpComp<2, false>(int argc, char** argv)
+{
+    return start<Properties::TTag::FlowExpCompProblem<2, false>>(argc, argv, false);
+}
+
+}

flowexperimental/comp/flowexp_comp3.cpp

@@ -28,9 +28,9 @@
 namespace Opm {
 
 template<>
-int dispatchFlowExpComp<3>(int argc, char** argv)
+int dispatchFlowExpComp<3, true>(int argc, char** argv)
 {
-    return start<Properties::TTag::FlowExpCompProblem<3>>(argc, argv, false);
+    return start<Properties::TTag::FlowExpCompProblem<3, true>>(argc, argv, false);
 }
 
 }

flowexperimental/comp/flowexp_comp3_2p.cpp

@@ -0,0 +1,36 @@
+/*
+  Copyright 2024, SINTEF Digital
+
+  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 <opm/models/utils/start.hh>
+
+#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
+
+#include "flowexp_comp.hpp"
+
+namespace Opm {
+
+template<>
+int dispatchFlowExpComp<3, false>(int argc, char** argv)
+{
+    return start<Properties::TTag::FlowExpCompProblem<3, false>>(argc, argv, false);
+}
+
+}

flowexperimental/comp/flowexp_comp4.cpp

@@ -28,9 +28,9 @@
 namespace Opm {
 
 template<>
-int dispatchFlowExpComp<4>(int argc, char** argv)
+int dispatchFlowExpComp<4, true>(int argc, char** argv)
 {
-    return start<Properties::TTag::FlowExpCompProblem<4>>(argc, argv, false);
+    return start<Properties::TTag::FlowExpCompProblem<4, true>>(argc, argv, false);
 }
 
 }

flowexperimental/comp/flowexp_comp4_2p.cpp

@@ -0,0 +1,36 @@
+/*
+  Copyright 2024, SINTEF Digital
+
+  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 <opm/models/utils/start.hh>
+
+#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
+
+#include "flowexp_comp.hpp"
+
+namespace Opm {
+
+template<>
+int dispatchFlowExpComp<4, false>(int argc, char** argv)
+{
+    return start<Properties::TTag::FlowExpCompProblem<4, false>>(argc, argv, false);
+}
+
+}

flowexperimental/comp/flowexp_comp5.cpp

@@ -28,9 +28,9 @@
 namespace Opm {
 
 template<>
-int dispatchFlowExpComp<5>(int argc, char** argv)
+int dispatchFlowExpComp<5, true>(int argc, char** argv)
 {
-    return start<Properties::TTag::FlowExpCompProblem<5>>(argc, argv, false);
+    return start<Properties::TTag::FlowExpCompProblem<5, true>>(argc, argv, false);
 }
 
 }

flowexperimental/comp/flowexp_comp5_2p.cpp

@@ -0,0 +1,36 @@
+/*
+  Copyright 2024, SINTEF Digital
+
+  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 <opm/models/utils/start.hh>
+
+#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
+
+#include "flowexp_comp.hpp"
+
+namespace Opm {
+
+template<>
+int dispatchFlowExpComp<5, false>(int argc, char** argv)
+{
+    return start<Properties::TTag::FlowExpCompProblem<5, false>>(argc, argv, false);
+}
+
+}

flowexperimental/comp/flowexp_comp6.cpp

@@ -28,9 +28,9 @@
 namespace Opm {
 
 template<>
-int dispatchFlowExpComp<6>(int argc, char** argv)
+int dispatchFlowExpComp<6, true>(int argc, char** argv)
 {
-    return start<Properties::TTag::FlowExpCompProblem<6>>(argc, argv, false);
+    return start<Properties::TTag::FlowExpCompProblem<6, true>>(argc, argv, false);
 }
 
 }

flowexperimental/comp/flowexp_comp6_2p.cpp

@@ -0,0 +1,36 @@
+/*
+  Copyright 2024, SINTEF Digital
+
+  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 <opm/models/utils/start.hh>
+
+#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
+
+#include "flowexp_comp.hpp"
+
+namespace Opm {
+
+template<>
+int dispatchFlowExpComp<6, false>(int argc, char** argv)
+{
+    return start<Properties::TTag::FlowExpCompProblem<6, false>>(argc, argv, false);
+}
+
+}

flowexperimental/comp/flowexp_comp7.cpp

@@ -28,9 +28,9 @@
 namespace Opm {
 
 template<>
-int dispatchFlowExpComp<7>(int argc, char** argv)
+int dispatchFlowExpComp<7, true>(int argc, char** argv)
 {
-    return start<Properties::TTag::FlowExpCompProblem<7>>(argc, argv, false);
+    return start<Properties::TTag::FlowExpCompProblem<7, true>>(argc, argv, false);
 }
 
 }

flowexperimental/comp/flowexp_comp7_2p.cpp

@@ -0,0 +1,36 @@
+/*
+  Copyright 2024, SINTEF Digital
+
+  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 <opm/models/utils/start.hh>
+
+#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
+
+#include "flowexp_comp.hpp"
+
+namespace Opm {
+
+template<>
+int dispatchFlowExpComp<7, false>(int argc, char** argv)
+{
+    return start<Properties::TTag::FlowExpCompProblem<7, false>>(argc, argv, false);
+}
+
+}

opm/models/common/multiphasebaseproperties.hh

@@ -83,6 +83,8 @@ struct EnableDispersion { using type = UndefinedProperty; };
 //! Enable convective mixing?
 template<class TypeTag, class MyTypeTag>
 struct EnableConvectiveMixing { using type = UndefinedProperty; };
+template <class TypeTag, class MyTypeTag>
+struct EnableWater { using type = UndefinedProperty; };
 
 } // namespace Opm::Properties
 

opm/models/ptflash/flashindices.hh

@@ -48,16 +48,20 @@ class FlashIndices
 {
     static constexpr int numComponents = getPropValue<TypeTag, Properties::NumComponents>();
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableWater = getPropValue<TypeTag, Properties::EnableWater>() };
     using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>;
 
 public:
-    static constexpr bool waterEnabled = false;
+    //! All phases active (note: immiscible/"dummy" water phase)
+    static constexpr bool waterEnabled = enableWater;
     static constexpr bool gasEnabled = true;
     static constexpr bool oilEnabled = true;
-    static constexpr int waterPhaseIdx = -1;
+
-    static constexpr int numPhases = 2;
+    //! number of active phases
+    static constexpr int numPhases = enableWater ? 3 : 2;
+
     //! number of equations/primary variables
-    static const int numEq = numComponents + EnergyIndices::numEq_;
+    static const int numEq = numComponents + EnergyIndices::numEq_ + (enableWater ? 1 : 0);
 
     // Primary variable indices
 
@@ -66,6 +70,9 @@ public:
 
     //! Index of the molefraction of the first component
     static constexpr int z0Idx = pressure0Idx + 1;
+
+    //! Index of water saturation
+    static constexpr int water0Idx = enableWater ? z0Idx + numComponents - 1 : -1000;
     
     // equation indices
 

opm/models/ptflash/flashintensivequantities.hh

@@ -76,6 +76,9 @@ class FlashIntensiveQuantities
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { dimWorld = GridView::dimensionworld };
     enum { pressure0Idx = Indices::pressure0Idx };
+    enum { water0Idx = Indices::water0Idx};
+
+    static constexpr bool waterEnabled = Indices::waterEnabled;
 
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
@@ -216,19 +219,26 @@ public:
 
 
         // Update saturation
-        // \Note: the current implementation assume oil-gas system.
+        Evaluation Sw = 0.0;
+        if constexpr (waterEnabled) {
+            Sw = priVars.makeEvaluation(water0Idx, timeIdx);
+        }
         Evaluation L = fluidState_.L();
-        Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
+        Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
-        Evaluation Sg = Opm::max(1 - So, 0.0);
+        Evaluation Sg = Opm::max(1 - So - Sw, 0.0);
-        Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg);
+        Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg) + Opm::getValue(Sw);
         So /= sumS;
         Sg /= sumS;
 
-        fluidState_.setSaturation(0, So);
+        fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So);
-        fluidState_.setSaturation(1, Sg);
+        fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg);
+        if constexpr (waterEnabled) {
+            Sw /= sumS;
+            fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw);
+        }
 
-        fluidState_.setCompressFactor(0, Z_L);
+        fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L);
-        fluidState_.setCompressFactor(1, Z_V);
+        fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V);
 
         // Print saturation
         if (flashVerbosity >= 5) {
@@ -250,7 +260,10 @@ public:
 
         // set the phase viscosity and density
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            paramCache.updatePhase(fluidState_, phaseIdx);
+            if (phaseIdx == static_cast<unsigned int>(FluidSystem::oilPhaseIdx) 
+             || phaseIdx == static_cast<unsigned int>(FluidSystem::gasPhaseIdx)) {
+                paramCache.updatePhase(fluidState_, phaseIdx);
+            }
 
             const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
 

opm/models/ptflash/flashlocalresidual.hh

@@ -49,18 +49,24 @@ class FlashLocalResidual: public GetPropType<TypeTag, Properties::DiscLocalResid
     using Indices = GetPropType<TypeTag, Properties::Indices>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
+    enum { water0Idx = Indices::water0Idx };
     enum { conti0EqIdx = Indices::conti0EqIdx };
 
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
     using DiffusionModule = Opm::DiffusionModule<TypeTag, enableDiffusion>;
 
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>;
 
+    static const bool waterEnabled = Indices::waterEnabled;
+
     using Toolbox = Opm::MathToolbox<Evaluation>;
 
 public:
@@ -77,15 +83,25 @@ public:
         const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
         const auto& fs = intQuants.fluidState();
 
-        // compute storage term of all components within all phases
+        // compute water storage term
-        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+        if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
-            unsigned eqIdx = conti0EqIdx + compIdx;
+            unsigned eqIdx = conti0EqIdx + numComponents;
-            storage[eqIdx] +=
+            storage[eqIdx] = 
-            Toolbox::template decay<LhsEval>(fs.massFraction(phaseIdx, compIdx))
+                Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
-                * Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
                 * Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
                 * Toolbox::template decay<LhsEval>(intQuants.porosity());
         }
+        else {
+            // compute storage term of all components within oil/gas phases
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                unsigned eqIdx = conti0EqIdx + compIdx;
+                storage[eqIdx] +=
+                Toolbox::template decay<LhsEval>(fs.massFraction(phaseIdx, compIdx))
+                    * Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
+                    * Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
+                    * Toolbox::template decay<LhsEval>(intQuants.porosity());
+            }
+        }
 
         EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
     }
@@ -146,19 +162,31 @@ public:
                     up.fluidState().density(phaseIdx)
                     * extQuants.volumeFlux(phaseIdx);
 
-                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
-                    flux[conti0EqIdx + compIdx] +=
+                    unsigned eqIdx = conti0EqIdx + numComponents;
-                        tmp*up.fluidState().massFraction(phaseIdx, compIdx);
+                    flux[eqIdx] = tmp;
+                }
+                else {
+                    for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                        flux[conti0EqIdx + compIdx] +=
+                            tmp*up.fluidState().massFraction(phaseIdx, compIdx);
+                    }
                 }
             }
             else {
                 Evaluation tmp =
                     Toolbox::value(up.fluidState().density(phaseIdx))
                      * extQuants.volumeFlux(phaseIdx);
-
+                
-                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
-                    flux[conti0EqIdx + compIdx] +=
+                    unsigned eqIdx = conti0EqIdx + numComponents;
-                        tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx));
+                    flux[eqIdx] = tmp;
+                }
+                else {
+                    for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                        flux[conti0EqIdx + compIdx] +=
+                            tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx));
+                    }
                 }
             }
         }

opm/models/ptflash/flashmodel.hh

@@ -135,6 +135,10 @@ template<class TypeTag>
 struct EnableEnergy<TypeTag, TTag::FlashModel>
 { static constexpr bool value = false; };
 
+template<class TypeTag>
+struct EnableWater<TypeTag, TTag::MultiPhaseBaseModel> 
+{ static constexpr int value = GetPropType<TypeTag, Properties::FluidSystem>::waterEnabled; };
+
 } // namespace Opm::Properties
 
 namespace Opm {

opm/models/ptflash/flashnewtonmethod.hh

@@ -63,6 +63,8 @@ class FlashNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMeth
     enum { z0Idx = Indices::z0Idx };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
 
+    static constexpr bool waterEnabled = Indices::waterEnabled;
+
 public:
     /*!
      * \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& )
@@ -125,6 +127,14 @@ protected:
         for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
             clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol);
         }
+
+        if constexpr (waterEnabled) {
+            // limit change in water saturation to 0.2
+            constexpr Scalar dSwMax = 0.2;
+            if (update[Indices::water0Idx] > dSwMax) {
+                nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax;
+            }
+        }
     }
 private:
     void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const

opm/models/ptflash/flashprimaryvariables.hh

@@ -61,6 +61,14 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     // primary variable indices
     enum { z0Idx = Indices::z0Idx };
     enum { pressure0Idx = Indices::pressure0Idx };
+    enum { water0Idx = Indices::water0Idx };
+
+    static constexpr bool waterEnabled = Indices::waterEnabled;
+
+    // phase indices
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
 
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
@@ -108,10 +116,17 @@ public:
         // the energy module
         EnergyModule::setPriVarTemperatures(*this, fluidState);
 
+        // assign components total fraction
         for (int i = 0; i < numComponents - 1; ++i)
             (*this)[z0Idx + i] = getValue(fluidState.moleFraction(i));
 
-        (*this)[pressure0Idx] = getValue(fluidState.pressure(0));
+        // assign pressure
+        (*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx));
+
+        // assign water saturation
+        if constexpr (waterEnabled) {
+            (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx));
+        }
     }
 
     /*!
@@ -121,12 +136,15 @@ public:
      */
     void print(std::ostream& os = std::cout) const
     {
-        os << "(p_" << FluidSystem::phaseName(0) << " = "
+        os << "(p_" << FluidSystem::phaseName(FluidSystem::oilPhaseIdx) << " = "
            << this->operator[](pressure0Idx);
         for (unsigned compIdx = 0; compIdx < numComponents - 2; ++compIdx) {
             os << ", z_" << FluidSystem::componentName(compIdx) << " = "
                << this->operator[](z0Idx + compIdx);
         }
+        if constexpr (waterEnabled) {
+            os << ", S_w = " << this->operator[](water0Idx);
+        }
         os << ")" << std::flush;
     }
 };

opm/simulators/flow/FIPContainer.cpp

@@ -27,7 +27,7 @@
 
 #include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
-#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
 
 #include <opm/output/data/Solution.hpp>
 
@@ -438,11 +438,19 @@ INSTANTIATE_TYPE(double)
 INSTANTIATE_TYPE(float)
 #endif
 
-#define INSTANTIATE_COMP(NUM) \
+#define INSTANTIATE_COMP_THREEPHASE(NUM) \
-    template<class T> using FS##NUM = GenericOilGasFluidSystem<T, NUM>; \
+    template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM, true>; \
     template class FIPContainer<FS##NUM<double>>;
 
-INSTANTIATE_COMP(0)
+#define INSTANTIATE_COMP_TWOPHASE(NUM) \
+    template<class T> using GFS##NUM = GenericOilGasWaterFluidSystem<T, NUM, false>; \
+    template class FIPContainer<GFS##NUM<double>>;
+
+#define INSTANTIATE_COMP(NUM) \
+    INSTANTIATE_COMP_THREEPHASE(NUM) \
+    INSTANTIATE_COMP_TWOPHASE(NUM)
+
+INSTANTIATE_COMP_THREEPHASE(0)  // \Note: to register the parameter ForceDisableFluidInPlaceOutput
 INSTANTIATE_COMP(2)
 INSTANTIATE_COMP(3)
 INSTANTIATE_COMP(4)

opm/simulators/flow/FlowProblemComp.hpp

@@ -362,6 +362,9 @@ public:
             Dune::FieldVector<Scalar, numComponents> z(0.0);
             Scalar sumMoles = 0.0;
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                if (Indices::waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)){ 
+                    continue;
+                }
                 const auto saturation = fs.saturation(phaseIdx);
                 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                     Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation;
@@ -491,7 +494,7 @@ protected:
         const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
         const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
 
-        if (water_active && Indices::numPhases > 1)
+        if (water_active && Indices::numPhases > 2)
             waterSaturationData = fp.get_double("SWAT");
         else
             waterSaturationData.resize(numDof);
@@ -527,6 +530,10 @@ protected:
                                             - waterSaturationData[dofIdx]
                                             - gasSaturationData[dofIdx]);
             }
+            if (water_active) {
+                dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
+                                            waterSaturationData[dofIdx]);
+            }
 
             //////
             // set phase pressures

opm/simulators/flow/FlowProblemCompProperties.hpp

@@ -69,17 +69,10 @@ private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
-    // using Traits = ThreePhaseMaterialTraits<Scalar,
-    //                                         /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
-    //                                         /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
-    //                                         /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
-
-    // TODO: We should be able to use FluidSystem here and using Indices to handle the active phases
-    // some more development is needed
     using Traits = ThreePhaseMaterialTraits<Scalar,
-                                            /*wettingPhaseIdx=*/ 0,
+                                            /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
-                                            /*nonWettingPhaseIdx=*/ 1,
+                                            /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
-                                            /*gasPhaseIdx=*/ 2>;
+                                            /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
 
 public:
     using EclMaterialLawManager = ::Opm::EclMaterialLawManager<Traits>;

opm/simulators/flow/GenericOutputBlackoilModule.cpp

@@ -30,7 +30,7 @@
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 #include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
 
-#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/input/eclipse/EclipseState/Runspec.hpp>
@@ -1522,12 +1522,19 @@ INSTANTIATE_TYPE(double)
 INSTANTIATE_TYPE(float)
 #endif
 
-#define INSTANTIATE_COMP(NUM) \
+#define INSTANTIATE_COMP_THREEPHASE(NUM) \
-    template<class T> using FS##NUM = GenericOilGasFluidSystem<T, NUM>; \
+    template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM, true>; \
     template class GenericOutputBlackoilModule<FS##NUM<double>>;
 
-INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput
+#define INSTANTIATE_COMP_TWOPHASE(NUM) \
+    template<class T> using GFS##NUM = GenericOilGasWaterFluidSystem<T, NUM, false>; \
+    template class GenericOutputBlackoilModule<GFS##NUM<double>>;
 
+#define INSTANTIATE_COMP(NUM) \
+    INSTANTIATE_COMP_THREEPHASE(NUM) \
+    INSTANTIATE_COMP_TWOPHASE(NUM)
+
+INSTANTIATE_COMP_THREEPHASE(0)  // \Note: to register the parameter ForceDisableFluidInPlaceOutput
 INSTANTIATE_COMP(2)
 INSTANTIATE_COMP(3)
 INSTANTIATE_COMP(4)