cleaning StandardWell and StandardWellV

to make them more comparable, so that we can do refactoring easier
later.
This commit is contained in:
Kai Bao 2018-12-12 15:58:15 +01:00
parent 167306169b
commit 5535406b67
4 changed files with 120 additions and 124 deletions

View File

@ -66,7 +66,9 @@ namespace Opm
static const int numWellConservationEq = numEq - numPolymerEq - numEnergyEq;
// number of the well control equations
static const int numWellControlEq = 1;
static const int numWellEq = numWellConservationEq + numWellControlEq;
// number of the well equations that will always be used
// based on the solution strategy, there might be other well equations be introduced
static const int numStaticWellEq = numWellConservationEq + numWellControlEq;
// the positions of the primary variables for StandardWell
// the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g,
@ -85,7 +87,11 @@ namespace Opm
// the index for Bhp in primary variables and also the index of well control equation
// they both will be the last one in their respective system.
// TODO: we should have indices for the well equations and well primary variables separately
static const int Bhp = numWellEq - numWellControlEq;
static const int Bhp = numStaticWellEq - numWellControlEq;
// total number of the welll equations and primary variables
// for StandardWell, there is no extra well equations will be used
static const int numWellEq = numStaticWellEq;
using typename Base::Scalar;
@ -109,7 +115,6 @@ namespace Opm
// the matrix type for the diagonal matrix D
typedef Dune::FieldMatrix<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
// the matrix type for the non-diagonal matrix B and C^T
@ -303,11 +308,11 @@ namespace Opm
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state);
// TODO: to check whether all the paramters are required
void computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob_perfcells_dense,
const double Tw, const EvalWell& bhp, const double& cdp,
const bool& allow_cf, std::vector<EvalWell>& cq_s,
const std::vector<EvalWell>& mob,
const EvalWell& bhp,
const int perf,
const bool allow_cf, std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate, double& perf_vap_oil_rate) const;
// TODO: maybe we should provide a light version of computePerfRate, which does not include the

View File

@ -71,7 +71,9 @@ namespace Opm
static const int numWellConservationEq = numEq - numPolymerEq - numEnergyEq;
// number of the well control equations
static const int numWellControlEq = 1;
int numWellEq = numWellConservationEq + numWellControlEq;
// number of the well equations that will always be used
// based on the solution strategy, there might be other well equations be introduced
static const int numStaticWellEq = numWellConservationEq + numWellControlEq;
// the positions of the primary variables for StandardWell
// the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g,
@ -90,7 +92,11 @@ namespace Opm
// the index for Bhp in primary variables and also the index of well control equation
// they both will be the last one in their respective system.
// TODO: we should have indices for the well equations and well primary variables separately
int Bhp = numWellEq - numWellControlEq;
static const int Bhp = numStaticWellEq - numWellControlEq;
// total number of the welll equations and primary variables
// there might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq = numStaticWellEq;
using typename Base::Scalar;
@ -112,11 +118,13 @@ namespace Opm
typedef Dune::DynamicVector<Scalar> VectorBlockWellType;
typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
// the matrix type for the matrices
// since we will resize all the matrices individually, one single type for the three
// matrices will be plenty
typedef Dune::DynamicMatrix<Scalar> MatrixBlockWellType;
typedef Dune::BCRSMatrix <MatrixBlockWellType> MatWell;
// the matrix type for the diagonal matrix D
typedef Dune::DynamicMatrix<Scalar> DiagMatrixBlockWellType;
typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
// the matrix type for the non-diagonal matrix B and C^T
typedef Dune::DynamicMatrix<Scalar> OffDiagMatrixBlockWellType;
typedef Dune::BCRSMatrix<OffDiagMatrixBlockWellType> OffDiagMatWell;
typedef DenseAd::DynamicEvaluation<Scalar> EvalWell;
@ -227,10 +235,10 @@ namespace Opm
BVectorWell resWell_;
// two off-diagonal matrices
MatWell duneB_;
MatWell duneC_;
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// diagonal matrix for the well
MatWell invDuneD_;
DiagMatWell invDuneD_;
// several vector used in the matrix calculation
mutable BVectorWell Bx_;
@ -305,11 +313,11 @@ namespace Opm
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state);
// TODO: to check whether all the paramters are required
void computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob_perfcells_dense, const int perf,
const double Tw, const EvalWell& bhp, const double& cdp,
const bool& allow_cf, std::vector<EvalWell>& cq_s,
const std::vector<EvalWell>& mob,
const EvalWell& bhp,
const int perf,
const bool allow_cf, std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate, double& perf_vap_oil_rate) const;
// TODO: maybe we should provide a light version of computePerfRate, which does not include the

View File

