22 #include <opm/autodiff/StandardWellsSolvent.hpp>
35 StandardWellsSolvent::StandardWellsSolvent(
const Wells* wells_arg, WellCollection* well_collection)
36 : Base(wells_arg, well_collection)
37 , solvent_props_(nullptr)
48 StandardWellsSolvent::initSolvent(
const SolventPropsAdFromDeck* solvent_props,
49 const int solvent_pos,
50 const bool has_solvent)
52 solvent_props_ = solvent_props;
53 solvent_pos_ = solvent_pos;
54 has_solvent_ = has_solvent;
61 template<
class SolutionState,
class WellState>
63 StandardWellsSolvent::
64 computePropertiesForWellConnectionPressures(
const SolutionState& state,
66 std::vector<double>& b_perf,
67 std::vector<double>& rsmax_perf,
68 std::vector<double>& rvmax_perf,
69 std::vector<double>& surf_dens_perf)
74 const int nperf = wells().well_connpos[wells().number_of_wells];
75 const int nw = wells().number_of_wells;
78 const Vector perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
79 Vector avg_press = perf_press*0;
80 for (
int w = 0; w < nw; ++w) {
81 for (
int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
82 const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
83 const double p_avg = (perf_press[perf] + p_above)/2;
84 avg_press[perf] = p_avg;
88 const std::vector<int>& well_cells = wellOps().well_cells;
91 const ADB perf_temp =
subset(state.temperature, well_cells);
97 std::vector<PhasePresence> perf_cond(nperf);
98 for (
int perf = 0; perf < nperf; ++perf) {
99 perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
103 DataBlock b(nperf, pu.num_phases);
105 const Vector bw = fluid_->
bWat(avg_press_ad, perf_temp, well_cells).
value();
106 if (pu.phase_used[BlackoilPhases::Aqua]) {
107 b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
110 assert((*active_)[Oil]);
111 assert((*active_)[Gas]);
112 const ADB perf_rv =
subset(state.rv, well_cells);
113 const ADB perf_rs =
subset(state.rs, well_cells);
114 const Vector perf_so =
subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
115 if (pu.phase_used[BlackoilPhases::Liquid]) {
116 const Vector bo = fluid_->
bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).
value();
118 b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
121 rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
123 rsmax_perf.assign(0.0, nperf);
125 V surf_dens_copy =
superset(fluid_->
surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
126 for (
int phase = 1; phase < pu.num_phases; ++phase) {
127 if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) {
130 surf_dens_copy +=
superset(fluid_->
surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
133 if (pu.phase_used[BlackoilPhases::Vapour]) {
137 Vector bg = fluid_->
bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).
value();
138 Vector rhog = fluid_->
surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
142 const Vector bs = solvent_props_->
bSolvent(avg_press_ad,well_cells).
value();
146 const int nc = state.pressure.size();
149 const ADB& ss = state.solvent_saturation;
150 const ADB& sg = ((*active_)[ Gas ]
151 ? state.saturation[ pu.phase_pos[ Gas ] ]
154 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
155 Vector F_solvent =
subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
157 Vector injectedSolventFraction = Eigen::Map<const Vector>(&xw.solventFraction()[0], nperf);
159 Vector isProducer = Vector::Zero(nperf);
160 Vector ones = Vector::Constant(nperf,1.0);
161 for (
int w = 0; w < nw; ++w) {
162 if(wells().type[w] == PRODUCER) {
163 for (
int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
164 isProducer[perf] = 1;
169 F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
171 bg = bg * (ones - F_solvent);
172 bg = bg + F_solvent * bs;
175 rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
177 b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
178 surf_dens_copy +=
superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases);
182 rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
184 rvmax_perf.assign(0.0, nperf);
188 b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
189 surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
192 template <
class SolutionState>
194 StandardWellsSolvent::
195 computeWellFlux(
const SolutionState& state,
196 const std::vector<ADB>& mob_perfcells,
197 const std::vector<ADB>& b_perfcells,
199 std::vector<ADB>& cq_s)
const
203 const int np = wells().number_of_phases;
204 const int nw = wells().number_of_wells;
205 const int nperf = wells().well_connpos[nw];
206 Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
207 const std::vector<int>& well_cells = wellOps().well_cells;
212 const ADB& p_perfcells =
subset(state.pressure, well_cells);
215 const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
218 const ADB drawdown = p_perfcells - perfpressure;
224 Vector selectInjectingPerforations = Vector::Zero(nperf);
226 Vector selectProducingPerforations = Vector::Zero(nperf);
227 for (
int c = 0; c < nperf; ++c){
228 if (drawdown.value()[c] < 0)
229 selectInjectingPerforations[c] = 1;
231 selectProducingPerforations[c] = 1;
235 const Vector numInjectingPerforations = (wellOps().p2w *
ADB::constant(selectInjectingPerforations)).value();
236 const Vector numProducingPerforations = (wellOps().p2w *
ADB::constant(selectProducingPerforations)).value();
237 for (
int w = 0; w < nw; ++w) {
238 if (!wells().allow_cf[w]) {
239 for (
int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
244 if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
245 selectProducingPerforations[perf] = 0.0;
246 }
else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
247 selectInjectingPerforations[perf] = 0.0;
257 for (
int phase = 0; phase < np; ++phase) {
258 cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
259 cq_ps[phase] = b_perfcells[phase] * cq_p[phase];
262 Vector ones = Vector::Constant(nperf,1.0);
265 const Opm::PhaseUsage& pu = fluid_->
phaseUsage();
266 if ((*active_)[Oil] && (*active_)[Gas]) {
267 const int oilpos = pu.phase_pos[Oil];
268 const int gaspos = pu.phase_pos[Gas];
269 const ADB cq_psOil = cq_ps[oilpos];
270 ADB cq_psGas = cq_ps[gaspos];
271 const ADB& rv_perfcells =
subset(state.rv, well_cells);
272 const ADB& rs_perfcells =
subset(state.rs, well_cells);
273 cq_ps[gaspos] += rs_perfcells * cq_psOil;
278 const ADB& ss = state.solvent_saturation;
279 const ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
281 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
282 F_gas -=
subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
283 cq_psGas = cq_psGas * F_gas;
285 cq_ps[oilpos] += rv_perfcells * cq_psGas;
290 ADB total_mob = mob_perfcells[0];
291 for (
int phase = 1; phase < np; ++phase) {
292 total_mob += mob_perfcells[phase];
295 const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
298 if (store_well_perforation_fluxes_) {
300 Vector& wf =
const_cast<Vector&
>(well_perforation_fluxes_);
302 for (
int phase = 0; phase < np; ++phase) {
303 wf += cq_p[phase].value();
312 const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
315 for (
int phase = 0; phase < np; ++phase) {
316 const ADB& q_ps = wellOps().p2w * cq_ps[phase];
317 const ADB& q_s =
subset(state.qs, Span(nw, 1, phase*nw));
318 Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
319 const int pos = pu.phase_pos[phase];
320 wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,
ADB::constant(Vector::Zero(nw)))) - q_ps;
324 Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
325 std::vector<ADB> cmix_s(np,
ADB::null());
326 for (
int phase = 0; phase < np; ++phase) {
327 const int pos = pu.phase_pos[phase];
328 cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(
ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
334 if ((*active_)[Water]) {
335 const int watpos = pu.phase_pos[Water];
336 volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
339 if ((*active_)[Oil] && (*active_)[Gas]) {
341 const ADB& rv_perfcells =
subset(state.rv, well_cells);
342 const ADB& rs_perfcells =
subset(state.rs, well_cells);
343 const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
345 const int oilpos = pu.phase_pos[Oil];
346 const int gaspos = pu.phase_pos[Gas];
348 const ADB tmp_oil = (cmix_s[oilpos] - rv_perfcells * F_gas * cmix_s[gaspos]) / d;
349 volumeRatio += tmp_oil / b_perfcells[oilpos];
351 const ADB tmp_gas = (cmix_s[gaspos] - rs_perfcells * cmix_s[oilpos]) / d;
352 volumeRatio += tmp_gas / b_perfcells[gaspos];
355 if ((*active_)[Oil]) {
356 const int oilpos = pu.phase_pos[Oil];
357 volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
359 if ((*active_)[Gas]) {
360 const int gaspos = pu.phase_pos[Gas];
361 volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
367 ADB cqt_is = cqt_i/volumeRatio;
371 for (
int phase = 0; phase < np; ++phase) {
372 cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
376 aliveWells = Vector::Constant(nw, 1.0);
377 for (
int w = 0; w < nw; ++w) {
378 if (wbqt.value()[w] == 0) {
389 template <
class SolutionState,
class WellState>
391 StandardWellsSolvent::
392 computeWellConnectionPressures(
const SolutionState& state,
399 std::vector<double> b_perf;
400 std::vector<double> rsmax_perf;
401 std::vector<double> rvmax_perf;
402 std::vector<double> surf_dens_perf;
403 computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
405 const Vector pdepth = perf_cell_depth_;
406 const int nperf = wells().well_connpos[wells().number_of_wells];
407 const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
409 computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
417 template <
class ReservoirRes
idualQuant,
class SolutionState>
419 StandardWellsSolvent::
420 extractWellPerfProperties(
const SolutionState& state,
421 const std::vector<ReservoirResidualQuant>& rq,
422 std::vector<ADB>& mob_perfcells,
423 std::vector<ADB>& b_perfcells)
const
425 Base::extractWellPerfProperties(state, rq, mob_perfcells, b_perfcells);
428 const Opm::PhaseUsage& pu = fluid_->
phaseUsage();
429 int gas_pos = pu.phase_pos[Gas];
430 const std::vector<int>& well_cells = wellOps().well_cells;
431 const int nperf = well_cells.size();
437 mob_perfcells[gas_pos] +=
subset(rq[solvent_pos_].mob, well_cells);
440 const int nc = rq[solvent_pos_].mob.size();
443 const ADB& ss = state.solvent_saturation;
444 const ADB& sg = ((*active_)[ Gas ]
445 ? state.saturation[ pu.phase_pos[ Gas ] ]
448 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
449 ADB F_solvent =
subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
450 Vector ones = Vector::Constant(nperf,1.0);
452 b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
453 b_perfcells[gas_pos] += (F_solvent *
subset(rq[solvent_pos_].b, well_cells));
ADB bWat(const ADB &pw, const ADB &T, const Cells &cells) const
Water formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:446
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
V rsSat(const V &po, const Cells &cells) const
Bubble point curve for Rs as function of oil pressure.
ADB rvSat(const ADB &po, const Cells &cells) const
Condensation curve for Rv as function of oil pressure.
Definition: BlackoilPropsAdFromDeck.cpp:688
bool localWellsActive() const
return true if wells are available on this process
Definition: StandardWells_impl.hpp:148
V solventSurfaceDensity(const Cells &cells) const
Solvent surface density.
Definition: SolventPropsAdFromDeck.cpp:447
Vector & wellPerforationPressureDiffs()
Diff to bhp for each well perforation.
Definition: StandardWells_impl.hpp:204
ADB bOil(const ADB &po, const ADB &T, const ADB &rs, const std::vector< PhasePresence > &cond, const Cells &cells) const
Oil formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:493
ADB bSolvent(const ADB &pg, const Cells &cells) const
Solvent formation volume factor.
Definition: SolventPropsAdFromDeck.cpp:339
ADB bGas(const ADB &pg, const ADB &T, const ADB &rv, const std::vector< PhasePresence > &cond, const Cells &cells) const
Gas formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:566
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:319
V surfaceDensity(const int phaseIdx, const Cells &cells) const
Densities of stock components at surface conditions.
Definition: BlackoilPropsAdFromDeck.cpp:242
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
PhaseUsage phaseUsage() const
Definition: BlackoilPropsAdFromDeck.cpp:231