Thermal msw refactoring

opm/simulators/timestepping/ConvergenceReport.cpp

@@ -76,6 +76,8 @@ namespace Opm
             return "MassBalance";
         case T::Pressure:
             return "Pressure";
+        case T::Energy:
+            return "Energy";
         case T::ControlBHP:
             return "ControlBHP";
         case T::ControlTHP:

opm/simulators/timestepping/ConvergenceReport.hpp

@@ -164,6 +164,7 @@ namespace Opm
                 Invalid,
                 MassBalance,
                 Pressure,
+                Energy,
                 ControlBHP,
                 ControlTHP,
                 ControlRate,

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -2092,6 +2092,12 @@ namespace Opm {
     int
     BlackoilWellModel<TypeTag>::numConservationQuantities() const
     {
+        // TODO: when the energy equation joins, here should we start
+        // of the refactoring related to the the usage of the numConservationQuantities()
+        // since energy equation is also a conservation equation
+        // at the current stage, we only have energy equations for MSW, so we should not
+        // refactor at this stage.
+
         // The numPhases() functions returns 1-3, depending on which
         // of the (oil, water, gas) phases are active. For each of those phases,
         // if the phase is active the corresponding component is present and

opm/simulators/wells/MSWellHelpers.cpp

@@ -376,7 +376,9 @@ using Mat = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,M,N>>;
     INSTANTIATE_PARALLELLMSWELLB(T, 3, 4)   \
     INSTANTIATE_PARALLELLMSWELLB(T, 4, 3)   \
     INSTANTIATE_PARALLELLMSWELLB(T, 4, 4)   \
-    INSTANTIATE_PARALLELLMSWELLB(T, 4, 5)
+    INSTANTIATE_PARALLELLMSWELLB(T, 4, 5)   \
+    INSTANTIATE_PARALLELLMSWELLB(T, 5, 4)   \
+    INSTANTIATE_PARALLELLMSWELLB(T, 5, 5)
 
 #define INSTANTIATE_UMF(T,Dim)                                             \
     template Vec<T,Dim> applyUMFPack(Dune::UMFPack<Mat<T,Dim>>&,           \
@@ -414,13 +416,15 @@ using Mat = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,M,N>>;
     INSTANTIATE_UMF(T,2)    \
     INSTANTIATE_UMF(T,3)    \
     INSTANTIATE_UMF(T,4)    \
+    INSTANTIATE_UMF(T,5)    \
     INSTANTIATE_EVAL(T,3)   \
     INSTANTIATE_EVAL(T,4)   \
     INSTANTIATE_EVAL(T,5)   \
     INSTANTIATE_EVAL(T,6)   \
     INSTANTIATE_EVAL(T,7)   \
     INSTANTIATE_EVAL(T,8)   \
     INSTANTIATE_EVAL(T,9)   \
+    INSTANTIATE_EVAL(T,10)  \
     INSTANTIATE_ALL_PARALLELLMSWELLB(T)
 
 INSTANTIATE_TYPE(double)

opm/simulators/wells/MultisegmentWell.hpp

@@ -27,6 +27,9 @@
 #include <opm/simulators/wells/WellInterface.hpp>
 #include <opm/simulators/wells/MultisegmentWellEval.hpp>
 
+#include <limits>
+#include <string_view>
+
 namespace Opm {
 
     class DeferredLogger;
@@ -56,12 +59,36 @@ namespace Opm {
 
         using Base::has_solvent;
         using Base::has_polymer;
+        using Base::has_energy;
+        using Base::has_brine;
         using Base::Water;
         using Base::Oil;
         using Base::Gas;
 
         using typename Base::Scalar;
 
+        // True when the composition switch primary variable is active, i.e. both oil and
+        // gas phases are present so that Rs/Rv are stored in the fluid state. Matches the
+        // fluid state's enableDissolution flag (see WellInterface::BlackOilFluidStateType).
+        static constexpr bool compositionSwitchEnabled =
+            Indices::compositionSwitchIdx != std::numeric_limits<unsigned>::max();
+
+        // True when the segment fluid state stores a temperature (any thermal mode, not just
+        // the fully implicit one). Matches the fluid state's enableTemperature flag (see
+        // WellInterface::BlackOilFluidStateType); used to decide whether createFluidState()
+        // must set the temperature explicitly.
+        static constexpr bool enable_temperature =
+            Base::energyModuleType != EnergyModules::NoTemperature;
+
+        // Scaling factor applied to the well-side energy equation so that its residual
+        // lands on the same numerical scale as the mass-balance equations (the raw
+        // energy flux/accumulation is ~1e4-1e9, while mass residuals are ~O(1)).
+        // We reuse the very same factor the reservoir energy equation uses (see
+        // BlackOilEnergyScalingFactor in blackoilmodel.hh) so the coupled well and
+        // reservoir energy equations live on the same scale.
+        static constexpr Scalar energy_scaling_factor_ =
+            getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+
         /// the matrix and vector types for the reservoir
         using typename Base::BVector;
         using typename Base::Eval;
@@ -73,6 +100,12 @@ namespace Opm {
         using typename Base::PressureMatrix;
         using FSInfo = std::tuple<Scalar, Scalar>;
 
+        // a fluid state to calculate the properties inside the wellbore for each segment
+        // it will be probably used for more things, but at the moment, it is for the enthalpy
+        // calculation in the wellbore.
+        template <typename ValueType>
+        using SegmentFluidState = Base::template BlackOilFluidStateType<ValueType>;
+
         MultisegmentWell(const Well& well,
                          const ParallelWellInfo<Scalar>& pw_info,
                          const int time_step,
@@ -173,6 +206,16 @@ namespace Opm {
 
         // the intial amount of fluids in each segment under surface condition
         std::vector<std::vector<Scalar> > segment_fluid_initial_;
+        // total energy inside the segments at the beginning of the time step
+        std::vector<Scalar> segment_initial_energy_;
+
+        // segment fluid state
+        std::vector<SegmentFluidState<EvalWell>> segment_fluid_state_;
+
+        // fluid state under the wellhead condition, it is used to calculate the enthalpy
+        // under operation condition for energy injection
+        // because BHP will be involved, we use EvalWell type here
+        SegmentFluidState<EvalWell> wellhead_fluid_state_;
 
         mutable int debug_cost_counter_ = 0;
 
@@ -184,7 +227,7 @@ namespace Opm {
                              const Scalar relaxation_factor = 1.0);
 
         // computing the accumulation term for later use in well mass equations
-        void computeInitialSegmentFluids(const Simulator& simulator, DeferredLogger& deferred_logger);
+        void computeInitialSegmentFluids(const FSInfo& info, DeferredLogger& deferred_logger);
 
         // compute the pressure difference between the perforation and cell center
         void computePerfCellPressDiffs(const Simulator& simulator);
@@ -289,8 +332,8 @@ namespace Opm {
 
         void updateWaterThroughput(const double dt, WellStateType& well_state) const override;
 
-        EvalWell getSegmentSurfaceVolume(const Simulator& simulator,
-                                         const int seg_idx,
+        EvalWell getSegmentSurfaceVolume(const int seg_idx,
+                                         const FSInfo& info,
                                          DeferredLogger& deferred_logger) const;
 
         // turn on crossflow to avoid singular well equations
@@ -332,9 +375,70 @@ namespace Opm {
                        DeferredLogger& deferred_logger) const override;
 
         FSInfo getFirstPerforationFluidStateInfo(const Simulator& simulator) const;
+
+        // this function can potentially be shared between multisegment wells and standard wells
+        // The optional @p volume_ratio output receives the segment volume ratio (reservoir
+        // volume per unit surface volume, i.e. the sum of the unnormalized phase saturations),
+        // so callers can cache and reuse it for the segment surface volume.
+        template <typename ValueType = EvalWell>
+        SegmentFluidState<ValueType>
+        createFluidState(const std::vector<ValueType>& fluid_composition,
+                         const ValueType& pressure,
+                         const ValueType& temperature,
+                         DeferredLogger& deferred_logger,
+                         ValueType* volume_ratio = nullptr) const;
+
+        SegmentFluidState<EvalWell>
+        createSegmentFluidState(int seg, const FSInfo& info, DeferredLogger& deferred_logger,
+                                EvalWell* volume_ratio = nullptr) const;
+
+        void computeInitialSegmentEnergy();
+
+        // assemble the energy equation contribution for a single perforation/connection
+        void assemblePerforationEnergyEq(const IntensiveQuantities& int_quants,
+                                         const std::vector<EvalWell>& cq_s,
+                                         const int seg,
+                                         const int local_perf_index,
+                                         DeferredLogger& deferred_logger);
+
+        void updateWellHeadCondition(const Simulator& simulator,
+                                     const Scalar first_perf_temperature,
+                                     DeferredLogger& deferred_logger);
+
+        void updateSegmentFluidState(const FSInfo& info, DeferredLogger& deferred_logger);
+
+        template <typename ValueType = EvalWell>
+        ValueType computeSegmentEnergy(int seg) const;
+
+        // Convert per-component surface volumetric rates to a phase reservoir
+        // volumetric rate using @p fs as the upwind fluid state. Handles the
+        // dissolved-gas/vaporized-oil coupling (rs/rv). When the determinant
+        // (1 - rs*rv) is non-positive, falls back to the no-dissolution case
+        // and logs a debug message tagged with @p context.
+        // @p fs may be either a wellbore SegmentFluidState (properties already
+        // EvalWell) or a reservoir-cell intensive-quantity fluid state (Eval,
+        // extended to EvalWell on the fly).
+        // @return the reservoir volumetric rate of @p phaseIdx.
+        template <typename FluidStateT>
+        EvalWell surfaceToReservoirRate(unsigned phaseIdx,
+                                        const FluidStateT& fs,
+                                        const std::vector<EvalWell>& surface_rates,
+                                        int seg,
+                                        std::string_view context,
+                                        DeferredLogger& deferred_logger) const;
+
+        // Compute the energy flux carried by fluid flowing from segment
+        // @p seg toward its outlet, using @p upwind_fs as the upwind
+        // fluid state. @p context is used in diagnostic messages.
+        // @return the energy flux as sum_phases(reservoir_rate * enthalpy * density).
+        EvalWell computeSegmentEnergyRate(int seg,
+                                          int upwind_seg,
+                                          const SegmentFluidState<EvalWell>& upwind_fs,
+                                          std::string_view context,
+                                          DeferredLogger& deferred_logger) const;
     };
 
-}
+} // namespace Opm
 
 #include "MultisegmentWell_impl.hpp"
 

opm/simulators/wells/MultisegmentWellAssemble.cpp

@@ -352,7 +352,17 @@ assembleOutflowTerm(const int seg,
     if constexpr (has_gfrac_variable) {
         eqns.D()[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
     }
-    // pressure derivative should be zero
+    // pressure derivative should be zero for non-thermal cases
+
+    if constexpr (enable_energy) {
+        // Only the energy equation depends on pressure and temperature; the mass
+        // surface rates carry no SPres/Temperature derivatives.
+        if (comp_idx == Temperature) {
+            // energy flux depends on the pressure and temperature
+            eqns.D()[seg][seg_upwind][comp_idx][SPres] -= segment_rate.derivative(SPres + Indices::numEq);
+            eqns.D()[seg][seg_upwind][comp_idx][Temperature] -= segment_rate.derivative(Temperature + Indices::numEq);
+        }
+    }
 }
 
 template<class FluidSystem, class Indices>
@@ -377,7 +387,17 @@ assembleInflowTerm(const int seg,
     if constexpr (has_gfrac_variable) {
         eqns.D()[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
     }
-    // pressure derivative should be zero
+    // pressure derivative should be zero for non-thermal cases
+
+    if constexpr (enable_energy) {
+        // Only the energy equation depends on pressure and temperature; the mass
+        // surface rates carry no SPres/Temperature derivatives.
+        if (comp_idx == Temperature) {
+            // energy flux depends on the pressure and temperature
+            eqns.D()[seg][inlet_upwind][comp_idx][SPres] += inlet_rate.derivative(SPres + Indices::numEq);
+            eqns.D()[seg][inlet_upwind][comp_idx][Temperature] += inlet_rate.derivative(Temperature + Indices::numEq);
+        }
+    }
 }
 
 template<class FluidSystem, class Indices>

opm/simulators/wells/MultisegmentWellAssemble.hpp

@@ -47,12 +47,14 @@ class MultisegmentWellAssemble
     static constexpr int WQTotal = PrimaryVariables::WQTotal;
     static constexpr bool has_wfrac_variable = PrimaryVariables::has_wfrac_variable;
     static constexpr bool has_gfrac_variable = PrimaryVariables::has_gfrac_variable;
+    static constexpr bool enable_energy = PrimaryVariables::enable_energy;
     static constexpr int WFrac = PrimaryVariables::WFrac;
     static constexpr int GFrac = PrimaryVariables::GFrac;
+    static constexpr int Temperature = PrimaryVariables::Temperature;
     static constexpr int SPres = PrimaryVariables::SPres;
 
 public:
-    static constexpr int numWellEq = Indices::numPhases+1;
+    static constexpr int numWellEq = Indices::numPhases + 1 + enable_energy;
     using Scalar = typename FluidSystem::Scalar;
     using IndexTraits = typename FluidSystem::IndexTraitsType;
     using Equations = MultisegmentWellEquations<Scalar, IndexTraits, numWellEq,Indices::numEq>;

opm/simulators/wells/MultisegmentWellEquations.cpp

@@ -484,7 +484,9 @@ sumDistributed(Parallel::Communication comm)
     INSTANTIATE(T,3,4)      \
     INSTANTIATE(T,4,3)      \
     INSTANTIATE(T,4,4)      \
-    INSTANTIATE(T,4,5)
+    INSTANTIATE(T,4,5)      \
+    INSTANTIATE(T,5,4)      \
+    INSTANTIATE(T,5,5)
 
 INSTANTIATE_TYPE(double)
 

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -89,7 +89,11 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
                    const bool relax_tolerance,
                    const bool well_is_stopped) const
 {
-    assert(int(B_avg.size()) == baseif_.numConservationQuantities());
+    // TODO: numConservationQuantities() should be refactored to incorporate the energy equation
+    // TODO: when enable_energy is true, at the current stage, B_avg can be baseif_.numConservationQuantities() if calculated in BlackoilWellModel
+    // or baseif_.numConservationQuantities() + 1 if calculated in getReservoirConvergence() in BlackoilModel
+    // we are only accessing the first baseif_.numConservationQuantities() elements with current form of the code
+    assert(int(B_avg.size()) == baseif_.numConservationQuantities() || enable_energy);
 
     // checking if any residual is NaN or too large. The two large one is only handled for the well flux
     std::vector<std::vector<Scalar>> abs_residual(this->numberOfSegments(),
@@ -103,14 +107,19 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
     std::vector<Scalar> maximum_residual(numWellEq, 0.0);
 
     ConvergenceReport report;
-    // TODO: the following is a little complicated, maybe can be simplified in some way?
+    // TODO: the energy equation should not have a B_avg here to use, maybe we can use 1 for it
     for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
         for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
             if (eq_idx < baseif_.numConservationQuantities()) { // phase or component mass equations
                 const Scalar flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
                 if (flux_residual > maximum_residual[eq_idx]) {
                     maximum_residual[eq_idx] = flux_residual;
                 }
+            } else if (enable_energy && eq_idx == Temperature) { // energy equation
+                const Scalar energy_residual = abs_residual[seg][eq_idx];
+                if (energy_residual > maximum_residual[eq_idx]) {
+                    maximum_residual[eq_idx] = energy_residual;
+                }
             } else { // pressure or control equation
                 // for the top segment (seg == 0), it is control equation, will be checked later separately
                 if (seg > 0) {
@@ -139,7 +148,24 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
             } else if (flux_residual > relaxed_inner_tolerance_flow_ms_well) {
                 report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, baseif_.name()});
             }
-        } else { // pressure equation
+        } else if (enable_energy && eq_idx == Temperature) {
+            const Scalar energy_residual = maximum_residual[eq_idx];
+            // TODO: possibly the dummy_phase should be something else, while requires extension to WellConvergenceMetric
+            constexpr int dummy_phase = -1;
+            // The well-side energy equation is scaled (by MultisegmentWell::energy_scaling_factor_,
+            // the same factor the reservoir energy equation uses) so that its residual is on the
+            // same magnitude as the mass-balance equations. We can therefore reuse the mass-balance
+            // well tolerances rather than a hand-tuned energy-specific constant.
+            if (std::isnan(energy_residual)) {
+                report.setWellFailed({CR::WellFailure::Type::Energy, CR::Severity::NotANumber, dummy_phase, baseif_.name()});
+            } else if (energy_residual > max_residual_allowed) {
+                report.setWellFailed({CR::WellFailure::Type::Energy, CR::Severity::TooLarge, dummy_phase, baseif_.name()});
+            } else if (!relax_tolerance && energy_residual > tolerance_wells) {
+                report.setWellFailed({CR::WellFailure::Type::Energy, CR::Severity::Normal, dummy_phase, baseif_.name()});
+            } else if (energy_residual > relaxed_inner_tolerance_flow_ms_well) {
+                report.setWellFailed({CR::WellFailure::Type::Energy, CR::Severity::Normal, dummy_phase, baseif_.name()});
+            }
+        } else if (eq_idx == SPres) { // pressure equation
             const Scalar pressure_residual = maximum_residual[eq_idx];
             const int dummy_component = -1;
             if (std::isnan(pressure_residual)) {
@@ -437,7 +463,10 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
             Scalar residual = 0.;
             if (eq_idx < baseif_.numConservationQuantities()) {
                 residual = std::abs(linSys_.residual()[seg][eq_idx]) * B_avg[eq_idx];
-            } else {
+            } else if (eq_idx == Temperature) {
+                residual = std::abs(linSys_.residual()[seg][eq_idx]);
+            } else if (eq_idx == SPres) {
+                // skip seg 0 because it holds the control equation
                 if (seg > 0) {
                     residual = std::abs(linSys_.residual()[seg][eq_idx]);
                 }
@@ -548,7 +577,10 @@ getResidualMeasureValue(const WellState<Scalar, IndexTraits>& well_state,
     const Scalar rate_tolerance = tolerance_wells;
     [[maybe_unused]] int count = 0;
     Scalar sum = 0;
-    for (int eq_idx = 0; eq_idx < numWellEq - 1; ++eq_idx) {
+    // Loop over mass balance equations
+    // TODO: we should probably also do something related to the energy equation also
+    // in the residual measures
+    for (int eq_idx = 0; eq_idx < baseif_.numConservationQuantities(); ++eq_idx) {
         if (residuals[eq_idx] > rate_tolerance) {
             sum += residuals[eq_idx] / rate_tolerance;
             ++count;

opm/simulators/wells/MultisegmentWellEval.hpp

@@ -54,6 +54,9 @@ class MultisegmentWellEval : public MultisegmentWellGeneric<typename FluidSystem
     static constexpr int SPres = PrimaryVariables::SPres;
     static constexpr int WQTotal = PrimaryVariables::WQTotal;
 
+    static constexpr bool enable_energy = PrimaryVariables::enable_energy;
+    static constexpr int Temperature = PrimaryVariables::Temperature;
+
     using Equations = MultisegmentWellEquations<Scalar, IndexTraits, numWellEq, Indices::numEq>;
     using MSWSegments = MultisegmentWellSegments<FluidSystem,Indices>;