Merge branch 'OPM:master' into update-well-state-with-target-modifications

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

CMakeLists_files.cmake

@@ -90,6 +90,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/flow/BlackoilModelParameters.cpp
   opm/simulators/flow/BlackoilModelConvergenceMonitor.cpp
   opm/simulators/flow/CollectDataOnIORank.cpp
+  opm/simulators/flow/CompositionalContainer.cpp
   opm/simulators/flow/ConvergenceOutputConfiguration.cpp
   opm/simulators/flow/EclGenericWriter.cpp
   opm/simulators/flow/ExtboContainer.cpp
@@ -119,6 +120,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/flow/SimulatorReportBanners.cpp
   opm/simulators/flow/SimulatorSerializer.cpp
   opm/simulators/flow/SolutionContainers.cpp
+  opm/simulators/flow/TracerContainer.cpp
   opm/simulators/flow/Transmissibility.cpp
   opm/simulators/flow/ValidationFunctions.cpp
   opm/simulators/flow/equil/EquilibrationHelpers.cpp
@@ -821,6 +823,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/BlackoilModelProperties.hpp
   opm/simulators/flow/CollectDataOnIORank.hpp
   opm/simulators/flow/CollectDataOnIORank_impl.hpp
+  opm/simulators/flow/CompositionalContainer.hpp
   opm/simulators/flow/ConvergenceOutputConfiguration.hpp
   opm/simulators/flow/countGlobalCells.hpp
   opm/simulators/flow/CpGridVanguard.hpp
@@ -878,6 +881,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/SolutionContainers.hpp
   opm/simulators/flow/SubDomain.hpp
   opm/simulators/flow/TTagFlowProblemTPFA.hpp
+  opm/simulators/flow/TracerContainer.hpp
   opm/simulators/flow/TracerModel.hpp
   opm/simulators/flow/Transmissibility.hpp
   opm/simulators/flow/Transmissibility_impl.hpp

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/blackoil/blackoilindices.hh

@@ -94,10 +94,10 @@ struct BlackOilIndices
                                  numEnergy + numFoam + numBrine + numMICPs;
 
     //! \brief returns the index of "active" component
-    static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx)
+    static constexpr int canonicalToActiveComponentIndex(const int compIdx)
     { return compIdx; }
 
-    static constexpr unsigned activeToCanonicalComponentIndex(unsigned compIdx)
+    static constexpr int activeToCanonicalComponentIndex(const int compIdx)
     { return compIdx; }
 
     ////////

opm/models/blackoil/blackoilonephaseindices.hh

@@ -30,6 +30,8 @@
 
 #include <cassert>
 
+#include <opm/common/utility/ConstexprAssert.hpp>
+
 namespace Opm {
 
 /*!
@@ -179,15 +181,15 @@ struct BlackOilOnePhaseIndices
     //////////////////////
 
     //! \brief returns the index of "active" component
-    static constexpr unsigned canonicalToActiveComponentIndex(unsigned /*compIdx*/)
+    static constexpr int canonicalToActiveComponentIndex(const int /*compIdx*/)
     {
         return 0;
     }
 
-    static unsigned activeToCanonicalComponentIndex([[maybe_unused]] unsigned compIdx)
+    static constexpr int activeToCanonicalComponentIndex([[maybe_unused]] const int compIdx)
     {
         // assumes canonical oil = 0, water = 1, gas = 2;
-        assert(compIdx == 0);
+        constexpr_assert(compIdx == 0);
         if (gasEnabled) {
             return 2;
         } else if (waterEnabled) {

opm/models/blackoil/blackoiltwophaseindices.hh

@@ -181,7 +181,7 @@ struct BlackOilTwoPhaseIndices
     //////////////////////
 
     //! \brief returns the index of "active" component
-    static constexpr unsigned canonicalToActiveComponentIndex(const unsigned compIdx)
+    static constexpr int canonicalToActiveComponentIndex(const int compIdx)
     {
         // assumes canonical oil = 0, water = 1, gas = 2;
         if (!gasEnabled) {
@@ -201,10 +201,10 @@ struct BlackOilTwoPhaseIndices
         return compIdx - 1;
     }
 
-    static unsigned activeToCanonicalComponentIndex(unsigned compIdx)
+    static constexpr int activeToCanonicalComponentIndex(const int compIdx)
     {
         // assumes canonical oil = 0, water = 1, gas = 2;
-        assert(compIdx < 2);
+        constexpr_assert(compIdx < 2);
         if (!gasEnabled) {
             // oil = 0, water = 1
             return compIdx;
@@ -212,7 +212,7 @@ struct BlackOilTwoPhaseIndices
             // oil = 0, gas = 1
             return compIdx * 2;
         } else {
-            assert(!oilEnabled);
+            constexpr_assert(!oilEnabled);
         }
 
         // water = 0, gas = 1;

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
 
@@ -67,6 +71,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
 
     //! Index of the mass conservation equation for the first

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) {
+            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,7 +83,16 @@ 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
+        if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
+            unsigned eqIdx = conti0EqIdx + numComponents;
+            storage[eqIdx] = 
+                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] +=
@@ -86,6 +101,7 @@ public:
                     * Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
                     * Toolbox::template decay<LhsEval>(intQuants.porosity());
             }
+        }
 
         EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
     }
@@ -146,22 +162,34 @@ public:
                     up.fluidState().density(phaseIdx)
                     * extQuants.volumeFlux(phaseIdx);
 
+                if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
+                    unsigned eqIdx = conti0EqIdx + numComponents;
+                    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);
                 
+                if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
+                    unsigned eqIdx = conti0EqIdx + numComponents;
+                    flux[eqIdx] = tmp;
+                }
+                else {
                     for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                         flux[conti0EqIdx + compIdx] +=
                             tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx));
                     }
                 }
             }
+        }
 
         EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
     }

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/CompositionalContainer.cpp

@@ -0,0 +1,191 @@
+// -*- 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.
+*/
+
+#include <config.h>
+#include <opm/simulators/flow/CompositionalContainer.hpp>
+
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
+
+#include <opm/output/data/Solution.hpp>
+
+#include <algorithm>
+#include <tuple>
+
+#include <fmt/format.h>
+
+namespace Opm {
+
+template<class FluidSystem>
+void CompositionalContainer<FluidSystem>::
+allocate(const unsigned bufferSize,
+         std::map<std::string, int>& rstKeywords)
+{
+    if (auto& zmf = rstKeywords["ZMF"]; zmf > 0) {
+        this->allocated_ = true;
+        zmf = 0;
+        for (int i = 0; i < numComponents; ++i) {
+            moleFractions_[i].resize(bufferSize, 0.0);
+        }
+    }
+
+    if (auto& xmf = rstKeywords["XMF"]; xmf > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) {
+        this->allocated_ = true;
+        xmf = 0;
+        for (int i = 0; i < numComponents; ++i) {
+            phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0);
+        }
+    }
+
+    if (auto& ymf = rstKeywords["YMF"]; ymf > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
+        this->allocated_ = true;
+        ymf = 0;
+        for (int i = 0; i < numComponents; ++i) {
+            phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0);
+        }
+    }
+}
+
+template<class FluidSystem>
+void CompositionalContainer<FluidSystem>::
+assignGasFractions(const unsigned globalDofIdx,
+                   const AssignFunction& fractions)
+{
+    if (phaseMoleFractions_[gasPhaseIdx][0].empty()) {
+        return;
+    }
+
+    std::for_each(phaseMoleFractions_[gasPhaseIdx].begin(),
+                  phaseMoleFractions_[gasPhaseIdx].end(),
+                  [globalDofIdx, &fractions, c = 0](auto& comp) mutable
+                  { comp[globalDofIdx] = fractions(c++); });
+}
+
+template<class FluidSystem>
+void CompositionalContainer<FluidSystem>::
+assignMoleFractions(const unsigned globalDofIdx,
+                    const AssignFunction& fractions)
+{
+    if (moleFractions_.empty()) {
+        return;
+    }
+
+    std::for_each(moleFractions_.begin(), moleFractions_.end(),
+                  [&fractions, globalDofIdx, c = 0](auto& comp) mutable
+                  { comp[globalDofIdx] = fractions(c++); });
+}
+
+template<class FluidSystem>
+void CompositionalContainer<FluidSystem>::
+assignOilFractions(const unsigned globalDofIdx,
+                   const AssignFunction& fractions)
+{
+    if (phaseMoleFractions_[oilPhaseIdx][0].empty()) {
+        return;
+    }
+
+    std::for_each(phaseMoleFractions_[oilPhaseIdx].begin(),
+                  phaseMoleFractions_[oilPhaseIdx].end(),
+                  [globalDofIdx, &fractions, c = 0](auto& comp) mutable
+                  { comp[globalDofIdx] = fractions(c++); });
+}
+
+template<class FluidSystem>
+void CompositionalContainer<FluidSystem>::
+outputRestart(data::Solution& sol,
+              ScalarBuffer& oil_saturation)
+{
+    using DataEntry =
+        std::tuple<std::string, UnitSystem::measure, std::vector<Scalar>&>;
+
+    auto doInsert = [&sol](DataEntry&       entry,
+                           const data::TargetType target)
+    {
+        if (std::get<2>(entry).empty()) {
+            return;
+        }
+
+        sol.insert(std::get<std::string>(entry),
+                   std::get<UnitSystem::measure>(entry),
+                   std::move(std::get<2>(entry)),
+                   target);
+    };
+
+    auto entries = std::vector<DataEntry>{};
+
+    // ZMF
+    if (!moleFractions_[0].empty()) {
+        for (int i = 0; i < numComponents; ++i) {
+            const auto name = fmt::format("ZMF{}", i + 1);  // Generate ZMF1, ZMF2, ...
+            entries.emplace_back(name, UnitSystem::measure::identity, moleFractions_[i]);
+        }
+    }
+
+    // XMF
+    if (!phaseMoleFractions_[oilPhaseIdx][0].empty()) {
+        for (int i = 0; i < numComponents; ++i) {
+            const auto name = fmt::format("XMF{}", i + 1);  // Generate XMF1, XMF2, ...
+            entries.emplace_back(name, UnitSystem::measure::identity,
+                                 phaseMoleFractions_[oilPhaseIdx][i]);
+        }
+    }
+
+    // YMF
+    if (!phaseMoleFractions_[gasPhaseIdx][0].empty()) {
+        for (int i = 0; i < numComponents; ++i) {
+            const auto name = fmt::format("YMF{}", i + 1);  // Generate YMF1, YMF2, ...
+            entries.emplace_back(name, UnitSystem::measure::identity,
+                                 phaseMoleFractions_[gasPhaseIdx][i]);
+        }
+    }
+
+    if (!oil_saturation.empty()) {
+        entries.emplace_back("SOIL", UnitSystem::measure::identity, oil_saturation);
+    }
+
+    std::for_each(entries.begin(), entries.end(),
+                  [&doInsert](auto& array)
+                  { doInsert(array, data::TargetType::RESTART_SOLUTION); });
+
+    this->allocated_ = false;
+}
+
+#define INSTANTIATE_COMP_THREEPHASE(NUM) \
+    template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM, true>; \
+    template class CompositionalContainer<FS##NUM<double>>;
+
+#define INSTANTIATE_COMP_TWOPHASE(NUM) \
+    template<class T> using GFS##NUM = GenericOilGasWaterFluidSystem<T, NUM, false>; \
+    template class CompositionalContainer<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)
+INSTANTIATE_COMP(5)
+INSTANTIATE_COMP(6)
+INSTANTIATE_COMP(7)
+
+} // namespace Opm

opm/simulators/flow/CompositionalContainer.hpp

