From dec60a8bd8c88081b492e3e8214fb98613547cbc Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 12 Aug 2016 14:34:27 +0200 Subject: [PATCH] The well model no uses the fluidState and fluidsystem from ebos The well model is modified to use the fluidState and fluidsystem from ebos. UpdateLegacyState() is no longer needed. --- opm/autodiff/BlackoilModelEbos.hpp | 409 +++++++++++++++++++++-- opm/autodiff/StandardWellsDense.hpp | 27 +- opm/autodiff/StandardWellsDense_impl.hpp | 60 ++-- 3 files changed, 433 insertions(+), 63 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index f2a0f0816..01aa91268 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -300,7 +300,7 @@ namespace Opm { makeConstantState(state0); // Compute initial accumulation contributions // and well connection pressures. - wellModel().computeWellConnectionPressures(state0, well_state); + computeWellConnectionPressures(state0, well_state); wellModel().computeAccumWells(state0); } @@ -309,17 +309,18 @@ namespace Opm { return iter_report; } + double dt = timer.currentStepLength(); std::vector mob_perfcells; std::vector b_perfcells; - wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); + //wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && iterationIdx == 0 ) { // solve the well equations as a pre-processing step iter_report = solveWellEq(mob_perfcells, b_perfcells, dt, state, well_state); } V aliveWells; - std::vector cq_s; - wellModel().computeWellFluxDense(state, mob_perfcells, b_perfcells, cq_s); + std::vector cq_s; + computeWellFluxDense(state, cq_s); wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); wellModel().addWellFluxEq(cq_s, state, dt, residual_); addWellContributionToMassBalanceEq(cq_s, state, well_state); @@ -333,6 +334,341 @@ namespace Opm { return iter_report; } + typedef DenseAd::Evaluation Eval; + EvalWell extendEval(Eval in) const { + EvalWell out = 0.0; + out.value = in.value; + for(int i = 0;i<3;++i) { + out.derivatives[i] = in.derivatives[flowToEbosPvIdx(i)]; + } + return out; + } + + template + void + computeWellFluxDense(const SolutionState& state, + std::vector& cq_s) const + { + if( ! localWellsActive() ) return ; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + V Tw = Eigen::Map(wellModel().wells().WI, nperf); + const std::vector& well_cells = wellModel().wellOps().well_cells; + std::vector well_id(nperf); + + // pressure diffs computed already (once per step, not changing per iteration) + const V& cdp = wellModel().wellPerforationPressureDiffs(); + + std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + std::vector cmix_s_ADB = wellModel().wellVolumeFractions(state); + + for (int w = 0; w < nw; ++w) { + + EvalWell bhp = wellModel().extractDenseADWell(state.bhp,w); + + + // TODO: fix for 2-phase case + std::vector cmix_s(np,0.0); + for (int phase = 0; phase < np; ++phase) { + cmix_s[phase] = wellModel().extractDenseADWell(cmix_s_ADB[phase],w); + } + + //std::cout <<"cmix gas "<< w<< " "< b_perfcells_dense(np, 0.0); + std::vector mob_perfcells_dense(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); + b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); + mob_perfcells_dense[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); + } + + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + cdp[perf]; + EvalWell drawdown = pressure - well_pressure; + + // injection perforations + if ( drawdown.value > 0 ) { + + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) + continue; + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + const EvalWell cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); + cq_ps[phase] = b_perfcells_dense[phase] * cq_p; + } + + + if ((active_)[Oil] && (active_)[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const EvalWell cq_psOil = cq_ps[oilpos]; + const EvalWell cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs * cq_psOil; + cq_ps[oilpos] += rv * cq_psGas; + } + + + // map to ADB + for (int phase = 0; phase < np; ++phase) { + cq_s_dense[phase][perf] = cq_ps[phase]; + } + + } else { + //Do nothing if crossflow is not allowed + if (!wells().allow_cf[w] && wells().type[w] == PRODUCER) + continue; + + // Using total mobilities + EvalWell total_mob_dense = mob_perfcells_dense[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob_dense += mob_perfcells_dense[phase]; + } + // injection perforations total volume rates + const EvalWell cqt_i = - Tw[perf] * (total_mob_dense * drawdown); + + // compute volume ratio between connection at standard conditions + EvalWell volumeRatio = 0.0; + if ((active_)[Water]) { + const int watpos = pu.phase_pos[Water]; + volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; + } + + if ((active_)[Oil] && (active_)[Gas]) { + EvalWell well_temperature = extendEval(fs.temperature(FluidSystem::oilPhaseIdx)); + EvalWell rsSatEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + EvalWell rvSatEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); + + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + EvalWell rvPerf = 0.0; + if (cmix_s[gaspos] > 0) + rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; + + if (rvPerf.value > rvSatEval.value) { + rvPerf = rvSatEval; + //rvPerf.value = rvSatEval.value; + } + + EvalWell rsPerf = 0.0; + if (cmix_s[oilpos] > 0) + rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; + + if (rsPerf.value > rsSatEval.value) { + //rsPerf = 0.0; + rsPerf= rsSatEval; + } + + // Incorporate RS/RV factors if both oil and gas active + const EvalWell d = 1.0 - rvPerf * rsPerf; + + const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; + //std::cout << "tmp_oil " < cq_s2; + //Vector aliveWells; + //computeWellFlux(state,mob_perfcells,b_perfcells, aliveWells,cq_s2); + + //for (int phase = 0; phase < np; ++phase) { + //if( !(((cq_s[phase].value() - cq_s2[phase].value()).abs()<1e-10).all()) ) { + //std::cout << "phase " << phase << std::endl; + //std::cout << cq_s2[phase].value() << std::endl; + //std::cout << cq_s[phase].value() << std::endl; + //} + //} + } + void + computeWellConnectionPressures(const SolutionState& state, + const WellState& xw) + { + if( ! localWellsActive() ) return ; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); + const V& pdepth = subset(depth, wellModel().wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + wellModel().computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity); + + } + + template + void + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + const std::vector& well_cells = wellModel().wellOps().well_cells; + const PhaseUsage& pu = fluid_.phaseUsage(); + const int np = fluid_.numPhases(); + b_perf.resize(nperf*np); + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + + std::vector perf_cond(nperf); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = phaseCondition_[well_cells[perf]]; + } + + // Compute the average pressure in each well block + const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); + //V avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_idx = well_cells[perf]; + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * pu.num_phases] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * pu.num_phases; + if (perf_cond[perf].hasFreeOil()) { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + else { + double rv = fs.Rv().value; + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * pu.num_phases; + if (perf_cond[perf].hasFreeGas()) { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + else { + double rs = fs.Rs().value; + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } + } + + +// // Use cell values for the temperature as the wells don't knows its temperature yet. +// const ADB perf_temp = subset(state.temperature, well_cells); + +// // Compute b, rsmax, rvmax values for perforations. +// // Evaluate the properties using average well block pressures +// // and cell values for rs, rv, phase condition and temperature. +// const ADB avg_press_ad = ADB::constant(avg_press); +// // const std::vector& pc = phaseCondition(); + +// DataBlock b(nperf, pu.num_phases); +// if (pu.phase_used[BlackoilPhases::Aqua]) { +// const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); +// b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; +// } +// assert((*active_)[Oil]); +// const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); +// if (pu.phase_used[BlackoilPhases::Liquid]) { +// const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null(); +// const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); +// b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; +// } +// if (pu.phase_used[BlackoilPhases::Vapour]) { +// const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null(); +// const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); +// b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; +// } +// if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) { +// const V rssat = fluid_.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); +// //rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + +// const V rvsat = fluid_.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); +// //rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); +// } + +// // b is row major, so can just copy data. +// //b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + + // Surface density. + // The compute density segment wants the surface densities as + // an np * number of wells cells array + V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); + + } + + + /// \brief Compute the residual norms of the mass balance for each phase, /// the well flux, and the well equation. @@ -710,8 +1046,23 @@ namespace Opm { for ( int idx = 0; idx < np; ++idx ) { - const ADB& tempB = rq_[idx].b; - B.col(idx) = 1./tempB.value(); + V b(nc); + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + + b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; + } + B.col(idx) = b; + } + + + for ( int idx = 0; idx < np; ++idx ) + { + //const ADB& tempB = rq_[idx].b; + //B.col(idx) = 1./tempB.value(); R.col(idx) = residual_.material_balance_eq[idx].value(); tempV.col(idx) = R.col(idx).abs()/pv; } @@ -1460,7 +1811,7 @@ namespace Opm { std::move(pAdbJacs[phaseIdx])); rq_[phaseIdx].mob = ADB::function(std::move(mobVal[phaseIdx]), - std::move(mobAdbJacs[phaseIdx])); + std::move(mobAdbJacs[phaseIdx])); rq_[phaseIdx].b = ADB::function(std::move(bVal[phaseIdx]), std::move(bAdbJacs[phaseIdx])); @@ -1522,7 +1873,7 @@ namespace Opm { prevEpisodeIdx = ebosSimulator_.episodeIndex(); convertResults(ebosSimulator_, /*sparsityPattern=*/state.saturation[0]); - updateLegacyState(ebosSimulator_, state); + //updateLegacyState(ebosSimulator_, state); if (param_.update_equations_scaling_) { updateEquationsScaling(); @@ -1536,7 +1887,6 @@ namespace Opm { SolutionState& state, WellState& well_state) { - V aliveWells; const int np = wells().number_of_phases; std::vector cq_s(np, ADB::null()); std::vector indices = wellModel().variableWellStateIndices(); @@ -1550,10 +1900,10 @@ namespace Opm { if (localWellsActive() ){ // If there are non well in the sudomain of the process // thene mob_perfcells_const and b_perfcells_const would be empty - for (int phase = 0; phase < np; ++phase) { - mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); - b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); - } +// for (int phase = 0; phase < np; ++phase) { +// mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); +// b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); +// } } int it = 0; @@ -1567,7 +1917,7 @@ namespace Opm { SolutionState wellSolutionState = state0; wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState, well_state); - wellModel().computeWellFluxDense(wellSolutionState, mob_perfcells_const, b_perfcells_const, cq_s); + computeWellFluxDense(wellSolutionState, cq_s); wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); wellModel().addWellFluxEq(cq_s, wellSolutionState, dt, residual_); //wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_); @@ -1672,10 +2022,25 @@ namespace Opm { Eigen::Array B(nc, np); Eigen::Array R(nc, np); Eigen::Array tempV(nc, np); + for ( int idx = 0; idx < np; ++idx ) { - const ADB& tempB = rq_[idx].b; - B.col(idx) = 1./tempB.value(); + V b(nc); + for (int cell_idx = 0; cell_idx < nc; ++cell_idx) { + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx); + + b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value; + } + B.col(idx) = b; + } + + for ( int idx = 0; idx < np; ++idx ) + { + //const ADB& tempB = rq_[idx].b; + //B.col(idx) = 1./tempB.value(); R.col(idx) = residual_.material_balance_eq[idx].value(); tempV.col(idx) = R.col(idx).abs()/pv; } @@ -1835,12 +2200,12 @@ namespace Opm { // TODO: added since the interfaces of the function are different // TODO: for StandardWells and MultisegmentWells - void - computeWellConnectionPressures(const SolutionState& state, - const WellState& well_state) - { - wellModel().computeWellConnectionPressures(state, well_state); - } +// void +// computeWellConnectionPressures(const SolutionState& state, +// const WellState& well_state) +// { +// wellModel().computeWellConnectionPressures(state, well_state); +// } /// \brief Compute the reduction within the convergence check. diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 601760fdb..6ffb2c8c8 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -60,11 +60,12 @@ namespace Opm { // --------- Types --------- using ADB = AutoDiffBlock; - typedef DenseAd::Evaluation Eval; + typedef DenseAd::Evaluation EvalWell; //typedef AutoDiffBlock ADB; using Vector = ADB::V; using V = ADB::V; + // copied from BlackoilModelBase // should put to somewhere better using DataBlock = Eigen::Array& b_perfcells, std::vector& cq_s) const; - Eval extractDenseAD(const ADB& data, int i, int j) const; - Eval extractDenseADWell(const ADB& data, int i) const; + EvalWell extractDenseAD(const ADB& data, int i, int j) const; + EvalWell extractDenseADWell(const ADB& data, int i) const; - const ADB convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const; + const ADB convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const; template @@ -215,6 +216,15 @@ namespace Opm { template void computeAccumWells(const SolutionState& state); + template + void computeWellConnectionDensitesPressures(const WellState& xw, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav); + protected: bool wells_active_; const Wells* wells_; @@ -246,14 +256,7 @@ namespace Opm { std::vector& rvmax_perf, std::vector& surf_dens_perf); - template - void computeWellConnectionDensitesPressures(const WellState& xw, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& surf_dens_perf, - const std::vector& depth_perf, - const double grav); + template diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 46fd6b85d..3374fe846 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -397,6 +397,7 @@ namespace Opm std::vector& cq_s) const { if( ! localWellsActive() ) return ; + const int np = wells().number_of_phases; const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; @@ -408,7 +409,7 @@ namespace Opm // pressure diffs computed already (once per step, not changing per iteration) const Vector& cdp = wellPerforationPressureDiffs(); - std::vector> cq_s_dense(np, std::vector(nperf,0.0)); + std::vector> cq_s_dense(np, std::vector(nperf,0.0)); std::vector cmix_s_ADB = wellVolumeFractions(state); const int oilpos = pu.phase_pos[Oil]; @@ -417,11 +418,12 @@ namespace Opm const ADB rvSat = fluid_->rvSat(perfpressure, cmix_s_ADB[oilpos], well_cells); for (int w = 0; w < nw; ++w) { - Eval bhp = extractDenseADWell(state.bhp,w); + + EvalWell bhp = extractDenseADWell(state.bhp,w); // TODO: fix for 2-phase case - std::vector cmix_s(np,0.0); + std::vector cmix_s(np,0.0); for (int phase = 0; phase < np; ++phase) { cmix_s[phase] = extractDenseADWell(cmix_s_ADB[phase],w); } @@ -431,11 +433,11 @@ namespace Opm for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { const int cell_idx = well_cells[perf]; well_id[perf] = w; - Eval pressure = extractDenseAD(state.pressure, cell_idx, cell_idx); - Eval rs = extractDenseAD(state.rs, cell_idx, cell_idx); - Eval rv = extractDenseAD(state.rv, cell_idx, cell_idx); - std::vector b_perfcells_dense(np, 0.0); - std::vector mob_perfcells_dense(np, 0.0); + EvalWell pressure = extractDenseAD(state.pressure, cell_idx, cell_idx); + EvalWell rs = extractDenseAD(state.rs, cell_idx, cell_idx); + EvalWell rv = extractDenseAD(state.rv, cell_idx, cell_idx); + std::vector b_perfcells_dense(np, 0.0); + std::vector mob_perfcells_dense(np, 0.0); for (int phase = 0; phase < np; ++phase) { b_perfcells_dense[phase] = extractDenseAD(b_perfcells[phase], perf, cell_idx); mob_perfcells_dense[phase] = extractDenseAD(mob_perfcells[phase], perf, cell_idx); @@ -443,7 +445,7 @@ namespace Opm } // Pressure drawdown (also used to determine direction of flow) - Eval drawdown = pressure - bhp - cdp[perf]; + EvalWell drawdown = pressure - bhp - cdp[perf]; // injection perforations if ( drawdown.value > 0 ) { @@ -452,9 +454,9 @@ namespace Opm if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) continue; // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, 0.0); + std::vector cq_ps(np, 0.0); for (int phase = 0; phase < np; ++phase) { - const Eval cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); + const EvalWell cq_p = - Tw[perf] * (mob_perfcells_dense[phase] * drawdown); cq_ps[phase] = b_perfcells_dense[phase] * cq_p; } @@ -462,8 +464,8 @@ namespace Opm if ((*active_)[Oil] && (*active_)[Gas]) { const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; - const Eval cq_psOil = cq_ps[oilpos]; - const Eval cq_psGas = cq_ps[gaspos]; + const EvalWell cq_psOil = cq_ps[oilpos]; + const EvalWell cq_psGas = cq_ps[gaspos]; cq_ps[gaspos] += rs * cq_psOil; cq_ps[oilpos] += rv * cq_psGas; } @@ -480,15 +482,15 @@ namespace Opm continue; // Using total mobilities - Eval total_mob_dense = mob_perfcells_dense[0]; + EvalWell total_mob_dense = mob_perfcells_dense[0]; for (int phase = 1; phase < np; ++phase) { total_mob_dense += mob_perfcells_dense[phase]; } // injection perforations total volume rates - const Eval cqt_i = - Tw[perf] * (total_mob_dense * drawdown); + const EvalWell cqt_i = - Tw[perf] * (total_mob_dense * drawdown); // compute volume ratio between connection at standard conditions - Eval volumeRatio = 0.0; + EvalWell volumeRatio = 0.0; if ((*active_)[Water]) { const int watpos = pu.phase_pos[Water]; volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; @@ -498,7 +500,7 @@ namespace Opm const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; - Eval rvPerf = 0.0; + EvalWell rvPerf = 0.0; if (cmix_s[gaspos] > 0) rvPerf = cmix_s[oilpos] / cmix_s[gaspos]; @@ -507,7 +509,7 @@ namespace Opm rvPerf.value = rvSat.value()[w]; } - Eval rsPerf = 0.0; + EvalWell rsPerf = 0.0; if (cmix_s[oilpos] > 0) rsPerf = cmix_s[gaspos] / cmix_s[oilpos]; @@ -517,13 +519,13 @@ namespace Opm } // Incorporate RS/RV factors if both oil and gas active - const Eval d = 1.0 - rvPerf * rsPerf; + const EvalWell d = 1.0 - rvPerf * rsPerf; - const Eval tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; + const EvalWell tmp_oil = (cmix_s[oilpos] - rvPerf * cmix_s[gaspos]) / d; //std::cout << "tmp_oil " < Eval; - Eval + typedef DenseAd::Evaluation EvalWell; + EvalWell StandardWellsDense:: extractDenseAD(const ADB& data, int i, int j) const { - Eval output = 0.0; + EvalWell output = 0.0; output.value = data.value()[i]; const int np = wells().number_of_phases; const std::vector& jac = data.derivative(); @@ -831,12 +833,12 @@ namespace Opm return output; } - typedef DenseAd::Evaluation Eval; - Eval + typedef DenseAd::Evaluation EvalWell; + EvalWell StandardWellsDense:: extractDenseADWell(const ADB& data, int i) const { - Eval output = 0.0; + EvalWell output = 0.0; output.value = data.value()[i]; const int nw = wells().number_of_wells; const int np = wells().number_of_phases; @@ -849,7 +851,7 @@ namespace Opm return output; } - const AutoDiffBlock StandardWellsDense::convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const + const AutoDiffBlock StandardWellsDense::convertToADB(const std::vector& local, const std::vector& well_cells, const int nc, const std::vector& well_id, const int nw, const int numVars) const { typedef typename ADB::M M; const int nLocal = local.size();