Skip to content

Limit temperature Newton update for multisegment wells #6876

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
Copilot wants to merge 2 commits into from
Closed

Limit temperature Newton update for multisegment wells #6876

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
294466e
Initial plan
Copilot Feb 24, 2026
60a7b1d
Add temperature primary variable with Newton update limiting for MSW …
Copilot Feb 24, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions opm/simulators/flow/BlackoilModelParameters.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
relaxed_tolerance_flow_well_ = Parameters::Get<Parameters::RelaxedWellFlowTol<Scalar>>();
relaxed_tolerance_pressure_ms_well_ = Parameters::Get<Parameters::RelaxedPressureTolMsw<Scalar>>();
max_pressure_change_ms_wells_ = Parameters::Get<Parameters::MaxPressureChangeMsWells<Scalar>>();
max_temperature_change_ms_wells_ = Parameters::Get<Parameters::MaxTemperatureChangeMsWells<Scalar>>();
max_inner_iter_ms_wells_ = Parameters::Get<Parameters::MaxInnerIterMsWells>();
strict_inner_iter_wells_ = Parameters::Get<Parameters::StrictInnerIterWells>();
strict_outer_iter_wells_ = Parameters::Get<Parameters::StrictOuterIterWells>();
Expand Down Expand Up @@ -197,6 +198,9 @@ void BlackoilModelParameters<Scalar>::registerParameters()
Parameters::Register<Parameters::MaxPressureChangeMsWells<Scalar>>
("Maximum relative pressure change for a single iteration "
"of the multi-segment well model");
Parameters::Register<Parameters::MaxTemperatureChangeMsWells<Scalar>>
("Maximum temperature change for a single iteration "
"of the multi-segment well model");
Parameters::Register<Parameters::MaxInnerIterMsWells>
("Maximum number of inner iterations for multi-segment wells");
Parameters::Register<Parameters::StrictInnerIterWells>
Expand Down
6 changes: 6 additions & 0 deletions opm/simulators/flow/BlackoilModelParameters.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,9 @@ struct TolerancePressureMsWells { static constexpr Scalar value = 0.01*1e5; };
template<class Scalar>
struct MaxPressureChangeMsWells { static constexpr Scalar value = 10*1e5; };

template<class Scalar>
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; };
Expand Down Expand Up @@ -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_;

Expand Down
6 changes: 5 additions & 1 deletion opm/simulators/wells/MSWellHelpers.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>>&, \
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion opm/simulators/wells/MultisegmentWellAssemble.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Scalar, IndexTraits, numWellEq,Indices::numEq>;
Expand Down
4 changes: 3 additions & 1 deletion opm/simulators/wells/MultisegmentWellEquations.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
32 changes: 30 additions & 2 deletions opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,6 +157,13 @@ update(const WellState<Scalar, IndexTraits>& 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();
}

Expand All @@ -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<std::array<Scalar, numWellEq>> old_primary_variables = value_;

Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -300,13 +315,18 @@ copyToWellState(const MultisegmentWellGeneric<Scalar, IndexTraits>& 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.

Expand Down Expand Up @@ -648,6 +668,14 @@ getSegmentPressure(const int seg) const
return evaluation_[seg][SPres];
}

template<class FluidSystem, class Indices>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
getSegmentTemperature(const int seg) const
{
return evaluation_[seg][Temperature];
}

template<class FluidSystem, class Indices>
typename MultisegmentWellPrimaryVariables<FluidSystem,Indices>::EvalWell
MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
Expand Down
10 changes: 8 additions & 2 deletions opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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<Scalar, IndexTraits>& mswell,
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 3 additions & 1 deletion opm/simulators/wells/MultisegmentWell_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(),
Expand Down
3 changes: 3 additions & 0 deletions opm/simulators/wells/SegmentState.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ SegmentState<Scalar>::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())
Expand All @@ -79,6 +80,7 @@ SegmentState<Scalar> SegmentState<Scalar>::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};
Expand Down Expand Up @@ -136,6 +138,7 @@ bool SegmentState<Scalar>::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 &&
Expand Down
2 changes: 2 additions & 0 deletions opm/simulators/wells/SegmentState.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -100,6 +101,7 @@ class SegmentState
std::vector<Scalar> phase_density;

std::vector<Scalar> pressure;
std::vector<Scalar> temperature;
std::vector<Scalar> pressure_drop_friction;
std::vector<Scalar> pressure_drop_hydrostatic;
std::vector<Scalar> pressure_drop_accel;
Expand Down
1 change: 1 addition & 0 deletions opm/simulators/wells/WellAssemble.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -311,6 +311,7 @@ using FS = BlackOilFluidSystem<Scalar, BlackOilDefaultFluidSystemIndices>;
INSTANTIATE_METHODS(FS<T>, DenseAd::Evaluation<T,7,0u>) \
INSTANTIATE_METHODS(FS<T>, DenseAd::Evaluation<T,8,0u>) \
INSTANTIATE_METHODS(FS<T>, DenseAd::Evaluation<T,9,0u>) \
INSTANTIATE_METHODS(FS<T>, DenseAd::Evaluation<T,10,0u>) \
INSTANTIATE_METHODS(FS<T>, DenseAd::Evaluation<T,-1,4u>) \
INSTANTIATE_METHODS(FS<T>, DenseAd::Evaluation<T,-1,5u>) \
INSTANTIATE_METHODS(FS<T>, DenseAd::Evaluation<T,-1,6u>) \
Expand Down
3 changes: 2 additions & 1 deletion opm/simulators/wells/WellHelpers.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down