@@ -0,0 +1,83 @@
+// -*- 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
+ * \copydoc Opm::OutputBlackOilModule
+ */
+#ifndef OPM_COMPOSITIONAL_CONTAINER_HPP
+#define OPM_COMPOSITIONAL_CONTAINER_HPP
+
+#include <array>
+#include <functional>
+#include <map>
+#include <string>
+#include <vector>
+
+namespace Opm {
+
+namespace data { class Solution; }
+
+template<class FluidSystem>
+class CompositionalContainer
+{
+    using Scalar = typename FluidSystem::Scalar;
+    using ScalarBuffer = std::vector<Scalar>;
+
+    static constexpr int numComponents = FluidSystem::numComponents;
+
+    static constexpr int numPhases = FluidSystem::numPhases;
+    static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
+    static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
+    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
+
+public:
+    void allocate(const unsigned bufferSize,
+                  std::map<std::string, int>& rstKeywords);
+
+    using AssignFunction = std::function<Scalar(const unsigned)>;
+
+    void assignGasFractions(const unsigned globalDofIdx,
+                            const AssignFunction& fractions);
+
+    void assignMoleFractions(const unsigned globalDofIdx,
+                             const AssignFunction& fractions);
+
+    void assignOilFractions(const unsigned globalDofIdx,
+                            const AssignFunction& fractions);
+
+    void outputRestart(data::Solution& sol,
+                       ScalarBuffer& oil_saturation);
+
+    bool allocated() const
+    { return allocated_; }
+
+private:
+    bool allocated_ = false;
+    // total mole fractions for each component
+    std::array<ScalarBuffer, numComponents> moleFractions_;
+    // mole fractions for each component in each phase
+    std::array<std::array<ScalarBuffer, numComponents>, numPhases> phaseMoleFractions_;
+};
+
+} // namespace Opm
+
+#endif // OPM_COMPOSITIONAL_CONTAINER_HPP

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

@@ -27,10 +27,10 @@
 
 #include <opm/grid/common/CommunicationUtils.hpp>
 
+#include <opm/material/fluidmatrixinteractions/EclHysteresisConfig.hpp>
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 #include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
-
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
-#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/input/eclipse/EclipseState/Runspec.hpp>
@@ -55,7 +55,6 @@
 #include <cstddef>
 #include <functional>
 #include <initializer_list>
-#include <stdexcept>
 #include <string>
 #include <string_view>
 #include <tuple>
@@ -104,8 +103,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState,
                             bool enableBrine,
                             bool enableSaltPrecipitation,
                             bool enableExtbo,
-                            bool enableMICP,
+                            bool enableMICP)
-                            bool isCompositional)
     : eclState_(eclState)
     , schedule_(schedule)
     , summaryState_(summaryState)
@@ -124,7 +122,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState,
     , enableSaltPrecipitation_(enableSaltPrecipitation)
     , enableExtbo_(enableExtbo)
     , enableMICP_(enableMICP)
-    , isCompositional_(isCompositional)
+    , tracerC_(eclState_)
     , local_data_valid_(false)
 {
     const auto& fp = eclState_.fieldProps();
@@ -527,38 +525,6 @@ assignToSolution(data::Solution& sol)
         DataEntry{"TMULT_RC", UnitSystem::measure::identity,           rockCompTransMultiplier_},
     };
 
-    // basically, for compositional, we can not use std::array for this.  We need to generate the ZMF1, ZMF2, and so on
-    // and also, we need to map these values.
-    // TODO: the following should go to a function
-    if (this->isCompositional_) {
-        auto compositionalEntries = std::vector<DataEntry>{};
-        {
-            // ZMF
-            for (int i = 0; i < numComponents; ++i) {
-                const auto name = fmt::format("ZMF{}", i + 1);  // Generate ZMF1, ZMF2, ...
-                compositionalEntries.emplace_back(name, UnitSystem::measure::identity, moleFractions_[i]);
-            }
-
-            // XMF
-            for (int i = 0; i < numComponents; ++i) {
-                const auto name = fmt::format("XMF{}", i + 1);  // Generate XMF1, XMF2, ...
-                compositionalEntries.emplace_back(name, UnitSystem::measure::identity,
-                                                  phaseMoleFractions_[oilPhaseIdx][i]);
-            }
-
-            // YMF
-            for (int i = 0; i < numComponents; ++i) {
-                const auto name = fmt::format("YMF{}", i + 1);  // Generate YMF1, YMF2, ...
-                compositionalEntries.emplace_back(name, UnitSystem::measure::identity,
-                                                  phaseMoleFractions_[gasPhaseIdx][i]);
-            }
-        }
-
-        for (auto& array: compositionalEntries) {
-            doInsert(array, data::TargetType::RESTART_SOLUTION);
-        }
-    }
-
     for (auto& array : baseSolutionVector) {
         doInsert(array, data::TargetType::RESTART_SOLUTION);
     }
@@ -602,15 +568,6 @@ assignToSolution(data::Solution& sol)
                    data::TargetType::RESTART_SOLUTION);
     }
 
-    if (this->isCompositional_ && FluidSystem::phaseIsActive(oilPhaseIdx) &&
-        ! this->saturation_[oilPhaseIdx].empty())
-    {
-        sol.insert("SOIL", UnitSystem::measure::identity,
-                   std::move(this->saturation_[oilPhaseIdx]),
-                   data::TargetType::RESTART_SOLUTION);
-    }
-
-
     if ((eclState_.runspec().co2Storage() || eclState_.runspec().h2Storage()) && !rsw_.empty()) {
         auto mfrac = std::vector<double>(this->rsw_.size(), 0.0);
 
@@ -674,41 +631,7 @@ assignToSolution(data::Solution& sol)
     this->fipC_.outputRestart(sol);
 
     // Tracers
-    if (! this->freeTracerConcentrations_.empty()) {
+    this->tracerC_.outputRestart(sol);
-        const auto& tracers = this->eclState_.tracer();
-        for (auto tracerIdx = 0*tracers.size();
-             tracerIdx < tracers.size(); ++tracerIdx)
-        {
-            sol.insert(tracers[tracerIdx].fname(),
-                       UnitSystem::measure::identity,
-                       std::move(freeTracerConcentrations_[tracerIdx]),
-                       data::TargetType::RESTART_TRACER_SOLUTION);
-        }
-
-        // Put freeTracerConcentrations container into a valid state.  Otherwise
-        // we'll move from vectors that have already been moved from if we
-        // get here and it's not a restart step.
-        this->freeTracerConcentrations_.clear();
-    }
-    if (! this->solTracerConcentrations_.empty()) {
-        const auto& tracers = this->eclState_.tracer();
-        for (auto tracerIdx = 0*tracers.size();
-             tracerIdx < tracers.size(); ++tracerIdx)
-        {
-            if (solTracerConcentrations_[tracerIdx].empty())
-                continue;
-
-            sol.insert(tracers[tracerIdx].sname(),
-                       UnitSystem::measure::identity,
-                       std::move(solTracerConcentrations_[tracerIdx]),
-                       data::TargetType::RESTART_TRACER_SOLUTION);
-        }
-
-        // Put solTracerConcentrations container into a valid state.  Otherwise
-        // we'll move from vectors that have already been moved from if we
-        // get here and it's not a restart step.
-        this->solTracerConcentrations_.clear();
-    }
 }
 
 template<class FluidSystem>
