From 294466e4393a6e350383e6f8c30bef0b0ff68fbe Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 24 Feb 2026 12:43:06 +0000 Subject: [PATCH 1/2] Initial plan From 60a7b1db19c99d589c8bfaceaa3af3aee05409ba Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 24 Feb 2026 12:56:29 +0000 Subject: [PATCH 2/2] Add temperature primary variable with Newton update limiting for MSW wells Add enable_energy flag and Temperature index to MultisegmentWellPrimaryVariables. Limit temperature Newton update using max_temperature_change (default 5K), matching the reservoir model approach. Add temperature to SegmentState. Add necessary template instantiations for energy-enabled configurations. Co-authored-by: GitPaean <6847941+GitPaean@users.noreply.github.com> --- .../flow/BlackoilModelParameters.cpp | 4 +++ .../flow/BlackoilModelParameters.hpp | 6 ++++ opm/simulators/wells/MSWellHelpers.cpp | 6 +++- .../wells/MultisegmentWellAssemble.hpp | 2 +- .../wells/MultisegmentWellEquations.cpp | 4 ++- .../MultisegmentWellPrimaryVariables.cpp | 32 +++++++++++++++++-- .../MultisegmentWellPrimaryVariables.hpp | 10 ++++-- .../wells/MultisegmentWell_impl.hpp | 4 ++- opm/simulators/wells/SegmentState.cpp | 3 ++ opm/simulators/wells/SegmentState.hpp | 2 ++ opm/simulators/wells/WellAssemble.cpp | 1 + opm/simulators/wells/WellHelpers.cpp | 3 +- 12 files changed, 68 insertions(+), 9 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelParameters.cpp b/opm/simulators/flow/BlackoilModelParameters.cpp index c6b9c94bc11..60fba38aef1 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -61,6 +61,7 @@ BlackoilModelParameters::BlackoilModelParameters() relaxed_tolerance_flow_well_ = Parameters::Get>(); relaxed_tolerance_pressure_ms_well_ = Parameters::Get>(); max_pressure_change_ms_wells_ = Parameters::Get>(); + max_temperature_change_ms_wells_ = Parameters::Get>(); max_inner_iter_ms_wells_ = Parameters::Get(); strict_inner_iter_wells_ = Parameters::Get(); strict_outer_iter_wells_ = Parameters::Get(); @@ -197,6 +198,9 @@ void BlackoilModelParameters::registerParameters() Parameters::Register> ("Maximum relative pressure change for a single iteration " "of the multi-segment well model"); + Parameters::Register> + ("Maximum temperature change for a single iteration " + "of the multi-segment well model"); Parameters::Register ("Maximum number of inner iterations for multi-segment wells"); Parameters::Register diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index 9ac306abbcc..04fb84f5c9b 100644 --- a/opm/simulators/flow/BlackoilModelParameters.hpp +++ b/opm/simulators/flow/BlackoilModelParameters.hpp @@ -114,6 +114,9 @@ struct TolerancePressureMsWells { static constexpr Scalar value = 0.01*1e5; }; template struct MaxPressureChangeMsWells { static constexpr Scalar value = 10*1e5; }; +template +struct MaxTemperatureChangeMsWells { static constexpr Scalar value = 5.0; }; // Kelvin + struct MaxNewtonIterationsWithInnerWellIterations { static constexpr int value = 99; }; struct MaxInnerIterMsWells { static constexpr int value = 100; }; struct MaxInnerIterWells { static constexpr int value = 50; }; @@ -248,6 +251,9 @@ struct BlackoilModelParameters /// Maximum pressure change over an iteratio for ms wells Scalar max_pressure_change_ms_wells_; + /// Maximum temperature change over an iteration for ms wells + Scalar max_temperature_change_ms_wells_; + /// Maximum inner iteration number for ms wells int max_inner_iter_ms_wells_; diff --git a/opm/simulators/wells/MSWellHelpers.cpp b/opm/simulators/wells/MSWellHelpers.cpp index eb1eca9c636..29c6acb1520 100644 --- a/opm/simulators/wells/MSWellHelpers.cpp +++ b/opm/simulators/wells/MSWellHelpers.cpp @@ -376,7 +376,9 @@ using Mat = Dune::BCRSMatrix>; 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 applyUMFPack(Dune::UMFPack>&, \ @@ -414,6 +416,7 @@ using Mat = Dune::BCRSMatrix>; 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) \ @@ -421,6 +424,7 @@ using Mat = Dune::BCRSMatrix>; INSTANTIATE_EVAL(T,7) \ INSTANTIATE_EVAL(T,8) \ INSTANTIATE_EVAL(T,9) \ + INSTANTIATE_EVAL(T,10) \ INSTANTIATE_ALL_PARALLELLMSWELLB(T) INSTANTIATE_TYPE(double) diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index 7ae4594c89b..de60258d4c7 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp @@ -52,7 +52,7 @@ class MultisegmentWellAssemble static constexpr int SPres = PrimaryVariables::SPres; public: - static constexpr int numWellEq = Indices::numPhases+1; + static constexpr int numWellEq = PrimaryVariables::numWellEq; using Scalar = typename FluidSystem::Scalar; using IndexTraits = typename FluidSystem::IndexTraitsType; using Equations = MultisegmentWellEquations; diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index 1dee8d391d0..b59697a3bcb 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -487,7 +487,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) diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index b1b585d9126..ab4c052997c 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -157,6 +157,13 @@ update(const WellState& well_state, } } } + // if thermal is active, we set the temperature + if constexpr (enable_energy) { + const auto& segment_temperature = segments.temperature; + for (std::size_t seg = 0; seg < value_.size(); ++seg) { + value_[seg][Temperature] = segment_temperature[seg]; + } + } setEvaluationsFromValues(); } @@ -166,7 +173,8 @@ updateNewton(const BVectorWell& dwells, const Scalar relaxation_factor, const Scalar dFLimit, const bool stop_or_zero_rate_target, - const Scalar max_pressure_change) + const Scalar max_pressure_change, + const Scalar max_temperature_change) { const std::vector> old_primary_variables = value_; @@ -209,6 +217,13 @@ updateNewton(const BVectorWell& dwells, } } } + + // update the segment temperature + if constexpr (enable_energy) { + const int sign = dwells[seg][Temperature] > 0. ? 1 : -1; + const Scalar dx_limited = sign * std::min(std::abs(dwells[seg][Temperature]) * relaxation_factor, max_temperature_change); + value_[seg][Temperature] = std::max(old_primary_variables[seg][Temperature] - dx_limited, Scalar{0.0}); + } } if (stop_or_zero_rate_target) { @@ -300,13 +315,18 @@ copyToWellState(const MultisegmentWellGeneric& mswell, // update the segment pressure segment_pressure[seg] = value_[seg][SPres]; + // update the segment temperature + if constexpr (enable_energy) { + segments.temperature[seg] = value_[seg][Temperature]; + } + if (seg == 0) { // top segment ws.bhp = segment_pressure[seg]; } // Calculate other per-phase dynamic quantities. - const Scalar temperature = 0.0; // Ignore thermal effects + const Scalar temperature = enable_energy ? value_[seg][Temperature] : Scalar{0.0}; const Scalar saltConc = 0.0; // Ignore salt precipitation const Scalar Rvw = 0.0; // Ignore vaporised water. @@ -648,6 +668,14 @@ getSegmentPressure(const int seg) const return evaluation_[seg][SPres]; } +template +typename MultisegmentWellPrimaryVariables::EvalWell +MultisegmentWellPrimaryVariables:: +getSegmentTemperature(const int seg) const +{ + return evaluation_[seg][Temperature]; +} + template typename MultisegmentWellPrimaryVariables::EvalWell MultisegmentWellPrimaryVariables:: diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp index 590f87c7b83..9e2a89028dc 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp @@ -59,14 +59,16 @@ class MultisegmentWellPrimaryVariables static constexpr bool has_wfrac_variable = Indices::waterEnabled && Indices::oilEnabled; static constexpr bool has_gfrac_variable = Indices::gasEnabled && Indices::numPhases > 1; + static constexpr bool enable_energy = Indices::enableFullyImplicitThermal; static constexpr int WQTotal = 0; static constexpr int WFrac = has_wfrac_variable ? 1 : -1000; static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000; static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1; + static constexpr int Temperature = enable_energy ? SPres + 1 : -1000; // the number of well equations TODO: it should have a more general strategy for it - 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; @@ -91,7 +93,8 @@ class MultisegmentWellPrimaryVariables const Scalar relaxation_factor, const Scalar DFLimit, const bool stop_or_zero_rate_target, - const Scalar max_pressure_change); + const Scalar max_pressure_change, + const Scalar max_temperature_change); //! \brief Copy values to well state. void copyToWellState(const MultisegmentWellGeneric& mswell, @@ -121,6 +124,9 @@ class MultisegmentWellPrimaryVariables //! \brief Get pressure for a segment. EvalWell getSegmentPressure(const int seg) const; + //! \brief Get temperature for a segment. + EvalWell getSegmentTemperature(const int seg) const; + //! \brief Get rate for a component in a segment. EvalWell getSegmentRate(const int seg, const int comp_idx) const; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 89e305d1a17..8541e397436 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -728,13 +728,15 @@ namespace Opm const Scalar dFLimit = this->param_.dwell_fraction_max_; const Scalar max_pressure_change = this->param_.max_pressure_change_ms_wells_; + const Scalar max_temperature_change = this->param_.max_temperature_change_ms_wells_; const bool stop_or_zero_rate_target = this->stoppedOrZeroRateTarget(groupStateHelper); this->primary_variables_.updateNewton(dwells, relaxation_factor, dFLimit, stop_or_zero_rate_target, - max_pressure_change); + max_pressure_change, + max_temperature_change); const auto& summary_state = simulator.vanguard().summaryState(); this->primary_variables_.copyToWellState(*this, getRefDensity(), diff --git a/opm/simulators/wells/SegmentState.cpp b/opm/simulators/wells/SegmentState.cpp index 77d2f7d0c8a..bae23e27c1b 100644 --- a/opm/simulators/wells/SegmentState.cpp +++ b/opm/simulators/wells/SegmentState.cpp @@ -60,6 +60,7 @@ SegmentState::SegmentState(int num_phases, const WellSegments& segments) , phase_viscosity (segments.size() * num_phases) , phase_density (segments.size() * (num_phases + 2)) // +2 for mixture with and without exponents , pressure (segments.size()) + , temperature (segments.size()) , pressure_drop_friction (segments.size()) , pressure_drop_hydrostatic(segments.size()) , pressure_drop_accel (segments.size()) @@ -79,6 +80,7 @@ SegmentState SegmentState::serializationTestObject() result.phase_viscosity = {12.0, 12.5}; result.phase_density = {13.0, 13.5, 13.6, 13.75}; result.pressure = {14.0, 15.0}; + result.temperature = {14.5, 15.5}; result.pressure_drop_friction = {16.0}; result.pressure_drop_hydrostatic = {17.0, 18.0}; result.pressure_drop_accel = {19.0}; @@ -136,6 +138,7 @@ bool SegmentState::operator==(const SegmentState& rhs) const this->phase_viscosity == rhs.phase_viscosity && this->phase_density == rhs.phase_density && this->pressure == rhs.pressure && + this->temperature == rhs.temperature && this->pressure_drop_friction == rhs.pressure_drop_friction && this->pressure_drop_hydrostatic == rhs.pressure_drop_hydrostatic && this->pressure_drop_accel == rhs.pressure_drop_accel && diff --git a/opm/simulators/wells/SegmentState.hpp b/opm/simulators/wells/SegmentState.hpp index 7c6724bb1f1..78ab4a4d899 100644 --- a/opm/simulators/wells/SegmentState.hpp +++ b/opm/simulators/wells/SegmentState.hpp @@ -57,6 +57,7 @@ class SegmentState serializer(phase_viscosity); serializer(phase_density); serializer(pressure); + serializer(temperature); serializer(pressure_drop_friction); serializer(pressure_drop_hydrostatic); serializer(pressure_drop_accel); @@ -100,6 +101,7 @@ class SegmentState std::vector phase_density; std::vector pressure; + std::vector temperature; std::vector pressure_drop_friction; std::vector pressure_drop_hydrostatic; std::vector pressure_drop_accel; diff --git a/opm/simulators/wells/WellAssemble.cpp b/opm/simulators/wells/WellAssemble.cpp index 164f503cb8b..276850175b3 100644 --- a/opm/simulators/wells/WellAssemble.cpp +++ b/opm/simulators/wells/WellAssemble.cpp @@ -311,6 +311,7 @@ using FS = BlackOilFluidSystem; INSTANTIATE_METHODS(FS, DenseAd::Evaluation) \ INSTANTIATE_METHODS(FS, DenseAd::Evaluation) \ INSTANTIATE_METHODS(FS, DenseAd::Evaluation) \ + INSTANTIATE_METHODS(FS, DenseAd::Evaluation) \ INSTANTIATE_METHODS(FS, DenseAd::Evaluation) \ INSTANTIATE_METHODS(FS, DenseAd::Evaluation) \ INSTANTIATE_METHODS(FS, DenseAd::Evaluation) \ diff --git a/opm/simulators/wells/WellHelpers.cpp b/opm/simulators/wells/WellHelpers.cpp index ab2e5673c9b..a608286026c 100644 --- a/opm/simulators/wells/WellHelpers.cpp +++ b/opm/simulators/wells/WellHelpers.cpp @@ -252,7 +252,8 @@ using Comm = Parallel::Communication; INSTANTIATE(T,6) \ INSTANTIATE_WE(T,2) \ INSTANTIATE_WE(T,3) \ - INSTANTIATE_WE(T,4) + INSTANTIATE_WE(T,4) \ + INSTANTIATE_WE(T,5) INSTANTIATE_TYPE(double)