diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 209a1ed7c..fa8a01dc6 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -162,13 +162,13 @@ namespace Opm void computeConnLevelProdInd(const FluidState& fs, const std::function& connPICalc, - const std::vector& mobility, + const std::vector& mobility, double* connPI) const; void computeConnLevelInjInd(const FluidState& fs, const Phase preferred_phase, const std::function& connIICalc, - const std::vector& mobility, + const std::vector& mobility, double* connII, DeferredLogger& deferred_logger) const; @@ -198,27 +198,63 @@ namespace Opm // compute the pressure difference between the perforation and cell center void computePerfCellPressDiffs(const Simulator& ebosSimulator); - void computePerfRatePressure(const IntensiveQuantities& int_quants, - const std::vector& mob_perfcells, - const double Tw, - const int seg, - const int perf, - const EvalWell& segment_pressure, - const bool& allow_cf, - std::vector& cq_s, - EvalWell& perf_press, - double& perf_dis_gas_rate, - double& perf_vap_oil_rate, - DeferredLogger& deferred_logger) const; + void computePerfRateScalar(const IntensiveQuantities& int_quants, + const std::vector& mob_perfcells, + const double Tw, + const int seg, + const int perf, + const Scalar& segment_pressure, + const bool& allow_cf, + std::vector& cq_s, + Scalar& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const; + + void computePerfRateEval(const IntensiveQuantities& int_quants, + const std::vector& mob_perfcells, + const double Tw, + const int seg, + const int perf, + const EvalWell& segment_pressure, + const bool& allow_cf, + std::vector& cq_s, + EvalWell& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const; + + template + void computePerfRate(const Value& pressure_cell, + const Value& rs, + const Value& rv, + const std::vector& b_perfcells, + const std::vector& mob_perfcells, + const double Tw, + const int perf, + const Value& segment_pressure, + const Value& segment_density, + const bool& allow_cf, + const std::vector& cmix_s, + std::vector& cq_s, + Value& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const; // compute the fluid properties, such as densities, viscosities, and so on, in the segments // They will be treated implicitly, so they need to be of Evaluation type void computeSegmentFluidProperties(const Simulator& ebosSimulator); // get the mobility for specific perforation - void getMobility(const Simulator& ebosSimulator, - const int perf, - std::vector& mob) const; + void getMobilityEval(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const; + + // get the mobility for specific perforation + void getMobilityScalar(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const; void computeWellRatesAtBhpLimit(const Simulator& ebosSimulator, std::vector& well_flux, diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 2500afb2c..a769d7b00 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -1218,147 +1218,6 @@ getSegmentSurfaceVolume(const EvalWell& temperature, return volume / vol_ratio; } -template -void -MultisegmentWellEval:: -computePerfRatePressure(const EvalWell& pressure_cell, - const EvalWell& rs, - const EvalWell& rv, - const std::vector& b_perfcells, - const std::vector& mob_perfcells, - const double Tw, - const int seg, - const int perf, - const EvalWell& segment_pressure, - const bool& allow_cf, - std::vector& cq_s, - EvalWell& perf_press, - double& perf_dis_gas_rate, - double& perf_vap_oil_rate, - DeferredLogger& deferred_logger) const -{ - std::vector cmix_s(baseif_.numComponents(), 0.0); - - // the composition of the components inside wellbore - for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { - cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx); - } - - - // pressure difference between the segment and the perforation - const EvalWell perf_seg_press_diff = baseif_.gravity() * segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf]; - // pressure difference between the perforation and the grid cell - const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; - - perf_press = pressure_cell - cell_perf_press_diff; - // Pressure drawdown (also used to determine direction of flow) - // TODO: not 100% sure about the sign of the seg_perf_press_diff - const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff); - - // producing perforations - if ( drawdown > 0.0) { - // Do nothing is crossflow is not allowed - if (!allow_cf && baseif_.isInjector()) { - return; - } - - // compute component volumetric rates at standard conditions - for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { - const EvalWell cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown); - cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const EvalWell cq_s_oil = cq_s[oilCompIdx]; - const EvalWell cq_s_gas = cq_s[gasCompIdx]; - cq_s[gasCompIdx] += rs * cq_s_oil; - cq_s[oilCompIdx] += rv * cq_s_gas; - } - } else { // injecting perforations - // Do nothing if crossflow is not allowed - if (!allow_cf && baseif_.isProducer()) { - return; - } - - // for injecting perforations, we use total mobility - EvalWell total_mob = mob_perfcells[0]; - for (int comp_idx = 1; comp_idx < baseif_.numComponents(); ++comp_idx) { - total_mob += mob_perfcells[comp_idx]; - } - - // injection perforations total volume rates - const EvalWell cqt_i = - Tw * (total_mob * drawdown); - - // compute volume ratio between connection and at standard conditions - EvalWell volume_ratio = 0.0; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - - // Incorporate RS/RV factors if both oil and gas active - // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations - // basically, for injecting perforations, the wellbore is the upstreaming side. - const EvalWell d = 1.0 - rv * rs; - - if (d.value() == 0.0) { - OPM_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << baseif_.name() - << " during flux calculation" - << " with rs " << rs << " and rv " << rv, deferred_logger); - } - - const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d; - volume_ratio += tmp_oil / b_perfcells[oilCompIdx]; - - const EvalWell tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d; - volume_ratio += tmp_gas / b_perfcells[gasCompIdx]; - } else { // not having gas and oil at the same time - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx]; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx]; - } - } - // injecting connections total volumerates at standard conditions - EvalWell cqt_is = cqt_i / volume_ratio; - for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) { - cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; - } - } // end for injection perforations - - // calculating the perforation solution gas rate and solution oil rates - if (baseif_.isProducer()) { - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells - // s means standard condition, r means reservoir condition - // q_os = q_or * b_o + rv * q_gr * b_g - // q_gs = q_gr * g_g + rs * q_or * b_o - // d = 1.0 - rs * rv - // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) - // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) - - const double d = 1.0 - rv.value() * rs.value(); - // vaporized oil into gas - // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d - perf_vap_oil_rate = rv.value() * (cq_s[gasCompIdx].value() - rs.value() * cq_s[oilCompIdx].value()) / d; - // dissolved of gas in oil - // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d - perf_dis_gas_rate = rs.value() * (cq_s[oilCompIdx].value() - rv.value() * cq_s[gasCompIdx].value()) / d; - } - } -} - template void MultisegmentWellEval:: diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index ad2e69ce3..bebc19d04 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -150,22 +150,6 @@ protected: const double max_residual_allowed, DeferredLogger& deferred_logger) const; - void computePerfRatePressure(const EvalWell& pressure_cell, - const EvalWell& rs, - const EvalWell& rv, - const std::vector& b_perfcells, - const std::vector& mob_perfcells, - const double Tw, - const int seg, - const int perf, - const EvalWell& segment_pressure, - const bool& allow_cf, - std::vector& cq_s, - EvalWell& perf_press, - double& perf_dis_gas_rate, - double& perf_vap_oil_rate, - DeferredLogger& deferred_logger) const; - /// check whether the well equations get converged for this well ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 3c62a01f0..870730e19 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -615,8 +615,8 @@ namespace Opm return wellPICalc.connectionProdIndStandard(allPerfID, mobility); }; - std::vector mob(this->num_components_, 0.0); - getMobility(ebosSimulator, static_cast(subsetPerfID), mob); + std::vector mob(this->num_components_, 0.0); + getMobilityScalar(ebosSimulator, static_cast(subsetPerfID), mob); const auto& fs = fluidState(subsetPerfID); setToZero(connPI); @@ -677,23 +677,156 @@ namespace Opm + template + template + void + MultisegmentWell:: + computePerfRate(const Value& pressure_cell, + const Value& rs, + const Value& rv, + const std::vector& b_perfcells, + const std::vector& mob_perfcells, + const double Tw, + const int perf, + const Value& segment_pressure, + const Value& segment_density, + const bool& allow_cf, + const std::vector& cmix_s, + std::vector& cq_s, + Value& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const + { + // pressure difference between the segment and the perforation + const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf]; + // pressure difference between the perforation and the grid cell + const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; + perf_press = pressure_cell - cell_perf_press_diff; + // Pressure drawdown (also used to determine direction of flow) + // TODO: not 100% sure about the sign of the seg_perf_press_diff + const Value drawdown = perf_press - (segment_pressure + perf_seg_press_diff); + + // producing perforations + if ( drawdown > 0.0) { + // Do nothing is crossflow is not allowed + if (!allow_cf && this->isInjector()) { + return; + } + + // compute component volumetric rates at standard conditions + for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) { + const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown); + cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const Value cq_s_oil = cq_s[oilCompIdx]; + const Value cq_s_gas = cq_s[gasCompIdx]; + cq_s[gasCompIdx] += rs * cq_s_oil; + cq_s[oilCompIdx] += rv * cq_s_gas; + } + } else { // injecting perforations + // Do nothing if crossflow is not allowed + if (!allow_cf && this->isProducer()) { + return; + } + + // for injecting perforations, we use total mobility + Value total_mob = mob_perfcells[0]; + for (int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) { + total_mob += mob_perfcells[comp_idx]; + } + + // injection perforations total volume rates + const Value cqt_i = - Tw * (total_mob * drawdown); + + // compute volume ratio between connection and at standard conditions + Value volume_ratio = 0.0; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + + // Incorporate RS/RV factors if both oil and gas active + // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations + // basically, for injecting perforations, the wellbore is the upstreaming side. + const Value d = 1.0 - rv * rs; + + if (getValue(d) == 0.0) { + OPM_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << this->name() + << " during flux calculation" + << " with rs " << rs << " and rv " << rv, deferred_logger); + } + + const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d; + volume_ratio += tmp_oil / b_perfcells[oilCompIdx]; + + const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d; + volume_ratio += tmp_gas / b_perfcells[gasCompIdx]; + } else { // not having gas and oil at the same time + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx]; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx]; + } + } + // injecting connections total volumerates at standard conditions + Value cqt_is = cqt_i / volume_ratio; + for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) { + cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; + } + } // end for injection perforations + + // calculating the perforation solution gas rate and solution oil rates + if (this->isProducer()) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells + // s means standard condition, r means reservoir condition + // q_os = q_or * b_o + rv * q_gr * b_g + // q_gs = q_gr * g_g + rs * q_or * b_o + // d = 1.0 - rs * rv + // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) + // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) + + const double d = 1.0 - getValue(rv) * getValue(rs); + // vaporized oil into gas + // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d + perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d; + // dissolved of gas in oil + // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d + perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d; + } + } + } template void MultisegmentWell:: - computePerfRatePressure(const IntensiveQuantities& int_quants, - const std::vector& mob_perfcells, - const double Tw, - const int seg, - const int perf, - const EvalWell& segment_pressure, - const bool& allow_cf, - std::vector& cq_s, - EvalWell& perf_press, - double& perf_dis_gas_rate, - double& perf_vap_oil_rate, - DeferredLogger& deferred_logger) const + computePerfRateEval(const IntensiveQuantities& int_quants, + const std::vector& mob_perfcells, + const double Tw, + const int seg, + const int perf, + const EvalWell& segment_pressure, + const bool& allow_cf, + std::vector& cq_s, + EvalWell& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const { const auto& fs = int_quants.fluidState(); @@ -714,26 +847,88 @@ namespace Opm b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx)); } - this->MSWEval::computePerfRatePressure(pressure_cell, - rs, - rv, - b_perfcells, - mob_perfcells, - Tw, - seg, - perf, - segment_pressure, - allow_cf, - cq_s, - perf_press, - perf_dis_gas_rate, - perf_vap_oil_rate, - deferred_logger); + std::vector cmix_s(this->numComponents(), 0.0); + for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) { + cmix_s[comp_idx] = this->surfaceVolumeFraction(seg, comp_idx); + } + + this->computePerfRate(pressure_cell, + rs, + rv, + b_perfcells, + mob_perfcells, + Tw, + perf, + segment_pressure, + this->segment_densities_[seg], + allow_cf, + cmix_s, + cq_s, + perf_press, + perf_dis_gas_rate, + perf_vap_oil_rate, + deferred_logger); } + template + void + MultisegmentWell:: + computePerfRateScalar(const IntensiveQuantities& int_quants, + const std::vector& mob_perfcells, + const double Tw, + const int seg, + const int perf, + const Scalar& segment_pressure, + const bool& allow_cf, + std::vector& cq_s, + Scalar& perf_press, + double& perf_dis_gas_rate, + double& perf_vap_oil_rate, + DeferredLogger& deferred_logger) const + { + const auto& fs = int_quants.fluidState(); + + const Scalar pressure_cell = getValue(fs.pressure(FluidSystem::oilPhaseIdx)); + const Scalar rs = getValue(fs.Rs()); + const Scalar rv = getValue(fs.Rv()); + + // not using number_of_phases_ because of solvent + std::vector b_perfcells(this->num_components_, 0.0); + + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + b_perfcells[compIdx] = getValue(fs.invB(phaseIdx)); + } + + std::vector cmix_s(this->numComponents(), 0.0); + for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) { + cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx)); + } + + this->computePerfRate(pressure_cell, + rs, + rv, + b_perfcells, + mob_perfcells, + Tw, + perf, + segment_pressure, + getValue(this->segment_densities_[seg]), + allow_cf, + cmix_s, + cq_s, + perf_press, + perf_dis_gas_rate, + perf_vap_oil_rate, + deferred_logger); + } template void @@ -777,9 +972,9 @@ namespace Opm template void MultisegmentWell:: - getMobility(const Simulator& ebosSimulator, - const int perf, - std::vector& mob) const + getMobilityEval(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const { // TODO: most of this function, if not the whole function, can be moved to the base class const int cell_idx = this->well_cells_[perf]; @@ -826,6 +1021,58 @@ namespace Opm } + template + void + MultisegmentWell:: + getMobilityScalar(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const + { + // TODO: most of this function, if not the whole function, can be moved to the base class + const int cell_idx = this->well_cells_[perf]; + assert (int(mob.size()) == this->num_components_); + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); + + // either use mobility of the perforation cell or calcualte its own + // based on passing the saturation table index + const int satid = this->saturation_table_number_[perf] - 1; + const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); + if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell + + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx)); + } + // if (has_solvent) { + // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); + // } + } else { + + const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); + Scalar relativePerms[3] = { 0.0, 0.0, 0.0 }; + MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); + + // reset the satnumvalue back to original + materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); + + // compute the mobility + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx)); + } + } + } + + template @@ -916,11 +1163,10 @@ namespace Opm ref_depth = segment_depth; seg_bhp_press_diff += dp; for (const int perf : this->segment_perforations_[seg]) { - //std::vector mob(this->num_components_, {numWellEq_ + numEq, 0.0}); - std::vector mob(this->num_components_, 0.0); + std::vector mob(this->num_components_, 0.0); // TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration - getMobility(ebos_simulator, perf, mob); + getMobilityScalar(ebos_simulator, perf, mob); const int cell_idx = this->well_cells_[perf]; const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); @@ -963,7 +1209,7 @@ namespace Opm std::vector ipr_a_perf(this->ipr_a_.size()); std::vector ipr_b_perf(this->ipr_b_.size()); for (int p = 0; p < this->number_of_phases_; ++p) { - const double tw_mob = tw_perf * mob[p].value() * b_perf[p]; + const double tw_mob = tw_perf * mob[p] * b_perf[p]; ipr_a_perf[p] += tw_mob * pressure_diff; ipr_b_perf[p] += tw_mob; } @@ -1274,14 +1520,14 @@ namespace Opm const int cell_idx = this->well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); std::vector mob(this->num_components_, 0.0); - getMobility(ebosSimulator, perf, mob); + getMobilityEval(ebosSimulator, perf, mob); const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; std::vector cq_s(this->num_components_, 0.0); EvalWell perf_press; double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - computePerfRatePressure(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); + computePerfRateEval(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); // updating the solution gas rate and solution oil rate if (this->isProducer()) { @@ -1523,34 +1769,30 @@ namespace Opm DeferredLogger& deferred_logger) const { // Calculate the rates that follow from the current primary variables. - std::vector well_q_s(this->num_components_, 0.0); + std::vector well_q_s(this->num_components_, 0.0); const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { // calculating the perforation rate for each perforation that belongs to this segment - const EvalWell seg_pressure = this->getSegmentPressure(seg); + const Scalar seg_pressure = getValue(this->getSegmentPressure(seg)); for (const int perf : this->segment_perforations_[seg]) { const int cell_idx = this->well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - std::vector mob(this->num_components_, 0.0); - getMobility(ebosSimulator, perf, mob); + std::vector mob(this->num_components_, 0.0); + getMobilityScalar(ebosSimulator, perf, mob); const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; - std::vector cq_s(this->num_components_, 0.0); - EvalWell perf_press; + std::vector cq_s(this->num_components_, 0.0); + Scalar perf_press; double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - computePerfRatePressure(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); + computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); for (int comp = 0; comp < this->num_components_; ++comp) { well_q_s[comp] += cq_s[comp]; } } } - std::vector well_q_s_noderiv(well_q_s.size()); - for (int comp = 0; comp < this->num_components_; ++comp) { - well_q_s_noderiv[comp] = well_q_s[comp].value(); - } - return well_q_s_noderiv; + return well_q_s; } @@ -1562,7 +1804,7 @@ namespace Opm MultisegmentWell:: computeConnLevelProdInd(const typename MultisegmentWell::FluidState& fs, const std::function& connPICalc, - const std::vector& mobility, + const std::vector& mobility, double* connPI) const { const auto& pu = this->phaseUsage(); @@ -1571,7 +1813,7 @@ namespace Opm // Note: E100's notion of PI value phase mobility includes // the reciprocal FVF. const auto connMob = - mobility[ this->flowPhaseToEbosCompIdx(p) ].value() + mobility[ this->flowPhaseToEbosCompIdx(p) ] * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value(); connPI[p] = connPICalc(connMob); @@ -1601,7 +1843,7 @@ namespace Opm computeConnLevelInjInd(const typename MultisegmentWell::FluidState& fs, const Phase preferred_phase, const std::function& connIICalc, - const std::vector& mobility, + const std::vector& mobility, double* connII, DeferredLogger& deferred_logger) const { @@ -1627,9 +1869,8 @@ namespace Opm deferred_logger); } - const auto zero = EvalWell { 0.0 }; - const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero); - connII[phase_pos] = connIICalc(mt.value() * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value()); + const Scalar mt = std::accumulate(mobility.begin(), mobility.end(), 0.0); + connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value()); } } // namespace Opm