@@ -828,16 +751,15 @@ doAllocBuffers(const unsigned bufferSize,
                const bool     substep,
                const bool     log,
                const bool     isRestart,
-               const bool     vapparsActive,
+               const EclHysteresisConfig* hysteresisConfig,
-               const bool     enablePCHysteresis,
+               const unsigned numOutputNnc,
-               const bool     enableNonWettingHysteresis,
+               std::map<std::string, int> rstKeywords)
-               const bool     enableWettingHysteresis,
-               const unsigned numTracers,
-               const std::vector<bool>& enableSolTracers,
-               const unsigned numOutputNnc)
 {
+    if (rstKeywords.empty()) {
+        rstKeywords = schedule_.rst_keywords(reportStepNum);
+    }
+
     // Output RESTART_OPM_EXTENDED only when explicitly requested by user.
-    std::map<std::string, int> rstKeywords = schedule_.rst_keywords(reportStepNum);
     for (auto& [keyword, should_write] : rstKeywords) {
         if (this->isOutputCreationDirective_(keyword)) {
             // 'BASIC', 'FREQ' and similar.  Don't attempt to create
@@ -1027,11 +949,13 @@ doAllocBuffers(const unsigned bufferSize,
         this->micpC_.allocate(bufferSize);
     }
 
+    const bool vapparsActive = schedule_[std::max(reportStepNum, 0u)].oilvap().getType() ==
+                                  OilVaporizationProperties::OilVaporization::VAPPARS;
     if (vapparsActive) {
         soMax_.resize(bufferSize, 0.0);
     }
 
-    if (enableNonWettingHysteresis) {
+    if (hysteresisConfig && hysteresisConfig->enableNonWettingHysteresis()) {
         if (FluidSystem::phaseIsActive(oilPhaseIdx)){
             if (FluidSystem::phaseIsActive(waterPhaseIdx)){
                 soMax_.resize(bufferSize, 0.0);
@@ -1043,7 +967,7 @@ doAllocBuffers(const unsigned bufferSize,
             //TODO add support for gas-water 
         }
     }
-    if (enableWettingHysteresis) {
+    if (hysteresisConfig && hysteresisConfig->enableWettingHysteresis()) {
         if (FluidSystem::phaseIsActive(oilPhaseIdx)){
             if (FluidSystem::phaseIsActive(waterPhaseIdx)){
                 swMax_.resize(bufferSize, 0.0);
@@ -1055,7 +979,7 @@ doAllocBuffers(const unsigned bufferSize,
             //TODO add support for gas-water 
         }
     }
-    if (enablePCHysteresis) {
+    if (hysteresisConfig && hysteresisConfig->enablePCHysteresis()) {
         if (FluidSystem::phaseIsActive(oilPhaseIdx)){
             if (FluidSystem::phaseIsActive(waterPhaseIdx)){
                 swmin_.resize(bufferSize, 0.0);
@@ -1259,19 +1183,7 @@ doAllocBuffers(const unsigned bufferSize,
     }
 
     // tracers
-    if (numTracers > 0) {
+    this->tracerC_.allocate(bufferSize);
-        freeTracerConcentrations_.resize(numTracers);
-        for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
-        {
-            freeTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
-        }
-        solTracerConcentrations_.resize(numTracers);
-        for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
-        {
-            if (enableSolTracers[tracerIdx])
-                solTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
-        }
-    }
 
     if (rstKeywords["RESIDUAL"] > 0) {
         rstKeywords["RESIDUAL"] = 0;
@@ -1293,30 +1205,6 @@ doAllocBuffers(const unsigned bufferSize,
         overburdenPressure_.resize(bufferSize, 0.0);
     }
 
-    if (this->isCompositional_) {
-        if (rstKeywords["ZMF"] > 0) {
-            rstKeywords["ZMF"] = 0;
-            for (int i = 0; i < numComponents; ++i) {
-                moleFractions_[i].resize(bufferSize, 0.0);
-            }
-        }
-
-        if (rstKeywords["XMF"] > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) {
-            rstKeywords["XMF"] = 0;
-            for (int i = 0; i < numComponents; ++i) {
-                phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0);
-            }
-        }
-
-        if (rstKeywords["YMF"] > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
-            rstKeywords["YMF"] = 0;
-            for (int i = 0; i < numComponents; ++i) {
-                phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0);
-            }
-        }
-    }
-
-
     //Warn for any unhandled keyword
     if (log) {
         for (auto& keyValue: rstKeywords) {
@@ -1585,12 +1473,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)

opm/simulators/flow/GenericOutputBlackoilModule.hpp

@@ -40,6 +40,7 @@
 #include <opm/simulators/flow/MechContainer.hpp>
 #include <opm/simulators/flow/MICPContainer.hpp>
 #include <opm/simulators/flow/RegionPhasePVAverage.hpp>
+#include <opm/simulators/flow/TracerContainer.hpp>
 
 #include <opm/simulators/utils/ParallelCommunication.hpp>
 
@@ -62,6 +63,7 @@ struct ForceDisableResvFluidInPlaceOutput { static constexpr bool value = false;
 namespace Opm {
 
 namespace data { class Solution; }
+class EclHysteresisConfig;
 class EclipseState;
 class Schedule;
 class SummaryConfig;
@@ -319,21 +321,16 @@ protected:
                                 bool enableBrine,
                                 bool enableSaltPrecipitation,
                                 bool enableExtbo,
-                                bool enableMICP,
+                                bool enableMICP);
-                                bool isCompositional = false);
 
     void doAllocBuffers(unsigned bufferSize,
                         unsigned reportStepNum,
                         const bool substep,
                         const bool log,
                         const bool isRestart,
-                        const bool vapparsActive = false,
+                        const EclHysteresisConfig* hysteresisConfig,
-                        const bool enablePCHysteresis = false,
+                        unsigned numOutputNnc = 0,
-                        const bool enableNonWettingHysteresis =false,
+                        std::map<std::string, int> rstKeywords = {});
-                        const bool enableWettingHysteresis = false,
-                        unsigned numTracers = 0,
-                        const std::vector<bool>& enableSolTracers = {},
-                        unsigned numOutputNnc = 0);
 
     void makeRegionSum(Inplace& inplace,
                        const std::string& region_name,
@@ -390,7 +387,6 @@ protected:
     bool enableSaltPrecipitation_{false};
     bool enableExtbo_{false};
     bool enableMICP_{false};
-    bool isCompositional_{false};
 
     bool forceDisableFipOutput_{false};
     bool forceDisableFipresvOutput_{false};
@@ -468,12 +464,7 @@ protected:
     std::array<ScalarBuffer, numPhases> viscosity_;
     std::array<ScalarBuffer, numPhases> relativePermeability_;
 
-    // total mole fractions for each component
+    TracerContainer<FluidSystem> tracerC_;
-    std::array<ScalarBuffer, numComponents> moleFractions_;
-    // mole fractions for each component in each phase
-    std::array<std::array<ScalarBuffer, numComponents>, numPhases> phaseMoleFractions_;
-    std::vector<ScalarBuffer> freeTracerConcentrations_;
-    std::vector<ScalarBuffer> solTracerConcentrations_;
 
     std::array<ScalarBuffer, numPhases> residual_;
 

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -184,12 +184,7 @@ public:
                              substep,
                              log,
                              isRestart,
-                             problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)),
+                             &problem.materialLawManager()->hysteresisConfig(),
-                             problem.materialLawManager()->enablePCHysteresis(),
-                             problem.materialLawManager()->enableNonWettingHysteresis(),
-                             problem.materialLawManager()->enableWettingHysteresis(),
-                             problem.tracerModel().numTracers(),
-                             problem.tracerModel().enableSolTracers(),
                              problem.eclWriter()->getOutputNnc().size());
     }
 
@@ -647,25 +642,14 @@ public:
 
             // tracers
             const auto& tracerModel = simulator_.problem().tracerModel();
-            if (! this->freeTracerConcentrations_.empty()) {
+            this->tracerC_.assignFreeConcentrations(globalDofIdx,
-                for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
+                                                    [globalDofIdx, &tracerModel](const unsigned tracerIdx)
-                    if (this->freeTracerConcentrations_[tracerIdx].empty()) {
+                                                    { return tracerModel.freeTracerConcentration(tracerIdx,
-                        continue;
+                                                                                                 globalDofIdx); });
-                    }
+            this->tracerC_.assignSolConcentrations(globalDofIdx,
-                    this->freeTracerConcentrations_[tracerIdx][globalDofIdx] =
+                                                   [globalDofIdx, &tracerModel](const unsigned tracerIdx)
-                        tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx);
+                                                   { return tracerModel.solTracerConcentration(tracerIdx,
-                }
+                                                                                               globalDofIdx); });
-            }
-            if (! this->solTracerConcentrations_.empty()) {
-                for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
-                    if (this->solTracerConcentrations_[tracerIdx].empty()) {
-                        continue;
-                    }
-                    this->solTracerConcentrations_[tracerIdx][globalDofIdx] =
-                        tracerModel.solTracerConcentration(tracerIdx, globalDofIdx);
-                    
-                }
-            }
 
             // output residual
             for ( int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx )

opm/simulators/flow/OutputCompositionalModule.hpp

@@ -32,15 +32,20 @@
 #include <opm/simulators/utils/moduleVersion.hpp>
 
 #include <opm/common/Exceptions.hpp>
+#include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 
 #include <opm/material/common/Valgrind.hpp>
+
+#include <opm/models/blackoil/blackoilproperties.hh>
+#include <opm/models/common/multiphasebaseproperties.hh>
 #include <opm/models/utils/parametersystem.hpp>
 #include <opm/models/utils/propertysystem.hh>
 
+#include <opm/simulators/flow/CompositionalContainer.hpp>
 #include <opm/simulators/flow/FlowBaseVanguard.hpp>
 #include <opm/simulators/flow/GenericOutputBlackoilModule.hpp>
 
@@ -53,8 +58,7 @@
 #include <vector>
 
 
-namespace Opm
+namespace Opm {
-{
 
 // forward declaration
 template <class TypeTag>
@@ -102,8 +106,7 @@ public:
                    getPropValue<TypeTag, Properties::EnableBrine>(),
                    getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
                    getPropValue<TypeTag, Properties::EnableExtbo>(),
-                   getPropValue<TypeTag, Properties::EnableMICP>(),
+                   getPropValue<TypeTag, Properties::EnableMICP>())
-                   true)
         , simulator_(simulator)
     {
         for (auto& region_pair : this->regions_) {
@@ -155,11 +158,19 @@ public:
             return;
         }
 
-        this->doAllocBuffers(bufferSize,
+        auto rstKeywords = this->schedule_.rst_keywords(reportStepNum);
-                              reportStepNum,
+        this->compC_.allocate(bufferSize, rstKeywords);
-                              substep,
+
-                              log,
+        this->doAllocBuffers(bufferSize, reportStepNum, substep, log, isRestart,
-                              isRestart);
+                             /* hysteresisConfig = */ nullptr,
+                             /* numOutputNnc =*/ 0,
+                             std::move(rstKeywords));
+    }
+
+    void assignToSolution(data::Solution& sol)
+    {
+        this->compC_.outputRestart(sol, this->saturation_[oilPhaseIdx]);
+        BaseType::assignToSolution(sol);
     }
 
     /*!
@@ -187,20 +198,21 @@ public:
                 Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]);
             }
 
-            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+            if (this->compC_.allocated()) {
-                if (this->moleFractions_[compIdx].empty()) continue;
+                this->compC_.assignMoleFractions(globalDofIdx,
+                                                 [&fs](const unsigned compIdx)
+                                                 { return getValue(fs.moleFraction(compIdx)); });
 
-                this->moleFractions_[compIdx][globalDofIdx] = getValue(fs.moleFraction(compIdx));
-            }
-            // XMF and YMF
-            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-                if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
-                    if (this->phaseMoleFractions_[oilPhaseIdx][compIdx].empty()) continue;
-                    this->phaseMoleFractions_[oilPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(oilPhaseIdx, compIdx));
-                }
                 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
-                    if (this->phaseMoleFractions_[gasPhaseIdx][compIdx].empty()) continue;
+                    this->compC_.assignGasFractions(globalDofIdx,
-                    this->phaseMoleFractions_[gasPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(gasPhaseIdx, compIdx));
+                                                    [&fs](const unsigned compIdx)
+                                                    { return getValue(fs.moleFraction(gasPhaseIdx, compIdx)); });
+                }
+
+                if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                    this->compC_.assignOilFractions(globalDofIdx,
+                                                    [&fs](const unsigned compIdx)
+                                                    { return getValue(fs.moleFraction(oilPhaseIdx, compIdx)); });
                 }
             }
 
@@ -332,6 +344,7 @@ private:
     }
 
     const Simulator& simulator_;
+    CompositionalContainer<FluidSystem> compC_;
 };
 
 } // namespace Opm

opm/simulators/flow/TracerContainer.cpp

@@ -0,0 +1,154 @@
+// -*- 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.
+*/
+
+#include <config.h>
+#include <opm/simulators/flow/TracerContainer.hpp>
+
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+
+#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
+
+#include <opm/output/data/Solution.hpp>
+
+#include <algorithm>
+#include <utility>
+
+namespace Opm {
+
+template<class FluidSystem>
+void TracerContainer<FluidSystem>::
+allocate(const unsigned bufferSize)
+{
+    const auto& tracers = eclState_.tracer();
+    if (!tracers.empty()) {
+        allocated_ = true;
+        freeConcentrations_.resize(tracers.size());
+        solConcentrations_.resize(tracers.size());
+        std::for_each(tracers.begin(), tracers.end(),
+                      [idx = 0, bufferSize, this](const auto& tracer) mutable
+                      {
+                          freeConcentrations_[idx].resize(bufferSize, 0.0);
+                          if (((tracer.phase == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
+                               (tracer.phase == Phase::OIL && FluidSystem::enableVaporizedOil())) &&
+                              (tracer.solution_concentration.has_value() ||
+                               tracer.solution_tvdp.has_value()))
+                          {
+                              solConcentrations_[idx].resize(bufferSize, 0.0);
+                          }
+                          ++idx;
+                      });
+    }
+}
+
+template<class FluidSystem>
+void TracerContainer<FluidSystem>::
+assignFreeConcentrations(const unsigned globalDofIdx,
+                         const AssignFunction& concentration)
+{
+    std::for_each(freeConcentrations_.begin(), freeConcentrations_.end(),
+                  [globalDofIdx, idx = 0, &concentration](auto& tracer) mutable
+                  {
+                      if (!tracer.empty()) {
+                          tracer[globalDofIdx] = concentration(idx);
+                      }
+                      ++idx;
+                  });
+}
+
+template<class FluidSystem>
+void TracerContainer<FluidSystem>::
+assignSolConcentrations(const unsigned globalDofIdx,
+                         const AssignFunction& concentration)
+{
+    std::for_each(solConcentrations_.begin(), solConcentrations_.end(),
+                  [globalDofIdx, idx = 0, &concentration](auto& tracer) mutable
+                  {
+                      if (!tracer.empty()) {
+                          tracer[globalDofIdx] = concentration(idx);
+                      }
+                      ++idx;
+                  });
+}
+
+template<class FluidSystem>
+void TracerContainer<FluidSystem>::
+outputRestart(data::Solution& sol)
+{
+    if (!this->allocated_) {
+        return;
+    }
+
+    const auto& tracers = this->eclState_.tracer();
+    std::for_each(tracers.begin(), tracers.end(),
+                    [idx = 0, &sol, this](const auto& tracer) mutable
+                    {
+                        sol.insert(tracer.fname(),
+                                   UnitSystem::measure::identity,
+                                   std::move(freeConcentrations_[idx]),
+                                   data::TargetType::RESTART_TRACER_SOLUTION);
+
+                        if (!solConcentrations_[idx].empty()) {
+                            sol.insert(tracer.sname(),
+                                       UnitSystem::measure::identity,
+                                       std::move(solConcentrations_[idx]),
+                                       data::TargetType::RESTART_TRACER_SOLUTION);
+                        }
+                        ++idx;
+                    });
+
+    this->allocated_ = false;
+}
+
+template<class T> using FS = BlackOilFluidSystem<T,BlackOilDefaultIndexTraits>;
+
+#define INSTANTIATE_TYPE(T) \
+    template class TracerContainer<FS<T>>;
+
+INSTANTIATE_TYPE(double)
+
+#if FLOW_INSTANTIATE_FLOAT
+INSTANTIATE_TYPE(float)
+#endif
+
+#define INSTANTIATE_COMP_THREEPHASE(NUM) \
+    template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM, true>; \
+    template class TracerContainer<FS##NUM<double>>;
+
+#define INSTANTIATE_COMP_TWOPHASE(NUM) \
+    template<class T> using GFS##NUM = GenericOilGasWaterFluidSystem<T, NUM, false>; \
+    template class TracerContainer<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)
+INSTANTIATE_COMP(5)
+INSTANTIATE_COMP(6)
+INSTANTIATE_COMP(7)
+
+} // namespace Opm

opm/simulators/flow/TracerContainer.hpp

@@ -0,0 +1,70 @@
+// -*- 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
+ * \copydoc Opm::OutputBlackOilModule
+ */
+#ifndef OPM_TRACER_CONTAINER_HPP
+#define OPM_TRACER_CONTAINER_HPP
+
+#include <functional>
+#include <vector>
+
+namespace Opm {
+
+namespace data { class Solution; }
+class EclipseState;
+
+template<class FluidSystem>
+class TracerContainer
+{
+    using Scalar = typename FluidSystem::Scalar;
+    using ScalarBuffer = std::vector<Scalar>;
+
+public:
+    TracerContainer(const EclipseState& eclState)
+        : eclState_(eclState)
+    {}
+
+    void allocate(const unsigned bufferSize);
+
+    using AssignFunction = std::function<Scalar(const unsigned)>;
+
+    void assignFreeConcentrations(const unsigned globalDofIdx,
+                                  const AssignFunction& concentration);
+
+    void assignSolConcentrations(const unsigned globalDofIdx,
+                                 const AssignFunction& concentration);
+
+    void outputRestart(data::Solution& sol);
+
+private:
+    const EclipseState& eclState_;
+
+    std::vector<ScalarBuffer> freeConcentrations_{};
+    std::vector<ScalarBuffer> solConcentrations_{};
+    bool allocated_{false};
+};
+
+} // namespace Opm
+
+#endif // OPM_TRACER_CONTAINER_HPP

opm/simulators/linalg/ISTLSolver.hpp

@@ -25,6 +25,7 @@
 #include <dune/istl/owneroverlapcopy.hh>
 #include <dune/istl/solver.hh>
 
+#include <opm/common/CriticalError.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/Exceptions.hpp>
 #include <opm/common/TimingMacros.hpp>
@@ -377,10 +378,11 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         void prepare(const Matrix& M, Vector& b)
         {
             OPM_TIMEBLOCK(istlSolverPrepare);
-
+            try {
                 initPrepare(M,b);
 
                 prepareFlexibleSolver();
+            } OPM_CATCH_AND_RETHROW_AS_CRITICAL_ERROR("This is likely due to a faulty linear solver JSON specification. Check for errors related to missing nodes.");
         }
 
 

opm/simulators/linalg/gpuistl/GpuBuffer.cpp

@@ -188,6 +188,7 @@ GpuBuffer<T>::copyFromHost(const std::vector<T>& data)
 {
     copyFromHost(data.data(), data.size());
 }
+
 template <class T>
 void
 GpuBuffer<T>::copyToHost(std::vector<T>& data) const

opm/simulators/linalg/gpuistl/GpuBuffer.hpp

@@ -27,6 +27,7 @@
 #include <opm/simulators/linalg/gpuistl/GpuView.hpp>
 #include <vector>
 #include <string>
+#include <cuda_runtime.h>
 
 
 namespace Opm::gpuistl

opm/simulators/linalg/gpuistl/GpuDILU.cpp

@@ -37,6 +37,7 @@
 #include <functional>
 #include <utility>
 #include <string>
+
 namespace Opm::gpuistl
 {
 
@@ -93,7 +94,8 @@ GpuDILU<M, X, Y, l>::GpuDILU(const M& A, bool splitMatrix, bool tuneKernels, int
         }
     }
 
-    computeDiagAndMoveReorderedData(m_moveThreadBlockSize, m_DILUFactorizationThreadBlockSize);
+    reorderAndSplitMatrix(m_moveThreadBlockSize);
+    computeDiagonal(m_DILUFactorizationThreadBlockSize);
 
     if (m_tuneThreadBlockSizes) {
         tuneThreadBlockSizes();
@@ -110,10 +112,31 @@ template <class M, class X, class Y, int l>
 void
 GpuDILU<M, X, Y, l>::apply(X& v, const Y& d)
 {
+    // ensure that this stream only starts doing work when main stream is completed up to this point
+    OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0));
+    OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0));
+
     OPM_TIMEBLOCK(prec_apply);
     {
+        const auto ptrs = std::make_pair(v.data(), d.data());
+
+        auto it = m_apply_graphs.find(ptrs);
+
+        if (it == m_apply_graphs.end()) {
+            OPM_GPU_SAFE_CALL(cudaStreamBeginCapture(m_stream.get(), cudaStreamCaptureModeGlobal));
+
+            // The apply functions contains lots of small function calls which call a kernel each
             apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize);
+
+            OPM_GPU_SAFE_CALL(cudaStreamEndCapture(m_stream.get(), &m_apply_graphs[ptrs].get()));
+            OPM_GPU_SAFE_CALL(cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0));
         }
+        OPM_GPU_SAFE_CALL(cudaGraphLaunch(m_executableGraphs[ptrs].get(), 0));
+    }
+
+    // ensure that main stream only continues after this stream is completed
+    OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get()));
+    OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0));
 }
 
 template <class M, class X, class Y, int l>
