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

Merge pull request #1358 from GitPaean/cleaningup_well_models

Browse Source 
Cleaningup well models
This commit is contained in:
Atgeirr Flø Rasmussen 2017-12-01 10:51:30 +01:00 committed by GitHub
commit 95c8a0e673
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 142 additions and 176 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
opm/autodiff/BlackoilWellModel.hpp
View File

@ -217,7 +217,6 @@ namespace Opm {
std::vector<double>& 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<double>& well_potentials);
const std::vector<double>& wellPerfEfficiencyFactors() const;

20
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]);
}
}
@ -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_);
}
}
}

8
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,
@ -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<double>& 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;

60
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)
@ -251,11 +252,11 @@ namespace Opm
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
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<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_;

21
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,
@ -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<double>& 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<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
@ -270,7 +267,7 @@ namespace Opm
void computeConnectionPressureDelta();
void computeWellConnectionDensitesPressures(const WellState& xw,
void computeWellConnectionDensitesPressures(const WellState& well_state,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& 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,

147
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_;
@ -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<typename TypeTag>
void
StandardWell<TypeTag>::
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<double> 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<TypeTag>::
computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw,
const WellState& well_state,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
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
@ -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<TypeTag>::
getWellConvergence(const std::vector<double>& 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<Scalar> res(numComp);
for (int comp = 0; comp < numComp; ++comp) {
// using num_components_ here for flow_ebos running 2p case.
std::vector<double> 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<double> 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<typename TypeTag>
void
StandardWell<TypeTag>::
computeWellConnectionDensitesPressures(const WellState& xw,
computeWellConnectionDensitesPressures(const WellState& well_state,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& 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<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] = 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<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 +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<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 +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<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];

17
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() {}
@ -192,16 +193,15 @@ namespace Opm
const WellState& well_state,
std::vector<double>& 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<bool>& 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;

27
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>::
@ -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);
}
}

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) );
}
}