Skip to content

Store primary variables in well state, and improve well initialization #6114

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
GitPaean merged 7 commits into from
May 8, 2025
Merged

Store primary variables in well state, and improve well initialization #6114

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
2f1e445
compiling version for blackoil
hnil Feb 14, 2025
2e81e67
Updates to clean up.
Mar 28, 2025
e77484b
Refine producer initialization.
Apr 2, 2025
a828ce3
Ensure that the initialization is not producing above-bhp-limit rates.
Apr 29, 2025
463dfcf
Cleanup in response to review comments.
May 6, 2025
272fb20
Do not re-initialize stopped wells.
May 7, 2025
1e6c89c
Improve treatment of primary variables.
May 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 4 additions & 6 deletions opm/simulators/wells/BlackoilWellModel_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -429,16 +429,14 @@ namespace Opm {
// we need the inj_multiplier from the previous time step
this->initInjMult();

const auto& summaryState = simulator_.vanguard().summaryState();
if (alternative_well_rate_init_) {
// Update the well rates of well_state_, if only single-phase rates, to
// have proper multi-phase rates proportional to rates at bhp zero.
// This is done only for producers, as injectors will only have a single
// nonzero phase anyway.
for (const auto& well : well_container_) {
const bool zero_target = well->stoppedOrZeroRateTarget(simulator_, this->wellState(), local_deferredLogger);
if (well->isProducer() && !zero_target) {
well->updateWellStateRates(simulator_, this->wellState(), local_deferredLogger);
if (well->isProducer()) {
well->initializeProducerWellState(simulator_, this->wellState(), local_deferredLogger);

@GitPaean GitPaean May 7, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removing && !zero_target caused the regression for the case, it is part of the difficulties that difficult to describe the flow fractions for zero rate situation when we use only surface rates to initialize the primary variables. We should have fraction information, basically zero total rate and proper fraction for the phases.

11_udt-1d-03
some bumps at the end
image

@GitPaean GitPaean May 7, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not saying the master branch is well designed, I believe the initialized rates play difference in the following logic with nonzero_rate_original.

            if (converged) {
                const bool zero_target = this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
                if (this->wellIsStopped() && !zero_target && nonzero_rate_original) {
                    // Well had non-zero rate, but was stopped during local well-solve. We re-open the well 
                    // for the next global iteration, but if the zero rate persists, it will be stopped.
                    // This logic is introduced to prevent/ameliorate stopped/revived oscillations  
                    this->operability_status_.resetOperability();
                    this->openWell();
                    deferred_logger.debug("    " + this->name() + " is re-opened after being stopped during local solve");
                }
            }

It will open a well previously stopped due to not solvable or operable.

}
}
}
Expand Down Expand Up @@ -466,7 +464,7 @@ namespace Opm {
const Group& fieldGroup = this->schedule().getGroup("FIELD", reportStepIdx);
WellGroupHelpers<Scalar>::updateGuideRates(fieldGroup,
this->schedule(),
summaryState,
this->summaryState(),
this->phase_usage_,
reportStepIdx,
simulationTime,
Expand Down Expand Up @@ -559,7 +557,7 @@ namespace Opm {

// initialize rates/previous rates to prevent zero fractions in vfp-interpolation
if (well->isProducer()) {
well->updateWellStateRates(simulator_, this->wellState(), deferred_logger);
well->initializeProducerWellState(simulator_, this->wellState(), deferred_logger);
}
if (well->isVFPActive(deferred_logger)) {
well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
Expand Down
20 changes: 19 additions & 1 deletion opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,16 @@ update(const WellState<Scalar>& well_state,
if (stop_or_zero_rate_target && seg == 0) {
value_[seg][WQTotal] = 0;
}
if (std::abs(total_seg_rate) > 0.) {

// todo to map old fraction to new perforations for now start from scratch.
if (ws.primaryvar.size() == value_.size() * numWellEq) {
if (has_wfrac_variable) {
value_[seg][WFrac] = ws.primaryvar[seg * numWellEq + WFrac];
}
if (has_gfrac_variable) {
value_[seg][GFrac] = ws.primaryvar[seg * numWellEq + GFrac];
}
} else if (std::abs(total_seg_rate) > 0.) {
if (has_wfrac_variable) {
const int water_pos = pu.phase_pos[Water];
value_[seg][WFrac] = well_.scalingFactor(water_pos) * segment_rates[well_.numPhases() * seg + water_pos] / total_seg_rate;
Expand Down Expand Up @@ -230,6 +239,15 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
const int oil_pos = pu.phase_pos[Oil];

auto& ws = well_state.well(well_.indexOfWell());

// Store primary variables
ws.primaryvar.resize(value_.size() * numWellEq);
for (std::size_t seg = 0; seg < value_.size(); ++seg) {
for (int ii = 0; ii < numWellEq; ++ii) {
ws.primaryvar[seg * numWellEq + ii] = value_[seg][ii];
}
}

auto& segments = ws.segments;
auto& segment_rates = segments.rates;
auto& disgas = segments.dissolved_gas_rate;
Expand Down
16 changes: 7 additions & 9 deletions opm/simulators/wells/SingleWellState.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,6 +246,7 @@ update_producer_targets(const Well& ecl_well, const SummaryState& st)
return;
}

std::fill(this->surface_rates.begin(), this->surface_rates.end(), 0.0);
switch (prod_controls.cmode) {
case Well::ProducerCMode::ORAT:
assert(this->pu.phase_used[BlackoilPhases::Liquid]);
Expand All @@ -262,15 +263,12 @@ update_producer_targets(const Well& ecl_well, const SummaryState& st)
case Well::ProducerCMode::GRUP:
case Well::ProducerCMode::THP:
case Well::ProducerCMode::BHP:
if (this->pu.phase_used[BlackoilPhases::Liquid]) {
this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
}
if (this->pu.phase_used[BlackoilPhases::Aqua]) {
this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
}
if (this->pu.phase_used[BlackoilPhases::Vapour]){
this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -100000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
}
// Keeping all rates at zero, they will be initialized properly in
// a call to WellInterface::initializeProducerWellState() later, which will
// use the reservoir state to find a better initial value.
// This also applies to the ORAT/WRAT/GRAT above, but then only the
// rate not set in the above will be modified in
// WellInterface::initializeProducerWellState().
break;

default:
Expand Down
5 changes: 3 additions & 2 deletions opm/simulators/wells/SingleWellState.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ class SingleWellState {
serializer(well_potentials);
serializer(productivity_index);
serializer(implicit_ipr_a);
serializer(implicit_ipr_b);
serializer(implicit_ipr_b);
serializer(surface_rates);
serializer(reservoir_rates);
serializer(prev_surface_rates);
Expand All @@ -76,6 +76,7 @@ class SingleWellState {
serializer(production_cmode);
serializer(filtrate_conc);
serializer(perf_data);
serializer(primaryvar);
}

bool operator==(const SingleWellState&) const;
Expand Down Expand Up @@ -115,7 +116,7 @@ class SingleWellState {
Events events;
WellInjectorCMode injection_cmode{WellInjectorCMode::CMODE_UNDEFINED};
WellProducerCMode production_cmode{WellProducerCMode::CMODE_UNDEFINED};

std::vector<Scalar> primaryvar;

/// Special purpose method to support dynamically rescaling a well's
/// CTFs through WELPI.
Expand Down
14 changes: 12 additions & 2 deletions opm/simulators/wells/StandardWellPrimaryVariables.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,17 @@ update(const WellState<Scalar>& well_state,
}
}

if (std::abs(total_well_rate) > 0.) {
if (ws.primaryvar.size() > 0) {
if constexpr (has_wfrac_variable) {
value_[WFrac] = ws.primaryvar[WFrac];
}
if constexpr (has_gfrac_variable) {
value_[GFrac] = ws.primaryvar[GFrac];
}
if constexpr (Indices::enableSolvent) {
value_[SFrac] = ws.primaryvar[SFrac];
}
} else if (std::abs(total_well_rate) > 0.) {
if constexpr (has_wfrac_variable) {
value_[WFrac] = well_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
}
Expand All @@ -175,7 +185,6 @@ update(const WellState<Scalar>& well_state,
if constexpr (Indices::enableSolvent) {
value_[SFrac] = well_.scalingFactor(Indices::contiSolventEqIdx) * ws.sum_solvent_rates() / total_well_rate ;
}

} else { // total_well_rate == 0
if (well_.isInjector()) {
// only single phase injection handled
Expand Down Expand Up @@ -399,6 +408,7 @@ copyToWellState(WellState<Scalar>& well_state,
}

auto& ws = well_state.well(well_.indexOfWell());
ws.primaryvar = value_;
ws.bhp = value_[Bhp];

// calculate the phase rates based on the primary variables
Expand Down
8 changes: 4 additions & 4 deletions opm/simulators/wells/WellInterface.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -327,10 +327,10 @@ class WellInterface : public WellInterfaceIndices<GetPropType<TypeTag, Propertie

/// Modify the well_state's rates if there is only one nonzero rate.
/// If so, that rate is kept as is, but the others are set proportionally
/// to the rates returned by computeCurrentWellRates().
void updateWellStateRates(const Simulator& simulator,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const;
/// to the rates at bhp limit or bhp 1 bar.
void initializeProducerWellState(const Simulator& simulator,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const;

void solveWellEquation(const Simulator& simulator,
WellState<Scalar>& well_state,
Expand Down
65 changes: 51 additions & 14 deletions opm/simulators/wells/WellInterface_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1603,10 +1603,11 @@ namespace Opm
template <typename TypeTag>
void
WellInterface<TypeTag>::
updateWellStateRates(const Simulator& simulator,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const
initializeProducerWellState(const Simulator& simulator,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const
{
assert(this->isProducer());
OPM_TIMEFUNCTION();
// Check if the rates of this well only are single-phase, do nothing
// if more than one nonzero rate.
Expand All @@ -1624,29 +1625,65 @@ namespace Opm
}
}

// Calculate the rates that follow from the current primary variables.
std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
// Calculate rates at bhp limit, or 1 bar if no limit.
std::vector<Scalar> well_q_s(this->number_of_phases_, 0.0);
bool rates_evaluated_at_1bar = false;
{
const auto& summary_state = simulator.vanguard().summaryState();
const auto& prod_controls = this->well_ecl_.productionControls(summary_state);
const double bhp_limit = std::max(prod_controls.bhp_limit, 1.0 * unit::barsa);
this->computeWellRatesWithBhp(simulator, bhp_limit, well_q_s, deferred_logger);
// Remember of we evaluated the rates at (approx.) 1 bar or not.
rates_evaluated_at_1bar = (bhp_limit < 1.1 * unit::barsa);
// Check that no rates are positive.
if (std::any_of(well_q_s.begin(), well_q_s.end(), [](Scalar q) { return q > 0.0; })) {
// Did we evaluate at 1 bar? If not, then we can try again at 1 bar.
if (!rates_evaluated_at_1bar) {
this->computeWellRatesWithBhp(simulator, 1.0 * unit::barsa, well_q_s, deferred_logger);
rates_evaluated_at_1bar = true;
}
// At this point we can only set the wrong-direction (if any) values to zero.
for (auto& q : well_q_s) {
q = std::min(q, Scalar{0.0});
}
}
}

if (nonzero_rate_index == -1) {
// No nonzero rates.
// Use the computed rate directly
// No nonzero rates on input.
// Use the computed rate directly, or scaled by a factor
// 0.5 (to avoid too high values) if it was evaluated at 1 bar.
const Scalar factor = rates_evaluated_at_1bar ? 0.5 : 1.0;
for (int p = 0; p < this->number_of_phases_; ++p) {
ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
ws.surface_rates[p] = factor * well_q_s[p];
}
return;
}

// Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
// If we are here, we had a single nonzero rate for the well,
// typically from a rate constraint. We must make sure it is
// respected, so if it was lower than the calculated rate for
// the same phase we scale all rates to match.
const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
const int comp_idx_nz = this->flowPhaseToModelCompIdx(nonzero_rate_index);
if (std::abs(well_q_s[comp_idx_nz]) > floating_point_error_epsilon) {
const Scalar computed_rate = well_q_s[nonzero_rate_index];
if (std::abs(initial_nonzero_rate) < std::abs(computed_rate)) {

@GitPaean GitPaean May 6, 2025 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we need to check std::abs(initial_nonzero_rate) < std::abs(computed_rate)? It might save some small computation potentially, while it adds more complication to the code. The results should be the same?

@atgeirr atgeirr May 6, 2025

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If initial_nonzero_rate is greater than computed_rate we should not scale everything to match that rate, because the computed_rate comes from a bhp limit and we should not go above that rate (i.e. factor would be > 1).

An example: we have a GRAT control set to a high number as an important constraint for the later parts of the simulation, but at the start of simulation the gas rate is low. Then we might scale the rate up to be much higher than the BHP-limit rate.

// Note that both rates below are negative. The factor should be < 1.0.
const Scalar factor = initial_nonzero_rate / computed_rate;
assert(factor < 1.0);
for (int p = 0; p < this->number_of_phases_; ++p) {
// We skip the nonzero_rate_index, as that should remain as it was.
if (p != nonzero_rate_index) {
const int comp_idx = this->flowPhaseToModelCompIdx(p);
Scalar& rate = ws.surface_rates[p];
rate = (initial_nonzero_rate / well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
ws.surface_rates[p] = factor * well_q_s[p];
}
}
return;
}

// If we are here, we had a single nonzero rate, but it was
// higher than the one calculated from the bhp limit, so we
// use the calculated rates.
for (int p = 0; p < this->number_of_phases_; ++p) {
ws.surface_rates[p] = well_q_s[p];
}
}

Expand Down
1 change: 1 addition & 0 deletions opm/simulators/wells/WellState.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -409,6 +409,7 @@ void WellState<Scalar>::init(const std::vector<Scalar>& cellPressures,
}

new_well.surface_rates = prev_well.surface_rates;
new_well.prev_surface_rates = prev_well.prev_surface_rates;
new_well.reservoir_rates = prev_well.reservoir_rates;
new_well.well_potentials = prev_well.well_potentials;

Expand Down