@@ -135,7 +158,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                     m_gpuDInvFloat->data(),
                     d.data(),
                     v.data(),
-                    lowerSolveThreadBlockSize);
+                    lowerSolveThreadBlockSize,
+                    m_stream.get());
             } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) {
                 detail::DILU::solveLowerLevelSetSplit<blocksize_, field_type, float, field_type>(
                     m_gpuMatrixReorderedLowerFloat->getNonZeroValues().data(),
@@ -147,7 +171,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                     m_gpuDInv.data(),
                     d.data(),
                     v.data(),
-                    lowerSolveThreadBlockSize);
+                    lowerSolveThreadBlockSize,
+                    m_stream.get());
             } else {
                 detail::DILU::solveLowerLevelSetSplit<blocksize_, field_type, field_type, field_type>(
                     m_gpuMatrixReorderedLower->getNonZeroValues().data(),
@@ -159,7 +184,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                     m_gpuDInv.data(),
                     d.data(),
                     v.data(),
-                    lowerSolveThreadBlockSize);
+                    lowerSolveThreadBlockSize,
+                    m_stream.get());
             }
         } else {
             detail::DILU::solveLowerLevelSet<field_type, blocksize_>(
@@ -172,7 +198,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                 m_gpuDInv.data(),
                 d.data(),
                 v.data(),
-                lowerSolveThreadBlockSize);
+                lowerSolveThreadBlockSize,
+                m_stream.get());
         }
         levelStartIdx += numOfRowsInLevel;
     }
@@ -193,7 +220,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                     numOfRowsInLevel,
                     m_gpuDInvFloat->data(),
                     v.data(),
-                    upperSolveThreadBlockSize);
+                    upperSolveThreadBlockSize,
+                    m_stream.get());
             } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG){
                 detail::DILU::solveUpperLevelSetSplit<blocksize_, field_type, float>(
                     m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(),
@@ -204,7 +232,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                     numOfRowsInLevel,
                     m_gpuDInv.data(),
                     v.data(),
-                    upperSolveThreadBlockSize);
+                    upperSolveThreadBlockSize,
+                    m_stream.get());
             } else {
                 detail::DILU::solveUpperLevelSetSplit<blocksize_, field_type, field_type>(
                     m_gpuMatrixReorderedUpper->getNonZeroValues().data(),
@@ -215,7 +244,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                     numOfRowsInLevel,
                     m_gpuDInv.data(),
                     v.data(),
-                    upperSolveThreadBlockSize);
+                    upperSolveThreadBlockSize,
+                    m_stream.get());
             }
         } else {
             detail::DILU::solveUpperLevelSet<field_type, blocksize_>(
@@ -227,7 +257,8 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
                 numOfRowsInLevel,
                 m_gpuDInv.data(),
                 v.data(),
-                upperSolveThreadBlockSize);
+                upperSolveThreadBlockSize,
+                m_stream.get());
         }
     }
 }
@@ -252,7 +283,9 @@ GpuDILU<M, X, Y, l>::update()
 {
     OPM_TIMEBLOCK(prec_update);
     {
-        update(m_moveThreadBlockSize, m_DILUFactorizationThreadBlockSize);
+        m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu
+        reorderAndSplitMatrix(m_moveThreadBlockSize);
+        computeDiagonal(m_DILUFactorizationThreadBlockSize);
     }
 }
 
@@ -260,13 +293,22 @@ template <class M, class X, class Y, int l>
 void
 GpuDILU<M, X, Y, l>::update(int moveThreadBlockSize, int factorizationBlockSize)
 {
-    m_gpuMatrix.updateNonzeroValues(m_cpuMatrix, true); // send updated matrix to the gpu
+    // ensure that this stream only starts doing work when main stream is completed up to this point
-    computeDiagAndMoveReorderedData(moveThreadBlockSize, factorizationBlockSize);
+    OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0));
+    OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0));
+
+    m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu
+    reorderAndSplitMatrix(moveThreadBlockSize);
+    computeDiagonal(factorizationBlockSize);
+
+    // ensure that main stream only continues after this stream is completed
+    OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get()));
+    OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0));
 }
 
 template <class M, class X, class Y, int l>
 void
-GpuDILU<M, X, Y, l>::computeDiagAndMoveReorderedData(int moveThreadBlockSize, int factorizationBlockSize)
+GpuDILU<M, X, Y, l>::reorderAndSplitMatrix(int moveThreadBlockSize)
 {
     if (m_splitMatrix) {
         detail::copyMatDataToReorderedSplit<field_type, blocksize_>(
@@ -290,7 +332,12 @@ GpuDILU<M, X, Y, l>::computeDiagAndMoveReorderedData(int moveThreadBlockSize, in
                                                                 m_gpuMatrixReordered->N(),
                                                                 moveThreadBlockSize);
     }
+}
 
+template <class M, class X, class Y, int l>
+void
+GpuDILU<M, X, Y, l>::computeDiagonal(int factorizationBlockSize)
+{
     int levelStartIdx = 0;
     for (int level = 0; level < m_levelSets.size(); ++level) {
         const int numOfRowsInLevel = m_levelSets[level].size();

opm/simulators/linalg/gpuistl/GpuDILU.hpp

@@ -24,8 +24,10 @@
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/kernel_enums.hpp>
+#include <opm/simulators/linalg/gpuistl/gpu_resources.hpp>
 #include <vector>
-
+#include <map>
+#include <utility>
 
 
 namespace Opm::gpuistl
@@ -82,8 +84,11 @@ public:
     //! \brief Updates the matrix data.
     void update() final;
 
+    //! \brief perform matrix splitting and reordering
+    void reorderAndSplitMatrix(int moveThreadBlockSize);
+
     //! \brief Compute the diagonal of the DILU, and update the data of the reordered matrix
-    void computeDiagAndMoveReorderedData(int moveThreadBlockSize, int factorizationThreadBlockSize);
+    void computeDiagonal(int factorizationThreadBlockSize);
 
     //! \brief function that will experimentally tune the thread block sizes of the important cuda kernels
     void tuneThreadBlockSizes();
@@ -153,6 +158,16 @@ private:
     int m_lowerSolveThreadBlockSize = -1;
     int m_moveThreadBlockSize = -1;
     int m_DILUFactorizationThreadBlockSize = -1;
+
+    // Graphs for Apply
+    std::map<std::pair<field_type*, const field_type*>, GPUGraph> m_apply_graphs;
+    std::map<std::pair<field_type*, const field_type*>, GPUGraphExec> m_executableGraphs;
+
+    // Stream for the DILU operations on the GPU
+    GPUStream m_stream{};
+    // Events for synchronization with main stream
+    GPUEvent m_before{};
+    GPUEvent m_after{};
 };
 } // end namespace Opm::gpuistl
 

opm/simulators/linalg/gpuistl/GpuVector.cpp

@@ -266,6 +266,19 @@ GpuVector<T>::copyFromHost(const T* dataPointer, size_t numberOfElements)
     OPM_GPU_SAFE_CALL(cudaMemcpy(data(), dataPointer, numberOfElements * sizeof(T), cudaMemcpyHostToDevice));
 }
 
+template <class T>
+void
+GpuVector<T>::copyFromHost(const T* dataPointer, size_t numberOfElements, cudaStream_t stream)
+{
+    if (numberOfElements > dim()) {
+        OPM_THROW(std::runtime_error,
+                  fmt::format("Requesting to copy too many elements. Vector has {} elements, while {} was requested.",
+                              dim(),
+                              numberOfElements));
+    }
+    OPM_GPU_SAFE_CALL(cudaMemcpyAsync(data(), dataPointer, numberOfElements * sizeof(T), cudaMemcpyHostToDevice, stream));
+}
+
 template <class T>
 void
 GpuVector<T>::copyToHost(T* dataPointer, size_t numberOfElements) const

opm/simulators/linalg/gpuistl/GpuVector.hpp

@@ -203,6 +203,7 @@ public:
      * @note assumes that this vector has numberOfElements elements
      */
     void copyFromHost(const T* dataPointer, size_t numberOfElements);
