20 #ifndef OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED
21 #define OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED
32 #include <opm/core/grid.h>
33 #include <opm/core/linalg/LinearSolverInterface.hpp>
34 #include <opm/core/linalg/ParallelIstlInformation.hpp>
35 #include <opm/core/props/rock/RockCompressibility.hpp>
36 #include <opm/common/ErrorMacros.hpp>
37 #include <opm/common/Exceptions.hpp>
38 #include <opm/core/utility/Units.hpp>
39 #include <opm/core/well_controls.h>
40 #include <opm/core/utility/parameters/ParameterGroup.hpp>
59 for (
int phase = 0; phase < maxnp; ++phase) {
60 if (pu.phase_used[phase]) {
77 const RockCompressibility* rock_comp_props,
79 const Wells* wells_arg,
81 const EclipseStateConstPtr eclState,
82 const bool has_disgas,
83 const bool has_vapoil,
84 const bool terminal_output,
85 const bool has_solvent)
86 :
Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
87 eclState, has_disgas, has_vapoil, terminal_output),
88 has_solvent_(has_solvent),
89 solvent_pos_(detail::
solventPos(fluid.phaseUsage())),
90 solvent_props_(solvent_props)
100 OPM_THROW(std::runtime_error,
"Solvent option only works with dead gas\n");
112 template <
class Gr
id>
116 Base::makeConstantState(state);
117 state.solvent_saturation =
ADB::constant(state.solvent_saturation.value());
124 template <
class Gr
id>
129 std::vector<V> vars0 = Base::variableStateInitials(x, xw);
130 assert(
int(vars0.size()) == fluid_.numPhases() + 2);
134 assert (not x.solvent_saturation().empty());
135 const int nc = x.solvent_saturation().size();
136 const V ss = Eigen::Map<const V>(&x.solvent_saturation()[0], nc);
138 auto solvent_pos = vars0.begin() + fluid_.numPhases();
139 assert(solvent_pos == vars0.end() - 2);
140 vars0.insert(solvent_pos, ss);
149 template <
class Gr
id>
153 std::vector<int> ind = Base::variableStateIndices();
154 assert(ind.size() == 5);
158 ind[Solvent] = fluid_.numPhases();
169 template <
class Gr
id>
172 const std::vector<int>& indices,
173 std::vector<ADB>& vars)
const
175 SolutionState state = Base::variableStateExtractVars(x, indices, vars);
177 state.solvent_saturation = std::move(vars[indices[Solvent]]);
178 if (active_[
Oil ]) {
180 const Opm::PhaseUsage pu = fluid_.phaseUsage();
181 state.saturation[pu.phase_pos[
Oil ]] -= state.solvent_saturation;
191 template <
class Gr
id>
196 Base::computeAccum(state, aix);
200 const ADB& press = state.pressure;
201 const ADB& ss = state.solvent_saturation;
202 const ADB pv_mult = poroMult(press);
203 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
205 const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[
Gas]];
206 rq_[solvent_pos_].b = solvent_props_.bSolvent(pg,cells_);
207 rq_[solvent_pos_].accum[aix] = pv_mult * rq_[solvent_pos_].b * ss;
215 template <
class Gr
id>
221 Base::assembleMassBalanceEq(state);
224 residual_.material_balance_eq[ solvent_pos_ ] =
225 pvdt_ * (rq_[solvent_pos_].accum[1] - rq_[solvent_pos_].accum[0])
226 + ops_.div*rq_[solvent_pos_].mflux;
231 template <
class Gr
id>
235 Base::updateEquationsScaling();
236 assert(
MaxNumPhases + 1 == residual_.matbalscale.size());
238 const ADB& temp_b = rq_[solvent_pos_].b;
240 residual_.matbalscale[solvent_pos_] = B.mean();
244 template <
class Gr
id>
253 Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
256 const int nperf = wells().well_connpos[wells().number_of_wells];
257 const int nc = Opm::AutoDiffGrid::numCells(grid_);
259 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
261 const ADB& ss = state.solvent_saturation;
262 const ADB& sg = (active_[
Gas ]
263 ? state.saturation[ pu.phase_pos[
Gas ] ]
266 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
268 ADB F_solvent =
subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
270 const int nw = wells().number_of_wells;
271 V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
273 V isProducer = V::Zero(nperf);
274 V ones = V::Constant(nperf,1.0);
275 for (
int w = 0; w < nw; ++w) {
276 if(wells().type[w] == PRODUCER) {
277 for (
int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
278 isProducer[perf] = 1;
283 const ADB& rs_perfcells =
subset(state.rs, well_cells);
284 int gas_pos = fluid_.phaseUsage().phase_pos[
Gas];
285 int oil_pos = fluid_.phaseUsage().phase_pos[
Oil];
288 assert(!has_vapoil_);
289 const ADB cq_s_solvent = (isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction) * (cq_s[gas_pos] - rs_perfcells * cq_s[oil_pos]);
293 residual_.material_balance_eq[solvent_pos_] -=
superset(cq_s_solvent, well_cells, nc);
297 residual_.material_balance_eq[gas_pos] +=
superset(cq_s_solvent, well_cells, nc);
302 template <
class Gr
id>
306 if( ! Base::localWellsActive() ) return ;
312 const int nperf = wells().well_connpos[wells().number_of_wells];
313 const int nw = wells().number_of_wells;
314 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
317 const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
318 V avg_press = perf_press*0;
319 for (
int w = 0; w < nw; ++w) {
320 for (
int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
321 const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
322 const double p_avg = (perf_press[perf] + p_above)/2;
323 avg_press[perf] = p_avg;
328 const ADB perf_temp =
subset(state.temperature, well_cells);
331 const PhaseUsage& pu = fluid_.phaseUsage();
333 DataBlock surf_dens(nperf, pu.num_phases);
334 for (
int phase = 0; phase < pu.num_phases; ++ phase) {
335 surf_dens.col(phase) = V::Constant(nperf, fluid_.surfaceDensity()[pu.phase_pos[phase]]);
342 std::vector<PhasePresence> perf_cond(nperf);
343 const std::vector<PhasePresence>& pc = phaseCondition();
344 for (
int perf = 0; perf < nperf; ++perf) {
345 perf_cond[perf] = pc[well_cells[perf]];
349 std::vector<double> rsmax_perf(nperf, 0.0);
350 std::vector<double> rvmax_perf(nperf, 0.0);
351 if (pu.phase_used[BlackoilPhases::Aqua]) {
352 const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
353 b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
355 assert(active_[
Oil]);
356 const V perf_so =
subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
357 if (pu.phase_used[BlackoilPhases::Liquid]) {
358 const ADB perf_rs =
subset(state.rs, well_cells);
359 const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
360 b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
361 const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
362 rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
364 if (pu.phase_used[BlackoilPhases::Vapour]) {
365 const ADB perf_rv =
subset(state.rv, well_cells);
366 V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
369 const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
371 const int nc = Opm::AutoDiffGrid::numCells(grid_);
374 const ADB& ss = state.solvent_saturation;
375 const ADB& sg = (active_[
Gas ]
376 ? state.saturation[ pu.phase_pos[
Gas ] ]
380 V F_solvent =
subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
382 V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
384 V isProducer = V::Zero(nperf);
385 V ones = V::Constant(nperf,1.0);
386 for (
int w = 0; w < nw; ++w) {
387 if(wells().type[w] == PRODUCER) {
388 for (
int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
389 isProducer[perf] = 1;
393 F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
395 bg = bg * (ones - F_solvent);
396 bg = bg + F_solvent * bs;
398 const V& rhog = surf_dens.col(pu.phase_pos[BlackoilPhases::Vapour]);
399 const V& rhos = solvent_props_.solventSurfaceDensity(well_cells);
400 surf_dens.col(pu.phase_pos[BlackoilPhases::Vapour]) = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
402 b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
404 const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
405 rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
409 std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
410 std::vector<double> surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases);
414 const V pdepth =
subset(depth, well_cells);
415 std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
421 std::vector<double> cd =
423 wells(), xw, fluid_.phaseUsage(),
424 b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
427 std::vector<double> cdp =
429 wells(), perf_depth, cd, grav);
432 Base::well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf);
433 Base::well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
441 template <
class Gr
id>
449 const int np = fluid_.numPhases();
450 const int nc = Opm::AutoDiffGrid::numCells(grid_);
451 const V zero = V::Zero(nc);
452 const int solvent_start = nc * np;
453 const V dss =
subset(dx,
Span(nc, 1, solvent_start));
456 V modified_dx = V::Zero(dx.size() - nc);
457 modified_dx.head(solvent_start) = dx.head(solvent_start);
458 const int tail_len = dx.size() - solvent_start - nc;
459 modified_dx.tail(tail_len) = dx.tail(tail_len);
462 Base::updateState(modified_dx, reservoir_state, well_state);
465 const V ss_old = Eigen::Map<const V>(&reservoir_state.solvent_saturation()[0], nc, 1);
466 const V ss = (ss_old - dss).max(zero);
467 std::copy(&ss[0], &ss[0] + nc, reservoir_state.solvent_saturation().begin());
470 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
471 const int oilpos = pu.phase_pos[
Oil ];
472 for (
int c = 0; c < nc; ++c) {
473 reservoir_state.saturation()[c*np + oilpos] = 1 - ss[c];
474 if (pu.phase_used[
Gas ]) {
475 const int gaspos = pu.phase_pos[
Gas ];
476 reservoir_state.saturation()[c*np + oilpos] -= reservoir_state.saturation()[c*np + gaspos];
478 if (pu.phase_used[
Water ]) {
479 const int waterpos = pu.phase_pos[
Water ];
480 reservoir_state.saturation()[c*np + oilpos] -= reservoir_state.saturation()[c*np + waterpos];
486 Base::updateState(dx, reservoir_state, well_state);
494 template <
class Gr
id>
499 const ADB& phasePressure,
502 Base::computeMassFlux(actph, transi, kr, phasePressure, state);
504 const int canonicalPhaseIdx = canph_[ actph ];
505 if (canonicalPhaseIdx ==
Gas) {
507 const int nc = Opm::UgGridHelpers::numCells(grid_);
509 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
511 const ADB& ss = state.solvent_saturation;
512 const ADB& sg = (active_[
Gas ]
513 ? state.saturation[ pu.phase_pos[
Gas ] ]
517 ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
518 V ones = V::Constant(nc, 1.0);
520 const ADB tr_mult = transMult(state.pressure);
521 const ADB mu = solvent_props_.muSolvent(phasePressure,cells_);
523 rq_[solvent_pos_].mob = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * tr_mult * kr / mu;
525 const ADB rho_solvent = solvent_props_.solventSurfaceDensity(cells_) * rq_[solvent_pos_].b;
526 const ADB rhoavg_solvent = ops_.caver * rho_solvent;
527 rq_[ solvent_pos_ ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg_solvent * (ops_.ngrad * geo_.z().matrix()));
531 rq_[solvent_pos_].mflux = upwind_solvent.
select(rq_[solvent_pos_].b * rq_[solvent_pos_].mob) * (transi * rq_[solvent_pos_].dh);
534 rq_[actph].mob = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * rq_[actph].mob;
536 const ADB& b = rq_[ actph ].b;
537 const ADB& mob = rq_[ actph ].mob;
538 const ADB& dh = rq_[ actph ].dh;
540 rq_[ actph ].mflux = upwind_gas.
select(b * mob) * (transi * dh);
546 template <
class Gr
id>
551 const int nc = numCells(grid_);
555 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
557 ? state.saturation[ pu.phase_pos[
Water ] ]
560 const ADB& so = (active_[
Oil ]
561 ? state.saturation[ pu.phase_pos[
Oil ] ]
564 const ADB& sg = (active_[
Gas ]
565 ? state.saturation[ pu.phase_pos[
Gas ] ]
569 const ADB& ss = state.solvent_saturation;
570 return fluid_.relperm(sw, so, sg+ss, cells_);
572 return fluid_.relperm(sw, so, sg, cells_);
578 template <
class Gr
id>
582 const bool initial_assembly)
589 updateWellControls(well_state);
592 SolutionState state = variableState(reservoir_state, well_state);
594 if (initial_assembly) {
597 makeConstantState(state0);
600 computeAccum(state0, 0);
601 computeWellConnectionPressures(state0, well_state);
605 assembleMassBalanceEq(state);
608 if ( ! wellsActive() ) {
614 const int np = wells().number_of_phases;
617 const int nw = wells().number_of_wells;
618 const int nperf = wells().well_connpos[nw];
619 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
621 std::vector<ADB> mob_perfcells(np,
ADB::null());
622 std::vector<ADB> b_perfcells(np,
ADB::null());
623 for (
int phase = 0; phase < np; ++phase) {
624 mob_perfcells[phase] =
subset(rq_[phase].mob, well_cells);
625 b_perfcells[phase] =
subset(rq_[phase].b, well_cells);
629 int gas_pos = fluid_.phaseUsage().phase_pos[
Gas];
635 mob_perfcells[gas_pos] +=
subset(rq_[solvent_pos_].mob, well_cells);
638 const int nc = Opm::AutoDiffGrid::numCells(grid_);
640 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
642 const ADB& ss = state.solvent_saturation;
643 const ADB& sg = (active_[
Gas ]
644 ? state.saturation[ pu.phase_pos[
Gas ] ]
648 ADB F_solvent =
subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
649 V ones = V::Constant(nperf,1.0);
650 b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
651 b_perfcells[gas_pos] += (F_solvent *
subset(rq_[solvent_pos_].b, well_cells));
653 if (param_.solve_welleq_initially_ && initial_assembly) {
655 Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
658 Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
659 Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
660 Base::addWellFluxEq(cq_s, state);
661 addWellContributionToMassBalanceEq(cq_s, state, well_state);
662 Base::addWellControlEq(state, well_state, aliveWells);
670 #endif // OPM_BLACKOILSOLVENT_IMPL_HEADER_INCLUDED
Definition: GeoProps.hpp:53
Base::SolutionState SolutionState
Definition: BlackoilSolventModel.hpp:101
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:287
std::vector< ADB > material_balance_eq
Definition: LinearisedBlackoilResidual.hpp:54
LinearisedBlackoilResidual residual_
Definition: BlackoilModelBase.hpp:279
Definition: BlackoilPropsAdInterface.hpp:38
std::vector< double > matbalscale
Definition: LinearisedBlackoilResidual.hpp:66
Eigen::Array< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:98
Definition: BlackoilModelEnums.hpp:47
Definition: AdditionalObjectDeleter.hpp:22
ModelTraits< BlackoilSolventModel< Grid > >::ModelParameters ModelParameters
Definition: BlackoilModelBase.hpp:108
void computeAccum(const SolutionState &state, const int aix)
Definition: BlackoilSolventModel_impl.hpp:193
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:103
std::vector< std::string > material_name_
Definition: BlackoilModelBase.hpp:286
std::vector< V > variableStateInitials(const ReservoirState &x, const WellState &xw) const
Definition: BlackoilSolventModel_impl.hpp:126
int solventPos(const PU &pu)
Definition: BlackoilSolventModel_impl.hpp:55
Definition: BlackoilModelEnums.hpp:32
void makeConstantState(SolutionState &state) const
Definition: BlackoilSolventModel_impl.hpp:114
Definition: AutoDiffHelpers.hpp:617
const int solvent_pos_
Definition: BlackoilSolventModel.hpp:107
virtual int numPhases() const =0
Definition: BlackoilModelEnums.hpp:46
Definition: GridHelpers.hpp:47
Selection. Choose first of two elements if selection basis element is nonnegative.
Definition: AutoDiffHelpers.hpp:359
static std::vector< double > computeConnectionPressureDelta(const Wells &wells, const std::vector< double > &z_perf, const std::vector< double > &dens_perf, const double gravity)
const BlackoilPropsAdInterface & fluid_
Definition: BlackoilModelBase.hpp:250
void assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Definition: BlackoilSolventModel_impl.hpp:580
void addWellContributionToMassBalanceEq(const std::vector< ADB > &cq_s, const SolutionState &state, WellState &xw)
Definition: BlackoilSolventModel_impl.hpp:245
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:314
static AutoDiffBlock constant(V &&val)
Definition: AutoDiffBlock.hpp:110
Definition: BlackoilModelEnums.hpp:30
SolutionState variableStateExtractVars(const ReservoirState &x, const std::vector< int > &indices, std::vector< ADB > &vars) const
Definition: BlackoilSolventModel_impl.hpp:171
Base::WellState WellState
Definition: BlackoilSolventModel.hpp:46
void computeMassFlux(const int actph, const V &transi, const ADB &kr, const ADB &p, const SolutionState &state)
Definition: BlackoilSolventModel_impl.hpp:496
Definition: BlackoilModelEnums.hpp:29
Eigen::Array< double, Eigen::Dynamic, 1 > cellCentroidsZToEigen(const UnstructuredGrid &grid)
Get the z coordinates of the cell centroids of a grid.
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
std::vector< ADB > select(const std::vector< ADB > &xc) const
Apply selector to multiple per-cell quantities.
Definition: AutoDiffHelpers.hpp:225
const bool has_vapoil_
Definition: BlackoilModelBase.hpp:264
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Definition: BlackoilSolventModel_impl.hpp:442
void computeWellConnectionPressures(const SolutionState &state, const WellState &xw)
Definition: BlackoilSolventModel_impl.hpp:303
void updateEquationsScaling()
Definition: BlackoilSolventModel_impl.hpp:233
Definition: AutoDiffHelpers.hpp:176
std::vector< ReservoirResidualQuant > rq_
Definition: BlackoilModelBase.hpp:271
double getGravity(const double *g, const int dim)
Definition: BlackoilModelBase_impl.hpp:698
BlackoilSolventModel(const typename Base::ModelParameters ¶m, const Grid &grid, const BlackoilPropsAdInterface &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const SolventPropsAdFromDeck &solvent_props, const Wells *wells, const NewtonIterationBlackoilInterface &linsolver, const EclipseStateConstPtr eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output, const bool has_solvent)
Definition: BlackoilSolventModel_impl.hpp:73
Definition: SolventPropsAdFromDeck.hpp:38
void assembleMassBalanceEq(const SolutionState &state)
Definition: BlackoilSolventModel_impl.hpp:218
const bool has_solvent_
Definition: BlackoilSolventModel.hpp:106
Base::ReservoirState ReservoirState
Definition: BlackoilSolventModel.hpp:45
std::vector< ADB > computeRelPerm(const SolutionState &state) const
Definition: BlackoilSolventModel_impl.hpp:548
std::vector< int > variableStateIndices() const
Definition: BlackoilSolventModel_impl.hpp:151
ADB::V V
Definition: BlackoilModelBase.hpp:103
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:436
Base::DataBlock DataBlock
Definition: BlackoilSolventModel.hpp:102
static std::vector< double > computeConnectionDensities(const Wells &wells, const WellStateFullyImplicitBlackoil &wstate, const PhaseUsage &phase_usage, const std::vector< double > &b_perf, const std::vector< double > &rsmax_perf, const std::vector< double > &rvmax_perf, const std::vector< double > &surf_dens_perf)
Definition: BlackoilModelEnums.hpp:31