rename to implicit_ipr

This commit is contained in:
Stein Krogstad 2023-11-07 13:31:12 +01:00
parent 746e05db5d
commit 3e6732a67b
5 changed files with 62 additions and 71 deletions

View File

@ -1326,16 +1326,16 @@ namespace Opm
updateIPR(ebos_simulator, deferred_logger);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
ws.ipr_a[idx] = this->ipr_a_[comp_idx];
ws.ipr_b[idx] = this->ipr_b_[comp_idx];
ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
}
return;
*/
}
const auto& group_state = ebos_simulator.problem().wellModel().groupState();
std::fill(ws.ipr_a.begin(), ws.ipr_a.end(), 0.);
std::fill(ws.ipr_b.begin(), ws.ipr_b.end(), 0.);
std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
//WellState well_state_copy = well_state;
auto inj_controls = Well::InjectionControls(0);
auto prod_controls = Well::ProductionControls(0);
@ -1358,9 +1358,9 @@ namespace Opm
const EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
for (size_t pvIdx = 0; pvIdx < num_eq; ++pvIdx) {
ws.ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
}
ws.ipr_a[idx] = ws.ipr_b[idx]*ws.bhp - comp_rate.value();
ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
}
// reset cmode
ws.production_cmode = cmode;

View File

@ -42,8 +42,8 @@ SingleWellState::SingleWellState(const std::string& name_,
, temperature(temp)
, well_potentials(pu_.num_phases)
, productivity_index(pu_.num_phases)
, ipr_a(pu_.num_phases)
, ipr_b(pu_.num_phases)
, implicit_ipr_a(pu_.num_phases)
, implicit_ipr_b(pu_.num_phases)
, surface_rates(pu_.num_phases)
, reservoir_rates(pu_.num_phases)
, prev_surface_rates(pu_.num_phases)
@ -91,8 +91,8 @@ void SingleWellState::shut() {
std::fill(this->prev_surface_rates.begin(), this->prev_surface_rates.end(), 0);
std::fill(this->reservoir_rates.begin(), this->reservoir_rates.end(), 0);
std::fill(this->productivity_index.begin(), this->productivity_index.end(), 0);
std::fill(this->ipr_a.begin(), this->ipr_a.end(), 0);
std::fill(this->ipr_b.begin(), this->ipr_b.end(), 0);
std::fill(this->implicit_ipr_a.begin(), this->implicit_ipr_a.end(), 0);
std::fill(this->implicit_ipr_b.begin(), this->implicit_ipr_b.end(), 0);
auto& connpi = this->perf_data.prod_index;
connpi.assign(connpi.size(), 0);
@ -309,8 +309,8 @@ bool SingleWellState::operator==(const SingleWellState& rhs) const
this->phase_mixing_rates == rhs.phase_mixing_rates &&
this->well_potentials == rhs.well_potentials &&
this->productivity_index == rhs.productivity_index &&
this->ipr_a == rhs.ipr_a &&
this->ipr_b == rhs.ipr_b &&
this->implicit_ipr_a == rhs.implicit_ipr_a &&
this->implicit_ipr_b == rhs.implicit_ipr_b &&
this->surface_rates == rhs.surface_rates &&
this->reservoir_rates == rhs.reservoir_rates &&
this->prev_surface_rates == rhs.prev_surface_rates &&

View File

@ -61,8 +61,8 @@ public:
serializer(phase_mixing_rates);
serializer(well_potentials);
serializer(productivity_index);
serializer(ipr_a);
serializer(ipr_b);
serializer(implicit_ipr_a);
serializer(implicit_ipr_b);
serializer(surface_rates);
serializer(reservoir_rates);
serializer(prev_surface_rates);
@ -100,8 +100,8 @@ public:
std::vector<double> well_potentials;
std::vector<double> productivity_index;
std::vector<double> ipr_a;
std::vector<double> ipr_b;
std::vector<double> implicit_ipr_a;
std::vector<double> implicit_ipr_b;
std::vector<double> surface_rates;
std::vector<double> reservoir_rates;
std::vector<double> prev_surface_rates;

View File

@ -843,7 +843,6 @@ namespace Opm
// For a component rate r the derivative dr/dbhp is obtained by
// dr/dbhp = - (partial r/partial x) * inv(partial Eq/partial x) * (partial Eq/partial control_value)
// where Eq(x)=0 is the well equation setup with bhp control and primary varables x
//StandardWell<TypeTag> well_copy(*this);
// We shouldn't have zero rates at this stage, but check
bool zero_rates;
@ -854,66 +853,58 @@ namespace Opm
}
auto& ws = well_state.well(this->index_of_well_);
if (zero_rates) {
const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, reverting to explicit IPR-calulations", this->name());
const auto msg = fmt::format("updateIPRImplicit: Well {} has zero rate, IPRs might be probelmatic", this->name());
deferred_logger.debug(msg);
updateIPR(ebosSimulator, deferred_logger);
/*
// could revert to standard approach here
updateIPR(ebos_simulator, deferred_logger);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
ws.ipr_a[idx] = this->ipr_a_[comp_idx];
ws.ipr_b[idx] = this->ipr_b_[comp_idx];
ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
}
} else {
const auto& group_state = ebosSimulator.problem().wellModel().groupState();
return;
*/
}
const auto& group_state = ebosSimulator.problem().wellModel().groupState();
// XXX maybe don't update this
std::fill(ws.ipr_a.begin(), ws.ipr_a.end(), 0.);
std::fill(ws.ipr_b.begin(), ws.ipr_b.end(), 0.);
//WellState well_state_copy = well_state;
auto inj_controls = Well::InjectionControls(0);
auto prod_controls = Well::ProductionControls(0);
prod_controls.addControl(Well::ProducerCMode::BHP);
prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
std::fill(ws.implicit_ipr_a.begin(), ws.implicit_ipr_a.end(), 0.);
std::fill(ws.implicit_ipr_b.begin(), ws.implicit_ipr_b.end(), 0.);
auto inj_controls = Well::InjectionControls(0);
auto prod_controls = Well::ProductionControls(0);
prod_controls.addControl(Well::ProducerCMode::BHP);
prod_controls.bhp_limit = well_state.well(this->index_of_well_).bhp;
// Set current control to bhp, and bhp value in state, modify bhp limit in control object.
const auto cmode = ws.production_cmode;
ws.production_cmode = Well::ProducerCMode::BHP;
const double dt = ebosSimulator.timeStepSize();
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
// Set current control to bhp, and bhp value in state, modify bhp limit in control object.
const auto cmode = ws.production_cmode;
ws.production_cmode = Well::ProducerCMode::BHP;
const double dt = ebosSimulator.timeStepSize();
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
const double nEq = this->primary_variables_.numWellEq();
BVectorWell rhs(1);
rhs[0].resize(nEq);
// rhs = 0 except -1 for control eq
for (size_t i=0; i < nEq; ++i){
rhs[0][i] = 0.0;
}
rhs[0][Bhp] = -1.0;
BVectorWell x_well(1);
x_well[0].resize(nEq);
this->linSys_.solve(rhs, x_well);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
ws.ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
}
ws.ipr_a[idx] = ws.ipr_b[idx]*ws.bhp - comp_rate.value();
//for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
// this->ipr_b_[comp_idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
//}
// XXX maybe don't update this
//this->ipr_a_[comp_idx] = this->ipr_b_[comp_idx]*ws.bhp - comp_rate.value();
// For ipr in well_state use same ordering as potentials etc.
//const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
//ws.ipr_a[idx] = this->ipr_a_[comp_idx];
//ws.ipr_b[idx] = this->ipr_b_[comp_idx];
}
// reset cmode
ws.production_cmode = cmode;
const double nEq = this->primary_variables_.numWellEq();
BVectorWell rhs(1);
rhs[0].resize(nEq);
// rhs = 0 except -1 for control eq
for (size_t i=0; i < nEq; ++i){
rhs[0][i] = 0.0;
}
rhs[0][Bhp] = -1.0;
BVectorWell x_well(1);
x_well[0].resize(nEq);
this->linSys_.solve(rhs, x_well);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
const int idx = this->ebosCompIdxToFlowCompIdx(comp_idx);
for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
ws.implicit_ipr_b[idx] -= x_well[0][pvIdx]*comp_rate.derivative(pvIdx+Indices::numEq);
}
ws.implicit_ipr_a[idx] = ws.implicit_ipr_b[idx]*ws.bhp - comp_rate.value();
}
// reset cmode
ws.production_cmode = cmode;
}
template<typename TypeTag>

View File

@ -97,11 +97,11 @@ public:
std::pair<double,double>retval(0.0, 0.0);
const VFPProdTable& table = this->m_prod.getTable(table_id);
const auto& pu = well_state_.phaseUsage();
const auto& ipr_a= well_state_.well(well_index).ipr_a;
const auto& ipr_a= well_state_.well(well_index).implicit_ipr_a;
const auto& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
const auto& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]]:0.0;
const auto& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]]:0.0;
const auto& ipr_b= well_state_.well(well_index).ipr_b;
const auto& ipr_b= well_state_.well(well_index).implicit_ipr_b;
const auto& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
const auto& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]]:0.0;
const auto& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]]:0.0;