Skip to content
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

Safer well state #5990

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

Safer well state #5990

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3027dc6
compiling version for blackoil
hnil Feb 14, 2025
3e3b72d
hopefully better initialization of wells
hnil Feb 14, 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: 5 additions & 5 deletions opm/simulators/wells/BlackoilWellModel_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -436,10 +436,10 @@ namespace Opm {
// 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);
}
//const bool zero_target = well->stoppedOrZeroRateTarget(simulator_, this->wellState(), local_deferredLogger);
//if (well->isProducer()){// && !zero_target) {
well->initializeWellState(simulator_, this->groupState(), this->wellState(), local_deferredLogger);
//}
}
}

Expand Down Expand Up @@ -559,7 +559,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->initializeWellState(simulator_, this->groupState(), this->wellState(), deferred_logger);
}
if (well->isVFPActive(deferred_logger)) {
well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
Expand Down
23 changes: 20 additions & 3 deletions opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,16 +111,29 @@ update(const WellState<Scalar>& well_state,
if (stop_or_zero_rate_target && seg == 0) {
value_[seg][WQTotal] = 0;
}
assert(ws.initializedFromReservoir());
// tot to map old fraction to new perforations for now start from scratch.
bool prim_set = ws.stw_primaryvar.size() == value_.size();
if(prim_set){
//if (std::abs(total_seg_rate) > 0.) {
if (has_wfrac_variable) {
value_[seg][WFrac] = ws.multiseg_primaryvar[seg][WFrac];
}
if (has_gfrac_variable) {
value_[seg][GFrac] = ws.multiseg_primaryvar[seg][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;
}
if (has_gfrac_variable) {
const int gas_pos = pu.phase_pos[Gas];
value_[seg][GFrac] = well_.scalingFactor(gas_pos) * segment_rates[well_.numPhases() * seg + gas_pos] / total_seg_rate;
value_[seg][GFrac] = well_.scalingFactor(gas_pos) * segment_rates[well_.numPhases() * seg + gas_pos] / total_seg_rate;
}
} else { // total_seg_rate == 0
// what about water and gas injection?
} else { // total_seg_rate == 0
if (well_.isInjector()) {
// only single phase injection handled
auto phase = well.getInjectionProperties().injectorType;
Expand Down Expand Up @@ -151,6 +164,7 @@ update(const WellState<Scalar>& well_state,
}
}
}
}
}
}

Expand Down Expand Up @@ -221,13 +235,16 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;

const auto pvtReg = std::max(well_.wellEcl().pvt_table_number() - 1, 0);

const PhaseUsage& pu = well_.phaseUsage();
const int oil_pos = pu.phase_pos[Oil];

auto& ws = well_state.well(well_.indexOfWell());
if constexpr ( numWellEq == 4) {
ws.multiseg_primaryvar = value_;
}
auto& segments = ws.segments;
auto& segment_rates = segments.rates;
auto& disgas = segments.dissolved_gas_rate;
Expand Down
1 change: 1 addition & 0 deletions opm/simulators/wells/SingleWellState.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ SingleWellState(const std::string& name_,
, prev_surface_rates(pu_.num_phases)
, perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu_.num_phases)
, trivial_target(false)
, initialized_from_reservoir_(false)
{
for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
Expand Down
10 changes: 9 additions & 1 deletion opm/simulators/wells/SingleWellState.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,9 @@ class SingleWellState {
serializer(production_cmode);
serializer(filtrate_conc);
serializer(perf_data);
serializer(stw_primaryvar);
serializer(multiseg_primaryvar);
serializer(initialized_from_reservoir_);
}

bool operator==(const SingleWellState&) const;
Expand Down Expand Up @@ -142,8 +145,13 @@ class SingleWellState {

Scalar sum_filtrate_rate() const;
Scalar sum_filtrate_total() const;

std::vector<Scalar> stw_primaryvar;
std::vector<std::array<Scalar, 4>> multiseg_primaryvar;
bool initializedFromReservoir() const { return initialized_from_reservoir_; }
void setInitializedFromReservoir(bool value) { initialized_from_reservoir_ = value; }
private:
bool initialized_from_reservoir_ = false;

Scalar sum_connection_rates(const std::vector<Scalar>& connection_rates) const;
};

Expand Down
32 changes: 23 additions & 9 deletions opm/simulators/wells/StandardWellPrimaryVariables.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,20 +162,34 @@ update(const WellState<Scalar>& well_state,
value_[WQTotal] = 0.;
}
}