+    void copyFromHost(const T* dataPointer, size_t numberOfElements, cudaStream_t stream);
 
     /**
      * @brief copyFromHost copies numberOfElements to the CPU memory dataPointer

opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp

@@ -37,6 +37,7 @@
 #include <string>
 #include <tuple>
 #include <utility>
+
 namespace Opm::gpuistl
 {
 
@@ -58,7 +59,6 @@ OpmGpuILU0<M, X, Y, l>::OpmGpuILU0(const M& A, bool splitMatrix, bool tuneKernel
     , m_tuneThreadBlockSizes(tuneKernels)
     , m_mixedPrecisionScheme(makeMatrixStorageMPScheme(mixedPrecisionScheme))
 {
-
     // TODO: Should in some way verify that this matrix is symmetric, only do it debug mode?
     // Some sanity check
     OPM_ERROR_IF(A.N() != m_gpuMatrix.N(),
@@ -98,7 +98,8 @@ OpmGpuILU0<M, X, Y, l>::OpmGpuILU0(const M& A, bool splitMatrix, bool tuneKernel
         }
     }
 
-    LUFactorizeAndMoveData(m_moveThreadBlockSize, m_ILU0FactorizationThreadBlockSize);
+    reorderAndSplitMatrix(m_moveThreadBlockSize);
+    LUFactorizeMatrix(m_ILU0FactorizationThreadBlockSize);
 
     if (m_tuneThreadBlockSizes) {
         tuneThreadBlockSizes();
@@ -117,7 +118,29 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d)
 {
     OPM_TIMEBLOCK(prec_apply);
     {
+        // ensure that this stream only starts doing work when main stream is completed up to this point
+        OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0));
+        OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0));
+
+        const auto ptrs = std::make_pair(v.data(), d.data());
+
+        auto it = m_apply_graphs.find(ptrs);
+
+        if (it == m_apply_graphs.end()) {
+            OPM_GPU_SAFE_CALL(cudaStreamBeginCapture(m_stream.get(), cudaStreamCaptureModeGlobal));
+
+            // The apply functions contains lots of small function calls which call a kernel each
             apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize);
+
+            OPM_GPU_SAFE_CALL(cudaStreamEndCapture(m_stream.get(), &m_apply_graphs[ptrs].get()));
+            OPM_GPU_SAFE_CALL(cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0));
+        }
+        OPM_GPU_SAFE_CALL(cudaGraphLaunch(m_executableGraphs[ptrs].get(), 0));
+
+
+        // ensure that main stream only continues after this stream is completed
+        OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get()));
+        OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0));
     }
 }
 
@@ -142,7 +165,8 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i
                     numOfRowsInLevel,
                     d.data(),
                     v.data(),
-                    lowerSolveThreadBlockSize);
+                    lowerSolveThreadBlockSize,
+                    m_stream.get());
             }
             else{
                 detail::ILU0::solveLowerLevelSetSplit<blocksize_, field_type, field_type>(
@@ -154,7 +178,8 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i
                     numOfRowsInLevel,
                     d.data(),
                     v.data(),
-                    lowerSolveThreadBlockSize);
+                    lowerSolveThreadBlockSize,
+                    m_stream.get());
             }
         } else {
             detail::ILU0::solveLowerLevelSet<field_type, blocksize_>(m_gpuReorderedLU->getNonZeroValues().data(),
@@ -165,7 +190,8 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i
                                                                      numOfRowsInLevel,
                                                                      d.data(),
                                                                      v.data(),
-                                                                     lowerSolveThreadBlockSize);
+                                                                     lowerSolveThreadBlockSize,
+                                                                     m_stream.get());
         }
         levelStartIdx += numOfRowsInLevel;
     }
@@ -185,7 +211,8 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i
                         numOfRowsInLevel,
                         m_gpuMatrixReorderedDiagFloat.value().data(),
                         v.data(),
-                        upperSolveThreadBlockSize);
+                        upperSolveThreadBlockSize,
+                        m_stream.get());
             }
             else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) {
                     detail::ILU0::solveUpperLevelSetSplit<blocksize_, field_type, float, field_type>(
@@ -197,7 +224,8 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i
                         numOfRowsInLevel,
                         m_gpuMatrixReorderedDiag.value().data(),
                         v.data(),
-                        upperSolveThreadBlockSize);
+                        upperSolveThreadBlockSize,
+                        m_stream.get());
             }
             else{
                 detail::ILU0::solveUpperLevelSetSplit<blocksize_, field_type, field_type, field_type>(
@@ -209,7 +237,8 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i
                     numOfRowsInLevel,
                     m_gpuMatrixReorderedDiag.value().data(),
                     v.data(),
-                    upperSolveThreadBlockSize);
+                    upperSolveThreadBlockSize,
+                    m_stream.get());
             }
         } else {
             detail::ILU0::solveUpperLevelSet<field_type, blocksize_>(m_gpuReorderedLU->getNonZeroValues().data(),
@@ -219,7 +248,8 @@ OpmGpuILU0<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i
                                                                      levelStartIdx,
                                                                      numOfRowsInLevel,
                                                                      v.data(),
-                                                                     upperSolveThreadBlockSize);
+                                                                     upperSolveThreadBlockSize,
+                                                                     m_stream.get());
         }
     }
 }
@@ -241,10 +271,9 @@ template <class M, class X, class Y, int l>
 void
 OpmGpuILU0<M, X, Y, l>::update()
 {
-    OPM_TIMEBLOCK(prec_update);
+    m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu
-    {
+    reorderAndSplitMatrix(m_moveThreadBlockSize);
-        update(m_moveThreadBlockSize, m_ILU0FactorizationThreadBlockSize);
+    LUFactorizeMatrix(m_ILU0FactorizationThreadBlockSize);
-    }
 }
 
 template <class M, class X, class Y, int l>
@@ -253,13 +282,23 @@ OpmGpuILU0<M, X, Y, l>::update(int moveThreadBlockSize, int factorizationThreadB
 {
     OPM_TIMEBLOCK(prec_update);
     {
+        // ensure that this stream only starts doing work when main stream is completed up to this point
+        OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0));
+        OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0));
+
         m_gpuMatrix.updateNonzeroValues(m_cpuMatrix, true); // send updated matrix to the gpu
-        LUFactorizeAndMoveData(moveThreadBlockSize, factorizationThreadBlockSize);
+        reorderAndSplitMatrix(moveThreadBlockSize);
+        LUFactorizeMatrix(factorizationThreadBlockSize);
+
+        // ensure that main stream only continues after this stream is completed
+        OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get()));
+        OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0));
     }
 }
+
 template <class M, class X, class Y, int l>
 void
-OpmGpuILU0<M, X, Y, l>::LUFactorizeAndMoveData(int moveThreadBlockSize, int factorizationThreadBlockSize)
+OpmGpuILU0<M, X, Y, l>::reorderAndSplitMatrix(int moveThreadBlockSize)
 {
     if (m_splitMatrix) {
         detail::copyMatDataToReorderedSplit<field_type, blocksize_>(
@@ -283,6 +322,12 @@ OpmGpuILU0<M, X, Y, l>::LUFactorizeAndMoveData(int moveThreadBlockSize, int fact
                                                                m_gpuReorderedLU->N(),
                                                                moveThreadBlockSize);
     }
+}
+
+template <class M, class X, class Y, int l>
+void
+OpmGpuILU0<M, X, Y, l>::LUFactorizeMatrix(int factorizationThreadBlockSize)
+{
     int levelStartIdx = 0;
     for (int level = 0; level < m_levelSets.size(); ++level) {
         const int numOfRowsInLevel = m_levelSets[level].size();

opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp

@@ -24,6 +24,7 @@
 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
 #include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
+#include <opm/simulators/linalg/gpuistl/gpu_resources.hpp>
 #include <opm/simulators/linalg/gpuistl/detail/kernel_enums.hpp>
 #include <optional>
 #include <type_traits>
@@ -84,8 +85,11 @@ public:
     //! \brief Updates the matrix data.
     void update() final;
 
+    //! \brief perform matrix splitting and reordering
+    void reorderAndSplitMatrix(int moveThreadBlockSize);
+
     //! \brief Compute LU factorization, and update the data of the reordered matrix
-    void LUFactorizeAndMoveData(int moveThreadBlockSize, int factorizationThreadBlockSize);
+    void LUFactorizeMatrix(int factorizationThreadBlockSize);
 
     //! \brief function that will experimentally tune the thread block sizes of the important cuda kernels
     void tuneThreadBlockSizes();
@@ -152,6 +156,16 @@ private:
     int m_lowerSolveThreadBlockSize = -1;
     int m_moveThreadBlockSize = -1;
     int m_ILU0FactorizationThreadBlockSize = -1;
+
+    // Graphs for Apply
+    std::map<std::pair<field_type*, const field_type*>, GPUGraph> m_apply_graphs;
+    std::map<std::pair<field_type*, const field_type*>, GPUGraphExec> m_executableGraphs;
+
+    // Stream for the DILU operations on the GPU
+    GPUStream m_stream{};
+    // Events for synchronization with main stream
+    GPUEvent m_before{};
+    GPUEvent m_after{};
 };
 } // end namespace Opm::gpuistl
 

opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu

@@ -85,10 +85,12 @@ namespace
             // TODO: removce the first condition in the for loop
             for (int block = nnzIdx; block < nnzIdxLim; ++block) {
                 const int col = colIndices[block];
-                mmvMixedGeneral<blocksize, MatrixScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(&mat[block * blocksize * blocksize], &v[col * blocksize], rhs);
+                mmvMixedGeneral<blocksize, MatrixScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(
+                    &mat[block * blocksize * blocksize], &v[col * blocksize], rhs);
             }
 
-            mvMixedGeneral<blocksize, DiagonalScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(&dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]);
+            mvMixedGeneral<blocksize, DiagonalScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(
+                &dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]);
         }
     }
 
@@ -137,10 +139,12 @@ namespace
             LinearSolverScalar rhs[blocksize] = {0};
             for (int block = nnzIdx; block < nnzIdxLim; ++block) {
                 const int col = colIndices[block];
-                umvMixedGeneral<blocksize, MatrixScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(&mat[block * blocksize * blocksize], &v[col * blocksize], rhs);
+                umvMixedGeneral<blocksize, MatrixScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(
+                    &mat[block * blocksize * blocksize], &v[col * blocksize], rhs);
             }
 
-            mmvMixedGeneral<blocksize, DiagonalScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(&dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]);
+            mmvMixedGeneral<blocksize, DiagonalScalar, LinearSolverScalar, LinearSolverScalar, LinearSolverScalar>(
+                &dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]);
         }
     }
 
@@ -211,8 +215,8 @@ namespace
         }
     }
 
-    // TODO: rewrite such that during the factorization there is a dInv of InputScalar type that stores intermediate results
+    // TODO: rewrite such that during the factorization there is a dInv of InputScalar type that stores intermediate
-    // TOOD: The important part is to only cast after that is fully computed
+    // results TOOD: The important part is to only cast after that is fully computed
     template <int blocksize, class InputScalar, class OutputScalar, MatrixStorageMPScheme mixedPrecisionScheme>
     __global__ void cuComputeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat,
                                                int* lowerRowIndices,
@@ -239,7 +243,8 @@ namespace
             InputScalar dInvTmp[blocksize * blocksize];
             for (int i = 0; i < blocksize; ++i) {
                 for (int j = 0; j < blocksize; ++j) {
-                    dInvTmp[i * blocksize + j] = srcDiagonal[reorderedRowIdx * blocksize * blocksize + i * blocksize + j];
+                    dInvTmp[i * blocksize + j]
+                        = srcDiagonal[reorderedRowIdx * blocksize * blocksize + i * blocksize + j];
                 }
             }
 
@@ -257,8 +262,12 @@ namespace
 
                 if constexpr (detail::storeOffDiagonalAsFloat(mixedPrecisionScheme)) {
                     // TODO: think long and hard about whether this performs only the wanted memory transfers
-                    moveBlock<blocksize, InputScalar, OutputScalar>(&srcReorderedLowerMat[block * blocksize * blocksize], &dstLowerMat[block * blocksize * blocksize]);
+                    moveBlock<blocksize, InputScalar, OutputScalar>(
-                    moveBlock<blocksize, InputScalar, OutputScalar>(&srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], &dstUpperMat[symOppositeBlock * blocksize * blocksize]);
+                        &srcReorderedLowerMat[block * blocksize * blocksize],
+                        &dstLowerMat[block * blocksize * blocksize]);
+                    moveBlock<blocksize, InputScalar, OutputScalar>(
+                        &srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize],
+                        &dstUpperMat[symOppositeBlock * blocksize * blocksize]);
                 }
 
                 mmx2Subtraction<InputScalar, blocksize>(&srcReorderedLowerMat[block * blocksize * blocksize],
@@ -271,7 +280,8 @@ namespace
             moveBlock<blocksize, InputScalar, InputScalar>(dInvTmp, &dInv[reorderedRowIdx * blocksize * blocksize]);
 
             if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)) {
-                moveBlock<blocksize, InputScalar, OutputScalar>(dInvTmp, &dstDiag[reorderedRowIdx * blocksize * blocksize]); // important!
+                moveBlock<blocksize, InputScalar, OutputScalar>(
+                    dInvTmp, &dstDiag[reorderedRowIdx * blocksize * blocksize]); // important!
             }
         }
     }
@@ -289,12 +299,13 @@ solveLowerLevelSet(T* reorderedMat,
                    const T* dInv,
                    const T* d,
                    T* v,
-                   int thrBlockSize)
+                   int thrBlockSize,
+                   cudaStream_t stream)
 {
     int threadBlockSize
         = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveLowerLevelSet<T, blocksize>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveLowerLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveLowerLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize, 0, stream>>>(
         reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v);
 }
 
@@ -310,12 +321,14 @@ solveLowerLevelSetSplit(MatrixScalar* reorderedMat,
                         const DiagonalScalar* dInv,
                         const LinearSolverScalar* d,
                         LinearSolverScalar* v,
-                        int thrBlockSize)
+                        int thrBlockSize,
+                        cudaStream_t stream)
 {
     int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(
         cuSolveLowerLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveLowerLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveLowerLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>
+        <<<nThreadBlocks, threadBlockSize, 0, stream>>>(
             reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v);
 }
 // perform the upper solve for all rows in the same level set
@@ -329,12 +342,13 @@ solveUpperLevelSet(T* reorderedMat,
                    int rowsInLevelSet,
                    const T* dInv,
                    T* v,
-                   int thrBlockSize)
+                   int thrBlockSize,
+                   cudaStream_t stream)
 {
     int threadBlockSize
         = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveUpperLevelSet<T, blocksize>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveUpperLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveUpperLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize, 0, stream>>>(
         reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v);
 }
 
@@ -348,12 +362,14 @@ solveUpperLevelSetSplit(MatrixScalar* reorderedMat,
                         int rowsInLevelSet,
                         const DiagonalScalar* dInv,
                         LinearSolverScalar* v,
-                        int thrBlockSize)
+                        int thrBlockSize,
+                        cudaStream_t stream)
 {
     int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(
         cuSolveUpperLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveUpperLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveUpperLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>
+        <<<nThreadBlocks, threadBlockSize, 0, stream>>>(
             reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v);
 }
 
@@ -409,7 +425,8 @@ computeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat,
         int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(
             cuComputeDiluDiagonalSplit<blocksize, InputScalar, OutputScalar, scheme>, thrBlockSize);
         int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-        cuComputeDiluDiagonalSplit<blocksize, InputScalar, OutputScalar, scheme><<<nThreadBlocks, threadBlockSize>>>(srcReorderedLowerMat,
+        cuComputeDiluDiagonalSplit<blocksize, InputScalar, OutputScalar, scheme>
+            <<<nThreadBlocks, threadBlockSize>>>(srcReorderedLowerMat,
                                                  lowerRowIndices,
                                                  lowerColIndices,
                                                  srcReorderedUpperMat,
@@ -429,23 +446,113 @@ computeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat,
     }
 }
 
-// TODO: format
 #define INSTANTIATE_KERNEL_WRAPPERS(T, blocksize)                                                                      \
     template void computeDiluDiagonal<T, blocksize>(T*, int*, int*, int*, int*, const int, int, T*, int);              \
     template void computeDiluDiagonalSplit<blocksize, T, double, MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG>(   \
-        const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int);                                      \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const int,                                                                                                     \
+        int,                                                                                                           \
+        T*,                                                                                                            \
+        double*,                                                                                                       \
+        double*,                                                                                                       \
+        double*,                                                                                                       \
+        int);                                                                                                          \
     template void computeDiluDiagonalSplit<blocksize, T, float, MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG>(    \
-        const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int);                                      \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const int,                                                                                                     \
+        int,                                                                                                           \
+        T*,                                                                                                            \
+        float*,                                                                                                        \
+        float*,                                                                                                        \
+        float*,                                                                                                        \
+        int);                                                                                                          \
     template void computeDiluDiagonalSplit<blocksize, T, float, MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG>(      \
-        const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int);                                      \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const int,                                                                                                     \
+        int,                                                                                                           \
+        T*,                                                                                                            \
+        float*,                                                                                                        \
+        float*,                                                                                                        \
+        float*,                                                                                                        \
+        int);                                                                                                          \
     template void computeDiluDiagonalSplit<blocksize, T, double, MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG>(     \
-        const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int);                                      \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const int,                                                                                                     \
+        int,                                                                                                           \
+        T*,                                                                                                            \
+        double*,                                                                                                       \
+        double*,                                                                                                       \
+        double*,                                                                                                       \
+        int);                                                                                                          \
     template void computeDiluDiagonalSplit<blocksize, T, float, MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG>(     \
-        const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int);                                      \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const int,                                                                                                     \
+        int,                                                                                                           \
+        T*,                                                                                                            \
+        float*,                                                                                                        \
+        float*,                                                                                                        \
+        float*,                                                                                                        \
+        int);                                                                                                          \
     template void computeDiluDiagonalSplit<blocksize, T, double, MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG>(    \
-        const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int);                                      \
+        const T*,                                                                                                      \
-    template void solveUpperLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, const T*, T*, int);                 \
+        int*,                                                                                                          \
-    template void solveLowerLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, const T*, const T*, T*, int);
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const T*,                                                                                                      \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        const int,                                                                                                     \
+        int,                                                                                                           \
+        T*,                                                                                                            \
+        double*,                                                                                                       \
+        double*,                                                                                                       \
+        double*,                                                                                                       \
+        int);                                                                                                          \
+    template void solveUpperLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t);   \
+    template void solveLowerLevelSet<T, blocksize>(                                                                    \
+        T*, int*, int*, int*, int, int, const T*, const T*, T*, int, cudaStream_t);
 
 INSTANTIATE_KERNEL_WRAPPERS(float, 1);
 INSTANTIATE_KERNEL_WRAPPERS(float, 2);
@@ -462,9 +569,19 @@ INSTANTIATE_KERNEL_WRAPPERS(double, 6);
 
 #define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar)                 \
     template void solveUpperLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>(                \
-        MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, LinearSolverScalar*, int);                     \
+        MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, LinearSolverScalar*, int, cudaStream_t);     \
     template void solveLowerLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>(                \
-        MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, const LinearSolverScalar*, LinearSolverScalar*, int);
+        MatrixScalar*,                                                                                                 \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        int*,                                                                                                          \
+        int,                                                                                                           \
+        int,                                                                                                           \
+        const DiagonalScalar*,                                                                                         \
+        const LinearSolverScalar*,                                                                                     \
+        LinearSolverScalar*,                                                                                           \
+        int,                                                                                                           \
+        cudaStream_t);
 
 // TODO: be smarter about this... Surely this instantiates many more combinations that are actually needed
 #define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT_ALL(blocksize)                                                               \

opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp

@@ -54,7 +54,8 @@ void solveLowerLevelSet(T* reorderedMat,
                         const T* dInv,
                         const T* d,
                         T* v,
-                        int threadBlockSize);
+                        int threadBlockSize,
+                        cudaStream_t stream);
 
 /**
  * @brief Perform a lower solve on certain rows in a matrix that can safely be computed in parallel
@@ -82,7 +83,8 @@ void solveLowerLevelSetSplit(MatrixScalar* reorderedUpperMat,
                              const DiagonalScalar* dInv,
                              const LinearSolverScalar* d,
                              LinearSolverScalar* v,
-                             int threadBlockSize);
+                             int threadBlockSize,
+                             cudaStream_t stream);
 
 /**
  * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel
@@ -108,7 +110,8 @@ void solveUpperLevelSet(T* reorderedMat,
                         int rowsInLevelSet,
                         const T* dInv,
                         T* v,
-                        int threadBlockSize);
+                        int threadBlockSize,
+                        cudaStream_t stream);
 
 /**
  * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel
@@ -134,7 +137,8 @@ void solveUpperLevelSetSplit(MatrixScalar* reorderedUpperMat,
                              int rowsInLevelSet,
                              const DiagonalScalar* dInv,
                              LinearSolverScalar* v,
-                             int threadBlockSize);
+                             int threadBlockSize,
+                             cudaStream_t stream);
 
 /**
  * @brief Computes the ILU0 of the diagonal elements of the reordered matrix and stores it in a reordered vector

opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu

@@ -215,9 +215,9 @@ namespace
             // as of now, if we are using mixed precision, then we are always storing the off-diagonals as floats,
             // and sometimes also the diagonal.
             if constexpr (detail::usingMixedPrecision(mixedPrecisionScheme)) {
-                // if we are want to store the entire matrix as a float then we must also move the diagonal block from double to float
+                // if we are want to store the entire matrix as a float then we must also move the diagonal block from
-                // if not then we just use the double diagonal that is already now stored in srcDiagonal
+                // double to float if not then we just use the double diagonal that is already now stored in srcDiagonal
-                if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)){
+                if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)) {
                     moveBlock<blocksize, InputScalar, OutputScalar>(&srcDiagonal[reorderedIdx * scalarsInBlock],
                                                                     &dstDiagonal[reorderedIdx * scalarsInBlock]);
                 }
@@ -362,12 +362,13 @@ solveLowerLevelSet(T* reorderedMat,
                    int rowsInLevelSet,
                    const T* d,
                    T* v,
-                   int thrBlockSize)
+                   int thrBlockSize,
+                   cudaStream_t stream)
 {
     int threadBlockSize
         = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveLowerLevelSet<T, blocksize>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveLowerLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveLowerLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize, 0, stream>>>(
         reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v);
 }
 // perform the upper solve for all rows in the same level set
@@ -380,12 +381,13 @@ solveUpperLevelSet(T* reorderedMat,
                    int startIdx,
                    int rowsInLevelSet,
                    T* v,
-                   int thrBlockSize)
+                   int thrBlockSize,
+                   cudaStream_t stream)
 {
     int threadBlockSize
         = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveUpperLevelSet<T, blocksize>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveUpperLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveUpperLevelSet<T, blocksize><<<nThreadBlocks, threadBlockSize, 0, stream>>>(
         reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, v);
 }
 
@@ -399,12 +401,14 @@ solveLowerLevelSetSplit(MatrixScalar* reorderedMat,
                         int rowsInLevelSet,
                         const LinearSolverScalar* d,
                         LinearSolverScalar* v,
-                        int thrBlockSize)
+                        int thrBlockSize,
+                        cudaStream_t stream)
 {
     int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(
         cuSolveLowerLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveLowerLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveLowerLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar>
+        <<<nThreadBlocks, threadBlockSize, 0, stream>>>(
             reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v);
 }
 // perform the upper solve for all rows in the same level set
@@ -418,12 +422,14 @@ solveUpperLevelSetSplit(MatrixScalar* reorderedMat,
                         int rowsInLevelSet,
                         const DiagonalScalar* dInv,
                         LinearSolverScalar* v,
-                        int thrBlockSize)
+                        int thrBlockSize,
+                        cudaStream_t stream)
 {
     int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(
         cuSolveUpperLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>, thrBlockSize);
     int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
-    cuSolveUpperLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar><<<nThreadBlocks, threadBlockSize>>>(
+    cuSolveUpperLevelSetSplit<blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar>
+        <<<nThreadBlocks, threadBlockSize, 0, stream>>>(
             reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v);
 }
 
@@ -484,8 +490,8 @@ LUFactorizationSplit(InputScalar* srcReorderedLowerMat,
 }
 
 #define INSTANTIATE_KERNEL_WRAPPERS(T, blocksize)                                                                      \
-    template void solveUpperLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, T*, int);                           \
+    template void solveUpperLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, T*, int, cudaStream_t);             \
-    template void solveLowerLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, const T*, T*, int);                 \
+    template void solveLowerLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t);   \
     template void LUFactorization<T, blocksize>(T*, int*, int*, int*, int*, size_t, int, int);                         \
     template void LUFactorizationSplit<blocksize, T, float, MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG>(        \
         T*, int*, int*, T*, int*, int*, T*, float*, float*, float*, int*, int*, const int, int, int);                  \
@@ -514,32 +520,32 @@ INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(double)
 #define INSTANTIATE_MIXED_PRECISION_KERNEL_WRAPPERS(blocksize)                                                         \
     /* double preconditioner */                                                                                        \
     template void solveLowerLevelSetSplit<blocksize, double, double>(                                                  \
-        double*, int*, int*, int*, int, int, const double*, double*, int);                                             \
+        double*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t);                               \
     /* float matrix, double compute preconditioner */                                                                  \
     template void solveLowerLevelSetSplit<blocksize, double, float>(                                                   \
-        float*, int*, int*, int*, int, int, const double*, double*, int);                                              \
+        float*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t);                                \
     /* float preconditioner */                                                                                         \
     template void solveLowerLevelSetSplit<blocksize, float, float>(                                                    \
-        float*, int*, int*, int*, int, int, const float*, float*, int);                                                \
+        float*, int*, int*, int*, int, int, const float*, float*, int, cudaStream_t);                                  \
                                                                                                                        \
     /* double preconditioner */                                                                                        \
     template void solveUpperLevelSetSplit<blocksize, double, double, double>(                                          \
-        double*, int*, int*, int*, int, int, const double*, double*, int);                                             \
+        double*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t);                               \
     /* float matrix, double compute preconditioner */                                                                  \
     template void solveUpperLevelSetSplit<blocksize, double, float, double>(                                           \
-        float*, int*, int*, int*, int, int, const double*, double*, int);                                               \
+        float*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t);                                \
     /* float preconditioner */                                                                                         \
     template void solveUpperLevelSetSplit<blocksize, float, float, double>(                                            \
-        float*, int*, int*, int*, int, int, const double*, float*, int);                                                \
+        float*, int*, int*, int*, int, int, const double*, float*, int, cudaStream_t);                                 \
     /* double preconditioner */                                                                                        \
     template void solveUpperLevelSetSplit<blocksize, double, double, float>(                                           \
-        double*, int*, int*, int*, int, int, const float*, double*, int);                                             \
+        double*, int*, int*, int*, int, int, const float*, double*, int, cudaStream_t);                                \
     /* float matrix, double compute preconditioner */                                                                  \
     template void solveUpperLevelSetSplit<blocksize, double, float, float>(                                            \
-        float*, int*, int*, int*, int, int, const float*, double*, int);                                               \
+        float*, int*, int*, int*, int, int, const float*, double*, int, cudaStream_t);                                 \
     /* float preconditioner */                                                                                         \
     template void solveUpperLevelSetSplit<blocksize, float, float, float>(                                             \
-        float*, int*, int*, int*, int, int, const float*, float*, int);
+        float*, int*, int*, int*, int, int, const float*, float*, int, cudaStream_t);
 
 
 INSTANTIATE_MIXED_PRECISION_KERNEL_WRAPPERS(1);

opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp

@@ -46,7 +46,8 @@ void solveUpperLevelSet(T* reorderedMat,
                         int startIdx,
                         int rowsInLevelSet,
                         T* v,
-                        int threadBlockSize);
+                        int threadBlockSize,
+                        cudaStream_t stream);
 
 /**
  * @brief Perform a lower solve on certain rows in a matrix that can safely be computed in parallel
@@ -72,7 +73,8 @@ void solveLowerLevelSet(T* reorderedMat,
                         int rowsInLevelSet,
                         const T* d,
                         T* v,
-                        int threadBlockSize);
+                        int threadBlockSize,
+                        cudaStream_t stream);
 
 /**
  * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel
@@ -99,7 +101,8 @@ void solveUpperLevelSetSplit(MatrixScalar* reorderedMat,
                              int rowsInLevelSet,
                              const DiagonalScalar* dInv,
                              LinearSolverScalar* v,
-                             int threadBlockSize);
+                             int threadBlockSize,
+                             cudaStream_t stream);
 
 /**
  * @brief Perform an lower solve on certain rows in a matrix that can safely be computed in parallel
@@ -127,7 +130,8 @@ void solveLowerLevelSetSplit(MatrixScalar* reorderedLowerMat,
                              int rowsInLevelSet,
                              const LinearSolverScalar* d,
                              LinearSolverScalar* v,
-                             int threadBlockSize);
+                             int threadBlockSize,
+                             cudaStream_t stream);
 
 /**
  * @brief Computes the ILU Factorization of the input bcsr matrix, which is stored in a reordered way. The diagonal

opm/simulators/utils/FullySupportedFlowKeywords.cpp

@@ -32,6 +32,12 @@ const SupportedKeywordItems<std::string>&
 fullySupported()
 {
    static const SupportedKeywordItems<std::string> fully_supported_keywords_strings = {
+         {
+            "GEFAC",
+            {
+               {3,{true, is_bool_convertible {}, "GEFAC(GRPNETWK): String value must be convertible to bool."}}, // USE_GEFAC_IN_NETWORK
+            },
+         },
          {
             "NEXTSTEP",
             {
@@ -44,6 +50,12 @@ fullySupported()
                {3,{true, allow_values<std::string> {"ORAT", "WRAT", "GRAT", "LRAT", "RESV", "BHP"}, "WCONHIST(TARGET): should be set to ORAT/WRAT/GRAT/LRAT/RESV or BHP"}}, // CMODE
             },
          },
+         {
+            "WEFAC",
+            {
+               {3,{true, is_bool_convertible {}, "WEFAC(WELNETWK): String value must be convertible to bool."}}, // USE_WEFAC_IN_NETWORK
+            },
+         },
    };
 
    return fully_supported_keywords_strings;
@@ -66,6 +78,12 @@ const SupportedKeywordItems<double>&
 fullySupported()
 {
    static const SupportedKeywordItems<double> fully_supported_keywords_double = {
+         {
+            "NEFAC",
+            {
+               {2,{true, [](double x) { return x > 0 && x <= 1.0; }, "NEFAC(EFF_FACTOR: Efficiency must be in the range (0,1]"}}, // NETWORK_EFF_FACTOR
+            },
+         },
          {
             "WPIMULT",
             {

opm/simulators/utils/PartiallySupportedFlowKeywords.cpp

@@ -111,12 +111,6 @@ partiallySupported()
                {8,{true, allow_values<std::string> {"NO"}, "GECON(ENDRUN): End run not implemented"}},
             },
          },
-         {
-            "GEFAC",
-            {
-               {3,{true, allow_values<std::string> {"YES"}, "GEFAC(GRPNETWK): Extended Network Model efficiency NO option not implemented"}}, // TRANSFER_EXT_NET
-            },
-         },
          {
             "GRIDOPTS",
             {
@@ -273,12 +267,6 @@ partiallySupported()
                {5,{true, allow_values<std::string> {"NO"}, "WAGHYSTR(WATER_MODEL): only the NO option is supported – will STOP"}}, // WATER_MODEL
             },
          },
-         {
-            "WEFAC",
-            {
-               {3,{true, allow_values<std::string> {"YES"}, "WEFAC(WELNETWK): only the YES option is supported"}}, // EXTENDED_NETWORK_OPT
-            },
-         },
          {
             "WELSPECS",
             {

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -397,7 +397,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"NARROW", {true, std::nullopt}},
         {"NCONSUMP", {true, std::nullopt}},
         {"NCOMPS", {false, std::nullopt}},
-        {"NEFAC", {true, std::nullopt}},
         {"NETCOMPA", {true, std::nullopt}},
         {"NEXT", {false, std::nullopt}},
         {"NEXTSTPL", {true, std::nullopt}},

opm/simulators/wells/BlackoilWellModel.hpp

@@ -118,13 +118,13 @@ template<class Scalar> class WellContributions;
 
             // TODO: where we should put these types, WellInterface or Well Model?
             // or there is some other strategy, like TypeTag
-            typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
+            using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
-            typedef Dune::BlockVector<VectorBlockType> BVector;
+            using BVector = Dune::BlockVector<VectorBlockType>;
 
-            typedef BlackOilPolymerModule<TypeTag> PolymerModule;
+            using PolymerModule = BlackOilPolymerModule<TypeTag>;
-            typedef BlackOilMICPModule<TypeTag> MICPModule;
+            using MICPModule = BlackOilMICPModule<TypeTag>;
 
-            // For the conversion between the surface volume rate and resrevoir voidage rate
+            // For the conversion between the surface volume rate and reservoir voidage rate
             using RateConverterType = RateConverter::
                 SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
 

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1373,6 +1373,10 @@ updateAndCommunicateGroupData(const int reportStepIdx,
                                                                 reportStepIdx,
                                                                 well_state_nupcol,
                                                                 this->groupState());
+    WellGroupHelpers<Scalar>::updateNetworkLeafNodeProductionRates(schedule(),
+                                                                   reportStepIdx,
+                                                                   well_state_nupcol,
+                                                                   this->groupState());
 
     WellGroupHelpers<Scalar>::updateGroupProductionRates(fieldGroup,
                                                          schedule(),

opm/simulators/wells/GroupState.cpp

@@ -41,6 +41,7 @@ GroupState<Scalar> GroupState<Scalar>::serializationTestObject()
 {
     GroupState result(3);
     result.m_production_rates = {{"test1", {1.0, 2.0}}};
+    result.m_network_leaf_node_production_rates={{"test1", {1.0, 20}}};
     result.production_controls = {{"test2", Group::ProductionCMode::LRAT}};
     result.prod_red_rates = {{"test3", {3.0, 4.0, 5.0}}};
     result.inj_red_rates = {{"test4", {6.0, 7.0}}};
@@ -61,6 +62,7 @@ template<class Scalar>
 bool GroupState<Scalar>::operator==(const GroupState& other) const
 {
     return this->m_production_rates == other.m_production_rates &&
+           this->m_network_leaf_node_production_rates == other.m_network_leaf_node_production_rates &&
            this->production_controls == other.production_controls &&
            this->prod_red_rates == other.prod_red_rates &&
            this->inj_red_rates == other.inj_red_rates &&
@@ -93,6 +95,16 @@ void GroupState<Scalar>::update_production_rates(const std::string& gname,
     this->m_production_rates[gname] = rates;
 }
 
+template<class Scalar>
+void GroupState<Scalar>::update_network_leaf_node_production_rates(const std::string& gname,
+                                                 const std::vector<Scalar>& rates)
+{
+    if (rates.size() != this->num_phases)
+        throw std::logic_error("Wrong number of phases");
+
+    this->m_network_leaf_node_production_rates[gname] = rates;
+}
+
 template<class Scalar>
 const std::vector<Scalar>&
 GroupState<Scalar>::production_rates(const std::string& gname) const
@@ -104,6 +116,17 @@ GroupState<Scalar>::production_rates(const std::string& gname) const
     return group_iter->second;
 }
 
+template<class Scalar>
+const std::vector<Scalar>&
+GroupState<Scalar>::network_leaf_node_production_rates(const std::string& gname) const
+{
+    auto group_iter = this->m_network_leaf_node_production_rates.find(gname);
+    if (group_iter == this->m_network_leaf_node_production_rates.end())
+        throw std::logic_error("No such group");
+
+    return group_iter->second;
+}
+
 //-------------------------------------------------------------------------
 
 template<class Scalar>

opm/simulators/wells/GroupState.hpp

@@ -52,7 +52,10 @@ public:
     bool has_production_rates(const std::string& gname) const;
     void update_production_rates(const std::string& gname,
                                  const std::vector<Scalar>& rates);
+    void update_network_leaf_node_production_rates(const std::string& gname,
+                                 const std::vector<Scalar>& rates);
     const std::vector<Scalar>& production_rates(const std::string& gname) const;
+    const std::vector<Scalar>& network_leaf_node_production_rates(const std::string& gname) const;
 
     void update_well_group_thp(const std::string& gname, const double& thp);
     Scalar well_group_thp(const std::string& gname) const;
@@ -130,6 +133,7 @@ public:
 
         auto forAllGroupData = [&](auto& func) {
             iterateContainer(m_production_rates, func);
+            iterateContainer(m_network_leaf_node_production_rates, func);
             iterateContainer(prod_red_rates, func);
             iterateContainer(inj_red_rates, func);
             iterateContainer(inj_resv_rates, func);
@@ -187,6 +191,7 @@ public:
     {
         serializer(num_phases);
         serializer(m_production_rates);
+        serializer(m_network_leaf_node_production_rates);
         serializer(production_controls);
         serializer(group_thp);
         serializer(prod_red_rates);
@@ -205,6 +210,7 @@ public:
 private:
     std::size_t num_phases{};
     std::map<std::string, std::vector<Scalar>> m_production_rates;
+    std::map<std::string, std::vector<Scalar>> m_network_leaf_node_production_rates;
     std::map<std::string, Group::ProductionCMode> production_controls;
     std::map<std::string, std::vector<Scalar>> prod_red_rates;
     std::map<std::string, std::vector<Scalar>> inj_red_rates;

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -152,6 +152,7 @@ update(const WellState<Scalar>& well_state,
             }
         }
     }
+    init();
 }
 
 template<class FluidSystem, class Indices>
@@ -208,6 +209,7 @@ updateNewton(const BVectorWell& dwells,
     if (stop_or_zero_rate_target) {
         value_[0][WQTotal] = 0.;
     }
+    init();
 }
 
 template<class FluidSystem, class Indices>
@@ -515,7 +517,7 @@ template<typename FluidSystem, typename Indices>
 typename MultisegmentWellPrimaryVariables<FluidSystem,Indices>::EvalWell
 MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
 volumeFraction(const int seg,
-               const unsigned compIdx) const
+               const int compIdx) const
 {
     if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
         return evaluation_[seg][WFrac];
@@ -584,7 +586,7 @@ typename MultisegmentWellPrimaryVariables<FluidSystem,Indices>::EvalWell
 MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
 getSegmentRateUpwinding(const int seg,
                         const int seg_upwind,
-                        const std::size_t comp_idx) const
+                        const int comp_idx) const
 {
     // the result will contain the derivative with respect to WQTotal in segment seg,
     // and the derivatives with respect to WFrac GFrac in segment seg_upwind.

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -115,7 +115,7 @@ public:
     //! \brief Returns upwinding rate for a component in a segment.
     EvalWell getSegmentRateUpwinding(const int seg,
                                      const int seg_upwind,
-                                     const std::size_t comp_idx) const;
+                                     const int comp_idx) const;
 
     //! \brief Get bottomhole pressure.
     EvalWell getBhp() const;
@@ -154,7 +154,7 @@ private:
 
     //! \brief Returns volume fraction for component in a segment.
     EvalWell volumeFraction(const int seg,
-                            const unsigned compIdx) const;
+                            const int compIdx) const;
 
     //! \brief The values for the primary variables
     //! \details Based on different solution strategies, the wells can have different primary variables

opm/simulators/wells/RatioCalculator.cpp

@@ -51,9 +51,9 @@ namespace Opm {
 
 template<class Value>
 RatioCalculator<Value>::
-RatioCalculator(unsigned gasCompIdx,
+RatioCalculator(int gasCompIdx,
-                unsigned oilCompIdx,
+                int oilCompIdx,
-                unsigned waterCompIdx,
+                int waterCompIdx,
                 std::string_view name)
     : gasComp_{gasCompIdx}
     , oilComp_(oilCompIdx)

opm/simulators/wells/RatioCalculator.hpp

@@ -39,9 +39,9 @@ class RatioCalculator
 public:
     using Scalar = decltype(getValue(Value{}));
 
-    RatioCalculator(unsigned gasCompIdx,
+    RatioCalculator(int gasCompIdx,
-                    unsigned oilCompIdx,
+                    int oilCompIdx,
-                    unsigned waterCompIdx,
+                    int waterCompIdx,
                     std::string_view name);
 
     void disOilVapWatVolumeRatio(Value& volumeRatio,
@@ -91,9 +91,9 @@ public:
                               const bool isProducer) const;
 
 private:
-    unsigned gasComp_;
+    int gasComp_;
-    unsigned oilComp_;
+    int oilComp_;
-    unsigned waterComp_;
+    int waterComp_;
     std::string name_;
 };
 

opm/simulators/wells/StandardWellConnections.cpp

@@ -731,7 +731,7 @@ connectionRateFoam(const std::vector<EvalWell>& cq_s,
             }
             case Phase::SOLVENT: {
                 if constexpr (Indices::enableSolvent)
-                    return static_cast<unsigned>(Indices::contiSolventEqIdx);
+                    return Indices::contiSolventEqIdx;
                 else
                     OPM_DEFLOG_THROW(std::runtime_error, "Foam transport phase is SOLVENT but SOLVENT is not activated.", deferred_logger);
             }

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -223,6 +223,7 @@ update(const WellState<Scalar>& well_state,
 
     // BHP
     value_[Bhp] = ws.bhp;
+    init();
 }
 
 template<class FluidSystem, class Indices>
@@ -301,6 +302,7 @@ updateNewton(const BVectorWell& dwells,
     // so that bhp constaint can be an active control when needed.
     constexpr Scalar bhp_lower_limit = 1. * unit::barsa - 1. * unit::Pascal;
     value_[Bhp] = std::max(value_[Bhp] - dx1_limited, bhp_lower_limit);
+    init();
 }
 
 template<class FluidSystem, class Indices>
@@ -320,6 +322,7 @@ updateNewtonPolyMW(const BVectorWell& dwells)
             value_[pskin_index] -= relaxation_factor * dx_pskin;
         }
     }
+    init();
 }
 
 template<class FluidSystem, class Indices>
@@ -446,7 +449,7 @@ copyToWellStatePolyMW(WellState<Scalar>& well_state) const
 template<class FluidSystem, class Indices>
 typename StandardWellPrimaryVariables<FluidSystem,Indices>::EvalWell
 StandardWellPrimaryVariables<FluidSystem,Indices>::
-volumeFraction(const unsigned compIdx) const
+volumeFraction(const int compIdx) const
 {
     if (FluidSystem::numActivePhases() == 1) {
         return EvalWell(numWellEq_ + Indices::numEq, 1.0);
@@ -456,7 +459,7 @@ volumeFraction(const unsigned compIdx) const
         return evaluation_[GFrac];
     }
 
-    if (Indices::enableSolvent && compIdx == (unsigned)Indices::contiSolventEqIdx) {
+    if (Indices::enableSolvent && compIdx == Indices::contiSolventEqIdx) {
         return evaluation_[SFrac];
     }
 

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -157,7 +157,7 @@ private:
                                              DeferredLogger& deferred_logger) const;
 
     //! \brief Returns volume fraction for a component.
-    EvalWell volumeFraction(const unsigned compIdx) const;
+    EvalWell volumeFraction(const int compIdx) const;
 
     //! \brief Handle non-reasonable fractions due to numerical overshoot.
     void processFractions();

opm/simulators/wells/WellGroupHelpers.cpp

@@ -90,14 +90,15 @@ namespace Opm {
                       const Opm::WellState<Scalar>& wellState,
                       const int reportStepIdx,
                       const int phasePos,
-                      const bool injector)
+                      const bool injector,
+                      const bool network)
     {
 
         Scalar rate = 0.0;
         for (const std::string& groupName : group.groups()) {
             const auto& groupTmp = schedule.getGroup(groupName, reportStepIdx);
-            const auto& gefac = groupTmp.getGroupEfficiencyFactor();
+            const auto& gefac = groupTmp.getGroupEfficiencyFactor(network);
-            rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector);
+            rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector, network);
         }
 
         // only sum satelite production once
@@ -126,8 +127,7 @@ namespace Opm {
             if (wellEcl.getStatus() == Opm::Well::Status::SHUT)
                 continue;
 
-            const Scalar factor = wellEcl.getEfficiencyFactor() *
+            const Scalar factor = wellEcl.getEfficiencyFactor(network) * wellState[wellEcl.name()].efficiency_scaling_factor;
-                                  wellState[wellEcl.name()].efficiency_scaling_factor;
             const auto& ws = wellState.well(well_index.value());
             if (res_rates) {
                 const auto& well_rates = ws.reservoir_rates;
@@ -744,6 +744,30 @@ updateGroupProductionRates(const Group& group,
     group_state.update_production_rates(group.name(), rates);
 }
 
+template<class Scalar>
+void WellGroupHelpers<Scalar>::
+updateNetworkLeafNodeProductionRates(const Schedule& schedule,
+                                     const int reportStepIdx,
+                                     const WellState<Scalar>& wellState,
+                                     GroupState<Scalar>& group_state)
+{
+    const auto& network = schedule[reportStepIdx].network();
+    if (network.active()) {
+        const int np = wellState.numPhases();
+        for (const auto& group_name : network.leaf_nodes()) {
+            assert(schedule[reportStepIdx].groups.has(group_name));
+            const auto& group = schedule[reportStepIdx].groups.get(group_name);
+            std::vector<Scalar> network_rates(np, 0.0);
+            if (group.numWells() > 0) {
+                for (int phase = 0; phase < np; ++phase) {
+                    network_rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false, /*network*/ true);
+                }
+            }
+            group_state.update_network_leaf_node_production_rates(group_name, network_rates);
+        }
+    }
+}
+
 template<class Scalar>
 void WellGroupHelpers<Scalar>::
 updateREINForGroups(const Group& group,
@@ -928,23 +952,15 @@ computeNetworkPressures(const Network::ExtNetwork& network,
                 continue;
             }
 
-            node_inflows[node] = group_state.production_rates(node);
+            node_inflows[node] = group_state.network_leaf_node_production_rates(node);
             // Add the ALQ amounts to the gas rates if requested.
             if (network.node(node).add_gas_lift_gas()) {
                 const auto& group = schedule.getGroup(node, report_time_step);
                 for (const std::string& wellname : group.wells()) {
                     const Well& well = schedule.getWell(wellname, report_time_step);
+
                     if (well.isInjector() || !well_state.isOpen(wellname)) continue;
-                    // Here we use the efficiency unconditionally, but if WEFAC item 3
+                    const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) * well_state.getGlobalEfficiencyScalingFactor(wellname);
-                    // for the well is false (it defaults to true) then we should NOT use
-                    // the efficiency factor. Fixing this requires not only changing the
-                    // code here, but also:
-                    //    - Adding a member to the well for this flag, and setting it in Schedule::handleWEFAC().
-                    //    - Making the wells' maximum flows (i.e. not time-averaged by using a efficiency factor)
-                    //      available and using those (for wells with WEFAC(3) true only) when accumulating group
-                    //      rates, but ONLY for network calculations.
-                    const Scalar efficiency = well.getEfficiencyFactor() *
-                                              well_state.getGlobalEfficiencyScalingFactor(wellname);
                     node_inflows[node][BlackoilPhases::Vapour] += well_state.getALQ(wellname) * efficiency;
                 }
             }
@@ -962,13 +978,14 @@ computeNetworkPressures(const Network::ExtNetwork& network,
                 // Add downbranch rates to upbranch.
                 std::vector<Scalar>& up = node_inflows[(*upbranch).uptree_node()];
                 const std::vector<Scalar>& down = node_inflows[node];
+                // We now also support NEFAC
+                const Scalar efficiency = network.node(node).efficiency();
                 if (up.empty()) {
-                    up = down;
+                    up = std::vector<Scalar>(down.size(), 0.0);
-                } else {
+                }
                 assert (up.size() == down.size());
                 for (std::size_t ii = 0; ii < up.size(); ++ii) {
-                        up[ii] += down[ii];
+                    up[ii] += efficiency*down[ii];
-                    }
                 }
             }
         }

opm/simulators/wells/WellGroupHelpers.hpp

@@ -56,7 +56,8 @@ public:
                                     const Opm::WellState<Scalar>& wellState,
                                     const int reportStepIdx,
                                     const int phasePos,
-                                    const bool injector);
+                                    const bool injector,
+                                    const bool network = false);
 
     static Scalar satelliteProduction(const ScheduleState& sched,
                                       const std::vector<std::string>& groups,
@@ -199,6 +200,12 @@ public:
                                            const WellState<Scalar>& wellState,
                                            GroupState<Scalar>& group_state);
 
