diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp index 1a136bad7..f81fb6294 100644 --- a/opm/autodiff/BlackoilWellModel.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -217,7 +217,6 @@ namespace Opm { std::vector& voidage_conversion_coeffs) const; // Calculating well potentials for each well - // TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP. void computeWellPotentials(std::vector& well_potentials); const std::vector& wellPerfEfficiencyFactors() const; diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 1f4d68f97..f63d412f0 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -226,9 +226,11 @@ namespace Opm { const int pvtreg = pvt_region_idx_[well_cell_top]; if ( !well_ecl->isMultiSegment(time_step) || !param_.use_multisegment_well_) { - well_container.emplace_back(new StandardWell(well_ecl, time_step, wells(), param_, *rateConverter_, pvtreg ) ); + well_container.emplace_back(new StandardWell(well_ecl, time_step, wells(), + param_, *rateConverter_, pvtreg, numComponents() ) ); } else { - well_container.emplace_back(new MultisegmentWell(well_ecl, time_step, wells(), param_, *rateConverter_, pvtreg) ); + well_container.emplace_back(new MultisegmentWell(well_ecl, time_step, wells(), + param_, *rateConverter_, pvtreg, numComponents() ) ); } } } @@ -402,8 +404,11 @@ namespace Opm { resetWellControlFromState() const { const int nw = numWells(); + + assert(nw == int(well_container_.size()) ); + for (int w = 0; w < nw; ++w) { - WellControls* wc = wells()->ctrls[w]; + WellControls* wc = well_container_[w]->wellControls(); well_controls_set_current( wc, well_state_.currentControls()[w]); } } @@ -699,7 +704,7 @@ namespace Opm { well_state_.currentControls()[w] = control; // TODO: for VFP control, the perf_densities are still zero here, investigate better // way to handle it later. - well_container_[w]->updateWellStateWithTarget(control, well_state_); + well_container_[w]->updateWellStateWithTarget(well_state_); // The wells are not considered to be newly added // for next time step @@ -945,12 +950,7 @@ namespace Opm { wellCollection().updateWellTargets(well_state_.wellRates()); for (int w = 0; w < numWells(); ++w) { - // TODO: check whether we need current argument in updateWellStateWithTarget - // maybe there is some circumstances that the current is different from the one - // in the WellState. - // while probalby, the current argument can be removed - const int current = well_state_.currentControls()[w]; - well_container_[w]->updateWellStateWithTarget(current, well_state_); + well_container_[w]->updateWellStateWithTarget(well_state_); } } } diff --git a/opm/autodiff/MultisegmentWell.hpp b/opm/autodiff/MultisegmentWell.hpp index a4f9f7179..e14626884 100644 --- a/opm/autodiff/MultisegmentWell.hpp +++ b/opm/autodiff/MultisegmentWell.hpp @@ -102,7 +102,8 @@ namespace Opm MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param, const RateConverterType& rate_converter, - const int pvtRegionIdx); + const int pvtRegionIdx, + const int num_components); virtual void init(const PhaseUsage* phase_usage_arg, const std::vector* active_arg, @@ -120,8 +121,7 @@ namespace Opm /// updating the well state based the control mode specified with current // TODO: later will check wheter we need current - virtual void updateWellStateWithTarget(const int current, - WellState& well_state) const; + virtual void updateWellStateWithTarget(WellState& well_state) const; /// check whether the well equations get converged for this well virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const; @@ -184,12 +184,12 @@ namespace Opm using Base::gravity_; using Base::well_controls_; using Base::perf_depth_; + using Base::num_components_; // protected functions from the Base class using Base::active; using Base::phaseUsage; using Base::name; - using Base::numComponents; using Base::flowPhaseToEbosPhaseIdx; using Base::flowPhaseToEbosCompIdx; using Base::getAllowCrossFlow; diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index 320eec4e0..d78eed4a8 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -30,14 +30,15 @@ namespace Opm MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param, const RateConverterType& rate_converter, - const int pvtRegionIdx) - : Base(well, time_step, wells, param, rate_converter, pvtRegionIdx) + const int pvtRegionIdx, + const int num_components) + : Base(well, time_step, wells, param, rate_converter, pvtRegionIdx, num_components) , segment_perforations_(numberOfSegments()) , segment_inlets_(numberOfSegments()) , cell_perforation_depth_diffs_(number_of_perforations_, 0.0) , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0) , perforation_segment_depth_diffs_(number_of_perforations_, 0.0) - , segment_comp_initial_(numberOfSegments(), std::vector(numComponents(), 0.0)) + , segment_comp_initial_(numberOfSegments(), std::vector(num_components_, 0.0)) , segment_densities_(numberOfSegments(), 0.0) , segment_viscosities_(numberOfSegments(), 0.0) , segment_mass_rates_(numberOfSegments(), 0.0) @@ -251,11 +252,11 @@ namespace Opm template void MultisegmentWell:: - updateWellStateWithTarget(const int current, - WellState& well_state) const + updateWellStateWithTarget(WellState& well_state) const { // Updating well state bas on well control // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const int current = well_state.currentControls()[index_of_well_]; const double target = well_controls_iget_target(well_controls_, current); const double* distr = well_controls_iget_distr(well_controls_, current); switch (well_controls_iget_type(well_controls_, current)) { @@ -411,8 +412,7 @@ namespace Opm MultisegmentWell:: getWellConvergence(const std::vector& B_avg) const { - // assert((int(B_avg.size()) == numComponents()) || has_polymer); - assert( (int(B_avg.size()) == numComponents()) ); + assert(int(B_avg.size()) == num_components_); // checking if any residual is NaN or too large. The two large one is only handled for the well flux std::vector> abs_residual(numberOfSegments(), std::vector(numWellEq, 0.0)); @@ -428,7 +428,7 @@ namespace Opm // TODO: the following is a little complicated, maybe can be simplified in some way? for (int seg = 0; seg < numberOfSegments(); ++seg) { for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { - if (eq_idx < numComponents()) { // phase or component mass equations + if (eq_idx < num_components_) { // phase or component mass equations const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx]; // TODO: the report can not handle the segment number yet. if (std::isnan(flux_residual)) { @@ -470,7 +470,7 @@ namespace Opm if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found // check convergence for flux residuals - for ( int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) + for ( int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { report.converged = report.converged && (maximum_residual[comp_idx] < param_.tolerance_wells_); } @@ -717,7 +717,7 @@ namespace Opm for (int seg = 0; seg < numberOfSegments(); ++seg) { // TODO: probably it should be numWellEq -1 more accurately, // while by meaning it should be num_comp - for (int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { segment_comp_initial_[seg][comp_idx] = surfaceVolumeFraction(seg, comp_idx).value(); } } @@ -925,8 +925,7 @@ namespace Opm surfaceVolumeFraction(const int seg, const int comp_idx) const { EvalWell sum_volume_fraction_scaled = 0.; - const int num_comp = numComponents(); - for (int idx = 0; idx < num_comp; ++idx) { + for (int idx = 0; idx < num_components_; ++idx) { sum_volume_fraction_scaled += volumeFractionScaled(seg, idx); } @@ -950,11 +949,10 @@ namespace Opm const bool& allow_cf, std::vector& cq_s) const { - const int num_comp = numComponents(); - std::vector cmix_s(num_comp, 0.0); + std::vector cmix_s(num_components_, 0.0); // the composition of the components inside wellbore - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); } @@ -965,7 +963,7 @@ namespace Opm const EvalWell rv = extendEval(fs.Rv()); // not using number_of_phases_ because of solvent - std::vector b_perfcells(num_comp, 0.0); + std::vector b_perfcells(num_components_, 0.0); for (int phase = 0; phase < number_of_phases_; ++phase) { const int phase_idx_ebos = flowPhaseToEbosPhaseIdx(phase); @@ -991,7 +989,7 @@ namespace Opm } // compute component volumetric rates at standard conditions - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { const EvalWell cq_p = - well_index_[perf] * (mob_perfcells[comp_idx] * drawdown); cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; } @@ -1012,7 +1010,7 @@ namespace Opm // for injecting perforations, we use total mobility EvalWell total_mob = mob_perfcells[0]; - for (int comp_idx = 1; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 1; comp_idx < num_components_; ++comp_idx) { total_mob += mob_perfcells[comp_idx]; } @@ -1057,7 +1055,7 @@ namespace Opm } // injecting connections total volumerates at standard conditions EvalWell cqt_is = cqt_i / volume_ratio; - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; } } // end for injection perforations @@ -1119,16 +1117,15 @@ namespace Opm surf_dens[phase] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx(phase), pvt_region_index ); } - const int num_comp = numComponents(); const Opm::PhaseUsage& pu = phaseUsage(); for (int seg = 0; seg < numberOfSegments(); ++seg) { // the compostion of the components inside wellbore under surface condition - std::vector mix_s(num_comp, 0.0); - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + std::vector mix_s(num_components_, 0.0); + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); } - std::vector b(num_comp, 0.0); + std::vector b(num_components_, 0.0); // it is the phase viscosities asked for std::vector visc(number_of_phases_, 0.0); const EvalWell seg_pressure = getSegmentPressure(seg); @@ -1221,7 +1218,7 @@ namespace Opm } EvalWell volrat(0.0); - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { volrat += mix[comp_idx] / b[comp_idx]; } @@ -1233,7 +1230,7 @@ namespace Opm } EvalWell density(0.0); - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { density += surf_dens[comp_idx] * mix_s[comp_idx]; } segment_densities_[seg] = density / volrat; @@ -1298,7 +1295,7 @@ namespace Opm // TODO: most of this function, if not the whole function, can be moved to the base class const int np = number_of_phases_; const int cell_idx = well_cells_[perf]; - assert (int(mob.size()) == numComponents()); + assert (int(mob.size()) == num_components_); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); @@ -1769,7 +1766,6 @@ namespace Opm const bool allow_cf = getAllowCrossFlow(); const int nseg = numberOfSegments(); - const int num_comp = numComponents(); for (int seg = 0; seg < nseg; ++seg) { // calculating the accumulation term // TODO: without considering the efficiencty factor for now @@ -1777,7 +1773,7 @@ namespace Opm { const double volume = segmentSet()[seg].volume(); // for each component - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx]) + getSegmentRate(seg, comp_idx); @@ -1791,7 +1787,7 @@ namespace Opm // considering the contributions from the inlet segments { for (const int inlet : segment_inlets_[seg]) { - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx); resWell_[seg][comp_idx] -= inlet_rate.value(); for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { @@ -1806,12 +1802,12 @@ namespace Opm for (const int perf : segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - std::vector mob(num_comp, 0.0); + std::vector mob(num_components_, 0.0); getMobility(ebosSimulator, perf, mob); - std::vector cq_s(num_comp, 0.0); + std::vector cq_s(num_components_, 0.0); computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s); - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { // the cq_s entering mass balance equations need to consider the efficiency factors. const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_; diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 68df9dbb3..617346fbc 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -106,17 +106,15 @@ namespace Opm typedef DenseAd::Evaluation EvalWell; - // TODO: should these go to WellInterface? - static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; - static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; - static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; - static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx; + using Base::contiSolventEqIdx; + using Base::contiPolymerEqIdx; StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param, const RateConverterType& rate_converter, - const int pvtRegionIdx); + const int pvtRegionIdx, + const int num_components); virtual void init(const PhaseUsage* phase_usage_arg, const std::vector* active_arg, @@ -134,8 +132,7 @@ namespace Opm /// updating the well state based the control mode specified with current // TODO: later will check wheter we need current - virtual void updateWellStateWithTarget(const int current, - WellState& xw) const; + virtual void updateWellStateWithTarget(WellState& well_state) const; /// check whether the well equations get converged for this well virtual ConvergenceReport getWellConvergence(const std::vector& B_avg) const; @@ -169,7 +166,6 @@ namespace Opm using Base::active; using Base::flowPhaseToEbosPhaseIdx; using Base::flowPhaseToEbosCompIdx; - using Base::numComponents; using Base::wsolvent; using Base::wpolymer; using Base::wellHasTHPConstraints; @@ -193,6 +189,7 @@ namespace Opm using Base::index_of_well_; using Base::well_controls_; using Base::well_type_; + using Base::num_components_; using Base::perf_rep_radius_; using Base::perf_length_; @@ -254,7 +251,7 @@ namespace Opm // calculate the properties for the well connections // to calulate the pressure difference between well connections. void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw, + const WellState& well_state, std::vector& b_perf, std::vector& rsmax_perf, std::vector& rvmax_perf, @@ -270,7 +267,7 @@ namespace Opm void computeConnectionPressureDelta(); - void computeWellConnectionDensitesPressures(const WellState& xw, + void computeWellConnectionDensitesPressures(const WellState& well_state, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, @@ -280,7 +277,7 @@ namespace Opm void computeAccumWell(); void computeWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw); + const WellState& well_state); // TODO: to check whether all the paramters are required void computePerfRate(const IntensiveQuantities& intQuants, diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 7981176b6..00b5d1be3 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -27,8 +27,9 @@ namespace Opm StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param, const RateConverterType& rate_converter, - const int pvtRegionIdx) - : Base(well, time_step, wells, param, rate_converter, pvtRegionIdx) + const int pvtRegionIdx, + const int num_components) + : Base(well, time_step, wells, param, rate_converter, pvtRegionIdx, num_components) , perf_densities_(number_of_perforations_) , perf_pressure_diffs_(number_of_perforations_) , primary_variables_(numWellEq, 0.0) @@ -105,10 +106,10 @@ namespace Opm void StandardWell:: initPrimaryVariablesEvaluation() const { - // TODO: using numComp here is only to make the 2p + dummy phase work + // TODO: using num_components_ here is only to make the 2p + dummy phase work // TODO: in theory, we should use numWellEq here. // for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) { - for (int eqIdx = 0; eqIdx < numComponents(); ++eqIdx) { + for (int eqIdx = 0; eqIdx < num_components_; ++eqIdx) { assert( (size_t)eqIdx < primary_variables_.size() ); primary_variables_evaluation_[eqIdx] = 0.0; @@ -167,7 +168,7 @@ namespace Opm const int np = number_of_phases_; const double target_rate = well_controls_get_current_target(wc); - assert(comp_idx < numComponents()); + assert(comp_idx < num_components_); const auto pu = phaseUsage(); // TODO: the formulation for the injectors decides it only work with single phase @@ -357,8 +358,7 @@ namespace Opm wellSurfaceVolumeFraction(const int compIdx) const { EvalWell sum_volume_fraction_scaled = 0.; - const int numComp = numComponents(); - for (int idx = 0; idx < numComp; ++idx) { + for (int idx = 0; idx < num_components_; ++idx) { sum_volume_fraction_scaled += wellVolumeFractionScaled(idx); } @@ -398,9 +398,8 @@ namespace Opm { const Opm::PhaseUsage& pu = phaseUsage(); const int np = number_of_phases_; - const int numComp = numComponents(); - std::vector cmix_s(numComp,0.0); - for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + std::vector cmix_s(num_components_,0.0); + for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx); } const auto& fs = intQuants.fluidState(); @@ -408,7 +407,7 @@ namespace Opm const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); const EvalWell rs = extendEval(fs.Rs()); const EvalWell rv = extendEval(fs.Rv()); - std::vector b_perfcells_dense(numComp, 0.0); + std::vector b_perfcells_dense(num_components_, 0.0); for (int phase = 0; phase < np; ++phase) { const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); @@ -429,7 +428,7 @@ namespace Opm } // compute component volumetric rates at standard conditions - for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { const EvalWell cq_p = - Tw * (mob_perfcells_dense[componentIdx] * drawdown); cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p; } @@ -451,7 +450,7 @@ namespace Opm // Using total mobilities EvalWell total_mob_dense = mob_perfcells_dense[0]; - for (int componentIdx = 1; componentIdx < numComp; ++componentIdx) { + for (int componentIdx = 1; componentIdx < num_components_; ++componentIdx) { total_mob_dense += mob_perfcells_dense[componentIdx]; } @@ -503,7 +502,7 @@ namespace Opm // injecting connections total volumerates at standard conditions EvalWell cqt_is = cqt_i/volumeRatio; //std::cout << "volrat " << volumeRatio << " " << volrat_perf_[perf] << std::endl; - for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is; // * b_perfcells_dense[phase]; } } @@ -521,7 +520,6 @@ namespace Opm WellState& well_state, bool only_wells) { - const int numComp = numComponents(); const int np = number_of_phases_; // clear all entries @@ -546,12 +544,12 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - std::vector cq_s(numComp,0.0); - std::vector mob(numComp, 0.0); + std::vector cq_s(num_components_,0.0); + std::vector mob(num_components_, 0.0); getMobility(ebosSimulator, perf, mob); computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); - for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + 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_; @@ -590,14 +588,14 @@ namespace Opm } if (has_polymer) { - EvalWell cq_s_poly = cq_s[Water]; + // TODO: the application of well efficiency factor has not been tested with an example yet + EvalWell cq_s_poly = cq_s[Water] * well_efficiency_factor_; if (well_type_ == INJECTOR) { cq_s_poly *= wpolymer(); } else { cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); } if (!only_wells) { - // TODO: we need to consider the efficiency here. for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx); } @@ -610,7 +608,7 @@ namespace Opm } // add vol * dF/dt + Q to the well equations; - for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { + for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt; resWell_loc += getQs(componentIdx) * well_efficiency_factor_; for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { @@ -644,12 +642,12 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& fs = intQuants.fluidState(); - EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); - EvalWell bhp = getBhp(); + const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + const EvalWell bhp = getBhp(); // Pressure drawdown (also used to determine direction of flow) - EvalWell well_pressure = bhp + perf_pressure_diffs_[perf]; - EvalWell drawdown = pressure - well_pressure; + const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf]; + const EvalWell drawdown = pressure - well_pressure; if (drawdown.value() < 0 && well_type_ == INJECTOR) { return false; @@ -675,7 +673,7 @@ namespace Opm { const int np = number_of_phases_; const int cell_idx = well_cells_[perf]; - assert (int(mob.size()) == numComponents()); + assert (int(mob.size()) == num_components_); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); @@ -978,40 +976,40 @@ namespace Opm template void StandardWell:: - updateWellStateWithTarget(const int current, - WellState& xw) const + updateWellStateWithTarget(WellState& well_state) const { // number of phases const int np = number_of_phases_; const int well_index = index_of_well_; const WellControls* wc = well_controls_; + const int current = well_state.currentControls()[well_index]; // Updating well state and primary variables. // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE const double target = well_controls_iget_target(wc, current); const double* distr = well_controls_iget_distr(wc, current); switch (well_controls_iget_type(wc, current)) { case BHP: - xw.bhp()[well_index] = target; + well_state.bhp()[well_index] = target; // TODO: similar to the way below to handle THP // we should not something related to thp here when there is thp constraint break; case THP: { - xw.thp()[well_index] = target; + well_state.thp()[well_index] = target; const Opm::PhaseUsage& pu = phaseUsage(); std::vector rates(3, 0.0); if (active()[ Water ]) { - rates[ Water ] = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ]; + rates[ Water ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Water ] ]; } if (active()[ Oil ]) { - rates[ Oil ] = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ]; + rates[ Oil ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Oil ] ]; } if (active()[ Gas ]) { - rates[ Gas ] = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ]; + rates[ Gas ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Gas ] ]; } - xw.bhp()[well_index] = calculateBhpFromThp(rates, current); + well_state.bhp()[well_index] = calculateBhpFromThp(rates, current); break; } @@ -1034,9 +1032,9 @@ namespace Opm for (int phase = 0; phase < np; ++phase) { if (distr[phase] > 0.) { - xw.wellRates()[np*well_index + phase] = target / distr[phase]; + well_state.wellRates()[np*well_index + phase] = target / distr[phase]; } else { - xw.wellRates()[np * well_index + phase] = 0.; + well_state.wellRates()[np * well_index + phase] = 0.; } } } else if (well_type_ == PRODUCER) { @@ -1046,7 +1044,7 @@ namespace Opm double original_rates_under_phase_control = 0.0; for (int phase = 0; phase < np; ++phase) { if (distr[phase] > 0.0) { - original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase]; + original_rates_under_phase_control += well_state.wellRates()[np * well_index + phase] * distr[phase]; } } @@ -1054,17 +1052,17 @@ namespace Opm double scaling_factor = target / original_rates_under_phase_control; for (int phase = 0; phase < np; ++phase) { - xw.wellRates()[np * well_index + phase] *= scaling_factor; + well_state.wellRates()[np * well_index + phase] *= scaling_factor; } } else { // scaling factor is not well defied when original_rates_under_phase_control is zero // separating targets equally between phases under control const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl; for (int phase = 0; phase < np; ++phase) { if (distr[phase] > 0.0) { - xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase]; + well_state.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase]; } else { // this only happens for SURFACE_RATE control - xw.wellRates()[np * well_index + phase] = target_rate_divided; + well_state.wellRates()[np * well_index + phase] = target_rate_divided; } } } @@ -1075,7 +1073,7 @@ namespace Opm break; } // end of switch - updatePrimaryVariables(xw); + updatePrimaryVariables(well_state); } @@ -1086,18 +1084,16 @@ namespace Opm void StandardWell:: computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& xw, + const WellState& well_state, std::vector& b_perf, std::vector& rsmax_perf, std::vector& rvmax_perf, std::vector& surf_dens_perf) const { const int nperf = number_of_perforations_; - // TODO: can make this a member? - const int numComp = numComponents(); const PhaseUsage& pu = phaseUsage(); - b_perf.resize(nperf*numComp); - surf_dens_perf.resize(nperf*numComp); + b_perf.resize(nperf * num_components_); + surf_dens_perf.resize(nperf * num_components_); const int w = index_of_well_; //rs and rv are only used if both oil and gas is present @@ -1114,25 +1110,25 @@ namespace Opm // TODO: this is another place to show why WellState need to be a vector of WellState. // TODO: to check why should be perf - 1 - const double p_above = perf == 0 ? xw.bhp()[w] : xw.perfPress()[first_perf_ + perf - 1]; - const double p_avg = (xw.perfPress()[first_perf_ + perf] + p_above)/2; + const double p_above = perf == 0 ? well_state.bhp()[w] : well_state.perfPress()[first_perf_ + perf - 1]; + const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2; const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); if (pu.phase_used[Water]) { - b_perf[ pu.phase_pos[Water] + perf * numComp] = + b_perf[ pu.phase_pos[Water] + perf * num_components_] = FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); } if (pu.phase_used[Gas]) { - const int gaspos = pu.phase_pos[Gas] + perf * numComp; + const int gaspos = pu.phase_pos[Gas] + perf * num_components_; const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases; if (pu.phase_used[Oil]) { const int oilpos_well = pu.phase_pos[Oil] + w * pu.num_phases; - const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers + const double oilrate = std::abs(well_state.wellRates()[oilpos_well]); //in order to handle negative rates in producers rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); if (oilrate > 0) { - const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); + const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w); double rv = 0.0; if (gasrate > 0) { rv = oilrate / gasrate; @@ -1151,14 +1147,14 @@ namespace Opm } if (pu.phase_used[Oil]) { - const int oilpos = pu.phase_pos[Oil] + perf * numComp; + const int oilpos = pu.phase_pos[Oil] + perf * num_components_; const int oilpos_well = pu.phase_pos[Oil] + w * pu.num_phases; if (pu.phase_used[Gas]) { rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases; - const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); + const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w); if (gasrate > 0) { - const double oilrate = std::abs(xw.wellRates()[oilpos_well]); + const double oilrate = std::abs(well_state.wellRates()[oilpos_well]); double rs = 0.0; if (oilrate > 0) { rs = gasrate / oilrate; @@ -1175,13 +1171,13 @@ namespace Opm // Surface density. for (int p = 0; p < pu.num_phases; ++p) { - surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); + surf_dens_perf[num_components_ * perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); } // We use cell values for solvent injector if (has_solvent) { - b_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value(); - surf_dens_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventRefDensity(); + b_perf[num_components_ * perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value(); + surf_dens_perf[num_components_ * perf + contiSolventEqIdx] = intQuants.solventRefDensity(); } } } @@ -1202,7 +1198,7 @@ namespace Opm // Verify that we have consistent input. const int np = number_of_phases_; const int nperf = number_of_perforations_; - const int num_comp = numComponents(); + const int num_comp = num_components_; const PhaseUsage& phase_usage = phaseUsage(); // 1. Compute the flow (in surface volume units for each @@ -1333,31 +1329,27 @@ namespace Opm StandardWell:: getWellConvergence(const std::vector& B_avg) const { - typedef double Scalar; - typedef std::vector< Scalar > Vector; - const int np = number_of_phases_; - const int numComp = numComponents(); // the following implementation assume that the polymer is always after the w-o-g phases // For the polymer case, there is one more mass balance equations of reservoir than wells - assert((int(B_avg.size()) == numComp) || has_polymer); + assert((int(B_avg.size()) == num_components_) || has_polymer); const double tol_wells = param_.tolerance_wells_; const double maxResidualAllowed = param_.max_residual_allowed_; // TODO: it should be the number of numWellEq - // using numComp here for flow_ebos running 2p case. - std::vector res(numComp); - for (int comp = 0; comp < numComp; ++comp) { + // using num_components_ here for flow_ebos running 2p case. + std::vector res(num_components_); + for (int comp = 0; comp < num_components_; ++comp) { // magnitude of the residual matters res[comp] = std::abs(resWell_[0][comp]); } - Vector well_flux_residual(numComp); + std::vector well_flux_residual(num_components_); // Finish computation - for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + for ( int compIdx = 0; compIdx < num_components_; ++compIdx ) { well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx]; } @@ -1383,7 +1375,7 @@ namespace Opm if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found // check convergence - for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + for ( int compIdx = 0; compIdx < num_components_; ++compIdx ) { report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells); } @@ -1401,7 +1393,7 @@ namespace Opm template void StandardWell:: - computeWellConnectionDensitesPressures(const WellState& xw, + computeWellConnectionDensitesPressures(const WellState& well_state, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, @@ -1409,16 +1401,15 @@ namespace Opm { // Compute densities const int nperf = number_of_perforations_; - const int numComponent = numComponents(); const int np = number_of_phases_; std::vector perfRates(b_perf.size(),0.0); for (int perf = 0; perf < nperf; ++perf) { for (int phase = 0; phase < np; ++phase) { - perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[(first_perf_ + perf) * np + phase]; + perfRates[perf * num_components_ + phase] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + phase]; } if(has_solvent) { - perfRates[perf*numComponent + contiSolventEqIdx] = xw.perfRateSolvent()[first_perf_ + perf]; + perfRates[perf * num_components_ + contiSolventEqIdx] = well_state.perfRateSolvent()[first_perf_ + perf]; } } @@ -1579,7 +1570,6 @@ namespace Opm std::vector& well_flux) const { const int np = number_of_phases_; - const int numComp = numComponents(); well_flux.resize(np, 0.0); const bool allow_cf = crossFlowAllowed(ebosSimulator); @@ -1588,8 +1578,8 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); // flux for each perforation - std::vector cq_s(numComp, 0.0); - std::vector mob(numComp, 0.0); + std::vector cq_s(num_components_, 0.0); + std::vector mob(num_components_, 0.0); getMobility(ebosSimulator, perf, mob); computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); @@ -1973,10 +1963,9 @@ namespace Opm return; } // compute the well water velocity with out shear effects. - const int num_comp = numComponents(); const bool allow_cf = crossFlowAllowed(ebos_simulator); const EvalWell& bhp = getBhp(); - std::vector cq_s(num_comp,0.0); + std::vector cq_s(num_components_,0.0); computePerfRate(int_quant, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); // TODO: make area a member const double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf]; diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index 20a420a08..979fdf115 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -92,7 +92,7 @@ namespace Opm static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; - + static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; // For the conversion between the surface volume rate and resrevoir voidage rate using RateConverterType = RateConverter:: @@ -102,7 +102,8 @@ namespace Opm WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param, const RateConverterType& rate_converter, - const int pvtRegionIdx); + const int pvtRegionIdx, + const int num_components); /// Virutal destructor virtual ~WellInterface() {} @@ -192,16 +193,15 @@ namespace Opm const WellState& well_state, std::vector& well_potentials) = 0; - virtual void updateWellStateWithTarget(const int current, - WellState& xw) const = 0; + virtual void updateWellStateWithTarget(WellState& well_state) const = 0; - void updateWellControl(WellState& xw, + void updateWellControl(WellState& well_state, wellhelpers::WellSwitchingLogger& logger) const; virtual void updatePrimaryVariables(const WellState& well_state) const = 0; virtual void calculateExplicitQuantities(const Simulator& ebosSimulator, - const WellState& xw) = 0; // should be const? + const WellState& well_state) = 0; // should be const? protected: @@ -282,6 +282,8 @@ namespace Opm // We assume a well to not penetrate more than one pvt region. const int pvtRegionIdx_; + const int num_components_; + const std::vector& active() const; const PhaseUsage& phaseUsage() const; @@ -290,9 +292,6 @@ namespace Opm int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; - // TODO: it is dumplicated with BlackoilWellModel - int numComponents() const; - double wsolvent() const; double wpolymer() const; diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index 994d3d87c..bcc62dc22 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -28,12 +28,14 @@ namespace Opm WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param, const RateConverterType& rate_converter, - const int pvtRegionIdx) + const int pvtRegionIdx, + const int num_components) : well_ecl_(well) , current_step_(time_step) , param_(param) , rateConverter_(rate_converter) , pvtRegionIdx_(pvtRegionIdx) + , num_components_(num_components) { if (!well) { OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface"); @@ -269,27 +271,6 @@ namespace Opm - template - int - WellInterface:: - numComponents() const - { - // TODO: how about two phase polymer - if (number_of_phases_ == 2) { - return 2; - } - - int numComp = FluidSystem::numComponents; - - if (has_solvent) { - numComp++; - } - return numComp; - } - - - - template double WellInterface:: @@ -462,7 +443,7 @@ namespace Opm } if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) { - updateWellStateWithTarget(updated_control_index, well_state); + updateWellStateWithTarget(well_state); } } diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index 81af57adf..754812c92 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -126,7 +126,7 @@ struct SetupTest { BOOST_AUTO_TEST_CASE(TestStandardWellInput) { - SetupTest setup_test; + const SetupTest setup_test; const Wells* wells = setup_test.wells_manager->c_wells(); const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep); BOOST_CHECK_EQUAL( wells_ecl.size(), 2); @@ -145,14 +145,16 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) { std::vector(10, 0))); const int pvtIdx = 0; - BOOST_CHECK_THROW( StandardWell( well, -1, wells, param, *rateConverter, pvtIdx), std::invalid_argument); - BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param , *rateConverter, pvtIdx), std::invalid_argument); - BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param , *rateConverter, pvtIdx), std::invalid_argument); + const int num_comp = wells->number_of_phases; + + BOOST_CHECK_THROW( StandardWell( well, -1, wells, param, *rateConverter, pvtIdx, num_comp), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param , *rateConverter, pvtIdx, num_comp), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param , *rateConverter, pvtIdx, num_comp), std::invalid_argument); } BOOST_AUTO_TEST_CASE(TestBehavoir) { - SetupTest setup_test; + const SetupTest setup_test; const Wells* wells_struct = setup_test.wells_manager->c_wells(); const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep); const int current_timestep = setup_test.current_timestep; @@ -178,6 +180,8 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { using RateConverterType = Opm::RateConverter:: SurfaceToReservoirVoidage >; // Compute reservoir volumes for RESV controls. + // TODO: not sure why for this class the initlizer list does not work + // otherwise we should make a meaningful const PhaseUsage here. Opm::PhaseUsage phaseUsage; std::unique_ptr rateConverter; // Compute reservoir volumes for RESV controls. @@ -185,8 +189,9 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { std::vector(10, 0))); const int pvtIdx = 0; + const int num_comp = wells_struct->number_of_phases; - wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param, *rateConverter, pvtIdx) ); + wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param, *rateConverter, pvtIdx, num_comp) ); } }