From c398a6e424b7d78f936adadfeac43683feaf6f92 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Apr 2016 17:26:07 +0200 Subject: [PATCH] puting computeWellFlux in StandardWells --- opm/autodiff/BlackoilModelBase.hpp | 7 - opm/autodiff/BlackoilModelBase_impl.hpp | 156 +----------------- opm/autodiff/StandardWells.hpp | 10 ++ opm/autodiff/StandardWells_impl.hpp | 154 +++++++++++++++++ .../BlackoilPolymerModel_impl.hpp | 4 +- 5 files changed, 168 insertions(+), 163 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 6e27f4c20..0e8784a81 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -397,13 +397,6 @@ namespace Opm { SolutionState& state, WellState& well_state); - void - computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s) const; - void updatePerfPhaseRatesAndPressures(const std::vector& cq_s, const SolutionState& state, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 44a9c544d..079475ee0 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -858,7 +858,7 @@ namespace detail { } V aliveWells; std::vector cq_s; - asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + asImpl().stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); asImpl().addWellFluxEq(cq_s, state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); @@ -1024,158 +1024,6 @@ namespace detail { - - template - void - BlackoilModelBase::computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - 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(wells().WI, nperf); - const std::vector& well_cells = stdWells().wellOps().well_cells; - - // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = asImpl().stdWells().wellPerforationPressureDiffs(); - // Extract needed quantities for the perforation cells - const ADB& p_perfcells = subset(state.pressure, well_cells); - const ADB& rv_perfcells = subset(state.rv, well_cells); - const ADB& rs_perfcells = subset(state.rs, well_cells); - - // Perforation pressure - const ADB perfpressure = (stdWells().wellOps().w2p * state.bhp) + cdp; - - // Pressure drawdown (also used to determine direction of flow) - const ADB drawdown = p_perfcells - perfpressure; - - // Compute vectors with zero and ones that - // selects the wanted quantities. - - // selects injection perforations - V selectInjectingPerforations = V::Zero(nperf); - // selects producing perforations - V selectProducingPerforations = V::Zero(nperf); - for (int c = 0; c < nperf; ++c){ - if (drawdown.value()[c] < 0) - selectInjectingPerforations[c] = 1; - else - selectProducingPerforations[c] = 1; - } - - // Handle cross flow - const V numInjectingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); - const V numProducingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); - for (int w = 0; w < nw; ++w) { - if (!wells().allow_cf[w]) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - // Crossflow is not allowed; reverse flow is prevented. - // At least one of the perforation must be open in order to have a meeningful - // equation to solve. For the special case where all perforations have reverse flow, - // and the target rate is non-zero all of the perforations are keept open. - if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { - selectProducingPerforations[perf] = 0.0; - } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ - selectInjectingPerforations[perf] = 0.0; - } - } - } - } - - // HANDLE FLOW INTO WELLBORE - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); - cq_ps[phase] = b_perfcells[phase] * cq_p; - } - if (active_[Oil] && active_[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const ADB cq_psOil = cq_ps[oilpos]; - const ADB cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs_perfcells * cq_psOil; - cq_ps[oilpos] += rv_perfcells * cq_psGas; - } - - // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities - ADB total_mob = mob_perfcells[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob += mob_perfcells[phase]; - } - // injection perforations total volume rates - const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); - - // compute wellbore mixture for injecting perforations - // The wellbore mixture depends on the inflow from the reservoar - // and the well injection rates. - - // compute avg. and total wellbore phase volumetric rates at standard conds - const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); - std::vector wbq(np, ADB::null()); - ADB wbqt = ADB::constant(V::Zero(nw)); - for (int phase = 0; phase < np; ++phase) { - const ADB& q_ps = stdWells().wellOps().p2w * cq_ps[phase]; - const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); - Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); - const int pos = pu.phase_pos[phase]; - wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps; - wbqt += wbq[phase]; - } - // compute wellbore mixture at standard conditions. - Selector notDeadWells_selector(wbqt.value(), Selector::Zero); - std::vector cmix_s(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const int pos = pu.phase_pos[phase]; - cmix_s[phase] = stdWells().wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); - } - - // compute volume ratio between connection at standard conditions - ADB volumeRatio = ADB::constant(V::Zero(nperf)); - const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; - for (int phase = 0; phase < np; ++phase) { - ADB tmp = cmix_s[phase]; - - if (phase == Oil && active_[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - tmp -= rv_perfcells * cmix_s[gaspos] / d; - } - if (phase == Gas && active_[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - tmp -= rs_perfcells * cmix_s[oilpos] / d; - } - volumeRatio += tmp / b_perfcells[phase]; - } - - // injecting connections total volumerates at standard conditions - ADB cqt_is = cqt_i/volumeRatio; - - // connection phase volumerates at standard conditions - cq_s.resize(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; - } - - // check for dead wells (used in the well controll equations) - aliveWells = V::Constant(nw, 1.0); - for (int w = 0; w < nw; ++w) { - if (wbqt.value()[w] == 0) { - aliveWells[w] = 0.0; - } - } - } - - - - - template void BlackoilModelBase::updatePerfPhaseRatesAndPressures(const std::vector& cq_s, const SolutionState& state, @@ -1587,7 +1435,7 @@ namespace detail { SolutionState wellSolutionState = state0; asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState); - asImpl().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); + asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); asImpl().addWellFluxEq(cq_s, wellSolutionState); asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 3d2d85e8e..632ed183d 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -107,6 +107,16 @@ namespace Opm { std::vector& mob_perfcells, std::vector& b_perfcells) const; + template + void computeWellFlux(const SolutionState& state, + const Opm::PhaseUsage& phase_usage, + const std::vector& active, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + std::vector& cq_s) const; + + protected: bool wells_active_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 55de376f8..c154bf130 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -310,4 +310,158 @@ namespace Opm } } + + + + + + template + void + StandardWells:: + computeWellFlux(const SolutionState& state, + const Opm::PhaseUsage& pu, + const std::vector& active, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + 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]; + Vector Tw = Eigen::Map(wells().WI, nperf); + const std::vector& well_cells = wellOps().well_cells; + + // pressure diffs computed already (once per step, not changing per iteration) + const Vector& cdp = wellPerforationPressureDiffs(); + // Extract needed quantities for the perforation cells + const ADB& p_perfcells = subset(state.pressure, well_cells); + const ADB& rv_perfcells = subset(state.rv, well_cells); + const ADB& rs_perfcells = subset(state.rs, well_cells); + + // Perforation pressure + const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp; + + // Pressure drawdown (also used to determine direction of flow) + const ADB drawdown = p_perfcells - perfpressure; + + // Compute vectors with zero and ones that + // selects the wanted quantities. + + // selects injection perforations + Vector selectInjectingPerforations = Vector::Zero(nperf); + // selects producing perforations + Vector selectProducingPerforations = Vector::Zero(nperf); + for (int c = 0; c < nperf; ++c){ + if (drawdown.value()[c] < 0) + selectInjectingPerforations[c] = 1; + else + selectProducingPerforations[c] = 1; + } + + // Handle cross flow + const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); + const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); + for (int w = 0; w < nw; ++w) { + if (!wells().allow_cf[w]) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + // Crossflow is not allowed; reverse flow is prevented. + // At least one of the perforation must be open in order to have a meeningful + // equation to solve. For the special case where all perforations have reverse flow, + // and the target rate is non-zero all of the perforations are keept open. + if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { + selectProducingPerforations[perf] = 0.0; + } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ + selectInjectingPerforations[perf] = 0.0; + } + } + } + } + + // HANDLE FLOW INTO WELLBORE + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); + cq_ps[phase] = b_perfcells[phase] * cq_p; + } + if (active[Oil] && active[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const ADB cq_psOil = cq_ps[oilpos]; + const ADB cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs_perfcells * cq_psOil; + cq_ps[oilpos] += rv_perfcells * cq_psGas; + } + + // HANDLE FLOW OUT FROM WELLBORE + // Using total mobilities + ADB total_mob = mob_perfcells[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob += mob_perfcells[phase]; + } + // injection perforations total volume rates + const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); + + // compute wellbore mixture for injecting perforations + // The wellbore mixture depends on the inflow from the reservoar + // and the well injection rates. + + // compute avg. and total wellbore phase volumetric rates at standard conds + const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); + std::vector wbq(np, ADB::null()); + ADB wbqt = ADB::constant(Vector::Zero(nw)); + for (int phase = 0; phase < np; ++phase) { + const ADB& q_ps = wellOps().p2w * cq_ps[phase]; + const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); + Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); + const int pos = pu.phase_pos[phase]; + wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps; + wbqt += wbq[phase]; + } + // compute wellbore mixture at standard conditions. + Selector notDeadWells_selector(wbqt.value(), Selector::Zero); + std::vector cmix_s(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const int pos = pu.phase_pos[phase]; + cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); + } + + // compute volume ratio between connection at standard conditions + ADB volumeRatio = ADB::constant(Vector::Zero(nperf)); + const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; + for (int phase = 0; phase < np; ++phase) { + ADB tmp = cmix_s[phase]; + + if (phase == Oil && active[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + tmp -= rv_perfcells * cmix_s[gaspos] / d; + } + if (phase == Gas && active[Oil]) { + const int oilpos = pu.phase_pos[Oil]; + tmp -= rs_perfcells * cmix_s[oilpos] / d; + } + volumeRatio += tmp / b_perfcells[phase]; + } + + // injecting connections total volumerates at standard conditions + ADB cqt_is = cqt_i/volumeRatio; + + // connection phase volumerates at standard conditions + cq_s.resize(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; + } + + // check for dead wells (used in the well controll equations) + aliveWells = Vector::Constant(nw, 1.0); + for (int w = 0; w < nw; ++w) { + if (wbqt.value()[w] == 0) { + aliveWells[w] = 0.0; + } + } + } + } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index c50cf38d8..53a81315e 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -550,7 +550,7 @@ namespace Opm { Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); } - Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); if (has_plyshlog_) { std::vector water_vel_wells; @@ -569,7 +569,7 @@ namespace Opm { mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; } - Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); Base::addWellFluxEq(cq_s, state); addWellContributionToMassBalanceEq(cq_s, state, well_state);