diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index cc2176184..a9a034f1f 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -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& 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:: 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()}); + } } } } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 65329660f..605413d46 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -111,7 +111,7 @@ namespace Opm static const bool has_energy = getPropValue(); static const bool has_temperature = getPropValue(); // flag for polymer molecular weight related - static const bool has_polymermw = getPropValue(); + static constexpr bool has_polymermw = getPropValue(); static constexpr bool has_foam = getPropValue(); static constexpr bool has_brine = getPropValue(); static const int contiSolventEqIdx = Indices::contiSolventEqIdx;