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

drop aliases for Indices entries

Browse Source 
using Indices directly more clearly shows where the data comes
from without having to hop through hoops to do so.
This commit is contained in:
Arne Morten Kvarving
2021-09-06 12:47:21 +02:00
parent 87b618e93a
commit 2b1ac22c99
7 changed files with 67 additions and 93 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@@ -1050,7 +1050,7 @@ namespace Opm {
getWellContributions(WellContributions& wellContribs) const getWellContributions(WellContributions& wellContribs) const
{ {
// prepare for StandardWells // prepare for StandardWells
wellContribs.setBlockSize(StandardWell<TypeTag>::numEq, StandardWell<TypeTag>::numStaticWellEq); wellContribs.setBlockSize(StandardWell<TypeTag>::Indices::numEq, StandardWell<TypeTag>::numStaticWellEq);
for(unsigned int i = 0; i < well_container_.size(); i++){ for(unsigned int i = 0; i < well_container_.size(); i++){
auto& well = well_container_[i]; auto& well = well_container_[i];

4
opm/simulators/wells/MultisegmentWell.hpp
View File

@@ -58,10 +58,6 @@ namespace Opm
using typename Base::GLiftWellStateMap; using typename Base::GLiftWellStateMap;
using typename Base::GLiftSyncGroups; using typename Base::GLiftSyncGroups;
/// the number of reservior equations
using Base::numEq;
using Base::numPhases;
using Base::has_solvent; using Base::has_solvent;
using Base::has_polymer; using Base::has_polymer;
using Base::Water; using Base::Water;

20
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@@ -1219,7 +1219,7 @@ namespace Opm
this->resWell_[seg][comp_idx] += accumulation_term.value(); this->resWell_[seg][comp_idx] += accumulation_term.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq); this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
} }
} }
} }
@@ -1232,12 +1232,12 @@ namespace Opm
// segment_rate contains the derivatives with respect to GTotal in seg, // segment_rate contains the derivatives with respect to GTotal in seg,
// and WFrac and GFrac in seg_upwind // and WFrac and GFrac in seg_upwind
this->resWell_[seg][comp_idx] -= segment_rate.value(); this->resWell_[seg][comp_idx] -= segment_rate.value();
this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + numEq); this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq); this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
} }
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq); this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
} }
// pressure derivative should be zero // pressure derivative should be zero
} }
@@ -1253,12 +1253,12 @@ namespace Opm
// inlet_rate contains the derivatives with respect to GTotal in inlet, // inlet_rate contains the derivatives with respect to GTotal in inlet,
// and WFrac and GFrac in inlet_upwind // and WFrac and GFrac in inlet_upwind
this->resWell_[seg][comp_idx] += inlet_rate.value(); this->resWell_[seg][comp_idx] += inlet_rate.value();
this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + numEq); this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq); this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
} }
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq); this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
} }
// pressure derivative should be zero // pressure derivative should be zero
} }
@@ -1308,13 +1308,13 @@ namespace Opm
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // intput in transformed matrix
// the index name for the D should be eq_idx / pv_idx // the index name for the D should be eq_idx / pv_idx
this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + numEq); this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
} }
for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
} }

20
opm/simulators/wells/StandardWell.hpp
View File

