OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

make has_polymermw constexpr and use if constexpr to eliminate code

Browse Source
This commit is contained in:
Arne Morten Kvarving 2021-05-14 21:49:47 +02:00
parent dd04eb128a
commit 268201eeb1
2 changed files with 141 additions and 121 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

260
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -79,11 +79,13 @@ namespace Opm
}
// counting/updating primary variable numbers
if (this->has_polymermw && this->isInjector()) {
// adding a primary variable for water perforation rate per connection
numWellEq_ += number_of_perforations_;
// adding a primary variable for skin pressure per connection
numWellEq_ += number_of_perforations_;
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
// adding a primary variable for water perforation rate per connection
numWellEq_ += number_of_perforations_;
// adding a primary variable for skin pressure per connection
numWellEq_ += number_of_perforations_;
}
}
// with the updated numWellEq_, we can initialize the primary variables and matrices now
@ -406,10 +408,12 @@ namespace Opm
const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
EvalWell drawdown = pressure - well_pressure;
if (this->has_polymermw && this->isInjector()) {
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index];
drawdown += skin_pressure;
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index];
drawdown += skin_pressure;
}
}
// producing perforations
@ -596,8 +600,10 @@ namespace Opm
calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
// Equation assembly for this perforation.
if (has_polymer && this->has_polymermw && this->isInjector()) {
handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
if constexpr (has_polymer && Base::has_polymermw) {
if (this->isInjector()) {
handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
}
}
const int cell_idx = well_cells_[perf];
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
@ -698,14 +704,16 @@ namespace Opm
computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
if (has_polymer && this->has_polymermw && this->isInjector()) {
// Store the original water flux computed from the reservoir quantities.
// It will be required to assemble the injectivity equations.
const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
water_flux_s = cq_s[water_comp_idx];
// Modify the water flux for the rest of this function to depend directly on the
// local water velocity primary variable.
handleInjectivityRate(ebosSimulator, perf, cq_s);
if constexpr (has_polymer && Base::has_polymermw) {
if (this->isInjector()) {
// Store the original water flux computed from the reservoir quantities.
// It will be required to assemble the injectivity equations.
const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
water_flux_s = cq_s[water_comp_idx];
// Modify the water flux for the rest of this function to depend directly on the
// local water velocity primary variable.
handleInjectivityRate(ebosSimulator, perf, cq_s);
}
}
// updating the solution gas rate and solution oil rate
@ -790,7 +798,7 @@ namespace Opm
cq_s_poly *= well_efficiency_factor_;
connectionRates[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
if (this->has_polymermw) {
if constexpr (Base::has_polymermw) {
updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger);
}
}
@ -970,7 +978,11 @@ namespace Opm
OPM_DEFLOG_THROW(std::runtime_error, "Water is required when polymer is active", deferred_logger);
}
updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
// for the cases related to polymer molecular weight, we assume fully mixing
// as a result, the polymer and water share the same viscosity
if constexpr (!Base::has_polymermw) {
updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
}
}
}
@ -1068,17 +1080,19 @@ namespace Opm
updateExtraPrimaryVariables(const BVectorWell& dwells) const
{
// for the water velocity and skin pressure
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const int wat_vel_index = Bhp + 1 + perf;
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const int wat_vel_index = Bhp + 1 + perf;
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const double relaxation_factor = 0.9;
const double dx_wat_vel = dwells[0][wat_vel_index];
primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
const double relaxation_factor = 0.9;
const double dx_wat_vel = dwells[0][wat_vel_index];
primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
const double dx_pskin = dwells[0][pskin_index];
primary_variables_[pskin_index] -= relaxation_factor * dx_pskin;
const double dx_pskin = dwells[0][pskin_index];
primary_variables_[pskin_index] -= relaxation_factor * dx_pskin;
}
}
}
}
@ -1284,10 +1298,12 @@ namespace Opm
updateThp(well_state, deferred_logger);
// other primary variables related to polymer injectivity study
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
well_state.perfWaterVelocity()[first_perf_ + perf] = primary_variables_[Bhp + 1 + perf];
well_state.perfSkinPressure()[first_perf_ + perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf];
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
well_state.perfWaterVelocity()[first_perf_ + perf] = primary_variables_[Bhp + 1 + perf];
well_state.perfSkinPressure()[first_perf_ + perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf];
}
}
}
}
@ -2787,10 +2803,12 @@ namespace Opm
primary_variables_[Bhp] = well_state.bhp()[index_of_well_];
// other primary variables related to polymer injection
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
primary_variables_[Bhp + 1 + perf] = well_state.perfWaterVelocity()[first_perf_ + perf];
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
primary_variables_[Bhp + 1 + perf] = well_state.perfWaterVelocity()[first_perf_ + perf];
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
}
}
}
#ifndef NDEBUG
@ -2862,11 +2880,6 @@ namespace Opm
std::vector<EvalWell>& mob,
Opm::DeferredLogger& deferred_logger) const
{
// for the cases related to polymer molecular weight, we assume fully mixing
// as a result, the polymer and water share the same viscosity
if (this->has_polymermw) {
return;
}
const int cell_idx = well_cells_[perf];
const auto& int_quant = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const EvalWell polymer_concentration = extendEval(int_quant.polymerConcentration());
@ -3063,20 +3076,21 @@ namespace Opm
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const
{
if (!this->has_polymermw) {
if constexpr (Base::has_polymermw) {
const int water_table_id = well_ecl_.getPolymerProperties().m_skprwattable;
if (water_table_id <= 0) {
OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger);
}
const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_water(numWellEq_ + numEq, 0.0);
pskin_water = water_table_func.eval(throughput_eval, water_velocity);
return pskin_water;
} else {
OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
"while injecting skin pressure is requested for well " << name(), deferred_logger);
}
const int water_table_id = well_ecl_.getPolymerProperties().m_skprwattable;
if (water_table_id <= 0) {
OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger);
}
const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_water(numWellEq_ + numEq, 0.0);
pskin_water = water_table_func.eval(throughput_eval, water_velocity);
return pskin_water;
}
@ -3091,32 +3105,33 @@ namespace Opm
const EvalWell& poly_inj_conc,
Opm::DeferredLogger& deferred_logger) const
{
if (!this->has_polymermw) {
if constexpr (Base::has_polymermw) {
const double sign = water_velocity >= 0. ? 1.0 : -1.0;
const EvalWell water_velocity_abs = Opm::abs(water_velocity);
if (poly_inj_conc == 0.) {
return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
}
const int polymer_table_id = well_ecl_.getPolymerProperties().m_skprpolytable;
if (polymer_table_id <= 0) {
OPM_DEFLOG_THROW(std::runtime_error, "Unavailable SKPRPOLY table id used for well " << name(), deferred_logger);
}
const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
const double reference_concentration = skprpolytable.refConcentration;
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_poly(numWellEq_ + numEq, 0.0);
pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
if (poly_inj_conc == reference_concentration) {
return sign * pskin_poly;
}
// poly_inj_conc != reference concentration of the table, then some interpolation will be required
const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
return sign * pskin;
} else {
OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
"while injecting skin pressure is requested for well " << name(), deferred_logger);
}
const double sign = water_velocity >= 0. ? 1.0 : -1.0;
const EvalWell water_velocity_abs = Opm::abs(water_velocity);
if (poly_inj_conc == 0.) {
return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
}
const int polymer_table_id = well_ecl_.getPolymerProperties().m_skprpolytable;
if (polymer_table_id <= 0) {
OPM_DEFLOG_THROW(std::runtime_error, "Unavailable SKPRPOLY table id used for well " << name(), deferred_logger);
}
const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
const double reference_concentration = skprpolytable.refConcentration;
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_poly(numWellEq_ + numEq, 0.0);
pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
if (poly_inj_conc == reference_concentration) {
return sign * pskin_poly;
}
// poly_inj_conc != reference concentration of the table, then some interpolation will be required
const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
return sign * pskin;
}
@ -3130,19 +3145,20 @@ namespace Opm
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const
{
if (!this->has_polymermw) {
if constexpr (Base::has_polymermw) {
const int table_id = well_ecl_.getPolymerProperties().m_plymwinjtable;
const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
EvalWell molecular_weight(numWellEq_ + numEq, 0.);
if (wpolymer() == 0.) { // not injecting polymer
return molecular_weight;
}
molecular_weight = table_func.eval(throughput_eval, Opm::abs(water_velocity));
return molecular_weight;
} else {
OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
"while injecting polymer molecular weight is requested for well " << name(), deferred_logger);
}
const int table_id = well_ecl_.getPolymerProperties().m_plymwinjtable;
const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
EvalWell molecular_weight(numWellEq_ + numEq, 0.);
if (wpolymer() == 0.) { // not injecting polymer
return molecular_weight;
}
molecular_weight = table_func.eval(throughput_eval, Opm::abs(water_velocity));
return molecular_weight;
}
@ -3154,12 +3170,14 @@ namespace Opm
StandardWell<TypeTag>::
updateWaterThroughput(const double dt, WellState &well_state) const
{
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double perf_water_vel = primary_variables_[Bhp + 1 + perf];
// we do not consider the formation damage due to water flowing from reservoir into wellbore
if (perf_water_vel > 0.) {
well_state.perfThroughput()[first_perf_ + perf] += perf_water_vel * dt;
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double perf_water_vel = primary_variables_[Bhp + 1 + perf];
// we do not consider the formation damage due to water flowing from reservoir into wellbore
if (perf_water_vel > 0.) {
well_state.perfThroughput()[first_perf_ + perf] += perf_water_vel * dt;
}
}
}
}
@ -3336,35 +3354,37 @@ namespace Opm
// if different types of extra equations are involved, this function needs to be refactored further
// checking the convergence of the extra equations related to polymer injectivity
if (this->has_polymermw && this->isInjector()) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const int dummy_component = -1;
const double maxResidualAllowed = param_.max_residual_allowed_;
using CR = ConvergenceReport;
const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double wat_vel_residual = res[Bhp + 1 + perf];
if (std::isnan(wat_vel_residual)) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, name()});
} else if (wat_vel_residual > maxResidualAllowed * 10.) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::TooLarge, dummy_component, name()});
} else if (wat_vel_residual > wat_vel_tol) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::Normal, dummy_component, name()});
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const int dummy_component = -1;
const double maxResidualAllowed = param_.max_residual_allowed_;
using CR = ConvergenceReport;
const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double wat_vel_residual = res[Bhp + 1 + perf];
if (std::isnan(wat_vel_residual)) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, name()});
} else if (wat_vel_residual > maxResidualAllowed * 10.) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::TooLarge, dummy_component, name()});
} else if (wat_vel_residual > wat_vel_tol) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::Normal, dummy_component, name()});
}
}
}
// checking the convergence of the skin pressure
const double pskin_tol = 1000.; // 1000 pascal
const auto pskin_failure_type = CR::WellFailure::Type::Pressure;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double pskin_residual = res[Bhp + 1 + perf + number_of_perforations_];
if (std::isnan(pskin_residual)) {
report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, name()});
} else if (pskin_residual > maxResidualAllowed * 10.) {
report.setWellFailed({pskin_failure_type, CR::Severity::TooLarge, dummy_component, name()});
} else if (pskin_residual > pskin_tol) {
report.setWellFailed({pskin_failure_type, CR::Severity::Normal, dummy_component, name()});
// checking the convergence of the skin pressure
const double pskin_tol = 1000.; // 1000 pascal
const auto pskin_failure_type = CR::WellFailure::Type::Pressure;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double pskin_residual = res[Bhp + 1 + perf + number_of_perforations_];
if (std::isnan(pskin_residual)) {
report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, name()});
} else if (pskin_residual > maxResidualAllowed * 10.) {
report.setWellFailed({pskin_failure_type, CR::Severity::TooLarge, dummy_component, name()});
} else if (pskin_residual > pskin_tol) {
report.setWellFailed({pskin_failure_type, CR::Severity::Normal, dummy_component, name()});
}
}
}
}

2
opm/simulators/wells/WellInterface.hpp
View File

@ -111,7 +111,7 @@ namespace Opm
static const bool has_energy = getPropValue<TypeTag, Properties::EnableEnergy>();
static const bool has_temperature = getPropValue<TypeTag, Properties::EnableTemperature>();
// flag for polymer molecular weight related
static const bool has_polymermw = getPropValue<TypeTag, Properties::EnablePolymerMW>();
static constexpr bool has_polymermw = getPropValue<TypeTag, Properties::EnablePolymerMW>();
static constexpr bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>();
static constexpr bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>();
static const int contiSolventEqIdx = Indices::contiSolventEqIdx;