+    static void updateNetworkLeafNodeProductionRates(const Schedule& schedule,
+                                                     const int reportStepIdx,
+                                                     const WellState<Scalar>& wellState,
+                                                     GroupState<Scalar>& group_state);
+
+
     static void updateWellRatesFromGroupTargetScale(const Scalar scale,
                                                     const Group& group,
                                                     const Schedule& schedule,

opm/simulators/wells/WellInterfaceIndices.cpp

@@ -78,7 +78,7 @@ flowPhaseToModelCompIdx(const int phaseIdx) const
 template<class FluidSystem, class Indices>
 int
 WellInterfaceIndices<FluidSystem,Indices>::
-modelCompIdxToFlowCompIdx(const unsigned compIdx) const
+modelCompIdxToFlowCompIdx(const int compIdx) const
 {
     const auto& pu = this->phaseUsage();
     if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == compIdx)

opm/simulators/wells/WellInterfaceIndices.hpp

@@ -41,7 +41,7 @@ public:
     using ModelParameters = typename WellInterfaceFluidSystem<FluidSystem>::ModelParameters;
 
     int flowPhaseToModelCompIdx(const int phaseIdx) const;
-    int modelCompIdxToFlowCompIdx(const unsigned compIdx) const;
+    int modelCompIdxToFlowCompIdx(const int compIdx) const;
     Scalar scalingFactor(const int phaseIdx) const;
 
     template <class EvalWell>