@@ -88,9 +88,6 @@ namespace Opm
using typename Base::GLiftWellStateMap; using typename Base::GLiftWellStateMap;
using typename Base::GLiftSyncGroups; using typename Base::GLiftSyncGroups;
using Base::numEq;
using Base::numPhases;
using Base::has_solvent; using Base::has_solvent;
using Base::has_zFraction; using Base::has_zFraction;
using Base::has_polymer; using Base::has_polymer;
@@ -103,20 +100,18 @@ namespace Opm
using FoamModule = BlackOilFoamModule<TypeTag>; using FoamModule = BlackOilFoamModule<TypeTag>;
using BrineModule = BlackOilBrineModule<TypeTag>; using BrineModule = BlackOilBrineModule<TypeTag>;
static const int numSolventEq = Indices::numSolvents;
// number of the conservation equations // number of the conservation equations
static const int numWellConservationEq = numPhases + numSolventEq; static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
// number of the well control equations // number of the well control equations
static const int numWellControlEq = 1; static constexpr int numWellControlEq = 1;
// number of the well equations that will always be used // number of the well equations that will always be used
// based on the solution strategy, there might be other well equations be introduced // based on the solution strategy, there might be other well equations be introduced
static const int numStaticWellEq = numWellConservationEq + numWellControlEq; static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
// the index for Bhp in primary variables and also the index of well control equation // 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. // 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 // TODO: we should have indices for the well equations and well primary variables separately
static const int Bhp = numStaticWellEq - numWellControlEq; static constexpr int Bhp = numStaticWellEq - numWellControlEq;
using typename Base::Scalar; using typename Base::Scalar;
@@ -132,13 +127,6 @@ namespace Opm
using EvalWell = typename StdWellEval::EvalWell; using EvalWell = typename StdWellEval::EvalWell;
using BVectorWell = typename StdWellEval::BVectorWell; using BVectorWell = typename StdWellEval::BVectorWell;
using Base::contiSolventEqIdx;
using Base::contiZfracEqIdx;
using Base::contiPolymerEqIdx;
using Base::contiFoamEqIdx;
using Base::contiBrineEqIdx;
static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
StandardWell(const Well& well, StandardWell(const Well& well,
const ParallelWellInfo& pw_info, const ParallelWellInfo& pw_info,
const int time_step, const int time_step,

90
opm/simulators/wells/StandardWell_impl.hpp
View File

@@ -118,7 +118,7 @@ namespace Opm
const EvalWell pressure = this->extendEval(getPerfCellPressure(fs)); const EvalWell pressure = this->extendEval(getPerfCellPressure(fs));
const EvalWell rs = this->extendEval(fs.Rs()); const EvalWell rs = this->extendEval(fs.Rs());
const EvalWell rv = this->extendEval(fs.Rv()); const EvalWell rv = this->extendEval(fs.Rv());
std::vector<EvalWell> b_perfcells_dense(num_components_, EvalWell{this->numWellEq_ + numEq, 0.0}); std::vector<EvalWell> b_perfcells_dense(num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) { if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue; continue;
@@ -127,7 +127,7 @@ namespace Opm
b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx)); b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx));
} }
if constexpr (has_solvent) { if constexpr (has_solvent) {
b_perfcells_dense[contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor()); b_perfcells_dense[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor());
} }
if constexpr (has_zFraction) { if constexpr (has_zFraction) {
@@ -138,7 +138,7 @@ namespace Opm
} }
} }
EvalWell skin_pressure = EvalWell{this->numWellEq_ + numEq, 0.0}; EvalWell skin_pressure = EvalWell{this->numWellEq_ + Indices::numEq, 0.0};
if (has_polymermw) { if (has_polymermw) {
if (this->isInjector()) { if (this->isInjector()) {
const int pskin_index = Bhp + 1 + this->numPerfs() + perf; const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
@@ -147,7 +147,7 @@ namespace Opm
} }
// surface volume fraction of fluids within wellbore // surface volume fraction of fluids within wellbore
std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + numEq}); std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + Indices::numEq});
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) { for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx); cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
} }
@@ -194,7 +194,7 @@ namespace Opm
b_perfcells_dense[compIdx] = fs.invB(phaseIdx).value(); b_perfcells_dense[compIdx] = fs.invB(phaseIdx).value();
} }
if constexpr (has_solvent) { if constexpr (has_solvent) {
b_perfcells_dense[contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value(); b_perfcells_dense[Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
} }
if constexpr (has_zFraction) { if constexpr (has_zFraction) {
@@ -435,9 +435,9 @@ namespace Opm
auto& perf_rates = perf_data.phase_rates; auto& perf_rates = perf_data.phase_rates;
for (int perf = 0; perf < number_of_perforations_; ++perf) { for (int perf = 0; perf < number_of_perforations_; ++perf) {
// Calculate perforation quantities. // Calculate perforation quantities.
std::vector<EvalWell> cq_s(num_components_, {this->numWellEq_ + numEq, 0.0}); std::vector<EvalWell> cq_s(num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
EvalWell water_flux_s{this->numWellEq_ + numEq, 0.0}; EvalWell water_flux_s{this->numWellEq_ + Indices::numEq, 0.0};
EvalWell cq_s_zfrac_effective{this->numWellEq_ + numEq, 0.0}; EvalWell cq_s_zfrac_effective{this->numWellEq_ + Indices::numEq, 0.0};
calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger); calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
// Equation assembly for this perforation. // Equation assembly for this perforation.
@@ -459,16 +459,16 @@ namespace Opm
// assemble the jacobians // assemble the jacobians
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
this->invDuneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+numEq); this->invDuneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
} }
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx); this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
} }
// Store the perforation phase flux for later usage. // Store the perforation phase flux for later usage.
if (has_solvent && componentIdx == contiSolventEqIdx) { if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
auto& perf_rate_solvent = perf_data.solvent_rates; auto& perf_rate_solvent = perf_data.solvent_rates;
perf_rate_solvent[perf] = cq_s[componentIdx].value(); perf_rate_solvent[perf] = cq_s[componentIdx].value();
} else { } else {
@@ -478,7 +478,7 @@ namespace Opm
if constexpr (has_zFraction) { if constexpr (has_zFraction) {
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
this->duneC_[0][cell_idx][pvIdx][contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+numEq); this->duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
} }
} }
} }
@@ -492,14 +492,14 @@ namespace Opm
for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) { for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
// TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
// since all the rates are under surface condition // since all the rates are under surface condition
EvalWell resWell_loc(this->numWellEq_ + numEq, 0.0); EvalWell resWell_loc(this->numWellEq_ + Indices::numEq, 0.0);
if (FluidSystem::numActivePhases() > 1) { if (FluidSystem::numActivePhases() > 1) {
assert(dt > 0); assert(dt > 0);
resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt; resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
} }
resWell_loc -= this->getQs(componentIdx) * well_efficiency_factor_; resWell_loc -= this->getQs(componentIdx) * well_efficiency_factor_;
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
this->invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq); this->invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
} }
this->resWell_[0][componentIdx] += resWell_loc.value(); this->resWell_[0][componentIdx] += resWell_loc.value();
} }
@@ -536,7 +536,7 @@ namespace Opm
const EvalWell& bhp = this->getBhp(); const EvalWell& bhp = this->getBhp();
const int cell_idx = well_cells_[perf]; const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> mob(num_components_, {this->numWellEq_ + numEq, 0.}); std::vector<EvalWell> mob(num_components_, {this->numWellEq_ + Indices::numEq, 0.});
getMobilityEval(ebosSimulator, perf, mob, deferred_logger); getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
double perf_dis_gas_rate = 0.; double perf_dis_gas_rate = 0.;
@@ -567,7 +567,7 @@ namespace Opm
} }
if constexpr (has_energy) { if constexpr (has_energy) {
connectionRates[perf][contiEnergyEqIdx] = 0.0; connectionRates[perf][Indices::contiEnergyEqIdx] = 0.0;
} }
if constexpr (has_energy) { if constexpr (has_energy) {
@@ -580,7 +580,7 @@ namespace Opm
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
// convert to reservoar conditions // convert to reservoar conditions
EvalWell cq_r_thermal(this->numWellEq_ + numEq, 0.); EvalWell cq_r_thermal(this->numWellEq_ + Indices::numEq, 0.);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if(FluidSystem::waterPhaseIdx == phaseIdx) if(FluidSystem::waterPhaseIdx == phaseIdx)
@@ -623,7 +623,7 @@ namespace Opm
} }
// compute the thermal flux // compute the thermal flux
cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx)); cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
connectionRates[perf][contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal); connectionRates[perf][Indices::contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal);
} }
} }
@@ -641,7 +641,7 @@ namespace Opm
perf_rate_polymer[perf] = cq_s_poly.value(); perf_rate_polymer[perf] = cq_s_poly.value();
cq_s_poly *= well_efficiency_factor_; cq_s_poly *= well_efficiency_factor_;
connectionRates[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly); connectionRates[perf][Indices::contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
if constexpr (Base::has_polymermw) { if constexpr (Base::has_polymermw) {
updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger); updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger);
@@ -657,7 +657,7 @@ namespace Opm
} else { } else {
cq_s_foam *= this->extendEval(intQuants.foamConcentration()); cq_s_foam *= this->extendEval(intQuants.foamConcentration());
} }
connectionRates[perf][contiFoamEqIdx] = Base::restrictEval(cq_s_foam); connectionRates[perf][Indices::contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
} }
if constexpr (has_zFraction) { if constexpr (has_zFraction) {
@@ -674,7 +674,7 @@ namespace Opm
perf_rate_solvent[perf] = cq_s_zfrac_effective.value(); perf_rate_solvent[perf] = cq_s_zfrac_effective.value();
cq_s_zfrac_effective *= well_efficiency_factor_; cq_s_zfrac_effective *= well_efficiency_factor_;
connectionRates[perf][contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective); connectionRates[perf][Indices::contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective);
} }
if constexpr (has_brine) { if constexpr (has_brine) {
@@ -691,7 +691,7 @@ namespace Opm
perf_rate_brine[perf] = cq_s_sm.value(); perf_rate_brine[perf] = cq_s_sm.value();
cq_s_sm *= well_efficiency_factor_; cq_s_sm *= well_efficiency_factor_;
connectionRates[perf][contiBrineEqIdx] = Base::restrictEval(cq_s_sm); connectionRates[perf][Indices::contiBrineEqIdx] = Base::restrictEval(cq_s_sm);
} }
// Store the perforation pressure for later usage. // Store the perforation pressure for later usage.
@@ -729,7 +729,7 @@ namespace Opm
mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx)); mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
} }
if (has_solvent) { if (has_solvent) {
mob[contiSolventEqIdx] = this->extendEval(intQuants.solventMobility()); mob[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventMobility());
} }
} else { } else {
@@ -798,7 +798,7 @@ namespace Opm
mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx)); mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
} }
if (has_solvent) { if (has_solvent) {
mob[contiSolventEqIdx] = getValue(intQuants.solventMobility()); mob[Indices::contiSolventEqIdx] = getValue(intQuants.solventMobility());
} }
} else { } else {
@@ -834,7 +834,7 @@ namespace Opm
// for the cases related to polymer molecular weight, we assume fully mixing // for the cases related to polymer molecular weight, we assume fully mixing
// as a result, the polymer and water share the same viscosity // as a result, the polymer and water share the same viscosity
if constexpr (!Base::has_polymermw) { if constexpr (!Base::has_polymermw) {
std::vector<EvalWell> mob_eval(num_components_, {this->numWellEq_ + numEq, 0.}); std::vector<EvalWell> mob_eval(num_components_, {this->numWellEq_ + Indices::numEq, 0.});
updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger); updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
for (size_t i = 0; i < mob.size(); ++i) { for (size_t i = 0; i < mob.size(); ++i) {
mob[i] = getValue(mob_eval[i]); mob[i] = getValue(mob_eval[i]);
@@ -938,7 +938,7 @@ namespace Opm
std::fill(ipr_b_.begin(), ipr_b_.end(), 0.); std::fill(ipr_b_.begin(), ipr_b_.end(), 0.);
for (int perf = 0; perf < number_of_perforations_; ++perf) { for (int perf = 0; perf < number_of_perforations_; ++perf) {
std::vector<EvalWell> mob(num_components_, {this->numWellEq_ + numEq, 0.0}); std::vector<EvalWell> mob(num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
// TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration // TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobilityEval(ebos_simulator, perf, mob, deferred_logger); getMobilityEval(ebos_simulator, perf, mob, deferred_logger);
@@ -1307,8 +1307,8 @@ namespace Opm
// We use cell values for solvent injector // We use cell values for solvent injector
if constexpr (has_solvent) { if constexpr (has_solvent) {
b_perf[num_components_ * perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value(); b_perf[num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
surf_dens_perf[num_components_ * perf + contiSolventEqIdx] = intQuants.solventRefDensity(); surf_dens_perf[num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventRefDensity();
} }
} }
} }
@@ -1392,7 +1392,7 @@ namespace Opm
return wellPICalc.connectionProdIndStandard(allPerfID, mobility); return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
}; };
std::vector<EvalWell> mob(num_components_, {this->numWellEq_ + numEq, 0.0}); std::vector<EvalWell> mob(num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
getMobilityEval(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger); getMobilityEval(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
const auto& fs = fluidState(subsetPerfID); const auto& fs = fluidState(subsetPerfID);
@@ -1451,7 +1451,7 @@ namespace Opm
if constexpr (has_solvent) { if constexpr (has_solvent) {
const auto& solvent_perf_rates_state = perf_data.solvent_rates; const auto& solvent_perf_rates_state = perf_data.solvent_rates;
for (int perf = 0; perf < nperf; ++perf) { for (int perf = 0; perf < nperf; ++perf) {
perfRates[perf * num_components_ + contiSolventEqIdx] = solvent_perf_rates_state[perf]; perfRates[perf * num_components_ + Indices::contiSolventEqIdx] = solvent_perf_rates_state[perf];
} }
} }
@@ -1482,7 +1482,7 @@ namespace Opm
perfRates[perf * num_components_ + p] = well_tw_fraction * intQuants.mobility(ebosPhaseIdx).value() / total_mobility; perfRates[perf * num_components_ + p] = well_tw_fraction * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
} }
if constexpr (has_solvent) { if constexpr (has_solvent) {
perfRates[perf * num_components_ + contiSolventEqIdx] = well_tw_fraction * intQuants.solventInverseFormationVolumeFactor().value() / total_mobility; perfRates[perf * num_components_ + Indices::contiSolventEqIdx] = well_tw_fraction * intQuants.solventInverseFormationVolumeFactor().value() / total_mobility;
} }
} }
} }
@@ -1976,7 +1976,7 @@ namespace Opm
const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator); const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
const EvalWell& bhp = this->getBhp(); const EvalWell& bhp = this->getBhp();
std::vector<EvalWell> cq_s(num_components_, {this->numWellEq_ + numEq, 0.}); std::vector<EvalWell> cq_s(num_components_, {this->numWellEq_ + Indices::numEq, 0.});
double perf_dis_gas_rate = 0.; double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.; double perf_vap_oil_rate = 0.;
double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx); double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
@@ -2050,9 +2050,9 @@ namespace Opm
OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger); OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger);
} }
const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id); const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
const EvalWell throughput_eval(this->numWellEq_ + numEq, throughput); const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero // the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_water(this->numWellEq_ + numEq, 0.0); EvalWell pskin_water(this->numWellEq_ + Indices::numEq, 0.0);
pskin_water = water_table_func.eval(throughput_eval, water_velocity); pskin_water = water_table_func.eval(throughput_eval, water_velocity);
return pskin_water; return pskin_water;
} else { } else {
@@ -2085,9 +2085,9 @@ namespace Opm
} }
const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id); const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
const double reference_concentration = skprpolytable.refConcentration; const double reference_concentration = skprpolytable.refConcentration;
const EvalWell throughput_eval(this->numWellEq_ + numEq, throughput); const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero // the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_poly(this->numWellEq_ + numEq, 0.0); EvalWell pskin_poly(this->numWellEq_ + Indices::numEq, 0.0);
pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs); pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
if (poly_inj_conc == reference_concentration) { if (poly_inj_conc == reference_concentration) {
return sign * pskin_poly; return sign * pskin_poly;
@@ -2116,8 +2116,8 @@ namespace Opm
if constexpr (Base::has_polymermw) { if constexpr (Base::has_polymermw) {
const int table_id = well_ecl_.getPolymerProperties().m_plymwinjtable; const int table_id = well_ecl_.getPolymerProperties().m_plymwinjtable;
const auto& table_func = PolymerModule::getPlymwinjTable(table_id); const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
const EvalWell throughput_eval(this->numWellEq_ + numEq, throughput); const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
EvalWell molecular_weight(this->numWellEq_ + numEq, 0.); EvalWell molecular_weight(this->numWellEq_ + Indices::numEq, 0.);
if (wpolymer() == 0.) { // not injecting polymer if (wpolymer() == 0.) { // not injecting polymer
return molecular_weight; return molecular_weight;
} }
@@ -2208,7 +2208,7 @@ namespace Opm
const double throughput = perf_water_throughput[perf]; const double throughput = perf_water_throughput[perf];
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf; const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
EvalWell poly_conc(this->numWellEq_ + numEq, 0.0); EvalWell poly_conc(this->numWellEq_ + Indices::numEq, 0.0);
poly_conc.setValue(wpolymer()); poly_conc.setValue(wpolymer());
// equation for the skin pressure // equation for the skin pressure
@@ -2217,12 +2217,12 @@ namespace Opm
this->resWell_[0][pskin_index] = eq_pskin.value(); this->resWell_[0][pskin_index] = eq_pskin.value();
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) { for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
this->invDuneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+numEq); this->invDuneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
this->invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+numEq); this->invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
} }
// the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx); this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
} }
} }
@@ -2284,7 +2284,7 @@ namespace Opm
cq_s_polymw *= 0.; cq_s_polymw *= 0.;
} }
} }
connectionRates[perf][this->contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw); connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
} }
@@ -2507,7 +2507,7 @@ namespace Opm
deferred_logger); deferred_logger);
} }
const auto zero = EvalWell { this->numWellEq_ + this->numEq, 0.0 }; const auto zero = EvalWell { this->numWellEq_ + Indices::numEq, 0.0 };
const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero); const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
connII[phase_pos] = connIICalc(mt.value() * fs.invB(flowPhaseToEbosPhaseIdx(phase_pos)).value()); connII[phase_pos] = connIICalc(mt.value() * fs.invB(flowPhaseToEbosPhaseIdx(phase_pos)).value());
} }

