From cb928f90f05ec1e57c16d7686f6a52d2fff4d172 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 3 Nov 2020 11:11:00 +0100 Subject: [PATCH] Refactor per-perforation code into separate function. Also make some methods const. --- opm/simulators/wells/StandardWell.hpp | 32 +- opm/simulators/wells/StandardWell_impl.hpp | 412 ++++++++++++-------- opm/simulators/wells/WellInterface.hpp | 4 +- opm/simulators/wells/WellInterface_impl.hpp | 2 +- 4 files changed, 268 insertions(+), 182 deletions(-) diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 73383ebd2..c2f7ba32a 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -68,6 +68,7 @@ namespace Opm using typename Base::RateConverterType; using typename Base::SparseMatrixAdapter; using typename Base::FluidState; + using typename Base::RateVector; using GasLiftHandler = Opm::GasLiftRuntime; using Base::numEq; @@ -511,6 +512,19 @@ namespace Opm WellState& well_state, Opm::DeferredLogger& deferred_logger) override; + void assembleWellEqWithoutIterationImpl(const Simulator& ebosSimulator, + const double dt, + WellState& well_state, + Opm::DeferredLogger& deferred_logger); + + void calculateSinglePerf(const Simulator& ebosSimulator, + const int perf, + WellState& well_state, + std::vector& connectionRates, + std::vector& cq_s, + EvalWell& water_flux_s, + Opm::DeferredLogger& deferred_logger) const; + // check whether the well is operable under the current reservoir condition // mostly related to BHP limit and THP limit void updateWellOperability(const Simulator& ebos_simulator, @@ -576,12 +590,17 @@ namespace Opm const EvalWell& water_velocity, Opm::DeferredLogger& deferred_logger) const; + // modify the water rate for polymer injectivity study + void handleInjectivityRate(const Simulator& ebosSimulator, + const int perf, + std::vector& cq_s) const; + // handle the extra equations for polymer injectivity study - void handleInjectivityRateAndEquations(const IntensiveQuantities& int_quants, - const WellState& well_state, - const int perf, - std::vector& cq_s, - Opm::DeferredLogger& deferred_logger); + void handleInjectivityEquations(const Simulator& ebosSimulator, + const WellState& well_state, + const int perf, + const EvalWell& water_flux_s, + Opm::DeferredLogger& deferred_logger); virtual void updateWaterThroughput(const double dt, WellState& well_state) const override; @@ -599,7 +618,8 @@ namespace Opm const IntensiveQuantities& int_quants, const WellState& well_state, const int perf, - DeferredLogger& deferred_logger); + std::vector& connectionRates, + DeferredLogger& deferred_logger) const; std::optional computeBhpAtThpLimitProd(const WellState& well_state, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 7a63dc823..d2fb83305 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -524,6 +524,8 @@ namespace Opm } + + template void StandardWell:: @@ -548,7 +550,6 @@ namespace Opm - template void StandardWell:: @@ -572,13 +573,24 @@ namespace Opm invDuneD_ = 0.0; resWell_ = 0.0; + assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, deferred_logger); + } + + + + + template + void + StandardWell:: + assembleWellEqWithoutIterationImpl(const Simulator& ebosSimulator, + const double dt, + WellState& well_state, + Opm::DeferredLogger& deferred_logger) + { + // TODO: it probably can be static member for StandardWell const double volume = 0.002831684659200; // 0.1 cu ft; - const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); - - const EvalWell& bhp = getBhp(); - // the solution gas rate and solution oil rate needs to be reset to be zero for well_state. well_state.wellVaporizedOilRates()[index_of_well_] = 0.; well_state.wellDissolvedGasRates()[index_of_well_] = 0.; @@ -588,41 +600,23 @@ namespace Opm well_state.productivityIndex()[np*index_of_well_ + p] = 0.; } + std::vector connectionRates = connectionRates_; // Copy to get right size. for (int perf = 0; perf < number_of_perforations_; ++perf) { + // Calculate perforation quantities. + std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.0}); + EvalWell water_flux_s{numWellEq_ + numEq, 0.0}; + calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, deferred_logger); - const int cell_idx = well_cells_[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - std::vector mob(num_components_, {numWellEq_ + numEq, 0.}); - getMobility(ebosSimulator, perf, mob, deferred_logger); - - std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.}); - double perf_dis_gas_rate = 0.; - double perf_vap_oil_rate = 0.; - double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); - const double Tw = well_index_[perf] * trans_mult; - computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf, - cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); - - // better way to do here is that use the cq_s and then replace the cq_s_water here? + // Equation assembly for this perforation. if (has_polymer && this->has_polymermw && this->isInjector()) { - handleInjectivityRateAndEquations(intQuants, well_state, perf, cq_s, deferred_logger); + handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger); } - - // updating the solution gas rate and solution oil rate - if (this->isProducer()) { - well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate; - well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate; - } - - if (has_energy) { - connectionRates_[perf][contiEnergyEqIdx] = 0.0; - } - + const int cell_idx = well_cells_[perf]; for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { // the cq_s entering mass balance equations need to consider the efficiency factors. const EvalWell cq_s_effective = cq_s[componentIdx] * well_efficiency_factor_; - connectionRates_[perf][componentIdx] = Base::restrictEval(cq_s_effective); + connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective); // subtract sum of phase fluxes in the well equations. resWell_[0][componentIdx] += cq_s_effective.value(); @@ -645,135 +639,9 @@ namespace Opm well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value(); } } - if (has_energy) { - - auto fs = intQuants.fluidState(); - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - // convert to reservoar conditions - EvalWell cq_r_thermal(numWellEq_ + numEq, 0.); - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - - if(FluidSystem::waterPhaseIdx == phaseIdx) - cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx)); - - // remove dissolved gas and vapporized oil - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - // q_os = q_or * b_o + rv * q_gr * b_g - // q_gs = q_gr * g_g + rs * q_or * b_o - // d = 1.0 - rs * rv - const EvalWell d = extendEval(1.0 - fs.Rv() * fs.Rs()); - // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) - if(FluidSystem::gasPhaseIdx == phaseIdx) - cq_r_thermal = (cq_s[gasCompIdx] - extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) ); - // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) - if(FluidSystem::oilPhaseIdx == phaseIdx) - cq_r_thermal = (cq_s[oilCompIdx] - extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) ); - - } else { - cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx)); - } - - // change temperature for injecting fluids - if (this->isInjector() && cq_s[activeCompIdx] > 0.0){ - // only handles single phase injection now - assert(this->well_ecl_.injectorType() != InjectorType::MULTI); - fs.setTemperature(this->well_ecl_.temperature()); - typedef typename std::decay::type::Scalar FsScalar; - typename FluidSystem::template ParameterCache paramCache; - const unsigned pvtRegionIdx = intQuants.pvtRegionIndex(); - paramCache.setRegionIndex(pvtRegionIdx); - paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx)); - paramCache.updatePhase(fs, phaseIdx); - - const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx); - fs.setDensity(phaseIdx, rho); - const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx); - fs.setEnthalpy(phaseIdx, h); - } - // compute the thermal flux - cq_r_thermal *= extendEval(fs.enthalpy(phaseIdx)) * extendEval(fs.density(phaseIdx)); - connectionRates_[perf][contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal); - } - } - - if (has_polymer) { - // TODO: the application of well efficiency factor has not been tested with an example yet - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - EvalWell cq_s_poly = cq_s[waterCompIdx]; - if (this->isInjector()) { - cq_s_poly *= wpolymer(); - } else { - cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); - } - // Note. Efficiency factor is handled in the output layer - well_state.perfRatePolymer()[first_perf_ + perf] = cq_s_poly.value(); - - cq_s_poly *= well_efficiency_factor_; - connectionRates_[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly); - - if (this->has_polymermw) { - updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, deferred_logger); - } - } - - if (has_foam) { - // TODO: the application of well efficiency factor has not been tested with an example yet - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - EvalWell cq_s_foam = cq_s[gasCompIdx] * well_efficiency_factor_; - if (this->isInjector()) { - cq_s_foam *= wfoam(); - } else { - cq_s_foam *= extendEval(intQuants.foamConcentration()); - } - connectionRates_[perf][contiFoamEqIdx] = Base::restrictEval(cq_s_foam); - } - - if (has_brine) { - // TODO: the application of well efficiency factor has not been tested with an example yet - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - EvalWell cq_s_sm = cq_s[waterCompIdx]; - if (this->isInjector()) { - cq_s_sm *= wsalt(); - } else { - cq_s_sm *= extendEval(intQuants.fluidState().saltConcentration()); - } - // Note. Efficiency factor is handled in the output layer - well_state.perfRateBrine()[first_perf_ + perf] = cq_s_sm.value(); - - cq_s_sm *= well_efficiency_factor_; - connectionRates_[perf][contiBrineEqIdx] = Base::restrictEval(cq_s_sm); - } - - // Store the perforation pressure for later usage. - well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf]; - - // Compute Productivity index if asked for - const auto& pu = phaseUsage(); - const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig(); - const Opm::Schedule& schedule = ebosSimulator.vanguard().schedule(); - for (int p = 0; p < np; ++p) { - if ( (pu.phase_pos[Water] == p && (summaryConfig.hasSummaryKey("WPIW:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name()))) - || (pu.phase_pos[Oil] == p && (summaryConfig.hasSummaryKey("WPIO:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name()))) - || (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) { - - const unsigned int compIdx = flowPhaseToEbosCompIdx(p); - const auto& fs = intQuants.fluidState(); - Eval perf_pressure = getPerfCellPressure(fs); - const double drawdown = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value(); - const bool new_well = schedule.hasWellGroupEvent(name(), ScheduleEvents::NEW_WELL, current_step_); - double productivity_index = cq_s[compIdx].value() / drawdown; - scaleProductivityIndex(perf, productivity_index, new_well, deferred_logger); - well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index; - } - } - } + // Update the connection + connectionRates_ = connectionRates; // add vol * dF/dt + Q to the well equations; for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) { @@ -809,6 +677,185 @@ namespace Opm + template + void + StandardWell:: + calculateSinglePerf(const Simulator& ebosSimulator, + const int perf, + WellState& well_state, + std::vector& connectionRates, + std::vector& cq_s, + EvalWell& water_flux_s, + Opm::DeferredLogger& deferred_logger) const + { + const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); + const int np = number_of_phases_; + const EvalWell& bhp = getBhp(); + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + std::vector mob(num_components_, {numWellEq_ + numEq, 0.}); + getMobility(ebosSimulator, perf, mob, deferred_logger); + + double perf_dis_gas_rate = 0.; + double perf_vap_oil_rate = 0.; + double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); + const double Tw = well_index_[perf] * trans_mult; + 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); + } + + // updating the solution gas rate and solution oil rate + if (this->isProducer()) { + well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate; + well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate; + } + + if (has_energy) { + connectionRates[perf][contiEnergyEqIdx] = 0.0; + } + + if (has_energy) { + + auto fs = intQuants.fluidState(); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + // convert to reservoar conditions + EvalWell cq_r_thermal(numWellEq_ + numEq, 0.); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + + if(FluidSystem::waterPhaseIdx == phaseIdx) + cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx)); + + // remove dissolved gas and vapporized oil + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + // q_os = q_or * b_o + rv * q_gr * b_g + // q_gs = q_gr * g_g + rs * q_or * b_o + // d = 1.0 - rs * rv + const EvalWell d = extendEval(1.0 - fs.Rv() * fs.Rs()); + // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) + if(FluidSystem::gasPhaseIdx == phaseIdx) + cq_r_thermal = (cq_s[gasCompIdx] - extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) ); + // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) + if(FluidSystem::oilPhaseIdx == phaseIdx) + cq_r_thermal = (cq_s[oilCompIdx] - extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) ); + + } else { + cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx)); + } + + // change temperature for injecting fluids + if (this->isInjector() && cq_s[activeCompIdx] > 0.0){ + // only handles single phase injection now + assert(this->well_ecl_.injectorType() != InjectorType::MULTI); + fs.setTemperature(this->well_ecl_.temperature()); + typedef typename std::decay::type::Scalar FsScalar; + typename FluidSystem::template ParameterCache paramCache; + const unsigned pvtRegionIdx = intQuants.pvtRegionIndex(); + paramCache.setRegionIndex(pvtRegionIdx); + paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx)); + paramCache.updatePhase(fs, phaseIdx); + + const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx); + fs.setDensity(phaseIdx, rho); + const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx); + fs.setEnthalpy(phaseIdx, h); + } + // compute the thermal flux + cq_r_thermal *= extendEval(fs.enthalpy(phaseIdx)) * extendEval(fs.density(phaseIdx)); + connectionRates[perf][contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal); + } + } + + if (has_polymer) { + // TODO: the application of well efficiency factor has not been tested with an example yet + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + EvalWell cq_s_poly = cq_s[waterCompIdx]; + if (this->isInjector()) { + cq_s_poly *= wpolymer(); + } else { + cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); + } + // Note. Efficiency factor is handled in the output layer + well_state.perfRatePolymer()[first_perf_ + perf] = cq_s_poly.value(); + + cq_s_poly *= well_efficiency_factor_; + connectionRates[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly); + + if (this->has_polymermw) { + updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger); + } + } + + if (has_foam) { + // TODO: the application of well efficiency factor has not been tested with an example yet + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + EvalWell cq_s_foam = cq_s[gasCompIdx] * well_efficiency_factor_; + if (this->isInjector()) { + cq_s_foam *= wfoam(); + } else { + cq_s_foam *= extendEval(intQuants.foamConcentration()); + } + connectionRates[perf][contiFoamEqIdx] = Base::restrictEval(cq_s_foam); + } + + if (has_brine) { + // TODO: the application of well efficiency factor has not been tested with an example yet + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + EvalWell cq_s_sm = cq_s[waterCompIdx]; + if (this->isInjector()) { + cq_s_sm *= wsalt(); + } else { + cq_s_sm *= extendEval(intQuants.fluidState().saltConcentration()); + } + // Note. Efficiency factor is handled in the output layer + well_state.perfRateBrine()[first_perf_ + perf] = cq_s_sm.value(); + + cq_s_sm *= well_efficiency_factor_; + connectionRates[perf][contiBrineEqIdx] = Base::restrictEval(cq_s_sm); + } + + // Store the perforation pressure for later usage. + well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf]; + + // Compute Productivity index if asked for + const auto& pu = phaseUsage(); + const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig(); + const Opm::Schedule& schedule = ebosSimulator.vanguard().schedule(); + for (int p = 0; p < np; ++p) { + if ( (pu.phase_pos[Water] == p && (summaryConfig.hasSummaryKey("WPIW:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name()))) + || (pu.phase_pos[Oil] == p && (summaryConfig.hasSummaryKey("WPIO:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name()))) + || (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) { + + const unsigned int compIdx = flowPhaseToEbosCompIdx(p); + const auto& fs = intQuants.fluidState(); + Eval perf_pressure = getPerfCellPressure(fs); + const double drawdown = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value(); + const bool new_well = schedule.hasWellGroupEvent(name(), ScheduleEvents::NEW_WELL, current_step_); + double productivity_index = cq_s[compIdx].value() / drawdown; + scaleProductivityIndex(perf, productivity_index, new_well, deferred_logger); + well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index; + } + } + + } + + + + template void StandardWell::assembleControlEq(const WellState& well_state, @@ -3393,14 +3440,37 @@ namespace Opm template void StandardWell:: - handleInjectivityRateAndEquations(const IntensiveQuantities& int_quants, - const WellState& well_state, - const int perf, - std::vector& cq_s, - Opm::DeferredLogger& deferred_logger) + handleInjectivityRate(const Simulator& ebosSimulator, + const int perf, + std::vector& cq_s) const { + const int cell_idx = well_cells_[perf]; + const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + const auto& fs = int_quants.fluidState(); + const EvalWell b_w = extendEval(fs.invB(FluidSystem::waterPhaseIdx)); + const double area = M_PI * bore_diameters_[perf] * perf_length_[perf]; + const int wat_vel_index = Bhp + 1 + perf; const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - const EvalWell& water_flux_s = cq_s[water_comp_idx]; + + // water rate is update to use the form from water velocity, since water velocity is + // a primary variable now + cq_s[water_comp_idx] = area * primary_variables_evaluation_[wat_vel_index] * b_w; + } + + + + + template + void + StandardWell:: + handleInjectivityEquations(const Simulator& ebosSimulator, + const WellState& well_state, + const int perf, + const EvalWell& water_flux_s, + Opm::DeferredLogger& deferred_logger) + { + const int cell_idx = well_cells_[perf]; + const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = int_quants.fluidState(); const EvalWell b_w = extendEval(fs.invB(FluidSystem::waterPhaseIdx)); const EvalWell water_flux_r = water_flux_s / b_w; @@ -3428,12 +3498,7 @@ namespace Opm invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+numEq); } - // water rate is update to use the form from water velocity, since water velocity is - // a primary variable now - cq_s[water_comp_idx] = area * primary_variables_evaluation_[wat_vel_index] * b_w; - // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B - const int cell_idx = well_cells_[perf]; for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx); } @@ -3584,7 +3649,8 @@ namespace Opm const IntensiveQuantities& int_quants, const WellState& well_state, const int perf, - DeferredLogger& deferred_logger) + std::vector& connectionRates, + DeferredLogger& deferred_logger) const { // the source term related to transport of molecular weight EvalWell cq_s_polymw = cq_s_poly; @@ -3609,7 +3675,7 @@ namespace Opm cq_s_polymw *= 0.; } } - connectionRates_[perf][this->contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw); + connectionRates[perf][this->contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw); } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index cb9024841..8d8894e0c 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -511,13 +511,13 @@ namespace Opm const std::vector& B_avg, Opm::DeferredLogger& deferred_logger); - void scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger); + void scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) const; void initCompletions(); // count the number of times an output log message is created in the productivity // index calculations - int well_productivity_index_logger_counter_; + mutable int well_productivity_index_logger_counter_; bool checkConstraints(WellState& well_state, const Schedule& schedule, diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 52dc68621..67c4641c8 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1382,7 +1382,7 @@ namespace Opm template void - WellInterface::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) + WellInterface::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) const { const auto& connection = well_ecl_.getConnections()[(*perf_data_)[perfIdx].ecl_index]; if (well_ecl_.getDrainageRadius() < 0) {