From e5a1729d197a995f752e3e6c264a2f7022d86b2e Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 29 Apr 2026 21:46:41 +0200 Subject: [PATCH 01/11] adding temperature to SegmentState and PerfData --- opm/simulators/wells/PerfData.cpp | 4 +++ opm/simulators/wells/PerfData.hpp | 2 ++ opm/simulators/wells/SegmentState.cpp | 3 ++ opm/simulators/wells/SegmentState.hpp | 3 ++ opm/simulators/wells/WellState.cpp | 44 +++++++++++++++++---------- opm/simulators/wells/WellState.hpp | 1 + 6 files changed, 41 insertions(+), 16 deletions(-) diff --git a/opm/simulators/wells/PerfData.cpp b/opm/simulators/wells/PerfData.cpp index ff2765860e4..a49de4a257b 100644 --- a/opm/simulators/wells/PerfData.cpp +++ b/opm/simulators/wells/PerfData.cpp @@ -33,6 +33,7 @@ PerfData::PerfData(const std::size_t num_perf, const std::size_t num_phases) : injector(injector_) , pressure(num_perf) + , temperature(num_perf) , rates(num_perf) , phase_rates(num_perf * num_phases) , phase_mixing_rates(num_perf) @@ -73,6 +74,7 @@ PerfData PerfData::serializationTestObject() PerfData result; result.injector = true; result.pressure = {2.0, 3.0, 4.0}; + result.temperature = {345.}; result.rates = {5.0, 6.0}; result.phase_rates = {7.0}; result.phase_mixing_rates = { {1.0, 2.0, 3.0, 4.0}}; @@ -124,6 +126,7 @@ bool PerfData::try_assign(const PerfData& other) } this->pressure = other.pressure; + this->temperature = other.temperature; this->rates = other.rates; this->phase_rates = other.phase_rates; this->phase_mixing_rates = other.phase_mixing_rates; @@ -150,6 +153,7 @@ bool PerfData::operator==(const PerfData& rhs) const { return (this->injector == rhs.injector) && (this->pressure == rhs.pressure) + && (this->temperature == rhs.temperature) && (this->rates == rhs.rates) && (this->phase_rates == rhs.phase_rates) && (this->phase_mixing_rates == rhs.phase_mixing_rates) diff --git a/opm/simulators/wells/PerfData.hpp b/opm/simulators/wells/PerfData.hpp index 86e5f1c9b89..9ef799610e7 100644 --- a/opm/simulators/wells/PerfData.hpp +++ b/opm/simulators/wells/PerfData.hpp @@ -57,6 +57,7 @@ class PerfData { serializer(injector); serializer(pressure); + serializer(temperature); serializer(rates); serializer(phase_rates); serializer(phase_mixing_rates); @@ -91,6 +92,7 @@ class PerfData // result, e.g., a flow rate, then please update try_assign() as well. std::vector pressure{}; + std::vector temperature{}; std::vector rates{}; std::vector phase_rates{}; std::vector> phase_mixing_rates{}; diff --git a/opm/simulators/wells/SegmentState.cpp b/opm/simulators/wells/SegmentState.cpp index 77d2f7d0c8a..9f1e62e2e92 100644 --- a/opm/simulators/wells/SegmentState.cpp +++ b/opm/simulators/wells/SegmentState.cpp @@ -63,6 +63,7 @@ SegmentState::SegmentState(int num_phases, const WellSegments& segments) , pressure_drop_friction (segments.size()) , pressure_drop_hydrostatic(segments.size()) , pressure_drop_accel (segments.size()) + , temperature (segments.size()) , m_segment_number (make_segment_number(segments)) {} @@ -82,6 +83,7 @@ SegmentState SegmentState::serializationTestObject() result.pressure_drop_friction = {16.0}; result.pressure_drop_hydrostatic = {17.0, 18.0}; result.pressure_drop_accel = {19.0}; + result.temperature = {20.0, 21.0}; result.m_segment_number = {20, 21}; return result; @@ -139,6 +141,7 @@ bool SegmentState::operator==(const SegmentState& rhs) const this->pressure_drop_friction == rhs.pressure_drop_friction && this->pressure_drop_hydrostatic == rhs.pressure_drop_hydrostatic && this->pressure_drop_accel == rhs.pressure_drop_accel && + this->temperature == rhs.temperature && this->m_segment_number == rhs.m_segment_number; } diff --git a/opm/simulators/wells/SegmentState.hpp b/opm/simulators/wells/SegmentState.hpp index 7c6724bb1f1..2f833a479b5 100644 --- a/opm/simulators/wells/SegmentState.hpp +++ b/opm/simulators/wells/SegmentState.hpp @@ -60,6 +60,7 @@ class SegmentState serializer(pressure_drop_friction); serializer(pressure_drop_hydrostatic); serializer(pressure_drop_accel); + serializer(temperature); serializer(m_segment_number); } @@ -104,6 +105,8 @@ class SegmentState std::vector pressure_drop_hydrostatic; std::vector pressure_drop_accel; + std::vector temperature; + private: std::vector m_segment_number; }; diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index a22f5a50e40..029a227821f 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -188,11 +188,10 @@ void WellState:: initSingleProducer(const Well& well, const ParallelWellInfo& well_info, Scalar pressure_first_connection, + Scalar temperature_first_connection, const std::vector>& well_perf_data, const SummaryState& summary_state) { - const Scalar temp = 273.15 + 15.56; - auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, @@ -200,7 +199,7 @@ initSingleProducer(const Well& well, true, pressure_first_connection, well_perf_data, - temp}); + temperature_first_connection}); // the rest of the code needs to executed even if ws.perf_data is empty // as this does not say anything for the whole well if it is distributed. @@ -221,14 +220,14 @@ initSingleInjector(const Well& well, const std::vector>& well_perf_data, const SummaryState& summary_state) { - const Scalar temp = well.hasInjTemperature() ? well.inj_temperature() : temperature_first_connection; + const Scalar temperature = well.hasInjTemperature() ? well.inj_temperature() : temperature_first_connection; auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, this->phaseUsageInfo(), false, pressure_first_connection, well_perf_data, - temp}); + temperature}); // the rest of the code needs to executed even if ws.perf_data is empty // as this does not say anything for the whole well if it is distributed. @@ -250,23 +249,21 @@ initSingleWell(const std::vector& cellPressures, const SummaryState& summary_state) { Scalar pressure_first_connection = -1; + Scalar temperature_first_connection = -1; if (!well_perf_data.empty()) { - pressure_first_connection = cellPressures[well_perf_data[0].cell_index]; + const auto first_cell_index = well_perf_data[0].cell_index; + pressure_first_connection = cellPressures[first_cell_index]; + temperature_first_connection = cellTemperatures[first_cell_index]; } // The following call is necessary to ensure that processes that do not contain the first perforation get the correct value pressure_first_connection = well_info.broadcastFirstPerforationValue(pressure_first_connection); + temperature_first_connection = well_info.broadcastFirstPerforationValue(temperature_first_connection); if (well.isInjector()) { - Scalar temperature_first_connection = -1; - if (!well_perf_data.empty()) { - temperature_first_connection = cellTemperatures[well_perf_data[0].cell_index]; - } - // The following call is necessary to ensure that processes that do not contain the first perforation get the correct value - temperature_first_connection = well_info.broadcastFirstPerforationValue(temperature_first_connection); this->initSingleInjector(well, well_info, pressure_first_connection, temperature_first_connection, well_perf_data, summary_state); } else { - this->initSingleProducer(well, well_info, pressure_first_connection, + this->initSingleProducer(well, well_info, pressure_first_connection, temperature_first_connection, well_perf_data, summary_state); } } @@ -338,7 +335,9 @@ init(const std::vector& cellPressures, perf_data.phase_rates[this->numPhases()*perf + p] = ws.surface_rates[p] / Scalar(global_num_perf_this_well); } } - perf_data.pressure[perf] = cellPressures[well_perf_data[w][perf].cell_index]; + const auto cell_index = well_perf_data[w][perf].cell_index; + perf_data.pressure[perf] = cellPressures[cell_index]; + perf_data.temperature[perf] = cellTemperatures[cell_index]; } } @@ -822,6 +821,7 @@ initWellStateMSWell(const std::vector& wells_ecl, const auto& perf_rates = perf_data.phase_rates; const auto& perf_press = perf_data.pressure; + const auto& perf_temperature = perf_data.temperature; // The function calculateSegmentRates as well as the loop filling the segment_pressure work // with *global* containers. Now we create global vectors containing the phase_rates and // pressures of all processes. @@ -835,6 +835,7 @@ initWellStateMSWell(const std::vector& wells_ecl, std::vector perforation_rates(number_of_global_perfs * np, 0.0); std::vector perforation_pressures(number_of_global_perfs, 0.0); + std::vector perforation_temperature(number_of_global_perfs, 0.0); assert(perf_data.size() == perf_press.size()); assert(perf_data.size() * np == perf_rates.size()); @@ -842,6 +843,7 @@ initWellStateMSWell(const std::vector& wells_ecl, if (auto candidate = active_perf_index_local_to_global.find(perf); candidate != active_perf_index_local_to_global.end()) { const int global_active_perf_index = candidate->second; perforation_pressures[global_active_perf_index] = perf_press[perf]; + perforation_temperature[global_active_perf_index] = perf_temperature[perf]; for (int i = 0; i < np; i++) { perforation_rates[global_active_perf_index * np + i] = perf_rates[perf * np + i]; } @@ -863,21 +865,31 @@ initWellStateMSWell(const std::vector& wells_ecl, // for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment // if there is no perforation associated with this segment, it uses the pressure from the outlet segment - // which requres the ordering is successful + // which requires the ordering is successful // Not sure what is the best way to handle the initialization, hopefully, the bad initialization can be // improved during the solveWellEq process { // top segment is always the first one, and its pressure is the well bhp auto& segment_pressure = ws.segments.pressure; + auto& segment_temperature = ws.segments.temperature; segment_pressure[0] = ws.bhp; + // TODO: the temperature remains to be made to be more consistent + // TODO: we probably should distinguish whether it is a thermal case to avoid having garbage data in the segment temperature + segment_temperature[0] = ws.temperature; for (int seg = 1; seg < well_nseg; ++seg) { if (!segment_perforations[seg].empty()) { const int first_perf_global_index = segment_perforations[seg][0]; segment_pressure[seg] = perforation_pressures[first_perf_global_index]; + // TODO: the temperature remains to be made to be more consistent + const int outlet_seg = segment_set[seg].outletSegment(); + const int outlet_seg_index = segment_set.segmentNumberToIndex(outlet_seg); + segment_temperature[seg] = ws.producer ? perforation_temperature[first_perf_global_index] : segment_temperature[outlet_seg_index]; } else { // using the outlet segment pressure // it needs the ordering is correct const int outlet_seg = segment_set[seg].outletSegment(); - segment_pressure[seg] = segment_pressure[segment_set.segmentNumberToIndex(outlet_seg)]; + const int outlet_seg_index = segment_set.segmentNumberToIndex(outlet_seg); + segment_pressure[seg] = segment_pressure[outlet_seg_index]; + segment_temperature[seg] = segment_temperature[outlet_seg_index]; } } } diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 243e2a29092..3ac7a333ce1 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -408,6 +408,7 @@ class WellState void initSingleProducer(const Well& well, const ParallelWellInfo& well_info, Scalar pressure_first_connection, + Scalar temperature_first_connection, const std::vector>& well_perf_data, const SummaryState& summary_state); From 7cf7be1506ecd5ce8db7761d6c86e21d44ec4ddc Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 29 Apr 2026 21:48:07 +0200 Subject: [PATCH 02/11] adding Energy to WellFailure --- opm/simulators/timestepping/ConvergenceReport.cpp | 2 ++ opm/simulators/timestepping/ConvergenceReport.hpp | 1 + 2 files changed, 3 insertions(+) diff --git a/opm/simulators/timestepping/ConvergenceReport.cpp b/opm/simulators/timestepping/ConvergenceReport.cpp index ad439c90e5f..9bf81495b79 100644 --- a/opm/simulators/timestepping/ConvergenceReport.cpp +++ b/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: diff --git a/opm/simulators/timestepping/ConvergenceReport.hpp b/opm/simulators/timestepping/ConvergenceReport.hpp index 563905738ab..1fc96a70fdc 100644 --- a/opm/simulators/timestepping/ConvergenceReport.hpp +++ b/opm/simulators/timestepping/ConvergenceReport.hpp @@ -164,6 +164,7 @@ namespace Opm Invalid, MassBalance, Pressure, + Energy, ControlBHP, ControlTHP, ControlRate, From d4b7c820eaa963c29d176770e6d3146f8157d696 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 29 Apr 2026 21:49:16 +0200 Subject: [PATCH 03/11] adding energy equation for multisegment wells including adding tempearture as the primary variable for the MSW. to calculate the enthalpy, for each segment, we generate a SegmentFluidState. --- .../wells/BlackoilWellModel_impl.hpp | 6 + opm/simulators/wells/MSWellHelpers.cpp | 6 +- opm/simulators/wells/MultisegmentWell.hpp | 85 +++- .../wells/MultisegmentWellAssemble.cpp | 12 + .../wells/MultisegmentWellAssemble.hpp | 4 +- .../wells/MultisegmentWellEquations.cpp | 4 +- opm/simulators/wells/MultisegmentWellEval.cpp | 40 +- opm/simulators/wells/MultisegmentWellEval.hpp | 3 + .../MultisegmentWellPrimaryVariables.cpp | 34 ++ .../MultisegmentWellPrimaryVariables.hpp | 25 +- .../wells/MultisegmentWellSegments.cpp | 1 + .../wells/MultisegmentWell_impl.hpp | 462 +++++++++++++++++- opm/simulators/wells/WellAssemble.cpp | 1 + opm/simulators/wells/WellHelpers.cpp | 3 +- opm/simulators/wells/WellInterfaceGeneric.hpp | 2 + 15 files changed, 664 insertions(+), 24 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index b3a8364ae44..19256a45200 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -2159,6 +2159,12 @@ namespace Opm { int BlackoilWellModel::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 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/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 1e41a3b1d70..eeee33cbe5f 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -27,6 +27,8 @@ #include #include +#include + namespace Opm { class DeferredLogger; @@ -56,6 +58,7 @@ namespace Opm { using Base::has_solvent; using Base::has_polymer; + using Base::has_energy; using Base::Water; using Base::Oil; using Base::Gas; @@ -73,6 +76,25 @@ namespace Opm { using typename Base::PressureMatrix; using FSInfo = std::tuple; + // 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 + // templated on Value type because we want to use it for surface condition calculation injection + // while we shifted to use the bottom-hole/well-head condition instead, so there is only EvalWell type + // used here for now. + template + using SegmentFluidState = BlackOilFluidState= 0, + /*has_watVapor*/ false, + /*has_brine*/ false, + /*has_saltPrecip*/ false, + /*has_disgas_in_water*/ false, + has_solvent, + Indices::numPhases>; + MultisegmentWell(const Well& well, const ParallelWellInfo& pw_info, const int time_step, @@ -173,6 +195,16 @@ namespace Opm { // the intial amount of fluids in each segment under surface condition std::vector > segment_fluid_initial_; + // total energy inside the segments at the beginning of the time step + std::vector segment_initial_energy_; + + // segment fluid state + std::vector> 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 wellhead_fluid_state_; mutable int debug_cost_counter_ = 0; @@ -332,9 +364,60 @@ 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 + // TODO: this function largely overlaps with calculatePhaseProperties(), some refactoring/unification should be done + template + SegmentFluidState + createFluidState(const std::vector& fluid_composition, + const ValueType& pressure, + const ValueType& temperature, + DeferredLogger& deferred_logger) const; + + SegmentFluidState + createSegmentFluidState(int seg, DeferredLogger& deferred_logger) const; + + void computeInitialSegmentEnergy(); + + // assemble the energy equation contribution for a single perforation/connection + void assemblePerforationEnergyEq(const IntensiveQuantities& int_quants, + const std::vector& cq_s, + const int seg, + const int local_perf_index, + DeferredLogger& deferred_logger); + + void updateWellHeadCondition(const Simulator& simulator, DeferredLogger& deferred_logger); + + void updateSegmentFluidState(DeferredLogger& deferred_logger); + + template + ValueType computeSegmentEnergy(int seg) const; + + // Convert per-component surface volumetric rates to a phase reservoir + // volumetric rate using the upwind segment 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. + // @return the reservoir volumetric rate of @p phaseIdx. + EvalWell surfaceToReservoirRate(unsigned phaseIdx, + const SegmentFluidState& segment_fs, + const std::vector& 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& upwind_fs, + std::string_view context, + DeferredLogger& deferred_logger) const; }; -} +} // namespace Opm #include "MultisegmentWell_impl.hpp" diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp index a4f2f2844bb..b14bf9092a2 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.cpp +++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp @@ -352,6 +352,12 @@ assembleOutflowTerm(const int seg, eqns.D()[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq); } // pressure derivative should be zero + + if constexpr (enable_energy) { + // energy flux depends on the pressure + 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 @@ -377,6 +383,12 @@ assembleInflowTerm(const int seg, eqns.D()[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq); } // pressure derivative should be zero + + if constexpr (enable_energy) { + // energy flux depends on the pressure + 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 diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp index 7ae4594c89b..b57f046d39d 100644 --- a/opm/simulators/wells/MultisegmentWellAssemble.hpp +++ b/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; diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index def1efa40f2..52e2cdcb0dd 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/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) diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 580ca59f652..bb52b0263ad 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -89,7 +89,11 @@ getWellConvergence(const WellState& 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> abs_residual(this->numberOfSegments(), @@ -103,7 +107,7 @@ getWellConvergence(const WellState& well_state, std::vector 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 @@ -111,6 +115,11 @@ getWellConvergence(const WellState& well_state, 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,22 @@ getWellConvergence(const WellState& 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; + constexpr Scalar tolerance_energy_wells = 1.e-5; // TODO: need to determine a reasonable value for this + constexpr Scalar relaxed_tolerance_energy_wells = 1.e-4; // TODO: need to determine a reasonable value for this + if (std::isnan(energy_residual)) { + report.setWellFailed({CR::WellFailure::Type::Energy, CR::Severity::NotANumber, dummy_phase, baseif_.name()}); + } else if (std::isinf(energy_residual)) { + report.setWellFailed({CR::WellFailure::Type::Energy, CR::Severity::TooLarge, dummy_phase, baseif_.name()}); + } else if (!relax_tolerance && energy_residual > tolerance_energy_wells) { + report.setWellFailed({CR::WellFailure::Type::Energy, CR::Severity::Normal, dummy_phase, baseif_.name()}); + } else if (energy_residual > relaxed_tolerance_energy_wells) { + 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 +461,10 @@ getFiniteWellResiduals(const std::vector& 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 +575,10 @@ getResidualMeasureValue(const WellState& 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; diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index c5965b85d93..1d2be09f6da 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -54,6 +54,9 @@ class MultisegmentWellEval : public MultisegmentWellGeneric; using MSWSegments = MultisegmentWellSegments; diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 097f7fa819a..7b0f7f513e6 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -157,6 +157,14 @@ 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(); } @@ -209,6 +217,15 @@ updateNewton(const BVectorWell& dwells, } } } + + if constexpr (enable_energy) { + // TODO: how to regularize the temperature update remains to be investigated + const int sign = dwells[seg][Temperature] > 0. ? 1 : -1; + // currently we limit the temperature change to be 5 C/K in each iteration + constexpr Scalar max_temperature_change = 5.0; + 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) { @@ -242,6 +259,7 @@ copyToWellState(const MultisegmentWellGeneric& mswell, auto& disgas = segments.dissolved_gas_rate; auto& vapoil = segments.vaporized_oil_rate; auto& segment_pressure = segments.pressure; + auto& segment_temperature = segments.temperature; for (std::size_t seg = 0; seg < value_.size(); ++seg) { std::vector fractions(well_.numPhases(), 0.0); @@ -300,6 +318,14 @@ copyToWellState(const MultisegmentWellGeneric& mswell, // update the segment pressure segment_pressure[seg] = value_[seg][SPres]; + // update the segment temperature if thermal is active + if (enable_energy) { + segment_temperature[seg] = value_[seg][Temperature]; + if (seg == 0) { + ws.temperature = segment_temperature[seg]; + } + } + if (seg == 0) { // top segment ws.bhp = segment_pressure[seg]; } @@ -648,6 +674,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 bdd646d3b04..89bff6c1e8e 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp @@ -51,28 +51,34 @@ class MultisegmentWellPrimaryVariables // Table showing the primary variable indices, depending on what phases are present: // - // WOG OG WG WO W/O/G (single phase) - // WQTotal 0 0 0 0 0 - // WFrac 1 -1000 -1000 1 -1000 - // GFrac 2 1 1 -1000 -1000 - // Spres 3 2 2 2 1 + // WQTotal 0 0 0 0 0 + // WFrac 1 -1000 -1000 1 -1000 + // GFrac 2 1 1 -1000 -1000 + // Temperature(*) 3 2 2 2 1 + // SPres 3/4 2/3 2/3 2/3 1/2 + // (*) Temperature is only present when energy is enabled; + // SPres shifts by one when Temperature is present. 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; + // Temperature comes before SPres: it is a conservation equation consistent + // with reservoir equations, while SPres is well-specific. + static constexpr int Temperature = enable_energy ? has_wfrac_variable + has_gfrac_variable + 1 : -1000; + static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1 + enable_energy; // 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; using EvalWell = DenseAd::Evaluation; - using Equations = MultisegmentWellEquations; + using Equations = MultisegmentWellEquations; using BVectorWell = typename Equations::BVectorWell; explicit MultisegmentWellPrimaryVariables(const WellInterfaceIndices& well) @@ -121,6 +127,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/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 3cece323a68..407e0a5dcd7 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -196,6 +196,7 @@ computeFluidProperties(const EvalWell& temperature, } densities_[seg] = density / volrat; + // TODO: the enthalpy flux rate should be calculated in the similar manner. // calculate the mass rates mass_rates_[seg] = 0.; for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) { diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index a54f70f60b4..ff794ec1adf 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -47,6 +47,7 @@ #include #include +#include #include #if COMPILE_GPU_BRIDGE && (HAVE_CUDA || HAVE_OPENCL) @@ -73,6 +74,8 @@ namespace Opm , MSWEval(static_cast&>(*this), pw_info) , regularize_(false) , segment_fluid_initial_(this->numberOfSegments(), std::vector(this->num_conservation_quantities_, 0.0)) + , segment_initial_energy_(this->numberOfSegments(), 0.0) + , segment_fluid_state_(this->numberOfSegments(), SegmentFluidState{}) { // not handling solvent or polymer for now with multisegment well if constexpr (has_solvent) { @@ -83,10 +86,6 @@ namespace Opm OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet"); } - if constexpr (Base::has_energy) { - OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet"); - } - if constexpr (Base::has_foam) { OPM_THROW(std::runtime_error, "foam is not supported by multisegment well yet"); } @@ -766,7 +765,14 @@ namespace Opm auto& deferred_logger = groupStateHelper.deferredLogger(); updatePrimaryVariables(groupStateHelper); computePerfCellPressDiffs(simulator); + + // after updating the primary variables, we need to update the segment fluid state + updateSegmentFluidState(deferred_logger); + computeInitialSegmentFluids(simulator, deferred_logger); + if constexpr (has_energy) { + computeInitialSegmentEnergy(); + } } @@ -1785,6 +1791,7 @@ namespace Opm try{ const BVectorWell dx_well = this->linSys_.solve(); updateWellState(simulator, dx_well, groupStateHelper, well_state, relaxation_factor); + updateSegmentFluidState(deferred_logger); } catch(const NumericalProblem& exp) { // Add information about the well and log to deferred logger @@ -1921,6 +1928,11 @@ namespace Opm MultisegmentWellAssemble(*this). assemblePerforationEq(seg, local_perf_index, comp_idx, cq_s_effective, this->linSys_); } + + // assembling the energy equation for the perforation if needed + if constexpr (has_energy) { + assemblePerforationEnergyEq(int_quants, cq_s, seg, local_perf_index, deferred_logger); + } } } // Accumulate dissolved gas and vaporized oil flow rates across all ranks sharing this well. @@ -1933,9 +1945,9 @@ namespace Opm // accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed) this->linSys_.sumDistributed(this->parallel_well_info_.communication()); } + for (int seg = 0; seg < nseg; ++seg) { // calculating the accumulation term - // TODO: without considering the efficiency factor for now { const EvalWell segment_surface_volume = getSegmentSurfaceVolume(simulator, seg, deferred_logger); @@ -1950,6 +1962,13 @@ namespace Opm MultisegmentWellAssemble(*this). assembleAccumulationTerm(seg, comp_idx, accumulation_term, this->linSys_); } + + if constexpr (has_energy) { + const EvalWell segment_energy = this->computeSegmentEnergy(seg); + const EvalWell accumulation_term_energy = regularization_factor * (segment_energy - segment_initial_energy_[seg]) / dt; + MultisegmentWellAssemble(*this). + assembleAccumulationTerm(seg, MSWEval::PrimaryVariables::Temperature, accumulation_term_energy, this->linSys_); + } } // considering the contributions due to flowing out from the segment { @@ -1963,6 +1982,25 @@ namespace Opm MultisegmentWellAssemble(*this). assembleOutflowTerm(seg, seg_upwind, comp_idx, segment_rate, this->linSys_); } + if constexpr (has_energy) { + const bool top_injecting_segment = (seg == 0) && this->isInjector(); + if (top_injecting_segment) { + this->updateWellHeadCondition(simulator, deferred_logger); + } + + // Energy carried by fluid flowing out of this segment toward its outlet. + // Use upwind segment fluid properties for enthalpy and density. For the + // top injecting segment, use the wellhead fluid state instead. + const auto& upwind_fs = top_injecting_segment ? this->wellhead_fluid_state_ + : this->segment_fluid_state_[seg_upwind]; + assert((top_injecting_segment && seg_upwind == 0) || !top_injecting_segment); + + const EvalWell energy_rate = + this->computeSegmentEnergyRate(seg, seg_upwind, upwind_fs, + "energy outflow assembly", deferred_logger); + MultisegmentWellAssemble(*this). + assembleOutflowTerm(seg, seg_upwind, MSWEval::PrimaryVariables::Temperature, energy_rate, this->linSys_); + } } // considering the contributions from the inlet segments @@ -1979,6 +2017,19 @@ namespace Opm assembleInflowTerm(seg, inlet, inlet_upwind, comp_idx, inlet_rate, this->linSys_); } } + if constexpr (has_energy) { + for (const int inlet : this->segments_.inlets()[seg]) { + const int inlet_upwind = this->segments_.upwinding_segment(inlet); + // Energy carried by fluid flowing from inlet segment into this segment. + // Use upwind segment fluid properties for enthalpy and density. + const auto& upwind_fs = this->segment_fluid_state_[inlet_upwind]; + const EvalWell energy_rate = + this->computeSegmentEnergyRate(inlet, inlet_upwind, upwind_fs, + "energy inflow assembly", deferred_logger); + MultisegmentWellAssemble(*this). + assembleInflowTerm(seg, inlet, inlet_upwind, MSWEval::PrimaryVariables::Temperature, energy_rate, this->linSys_); + } + } } // the fourth equation, the pressure drop equation @@ -1986,7 +2037,7 @@ namespace Opm const bool stopped_or_zero_target = this->stoppedOrZeroRateTarget(groupStateHelper); // When solving with zero rate (well isolation), use empty group_state to isolate // from group constraints in assembly. - // Otherwise use real group state from groupStateHelper. + // Otherwise, use real group state from groupStateHelper. GroupState empty_group_state; // Note: Cannot use 'const auto&' here because pushGroupState() requires a // non-const reference. GroupStateHelper stores a non-const pointer to GroupState @@ -2392,6 +2443,405 @@ namespace Opm return this->parallel_well_info_.communication().size() == 1 ? info : this->parallel_well_info_.broadcastFirstPerforationValue(info); } + template + template + MultisegmentWell::SegmentFluidState + MultisegmentWell:: + createFluidState(const std::vector& fluid_composition, + const ValueType& pressure, + const ValueType& temperature, + DeferredLogger& deferred_logger) const + { + SegmentFluidState fluid_state; + if constexpr (has_energy) { + fluid_state.setTemperature(temperature); + } + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + // we assume there is no capillary pressure in the wellbore + fluid_state.setPressure(phaseIdx, pressure); + } + fluid_state.setPvtRegionIndex(this->pvtRegionIdx()); + + const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + + const ValueType zero_value {0.}; + // let us handle the dissolution first + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + constexpr Scalar epsilon = std::numeric_limits::epsilon(); + + switch (phaseIdx) { + case FluidSystem::oilPhaseIdx: { + if constexpr (Indices::compositionSwitchIdx >= 0) { + if (both_oil_gas) { + // starting with saturated rs value + ValueType rs = FluidSystem::saturatedDissolutionFactor(fluid_state, phaseIdx, fluid_state.pvtRegionIndex()); + if (fluid_composition[activeCompIdx] > epsilon) { + const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx); + const ValueType max_possible_rs = fluid_composition[gasCompIdx] / fluid_composition[activeCompIdx]; + rs = std::min(rs, max_possible_rs); + } + fluid_state.setRs(rs); + } else { + fluid_state.setRs(zero_value); + } + } + break; + } + case FluidSystem::gasPhaseIdx: { + if constexpr (Indices::compositionSwitchIdx >= 0) { + if (both_oil_gas) { + // starting with saturated rv value + ValueType rv = FluidSystem::saturatedVaporizationFactor(fluid_state, phaseIdx, fluid_state.pvtRegionIndex()); + const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx); + if (fluid_composition[activeCompIdx] > epsilon) { + const ValueType max_possible_rv = fluid_composition[oilCompIdx] / fluid_composition[activeCompIdx]; + rv = std::min(rv, max_possible_rv); + } + fluid_state.setRv(rv); + } else { + fluid_state.setRv(zero_value); + } + } + break; + } + case FluidSystem::waterPhaseIdx: { + // TODO: handle the water phase dissolution with gas later + break; + } + default: { + throw std::logic_error("Unhandled phase index " + std::to_string(phaseIdx)); + } + } + const auto& inv_b = FluidSystem::inverseFormationVolumeFactor(fluid_state, phaseIdx, fluid_state.pvtRegionIndex()); + fluid_state.setInvB(phaseIdx, inv_b); + } + + std::vector saturations (FluidSystem::numPhases, zero_value); + ValueType sum_saturation {0.0}; + // calculate the saturation for all the phases + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) { + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + saturations[phaseIdx] = fluid_composition[activeCompIdx] / fluid_state.invB(phaseIdx); + sum_saturation += saturations[phaseIdx]; + } else { + // remove dissolved gas and vaporized oil + // q_os = q_or * b_o + rv * q_gr * b_g + // q_gs = q_gr * g_g + rs * q_or * b_o + // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) + // d = 1.0 - rs * rv + const ValueType d = 1.0 - fluid_state.Rv() * fluid_state.Rs(); + if (d <= 0.0) { + deferred_logger.debug( + fmt::format("Problematic d value {} obtained for well {}" + " during createFluidState with rs {}" + ", rv {}. Continue as if no dissolution (rs = 0) and" + " vaporization (rv = 0)", + d, this->name(), fluid_state.Rs(), fluid_state.Rv()) ); + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + saturations[phaseIdx] = fluid_composition[activeCompIdx] / fluid_state.invB(phaseIdx); + } else { + const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx); + if (FluidSystem::gasPhaseIdx == phaseIdx) { + saturations[phaseIdx] = (fluid_composition[gasCompIdx] - + fluid_state.Rs() * fluid_composition[oilCompIdx]) / + (d * fluid_state.invB(phaseIdx)); + } else if (FluidSystem::oilPhaseIdx == phaseIdx) { + saturations[phaseIdx] = (fluid_composition[oilCompIdx] - + fluid_state.Rv() * fluid_composition[gasCompIdx]) / + (d * fluid_state.invB(phaseIdx)); + } + } + sum_saturation += saturations[phaseIdx]; + } + } + + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + fluid_state.setSaturation(phaseIdx, saturations[phaseIdx] / sum_saturation); + + typename FluidSystem::template ParameterCache paramCache; + paramCache.setRegionIndex(fluid_state.pvtRegionIndex()); + paramCache.updatePhase(fluid_state, phaseIdx); + fluid_state.setDensity(phaseIdx, FluidSystem::density(fluid_state, paramCache, phaseIdx)); + if constexpr (has_energy) { + fluid_state.setEnthalpy(phaseIdx, FluidSystem::enthalpy(fluid_state, paramCache, phaseIdx)); + } + } + return fluid_state; + } + + // it looks like these functions should go to MultisegmentWellSegments class + template + MultisegmentWell::template SegmentFluidState::EvalWell> + MultisegmentWell::createSegmentFluidState(const int seg, DeferredLogger& deferred_logger) const + { + const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg); + // TODO: 0 should be some other values fixed or we make it a std::optional + const EvalWell seg_temperature = has_energy ? this->primary_variables_.getSegmentTemperature(seg) : 0.; + + // TODO: with the energy equation joins, the num_conservation_quantities will be challenged + std::vector fluid_composition(this->numConservationQuantities(), 0.0); + for (int idx = 0; idx < this->numConservationQuantities(); ++idx) { + fluid_composition[idx] = this->primary_variables_.surfaceVolumeFraction(seg, idx); + } + + return createFluidState(fluid_composition, seg_pressure, seg_temperature, deferred_logger); + } + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + surfaceToReservoirRate(const unsigned phaseIdx, + const SegmentFluidState& segment_fs, + const std::vector& surface_rates, + const int seg, + const std::string_view context, + DeferredLogger& deferred_logger) const + { + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + const EvalWell& invB = segment_fs.invB(phaseIdx); + const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) + && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) { + return surface_rates[activeCompIdx] / invB; + } + + // remove dissolved gas and vaporized oil + const EvalWell rs = segment_fs.Rs(); + const EvalWell rv = segment_fs.Rv(); + const EvalWell d = 1. - rs * rv; + if (d <= 0.0) { + deferred_logger.debug( + fmt::format("Problematic d value {} obtained for well {}, segment {}" + " during {} with rs {}, rv {}. Continue as if no dissolution" + " (rs = 0) and vaporization (rv = 0) for this connection.", + d, this->name(), seg, context, rs, rv)); + return surface_rates[activeCompIdx] / invB; + } + const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx); + if (FluidSystem::gasPhaseIdx == phaseIdx) { + return (surface_rates[gasCompIdx] - rs * surface_rates[oilCompIdx]) / (d * invB); + } + if (FluidSystem::oilPhaseIdx == phaseIdx) { + return (surface_rates[oilCompIdx] - rv * surface_rates[gasCompIdx]) / (d * invB); + } + return EvalWell{0.0}; + } + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + computeSegmentEnergyRate(const int seg, + const int upwind_seg, + const SegmentFluidState& upwind_fs, + const std::string_view context, + DeferredLogger& deferred_logger) const + { + // surface volumetric rates per active phase, scaled by the well efficiency factor + std::vector surface_rates(this->num_conservation_quantities_, 0.0); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + surface_rates[activeCompIdx] = + this->primary_variables_.getSegmentRateUpwinding(seg, + upwind_seg, + activeCompIdx) * + this->well_efficiency_factor_; + } + + EvalWell energy_rate(0.0); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + const EvalWell reservoir_rate = + this->surfaceToReservoirRate(phaseIdx, upwind_fs, surface_rates, + seg, context, deferred_logger); + energy_rate += reservoir_rate * upwind_fs.enthalpy(phaseIdx) * upwind_fs.density(phaseIdx); + } + return energy_rate; + } + + template + void + MultisegmentWell:: + assemblePerforationEnergyEq(const IntensiveQuantities& int_quants, + const std::vector& cq_s, + const int seg, + const int local_perf_index, + DeferredLogger& deferred_logger) + { + // TODO: this part of the implementation remains to be debugged + const auto& fs = int_quants.fluidState(); + // segment fluid state for wellbore properties (used for injecting connections) + const auto& seg_fs = this->segment_fluid_state_[seg]; + + // TODO: should we only use the reservoir-cell properties for production cases? + // The reservoir cell fluid properties may be less accurate for injectors and + // could explain the abnormal results occasionally observed for them. + EvalWell energy_flux(0.0); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + // whether the connection is injecting (fluid flows from wellbore into reservoir) + const bool injecting = cq_s[activeCompIdx] > 0.0; + + EvalWell cq_r_thermal(0.0); + if (injecting) { + // use segment (wellbore) fluid properties for the upwind state + cq_r_thermal = this->surfaceToReservoirRate(phaseIdx, seg_fs, cq_s, + seg, "energy assembly (injecting)", + deferred_logger); + energy_flux += cq_r_thermal * seg_fs.enthalpy(phaseIdx) * seg_fs.density(phaseIdx); + } else { + // producing connection: use reservoir cell fluid properties + const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) + && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + const EvalWell invB = this->extendEval(fs.invB(phaseIdx)); + if (!both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx) { + cq_r_thermal = cq_s[activeCompIdx] / invB; + } else { + const EvalWell rs = this->extendEval(fs.Rs()); + const EvalWell rv = this->extendEval(fs.Rv()); + const EvalWell d = 1. - rs * rv; + if (d <= 0.0) { + deferred_logger.debug( + fmt::format("Problematic d value {} obtained for well {}" + " during energy assembly with rs {}, rv {}." + " Continue as if no dissolution (rs = 0) and" + " vaporization (rv = 0) for this connection.", + d, this->name(), rs, rv)); + cq_r_thermal = cq_s[activeCompIdx] / invB; + } else { + const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx); + if (FluidSystem::gasPhaseIdx == phaseIdx) { + cq_r_thermal = (cq_s[gasCompIdx] - rs * cq_s[oilCompIdx]) / (d * invB); + } else if (FluidSystem::oilPhaseIdx == phaseIdx) { + cq_r_thermal = (cq_s[oilCompIdx] - rv * cq_s[gasCompIdx]) / (d * invB); + } + } + } + energy_flux += cq_r_thermal * this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx)); + } + } + energy_flux *= this->well_efficiency_factor_; + this->connectionRates_[local_perf_index][Indices::contiEnergyEqIdx] = Base::restrictEval(energy_flux); + + MultisegmentWellAssemble(*this). + assemblePerforationEq(seg, local_perf_index, + MSWEval::PrimaryVariables::Temperature, + energy_flux, + this->linSys_); + } + + template + void + MultisegmentWell::computeInitialSegmentEnergy() + { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + segment_initial_energy_[seg] = computeSegmentEnergy(seg); + } + } + + template + void + MultisegmentWell::updateWellHeadCondition(const Simulator& simulator, DeferredLogger& deferred_logger) + { + if (!this->well_ecl_.isInjector()) return; + + std::vector fluid_composition(FluidSystem::numPhases, 0.0); + + // temperature should be the injecting temperature + // pressure should be the BHP + const EvalWell bhp = this->primary_variables_.getSegmentPressure(0); + const EvalWell inj_temperature = this->well_ecl_.inj_temperature(); + + const auto controls = this->well_ecl_.injectionControls(simulator.vanguard().summaryState()); + switch (controls.injector_type) { + case InjectorType::OIL: { + const unsigned oilActiveCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx); + fluid_composition[oilActiveCompIdx] = 1.0; + break; + } + case InjectorType::GAS: { + const unsigned gasActiveCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx); + fluid_composition[gasActiveCompIdx] = 1.0; + break; + } + case InjectorType::WATER: { + const unsigned waterActiveCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx); + fluid_composition[waterActiveCompIdx] = 1.0; + break; + } + default: { + throw std::logic_error("Unsupported injection type " + std::to_string(static_cast(controls.injector_type))); + } + } + + this->wellhead_fluid_state_ = createFluidState(fluid_composition, bhp, inj_temperature, deferred_logger); + } + + + template + template + ValueType + MultisegmentWell::computeSegmentEnergy(const int seg) const + { + auto obtain = [](const auto& val) { + if constexpr (std::is_same_v) { + return getValue(val); + } else { + return val; + } + }; + + ValueType result {0.}; + const auto& segment_fluid_state = this->segment_fluid_state_[seg]; + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + const auto u = obtain(segment_fluid_state.internalEnergy(phaseIdx)); + const auto s = obtain(segment_fluid_state.saturation(phaseIdx)); + const auto rho = obtain(segment_fluid_state.density(phaseIdx)); + const Scalar segment_volume = this->wellEcl().getSegments()[seg].volume(); + result += segment_volume * u * s * rho; + } + return result; + } + + + template + void + MultisegmentWell::updateSegmentFluidState(DeferredLogger& deferred_logger) + { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + segment_fluid_state_[seg] = this->createSegmentFluidState(seg, deferred_logger); + } + } + } // namespace Opm #endif 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) diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 54835a459b2..8f4f525b3ef 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -316,6 +316,8 @@ class WellInterfaceGeneric { // We assume a well to not penetrate more than one pvt region. const int pvtRegionIdx_; + // TODO: with energy conservation joins, we probably should change this to be num_mass_conservation_quantities_ + // or something like that. Energy term calculate somehow differently from the mass term const int num_conservation_quantities_; // number of phases From cb97e9d34b6692c5d403f17373e1f7050c200910 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 29 Apr 2026 21:54:57 +0200 Subject: [PATCH 04/11] adding SPE1CASE2_MSW_THERMAL for regressionTests to test the thermal equation solving for MSW --- regressionTests.cmake | 1 + 1 file changed, 1 insertion(+) diff --git a/regressionTests.cmake b/regressionTests.cmake index 5586df0ac37..5cd25270045 100644 --- a/regressionTests.cmake +++ b/regressionTests.cmake @@ -41,6 +41,7 @@ set(_spe1_tests SPE1CASE2_THERMAL_WATVISC SPE1CASE2_2P SPE1CASE2_TEMP + SPE1CASE2_MSW_THERMAL ) add_multiple_tests( From f3328511d27e3bc0c8bd89c067a2ff357b57515b Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 4 May 2026 17:41:58 +0200 Subject: [PATCH 05/11] computeFluidProperties and getSurfaceVolume uses dynamic temperature for thermal cases in MultisegmentWell. --- opm/simulators/wells/MultisegmentWell.hpp | 3 ++- .../wells/MultisegmentWellSegments.cpp | 18 +++++++++++---- .../wells/MultisegmentWellSegments.hpp | 10 +++++---- .../wells/MultisegmentWell_impl.hpp | 22 +++++++------------ 4 files changed, 30 insertions(+), 23 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index eeee33cbe5f..20699e56fcb 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -59,6 +59,7 @@ 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; @@ -89,7 +90,7 @@ namespace Opm { has_energy, Indices::compositionSwitchIdx >= 0, /*has_watVapor*/ false, - /*has_brine*/ false, + has_brine, /*has_saltPrecip*/ false, /*has_disgas_in_water*/ false, has_solvent, diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 407e0a5dcd7..9fd59898c42 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -21,6 +21,8 @@ #include #include +#include "examples/problems/co2injectionproblem.hh" + #include #include @@ -161,15 +163,21 @@ MultisegmentWellSegments(const int numSegments, template void MultisegmentWellSegments:: -computeFluidProperties(const EvalWell& temperature, - const EvalWell& saltConcentration, +computeFluidProperties(const Scalar firstPerfTemperature, + const Scalar firstPerfSaltConcentration, const PrimaryVariables& primary_variables, DeferredLogger& deferred_logger) { const int num_quantities = well_.numConservationQuantities(); PhaseCalcResult result(static_cast(num_quantities)); + // \Note: we do not have salt concentration supported as primary variable yet. + // \Note: this will change when we implement the brine equation in the multisegment well + const EvalWell saltConcentration = firstPerfSaltConcentration; + for (std::size_t seg = 0; seg < perforations_.size(); ++seg) { + const EvalWell temperature = enable_energy ? primary_variables.getSegmentTemperature(seg) : firstPerfTemperature; + calculatePhaseProperties(result, temperature, saltConcentration, primary_variables, seg, true, deferred_logger); @@ -255,12 +263,14 @@ getPressureDiffSegLocalPerf(const int seg, template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: -getSurfaceVolume(const EvalWell& temperature, - const EvalWell& saltConcentration, +getSurfaceVolume(const Scalar firstPerfTemperature, + const Scalar firstPerfSaltConcentration, const PrimaryVariables& primary_variables, const int seg_idx, DeferredLogger& deferred_logger) const { + const EvalWell saltConcentration {firstPerfSaltConcentration}; + const EvalWell temperature = enable_energy ? primary_variables.getSegmentTemperature(seg_idx) : firstPerfTemperature; PhaseCalcResult result(static_cast(well_.numConservationQuantities())); calculatePhaseProperties(result, temperature, saltConcentration, primary_variables, seg_idx, false, deferred_logger); diff --git a/opm/simulators/wells/MultisegmentWellSegments.hpp b/opm/simulators/wells/MultisegmentWellSegments.hpp index 7aa6450705c..6f2b199db65 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.hpp +++ b/opm/simulators/wells/MultisegmentWellSegments.hpp @@ -48,13 +48,15 @@ class MultisegmentWellSegments using EvalWell = typename PrimaryVariables::EvalWell; using IndexTraits = typename FluidSystem::IndexTraitsType; + static constexpr bool enable_energy = Indices::enableFullyImplicitThermal; + public: MultisegmentWellSegments(const int numSegments, const ParallelWellInfo& parallel_well_info, WellInterfaceGeneric& well); - void computeFluidProperties(const EvalWell& temperature, - const EvalWell& saltConcentration, + void computeFluidProperties(const Scalar firstPerfTemperature, + const Scalar firstPerfSaltConcentration, const PrimaryVariables& primary_variables, DeferredLogger& deferred_logger); @@ -68,8 +70,8 @@ class MultisegmentWellSegments Scalar getPressureDiffSegLocalPerf(const int seg, const int local_perf_index) const; - EvalWell getSurfaceVolume(const EvalWell& temperature, - const EvalWell& saltConcentration, + EvalWell getSurfaceVolume(const Scalar firstPerfTemperature, + const Scalar firstPerfSaltConcentration, const PrimaryVariables& primary_variables, const int seg_idx, DeferredLogger& deferred_logger) const; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index ff794ec1adf..7617c1faa85 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1149,9 +1149,6 @@ namespace Opm // get the temperature for later use. It is only useful when we are not handling // thermal related simulation // basically, it is a single value for all the segments - - EvalWell temperature; - EvalWell saltConcentration; // not sure how to handle the pvt region related to segment // for the current approach, we use the pvt region of the first perforated cell // although there are some text indicating using the pvt region of the lowest @@ -1159,11 +1156,11 @@ namespace Opm // TODO: later to investigate how to handle the pvt region auto info = this->getFirstPerforationFluidStateInfo(simulator); - temperature.setValue(std::get<0>(info)); - saltConcentration.setValue(std::get<1>(info)); + const Scalar firstPerfTemperature = std::get<0>(info); + const Scalar firstPerfSaltConcentration = std::get<1>(info); - this->segments_.computeFluidProperties(temperature, - saltConcentration, + this->segments_.computeFluidProperties(firstPerfTemperature, + firstPerfSaltConcentration, this->primary_variables_, deferred_logger); } @@ -2153,15 +2150,12 @@ namespace Opm const int seg_idx, DeferredLogger& deferred_logger) const { - EvalWell temperature; - EvalWell saltConcentration; - auto info = this->getFirstPerforationFluidStateInfo(simulator); - temperature.setValue(std::get<0>(info)); - saltConcentration.setValue(std::get<1>(info)); + const Scalar firstPerfTemperature = std::get<0>(info); + const Scalar firstPerfSaltConcentration = std::get<1>(info); - return this->segments_.getSurfaceVolume(temperature, - saltConcentration, + return this->segments_.getSurfaceVolume(firstPerfTemperature, + firstPerfSaltConcentration, this->primary_variables_, seg_idx, deferred_logger); From ea0f818821e57459f93f1075143154d1992f4826 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 6 May 2026 15:17:13 +0200 Subject: [PATCH 06/11] using perforation temperature for non-thermal createSegmentFluidState() instead of using 0, which is not physical --- opm/simulators/wells/MultisegmentWell.hpp | 10 ++--- .../wells/MultisegmentWell_impl.hpp | 37 +++++++++++-------- 2 files changed, 27 insertions(+), 20 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 20699e56fcb..2332d3d5b82 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -217,7 +217,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); @@ -322,8 +322,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 @@ -376,7 +376,7 @@ namespace Opm { DeferredLogger& deferred_logger) const; SegmentFluidState - createSegmentFluidState(int seg, DeferredLogger& deferred_logger) const; + createSegmentFluidState(int seg, const FSInfo& info, DeferredLogger& deferred_logger) const; void computeInitialSegmentEnergy(); @@ -389,7 +389,7 @@ namespace Opm { void updateWellHeadCondition(const Simulator& simulator, DeferredLogger& deferred_logger); - void updateSegmentFluidState(DeferredLogger& deferred_logger); + void updateSegmentFluidState(const FSInfo& info, DeferredLogger& deferred_logger); template ValueType computeSegmentEnergy(int seg) const; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 7617c1faa85..6c9e161b517 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -696,12 +696,11 @@ namespace Opm template void MultisegmentWell:: - computeInitialSegmentFluids(const Simulator& simulator, + computeInitialSegmentFluids(const FSInfo& info, DeferredLogger& deferred_logger) { for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - // TODO: trying to reduce the times for the surfaceVolumeFraction calculation - const Scalar surface_volume = getSegmentSurfaceVolume(simulator, seg, deferred_logger).value(); + const Scalar surface_volume = getSegmentSurfaceVolume(seg, info, deferred_logger).value(); for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) { segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value(); } @@ -766,10 +765,12 @@ namespace Opm updatePrimaryVariables(groupStateHelper); computePerfCellPressDiffs(simulator); + const auto info = this->getFirstPerforationFluidStateInfo(simulator); + // after updating the primary variables, we need to update the segment fluid state - updateSegmentFluidState(deferred_logger); + updateSegmentFluidState(info, deferred_logger); - computeInitialSegmentFluids(simulator, deferred_logger); + computeInitialSegmentFluids(info, deferred_logger); if constexpr (has_energy) { computeInitialSegmentEnergy(); } @@ -1788,7 +1789,8 @@ namespace Opm try{ const BVectorWell dx_well = this->linSys_.solve(); updateWellState(simulator, dx_well, groupStateHelper, well_state, relaxation_factor); - updateSegmentFluidState(deferred_logger); + const FSInfo info = this->getFirstPerforationFluidStateInfo(simulator); + updateSegmentFluidState(info, deferred_logger); } catch(const NumericalProblem& exp) { // Add information about the well and log to deferred logger @@ -1943,10 +1945,12 @@ namespace Opm this->linSys_.sumDistributed(this->parallel_well_info_.communication()); } + const FSInfo info = this->getFirstPerforationFluidStateInfo(simulator); + for (int seg = 0; seg < nseg; ++seg) { // calculating the accumulation term { - const EvalWell segment_surface_volume = getSegmentSurfaceVolume(simulator, seg, deferred_logger); + const EvalWell segment_surface_volume = getSegmentSurfaceVolume(seg, info, deferred_logger); // Add a regularization_factor to increase the accumulation term // This will make the system less stiff and help convergence for @@ -2146,11 +2150,10 @@ namespace Opm template typename MultisegmentWell::EvalWell MultisegmentWell:: - getSegmentSurfaceVolume(const Simulator& simulator, - const int seg_idx, + getSegmentSurfaceVolume(const int seg_idx, + const FSInfo& info, DeferredLogger& deferred_logger) const { - auto info = this->getFirstPerforationFluidStateInfo(simulator); const Scalar firstPerfTemperature = std::get<0>(info); const Scalar firstPerfSaltConcentration = std::get<1>(info); @@ -2582,11 +2585,14 @@ namespace Opm // it looks like these functions should go to MultisegmentWellSegments class template MultisegmentWell::template SegmentFluidState::EvalWell> - MultisegmentWell::createSegmentFluidState(const int seg, DeferredLogger& deferred_logger) const + MultisegmentWell::createSegmentFluidState(const int seg, const FSInfo& info, + DeferredLogger& deferred_logger) const { const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg); - // TODO: 0 should be some other values fixed or we make it a std::optional - const EvalWell seg_temperature = has_energy ? this->primary_variables_.getSegmentTemperature(seg) : 0.; + const Scalar firstPerfTemperature = std::get<0>(info); + // TODO: we can extend the salt concentration to the createFluidState + // const Scalar firstPerfSaltConcentration = std::get<1>(info); + const EvalWell seg_temperature = has_energy ? this->primary_variables_.getSegmentTemperature(seg) : firstPerfTemperature; // TODO: with the energy equation joins, the num_conservation_quantities will be challenged std::vector fluid_composition(this->numConservationQuantities(), 0.0); @@ -2829,10 +2835,11 @@ namespace Opm template void - MultisegmentWell::updateSegmentFluidState(DeferredLogger& deferred_logger) + MultisegmentWell::updateSegmentFluidState(const FSInfo& info, + DeferredLogger& deferred_logger) { for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - segment_fluid_state_[seg] = this->createSegmentFluidState(seg, deferred_logger); + segment_fluid_state_[seg] = this->createSegmentFluidState(seg, info, deferred_logger); } } From a611b604327e8f162dc2b1955e498969e8761df5 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 7 May 2026 10:36:01 +0200 Subject: [PATCH 07/11] using the FluidState definition from the WellInterface --- opm/simulators/wells/MultisegmentWell.hpp | 19 +++------------ opm/simulators/wells/WellInterface.hpp | 28 +++++++++++++---------- 2 files changed, 19 insertions(+), 28 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 2332d3d5b82..f5a6549d2db 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -79,22 +79,9 @@ namespace Opm { // 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 - // templated on Value type because we want to use it for surface condition calculation injection - // while we shifted to use the bottom-hole/well-head condition instead, so there is only EvalWell type - // used here for now. - template - using SegmentFluidState = BlackOilFluidState= 0, - /*has_watVapor*/ false, - has_brine, - /*has_saltPrecip*/ false, - /*has_disgas_in_water*/ false, - has_solvent, - Indices::numPhases>; + // calculation in the wellbore. + template + using SegmentFluidState = Base::template BlackOilFluidStateType; MultisegmentWell(const Well& well, const ParallelWellInfo& pw_info, diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index bd5e5ba459e..c5df4a98195 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -126,18 +126,22 @@ class WellInterface : public WellInterfaceIndices(); static constexpr bool has_micp = Indices::enableMICP; - // For the conversion between the surface volume rate and reservoir voidage rate - using FluidState = BlackOilFluidState= 0, - has_watVapor, - has_brine, - has_saltPrecip, - has_disgas_in_water, - has_solvent, - Indices::numPhases >; + template + using BlackOilFluidStateType = BlackOilFluidState= 0, + has_watVapor, + has_brine, + has_saltPrecip, + has_disgas_in_water, + has_solvent, + Indices::numPhases>; + + // fluid state for the reservoir fluid + using FluidState = BlackOilFluidStateType; + /// Constructor WellInterface(const Well& well, const ParallelWellInfo& pw_info, From 29aac30bd1e559c0beb163e3c741c08722760da5 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 7 May 2026 14:21:52 +0200 Subject: [PATCH 08/11] incoporating saltConcentration to SegmentFluidState and using the calculation in createSegmentFluidState to calculate the surface volume of the segments to reduce duplicated calculation. --- opm/simulators/wells/MultisegmentWell.hpp | 8 +++-- .../wells/MultisegmentWellSegments.cpp | 19 ++++++++++ .../wells/MultisegmentWellSegments.hpp | 11 ++++++ .../wells/MultisegmentWell_impl.hpp | 36 ++++++++++++++----- 4 files changed, 63 insertions(+), 11 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index f5a6549d2db..63aff6619a2 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -75,6 +75,8 @@ namespace Opm { using typename MSWEval::BVectorWell; using MSWEval::SPres; using typename Base::PressureMatrix; + // it contains the temperature and salt concentration of the first perforation + // TODO: kind of temporary measure, will evolve or be replaced based on future development and understanding using FSInfo = std::tuple; // a fluid state to calculate the properties inside the wellbore for each segment @@ -360,10 +362,12 @@ namespace Opm { createFluidState(const std::vector& fluid_composition, const ValueType& pressure, const ValueType& temperature, + const ValueType& saltConcentration, + ValueType& volume_ratio, DeferredLogger& deferred_logger) const; SegmentFluidState - createSegmentFluidState(int seg, const FSInfo& info, DeferredLogger& deferred_logger) const; + createSegmentFluidState(int seg, const FSInfo& info, DeferredLogger& deferred_logger); void computeInitialSegmentEnergy(); @@ -374,7 +378,7 @@ namespace Opm { const int local_perf_index, DeferredLogger& deferred_logger); - void updateWellHeadCondition(const Simulator& simulator, DeferredLogger& deferred_logger); + void updateWellHeadCondition(const Simulator& simulator, const FSInfo& info, DeferredLogger& deferred_logger); void updateSegmentFluidState(const FSInfo& info, DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 9fd59898c42..bf0c3d9d433 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -101,6 +101,7 @@ MultisegmentWellSegments(const int numSegments, , phase_densities_(numSegments, std::vector(well.numConservationQuantities(), 0.0)) // number of phase here? , phase_fractions_(numSegments, std::vector(well.numConservationQuantities(), 0.0)) // number of phase here? , phase_viscosities_(numSegments, std::vector(well.numConservationQuantities(), 0.0)) // number of phase here? + , volume_ratios_(numSegments, 0.0) , well_(well) { // since we decide to use the WellSegments from the well parser. we can reuse a lot from it. @@ -189,6 +190,8 @@ computeFluidProperties(const Scalar firstPerfTemperature, const auto& b = result.b; const auto& volrat = result.vol_ratio; + volume_ratios_[seg] = volrat; + viscosities_[seg] = 0.; // calculate the average viscosity for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { @@ -283,6 +286,22 @@ getSurfaceVolume(const Scalar firstPerfTemperature, return volume / vol_ratio; } +template +typename MultisegmentWellSegments::EvalWell +MultisegmentWellSegments:: +getSurfaceVolume(const int seg_idx, + DeferredLogger& deferred_logger) const +{ + const Scalar volume = well_.wellEcl().getSegments()[seg_idx].volume(); + const auto& vol_ratio = volume_ratios_[seg_idx]; + if (vol_ratio.value() == 0.) { + deferred_logger.warning(fmt::format("Volume ratio for segment {} for well {} is zero, " + "returning segment volume as surface volume.", seg_idx, this->well_.name())); + return volume; + } + return volume / vol_ratio; +} + template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: diff --git a/opm/simulators/wells/MultisegmentWellSegments.hpp b/opm/simulators/wells/MultisegmentWellSegments.hpp index 6f2b199db65..02aacbaf139 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.hpp +++ b/opm/simulators/wells/MultisegmentWellSegments.hpp @@ -76,6 +76,9 @@ class MultisegmentWellSegments const int seg_idx, DeferredLogger& deferred_logger) const; + EvalWell getSurfaceVolume(const int seg_idx, + DeferredLogger& deferred_logger) const; + EvalWell getFrictionPressureLoss(const int seg, const bool extra_reverse_flow_derivatives = false) const; @@ -128,6 +131,11 @@ class MultisegmentWellSegments return densities_[seg]; } + EvalWell& volumeRatio(const int seg) + { + return volume_ratios_[seg]; + } + Scalar local_perforation_depth_diff(const int local_perf_index) const { return local_perforation_depth_diffs_[local_perf_index]; @@ -179,6 +187,9 @@ class MultisegmentWellSegments std::vector> phase_fractions_; std::vector> phase_viscosities_; + // volume ratio between wellbore condition and surface condition for each segment + std::vector volume_ratios_; + WellInterfaceGeneric& well_; void copyPhaseDensities(const unsigned phaseIdx, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 6c9e161b517..29c344e73fb 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -700,7 +700,7 @@ namespace Opm DeferredLogger& deferred_logger) { for (int seg = 0; seg < this->numberOfSegments(); ++seg) { - const Scalar surface_volume = getSegmentSurfaceVolume(seg, info, deferred_logger).value(); + const Scalar surface_volume = this->segments_.getSurfaceVolume(seg, deferred_logger).value(); for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) { segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value(); } @@ -1950,7 +1950,7 @@ namespace Opm for (int seg = 0; seg < nseg; ++seg) { // calculating the accumulation term { - const EvalWell segment_surface_volume = getSegmentSurfaceVolume(seg, info, deferred_logger); + const EvalWell segment_surface_volume = this->segments_.getSurfaceVolume(seg, deferred_logger); // Add a regularization_factor to increase the accumulation term // This will make the system less stiff and help convergence for @@ -1986,7 +1986,7 @@ namespace Opm if constexpr (has_energy) { const bool top_injecting_segment = (seg == 0) && this->isInjector(); if (top_injecting_segment) { - this->updateWellHeadCondition(simulator, deferred_logger); + this->updateWellHeadCondition(simulator, info, deferred_logger); } // Energy carried by fluid flowing out of this segment toward its outlet. @@ -2447,12 +2447,17 @@ namespace Opm createFluidState(const std::vector& fluid_composition, const ValueType& pressure, const ValueType& temperature, + const ValueType& saltConcentration, + ValueType& volume_ratio, DeferredLogger& deferred_logger) const { SegmentFluidState fluid_state; if constexpr (has_energy) { fluid_state.setTemperature(temperature); } + if constexpr (has_brine) { + fluid_state.setSaltConcentration(saltConcentration); + } for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { continue; @@ -2483,7 +2488,7 @@ namespace Opm if (fluid_composition[activeCompIdx] > epsilon) { const unsigned gasCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx); const ValueType max_possible_rs = fluid_composition[gasCompIdx] / fluid_composition[activeCompIdx]; - rs = std::min(rs, max_possible_rs); + rs = std::clamp(rs, zero_value, max_possible_rs); } fluid_state.setRs(rs); } else { @@ -2500,7 +2505,7 @@ namespace Opm const unsigned oilCompIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx); if (fluid_composition[activeCompIdx] > epsilon) { const ValueType max_possible_rv = fluid_composition[oilCompIdx] / fluid_composition[activeCompIdx]; - rv = std::min(rv, max_possible_rv); + rv = std::clamp(rv, zero_value, max_possible_rv); } fluid_state.setRv(rv); } else { @@ -2564,6 +2569,7 @@ namespace Opm sum_saturation += saturations[phaseIdx]; } } + volume_ratio = sum_saturation; for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { @@ -2586,13 +2592,17 @@ namespace Opm template MultisegmentWell::template SegmentFluidState::EvalWell> MultisegmentWell::createSegmentFluidState(const int seg, const FSInfo& info, - DeferredLogger& deferred_logger) const + DeferredLogger& deferred_logger) { const EvalWell seg_pressure = this->primary_variables_.getSegmentPressure(seg); const Scalar firstPerfTemperature = std::get<0>(info); + const Scalar firstPerfSaltConcentration = std::get<1>(info); // TODO: we can extend the salt concentration to the createFluidState // const Scalar firstPerfSaltConcentration = std::get<1>(info); const EvalWell seg_temperature = has_energy ? this->primary_variables_.getSegmentTemperature(seg) : firstPerfTemperature; + // \Note: salt concentration is not yet a primary variable for the multisegment well; + // we use the value taken at the first perforation, matching what computeSegmentFluidProperties() did before. + const EvalWell seg_salt_concentration {firstPerfSaltConcentration}; // TODO: with the energy equation joins, the num_conservation_quantities will be challenged std::vector fluid_composition(this->numConservationQuantities(), 0.0); @@ -2600,7 +2610,11 @@ namespace Opm fluid_composition[idx] = this->primary_variables_.surfaceVolumeFraction(seg, idx); } - return createFluidState(fluid_composition, seg_pressure, seg_temperature, deferred_logger); + // we also update the volume ratio for this segment to calculate the surface volume later + auto& volume_ratio = this->segments_.volumeRatio(seg); + + return createFluidState(fluid_composition, seg_pressure, seg_temperature, + seg_salt_concentration, volume_ratio, deferred_logger); } template @@ -2767,7 +2781,9 @@ namespace Opm template void - MultisegmentWell::updateWellHeadCondition(const Simulator& simulator, DeferredLogger& deferred_logger) + MultisegmentWell::updateWellHeadCondition(const Simulator& simulator, + const FSInfo& info, + DeferredLogger& deferred_logger) { if (!this->well_ecl_.isInjector()) return; @@ -2800,7 +2816,9 @@ namespace Opm } } - this->wellhead_fluid_state_ = createFluidState(fluid_composition, bhp, inj_temperature, deferred_logger); + EvalWell volume_ratio {0.0}; // it is not needed, while for interface reason. + const EvalWell saltConcentration = std::get<1>(info); // saltConcentration + this->wellhead_fluid_state_ = createFluidState(fluid_composition, bhp, inj_temperature, saltConcentration, volume_ratio, deferred_logger); } From e483210dc8d0d67a77c4d6ca12ed7b2ad435f014 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sun, 10 May 2026 00:05:00 +0200 Subject: [PATCH 09/11] trying to use segment fluid state to calculate fluid properties remains to be cleaned up. --- .../wells/MultisegmentWellSegments.cpp | 29 ++++++++- .../wells/MultisegmentWellSegments.hpp | 17 ++++++ .../wells/MultisegmentWell_impl.hpp | 59 ++++++++++++++++--- 3 files changed, 96 insertions(+), 9 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index bf0c3d9d433..361d8a98e32 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -47,6 +47,7 @@ #include #include +#include #include #include #include @@ -58,7 +59,7 @@ namespace { template -std::vector surfaceDensities(const int pvt_region, const std::size_t num_quantities) +std::vector getSurfaceDensities(const int pvt_region, const std::size_t num_quantities) { using Scalar = typename FluidSystem::Scalar; auto surf_dens = std::vector(num_quantities); @@ -93,7 +94,7 @@ MultisegmentWellSegments(const int numSegments, // local information. This is an exception and intentionally, since here, we only need the local entries. , inlets_(well.wellEcl().getSegments().size()) , depth_diffs_(numSegments, 0.0) - , surface_densities_(surfaceDensities(well.pvtRegionIdx(), well.numConservationQuantities())) + , surface_densities_(getSurfaceDensities(well.pvtRegionIdx(), well.numConservationQuantities())) , densities_(numSegments, 0.0) , mass_rates_(numSegments, 0.0) , viscosities_(numSegments, 0.0) @@ -220,6 +221,30 @@ computeFluidProperties(const Scalar firstPerfTemperature, } } +template +void MultisegmentWellSegments:: +updateFluidProperties(const std::vector>& phase_densities, + const std::vector>& phase_viscosities, + const std::vector>& phase_fractions, + const std::vector& viscosities, + const std::vector& densities, + const std::vector& mass_rates) +{ + assert(phase_densities.size() == phase_densities_.size()); + assert(phase_viscosities.size() == phase_viscosities_.size()); + assert(phase_fractions.size() == phase_fractions_.size()); + assert(viscosities.size() == viscosities_.size()); + assert(densities.size() == densities_.size()); + assert(mass_rates.size() == mass_rates_.size()); + + phase_densities_ = phase_densities; + phase_viscosities_ = phase_viscosities; + phase_fractions_ = phase_fractions; + viscosities_ = viscosities; + densities_ = densities; + mass_rates_ = mass_rates; +} + template void MultisegmentWellSegments:: updateUpwindingSegments(const PrimaryVariables& primary_variables) diff --git a/opm/simulators/wells/MultisegmentWellSegments.hpp b/opm/simulators/wells/MultisegmentWellSegments.hpp index 02aacbaf139..ef4bba94890 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.hpp +++ b/opm/simulators/wells/MultisegmentWellSegments.hpp @@ -60,6 +60,13 @@ class MultisegmentWellSegments const PrimaryVariables& primary_variables, DeferredLogger& deferred_logger); + void updateFluidProperties(const std::vector>& phase_densities, + const std::vector>& phase_viscosities, + const std::vector>& phase_fractions, + const std::vector& viscosities, + const std::vector& densities, + const std::vector& mass_rates); + //! \brief Update upwinding segments. void updateUpwindingSegments(const PrimaryVariables& primary_variables); @@ -131,6 +138,11 @@ class MultisegmentWellSegments return densities_[seg]; } + const EvalWell& volumeRatio(const int seg) const + { + return volume_ratios_[seg]; + } + EvalWell& volumeRatio(const int seg) { return volume_ratios_[seg]; @@ -141,6 +153,11 @@ class MultisegmentWellSegments return local_perforation_depth_diffs_[local_perf_index]; } + const std::vector& surfaceDensities() const + { + return surface_densities_; + } + void copyPhaseDensities(SegmentState& segSol) const; private: diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 29c344e73fb..ef31f49c7c0 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1144,6 +1144,9 @@ namespace Opm MultisegmentWell:: computeSegmentFluidProperties(const Simulator& simulator, DeferredLogger& deferred_logger) { + static_cast(simulator); + static_cast(deferred_logger); + // TODO: the concept of phases and components are rather confusing in this function. // needs to be addressed sooner or later. @@ -1155,15 +1158,57 @@ namespace Opm // although there are some text indicating using the pvt region of the lowest // perforated cell // TODO: later to investigate how to handle the pvt region + // + // auto info = this->getFirstPerforationFluidStateInfo(simulator); + // const Scalar firstPerfTemperature = std::get<0>(info); + // const Scalar firstPerfSaltConcentration = std::get<1>(info); + // + // this->segments_.computeFluidProperties(firstPerfTemperature, + // firstPerfSaltConcentration, + // this->primary_variables_, + // deferred_logger); + const int nseg = this->numberOfSegments(); + const int nquantities = this->numConservationQuantities(); + const auto nseg_size = static_cast(nseg); + const auto nquantities_size = static_cast(nquantities); + std::vector> phase_densities(nseg_size, std::vector(nquantities_size, 0.0)); + std::vector> phase_viscosities(nseg_size, std::vector(nquantities_size, 0.0)); + std::vector> phase_fractions(nseg_size, std::vector(nquantities_size, 0.0)); + std::vector viscosities(nseg_size, 0.0); + std::vector densities(nseg_size, 0.0); + std::vector mass_rates(nseg_size, 0.0); + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { + const auto& fs = this->segment_fluid_state_[seg]; + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + const int compIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); - auto info = this->getFirstPerforationFluidStateInfo(simulator); - const Scalar firstPerfTemperature = std::get<0>(info); - const Scalar firstPerfSaltConcentration = std::get<1>(info); + phase_fractions[seg][compIdx] = fs.saturation(phaseIdx); + phase_viscosities[seg][compIdx] = fs.viscosity(phaseIdx); + viscosities[seg] += phase_fractions[seg][compIdx] * phase_viscosities[seg][compIdx]; + + phase_densities[seg][compIdx] = fs.density(phaseIdx); + densities[seg] += phase_fractions[seg][compIdx] * phase_densities[seg][compIdx]; + } + const auto& surface_densities = this->segments_.surfaceDensities(); + + for (int compIdx = 0; compIdx < nquantities; ++compIdx) { + const int upwind_seg = this->segments_.upwinding_segment(seg); + const EvalWell rate = this->primary_variables_.getSegmentRateUpwinding(seg, + upwind_seg, + compIdx); + mass_rates[seg] += rate * surface_densities[compIdx]; + } + } - this->segments_.computeFluidProperties(firstPerfTemperature, - firstPerfSaltConcentration, - this->primary_variables_, - deferred_logger); + this->segments_.updateFluidProperties(phase_densities, + phase_viscosities, + phase_fractions, + viscosities, + densities, + mass_rates); } template From 311b76b71aa9825d24d87ae4adf4e7fd0b844b9a Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Sun, 10 May 2026 11:31:15 +0200 Subject: [PATCH 10/11] WIP in cleaning up --- .../wells/MultisegmentWellSegments.cpp | 289 ++++-------------- .../wells/MultisegmentWellSegments.hpp | 39 --- .../wells/MultisegmentWell_impl.hpp | 1 + 3 files changed, 58 insertions(+), 271 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 361d8a98e32..ad4da37cdd2 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -163,63 +163,63 @@ MultisegmentWellSegments(const int numSegments, } } -template -void MultisegmentWellSegments:: -computeFluidProperties(const Scalar firstPerfTemperature, - const Scalar firstPerfSaltConcentration, - const PrimaryVariables& primary_variables, - DeferredLogger& deferred_logger) -{ - const int num_quantities = well_.numConservationQuantities(); - PhaseCalcResult result(static_cast(num_quantities)); - - // \Note: we do not have salt concentration supported as primary variable yet. - // \Note: this will change when we implement the brine equation in the multisegment well - const EvalWell saltConcentration = firstPerfSaltConcentration; - - for (std::size_t seg = 0; seg < perforations_.size(); ++seg) { - const EvalWell temperature = enable_energy ? primary_variables.getSegmentTemperature(seg) : firstPerfTemperature; - - calculatePhaseProperties(result, temperature, saltConcentration, - primary_variables, seg, true, deferred_logger); - - phase_densities_[seg] = result.phase_densities; - phase_viscosities_[seg] = result.phase_viscosities; - - const auto& mix = result.mix; - const auto& mix_s = result.mix_s; - const auto& b = result.b; - const auto& volrat = result.vol_ratio; - - volume_ratios_[seg] = volrat; - - viscosities_[seg] = 0.; - // calculate the average viscosity - for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { - const EvalWell fraction = mix[comp_idx] / b[comp_idx] / volrat; - // TODO: a little more work needs to be done to handle the negative fractions here - phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0; - viscosities_[seg] += phase_viscosities_[seg][comp_idx] * phase_fractions_[seg][comp_idx]; - } - - EvalWell density(0.0); - for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { - density += surface_densities_[comp_idx] * mix_s[comp_idx]; - } - densities_[seg] = density / volrat; - - // TODO: the enthalpy flux rate should be calculated in the similar manner. - // calculate the mass rates - mass_rates_[seg] = 0.; - for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) { - const int upwind_seg = upwinding_segments_[seg]; - const EvalWell rate = primary_variables.getSegmentRateUpwinding(seg, - upwind_seg, - comp_idx); - mass_rates_[seg] += rate * surface_densities_[comp_idx]; - } - } -} +// template +// void MultisegmentWellSegments:: +// computeFluidProperties(const Scalar firstPerfTemperature, +// const Scalar firstPerfSaltConcentration, +// const PrimaryVariables& primary_variables, +// DeferredLogger& deferred_logger) +// { +// const int num_quantities = well_.numConservationQuantities(); +// PhaseCalcResult result(static_cast(num_quantities)); +// +// // \Note: we do not have salt concentration supported as primary variable yet. +// // \Note: this will change when we implement the brine equation in the multisegment well +// const EvalWell saltConcentration = firstPerfSaltConcentration; +// +// for (std::size_t seg = 0; seg < perforations_.size(); ++seg) { +// const EvalWell temperature = enable_energy ? primary_variables.getSegmentTemperature(seg) : firstPerfTemperature; +// +// calculatePhaseProperties(result, temperature, saltConcentration, +// primary_variables, seg, true, deferred_logger); +// +// phase_densities_[seg] = result.phase_densities; +// phase_viscosities_[seg] = result.phase_viscosities; +// +// const auto& mix = result.mix; +// const auto& mix_s = result.mix_s; +// const auto& b = result.b; +// const auto& volrat = result.vol_ratio; +// +// volume_ratios_[seg] = volrat; +// +// viscosities_[seg] = 0.; +// // calculate the average viscosity +// for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { +// const EvalWell fraction = mix[comp_idx] / b[comp_idx] / volrat; +// // TODO: a little more work needs to be done to handle the negative fractions here +// phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0; +// viscosities_[seg] += phase_viscosities_[seg][comp_idx] * phase_fractions_[seg][comp_idx]; +// } +// +// EvalWell density(0.0); +// for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { +// density += surface_densities_[comp_idx] * mix_s[comp_idx]; +// } +// densities_[seg] = density / volrat; +// +// // TODO: the enthalpy flux rate should be calculated in the similar manner. +// // calculate the mass rates +// mass_rates_[seg] = 0.; +// for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) { +// const int upwind_seg = upwinding_segments_[seg]; +// const EvalWell rate = primary_variables.getSegmentRateUpwinding(seg, +// upwind_seg, +// comp_idx); +// mass_rates_[seg] += rate * surface_densities_[comp_idx]; +// } +// } +// } template void MultisegmentWellSegments:: @@ -288,29 +288,6 @@ getPressureDiffSegLocalPerf(const int seg, return well_.gravity() * densities_[seg].value() * local_perforation_depth_diffs_[local_perf_index]; } -template -typename MultisegmentWellSegments::EvalWell -MultisegmentWellSegments:: -getSurfaceVolume(const Scalar firstPerfTemperature, - const Scalar firstPerfSaltConcentration, - const PrimaryVariables& primary_variables, - const int seg_idx, - DeferredLogger& deferred_logger) const -{ - const EvalWell saltConcentration {firstPerfSaltConcentration}; - const EvalWell temperature = enable_energy ? primary_variables.getSegmentTemperature(seg_idx) : firstPerfTemperature; - PhaseCalcResult result(static_cast(well_.numConservationQuantities())); - calculatePhaseProperties(result, temperature, saltConcentration, - primary_variables, seg_idx, false, deferred_logger); - - const EvalWell& vol_ratio = result.vol_ratio; - - // We increase the segment volume with a factor 10 to stabilize the system. - const Scalar volume = well_.wellEcl().getSegments()[seg_idx].volume(); - - return volume / vol_ratio; -} - template typename MultisegmentWellSegments::EvalWell MultisegmentWellSegments:: @@ -821,158 +798,6 @@ mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const return mixDens; } -template -void -MultisegmentWellSegments:: -calculatePhaseProperties(PhaseCalcResult& result, - const EvalWell& temperature, - const EvalWell& saltConcentration, - const PrimaryVariables& primary_variables, - const int seg, - const bool update_visc_and_den, - DeferredLogger& deferred_logger) const -{ - result.clear(); - - auto& b = result.b; - auto& mix_s = result.mix_s; - auto& mix = result.mix; - auto& vol_ratio = result.vol_ratio; - auto& phase_viscosities = result.phase_viscosities; - auto& phase_densities = result.phase_densities; - - // we might use reference here - const EvalWell seg_pressure = primary_variables.getSegmentPressure(seg); - - const int num_quantities = well_.numConservationQuantities(); - for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { - mix_s[comp_idx] = primary_variables.surfaceVolumeFraction(seg, comp_idx); - } - - const bool waterActive = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); - const bool gasActive = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - const bool oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - - const int waterActiveCompIdx = waterActive ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::waterCompIdx) : -1; - const int gasActiveCompIdx = gasActive ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::gasCompIdx) : -1; - const int oilActiveCompIdx = oilActive ? FluidSystem::canonicalToActiveCompIdx(FluidSystem::oilCompIdx) : -1; - - const int pvt_region_index = well_.pvtRegionIdx(); - - // water phase - if (waterActive) { - // rsw is only for interface usage - const EvalWell rsw{0.}; - b[waterActiveCompIdx] = FluidSystem::waterPvt().inverseFormationVolumeFactor( - pvt_region_index, temperature, seg_pressure, rsw, saltConcentration); - if (update_visc_and_den) { - // TODO: should not we use phaseIndex here? - phase_viscosities[waterActiveCompIdx] = FluidSystem::waterPvt().viscosity( - pvt_region_index, temperature, seg_pressure, rsw, saltConcentration); - phase_densities[waterActiveCompIdx] = b[waterActiveCompIdx] * surface_densities_[waterActiveCompIdx]; - } - } - - EvalWell rv {0.}; - // gas phase - if (gasActive) { - // rvw is only for interface usage - const EvalWell rvw{0.}; - const bool oil_exist = oilActive && mix_s[oilActiveCompIdx] > 0.0; - if (oil_exist) { - const EvalWell rvmax = max( - FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure), - 0.); - if (mix_s[gasActiveCompIdx] > 0.0) { - rv = std::clamp(mix_s[oilActiveCompIdx] / mix_s[gasActiveCompIdx], EvalWell{0.}, rvmax); - } - b[gasActiveCompIdx] = FluidSystem::gasPvt().inverseFormationVolumeFactor( - pvt_region_index, temperature, seg_pressure, rv, rvw); - if (update_visc_and_den) { - phase_viscosities[gasActiveCompIdx] = FluidSystem::gasPvt().viscosity( - pvt_region_index, temperature, seg_pressure, rv, rvw); - phase_densities[gasActiveCompIdx] = b[gasActiveCompIdx] * surface_densities_[gasActiveCompIdx] - + rv * b[gasActiveCompIdx] * surface_densities_[oilActiveCompIdx]; - } - } else { // no oil here - b[gasActiveCompIdx] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor( - pvt_region_index, temperature, seg_pressure); - if (update_visc_and_den) { - phase_viscosities[gasActiveCompIdx] = FluidSystem::gasPvt().saturatedViscosity( - pvt_region_index, temperature, seg_pressure); - phase_densities[gasActiveCompIdx] = b[gasActiveCompIdx] * surface_densities_[gasActiveCompIdx]; - } - } - } - - EvalWell rs {0.}; - // oil phase - if (oilActive) { - const bool gas_exist = gasActive && mix_s[gasActiveCompIdx] > 0.0; - if (gas_exist) { - const EvalWell rsmax = max( - FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure), - 0.); - if (mix_s[oilActiveCompIdx] > 0.0) { - rs = std::clamp(mix_s[gasActiveCompIdx] / mix_s[oilActiveCompIdx], EvalWell{0.}, rsmax); - } - b[oilActiveCompIdx] = FluidSystem::oilPvt().inverseFormationVolumeFactor( - pvt_region_index, temperature, seg_pressure, rs); - if (update_visc_and_den) { - phase_viscosities[oilActiveCompIdx] = FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, - seg_pressure, rs); - phase_densities[oilActiveCompIdx] = b[oilActiveCompIdx] * surface_densities_[oilActiveCompIdx] - + rs * b[oilActiveCompIdx] * surface_densities_[gasActiveCompIdx]; - } - } else { // no gas phase - b[oilActiveCompIdx] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, - temperature, seg_pressure); - if (update_visc_and_den) { - phase_viscosities[oilActiveCompIdx] = FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, - temperature, seg_pressure); - phase_densities[oilActiveCompIdx] = b[oilActiveCompIdx] * surface_densities_[oilActiveCompIdx]; - } - } - } - - mix = mix_s; - if (oilActive && gasActive) { - const EvalWell d = 1.0 - rs * rv; - if (d <= 0.0) { - const std::string str = - fmt::format("Problematic d value {} obtained for well {} during segment density calculations " - "with rs {}, rv {} and pressure {}. Continue as if no dissolution (rs = 0) and " - "vaporization (rv = 0) for this connection.", - d, well_.name(), rs, rv, seg_pressure); - deferred_logger.debug(str); - } else { - if (rs > 0.0) { - mix[gasActiveCompIdx] = (mix_s[gasActiveCompIdx] - mix_s[oilActiveCompIdx] * rs) / d; - } - if (rv > 0.0) { - mix[oilActiveCompIdx] = (mix_s[oilActiveCompIdx] - mix_s[gasActiveCompIdx] * rv) / d; - } - } - } - - for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { - vol_ratio += mix[comp_idx] / b[comp_idx]; - } -} - -template -void -MultisegmentWellSegments:: -PhaseCalcResult::clear() -{ - std::ranges::fill(b, 0.0); - std::ranges::fill(mix, 0.0); - std::ranges::fill(mix_s, 0.0); - std::ranges::fill(phase_viscosities, 0.0); - std::ranges::fill(phase_densities, 0.0); - vol_ratio = 0.0; -} - #include INSTANTIATE_TYPE_INDICES(MultisegmentWellSegments, double) diff --git a/opm/simulators/wells/MultisegmentWellSegments.hpp b/opm/simulators/wells/MultisegmentWellSegments.hpp index ef4bba94890..a3b58731ece 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.hpp +++ b/opm/simulators/wells/MultisegmentWellSegments.hpp @@ -55,11 +55,6 @@ class MultisegmentWellSegments const ParallelWellInfo& parallel_well_info, WellInterfaceGeneric& well); - void computeFluidProperties(const Scalar firstPerfTemperature, - const Scalar firstPerfSaltConcentration, - const PrimaryVariables& primary_variables, - DeferredLogger& deferred_logger); - void updateFluidProperties(const std::vector>& phase_densities, const std::vector>& phase_viscosities, const std::vector>& phase_fractions, @@ -77,12 +72,6 @@ class MultisegmentWellSegments Scalar getPressureDiffSegLocalPerf(const int seg, const int local_perf_index) const; - EvalWell getSurfaceVolume(const Scalar firstPerfTemperature, - const Scalar firstPerfSaltConcentration, - const PrimaryVariables& primary_variables, - const int seg_idx, - DeferredLogger& deferred_logger) const; - EvalWell getSurfaceVolume(const int seg_idx, DeferredLogger& deferred_logger) const; @@ -216,34 +205,6 @@ class MultisegmentWellSegments Scalar mixtureDensity(const int seg) const; Scalar mixtureDensityWithExponents(const int seg) const; Scalar mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const; - - // this class is used to store the result of phase property calculation - struct PhaseCalcResult { - explicit PhaseCalcResult(const std::size_t num_quantities) - : b(num_quantities, 0.0) - , mix(num_quantities, 0.0) - , mix_s(num_quantities, 0.0) - , phase_viscosities(num_quantities, 0.0) - , phase_densities(num_quantities, 0.0) - {} - - void clear(); - - std::vector b; - std::vector mix; - std::vector mix_s; - std::vector phase_viscosities; - std::vector phase_densities; - EvalWell vol_ratio{0.}; - }; - - void calculatePhaseProperties(PhaseCalcResult& result, - const EvalWell& temperature, - const EvalWell& saltConcentration, - const PrimaryVariables& primary_variables, - int seg, - bool update_visc_and_den, - DeferredLogger& deferred_logger) const; }; } // namespace Opm diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index ef31f49c7c0..f61edb76af1 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1192,6 +1192,7 @@ namespace Opm phase_densities[seg][compIdx] = fs.density(phaseIdx); densities[seg] += phase_fractions[seg][compIdx] * phase_densities[seg][compIdx]; } + // TODO: surface densities might not be needed in segments anymore. const auto& surface_densities = this->segments_.surfaceDensities(); for (int compIdx = 0; compIdx < nquantities; ++compIdx) { From 59a1a372160b25b12fa2246da77a9aaf7e7718ae Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 11 May 2026 00:14:11 +0200 Subject: [PATCH 11/11] WIP in cleaning up --- opm/simulators/wells/MultisegmentWell.hpp | 6 +- .../wells/MultisegmentWellGeneric.hpp | 3 + .../wells/MultisegmentWellSegments.cpp | 58 ------------------- .../wells/MultisegmentWell_impl.hpp | 43 ++------------ 4 files changed, 8 insertions(+), 102 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 63aff6619a2..4b8537f3e5c 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -206,7 +206,7 @@ namespace Opm { const Scalar relaxation_factor = 1.0); // computing the accumulation term for later use in well mass equations - void computeInitialSegmentFluids(const FSInfo& info, DeferredLogger& deferred_logger); + void computeInitialSegmentFluids(DeferredLogger& deferred_logger); // compute the pressure difference between the perforation and cell center void computePerfCellPressDiffs(const Simulator& simulator); @@ -311,10 +311,6 @@ namespace Opm { void updateWaterThroughput(const double dt, WellStateType& well_state) const override; - EvalWell getSegmentSurfaceVolume(const int seg_idx, - const FSInfo& info, - DeferredLogger& deferred_logger) const; - // turn on crossflow to avoid singular well equations // when the well is banned from cross-flow and the BHP is not properly initialized, // we turn on crossflow to avoid singular well equations. It can result in wrong-signed diff --git a/opm/simulators/wells/MultisegmentWellGeneric.hpp b/opm/simulators/wells/MultisegmentWellGeneric.hpp index 06fa7dc1a18..eeaeedd94ec 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.hpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.hpp @@ -72,6 +72,9 @@ class MultisegmentWellGeneric const std::vector& seg_dp) const; const WellInterfaceGeneric& baseif_; + + // surface densities of the fluid (it is better in WellInterfaceGeneric, w) + std::vector surface_densities_; }; } diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index ad4da37cdd2..7bf1890acc3 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -163,64 +163,6 @@ MultisegmentWellSegments(const int numSegments, } } -// template -// void MultisegmentWellSegments:: -// computeFluidProperties(const Scalar firstPerfTemperature, -// const Scalar firstPerfSaltConcentration, -// const PrimaryVariables& primary_variables, -// DeferredLogger& deferred_logger) -// { -// const int num_quantities = well_.numConservationQuantities(); -// PhaseCalcResult result(static_cast(num_quantities)); -// -// // \Note: we do not have salt concentration supported as primary variable yet. -// // \Note: this will change when we implement the brine equation in the multisegment well -// const EvalWell saltConcentration = firstPerfSaltConcentration; -// -// for (std::size_t seg = 0; seg < perforations_.size(); ++seg) { -// const EvalWell temperature = enable_energy ? primary_variables.getSegmentTemperature(seg) : firstPerfTemperature; -// -// calculatePhaseProperties(result, temperature, saltConcentration, -// primary_variables, seg, true, deferred_logger); -// -// phase_densities_[seg] = result.phase_densities; -// phase_viscosities_[seg] = result.phase_viscosities; -// -// const auto& mix = result.mix; -// const auto& mix_s = result.mix_s; -// const auto& b = result.b; -// const auto& volrat = result.vol_ratio; -// -// volume_ratios_[seg] = volrat; -// -// viscosities_[seg] = 0.; -// // calculate the average viscosity -// for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { -// const EvalWell fraction = mix[comp_idx] / b[comp_idx] / volrat; -// // TODO: a little more work needs to be done to handle the negative fractions here -// phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0; -// viscosities_[seg] += phase_viscosities_[seg][comp_idx] * phase_fractions_[seg][comp_idx]; -// } -// -// EvalWell density(0.0); -// for (int comp_idx = 0; comp_idx < num_quantities; ++comp_idx) { -// density += surface_densities_[comp_idx] * mix_s[comp_idx]; -// } -// densities_[seg] = density / volrat; -// -// // TODO: the enthalpy flux rate should be calculated in the similar manner. -// // calculate the mass rates -// mass_rates_[seg] = 0.; -// for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) { -// const int upwind_seg = upwinding_segments_[seg]; -// const EvalWell rate = primary_variables.getSegmentRateUpwinding(seg, -// upwind_seg, -// comp_idx); -// mass_rates_[seg] += rate * surface_densities_[comp_idx]; -// } -// } -// } - template void MultisegmentWellSegments:: updateFluidProperties(const std::vector>& phase_densities, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index f61edb76af1..551c8ec3d08 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -696,8 +696,7 @@ namespace Opm template void MultisegmentWell:: - computeInitialSegmentFluids(const FSInfo& info, - DeferredLogger& deferred_logger) + computeInitialSegmentFluids(DeferredLogger& deferred_logger) { for (int seg = 0; seg < this->numberOfSegments(); ++seg) { const Scalar surface_volume = this->segments_.getSurfaceVolume(seg, deferred_logger).value(); @@ -770,7 +769,7 @@ namespace Opm // after updating the primary variables, we need to update the segment fluid state updateSegmentFluidState(info, deferred_logger); - computeInitialSegmentFluids(info, deferred_logger); + computeInitialSegmentFluids(deferred_logger); if constexpr (has_energy) { computeInitialSegmentEnergy(); } @@ -1150,23 +1149,6 @@ namespace Opm // TODO: the concept of phases and components are rather confusing in this function. // needs to be addressed sooner or later. - // get the temperature for later use. It is only useful when we are not handling - // thermal related simulation - // basically, it is a single value for all the segments - // not sure how to handle the pvt region related to segment - // for the current approach, we use the pvt region of the first perforated cell - // although there are some text indicating using the pvt region of the lowest - // perforated cell - // TODO: later to investigate how to handle the pvt region - // - // auto info = this->getFirstPerforationFluidStateInfo(simulator); - // const Scalar firstPerfTemperature = std::get<0>(info); - // const Scalar firstPerfSaltConcentration = std::get<1>(info); - // - // this->segments_.computeFluidProperties(firstPerfTemperature, - // firstPerfSaltConcentration, - // this->primary_variables_, - // deferred_logger); const int nseg = this->numberOfSegments(); const int nquantities = this->numConservationQuantities(); const auto nseg_size = static_cast(nseg); @@ -1185,6 +1167,7 @@ namespace Opm } const int compIdx = FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx)); + // there was some code to handle negative fractions before the refactoring of this part phase_fractions[seg][compIdx] = fs.saturation(phaseIdx); phase_viscosities[seg][compIdx] = fs.viscosity(phaseIdx); viscosities[seg] += phase_fractions[seg][compIdx] * phase_viscosities[seg][compIdx]; @@ -2193,25 +2176,7 @@ namespace Opm - template - typename MultisegmentWell::EvalWell - MultisegmentWell:: - getSegmentSurfaceVolume(const int seg_idx, - const FSInfo& info, - DeferredLogger& deferred_logger) const - { - const Scalar firstPerfTemperature = std::get<0>(info); - const Scalar firstPerfSaltConcentration = std::get<1>(info); - - return this->segments_.getSurfaceVolume(firstPerfTemperature, - firstPerfSaltConcentration, - this->primary_variables_, - seg_idx, - deferred_logger); - } - - - template + template std::optional::Scalar> MultisegmentWell:: computeBhpAtThpLimitProd(const WellStateType& well_state,