20
opm/simulators/wells/WellInterface.hpp
View File

@@ -88,14 +88,12 @@ public:
typename BlackoilWellModel<TypeTag>::GLiftWellStateMap; typename BlackoilWellModel<TypeTag>::GLiftWellStateMap;
using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups; using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups;
static const int numEq = Indices::numEq;
static const int numPhases = Indices::numPhases;
using Scalar = GetPropType<TypeTag, Properties::Scalar>; using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using VectorBlockType = Dune::FieldVector<Scalar, numEq>; using VectorBlockType = Dune::FieldVector<Scalar, Indices::numEq>;
using MatrixBlockType = Dune::FieldMatrix<Scalar, numEq, numEq>; using MatrixBlockType = Dune::FieldMatrix<Scalar, Indices::numEq, Indices::numEq>;
using BVector = Dune::BlockVector<VectorBlockType>; using BVector = Dune::BlockVector<VectorBlockType>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/numEq>; using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
using RateConverterType = using RateConverterType =
typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType; typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType;
@@ -114,21 +112,13 @@ public:
static constexpr bool has_polymermw = getPropValue<TypeTag, Properties::EnablePolymerMW>(); static constexpr bool has_polymermw = getPropValue<TypeTag, Properties::EnablePolymerMW>();
static constexpr bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>(); static constexpr bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>();
static constexpr bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>(); static constexpr bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>();
static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
static const int contiZfracEqIdx = Indices::contiZfracEqIdx;
static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
// index for the polymer molecular weight continuity equation
static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
static const int contiBrineEqIdx = Indices::contiBrineEqIdx;
static const bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
// For the conversion between the surface volume rate and reservoir voidage rate // For the conversion between the surface volume rate and reservoir voidage rate
using FluidState = BlackOilFluidState<Eval, using FluidState = BlackOilFluidState<Eval,
FluidSystem, FluidSystem,
has_temperature, has_temperature,
has_energy, has_energy,
compositionSwitchEnabled, Indices::compositionSwitchIdx >= 0,
has_brine, has_brine,
Indices::numPhases >; Indices::numPhases >;
/// Constructor /// Constructor
@@ -235,7 +225,7 @@ public:
{ {
Eval out = 0.0; Eval out = 0.0;
out.setValue(in.value()); out.setValue(in.value());
for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { for (int eqIdx = 0; eqIdx < Indices::numEq; ++eqIdx) {
out.setDerivative(eqIdx, in.derivative(eqIdx)); out.setDerivative(eqIdx, in.derivative(eqIdx));
} }
return out; return out;

4
tests/test_wellmodel.cpp
View File

@@ -177,7 +177,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
const auto& well = wells[0]; const auto& well = wells[0];
BOOST_CHECK_EQUAL(well->name(), "PROD1"); BOOST_CHECK_EQUAL(well->name(), "PROD1");
BOOST_CHECK(well->isProducer()); BOOST_CHECK(well->isProducer());
BOOST_CHECK(well->numEq == 3); BOOST_CHECK(StandardWell::Indices::numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4); BOOST_CHECK(well->numStaticWellEq== 4);
} }
@@ -186,7 +186,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
const auto& well = wells[1]; const auto& well = wells[1];
BOOST_CHECK_EQUAL(well->name(), "INJE1"); BOOST_CHECK_EQUAL(well->name(), "INJE1");
BOOST_CHECK(well->isInjector()); BOOST_CHECK(well->isInjector());
BOOST_CHECK(well->numEq == 3); BOOST_CHECK(StandardWell::Indices::numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4); BOOST_CHECK(well->numStaticWellEq== 4);
} }
} }