24 #ifndef OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
25 #define OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
39 #include <opm/core/grid.h>
40 #include <opm/core/linalg/LinearSolverInterface.hpp>
41 #include <opm/core/linalg/ParallelIstlInformation.hpp>
42 #include <opm/core/props/rock/RockCompressibility.hpp>
43 #include <opm/common/ErrorMacros.hpp>
44 #include <opm/common/Exceptions.hpp>
45 #include <opm/core/utility/Units.hpp>
46 #include <opm/core/well_controls.h>
47 #include <opm/core/utility/parameters/ParameterGroup.hpp>
48 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
59 #define OPM_AD_DUMP(foo) \
61 std::cout << "==========================================\n" \
63 << collapseJacs(foo) << std::endl; \
66 #define OPM_AD_DUMPVAL(foo) \
68 std::cout << "==========================================\n" \
70 << foo.value() << std::endl; \
73 #define OPM_AD_DISKVAL(foo) \
75 std::ofstream os(#foo); \
77 os << foo.value() << std::endl; \
86 typedef Eigen::Array<double,
99 std::vector<int> all_cells(nc);
101 for (
int c = 0; c < nc; ++c) { all_cells[c] = c; }
113 std::vector<bool> active(maxnp,
false);
115 for (
int p = 0; p < pu.MaxNumPhases; ++p) {
116 active[ p ] = pu.phase_used[ p ] != 0;
129 std::vector<int> act2can(maxnp, -1);
131 for (
int phase = 0; phase < maxnp; ++phase) {
132 if (pu.phase_used[ phase ]) {
133 act2can[ pu.phase_pos[ phase ] ] = phase;
143 template <
class Gr
id,
class Implementation>
149 const RockCompressibility* rock_comp_props,
150 const Wells* wells_arg,
152 Opm::EclipseStateConstPtr eclState,
153 const bool has_disgas,
154 const bool has_vapoil,
155 const bool terminal_output)
159 , rock_comp_props_(rock_comp_props)
161 , vfp_properties_(eclState->getTableManager()->getVFPInjTables(), eclState->getTableManager()->getVFPProdTables())
162 , linsolver_ (linsolver)
166 , ops_ (grid, eclState)
168 , has_disgas_(has_disgas)
169 , has_vapoil_(has_vapoil)
171 , use_threshold_pressure_(false)
172 , rq_ (fluid.numPhases())
173 , phaseCondition_(AutoDiffGrid::numCells(grid))
174 , isRs_(V::Zero(AutoDiffGrid::numCells(grid)))
175 , isRv_(V::Zero(AutoDiffGrid::numCells(grid)))
176 , isSg_(V::Zero(AutoDiffGrid::numCells(grid)))
180 { 1.1169, 1.0031, 0.0031 }} )
181 , terminal_output_ (terminal_output)
182 , material_name_{
"Water",
"Oil",
"Gas" }
184 assert(numMaterials() == 3);
186 if ( linsolver_.parallelInformation().type() ==
typeid(ParallelISTLInformation) )
188 const ParallelISTLInformation& info =
189 boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation());
190 if ( terminal_output_ ) {
192 terminal_output_ = (info.communicator().rank()==0);
194 int local_number_of_wells = wells_ ? wells_->number_of_wells : 0;
195 int global_number_of_wells = info.communicator().sum(local_number_of_wells);
196 wells_active_ = ( wells_ && global_number_of_wells > 0 );
199 wells_active_ = ( wells_ && wells_->number_of_wells > 0 );
202 wells_active_ = ( wells_ && wells_->number_of_wells > 0 );
210 template <
class Gr
id,
class Implementation>
217 pvdt_ = geo_.poreVolume() / dt;
219 updatePrimalVariableFromState(reservoir_state);
227 template <
class Gr
id,
class Implementation>
241 template <
class Gr
id,
class Implementation>
246 return residual_.sizeNonLinear();
253 template <
class Gr
id,
class Implementation>
258 return linsolver_.iterations();
265 template <
class Gr
id,
class Implementation>
270 return terminal_output_;
277 template <
class Gr
id,
class Implementation>
282 return fluid_.numPhases();
289 template <
class Gr
id,
class Implementation>
294 return material_name_.size();
301 template <
class Gr
id,
class Implementation>
306 assert(material_index < numMaterials());
307 return material_name_[material_index];
314 template <
class Gr
id,
class Implementation>
319 const int num_faces = AutoDiffGrid::numFaces(grid_);
320 if (
int(threshold_pressures_by_face.size()) != num_faces) {
321 OPM_THROW(std::runtime_error,
"Illegal size of threshold_pressures_by_face input, must be equal to number of faces.");
323 use_threshold_pressure_ =
true;
325 const int num_ifaces = ops_.internal_faces.size();
326 threshold_pressures_by_interior_face_.resize(num_ifaces);
327 for (
int ii = 0; ii < num_ifaces; ++ii) {
328 threshold_pressures_by_interior_face_[ii] = threshold_pressures_by_face[ops_.internal_faces[ii]];
336 template <
class Gr
id,
class Implementation>
338 : accum(2, ADB::null())
339 , mflux( ADB::null())
350 template <
class Gr
id,
class Implementation>
358 w2p = Eigen::SparseMatrix<double>(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells);
359 p2w = Eigen::SparseMatrix<double>(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]);
361 const int nw = wells->number_of_wells;
362 const int*
const wpos = wells->well_connpos;
364 typedef Eigen::Triplet<double> Tri;
366 std::vector<Tri> scatter, gather;
367 scatter.reserve(wpos[nw]);
368 gather .reserve(wpos[nw]);
370 for (
int w = 0, i = 0; w < nw; ++w) {
371 for (; i < wpos[ w + 1 ]; ++i) {
372 scatter.push_back(Tri(i, w, 1.0));
373 gather .push_back(Tri(w, i, 1.0));
377 w2p.setFromTriplets(scatter.begin(), scatter.end());
378 p2w.setFromTriplets(gather .begin(), gather .end());
386 template <
class Gr
id,
class Implementation>
395 state.temperature =
ADB::constant(state.temperature.value());
398 const int num_phases = state.saturation.
size();
399 for (
int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) {
400 state.saturation[phaseIdx] =
ADB::constant(state.saturation[phaseIdx].value());
406 ADB& pp = state.canonical_phase_pressures[canphase];
415 template <
class Gr
id,
class Implementation>
420 std::vector<V> vars0 =
asImpl().variableStateInitials(x, xw);
429 template <
class Gr
id,
class Implementation>
436 const int np = x.numPhases();
438 std::vector<V> vars0;
441 vars0.reserve(np + 1);
451 template <
class Gr
id,
class Implementation>
456 const int nc = numCells(
grid_);
457 const int np = x.numPhases();
459 assert (not x.pressure().empty());
460 const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
464 assert (not x.saturation().empty());
465 const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
470 const V sw = s.col(pu.phase_pos[ Water ]);
477 const V sg = s.col(pu.phase_pos[ Gas ]);
478 const V rs = Eigen::Map<const V>(& x.gasoilratio()[0], x.gasoilratio().size());
479 const V rv = Eigen::Map<const V>(& x.rv()[0], x.rv().size());
481 vars0.push_back(xvar);
489 template <
class Gr
id,
class Implementation>
498 const int nw =
wells().number_of_wells;
499 const int np =
wells().number_of_phases;
502 const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
503 const V qs = Eigen::Map<const V>(wrates.data(), nw*np);
507 assert (not xw.bhp().empty());
508 const V
bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
509 vars0.push_back(bhp);
514 vars0.push_back(
V());
515 vars0.push_back(
V());
523 template <
class Gr
id,
class Implementation>
528 std::vector<int> indices(5, -1);
532 indices[
Sw] = next++;
535 indices[
Xvar] = next++;
537 indices[
Qs] = next++;
538 indices[
Bhp] = next++;
546 template <
class Gr
id,
class Implementation>
552 std::vector<int> indices(5, -1);
554 indices[
Qs] = next++;
555 indices[
Bhp] = next++;
564 template <
class Gr
id,
class Implementation>
567 const std::vector<int>& indices,
568 std::vector<ADB>& vars)
const
571 const int nc = Opm::AutoDiffGrid::numCells(
grid_);
577 state.pressure = std::move(vars[indices[
Pressure]]);
580 const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
588 state.saturation[pu.phase_pos[
Water ]] = std::move(vars[indices[
Sw]]);
589 const ADB& sw = state.saturation[pu.phase_pos[
Water ]];
596 const ADB& xvar = vars[indices[
Xvar]];
597 ADB& sg = state.saturation[ pu.phase_pos[
Gas ] ];
604 ? state.saturation[ pu.phase_pos[
Water ] ]
606 state.canonical_phase_pressures =
computePressures(state.pressure, sw, so, sg);
607 const ADB rsSat =
fluidRsSat(state.canonical_phase_pressures[ Oil ], so ,
cells_);
613 const ADB rvSat =
fluidRvSat(state.canonical_phase_pressures[ Gas ], so ,
cells_);
624 state.saturation[pu.phase_pos[
Oil ]] = std::move(so);
636 template <
class Gr
id,
class Implementation>
639 std::vector<ADB>& vars,
643 state.qs = std::move(vars[indices[
Qs]]);
646 state.bhp = std::move(vars[indices[
Bhp]]);
653 template <
class Gr
id,
class Implementation>
660 const ADB& press = state.pressure;
661 const ADB& temp = state.temperature;
662 const std::vector<ADB>& sat = state.saturation;
663 const ADB& rs = state.rs;
664 const ADB& rv = state.rv;
668 const ADB pv_mult =
poroMult(press);
671 for (
int phase = 0; phase < maxnp; ++phase) {
673 const int pos = pu.phase_pos[ phase ];
674 rq_[pos].b =
fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
675 rq_[pos].accum[aix] = pv_mult *
rq_[pos].b * sat[pos];
683 const int po = pu.phase_pos[
Oil ];
684 const int pg = pu.phase_pos[
Gas ];
688 const ADB accum_gas_copy =
rq_[pg].accum[aix];
690 rq_[pg].accum[aix] += state.rs *
rq_[po].accum[aix];
691 rq_[po].accum[aix] += state.rv * accum_gas_copy;
702 for (
int dd = 0; dd < dim - 1; ++dd) {
703 assert(g[dd] == 0.0);
713 template <
class Gr
id,
class Implementation>
723 const int nperf =
wells().well_connpos[
wells().number_of_wells];
724 const int nw =
wells().number_of_wells;
725 const std::vector<int> well_cells(
wells().well_cells,
wells().well_cells + nperf);
728 const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
729 V avg_press = perf_press*0;
730 for (
int w = 0; w < nw; ++w) {
731 for (
int perf =
wells().well_connpos[w]; perf <
wells().well_connpos[w+1]; ++perf) {
732 const double p_above = perf ==
wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
733 const double p_avg = (perf_press[perf] + p_above)/2;
734 avg_press[perf] = p_avg;
739 const ADB perf_temp =
subset(state.temperature, well_cells);
745 std::vector<PhasePresence> perf_cond(nperf);
747 for (
int perf = 0; perf < nperf; ++perf) {
748 perf_cond[perf] = pc[well_cells[perf]];
751 DataBlock b(nperf, pu.num_phases);
752 std::vector<double> rsmax_perf(nperf, 0.0);
753 std::vector<double> rvmax_perf(nperf, 0.0);
754 if (pu.phase_used[BlackoilPhases::Aqua]) {
756 b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
759 const V perf_so =
subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
760 if (pu.phase_used[BlackoilPhases::Liquid]) {
761 const ADB perf_rs =
subset(state.rs, well_cells);
762 const V bo =
fluid_.
bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).
value();
763 b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
764 const V rssat =
fluidRsSat(avg_press, perf_so, well_cells);
765 rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
767 if (pu.phase_used[BlackoilPhases::Vapour]) {
768 const ADB perf_rv =
subset(state.rv, well_cells);
769 const V bg =
fluid_.
bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).
value();
770 b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
771 const V rvsat =
fluidRvSat(avg_press, perf_so, well_cells);
772 rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
775 std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
778 const V pdepth =
subset(depth, well_cells);
779 std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
782 DataBlock surf_dens(nperf, pu.num_phases);
783 for (
int phase = 0; phase < pu.num_phases; ++ phase) {
787 std::vector<double> surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases);
793 std::vector<double> cd =
796 b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
799 std::vector<double> cdp =
801 wells(), perf_depth, cd, grav);
812 template <
class Gr
id,
class Implementation>
817 const bool initial_assembly)
827 asImpl().makeConstantState(state0);
838 if (initial_assembly) {
841 asImpl().makeConstantState(state0);
844 asImpl().computeAccum(state0, 0);
858 asImpl().assembleMassBalanceEq(state);
867 const int np =
wells().number_of_phases;
870 const int nw =
wells().number_of_wells;
871 const int nperf =
wells().well_connpos[nw];
872 const std::vector<int> well_cells(
wells().well_cells,
wells().well_cells + nperf);
874 std::vector<ADB> mob_perfcells(np,
ADB::null());
875 std::vector<ADB> b_perfcells(np,
ADB::null());
876 for (
int phase = 0; phase < np; ++phase) {
877 mob_perfcells[phase] =
subset(
rq_[phase].mob, well_cells);
878 b_perfcells[phase] =
subset(
rq_[phase].b, well_cells);
880 if (
param_.solve_welleq_initially_ && initial_assembly) {
882 solveWellEq(mob_perfcells, b_perfcells, state, well_state);
885 asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
886 asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
887 asImpl().addWellFluxEq(cq_s, state);
888 asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
896 template <
class Gr
id,
class Implementation>
907 asImpl().computeAccum(state, 1);
913 V trans_all(transi.size() + trans_nnc.size());
914 trans_all << transi, trans_nnc;
916 const std::vector<ADB> kr =
asImpl().computeRelPerm(state);
917 #pragma omp parallel for schedule(static)
919 asImpl().computeMassFlux(phaseIdx, trans_all, kr[
canph_[phaseIdx]], state.canonical_phase_pressures[
canph_[phaseIdx]], state);
922 pvdt_ * (
rq_[phaseIdx].accum[1] -
rq_[phaseIdx].accum[0])
937 const ADB rs_face = upwindOil.
select(state.rs);
941 const ADB rv_face = upwindGas.
select(state.rv);
951 if (
param_.update_equations_scaling_) {
952 asImpl().updateEquationsScaling();
961 template <
class Gr
id,
class Implementation>
969 const int pos = pu.phase_pos[idx];
970 const ADB& temp_b =
rq_[pos].b;
971 B = 1. / temp_b.
value();
981 template <
class Gr
id,
class Implementation>
988 const int nc = Opm::AutoDiffGrid::numCells(
grid_);
989 const int nw =
wells().number_of_wells;
990 const int nperf =
wells().well_connpos[nw];
991 const int np =
wells().number_of_phases;
992 const std::vector<int> well_cells(
wells().well_cells,
wells().well_cells + nperf);
993 for (
int phase = 0; phase < np; ++phase) {
1001 template <
class Gr
id,
class Implementation>
1004 const std::vector<ADB>& mob_perfcells,
1005 const std::vector<ADB>& b_perfcells,
1007 std::vector<ADB>& cq_s)
1011 const int np =
wells().number_of_phases;
1012 const int nw =
wells().number_of_wells;
1013 const int nperf =
wells().well_connpos[nw];
1015 V Tw = Eigen::Map<const V>(
wells().WI, nperf);
1016 const std::vector<int> well_cells(
wells().well_cells,
wells().well_cells + nperf);
1021 const ADB& p_perfcells =
subset(state.pressure, well_cells);
1022 const ADB& rv_perfcells =
subset(state.rv, well_cells);
1023 const ADB& rs_perfcells =
subset(state.rs, well_cells);
1026 const ADB perfpressure = (
wops_.
w2p * state.bhp) + cdp;
1029 const ADB drawdown = p_perfcells - perfpressure;
1035 V selectInjectingPerforations = V::Zero(nperf);
1037 V selectProducingPerforations = V::Zero(nperf);
1038 for (
int c = 0; c < nperf; ++c){
1039 if (drawdown.
value()[c] < 0)
1040 selectInjectingPerforations[c] = 1;
1042 selectProducingPerforations[c] = 1;
1046 const V numInjectingPerforations = (
wops_.
p2w *
ADB::constant(selectInjectingPerforations)).value();
1047 const V numProducingPerforations = (
wops_.
p2w *
ADB::constant(selectProducingPerforations)).value();
1048 for (
int w = 0; w < nw; ++w) {
1049 if (!
wells().allow_cf[w]) {
1050 for (
int perf =
wells().well_connpos[w] ; perf <
wells().well_connpos[w+1]; ++perf) {
1055 if (
wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
1056 selectProducingPerforations[perf] = 0.0;
1057 }
else if (
wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
1058 selectInjectingPerforations[perf] = 0.0;
1066 std::vector<ADB> cq_ps(np,
ADB::null());
1067 for (
int phase = 0; phase < np; ++phase) {
1068 const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
1069 cq_ps[phase] = b_perfcells[phase] * cq_p;
1072 const int oilpos = pu.phase_pos[
Oil];
1073 const int gaspos = pu.phase_pos[
Gas];
1074 const ADB cq_psOil = cq_ps[oilpos];
1075 const ADB cq_psGas = cq_ps[gaspos];
1076 cq_ps[gaspos] += rs_perfcells * cq_psOil;
1077 cq_ps[oilpos] += rv_perfcells * cq_psGas;
1082 ADB total_mob = mob_perfcells[0];
1083 for (
int phase = 1; phase < np; ++phase) {
1084 total_mob += mob_perfcells[phase];
1087 const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
1094 const DataBlock compi = Eigen::Map<const DataBlock>(
wells().comp_frac, nw, np);
1097 for (
int phase = 0; phase < np; ++phase) {
1098 const ADB& q_ps =
wops_.
p2w * cq_ps[phase];
1099 const ADB& q_s =
subset(state.qs,
Span(nw, 1, phase*nw));
1101 const int pos = pu.phase_pos[phase];
1102 wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,
ADB::constant(V::Zero(nw)))) - q_ps;
1107 std::vector<ADB> cmix_s(np,
ADB::null());
1108 for (
int phase = 0; phase < np; ++phase) {
1109 const int pos = pu.phase_pos[phase];
1110 cmix_s[phase] =
wops_.
w2p * notDeadWells_selector.select(
ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
1115 const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
1116 for (
int phase = 0; phase < np; ++phase) {
1117 ADB tmp = cmix_s[phase];
1120 const int gaspos = pu.phase_pos[
Gas];
1121 tmp = tmp - rv_perfcells * cmix_s[gaspos] / d;
1124 const int oilpos = pu.phase_pos[
Oil];
1125 tmp = tmp - rs_perfcells * cmix_s[oilpos] / d;
1127 volumeRatio += tmp / b_perfcells[phase];
1131 ADB cqt_is = cqt_i/volumeRatio;
1134 for (
int phase = 0; phase < np; ++phase) {
1135 cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
1139 aliveWells = V::Constant(nw, 1.0);
1140 for (
int w = 0; w < nw; ++w) {
1141 if (wbqt.
value()[w] == 0) {
1142 aliveWells[w] = 0.0;
1151 template <
class Gr
id,
class Implementation>
1157 const int np =
wells().number_of_phases;
1158 const int nw =
wells().number_of_wells;
1159 const int nperf =
wells().well_connpos[nw];
1160 V cq =
superset(cq_s[0].value(),
Span(nperf, np, 0), nperf*np);
1161 for (
int phase = 1; phase < np; ++phase) {
1162 cq +=
superset(cq_s[phase].value(),
Span(nperf, np, phase), nperf*np);
1164 xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np);
1168 const V perfpressure = (
wops_.
w2p * state.bhp.value().matrix()).array() + cdp;
1169 xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf);
1176 template <
class Gr
id,
class Implementation>
1180 const int np =
wells().number_of_phases;
1181 const int nw =
wells().number_of_wells;
1183 for (
int phase = 0; phase < np; ++phase) {
1200 const int num_phases,
1201 const double* distr)
1204 for (
int phase = 0; phase < num_phases; ++phase) {
1207 rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase];
1214 const std::vector<double>& thp,
1215 const std::vector<double>& well_phase_flow_rate,
1217 const int num_phases,
1218 const WellType& well_type,
1219 const WellControls* wc,
1220 const int ctrl_index)
1222 const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
1223 const double target = well_controls_iget_target(wc, ctrl_index);
1224 const double* distr = well_controls_iget_distr(wc, ctrl_index);
1226 bool broken =
false;
1228 switch (well_type) {
1231 switch (ctrl_type) {
1233 broken = bhp[well] > target;
1237 broken = thp[well] > target;
1240 case RESERVOIR_RATE:
1243 well, num_phases, distr) > target;
1251 switch (ctrl_type) {
1253 broken = bhp[well] < target;
1257 broken = thp[well] < target;
1260 case RESERVOIR_RATE:
1266 well, num_phases, distr) < target;
1273 OPM_THROW(std::logic_error,
"Can only handle INJECTOR and PRODUCER wells.");
1292 const ADB::V& well_perforation_densities,
const double gravity) {
1293 if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
1301 const double well_ref_depth = wells.depth_ref[w];
1302 const double dh = vfp_ref_depth - well_ref_depth;
1303 const int perf = wells.well_connpos[w];
1304 const double rho = well_perforation_densities[perf];
1305 const double dp = rho*gravity*dh;
1312 const ADB::V& well_perforation_densities,
const double gravity) {
1313 const int nw = wells.number_of_wells;
1314 ADB::V retval = ADB::V::Zero(nw);
1316 #pragma omp parallel for schedule(static)
1317 for (
int i=0; i<nw; ++i) {
1327 template <
class Gr
id,
class Implementation>
1338 const int nw =
wells().number_of_wells;
1340 for (
int w = 0; w < nw; ++w) {
1341 const WellControls* wc =
wells().ctrls[w];
1343 const int nwc = well_controls_get_num(wc);
1346 for (
int c=0; c < nwc; ++c) {
1347 const WellControlType ctrl_type = well_controls_iget_type(wc, c);
1349 if (ctrl_type == THP) {
1359 template <
class Gr
id,
class Implementation>
1364 std::string modestring[4] = {
"BHP",
"THP",
"RESERVOIR_RATE",
"SURFACE_RATE" };
1367 const int np =
wells().number_of_phases;
1368 const int nw =
wells().number_of_wells;
1370 #pragma omp parallel for schedule(dynamic)
1371 for (
int w = 0; w < nw; ++w) {
1372 const WellControls* wc =
wells().ctrls[w];
1376 int current = xw.currentControls()[w];
1380 const int nwc = well_controls_get_num(wc);
1382 for (; ctrl_index < nwc; ++ctrl_index) {
1383 if (ctrl_index == current) {
1390 xw.bhp(), xw.thp(), xw.wellRates(),
1391 w, np,
wells().type[w], wc, ctrl_index)) {
1396 if (ctrl_index != nwc) {
1400 std::cout <<
"Switching control mode for well " <<
wells().name[w]
1401 <<
" from " << modestring[well_controls_iget_type(wc, current)]
1402 <<
" to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
1404 xw.currentControls()[w] = ctrl_index;
1405 current = xw.currentControls()[w];
1413 const double target = well_controls_iget_target(wc, current);
1414 const double* distr = well_controls_iget_distr(wc, current);
1415 switch (well_controls_iget_type(wc, current)) {
1417 xw.bhp()[w] = target;
1422 double liquid = 0.0;
1423 double vapour = 0.0;
1426 aqua = xw.wellRates()[w*np + pu.phase_pos[
Water ] ];
1429 liquid = xw.wellRates()[w*np + pu.phase_pos[
Oil ] ];
1432 vapour = xw.wellRates()[w*np + pu.phase_pos[
Gas ] ];
1435 const int vfp = well_controls_iget_vfp(wc, current);
1436 const double& thp = well_controls_iget_target(wc, current);
1437 const double& alq = well_controls_iget_alq(wc, current);
1440 const WellType& well_type =
wells().type[w];
1442 if (well_type == INJECTOR) {
1449 else if (well_type == PRODUCER) {
1457 OPM_THROW(std::logic_error,
"Expected PRODUCER or INJECTOR type of well");
1462 case RESERVOIR_RATE:
1470 for (
int phase = 0; phase < np; ++phase) {
1471 if (distr[phase] > 0.0) {
1472 xw.wellRates()[np*w + phase] = target * distr[phase];
1485 template <
class Gr
id,
class Implementation>
1487 const std::vector<ADB>& b_perfcells,
1492 const int np =
wells().number_of_phases;
1496 asImpl().makeConstantState(state0);
1498 std::vector<ADB> mob_perfcells_const(np,
ADB::null());
1499 std::vector<ADB> b_perfcells_const(np,
ADB::null());
1500 for (
int phase = 0; phase < np; ++phase) {
1501 mob_perfcells_const[phase] =
ADB::constant(mob_perfcells[phase].value());
1502 b_perfcells_const[phase] =
ADB::constant(b_perfcells[phase].value());
1509 std::vector<V> vars0;
1516 asImpl().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
1517 asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
1518 asImpl().addWellFluxEq(cq_s, wellSolutionState);
1529 std::vector<ADB> eqs;
1534 const std::vector<M>& Jn = total_residual.
derivative();
1535 typedef Eigen::SparseMatrix<double> Sp;
1537 Jn[0].toSparse(Jn0);
1538 const Eigen::SparseLU< Sp > solver(Jn0);
1540 const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
1541 assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1)));
1549 std::cout <<
"well converged iter: " << it << std::endl;
1551 const int nw =
wells().number_of_wells;
1555 ADB::V new_bhp = Eigen::Map<ADB::V>(well_state.bhp().data(), nw);
1558 std::vector<ADB::M> old_derivs = state.bhp.derivative();
1559 state.bhp =
ADB::function(std::move(new_bhp), std::move(old_derivs));
1565 const DataBlock wrates = Eigen::Map<const DataBlock>(well_state.wellRates().data(), nw, np).transpose();
1566 ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np);
1567 std::vector<ADB::M> old_derivs = state.qs.derivative();
1568 state.qs =
ADB::function(std::move(new_qs), std::move(old_derivs));
1570 asImpl().computeWellConnectionPressures(state, well_state);
1579 template <
class Gr
id,
class Implementation>
1582 const V& aliveWells)
1586 const int np =
wells().number_of_phases;
1587 const int nw =
wells().number_of_wells;
1594 aqua +=
subset(state.qs,
Span(nw, 1, BlackoilPhases::Aqua*nw));
1597 liquid +=
subset(state.qs,
Span(nw, 1, BlackoilPhases::Liquid*nw));
1600 vapour +=
subset(state.qs,
Span(nw, 1, BlackoilPhases::Vapour*nw));
1604 std::vector<int> inj_table_id(nw, -1);
1605 std::vector<int> prod_table_id(nw, -1);
1606 ADB::V thp_inj_target_v = ADB::V::Zero(nw);
1607 ADB::V thp_prod_target_v = ADB::V::Zero(nw);
1608 ADB::V alq_v = ADB::V::Zero(nw);
1611 ADB::V rho_v = ADB::V::Zero(nw);
1612 ADB::V vfp_ref_depth_v = ADB::V::Zero(nw);
1615 ADB::V bhp_targets = ADB::V::Zero(nw);
1616 ADB::V rate_targets = ADB::V::Zero(nw);
1617 Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
1620 std::vector<int> bhp_elems;
1621 std::vector<int> thp_inj_elems;
1622 std::vector<int> thp_prod_elems;
1623 std::vector<int> rate_elems;
1627 for (
int w = 0; w < nw; ++w) {
1628 auto wc =
wells().ctrls[w];
1633 const int current = xw.currentControls()[w];
1635 switch (well_controls_iget_type(wc, current)) {
1638 bhp_elems.push_back(w);
1639 bhp_targets(w) = well_controls_iget_target(wc, current);
1640 rate_targets(w) = -1e100;
1646 const int perf =
wells().well_connpos[w];
1649 const int table_id = well_controls_iget_vfp(wc, current);
1650 const double target = well_controls_iget_target(wc, current);
1652 const WellType& well_type =
wells().type[w];
1653 if (well_type == INJECTOR) {
1654 inj_table_id[w] = table_id;
1655 thp_inj_target_v[w] = target;
1660 thp_inj_elems.push_back(w);
1662 else if (well_type == PRODUCER) {
1663 prod_table_id[w] = table_id;
1664 thp_prod_target_v[w] = target;
1665 alq_v[w] = well_controls_iget_alq(wc, current);
1669 thp_prod_elems.push_back(w);
1672 OPM_THROW(std::logic_error,
"Expected INJECTOR or PRODUCER type well");
1674 bhp_targets(w) = -1e100;
1675 rate_targets(w) = -1e100;
1679 case RESERVOIR_RATE:
1682 rate_elems.push_back(w);
1687 const double*
const distr =
1688 well_controls_iget_distr(wc, current);
1690 for (
int p = 0; p < np; ++p) {
1691 rate_distr.insert(w, p*nw + w) = distr[p];
1694 bhp_targets(w) = -1.0e100;
1695 rate_targets(w) = well_controls_iget_target(wc, current);
1703 const ADB thp_prod_target =
ADB::constant(thp_prod_target_v);
1706 const ADB bhp_from_thp_prod =
vfp_properties_.
getProd()->
bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq);
1712 const ADB dp_inj =
superset(
subset(dp, thp_inj_elems), thp_inj_elems, nw);
1713 const ADB dp_prod =
superset(
subset(dp, thp_prod_elems), thp_prod_elems, nw);
1716 const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj;
1717 const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod;
1718 const ADB bhp_residual = state.bhp - bhp_targets;
1719 const ADB rate_residual = rate_distr * state.qs - rate_targets;
1723 superset(
subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) +
1724 superset(
subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) +
1730 Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
1731 for (
int w = 0; w < nw; ++w) {
1732 for (
int phase = 0; phase < np; ++phase) {
1733 rate_summer.insert(w, phase*nw + w) = 1.0;
1745 template <
class Gr
id,
class Implementation>
1763 double infinityNorm(
const ADB& a,
const boost::any& pinfo = boost::any() )
1765 static_cast<void>(pinfo);
1767 if ( pinfo.type() ==
typeid(ParallelISTLInformation) )
1769 const ParallelISTLInformation& real_info =
1770 boost::any_cast<
const ParallelISTLInformation&>(pinfo);
1772 real_info.computeReduction(a.
value(), Reduction::makeGlobalMaxFunctor<double>(), result);
1778 if( a.
value().size() > 0 ) {
1779 return a.
value().matrix().lpNorm<Eigen::Infinity> ();
1793 static_cast<void>(pinfo);
1795 if( a.
value().size() > 0 ) {
1796 result = a.
value().matrix().lpNorm<Eigen::Infinity> ();
1799 if ( pinfo.type() ==
typeid(ParallelISTLInformation) )
1801 const ParallelISTLInformation& real_info =
1802 boost::any_cast<
const ParallelISTLInformation&>(pinfo);
1803 result = real_info.communicator().max(result);
1815 template <
class Gr
id,
class Implementation>
1822 const int nc = numCells(
grid_);
1825 assert(null.size() == 0);
1826 const V zero = V::Zero(nc);
1832 varstart += dsw.size();
1835 varstart += dxvar.size();
1838 const V dwells =
subset(dx,
Span((np+1)*nw, 1, varstart));
1839 varstart += dwells.size();
1841 assert(varstart == dx.size());
1844 const double dpmaxrel =
dpMaxRel();
1845 const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
1846 const V absdpmax = dpmaxrel*p_old.abs();
1847 const V dp_limited =
sign(dp) * dp.abs().min(absdpmax);
1848 const V p = (p_old - dp_limited).max(zero);
1849 std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
1854 const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
1855 const double dsmax =
dsMax();
1865 maxVal = dsw.abs().max(maxVal);
1872 maxVal = dsg.abs().max(maxVal);
1876 maxVal = dso.abs().max(maxVal);
1878 V step = dsmax/maxVal;
1879 step = step.min(1.);
1882 const int pos = pu.phase_pos[
Water ];
1883 const V sw_old = s_old.col(pos);
1884 sw = sw_old - step * dsw;
1888 const int pos = pu.phase_pos[
Gas ];
1889 const V sg_old = s_old.col(pos);
1890 sg = sg_old - step * dsg;
1893 const int pos = pu.phase_pos[
Oil ];
1894 const V so_old = s_old.col(pos);
1895 so = so_old - step * dso;
1900 for (
int c = 0; c < nc; ++c) {
1902 sw[c] = sw[c] / (1-sg[c]);
1903 so[c] = so[c] / (1-sg[c]);
1910 for (
int c = 0; c < nc; ++c) {
1912 sw[c] = sw[c] / (1-so[c]);
1913 sg[c] = sg[c] / (1-so[c]);
1919 for (
int c = 0; c < nc; ++c) {
1921 so[c] = so[c] / (1-sw[c]);
1922 sg[c] = sg[c] / (1-sw[c]);
1933 for (
int c = 0; c < nc; ++c) {
1934 reservoir_state.saturation()[c*np + pu.phase_pos[
Water ]] = sw[c];
1937 for (
int c = 0; c < nc; ++c) {
1938 reservoir_state.saturation()[c*np + pu.phase_pos[
Gas ]] = sg[c];
1942 const int pos = pu.phase_pos[
Oil ];
1943 for (
int c = 0; c < nc; ++c) {
1944 reservoir_state.saturation()[c*np + pos] = so[c];
1949 const double drmaxrel =
drMaxRel();
1952 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
1953 const V drs =
isRs_ * dxvar;
1954 const V drs_limited =
sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
1955 rs = rs_old - drs_limited;
1959 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
1960 const V drv =
isRv_ * dxvar;
1961 const V drv_limited =
sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
1962 rv = rv_old - drv_limited;
1966 const double epsilon =
std::sqrt(std::numeric_limits<double>::epsilon());
1967 auto watOnly = sw > (1 - epsilon);
1976 auto hasGas = (sg > 0 &&
isRs_ == 0);
1979 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
1980 auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
1981 auto useSg = watOnly || hasGas || gasVaporized;
1982 for (
int c = 0; c < nc; ++c) {
1998 const V rvSat0 =
fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]),
cells_);
2002 auto hasOil = (so > 0 &&
isRv_ == 0);
2005 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
2006 auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
2007 auto useSg = watOnly || hasOil || oilCondensed;
2008 for (
int c = 0; c < nc; ++c) {
2020 std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
2024 std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
2038 template <
class Gr
id,
class Implementation>
2046 const int np =
wells().number_of_phases;
2047 const int nw =
wells().number_of_wells;
2051 const V dqs =
subset(dwells,
Span(np*nw, 1, varstart));
2052 varstart += dqs.size();
2053 const V dbhp =
subset(dwells,
Span(nw, 1, varstart));
2054 varstart += dbhp.size();
2055 assert(varstart == dwells.size());
2056 const double dpmaxrel =
dpMaxRel();
2063 const DataBlock wwr = Eigen::Map<const DataBlock>(dqs.data(), np, nw).transpose();
2064 const V dwr = Eigen::Map<const V>(wwr.data(), nw*np);
2065 const V wr_old = Eigen::Map<const V>(&well_state.wellRates()[0], nw*np);
2066 const V wr = wr_old - dwr;
2067 std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
2070 const V bhp_old = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
2071 const V dbhp_limited =
sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
2072 const V
bhp = bhp_old - dbhp_limited;
2073 std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
2081 #pragma omp parallel for schedule(static)
2082 for (
int w=0; w<nw; ++w) {
2083 const WellControls* wc =
wells().ctrls[w];
2084 const int nwc = well_controls_get_num(wc);
2088 for (
int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
2089 if (well_controls_iget_type(wc, ctrl_index) == THP) {
2091 double liquid = 0.0;
2092 double vapour = 0.0;
2095 aqua = wr[w*np + pu.phase_pos[
Water ] ];
2098 liquid = wr[w*np + pu.phase_pos[
Oil ] ];
2101 vapour = wr[w*np + pu.phase_pos[
Gas ] ];
2104 double alq = well_controls_iget_alq(wc, ctrl_index);
2105 int table_id = well_controls_iget_vfp(wc, ctrl_index);
2107 const WellType& well_type =
wells().type[w];
2108 if (well_type == INJECTOR) {
2115 else if (well_type == PRODUCER) {
2123 OPM_THROW(std::logic_error,
"Expected INJECTOR or PRODUCER well");
2138 template <
class Gr
id,
class Implementation>
2143 const int nc = numCells(
grid_);
2149 ? state.saturation[ pu.phase_pos[
Water ] ]
2153 ? state.saturation[ pu.phase_pos[
Oil ] ]
2157 ? state.saturation[ pu.phase_pos[
Gas ] ]
2167 template <
class Gr
id,
class Implementation>
2173 const ADB& sg)
const
2179 if (phaseIdx == BlackoilPhases::Liquid)
2181 pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
2190 if (phaseIdx == BlackoilPhases::Aqua) {
2191 pressure[phaseIdx] = po - pressure[phaseIdx];
2193 pressure[phaseIdx] += po;
2204 template <
class Gr
id,
class Implementation>
2216 return cp[
Gas].value() + po;
2223 template <
class Gr
id,
class Implementation>
2228 const ADB& phasePressure,
2232 const int canonicalPhaseIdx =
canph_[ actph ];
2234 const ADB tr_mult =
transMult(state.pressure);
2235 const ADB mu =
fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv, cond);
2236 rq_[ actph ].mob = tr_mult * kr / mu;
2239 const ADB rho =
fluidDensity(canonicalPhaseIdx,
rq_[actph].b, state.rs, state.rv);
2247 const ADB& b =
rq_[ actph ].b;
2248 const ADB& mob =
rq_[ actph ].mob;
2249 const ADB& dh =
rq_[ actph ].dh;
2251 rq_[ actph ].mflux = upwind.select(b * mob) * (transi * dh);
2258 template <
class Gr
id,
class Implementation>
2276 const M keep_high_potential(high_potential.matrix().asDiagonal());
2283 dp = keep_high_potential * (dp - threshold_modification);
2290 template <
class Gr
id,
class Implementation>
2294 std::vector<double> residualNorms;
2299 for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
2302 if (!std::isfinite(massBalanceResid)) {
2303 OPM_THROW(Opm::NumericalProblem,
2304 "Encountered a non-finite residual");
2306 residualNorms.push_back(massBalanceResid);
2312 if (!std::isfinite(wellFluxResid)) {
2313 OPM_THROW(Opm::NumericalProblem,
2314 "Encountered a non-finite residual");
2316 residualNorms.push_back(wellFluxResid);
2320 if (!std::isfinite(wellResid)) {
2321 OPM_THROW(Opm::NumericalProblem,
2322 "Encountered a non-finite residual");
2324 residualNorms.push_back(wellResid);
2326 return residualNorms;
2333 template <
class Gr
id,
class Implementation>
2336 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
2337 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
2338 std::vector<double>& R_sum,
2339 std::vector<double>& maxCoeff,
2340 std::vector<double>& B_avg,
2341 std::vector<double>& maxNormWell,
2345 const int np =
asImpl().numPhases();
2346 const int nm =
asImpl().numMaterials();
2352 const ParallelISTLInformation& info =
2356 std::vector<int> v(nc, 1);
2357 auto nc_and_pv = std::tuple<int, double>(0, 0.0);
2358 auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
2359 Opm::Reduction::makeGlobalSumFunctor<double>());
2360 auto nc_and_pv_containers = std::make_tuple(v,
geo_.
poreVolume());
2361 info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
2363 for (
int idx = 0; idx < nm; ++idx )
2365 auto values = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
2366 auto containers = std::make_tuple(B.col(idx),
2369 auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
2370 Opm::Reduction::makeGlobalMaxFunctor<double>(),
2371 Opm::Reduction::makeGlobalSumFunctor<double>());
2372 info.computeReduction(containers, operators, values);
2373 B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
2374 maxCoeff[idx] = std::get<1>(values);
2375 R_sum[idx] = std::get<2>(values);
2378 maxNormWell[idx] = 0.0;
2379 for (
int w = 0; w < nw; ++w ) {
2384 info.communicator().max(maxNormWell.data(), np);
2386 return std::get<1>(nc_and_pv);
2392 maxCoeff.resize(nm);
2394 maxNormWell.resize(np);
2395 for (
int idx = 0; idx < nm; ++idx )
2397 B_avg[idx] = B.col(idx).sum()/nc;
2398 maxCoeff[idx] = tempV.col(idx).maxCoeff();
2399 R_sum[idx] = R.col(idx).sum();
2403 maxNormWell[idx] = 0.0;
2404 for (
int w = 0; w < nw; ++w ) {
2418 template <
class Gr
id,
class Implementation>
2422 const double tol_mb =
param_.tolerance_mb_;
2423 const double tol_cnv =
param_.tolerance_cnv_;
2424 const double tol_wells =
param_.tolerance_wells_;
2426 const int nc = Opm::AutoDiffGrid::numCells(
grid_);
2428 const int np =
asImpl().numPhases();
2429 const int nm =
asImpl().numMaterials();
2430 assert(
int(
rq_.size()) == nm);
2434 std::vector<double> R_sum(nm);
2435 std::vector<double> B_avg(nm);
2436 std::vector<double> maxCoeff(nm);
2437 std::vector<double> maxNormWell(np);
2438 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
2439 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
2440 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
2442 for (
int idx = 0; idx < nm; ++idx )
2444 const ADB& tempB =
rq_[idx].b;
2445 B.col(idx) = 1./tempB.
value();
2447 tempV.col(idx) = R.col(idx).abs()/pv;
2451 R_sum, maxCoeff, B_avg, maxNormWell,
2454 std::vector<double> CNV(nm);
2455 std::vector<double> mass_balance_residual(nm);
2456 std::vector<double> well_flux_residual(np);
2458 bool converged_MB =
true;
2459 bool converged_CNV =
true;
2460 bool converged_Well =
true;
2462 for (
int idx = 0; idx < nm; ++idx )
2464 CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
2465 mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
2466 converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
2467 converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
2472 well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
2473 converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
2479 converged_Well = converged_Well && (residualWell < Opm::unit::barsa);
2480 const bool converged = converged_MB && converged_CNV && converged_Well;
2485 for (
int idx = 0; idx < nm; ++idx) {
2486 if (std::isnan(mass_balance_residual[idx])
2487 || std::isnan(CNV[idx])
2488 || (idx < np && std::isnan(well_flux_residual[idx]))) {
2489 OPM_THROW(Opm::NumericalProblem,
"NaN residual for phase " <<
materialName(idx));
2494 OPM_THROW(Opm::NumericalProblem,
"Too large residual for phase " <<
materialName(idx));
2497 if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) {
2498 OPM_THROW(Opm::NumericalProblem,
"NaN or too large residual for well control equation");
2504 if (iteration == 0) {
2505 std::cout <<
"\nIter";
2506 for (
int idx = 0; idx < nm; ++idx) {
2507 std::cout <<
" MB(" <<
materialName(idx).substr(0, 3) <<
") ";
2509 for (
int idx = 0; idx < nm; ++idx) {
2510 std::cout <<
" CNV(" <<
materialName(idx).substr(0, 1) <<
") ";
2512 for (
int idx = 0; idx < np; ++idx) {
2513 std::cout <<
" W-FLUX(" <<
materialName(idx).substr(0, 1) <<
")";
2517 const std::streamsize oprec = std::cout.precision(3);
2518 const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific);
2519 std::cout << std::setw(4) << iteration;
2520 for (
int idx = 0; idx < nm; ++idx) {
2521 std::cout << std::setw(11) << mass_balance_residual[idx];
2523 for (
int idx = 0; idx < nm; ++idx) {
2524 std::cout << std::setw(11) << CNV[idx];
2526 for (
int idx = 0; idx < np; ++idx) {
2527 std::cout << std::setw(11) << well_flux_residual[idx];
2529 std::cout << std::endl;
2530 std::cout.precision(oprec);
2531 std::cout.flags(oflags);
2540 template <
class Gr
id,
class Implementation>
2544 const double tol_wells =
param_.tolerance_wells_;
2546 const int nc = Opm::AutoDiffGrid::numCells(
grid_);
2548 const int np =
asImpl().numPhases();
2549 const int nm =
asImpl().numMaterials();
2552 std::vector<double> R_sum(nm);
2553 std::vector<double> B_avg(nm);
2554 std::vector<double> maxCoeff(nm);
2555 std::vector<double> maxNormWell(np);
2556 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
2557 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
2558 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
2559 for (
int idx = 0; idx < nm; ++idx )
2561 const ADB& tempB =
rq_[idx].b;
2562 B.col(idx) = 1./tempB.
value();
2564 tempV.col(idx) = R.col(idx).abs()/pv;
2569 std::vector<double> well_flux_residual(np);
2570 bool converged_Well =
true;
2572 for (
int idx = 0; idx < np; ++idx )
2574 well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
2575 converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
2580 converged_Well = converged_Well && (residualWell < Opm::unit::barsa);
2581 const bool converged = converged_Well;
2584 for (
int idx = 0; idx < np; ++idx) {
2585 if (std::isnan(well_flux_residual[idx])) {
2586 OPM_THROW(Opm::NumericalProblem,
"NaN residual for phase " <<
materialName(idx));
2589 OPM_THROW(Opm::NumericalProblem,
"Too large residual for phase " <<
materialName(idx));
2596 if (iteration == 0) {
2597 std::cout <<
"\nIter";
2598 for (
int idx = 0; idx < np; ++idx) {
2599 std::cout <<
" W-FLUX(" <<
materialName(idx).substr(0, 1) <<
")";
2603 const std::streamsize oprec = std::cout.precision(3);
2604 const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific);
2605 std::cout << std::setw(4) << iteration;
2606 for (
int idx = 0; idx < np; ++idx) {
2607 std::cout << std::setw(11) << well_flux_residual[idx];
2609 std::cout << std::endl;
2610 std::cout.precision(oprec);
2611 std::cout.flags(oflags);
2620 template <
class Gr
id,
class Implementation>
2627 const std::vector<PhasePresence>& cond)
const
2637 OPM_THROW(std::runtime_error,
"Unknown phase index " << phase);
2645 template <
class Gr
id,
class Implementation>
2652 const std::vector<PhasePresence>& cond)
const
2662 OPM_THROW(std::runtime_error,
"Unknown phase index " << phase);
2670 template <
class Gr
id,
class Implementation>
2675 const ADB& rv)
const
2678 ADB rho = rhos[phase] * b;
2681 rho += rhos[
Gas] * rs * b;
2685 rho += rhos[
Oil] * rv * b;
2694 template <
class Gr
id,
class Implementation>
2698 const std::vector<int>& cells)
const
2707 template <
class Gr
id,
class Implementation>
2711 const std::vector<int>& cells)
const
2720 template <
class Gr
id,
class Implementation>
2724 const std::vector<int>& cells)
const
2733 template <
class Gr
id,
class Implementation>
2737 const std::vector<int>& cells)
const
2746 template <
class Gr
id,
class Implementation>
2750 const int n = p.
size();
2754 #pragma omp parallel for schedule(static)
2755 for (
int i = 0; i < n; ++i) {
2759 ADB::M dpm_diag(dpm.matrix().asDiagonal());
2761 std::vector<ADB::M> jacs(num_blocks);
2762 #pragma omp parallel for schedule(dynamic)
2763 for (
int block = 0; block < num_blocks; ++block) {
2776 template <
class Gr
id,
class Implementation>
2780 const int n = p.
size();
2784 #pragma omp parallel for schedule(static)
2785 for (
int i = 0; i < n; ++i) {
2789 ADB::M dtm_diag(dtm.matrix().asDiagonal());
2791 std::vector<ADB::M> jacs(num_blocks);
2792 #pragma omp parallel for schedule(dynamic)
2793 for (
int block = 0; block < num_blocks; ++block) {
2806 template <
class Gr
id,
class Implementation>
2811 const int nc = numCells(
grid_);
2812 const int np = state.numPhases();
2815 const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
2818 const V so = s.col(pu.phase_pos[
Oil ]);
2819 const V sg = s.col(pu.phase_pos[ Gas ]);
2821 for (V::Index c = 0, e = sg.size(); c != e; ++c) {
2831 const V so = s.col(pu.phase_pos[
Oil ]);
2834 for (V::Index c = 0, e = so.size(); c != e; ++c) {
2846 template <
class Gr
id,
class Implementation>
2851 const int nc = numCells(
grid_);
2852 const int np = state.numPhases();
2855 const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
2864 const V sg = s.col(pu.phase_pos[ Gas ]);
2865 const V so = s.col(pu.phase_pos[
Oil ]);
2866 const V sw = s.col(pu.phase_pos[
Water ]);
2868 const double epsilon =
std::sqrt(std::numeric_limits<double>::epsilon());
2869 auto watOnly = sw > (1 - epsilon);
2870 auto hasOil = so > 0;
2871 auto hasGas = sg > 0;
2876 for (V::Index c = 0, e = sg.size(); c != e; ++c) {
2884 for (V::Index c = 0, e = so.size(); c != e; ++c) {
2896 template <
class Gr
id,
class Implementation>
2901 OPM_THROW(std::logic_error,
"updatePhaseCondFromPrimarVariable() logic requires active gas phase.");
2904 isRs_ = V::Zero(nc);
2905 isRv_ = V::Zero(nc);
2906 isSg_ = V::Zero(nc);
2907 for (
int c = 0; c < nc; ++c) {
2925 OPM_THROW(std::logic_error,
"Unknown primary variable enum value in cell " << c <<
": " <<
primalVariable_[c]);
2935 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
VFPEvaluation bhp(const VFPProdTable *table, const double &aqua, const double &liquid, const double &vapour, const double &thp, const double &alq)
Definition: VFPHelpers.hpp:517
Definition: GeoProps.hpp:53
std::vector< double > computeResidualNorms() const
Compute the residual norms of the mass balance for each phase, the well flux, and the well equation...
Definition: BlackoilModelBase_impl.hpp:2292
void makeConstantState(SolutionState &state) const
Definition: BlackoilModelBase_impl.hpp:388
Definition: BlackoilModelEnums.hpp:38
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
Definition: BlackoilModelEnums.hpp:43
double drMaxRel() const
Definition: BlackoilModelBase.hpp:518
std::vector< ADB > material_balance_eq
Definition: LinearisedBlackoilResidual.hpp:54
SolutionState variableStateExtractVars(const ReservoirState &x, const std::vector< int > &indices, std::vector< ADB > &vars) const
Definition: BlackoilModelBase_impl.hpp:566
virtual ADB muGas(const ADB &pg, const ADB &T, const ADB &rv, const std::vector< PhasePresence > &cond, const Cells &cells) const =0
LinearisedBlackoilResidual residual_
Definition: BlackoilModelBase.hpp:279
virtual ADB muWat(const ADB &pw, const ADB &T, const Cells &cells) const =0
Definition: BlackoilPropsAdInterface.hpp:38
double rateToCompare(const std::vector< double > &well_phase_flow_rate, const int well, const int num_phases, const double *distr)
Definition: BlackoilModelBase_impl.hpp:1198
std::vector< double > matbalscale
Definition: LinearisedBlackoilResidual.hpp:66
std::vector< int > primalVariable_
Definition: BlackoilModelBase.hpp:284
const DerivedGeology & geo_
Definition: BlackoilModelBase.hpp:251
M caver
Extract for each face the average of its adjacent cells' values.
Definition: AutoDiffHelpers.hpp:54
AutoDiffBlock< double > ADB
Definition: BlackoilModelBase_impl.hpp:83
Eigen::Array< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:98
V fluidRvSat(const V &p, const V &so, const std::vector< int > &cells) const
Definition: BlackoilModelBase_impl.hpp:2722
const VFPInjProperties * getInj() const
Definition: VFPProperties.hpp:61
static std::vector< AutoDiffBlock > variables(const std::vector< V > &initial_values)
Definition: AutoDiffBlock.hpp:202
Definition: BlackoilModelEnums.hpp:47
const Vector & transmissibility() const
Definition: GeoProps.hpp:177
Definition: AdditionalObjectDeleter.hpp:22
V isRs_
Definition: BlackoilModelBase.hpp:273
bool getConvergence(const double dt, const int iteration)
Definition: BlackoilModelBase_impl.hpp:2420
void afterStep(const double dt, ReservoirState &reservoir_state, WellState &well_state)
Definition: BlackoilModelBase_impl.hpp:230
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Definition: BlackoilModelBase_impl.hpp:1816
double infinityNorm(const ADB &a, const boost::any &pinfo=boost::any())
Compute the L-infinity norm of a vector This function is not suitable to compute on the well equatio...
Definition: BlackoilModelBase_impl.hpp:1763
ModelTraits< Implementation >::ModelParameters ModelParameters
Definition: BlackoilModelBase.hpp:108
double thp(int table_id, const double &aqua, const double &liquid, const double &vapour, const double &bhp, const double &alq) const
const VFPInjTable * getTable(const int table_id) const
const Wells & wells() const
Definition: BlackoilModelBase.hpp:309
const std::vector< bool > active_
Definition: BlackoilModelBase.hpp:257
double convergenceReduction(const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > &B, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > &tempV, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > &R, std::vector< double > &R_sum, std::vector< double > &maxCoeff, std::vector< double > &B_avg, std::vector< double > &maxNormWell, int nc, int nw) const
Compute the reduction within the convergence check.
Definition: BlackoilModelBase_impl.hpp:2335
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:103
void variableStateExtractWellsVars(const std::vector< int > &indices, std::vector< ADB > &vars, SolutionState &state) const
Definition: BlackoilModelBase_impl.hpp:638
const RockCompressibility * rock_comp_props_
Definition: BlackoilModelBase.hpp:252
int numBlocks() const
Number of Jacobian blocks.
Definition: AutoDiffBlock.hpp:419
Definition: BlackoilModelEnums.hpp:32
void computeWellConnectionPressures(const SolutionState &state, const WellState &xw)
Definition: BlackoilModelBase_impl.hpp:714
double computeHydrostaticCorrection(const Wells &wells, const int w, double vfp_ref_depth, const ADB::V &well_perforation_densities, const double gravity)
Definition: BlackoilModelBase_impl.hpp:1291
ModelTraits< Implementation >::ReservoirState ReservoirState
Definition: BlackoilModelBase.hpp:106
void updateWellState(const V &dwells, WellState &well_state)
Definition: BlackoilModelBase_impl.hpp:2040
const std::vector< int > canph_
Definition: BlackoilModelBase.hpp:259
bool wellsActive() const
Definition: BlackoilModelBase.hpp:305
Definition: AutoDiffHelpers.hpp:617
virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual &residual) const =0
Implementation & asImpl()
Definition: BlackoilModelBase.hpp:292
virtual ADB rsSat(const ADB &po, const Cells &cells) const =0
virtual std::vector< ADB > capPress(const ADB &sw, const ADB &so, const ADB &sg, const Cells &cells) const =0
ADB transMult(const ADB &p) const
Definition: BlackoilModelBase_impl.hpp:2778
double maxResidualAllowed() const
Definition: BlackoilModelBase.hpp:519
Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > DataBlock
Definition: BlackoilModelBase_impl.hpp:89
virtual int numPhases() const =0
int numMaterials() const
Definition: BlackoilModelBase_impl.hpp:292
Definition: BlackoilModelEnums.hpp:46
Definition: GridHelpers.hpp:47
Selection. Choose first of two elements if selection basis element is nonnegative.
Definition: AutoDiffHelpers.hpp:359
Definition: BlackoilModelEnums.hpp:37
bool getWellConvergence(const int iteration)
Definition: BlackoilModelBase_impl.hpp:2542
void computeMassFlux(const int actph, const V &transi, const ADB &kr, const ADB &p, const SolutionState &state)
Definition: BlackoilModelBase_impl.hpp:2225
const Vector & poreVolume() const
Definition: GeoProps.hpp:176
Definition: BlackoilModelEnums.hpp:36
V pvdt_
Definition: BlackoilModelBase.hpp:285
bool localWellsActive() const
Definition: BlackoilModelBase.hpp:307
V well_perforation_pressure_diffs_
Definition: BlackoilModelBase.hpp:277
static std::vector< double > computeConnectionPressureDelta(const Wells &wells, const std::vector< double > &z_perf, const std::vector< double > &dens_perf, const double gravity)
void classifyCondition(const ReservoirState &state)
Definition: BlackoilModelBase_impl.hpp:2808
const WellOps wops_
Definition: BlackoilModelBase.hpp:262
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Definition: AutoDiffMatrix.hpp:696
ADB select(const ADB &x1, const ADB &x2) const
Apply selector to ADB quantities.
Definition: AutoDiffHelpers.hpp:409
const BlackoilPropsAdInterface & fluid_
Definition: BlackoilModelBase.hpp:250
ModelTraits< Implementation >::SolutionState SolutionState
Definition: BlackoilModelBase.hpp:109
const VFPProdTable * getTable(const int table_id) const
V threshold_pressures_by_interior_face_
Definition: BlackoilModelBase.hpp:269
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:413
ModelParameters param_
Definition: BlackoilModelBase.hpp:266
void setThresholdPressures(const std::vector< double > &threshold_pressures_by_face)
Set threshold pressures that prevent or reduce flow. This prevents flow across faces if the potential...
Definition: BlackoilModelBase_impl.hpp:317
Eigen::SparseMatrix< double > p2w
Definition: BlackoilModelBase.hpp:244
ADB well_flux_eq
Definition: LinearisedBlackoilResidual.hpp:59
void updatePhaseCondFromPrimalVariable()
Update the phaseCondition_ member based on the primalVariable_ member.
Definition: BlackoilModelBase_impl.hpp:2898
std::vector< int > variableStateIndices() const
Definition: BlackoilModelBase_impl.hpp:525
BlackoilModelBase(const ModelParameters ¶m, const Grid &grid, const BlackoilPropsAdInterface &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const Wells *wells, const NewtonIterationBlackoilInterface &linsolver, Opm::EclipseStateConstPtr eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output)
Definition: BlackoilModelBase_impl.hpp:145
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
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:713
virtual ADB bGas(const ADB &pg, const ADB &T, const ADB &rv, const std::vector< PhasePresence > &cond, const Cells &cells) const =0
Definition: BlackoilModelEnums.hpp:30
const std::vector< int > cells_
Definition: BlackoilModelBase.hpp:260
virtual const boost::any & parallelInformation() const =0
Get the information about the parallelization of the grid.
void solveWellEq(const std::vector< ADB > &mob_perfcells, const std::vector< ADB > &b_perfcells, SolutionState &state, WellState &well_state)
Definition: BlackoilModelBase_impl.hpp:1486
std::vector< bool > activePhases(const PU &pu)
Definition: BlackoilModelBase_impl.hpp:110
SolutionState variableState(const ReservoirState &x, const WellState &xw) const
Definition: BlackoilModelBase_impl.hpp:417
const std::vector< PhasePresence > phaseCondition() const
Definition: BlackoilModelBase.hpp:468
VFPProperties vfp_properties_
Definition: BlackoilModelBase.hpp:254
Definition: AutoDiffMatrix.hpp:43
M ngrad
Extract for each internal face the difference of its adjacent cells' values (first - second)...
Definition: AutoDiffHelpers.hpp:50
void variableReservoirStateInitials(const ReservoirState &x, std::vector< V > &vars0) const
Definition: BlackoilModelBase_impl.hpp:453
ADB::V computeHydrostaticCorrection(const Wells &wells, const ADB::V vfp_ref_depth, const ADB::V &well_perforation_densities, const double gravity)
Definition: BlackoilModelBase_impl.hpp:1311
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.
bool use_threshold_pressure_
Definition: BlackoilModelBase.hpp:267
V nnc_trans
The NNC transmissibilities.
Definition: AutoDiffHelpers.hpp:68
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition: BlackoilModelBase_impl.hpp:256
ADB::V V
Definition: BlackoilModelBase_impl.hpp:84
const Grid & grid_
Definition: BlackoilModelBase.hpp:249
bool empty() const
Definition: VFPProdProperties.hpp:149
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelBase_impl.hpp:280
ADB::M M
Definition: BlackoilModelBase_impl.hpp:85
Definition: BlackoilModelEnums.hpp:44
M div
Extract for each cell the sum of its adjacent interior faces' (signed) values.
Definition: AutoDiffHelpers.hpp:56
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
V well_perforation_densities_
Definition: BlackoilModelBase.hpp:276
V isRv_
Definition: BlackoilModelBase.hpp:274
const bool has_vapoil_
Definition: BlackoilModelBase.hpp:264
ADB fluidDensity(const int phase, const ADB &b, const ADB &rs, const ADB &rv) const
Definition: BlackoilModelBase_impl.hpp:2672
std::vector< PhasePresence > phaseCondition_
Definition: BlackoilModelBase.hpp:272
HelperOps ops_
Definition: BlackoilModelBase.hpp:261
std::vector< V > variableStateInitials(const ReservoirState &x, const WellState &xw) const
Definition: BlackoilModelBase_impl.hpp:431
double dsMax() const
Definition: BlackoilModelBase.hpp:517
Eigen::SparseMatrix< double > w2p
Definition: BlackoilModelBase.hpp:243
const bool has_disgas_
Definition: BlackoilModelBase.hpp:263
ADB poroMult(const ADB &p) const
Definition: BlackoilModelBase_impl.hpp:2748
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Definition: AutoDiffBlock.hpp:183
void assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Definition: BlackoilModelBase_impl.hpp:815
IFaces internal_faces
Definition: AutoDiffHelpers.hpp:47
virtual std::vector< ADB > relperm(const ADB &sw, const ADB &so, const ADB &sg, const Cells &cells) const =0
void computeAccum(const SolutionState &state, const int aix)
Definition: BlackoilModelBase_impl.hpp:655
V solveJacobianSystem() const
Definition: BlackoilModelBase_impl.hpp:1746
Definition: BlackoilModelEnums.hpp:45
V isSg_
Definition: BlackoilModelBase.hpp:275
const double * gravity() const
Definition: GeoProps.hpp:180
void updateEquationsScaling()
Update the scaling factors for mass balance equations.
Definition: BlackoilModelBase_impl.hpp:963
AutoDiffBlock< double > vertcatCollapseJacs(const std::vector< AutoDiffBlock< double > > &x)
Definition: AutoDiffHelpers.hpp:524
Definition: AutoDiffHelpers.hpp:176
WellOps(const Wells *wells)
Definition: BlackoilModelBase_impl.hpp:352
std::vector< ReservoirResidualQuant > rq_
Definition: BlackoilModelBase.hpp:271
double getGravity(const double *g, const int dim)
Definition: BlackoilModelBase_impl.hpp:698
bool terminal_output_
Whether we print something to std::cout.
Definition: BlackoilModelBase.hpp:282
std::vector< int > active2Canonical(const PU &pu)
Definition: BlackoilModelBase_impl.hpp:126
ADB bhp(const std::vector< int > &table_id, const Wells &wells, const ADB &qs, const ADB &thp, const ADB &alq) const
int sizeNonLinear() const
The size (number of unknowns) of the nonlinear system of equations.
Definition: BlackoilModelBase_impl.hpp:244
virtual ADB muOil(const ADB &po, const ADB &T, const ADB &rs, const std::vector< PhasePresence > &cond, const Cells &cells) const =0
V fluidRsSat(const V &p, const V &so, const std::vector< int > &cells) const
Definition: BlackoilModelBase_impl.hpp:2696
void prepareStep(const double dt, ReservoirState &reservoir_state, WellState &well_state)
Definition: BlackoilModelBase_impl.hpp:213
void updatePrimalVariableFromState(const ReservoirState &state)
Definition: BlackoilModelBase_impl.hpp:2848
double dpMaxRel() const
Definition: BlackoilModelBase.hpp:516
ADB fluidReciprocFVF(const int phase, const ADB &p, const ADB &temp, const ADB &rs, const ADB &rv, const std::vector< PhasePresence > &cond) const
Definition: BlackoilModelBase_impl.hpp:2647
void addWellFluxEq(const std::vector< ADB > &cq_s, const SolutionState &state)
Definition: BlackoilModelBase_impl.hpp:1177
const VFPProdProperties * getProd() const
Definition: VFPProperties.hpp:68
std::vector< int > buildAllCells(const int nc)
Definition: BlackoilModelBase_impl.hpp:97
void updateWellControls(WellState &xw) const
Definition: BlackoilModelBase_impl.hpp:1360
std::vector< ADB > computeRelPerm(const SolutionState &state) const
Definition: BlackoilModelBase_impl.hpp:2140
void computeWellFlux(const SolutionState &state, const std::vector< ADB > &mob_perfcells, const std::vector< ADB > &b_perfcells, V &aliveWells, std::vector< ADB > &cq_s)
Definition: BlackoilModelBase_impl.hpp:1003
void assembleMassBalanceEq(const SolutionState &state)
Definition: BlackoilModelBase_impl.hpp:899
ADB well_eq
Definition: LinearisedBlackoilResidual.hpp:64
std::vector< int > variableWellStateIndices() const
Definition: BlackoilModelBase_impl.hpp:548
virtual const double * surfaceDensity(int regionIdx=0) const =0
V computeGasPressure(const V &po, const V &sw, const V &so, const V &sg) const
Definition: BlackoilModelBase_impl.hpp:2206
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModelBase_impl.hpp:268
void addWellContributionToMassBalanceEq(const std::vector< ADB > &cq_s, const SolutionState &state, const WellState &xw)
Definition: BlackoilModelBase_impl.hpp:983
bool empty() const
Definition: VFPInjProperties.hpp:138
ADB fluidViscosity(const int phase, const ADB &p, const ADB &temp, const ADB &rs, const ADB &rv, const std::vector< PhasePresence > &cond) const
Definition: BlackoilModelBase_impl.hpp:2622
ReservoirResidualQuant()
Definition: BlackoilModelBase_impl.hpp:337
const std::string & materialName(int material_index) const
Definition: BlackoilModelBase_impl.hpp:304
ADB::V V
Definition: BlackoilModelBase.hpp:103
void updatePerfPhaseRatesAndPressures(const std::vector< ADB > &cq_s, const SolutionState &state, WellState &xw)
Definition: BlackoilModelBase_impl.hpp:1152
virtual ADB bWat(const ADB &pw, const ADB &T, const Cells &cells) const =0
AutoDiff< Scalar > sqrt(const AutoDiff< Scalar > &x)
Definition: AutoDiff.hpp:287
bool constraintBroken(const std::vector< double > &bhp, const std::vector< double > &thp, const std::vector< double > &well_phase_flow_rate, const int well, const int num_phases, const WellType &well_type, const WellControls *wc, const int ctrl_index)
Definition: BlackoilModelBase_impl.hpp:1213
bool isVFPActive() const
Definition: BlackoilModelBase_impl.hpp:1328
double infinityNormWell(const ADB &a, const boost::any &pinfo)
Compute the L-infinity norm of a vector representing a well equation.
Definition: BlackoilModelBase_impl.hpp:1791
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:436
void applyThresholdPressures(ADB &dp)
Definition: BlackoilModelBase_impl.hpp:2260
const Vector & z() const
Definition: GeoProps.hpp:179
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)
virtual ADB rvSat(const ADB &po, const Cells &cells) const =0
ADB bhp(const std::vector< int > &table_id, const Wells &wells, const ADB &qs, const ADB &thp) const
const std::vector< M > & derivative() const
Function derivatives.
Definition: AutoDiffBlock.hpp:442
void addWellControlEq(const SolutionState &state, const WellState &xw, const V &aliveWells)
Definition: BlackoilModelBase_impl.hpp:1580
ModelTraits< Implementation >::WellState WellState
Definition: BlackoilModelBase.hpp:107
virtual ADB bOil(const ADB &po, const ADB &T, const ADB &rs, const std::vector< PhasePresence > &cond, const Cells &cells) const =0
double thp(int table_id, const double &aqua, const double &liquid, const double &vapour, const double &bhp) const
void variableWellStateInitials(const WellState &xw, std::vector< V > &vars0) const
Definition: BlackoilModelBase_impl.hpp:491
Definition: BlackoilModelEnums.hpp:31
const NewtonIterationBlackoilInterface & linsolver_
Definition: BlackoilModelBase.hpp:255
std::vector< ADB > computePressures(const ADB &po, const ADB &sw, const ADB &so, const ADB &sg) const
Definition: BlackoilModelBase_impl.hpp:2170
virtual PhaseUsage phaseUsage() const =0