21 #ifndef OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED 22 #define OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED 24 #include <opm/autodiff/BlackoilModelBase.hpp> 25 #include <opm/core/simulator/BlackoilState.hpp> 26 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp> 27 #include <opm/autodiff/BlackoilModelParameters.hpp> 28 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp> 29 #include <opm/autodiff/multiPhaseUpwind.hpp> 34 template<
class Gr
id,
class WellModel>
41 typedef typename Base::ReservoirState ReservoirState;
42 typedef typename Base::WellState WellState;
43 typedef typename Base::SolutionState SolutionState;
44 typedef typename Base::V V;
65 const RockCompressibility* rock_comp_props,
68 std::shared_ptr<const EclipseState> eclState,
69 const bool has_disgas,
70 const bool has_vapoil,
71 const bool terminal_output)
72 :
Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
73 eclState, has_disgas, has_vapoil, terminal_output)
82 const ReservoirState& reservoir_state,
83 const WellState& well_state)
86 Base::param_.solve_welleq_initially_ =
false;
87 SolutionState state0 = variableState(reservoir_state, well_state);
88 asImpl().makeConstantState(state0);
89 asImpl().computeAccum(state0, 0);
97 assemble(
const ReservoirState& reservoir_state,
98 WellState& well_state,
99 const bool initial_assembly)
103 SimulatorReport report;
109 SolutionState state =
asImpl().variableState(reservoir_state, well_state);
110 SolutionState state0 = state;
111 asImpl().makeConstantState(state0);
112 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
117 asImpl().wellModel().updateWellControls(well_state);
119 if (initial_assembly) {
121 const_cast<V&
>(total_flux_)
122 = Eigen::Map<const V>(reservoir_state.faceflux().data(), reservoir_state.faceflux().size());
123 const_cast<V&
>(total_wellperf_flux_)
124 = Eigen::Map<const V>(well_state.perfRates().data(), well_state.perfRates().size());
125 const_cast<DataBlock&
>(comp_wellperf_flux_)
126 = Eigen::Map<const DataBlock>(well_state.perfPhaseRates().data(), well_state.perfRates().size(),
numPhases());
127 assert(
numPhases() * well_state.perfRates().size() == well_state.perfPhaseRates().size());
128 asImpl().updatePhaseCondFromPrimalVariable(reservoir_state);
132 SolutionState state =
asImpl().variableState(reservoir_state, well_state);
134 if (initial_assembly) {
135 SolutionState state0 = state;
136 asImpl().makeConstantState(state0);
137 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
142 asImpl().assembleMassBalanceEq(state);
149 std::vector<ADB> mob_perfcells;
150 std::vector<ADB> b_perfcells;
151 asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
152 if (param_.solve_welleq_initially_ && initial_assembly) {
154 report +=
asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
157 std::vector<ADB> cq_s;
161 asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
163 asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
164 asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
165 asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
166 asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
184 ADB::function(mb[1].value(), { mb[1].derivative()[1], mb[1].derivative()[2] }),
185 ADB::function(mb[2].value(), { mb[2].derivative()[1], mb[2].derivative()[2] })
189 residual_.matbalscale,
190 residual_.singlePrecision
194 assert(dx_transport.size() == 2*n_transport);
195 V dx_full = V::Zero(n_full);
196 for (
int i = 0; i < 2*n_transport; ++i) {
197 dx_full(n_transport + i) = dx_transport(i);
213 using Base::maxResidualAllowed;
215 using Base::linsolver_;
216 using Base::residual_;
221 using Base::use_threshold_pressure_;
229 using Base::isVFPActive;
230 using Base::phaseCondition;
231 using Base::vfp_properties_;
235 V total_wellperf_flux_;
236 DataBlock comp_wellperf_flux_;
237 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind_flags_;
244 variableState(
const ReservoirState& x,
245 const WellState& xw)
const 248 std::vector<V> vars0 =
asImpl().variableStateInitials(x, xw);
250 const std::vector<int> indices =
asImpl().variableStateIndices();
251 vars[indices[Pressure]] =
ADB::constant(vars[indices[Pressure]].value());
252 vars[indices[Qs]] =
ADB::constant(vars[indices[Qs]].value());
253 vars[indices[Bhp]] =
ADB::constant(vars[indices[Bhp]].value());
254 return asImpl().variableStateExtractVars(x, indices, vars);
263 void assembleMassBalanceEq(
const SolutionState& state)
271 asImpl().computeAccum(state, 1);
275 const V transi =
subset(geo_.transmissibility(), ops_.internal_faces);
277 V trans_all(transi.size() + trans_nnc.size());
278 trans_all << transi, trans_nnc;
279 const ADB tr_mult =
asImpl().transMult(state.pressure);
280 const V gdz = geo_.gravity()[2] * (ops_.
grad * geo_.z().matrix());
283 const std::vector<PhasePresence>& cond =
asImpl().phaseCondition();
284 const std::vector<ADB> kr =
asImpl().computeRelPerm(state);
285 #pragma omp parallel for schedule(static) 286 for (
int phase_idx = 0; phase_idx <
numPhases(); ++phase_idx) {
288 const int canonical_phase_idx = canph_[ phase_idx ];
289 const ADB& phase_pressure = state.canonical_phase_pressures[canonical_phase_idx];
290 sd_.rq[phase_idx].mu =
asImpl().fluidViscosity(canonical_phase_idx, phase_pressure, state.temperature, state.rs, state.rv, cond);
291 sd_.rq[phase_idx].kr = kr[canonical_phase_idx];
294 sd_.rq[ phase_idx ].mob = tr_mult * sd_.rq[phase_idx].kr / sd_.rq[phase_idx].mu;
297 sd_.rq[phase_idx].rho =
asImpl().fluidDensity(canonical_phase_idx, sd_.rq[phase_idx].b, state.rs, state.rv);
298 const ADB rhoavg = ops_.
caver * sd_.rq[phase_idx].rho;
299 sd_.rq[ phase_idx ].dh = ops_.
grad * phase_pressure - rhoavg * gdz;
301 if (use_threshold_pressure_) {
302 asImpl().applyThresholdPressures(sd_.rq[ phase_idx ].dh);
307 const ADB gradp = ops_.
grad * state.pressure;
309 for (
int phase_idx = 0; phase_idx <
numPhases(); ++phase_idx) {
310 dh_sat[phase_idx] = gradp - sd_.rq[phase_idx].dh;
314 upwind_flags_ = multiPhaseUpwind(dh_sat, trans_all);
323 for (
int phase_idx = 0; phase_idx <
numPhases(); ++phase_idx) {
324 UpwindSelector<double> upwind(grid_, ops_, upwind_flags_.col(phase_idx));
325 mob[phase_idx] = upwind.select(sd_.rq[phase_idx].mob);
326 tot_mob += mob[phase_idx];
327 b[phase_idx] = upwind.select(sd_.rq[phase_idx].b);
328 if (canph_[phase_idx] == Oil) {
329 rs = upwind.select(state.rs);
331 if (canph_[phase_idx] == Gas) {
332 rv = upwind.select(state.rv);
337 for (
int phase_idx = 0; phase_idx <
numPhases(); ++phase_idx) {
339 for (
int other_phase = 0; other_phase <
numPhases(); ++other_phase) {
340 if (phase_idx != other_phase) {
341 gflux += mob[other_phase] * (dh_sat[phase_idx] - dh_sat[other_phase]);
344 sd_.rq[phase_idx].mflux = b[phase_idx] * (mob[phase_idx] / tot_mob) * (total_flux_ + trans_all * gflux);
347 #pragma omp parallel for schedule(static) 348 for (
int phase_idx = 0; phase_idx <
numPhases(); ++phase_idx) {
355 pvdt_ * (sd_.rq[phase_idx].accum[1] - sd_.rq[phase_idx].accum[0])
356 + ops_.
div*sd_.rq[phase_idx].mflux;
364 if (active_[ Oil ] && active_[ Gas ]) {
365 const int po = fluid_.
phaseUsage().phase_pos[ Oil ];
366 const int pg = fluid_.
phaseUsage().phase_pos[ Gas ];
371 if (param_.update_equations_scaling_) {
372 asImpl().updateEquationsScaling();
381 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> multiPhaseUpwind(
const std::vector<ADB>& head_diff,
382 const V& transmissibility)
385 const int num_connections = head_diff[0].size();
386 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind(num_connections,
numPhases());
387 for (
int conn = 0; conn < num_connections; ++conn) {
388 const double q = total_flux_[conn];
389 const double t = transmissibility[conn];
393 {{ head_diff[0].value()[conn], head_diff[1].value()[conn], head_diff[2].value()[conn] }},
394 {{ sd_.rq[0].mob.value()[a], sd_.rq[1].mob.value()[a], sd_.rq[2].mob.value()[a]}},
395 {{ sd_.rq[0].mob.value()[b], sd_.rq[1].mob.value()[b], sd_.rq[2].mob.value()[b]}},
398 for (
int ii = 0; ii <
numPhases(); ++ii) {
399 upwind(conn, ii) = up[ii];
409 void computeWellFlux(
const SolutionState& state,
410 const std::vector<ADB>& mob_perfcells,
411 const std::vector<ADB>& b_perfcells,
413 std::vector<ADB>& cq_s)
const 419 const int np =
asImpl().wells().number_of_phases;
420 const int nw =
asImpl().wells().number_of_wells;
421 const int nperf =
asImpl().wells().well_connpos[nw];
422 const Opm::PhaseUsage& pu =
asImpl().fluid_.phaseUsage();
426 for (
int phase = 0; phase <
numPhases(); ++phase) {
427 totmob_perfcells += mob_perfcells[phase];
431 std::vector<ADB> frac_flow(np,
ADB::null());
432 for (
int phase = 0; phase < np; ++phase) {
433 frac_flow[phase] = mob_perfcells[phase] / totmob_perfcells;
437 V is_inj = V::Zero(nperf);
438 V is_prod = V::Zero(nperf);
439 for (
int c = 0; c < nperf; ++c){
440 if (total_wellperf_flux_[c] > 0.0) {
448 std::vector<ADB> cq_s_prod(3,
ADB::null());
449 for (
int phase = 0; phase < np; ++phase) {
451 cq_s_prod[phase] = b_perfcells[phase] * frac_flow[phase] * total_wellperf_flux_;
454 const int oilpos = pu.phase_pos[Oil];
455 const int gaspos = pu.phase_pos[Gas];
456 const ADB cq_s_prod_oil = cq_s_prod[oilpos];
457 const ADB cq_s_prod_gas = cq_s_prod[gaspos];
458 cq_s_prod[gaspos] +=
subset(state.rs, Base::well_model_.wellOps().well_cells) * cq_s_prod_oil;
459 cq_s_prod[oilpos] +=
subset(state.rv, Base::well_model_.wellOps().well_cells) * cq_s_prod_gas;
464 for (
int phase = 0; phase < np; ++phase) {
465 const int pos = pu.phase_pos[phase];
467 const V cq_s_inj = comp_wellperf_flux_.col(pos);
468 cq_s[phase] = is_prod * cq_s_prod[phase] + is_inj * cq_s_inj;
476 bool getConvergence(
const SimulatorTimerInterface& timer,
const int iteration)
478 const double dt = timer.currentStepLength();
479 const double tol_mb = param_.tolerance_mb_;
480 const double tol_cnv = param_.tolerance_cnv_;
482 const int nc = Opm::AutoDiffGrid::numCells(grid_);
483 const int np =
asImpl().numPhases();
484 const int nm =
asImpl().numMaterials();
485 assert(
int(sd_.rq.size()) == nm);
487 const V& pv = geo_.poreVolume();
489 std::vector<double> R_sum(nm);
490 std::vector<double> B_avg(nm);
491 std::vector<double> maxCoeff(nm);
492 std::vector<double> maxNormWell(np);
493 Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
494 Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
495 Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
497 for (
int idx = 0; idx < nm; ++idx )
499 const ADB& tempB = sd_.rq[idx].b;
500 B.col(idx) = 1./tempB.value();
502 tempV.col(idx) = R.col(idx).abs()/pv;
506 R_sum, maxCoeff, B_avg, maxNormWell,
509 std::vector<double> CNV(nm);
510 std::vector<double> mass_balance_residual(nm);
511 std::vector<double> well_flux_residual(np);
513 bool converged_MB =
true;
514 bool converged_CNV =
true;
516 for (
int idx = 1; idx < nm; ++idx ) {
517 CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
518 mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
519 converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
520 converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
524 const bool converged = converged_MB && converged_CNV;
526 for (
int idx = 0; idx < nm; ++idx) {
527 if (std::isnan(mass_balance_residual[idx])
528 || std::isnan(CNV[idx])
529 || (idx < np && std::isnan(well_flux_residual[idx]))) {
530 OPM_THROW(Opm::NumericalProblem,
"NaN residual for phase " <<
materialName(idx));
532 if (mass_balance_residual[idx] > maxResidualAllowed()
533 || CNV[idx] > maxResidualAllowed()
534 || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
535 OPM_THROW(Opm::NumericalProblem,
"Too large residual for phase " <<
materialName(idx));
541 std::ostringstream os;
542 if (iteration == 0) {
544 for (
int idx = 1; idx < nm; ++idx) {
547 for (
int idx = 1; idx < nm; ++idx) {
548 os <<
" CNV(" <<
materialName(idx).substr(0, 1) <<
") ";
553 os.setf(std::ios::scientific);
554 os << std::setw(4) << iteration;
555 for (
int idx = 1; idx < nm; ++idx) {
556 os << std::setw(11) << mass_balance_residual[idx];
558 for (
int idx = 1; idx < nm; ++idx) {
559 os << std::setw(11) << CNV[idx];
561 OpmLog::info(os.str());
569 template <
class Gr
id,
class WellModel>
572 typedef BlackoilState ReservoirState;
583 #endif // OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
V nnc_trans
The NNC transmissibilities.
Definition: AutoDiffHelpers.hpp:72
V solveJacobianSystem() const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilTransportModel.hpp:176
M div
Extract for each cell the sum of its adjacent interior faces' (signed) values.
Definition: AutoDiffHelpers.hpp:60
bool localWellsActive() const
return true if wells are available on this process
Definition: BlackoilModelBase.hpp:389
virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual &residual) const =0
Solve the linear system Ax = b, with A being the combined derivative matrix of the residual and b bei...
int sizeNonLinear() const
The size of the non-linear system.
Definition: LinearisedBlackoilResidual.cpp:22
BlackoilTransportModel(const typename Base::ModelParameters ¶m, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const StandardWells &std_wells, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output)
Construct the model.
Definition: BlackoilTransportModel.hpp:61
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
std::array< double, 3 > connectionMultiPhaseUpwind(const std::array< double, 3 > &head_diff, const std::array< double, 3 > &mob1, const std::array< double, 3 > &mob2, const double transmissibility, const double flux)
Compute upwind directions for three-phase flow across a connection.
Definition: multiPhaseUpwind.cpp:31
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
Class for handling the standard well model.
Definition: StandardWells.hpp:51
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilModelBase_impl.hpp:220
const std::string & materialName(int material_index) const
The name of an active material in the model.
Definition: BlackoilModelBase_impl.hpp:406
TwoColInt connection_cells
The set of all connections' cells (face or nnc).
Definition: AutoDiffHelpers.hpp:75
M caver
Extract for each face the average of its adjacent cells' values.
Definition: AutoDiffHelpers.hpp:58
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
PhaseUsage phaseUsage() const
Definition: BlackoilPropsAdFromDeck.cpp:231
Traits to encapsulate the types used by classes using or extending this model.
Definition: BlackoilModelBase.hpp:60
Residual structure of the fully implicit solver.
Definition: LinearisedBlackoilResidual.hpp:47
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
bool terminal_output_
Whether we print something to std::cout.
Definition: BlackoilModelBase.hpp:353
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
Struct for containing iteration variables.
Definition: DefaultBlackoilSolutionState.hpp:29
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
BlackoilTransportModel< Grid, WellModel > & asImpl()
Access the most-derived class used for static polymorphism (CRTP).
Definition: BlackoilModelBase.hpp:370
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelBase_impl.hpp:383
bool wellsActive() const
return true if wells are available in the reservoir
Definition: BlackoilModelBase.hpp:386
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
M grad
Extract for each face the difference of its adjacent cells' values (second - first).
Definition: AutoDiffHelpers.hpp:56
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
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) const
Compute the reduction within the convergence check.
Definition: BlackoilModelBase_impl.hpp:1634
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
static std::vector< AutoDiffBlock > variables(const std::vector< V > &initial_values)
Construct a set of primary variables, each initialized to a given vector.
Definition: AutoDiffBlock.hpp:203
int numMaterials() const
The number of active materials in the model.
Definition: BlackoilModelBase_impl.hpp:394
A model implementation for the transport equation in three-phase black oil.
Definition: BlackoilTransportModel.hpp:35