24 #ifndef OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED 25 #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED 27 #include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp> 29 #include <opm/autodiff/AutoDiffBlock.hpp> 30 #include <opm/autodiff/AutoDiffHelpers.hpp> 31 #include <opm/autodiff/GridHelpers.hpp> 32 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp> 33 #include <opm/autodiff/GeoProps.hpp> 34 #include <opm/autodiff/WellDensitySegmented.hpp> 36 #include <opm/core/grid.h> 37 #include <opm/core/linalg/LinearSolverInterface.hpp> 38 #include <opm/core/linalg/ParallelIstlInformation.hpp> 39 #include <opm/core/props/rock/RockCompressibility.hpp> 40 #include <opm/common/ErrorMacros.hpp> 41 #include <opm/common/Exceptions.hpp> 42 #include <opm/parser/eclipse/Units/Units.hpp> 43 #include <opm/core/well_controls.h> 44 #include <opm/core/utility/parameters/ParameterGroup.hpp> 59 int polymerPos(
const PU& pu)
61 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
63 for (
int phase = 0; phase < maxnp; ++phase) {
64 if (pu.phase_used[phase]) {
81 const RockCompressibility* rock_comp_props,
85 std::shared_ptr< const EclipseState > eclipse_state,
86 const bool has_disgas,
87 const bool has_vapoil,
88 const bool has_polymer,
89 const bool has_plyshlog,
90 const bool has_shrate,
91 const std::vector<double>& wells_rep_radius,
92 const std::vector<double>& wells_perf_length,
93 const std::vector<double>& wells_bore_diameter,
94 const bool terminal_output)
95 :
Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver, eclipse_state,
96 has_disgas, has_vapoil, terminal_output),
97 polymer_props_ad_(polymer_props_ad),
98 has_polymer_(has_polymer),
99 has_plyshlog_(has_plyshlog),
100 has_shrate_(has_shrate),
101 poly_pos_(detail::polymerPos(fluid.phaseUsage())),
102 wells_rep_radius_(wells_rep_radius),
103 wells_perf_length_(wells_perf_length),
104 wells_bore_diameter_(wells_bore_diameter)
107 if (!active_[Water]) {
108 OPM_THROW(std::logic_error,
"Polymer must solved in water!\n");
110 residual_.matbalscale.resize(fluid_.
numPhases() + 1, 1.1169);
114 Base::material_name_.push_back(
"Polymer");
121 template <
class Gr
id>
125 const ReservoirState& reservoir_state,
126 const WellState& well_state)
128 Base::prepareStep(timer, reservoir_state, well_state);
129 auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX );
132 cmax_ = Eigen::Map<const V>(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_));
138 template <
class Gr
id>
142 ReservoirState& reservoir_state,
145 computeCmax(reservoir_state);
152 template <
class Gr
id>
156 Base::makeConstantState(state);
157 state.concentration =
ADB::constant(state.concentration.value());
164 template <
class Gr
id>
166 BlackoilPolymerModel<Grid>::variableStateInitials(
const ReservoirState& x,
167 const WellState& xw)
const 169 std::vector<V> vars0 = Base::variableStateInitials(x, xw);
170 assert(
int(vars0.size()) == fluid_.numPhases() + 2);
174 const auto& concentration = x.getCellData( x.CONCENTRATION );
175 assert (not concentration.empty());
176 const int nc = concentration.size();
177 const V c = Eigen::Map<const V>(concentration.data() , nc);
179 auto concentration_pos = vars0.begin() + fluid_.numPhases();
180 assert(concentration_pos == vars0.end() - 2);
181 vars0.insert(concentration_pos, c);
190 template <
class Gr
id>
192 BlackoilPolymerModel<Grid>::variableStateIndices()
const 194 std::vector<int> ind = Base::variableStateIndices();
195 assert(ind.size() == 5);
199 ind[Concentration] = fluid_.numPhases();
210 template <
class Gr
id>
211 typename BlackoilPolymerModel<Grid>::SolutionState
212 BlackoilPolymerModel<Grid>::variableStateExtractVars(
const ReservoirState& x,
213 const std::vector<int>& indices,
214 std::vector<ADB>& vars)
const 216 SolutionState state = Base::variableStateExtractVars(x, indices, vars);
218 state.concentration = std::move(vars[indices[Concentration]]);
227 template <
class Gr
id>
229 BlackoilPolymerModel<Grid>::computeAccum(
const SolutionState& state,
232 Base::computeAccum(state, aix);
236 const ADB& press = state.pressure;
237 const std::vector<ADB>& sat = state.saturation;
238 const ADB& c = state.concentration;
239 const ADB pv_mult = poroMult(press);
240 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
242 const ADB cmax =
ADB::constant(cmax_, state.concentration.blockPattern());
243 const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
244 const double rho_rock = polymer_props_ad_.rockDensity();
245 const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], AutoDiffGrid::numCells(grid_));
246 const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
248 sd_.rq[poly_pos_].accum[aix] = pv_mult * sd_.rq[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
249 + pv_mult * rho_rock * (1. - phi) / phi * ads;
259 template <
class Gr
id>
260 void BlackoilPolymerModel<Grid>::computeCmax(ReservoirState& state)
262 auto& max_concentration = state.getCellData( state.CMAX );
263 const auto& concentration = state.getCellData( state.CONCENTRATION );
264 std::transform( max_concentration.begin() ,
265 max_concentration.end() ,
266 concentration.begin() ,
267 max_concentration.begin() ,
268 [](
double c_max ,
double c) {
return std::max( c_max , c ); });
275 template <
class Gr
id>
277 BlackoilPolymerModel<Grid>::
278 assembleMassBalanceEq(
const SolutionState& state)
288 computeAccum(state, 1);
292 const V transi =
subset(geo_.transmissibility(), ops_.internal_faces);
294 const std::vector<ADB> kr = computeRelPerm(state);
295 for (
int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
296 sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
302 std::vector<double> water_vel;
303 std::vector<double> visc_mult;
305 computeWaterShearVelocityFaces(transi, state.canonical_phase_pressures, state, water_vel, visc_mult);
306 if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) {
308 OPM_THROW(std::runtime_error,
" failed in calculating the shear-multiplier. ");
312 for (
int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
313 const std::vector<PhasePresence>& cond = phaseCondition();
314 sd_.rq[phaseIdx].mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
315 sd_.rq[phaseIdx].rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
316 computeMassFlux(phaseIdx, transi, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
318 residual_.material_balance_eq[ phaseIdx ] =
319 pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
320 + ops_.div*sd_.rq[phaseIdx].mflux;
328 if (active_[ Oil ] && active_[ Gas ]) {
329 const int po = fluid_.phaseUsage().phase_pos[ Oil ];
330 const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
332 const UpwindSelector<double> upwindOil(grid_, ops_,
333 sd_.rq[po].dh.value());
334 const ADB rs_face = upwindOil.select(state.rs);
336 const UpwindSelector<double> upwindGas(grid_, ops_,
337 sd_.rq[pg].dh.value());
338 const ADB rv_face = upwindGas.select(state.rv);
340 residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
341 residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
349 residual_.material_balance_eq[poly_pos_] = pvdt_ * (sd_.rq[poly_pos_].accum[1] - sd_.rq[poly_pos_].accum[0])
350 + ops_.div*sd_.rq[poly_pos_].mflux;
354 if (param_.update_equations_scaling_) {
355 updateEquationsScaling();
364 template <
class Gr
id>
365 void BlackoilPolymerModel<Grid>::updateEquationsScaling()
367 Base::updateEquationsScaling();
369 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
370 residual_.matbalscale[poly_pos_] = residual_.matbalscale[water_pos];
378 template <
class Gr
id>
379 void BlackoilPolymerModel<Grid>::addWellContributionToMassBalanceEq(
const std::vector<ADB>& cq_s,
380 const SolutionState& state,
384 Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
388 const ADB mc = computeMc(state);
389 const int nc = xw.polymerInflow().size();
390 const V polyin = Eigen::Map<const V>(xw.polymerInflow().data(), nc);
391 const int nperf = wells().well_connpos[wells().number_of_wells];
392 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
393 const V poly_in_perf =
subset(polyin, well_cells);
394 const V poly_mc_perf =
subset(mc.value(), well_cells);
395 const ADB& cq_s_water = cq_s[fluid_.phaseUsage().phase_pos[Water]];
396 Selector<double> injector_selector(cq_s_water.value());
397 const V poly_perf = injector_selector.select(poly_in_perf, poly_mc_perf);
398 const ADB cq_s_poly = cq_s_water * poly_perf;
399 residual_.material_balance_eq[poly_pos_] -=
superset(cq_s_poly, well_cells, nc);
408 template <
class Gr
id>
410 ReservoirState& reservoir_state,
411 WellState& well_state)
415 const int np = fluid_.numPhases();
416 const int nc = Opm::AutoDiffGrid::numCells(grid_);
417 const V zero = V::Zero(nc);
418 const int concentration_start = nc * np;
419 const V dc =
subset(dx,
Span(nc, 1, concentration_start));
422 V modified_dx = V::Zero(dx.size() - nc);
423 modified_dx.head(concentration_start) = dx.head(concentration_start);
424 const int tail_len = dx.size() - concentration_start - nc;
425 modified_dx.tail(tail_len) = dx.tail(tail_len);
428 Base::updateState(modified_dx, reservoir_state, well_state);
431 auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION );
433 const V c_old = Eigen::Map<const V>(concentration.data(), nc, 1);
434 const V c = (c_old - dc).max(zero);
435 std::copy(&c[0], &c[0] + nc, concentration.begin());
439 Base::updateState(dx, reservoir_state, well_state);
447 template <
class Gr
id>
454 const ADB& phasePressure,
455 const SolutionState& state)
457 Base::computeMassFlux(actph, transi, kr, mu, rho, phasePressure, state);
460 const int canonicalPhaseIdx = canph_[ actph ];
461 if (canonicalPhaseIdx == Water) {
463 const ADB tr_mult = transMult(state.pressure);
465 const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
466 const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.
value());
468 sd_.rq[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
470 const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.
value());
471 sd_.rq[poly_pos_].mob = tr_mult * krw_eff * state.concentration * inv_poly_eff_visc;
472 sd_.rq[poly_pos_].b = sd_.rq[actph].b;
473 sd_.rq[poly_pos_].dh = sd_.rq[actph].dh;
476 sd_.rq[poly_pos_].mflux = upwind.
select(sd_.rq[poly_pos_].b * sd_.rq[poly_pos_].mob) * (transi * sd_.rq[poly_pos_].dh);
478 sd_.rq[ actph ].mflux = upwind.select(sd_.rq[actph].b * sd_.rq[actph].mob) * (transi * sd_.rq[actph].dh);
482 V shear_mult_faces_v = Eigen::Map<V>(shear_mult_faces_.data(), shear_mult_faces_.size());
484 sd_.rq[poly_pos_].mflux = sd_.rq[poly_pos_].mflux / shear_mult_faces_adb;
485 sd_.rq[actph].mflux = sd_.rq[actph].mflux / shear_mult_faces_adb;
494 template <
class Gr
id>
497 WellState& well_state,
498 const bool initial_assembly)
502 SimulatorReport report;
507 wellModel().updateWellControls(well_state);
510 SolutionState state = variableState(reservoir_state, well_state);
512 if (initial_assembly) {
514 SolutionState state0 = state;
515 makeConstantState(state0);
518 computeAccum(state0, 0);
520 wellModel().computeWellConnectionPressures(state0, well_state);
533 assembleMassBalanceEq(state);
535 if ( ! wellsActive() ) {
539 std::vector<ADB> mob_perfcells;
540 std::vector<ADB> b_perfcells;
541 wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
545 computeInjectionMobility(state, mob_perfcells);
547 if (param_.solve_welleq_initially_ && initial_assembly) {
549 Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
553 std::vector<ADB> cq_s;
555 wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
558 std::vector<double> water_vel_wells;
559 std::vector<double> visc_mult_wells;
561 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
562 computeWaterShearVelocityWells(state, well_state, cq_s[water_pos], water_vel_wells, visc_mult_wells);
564 if ( !polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_) ) {
565 OPM_THROW(std::runtime_error,
" failed in calculating the shear factors for wells ");
569 V shear_mult_wells_v = Eigen::Map<V>(shear_mult_wells_.data(), shear_mult_wells_.size());
571 mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
572 wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
575 wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
576 wellModel().addWellFluxEq(cq_s, state, residual_);
577 addWellContributionToMassBalanceEq(cq_s, state, well_state);
578 wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
580 report.converged =
true;
587 template <
class Gr
id>
591 return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
600 const std::vector<ADB>& phasePressure,
const SolutionState& state,
601 std::vector<double>& water_vel, std::vector<double>& visc_mult)
604 const int phase = fluid_.phaseUsage().phase_pos[Water];
606 const int canonicalPhaseIdx = canph_[phase];
608 const std::vector<PhasePresence> cond = phaseCondition();
610 const ADB tr_mult = transMult(state.pressure);
611 const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond);
612 sd_.rq[phase].mob = tr_mult * sd_.rq[phase].kr / mu;
615 const ADB rho = fluidDensity(canonicalPhaseIdx, sd_.rq[phase].b, state.rs, state.rv);
616 const ADB rhoavg = ops_.caver * rho;
617 sd_.rq[ phase ].dh = ops_.ngrad * phasePressure[ canonicalPhaseIdx ] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
618 if (use_threshold_pressure_) {
619 applyThresholdPressures(sd_.rq[ phase ].dh);
622 const ADB& b = sd_.rq[ phase ].b;
623 const ADB& mob = sd_.rq[ phase ].mob;
624 const ADB& dh = sd_.rq[ phase ].dh;
628 ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
631 ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.
value());
632 sd_.rq[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
634 const V& polymer_conc = state.concentration.
value();
635 V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
636 V visc_mult_faces = upwind.select(visc_mult_cells);
638 size_t nface = visc_mult_faces.size();
639 visc_mult.resize(nface);
640 std::copy(visc_mult_faces.data(), visc_mult_faces.data() + nface, visc_mult.begin());
642 sd_.rq[ phase ].mflux = (transi * upwind.select(b * mob)) * dh;
645 const auto& b_faces_adb = upwind.select(b);
646 std::vector<double> b_faces(b_faces_adb.value().data(), b_faces_adb.value().data() + b_faces_adb.size());
648 const auto& internal_faces = ops_.internal_faces;
650 std::vector<double> internal_face_areas;
651 internal_face_areas.resize(internal_faces.size());
653 for (
int i = 0; i < internal_faces.size(); ++i) {
654 internal_face_areas[i] = grid_.face_areas[internal_faces[i]];
658 const ADB phiavg_adb = ops_.caver * phi;
660 std::vector<double> phiavg(phiavg_adb.
value().data(), phiavg_adb.
value().data() + phiavg_adb.
size());
662 water_vel.resize(nface);
663 std::copy(sd_.rq[0].mflux.value().data(), sd_.rq[0].mflux.value().data() + nface, water_vel.begin());
665 for (
size_t i = 0; i < nface; ++i) {
666 water_vel[i] = water_vel[i] / (b_faces[i] * phiavg[i] * internal_face_areas[i]);
673 const Opm::PhaseUsage pu = fluid_.phaseUsage();
674 const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
675 const ADB& sw_upwind_adb = upwind.select(sw);
676 std::vector<double> sw_upwind(sw_upwind_adb.
value().data(), sw_upwind_adb.
value().data() + sw_upwind_adb.
size());
679 std::vector<double> perm;
680 perm.resize(transi.size());
682 for (
int i = 0; i < transi.size(); ++i) {
683 perm[i] = transi[i] / internal_faces[i];
687 const ADB& krw_adb = upwind.select(krw_eff);
688 std::vector<double> krw_upwind(krw_adb.
value().data(), krw_adb.
value().data() + krw_adb.
size());
690 const double& shrate_const = polymer_props_ad_.shrate();
692 const double epsilon = std::numeric_limits<double>::epsilon();
696 for (
size_t i = 0; i < water_vel.size(); ++i) {
699 if (std::abs(water_vel[i]) < epsilon) {
703 water_vel[i] *= shrate_const * std::sqrt(phiavg[i] / (perm[i] * sw_upwind[i] * krw_upwind[i]));
716 std::vector<double>& water_vel_wells, std::vector<double>& visc_mult_wells)
718 if( ! wellsActive() ) return ;
720 const int nw = wells().number_of_wells;
721 const int nperf = wells().well_connpos[nw];
722 const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
724 water_vel_wells.resize(cq_sw.
size());
725 std::copy(cq_sw.
value().data(), cq_sw.
value().data() + cq_sw.
size(), water_vel_wells.begin());
727 const V& polymer_conc = state.concentration.value();
729 V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
730 V visc_mult_wells_v =
subset(visc_mult_cells, well_cells);
732 visc_mult_wells.resize(visc_mult_wells_v.size());
733 std::copy(visc_mult_wells_v.data(), visc_mult_wells_v.data() + visc_mult_wells_v.size(), visc_mult_wells.begin());
735 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
736 ADB b_perfcells =
subset(sd_.rq[water_pos].b, well_cells);
738 const ADB& p_perfcells =
subset(state.pressure, well_cells);
739 const V& cdp = wellModel().wellPerforationPressureDiffs();
740 const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
742 const ADB drawdown = p_perfcells - perfpressure;
745 V selectInjectingPerforations = V::Zero(nperf);
746 for (
int c = 0; c < nperf; ++c) {
747 if (drawdown.
value()[c] < 0) {
748 selectInjectingPerforations[c] = 1;
753 for (
size_t i = 0; i < well_cells.size(); ++i) {
754 if (xw.polymerInflow()[well_cells[i]] == 0. && selectInjectingPerforations[i] == 1) {
755 visc_mult_wells[i] = 1.;
757 if (selectInjectingPerforations[i] == 1) {
758 if (xw.polymerInflow()[well_cells[i]] == 0.) {
759 visc_mult_wells[i] = 1.;
762 const double c_perf = state.concentration.value()[well_cells[i]];
763 visc_mult_wells[i] = polymer_props_ad_.viscMult(c_perf);
769 const ADB phi_wells_adb =
subset(phi, well_cells);
771 std::vector<double> phi_wells(phi_wells_adb.
value().data(), phi_wells_adb.
value().data() + phi_wells_adb.
size());
773 std::vector<double> b_wells(b_perfcells.
value().data(), b_perfcells.
value().data() + b_perfcells.
size());
775 for (
size_t i = 0; i < water_vel_wells.size(); ++i) {
776 water_vel_wells[i] = b_wells[i] * water_vel_wells[i] / (phi_wells[i] * 2. * M_PI * wells_rep_radius_[i] * wells_perf_length_[i]);
783 const double& shrate_const = polymer_props_ad_.shrate();
784 for (
size_t i = 0; i < water_vel_wells.size(); ++i) {
785 water_vel_wells[i] = shrate_const * water_vel_wells[i] / wells_bore_diameter_[i];
800 std::vector<ADB>& mob_perfcells)
803 if (has_polymer_ && wellModel().localWellsActive()) {
804 const std::vector<int> well_cells = wellModel().wellOps().well_cells;
805 const int nperf = well_cells.size();
808 const ADB& p_perfcells =
subset(state.pressure, well_cells);
809 const V& cdp = wellModel().wellPerforationPressureDiffs();
810 const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
812 const ADB drawdown = p_perfcells - perfpressure;
815 const ADB c_perfcells =
subset(state.concentration, well_cells);
819 std::vector<int> polymer_inj_cells;
820 std::vector<int> other_well_cells;
822 polymer_inj_cells.reserve(nperf);
823 other_well_cells.reserve(nperf);
825 for (
int c = 0; c < nperf; ++c) {
827 if (drawdown.
value()[c] < 0.0 && c_perfcells.
value()[c] > 0.0) {
828 polymer_inj_cells.push_back(c);
830 other_well_cells.push_back(c);
835 if ( !polymer_inj_cells.empty() ) {
837 const int water_pos = fluid_.phaseUsage().phase_pos[Water];
838 const ADB mu_perfcells =
subset(sd_.rq[water_pos].mu, well_cells);
840 const ADB c_poly_inj_cells =
subset(c_perfcells, polymer_inj_cells);
841 const ADB mu_poly_inj_cells =
subset(mu_perfcells, polymer_inj_cells);
843 const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(c_poly_inj_cells, mu_poly_inj_cells.value());
844 const ADB fully_mixing_visc = polymer_props_ad_.viscMult(c_poly_inj_cells) * mu_poly_inj_cells;
847 ADB mob_polymer_inj =
subset(mob_perfcells[water_pos], polymer_inj_cells);
848 const ADB mob_others =
subset(mob_perfcells[water_pos], other_well_cells);
850 mob_polymer_inj = mob_polymer_inj / inv_wat_eff_visc / fully_mixing_visc;
851 mob_perfcells[water_pos] =
superset(mob_polymer_inj, polymer_inj_cells, nperf) +
852 superset(mob_others, other_well_cells, nperf);
859 #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
void computeWaterShearVelocityFaces(const V &transi, const std::vector< ADB > &phasePressure, const SolutionState &state, std::vector< double > &water_vel, std::vector< double > &visc_mult)
Computing the water velocity without shear-thinning for the cell faces.
Definition: BlackoilPolymerModel_impl.hpp:599
Definition: PolymerPropsAd.hpp:32
SimulatorReport assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilPolymerModel_impl.hpp:496
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
std::vector< ADB > select(const std::vector< ADB > &xc) const
Apply selector to multiple per-cell quantities.
Definition: AutoDiffHelpers.hpp:230
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
Class for handling the standard well model.
Definition: StandardWells.hpp:51
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilPolymerModel_impl.hpp:409
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilPolymerModel_impl.hpp:124
void afterStep(const SimulatorTimerInterface &timer, ReservoirState &reservoir_state, WellState &well_state)
Called once after each time step.
Definition: BlackoilPolymerModel_impl.hpp:141
std::vector< ADB > material_balance_eq
The material_balance_eq vector has one element for each active phase, each of which has size equal to...
Definition: LinearisedBlackoilResidual.hpp:54
Definition: GridHelpers.cpp:27
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
BlackoilPolymerModel(const typename Base::ModelParameters ¶m, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const PolymerPropsAd &polymer_props_ad, const StandardWells &well_model, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclipseState, const bool has_disgas, const bool has_vapoil, const bool has_polymer, const bool has_plyshlog, const bool has_shrate, const std::vector< double > &wells_rep_radius, const std::vector< double > &wells_perf_length, const std::vector< double > &wells_bore_diameter, const bool terminal_output)
Construct the model.
Definition: BlackoilPolymerModel_impl.hpp:77
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
A model implementation for three-phase black oil with polymer.
Definition: BlackoilPolymerModel.hpp:45
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
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
Definition: AutoDiffHelpers.hpp:623
Upwind selection in absence of counter-current flow (i.e., without effects of gravity and/or capillar...
Definition: AutoDiffHelpers.hpp:181
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
int numPhases() const
Definition: BlackoilPropsAdFromDeck.cpp:225
void computeWaterShearVelocityWells(const SolutionState &state, WellState &xw, const ADB &cq_sw, std::vector< double > &water_vel_wells, std::vector< double > &visc_mult_wells)
Computing the water velocity without shear-thinning for the well perforations based on the water flux...
Definition: BlackoilPolymerModel_impl.hpp:715
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415