Compositional well

opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.cpp

@@ -394,7 +394,7 @@ struct SummaryConfigContext {
     bool is_well_comp(const std::string& keyword)
     {
         static const auto well_comp_kw = keyword_set {
-             "WXMF", "WYMF", "WZMF", "WCGMR", "WCOMR"
+             "WAMF", "WXMF", "WYMF", "WZMF", "WCGMR", "WCOMR",
         };
 
         return is_in_set(well_comp_kw, keyword);

opm/input/eclipse/Schedule/ScheduleState.cpp

@@ -311,6 +311,7 @@ bool ScheduleState::operator==(const ScheduleState& other) const {
         && this->bhp_defaults.get() == other.bhp_defaults.get()
         && this->source.get() == other.source.get()
         && this->wells == other.wells
+        && this->inj_streams == other.inj_streams
         && this->groups == other.groups
         && this->vfpprod == other.vfpprod
         && this->vfpinj == other.vfpinj

opm/input/eclipse/Schedule/ScheduleState.hpp

@@ -201,6 +201,9 @@ namespace Opm {
                 return (ptr != nullptr);
             }
 
+            void update(const K& key, std::shared_ptr<T> value) {
+                this->m_data.insert_or_assign(key, std::move(value));
+            }
 
             void update(T object) {
                 auto key = object.name();
@@ -509,6 +512,9 @@ namespace Opm {
         // constant flux aquifers
         std::unordered_map<int, SingleAquiferFlux> aqufluxs;
         BCProp bcprop;
+        // injection streams for compostional STREAM injection using WINJGAS 
+        map_member<std::string, std::vector<double>> inj_streams;
+
         std::unordered_map<std::string, double> target_wellpi;
         std::optional<NextStep> next_tstep;
 
@@ -546,6 +552,7 @@ namespace Opm {
             serializer(wells);
             serializer(aqufluxs);
             serializer(bcprop);
+            serializer(inj_streams);
             serializer(target_wellpi);
             serializer(this->next_tstep);
             serializer(m_start_time);

opm/input/eclipse/Schedule/Well/Well.hpp

@@ -164,6 +164,9 @@ class Well {
 
         double rsRvInj;
 
+        // injection stream compostion for compositional simulation
+        std::optional<std::vector<double>> gas_inj_composition{};
+
         bool operator==(const WellInjectionProperties& other) const;
         bool operator!=(const WellInjectionProperties& other) const;
 
@@ -230,6 +233,9 @@ class Well {
         void update_uda(const UDQConfig& udq_config, UDQActive& udq_active, UDAControl control, const UDAValue& value);
         void handleWTMULT(Well::WELTARGCMode cmode, double factor);
 
+        void setGasInjComposition(const std::vector<double>& composition);
+        const std::vector<double>& gasInjComposition() const;
+
         template<class Serializer>
         void serializeOp(Serializer& serializer)
         {
@@ -248,6 +254,7 @@ class Well {
             serializer(injectorType);
             serializer(controlMode);
             serializer(rsRvInj);
+            serializer(gas_inj_composition);
         }
     };
 

opm/input/eclipse/Schedule/Well/WellInjectionProperties.cpp

@@ -60,7 +60,8 @@ namespace Opm {
           injectionControls(0),
           injectorType(InjectorType::WATER),
           controlMode(InjectorCMode::CMODE_UNDEFINED),
-          rsRvInj(0.0)
+          rsRvInj(0.0),
+          gas_inj_composition(std::nullopt)
     {
     }
 
@@ -82,6 +83,7 @@ namespace Opm {
         result.injectorType = InjectorType::OIL;
         result.controlMode = InjectorCMode::BHP;
         result.rsRvInj = 11;
+        result.gas_inj_composition = std::vector<double>{1.0, 2.0, 3.0};
 
         return result;
     }
@@ -284,7 +286,8 @@ namespace Opm {
             (injectionControls == other.injectionControls) &&
             (injectorType == other.injectorType) &&
             (controlMode == other.controlMode) &&
-            (rsRvInj == other.rsRvInj))
+            (rsRvInj == other.rsRvInj) &&
+            (gas_inj_composition == other.gas_inj_composition))
             return true;
         else
             return false;
@@ -500,5 +503,15 @@ namespace Opm {
         }
     }
 
+    void Well::WellInjectionProperties::setGasInjComposition(const std::vector<double>& composition) {
+        gas_inj_composition = composition;
+    }
+
+    const std::vector<double>& Well::WellInjectionProperties::gasInjComposition() const {
+        if (!gas_inj_composition.has_value()) {
+            throw std::invalid_argument("Gas injection composition not set");
+        }
+        return gas_inj_composition.value();
+    }
 
 }

opm/input/eclipse/Schedule/Well/WellKeywordHandlers.cpp

@@ -445,6 +445,29 @@ void handleWCYCLE(HandlerContext& handlerContext)
     handlerContext.state().wcycle.update(std::move(new_config));
 }
 
+void handleWELLSTRE(HandlerContext& handlerContext)
+{
+    auto& inj_streams = handlerContext.state().inj_streams;
+    for (const auto& record : handlerContext.keyword) {
+        const auto stream_name = record.getItem<ParserKeywords::WELLSTRE::STREAM>().getTrimmedString(0);
+        const auto& composition = record.getItem<ParserKeywords::WELLSTRE::COMPOSITIONS>().getSIDoubleData();
+        const std::size_t num_comps = handlerContext.static_schedule().m_runspec.numComps();
+        if (composition.size() != num_comps) {
+            const std::string msg = fmt::format("The number of the composition values for stream '{}' is not the same as the number of components.", stream_name);
+            throw OpmInputError(msg, handlerContext.keyword.location());
+        }
+
+        const double sum = std::accumulate(composition.begin(), composition.end(), 0.0);
+        if (std::abs(sum - 1.0) > std::numeric_limits<double>::epsilon()) {
+            const std::string msg = fmt::format("The sum of the composition values for stream '{}' is not 1.0, but {}.", stream_name, sum);
+            throw OpmInputError(msg, handlerContext.keyword.location());
+        }
+        auto composition_ptr = std::make_shared<std::vector<double>>(composition);
+        inj_streams.update(stream_name, std::move(composition_ptr));
+    }
+
+}
+
 void handleWELOPEN(HandlerContext& handlerContext)
 {
     const auto& keyword = handlerContext.keyword;
@@ -527,6 +550,44 @@ void handleWELOPEN(HandlerContext& handlerContext)
     }
 }
 
+void handleWINJGAS(HandlerContext& handlerContext)
+{
+    // \Note: we do not support the item 4 MAKEUPGAS and item 5 STAGE in WINJGAS keyword yet
+    for (const auto& record : handlerContext.keyword) {
+        const std::string fluid_nature = record.getItem<ParserKeywords::WINJGAS::FLUID>().getTrimmedString(0);
+
+        // \Note: technically, only the first two characters are significant
+        // with some testing, we can determine whether we want to enforce this.
+        // at the moment, we only support full string STREAM for fluid nature
+        if (fluid_nature != "STREAM") {
+            const std::string msg = fmt::format("The fluid nature '{}' is not supported in WINJGAS keyword.", fluid_nature);
+            throw OpmInputError(msg, handlerContext.keyword.location());
+        }
+
+        const std::string stream_name = record.getItem<ParserKeywords::WINJGAS::STREAM>().getTrimmedString(0);
+        // we make sure the stream is defined in WELLSTRE keyword
+        const auto& inj_streams = handlerContext.state().inj_streams;
+        if (!inj_streams.has(stream_name)) {
+            const std::string msg = fmt::format("The stream '{}' is not defined in WELLSTRE keyword.", stream_name);
+            throw OpmInputError(msg, handlerContext.keyword.location());
+        }
+
+        const std::string wellNamePattern = record.getItem<ParserKeywords::WINJGAS::WELL>().getTrimmedString(0);
+        const auto well_names = handlerContext.wellNames(wellNamePattern, false);
+        for (const auto& well_name : well_names) {
+            auto well2 = handlerContext.state().wells.get(well_name);
+            auto injection = std::make_shared<Well::WellInjectionProperties>(well2.getInjectionProperties());
+
+            const auto& inj_stream = inj_streams.get(stream_name);
+            injection->setGasInjComposition(inj_stream);
+
+            if (well2.updateInjection(injection)) {
+                handlerContext.state().wells.update(std::move(well2));
+            }
+        }
+    }
+}
+
 void handleWELSPECS(HandlerContext& handlerContext)
 {
     using Kw = ParserKeywords::WELSPECS;
@@ -945,11 +1006,13 @@ getWellHandlers()
         { "WCONPROD", &handleWCONPROD },
         { "WCYCLE",   &handleWCYCLE   },
         { "WELOPEN" , &handleWELOPEN  },
+        { "WELLSTRE", &handleWELLSTRE },
         { "WELSPECS", &handleWELSPECS },
         { "WELSPECL", &handleWELSPECL },
         { "WELTARG" , &handleWELTARG  },
         { "WELTRAJ" , &handleWELTRAJ  },
         { "WHISTCTL", &handleWHISTCTL },
+        { "WINJGAS",  &handleWINJGAS  },
         { "WLIST"   , &handleWLIST    },
         { "WPAVE"   , &handleWPAVE    },
         { "WPAVEDEP", &handleWPAVEDEP },

opm/input/eclipse/share/keywords/001_Eclipse300/W/WELLSTRE

@@ -10,16 +10,11 @@
       "value_type": "STRING"
     },
     {
-        "item": 2,
-        "name": "XMOLE1",
-        "value_type": "DOUBLE",
-        "dimension": "1"
-      },
-      {
-        "item": 3,
-        "name": "XMOLE2",
-        "value_type": "DOUBLE",
-        "dimension": "1"
-      }
+      "item": 2,
+      "name": "COMPOSITIONS",
+      "value_type": "DOUBLE",
+      "dimension": "1",
+      "size_type": "ALL"
+    }
   ]
 }

opm/input/eclipse/share/keywords/001_Eclipse300/W/WELL_PROBE_COMP

@@ -4,6 +4,7 @@
     "SUMMARY"
   ],
   "deck_names": [
+    "WAMF",
     "WXMF",
     "WYMF",
     "WZMF",

opm/material/constraintsolvers/PTFlash.hpp

@@ -80,7 +80,7 @@ class PTFlash
      *
      */
     template <class FluidState>
-    static void solve(FluidState& fluid_state,
+    static bool solve(FluidState& fluid_state,
                       const std::string& twoPhaseMethod,
                       Scalar flash_tolerance,
                       const EOSType& eos_type,
@@ -116,6 +116,8 @@ class PTFlash
 
         // we update the derivatives in fluid_state
         updateDerivatives_(fluid_state_scalar, fluid_state, eos_type, is_single_phase);
+
+        return is_single_phase;
     } //end solve
 
     /*!
@@ -992,7 +994,7 @@ class PTFlash
         // p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary
         // pressure joins
 
-        constexpr size_t num_deri = numComponents;
+        constexpr size_t num_deri = InputEval::numVars;
         for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
             std::vector<double> deri(num_deri, 0.);
             // derivatives from P