if (std::abs(total_well_rate) > 0.) {
assert(ws.initializedFromReservoir());
//if (std::abs(total_well_rate) > 0.) {

bool prim_set = ws.stw_primaryvar.size()>0;
if(prim_set){
if constexpr (has_wfrac_variable) {
value_[WFrac] = ws.stw_primaryvar[WFrac];
}
if constexpr (has_gfrac_variable) {
value_[GFrac] = ws.stw_primaryvar[GFrac];
}
if constexpr (Indices::enableSolvent) {
value_[SFrac] = ws.stw_primaryvar[SFrac];
}
}else{
assert(prim_set==false);
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;
}
if constexpr (has_gfrac_variable) {
value_[GFrac] = well_.scalingFactor(pu.phase_pos[Gas]) *
(ws.surface_rates[pu.phase_pos[Gas]] -
(Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ;
(ws.surface_rates[pu.phase_pos[Gas]] -
(Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ;
}
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 @@ -220,8 +234,8 @@ update(const WellState<Scalar>& well_state,
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
}
}

// BHP
}
//BHP
value_[Bhp] = ws.bhp;
}

Expand Down Expand Up @@ -330,7 +344,6 @@ copyToWellState(WellState<Scalar>& well_state,
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;

const PhaseUsage& pu = well_.phaseUsage();
std::vector<Scalar> F(well_.numPhases(), 0.0);
[[maybe_unused]] Scalar F_solvent = 0.0;
Expand Down Expand Up @@ -393,8 +406,9 @@ copyToWellState(WellState<Scalar>& well_state,
F_solvent /= well_.scalingFactor(Indices::contiSolventEqIdx);
F[pu.phase_pos[Gas]] += F_solvent;
}

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

// calculate the phase rates based on the primary variables
Expand Down
11 changes: 8 additions & 3 deletions opm/simulators/wells/WellInterface.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -342,9 +342,14 @@ 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;
void initializeWellState(const Simulator& simulator,
const GroupState<Scalar>& group_state,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const;

// void updateWellStateRates(const Simulator& simulator,
// WellState<Scalar>& well_state,
// DeferredLogger& deferred_logger) const;

void solveWellEquation(const Simulator& simulator,
WellState<Scalar>& well_state,
Expand Down
117 changes: 82 additions & 35 deletions opm/simulators/wells/WellInterface_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1622,53 +1622,100 @@ namespace Opm
template <typename TypeTag>
void
WellInterface<TypeTag>::
updateWellStateRates(const Simulator& simulator,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const
initializeWellState(const Simulator& simulator,
const GroupState<Scalar>& group_state,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const
{
OPM_TIMEFUNCTION();
// Check if the rates of this well only are single-phase, do nothing
// if more than one nonzero rate.
auto& ws = well_state.well(this->index_of_well_);
int nonzero_rate_index = -1;
const Scalar floating_point_error_epsilon = 1e-14;
for (int p = 0; p < this->number_of_phases_; ++p) {
if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
if (nonzero_rate_index == -1) {
nonzero_rate_index = p;
} else {
// More than one nonzero rate.
return;
}
}
}

// int nonzero_rate_index = -1;
// const Scalar floating_point_error_epsilon = 1e-14;
// for (int p = 0; p < this->number_of_phases_; ++p) {
// if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
// if (nonzero_rate_index == -1) {
// nonzero_rate_index = p;
// } else {
// // More than one nonzero rate.
// return;
// }
// }
// }
// Calculate the rates that follow from the current primary variables.
std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);

if (nonzero_rate_index == -1) {
std::vector<Scalar> well_q_s;
Scalar bhp_tmp = 0.0;
computeWellRatesWithBhp(simulator,
bhp_tmp,
well_q_s,
deferred_logger);
// constcomputeCurrentWellRates(simulator, deferred_logger);

//if (nonzero_rate_index == -1) {
// No nonzero rates.
// Use the computed rate directly
for (int p = 0; p < this->number_of_phases_; ++p) {
ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
}
return;
for (int p = 0; p < this->number_of_phases_; ++p) {
ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
}
// set fractions

// Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
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) {
for (int p = 0; p < this->number_of_phases_; ++p) {
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.setInitializedFromReservoir(true);
this->updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);

}


// template <typename TypeTag>
// void
// WellInterface<TypeTag>::
// updateWellStateRates(const Simulator& simulator,
// WellState<Scalar>& well_state,
// DeferredLogger& deferred_logger) const
// {
// OPM_TIMEFUNCTION();
// // Check if the rates of this well only are single-phase, do nothing
// // if more than one nonzero rate.
// auto& ws = well_state.well(this->index_of_well_);
// int nonzero_rate_index = -1;
// const Scalar floating_point_error_epsilon = 1e-14;
// for (int p = 0; p < this->number_of_phases_; ++p) {
// if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
// if (nonzero_rate_index == -1) {
// nonzero_rate_index = p;
// } else {
// // More than one nonzero rate.
// return;
// }
// }
// }

// // Calculate the rates that follow from the current primary variables.
// std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);

// if (nonzero_rate_index == -1) {
// // No nonzero rates.
// // Use the computed rate directly
// for (int p = 0; p < this->number_of_phases_; ++p) {
// ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
// }
// return;
// }
// // set fractions
// // Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
// 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) {
// for (int p = 0; p < this->number_of_phases_; ++p) {
// 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]);
// }
// }
// }
// }

template <typename TypeTag>
std::vector<typename WellInterface<TypeTag>::Scalar>
WellInterface<TypeTag>::
Expand Down