OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-28 17:04:23 -06:00
Code Issues Packages Projects Releases Wiki Activity

removing numComponents() from WellInterface

Browse Source 
which is dumplicated from BlackoilWellModel.
This commit is contained in:
Kai Bao 2017-11-30 16:31:48 +01:00
parent 277a7809ac
commit ea3cbd1fe8
8 changed files with 100 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@ -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<TypeTag>(well_ecl, time_step, wells(), param_, *rateConverter_, pvtreg ) );
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells(),
param_, *rateConverter_, pvtreg, numComponents() ) );
} else {
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells(), param_, *rateConverter_, pvtreg) );
well_container.emplace_back(new MultisegmentWell<TypeTag>(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]);
}
}

5
opm/autodiff/MultisegmentWell.hpp
View File

@ -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<bool>* active_arg,
@ -184,12 +185,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;

56
opm/autodiff/MultisegmentWell_impl.hpp
View File

@ -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<double>(numComponents(), 0.0))
, segment_comp_initial_(numberOfSegments(), std::vector<double>(num_components_, 0.0))
, segment_densities_(numberOfSegments(), 0.0)
, segment_viscosities_(numberOfSegments(), 0.0)
, segment_mass_rates_(numberOfSegments(), 0.0)
@ -411,8 +412,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
getWellConvergence(const std::vector<double>& 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<std::vector<double>> abs_residual(numberOfSegments(), std::vector<double>(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<EvalWell>& cq_s) const
{
const int num_comp = numComponents();
std::vector<EvalWell> cmix_s(num_comp, 0.0);
std::vector<EvalWell> 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<EvalWell> b_perfcells(num_comp, 0.0);
std::vector<EvalWell> 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<EvalWell> mix_s(num_comp, 0.0);
for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
std::vector<EvalWell> 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<EvalWell> b(num_comp, 0.0);
std::vector<EvalWell> b(num_components_, 0.0);
// it is the phase viscosities asked for
std::vector<EvalWell> 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<EvalWell> mob(num_comp, 0.0);
std::vector<EvalWell> mob(num_components_, 0.0);
getMobility(ebosSimulator, perf, mob);
std::vector<EvalWell> cq_s(num_comp, 0.0);
std::vector<EvalWell> 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_;

12
opm/autodiff/StandardWell.hpp
View File

@ -106,17 +106,15 @@ namespace Opm
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> 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<bool>* active_arg,
@ -169,7 +167,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 +190,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_;

86
opm/autodiff/StandardWell_impl.hpp
View File

@ -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<TypeTag>::
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<EvalWell> cmix_s(numComp,0.0);
for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
std::vector<EvalWell> 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<EvalWell> b_perfcells_dense(numComp, 0.0);
std::vector<EvalWell> 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<EvalWell> cq_s(numComp,0.0);
std::vector<EvalWell> mob(numComp, 0.0);
std::vector<EvalWell> cq_s(num_components_,0.0);
std::vector<EvalWell> 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_;
@ -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) {
@ -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();
@ -1093,11 +1091,9 @@ namespace Opm
std::vector<double>& 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
@ -1119,12 +1115,12 @@ namespace Opm
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]) {
@ -1151,7 +1147,7 @@ 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);
@ -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
@ -1337,27 +1333,26 @@ namespace Opm
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<Scalar> res(numComp);
for (int comp = 0; comp < numComp; ++comp) {
// using num_components_ here for flow_ebos running 2p case.
std::vector<Scalar> 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);
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 +1378,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);
}
@ -1409,16 +1404,15 @@ namespace Opm
{
// Compute densities
const int nperf = number_of_perforations_;
const int numComponent = numComponents();
const int np = number_of_phases_;
std::vector<double> 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] = xw.perfPhaseRates()[(first_perf_ + perf) * np + phase];
}
if(has_solvent) {
perfRates[perf*numComponent + contiSolventEqIdx] = xw.perfRateSolvent()[first_perf_ + perf];
perfRates[perf * num_components_ + contiSolventEqIdx] = xw.perfRateSolvent()[first_perf_ + perf];
}
}
@ -1579,7 +1573,6 @@ namespace Opm
std::vector<double>& 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 +1581,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<EvalWell> cq_s(numComp, 0.0);
std::vector<EvalWell> mob(numComp, 0.0);
std::vector<EvalWell> cq_s(num_components_, 0.0);
std::vector<EvalWell> 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 +1966,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<EvalWell> cq_s(num_comp,0.0);
std::vector<EvalWell> 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];

10
opm/autodiff/WellInterface.hpp
View File

@ -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() {}
@ -282,6 +283,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<bool>& active() const;
const PhaseUsage& phaseUsage() const;
@ -290,9 +293,6 @@ namespace Opm
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
// TODO: it is dumplicated with BlackoilWellModel
int numComponents() const;
double wsolvent() const;
double wpolymer() const;

25
opm/autodiff/WellInterface_impl.hpp
View File

@ -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<typename TypeTag>
int
WellInterface<TypeTag>::
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<typename TypeTag>
double
WellInterface<TypeTag>::

17
tests/test_wellmodel.cpp
View File

@ -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<int>(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<FluidSystem, std::vector<int> >;
// 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<RateConverterType> rateConverter;
// Compute reservoir volumes for RESV controls.
@ -185,8 +189,9 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
std::vector<int>(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) );
}
}