compiling version for blackoil

This commit is contained in:
Halvor M Nilsen 2025-02-14 14:10:44 +01:00
parent 5cdde7230e
commit 3027dc6306
7 changed files with 198 additions and 125 deletions

View File

@ -436,9 +436,9 @@ 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);
}
}
}
@ -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());

View File

@ -111,46 +111,51 @@ 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.) {
assert(ws.initializedFromReservoir());
//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;
//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;
value_[seg][WFrac] = ws.multiseg_primaryvar[seg][WFrac];
}
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;
//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] = ws.multiseg_primaryvar[seg][GFrac];
}
} else { // total_seg_rate == 0
if (well_.isInjector()) {
// only single phase injection handled
auto phase = well.getInjectionProperties().injectorType;
if (has_wfrac_variable) {
if (phase == InjectorType::WATER) {
value_[seg][WFrac] = 1.0;
} else {
value_[seg][WFrac] = 0.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;
if (has_gfrac_variable) {
if (phase == InjectorType::GAS) {
value_[seg][GFrac] = 1.0;
} else {
value_[seg][GFrac] = 0.0;
}
}
// if (has_wfrac_variable) {
// if (phase == InjectorType::WATER) {
// value_[seg][WFrac] = 1.0;
// } else {
// value_[seg][WFrac] = 0.0;
// }
// }
} else if (well_.isProducer()) { // producers
if (has_wfrac_variable) {
value_[seg][WFrac] = 1.0 / well_.numPhases();
}
// if (has_gfrac_variable) {
// if (phase == InjectorType::GAS) {
// value_[seg][GFrac] = 1.0;
// } else {
// value_[seg][GFrac] = 0.0;
// }
// }
if (has_gfrac_variable) {
value_[seg][GFrac] = 1.0 / well_.numPhases();
}
}
}
// } else if (well_.isProducer()) { // producers
// if (has_wfrac_variable) {
// value_[seg][WFrac] = 1.0 / well_.numPhases();
// }
// if (has_gfrac_variable) {
// value_[seg][GFrac] = 1.0 / well_.numPhases();
// }
// }
// }
}
}
@ -221,13 +226,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;

View File

@ -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;

View File

@ -76,6 +76,9 @@ public:
serializer(production_cmode);
serializer(filtrate_conc);
serializer(perf_data);
serializer(stw_primaryvar);
serializer(multiseg_primaryvar);
serializer(initialized_from_reservoir_);
}
bool operator==(const SingleWellState&) const;
@ -142,8 +145,13 @@ public:
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;
};

View File

@ -163,63 +163,67 @@ update(const WellState<Scalar>& well_state,
}
}
if (std::abs(total_well_rate) > 0.) {
//if (std::abs(total_well_rate) > 0.) {
assert(ws.initializedFromReservoir());
if constexpr (has_wfrac_variable) {
value_[WFrac] = well_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
//value_[WFrac] = well_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
value_[WFrac] = ws.stw_primaryvar[WFrac];
}
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 ;
// 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 ;
value_[GFrac] = ws.stw_primaryvar[GFrac];
}
if constexpr (Indices::enableSolvent) {
value_[SFrac] = well_.scalingFactor(Indices::contiSolventEqIdx) * ws.sum_solvent_rates() / total_well_rate ;
value_[SFrac] = ws.stw_primaryvar[SFrac];
//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
if constexpr (has_wfrac_variable) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
auto phase = well_.wellEcl().getInjectionProperties().injectorType;
if (phase == InjectorType::WATER) {
value_[WFrac] = 1.0;
} else {
value_[WFrac] = 0.0;
}
}
}
if constexpr (has_gfrac_variable) {
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
auto phase = well_.wellEcl().getInjectionProperties().injectorType;
if (phase == InjectorType::GAS) {
value_[GFrac] = (1.0 - well_.rsRvInj());
if constexpr (Indices::enableSolvent) {
value_[GFrac] = 1.0 - well_.rsRvInj() - well_.wsolvent();
value_[SFrac] = well_.wsolvent();
}
} else {
value_[GFrac] = 0.0;
}
}
}
// } else { // total_well_rate == 0
// if (well_.isInjector()) {
// // only single phase injection handled
// if constexpr (has_wfrac_variable) {
// if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
// auto phase = well_.wellEcl().getInjectionProperties().injectorType;
// if (phase == InjectorType::WATER) {
// value_[WFrac] = 1.0;
// } else {
// value_[WFrac] = 0.0;
// }
// }
// }
// if constexpr (has_gfrac_variable) {
// if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// auto phase = well_.wellEcl().getInjectionProperties().injectorType;
// if (phase == InjectorType::GAS) {
// value_[GFrac] = (1.0 - well_.rsRvInj());
// if constexpr (Indices::enableSolvent) {
// value_[GFrac] = 1.0 - well_.rsRvInj() - well_.wsolvent();
// value_[SFrac] = well_.wsolvent();
// }
// } else {
// value_[GFrac] = 0.0;
// }
// }
// }
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_.isProducer()) { // producers
// TODO: the following are not addressed for the solvent case yet
if constexpr (has_wfrac_variable) {
value_[WFrac] = 1.0 / np;
}
// // TODO: it is possible to leave injector as a oil well,
// // when F_w and F_g both equals to zero, not sure under what kind of circumstance
// // this will happen.
// } else if (well_.isProducer()) { // producers
// // TODO: the following are not addressed for the solvent case yet
// if constexpr (has_wfrac_variable) {
// value_[WFrac] = 1.0 / np;
// }
if constexpr (has_gfrac_variable) {
value_[GFrac] = 1.0 / np;
}
} else {
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
}
}
// if constexpr (has_gfrac_variable) {
// value_[GFrac] = 1.0 / np;
// }
// } else {
// OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
// }
// }
// BHP
value_[Bhp] = ws.bhp;
@ -330,7 +334,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;
@ -393,8 +396,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

View File

@ -342,9 +342,14 @@ public:
/// 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,

View File

@ -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);
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) {
//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>::