@ -39,9 +39,9 @@ namespace Opm
{
assert(num_components_ == numWellConservationEq);
duneB_.setBuildMode( MatWell::row_wise );
duneC_.setBuildMode( MatWell::row_wise );
invDuneD_.setBuildMode( MatWell::row_wise );
duneB_.setBuildMode( OffDiagMatWell::row_wise );
duneC_.setBuildMode( OffDiagMatWell::row_wise );
invDuneD_.setBuildMode( DiagMatWell::row_wise );
}
@ -65,7 +65,7 @@ namespace Opm
}
// counting/updating primary variable numbers
if (this->has_polymermw) {
if (this->has_polymermw && well_type_ == INJECTOR) {
// adding a primary variable for water perforation rate per connection
numWellEq += number_of_perforations_;
// adding a primary variable for skin pressure per connection
@ -250,7 +250,6 @@ namespace Opm
// Oil fraction
EvalWell well_fraction(numWellEq + numEq, 1.0);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
well_fraction -= primary_variables_evaluation_[WFrac];
}
@ -308,17 +307,14 @@ namespace Opm
void
StandardWellV<TypeTag>::
computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob_perfcells_dense, const int perf,
const double Tw, const EvalWell& bhp, const double& cdp,
const bool& allow_cf, std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate, double& perf_vap_oil_rate) const
const std::vector<EvalWell>& mob,
const EvalWell& bhp,
const int perf,
const bool allow_cf,
std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate) const
{
// TODO: this should only happen once per well
// TODO: this can be moved closer where it is used
std::vector<EvalWell> cmix_s(num_components_, EvalWell{numWellEq + numEq});
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
}
const auto& fs = intQuants.fluidState();
const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
@ -338,23 +334,15 @@ namespace Opm
}
// Pressure drawdown (also used to determine direction of flow)
const EvalWell well_pressure = bhp + cdp;
const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
EvalWell drawdown = pressure - well_pressure;
// TODO: perf is introduced to handle the skin pressure only
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index];
/* if (drawdown.value() * skin_pressure.value() > 0.) {
std::cout << " the drawdown and skin_pressure has the same sign for well " << name() << " perf " << perf << std::endl;
if (this->has_polymermw && well_type_ == INJECTOR) {
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index];
drawdown += skin_pressure;
}
if (drawdown.value() > 0.) {
std::cout << " perf " << perf << " of well " << name() << " has crossflow " << std::endl;
} */
drawdown += skin_pressure;
// producing perforations
if ( drawdown.value() > 0 ) {
//Do nothing if crossflow is not allowed
@ -362,9 +350,10 @@ namespace Opm
return;
}
const double Tw = well_index_[perf];
// compute component volumetric rates at standard conditions
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
const EvalWell cq_p = - Tw * (mob_perfcells_dense[componentIdx] * drawdown);
const EvalWell cq_p = - Tw * (mob[componentIdx] * drawdown);
cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
}
@ -393,14 +382,21 @@ namespace Opm
}
// Using total mobilities
EvalWell total_mob_dense = mob_perfcells_dense[0];
EvalWell total_mob_dense = mob[0];
for (int componentIdx = 1; componentIdx < num_components_; ++componentIdx) {
total_mob_dense += mob_perfcells_dense[componentIdx];
total_mob_dense += mob[componentIdx];
}
// injection perforations total volume rates
const double Tw = well_index_[perf];
const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown);
// surface volume fraction of fluids within wellbore
std::vector<EvalWell> cmix_s(num_components_, EvalWell{numWellEq + numEq});
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
}
// compute volume ratio between connection at standard conditions
EvalWell volumeRatio(numWellEq + numEq, 0.);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
@ -517,12 +513,13 @@ namespace Opm
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> cq_s(num_components_, {numWellEq + numEq, 0.});
std::vector<EvalWell> mob(num_components_, {numWellEq + numEq, 0.});
getMobility(ebosSimulator, perf, mob);
std::vector<EvalWell> cq_s(num_components_, {numWellEq + numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(intQuants, mob, perf, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
computePerfRate(intQuants, mob, bhp, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
// better way to do here is that use the cq_s and then replace the cq_s_water here?
@ -636,9 +633,9 @@ namespace Opm
}
connectionRates_[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
if (this->has_polymermw) {
// TODO: should molecular weight be a Evalution type?
if (this->has_polymermw && well_type_ == INJECTOR) {
// the source term related to transport of molecular weight
// TODO: should molecular weight be a Evalution type?
EvalWell cq_s_polymw = cq_s_poly;
if (well_type_ == INJECTOR) {
const int wat_vel_index = Bhp + 1 + perf;
@ -646,8 +643,6 @@ namespace Opm
if (water_velocity > 0.) { // injecting
const double throughput = well_state.perfThroughput()[first_perf_ + perf];
const EvalWell molecular_weight = wpolymermw(throughput, water_velocity);
// std::cout << " well " << name() << " perf " << perf << " throughput " << throughput << " water_velocity " << water_velocity
// << " molecular_weight " << molecular_weight.value() << std::endl;
cq_s_polymw *= molecular_weight;
} else {
cq_s_polymw *= 0.;
@ -958,36 +953,21 @@ namespace Opm
}
// for the water velocity and skin pressure
// TODO: it will be some part need to come back to make more sophiscated.
{
if (this->has_polymermw && well_type_ == INJECTOR) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const int wat_vel_index = Bhp + 1 + perf;
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
// we can add a relaxation factor here
const double relaxation_factor = 0.9;
const double dx_wat_vel = dwells[0][wat_vel_index];
// const double wat_vel_current = primary_variables_[wat_vel_index];
primary_variables_[wat_vel_index] -= 0.9 * dx_wat_vel;
/*
// the following is the techniques implemented in MRST for better convergence rate, while it did not help
const double updated_water_vel = primary_variables_[wat_vel_index] - 0.9 * dx_wat_vel;
primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
if (updated_water_vel * primary_variables_[wat_vel_index] >= 0.) { // no changing sign
primary_variables_[wat_vel_index] = updated_water_vel;
} else {
if (perf_flow_switch_attemp_[perf]) {
primary_variables_[wat_vel_index] = updated_water_vel;
perf_flow_switch_attemp_[perf] = false;
} else {
primary_variables_[wat_vel_index] *= 1.e-8; // chopping it and waiting for next iteration to decide
perf_flow_switch_attemp_[perf] = true;
}
} */
const double dx_pskin = dwells[0][pskin_index];
primary_variables_[pskin_index] -= 0.9 * dx_pskin;
primary_variables_[pskin_index] -= relaxation_factor * dx_pskin;
}
}
}
@ -1145,10 +1125,12 @@ namespace Opm
updateThp(well_state);
// other primary variables related to polymer injection
for (int perf = 0; perf < number_of_perforations_; ++perf) {
well_state.perfWaterVelocity()[first_perf_ + perf] = primary_variables_[Bhp + 1 + perf];
well_state.perfSkinPressure()[first_perf_ + perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf];
// other primary variables related to polymer injectivity study
if (this->has_polymermw && well_type_ == INJECTOR) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
well_state.perfWaterVelocity()[first_perf_ + perf] = primary_variables_[Bhp + 1 + perf];
well_state.perfSkinPressure()[first_perf_ + perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf];
}
}
}
@ -2058,7 +2040,7 @@ namespace Opm
report.setWellFailed({type, CR::Severity::Normal, dummy_component, name()});
}
if (this->has_polymermw) {
if (this->has_polymermw && well_type_ == INJECTOR) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
@ -2297,12 +2279,13 @@ 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(num_components_, {numWellEq + numEq, 0.});
std::vector<EvalWell> mob(num_components_, {numWellEq + numEq, 0.});
getMobility(ebosSimulator, perf, mob);
std::vector<EvalWell> cq_s(num_components_, {numWellEq + numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(intQuants, mob, perf, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
computePerfRate(intQuants, mob, bhp, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
for(int p = 0; p < np; ++p) {
@ -2559,11 +2542,13 @@ namespace Opm
// BHP
primary_variables_[Bhp] = well_state.bhp()[index_of_well_];
// other primary variables related to polymer injection
for (int perf = 0; perf < number_of_perforations_; ++perf) {
primary_variables_[Bhp + 1 + perf] = well_state.perfWaterVelocity()[first_perf_ + perf];
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
}
// other primary variables related to polymer injectio
if (this->has_polymermw && well_type_ == INJECTOR) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
primary_variables_[Bhp + 1 + perf] = well_state.perfWaterVelocity()[first_perf_ + perf];
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
}
}
}
@ -2699,10 +2684,11 @@ namespace Opm
// TODO: do we need to turn on crossflow here?
const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
const EvalWell& bhp = getBhp();
std::vector<EvalWell> cq_s(num_components_, {numWellEq + numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(int_quant, mob, perf, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
computePerfRate(int_quant, mob, bhp, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
// TODO: make area a member
const double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf];
@ -2740,7 +2726,6 @@ namespace Opm
// B and C have 1 row, nc colums and nonzero
// at (0,j) only if this well has a perforation at cell j.
// TODO: we have some matrix type problem here.
for ( auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC )
{
const auto row_index = colC.index();
@ -2971,36 +2956,27 @@ namespace Opm
OPM_THROW(std::runtime_error, "Polymermw is not activated, "
"while injecting skin pressure is requested" << name());
}
// TODO: for the current implementation, without calculating the wellbore polymer concentration, poly_inj_conc can just be a scalar
// TODO: should we ignore the skin pressure when flowing from the reservoir into wellbore.
const double sign = water_velocity >= 0. ? 1.0 : -1.0;
const EvalWell water_velocity_abs = Opm::abs(water_velocity);
// TODO: we need to consider the direction of the water velocity here
if (poly_inj_conc == 0.) {
return sign * pskinwater(throughput, water_velocity_abs);
}
// otherwise, we need to use the skin pressure from polymer injection
// currently, we do not consider the reference concentration table. When there is polymer
// injection, we use SKPRPOLY table provided.
// TODO: we might need to do the interpolation when the polymer injection concentration
// is different from the reference concentration of the table.
const int polymer_table_id = well_ecl_->getPolymerProperties(current_step_).m_skprpolytable;
if (polymer_table_id <= 0) {
OPM_THROW(std::runtime_error, "Unavailable SKPRPOLY table id used for well " << name());
}
const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
// TODO: eventually it should be used
const double referene_concentration = skprpolytable.refConcentration;
const double reference_concentration = skprpolytable.refConcentration;
const EvalWell throughput_eval(numWellEq + numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_poly(numWellEq + numEq, 0.0);
skprpolytable.table_func.eval(throughput_eval, water_velocity_abs, pskin_poly);
if (poly_inj_conc == referene_concentration) {
if (poly_inj_conc == reference_concentration) {
return sign * pskin_poly;
}
// poly_inj_conc != reference concentration of the table, then some interpolation will be required
const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs);
const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / referene_concentration * poly_inj_conc;
const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
return sign * pskin;
}
@ -3020,9 +2996,6 @@ namespace Opm
}
const int table_id = well_ecl_->getPolymerProperties(current_step_).m_plymwinjtable;
const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
// TODO: more likely we can use extrapolation for the throughput while not for the velocity
// TODO: basically for velocity, we need to request a new table when it exceeds the velocity range
// TODO: from the table, we should give a clear message for that
const EvalWell throughput_eval(numWellEq + numEq, throughput);
EvalWell molecular_weight(numWellEq + numEq, 0.);
if (wpolymer() == 0.) { // not injecting polymer

View File

@ -279,15 +279,14 @@ namespace Opm
void
StandardWell<TypeTag>::
computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob_perfcells_dense,
const double Tw, const EvalWell& bhp, const double& cdp,
const bool& allow_cf, std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate, double& perf_vap_oil_rate) const
const std::vector<EvalWell>& mob,
const EvalWell& bhp,
const int perf,
const bool allow_cf,
std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate) const
{
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();
const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
const EvalWell rs = extendEval(fs.Rs());
@ -306,7 +305,7 @@ namespace Opm
}
// Pressure drawdown (also used to determine direction of flow)
const EvalWell well_pressure = bhp + cdp;
const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
const EvalWell drawdown = pressure - well_pressure;
// producing perforations
@ -316,9 +315,10 @@ namespace Opm
return;
}
const double Tw = well_index_[perf];
// compute component volumetric rates at standard conditions
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
const EvalWell cq_p = - Tw * (mob_perfcells_dense[componentIdx] * drawdown);
const EvalWell cq_p = - Tw * (mob[componentIdx] * drawdown);
cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
}
@ -347,14 +347,21 @@ namespace Opm
}
// Using total mobilities
EvalWell total_mob_dense = mob_perfcells_dense[0];
EvalWell total_mob_dense = mob[0];
for (int componentIdx = 1; componentIdx < num_components_; ++componentIdx) {
total_mob_dense += mob_perfcells_dense[componentIdx];
total_mob_dense += mob[componentIdx];
}
// injection perforations total volume rates
const double Tw = well_index_[perf];
const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown);
// surface volume fraction of fluids within wellbore
std::vector<EvalWell> cmix_s(num_components_, EvalWell{numWellEq + numEq});
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
}
// compute volume ratio between connection at standard conditions
EvalWell volumeRatio = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
@ -470,12 +477,13 @@ namespace Opm
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> cq_s(num_components_,0.0);
std::vector<EvalWell> mob(num_components_, 0.0);
getMobility(ebosSimulator, perf, mob);
std::vector<EvalWell> cq_s(num_components_, 0.0);
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
computePerfRate(intQuants, mob, bhp, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
// updating the solution gas rate and solution oil rate
@ -2166,12 +2174,13 @@ 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(num_components_, 0.0);
std::vector<EvalWell> mob(num_components_, 0.0);
getMobility(ebosSimulator, perf, mob);
std::vector<EvalWell> cq_s(num_components_, 0.0);
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
computePerfRate(intQuants, mob, bhp, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
for(int p = 0; p < np; ++p) {
@ -2560,10 +2569,11 @@ namespace Opm
// TODO: do we need to turn on crossflow here?
const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
const EvalWell& bhp = getBhp();
std::vector<EvalWell> cq_s(num_components_,0.0);
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(int_quant, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf,
computePerfRate(int_quant, mob, bhp, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate);
// TODO: make area a member
const double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf];