All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
BlackoilTransportModel.hpp
1 /*
2  Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
3  Copyright 2016 Statoil AS.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED
22 #define OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED
23 
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>
30 
31 namespace Opm {
32 
34  template<class Grid, class WellModel>
35  class BlackoilTransportModel : public BlackoilModelBase<Grid, WellModel, BlackoilTransportModel<Grid, WellModel> >
36  {
37  public:
39  friend Base;
40 
41  typedef typename Base::ReservoirState ReservoirState;
42  typedef typename Base::WellState WellState;
43  typedef typename Base::SolutionState SolutionState;
44  typedef typename Base::V V;
45 
46 
61  BlackoilTransportModel(const typename Base::ModelParameters& param,
62  const Grid& grid,
63  const BlackoilPropsAdFromDeck& fluid,
64  const DerivedGeology& geo,
65  const RockCompressibility* rock_comp_props,
66  const StandardWells& std_wells,
67  const NewtonIterationBlackoilInterface& linsolver,
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)
74  {
75  }
76 
77 
78 
79 
80 
81  void prepareStep(const SimulatorTimerInterface& timer,
82  const ReservoirState& reservoir_state,
83  const WellState& well_state)
84  {
85  Base::prepareStep(timer, reservoir_state, 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);
90  }
91 
92 
93 
94 
95 
96  SimulatorReport
97  assemble(const ReservoirState& reservoir_state,
98  WellState& well_state,
99  const bool initial_assembly)
100  {
101  using namespace Opm::AutoDiffGrid;
102 
103  SimulatorReport report;
104 
105  // If we have VFP tables, we need the well connection
106  // pressures for the "simple" hydrostatic correction
107  // between well depth and vfp table depth.
108  if (isVFPActive()) {
109  SolutionState state = asImpl().variableState(reservoir_state, well_state);
110  SolutionState state0 = state;
111  asImpl().makeConstantState(state0);
112  asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
113  }
114 
115  // Possibly switch well controls and updating well state to
116  // get reasonable initial conditions for the wells
117  asImpl().wellModel().updateWellControls(well_state);
118 
119  if (initial_assembly) {
120  // HACK
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);
129  }
130 
131  // Create the primary variables.
132  SolutionState state = asImpl().variableState(reservoir_state, well_state);
133 
134  if (initial_assembly) {
135  SolutionState state0 = state;
136  asImpl().makeConstantState(state0);
137  asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
138  }
139 
140 
141  // -------- Mass balance equations --------
142  asImpl().assembleMassBalanceEq(state);
143 
144  // -------- Well equations ----------
145  if ( ! wellsActive() ) {
146  return report;
147  }
148 
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) {
153  // solve the well equations as a pre-processing step
154  report += asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
155  }
156  V aliveWells;
157  std::vector<ADB> cq_s;
158 
159  // @afr changed
160  // asImpl().wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
161  asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
162  // end of changed
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_);
167 
168  return report;
169  }
170 
171 
172 
173 
177  {
178  const int n_transport = residual_.material_balance_eq[1].size();
179  const int n_full = residual_.sizeNonLinear();
180  const auto& mb = residual_.material_balance_eq;
181  LinearisedBlackoilResidual transport_res = {
182  {
183  // TODO: handle general 2-phase etc.
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] })
186  },
187  ADB::null(),
188  ADB::null(),
189  residual_.matbalscale,
190  residual_.singlePrecision
191  };
192  assert(transport_res.sizeNonLinear() == 2*n_transport);
193  V dx_transport = linsolver_.computeNewtonIncrement(transport_res);
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);
198  }
199  return dx_full;
200  }
201 
202 
203 
204 
205 
206  using Base::numPhases;
207  using Base::numMaterials;
208 
209  protected:
210  using Base::asImpl;
211  using Base::materialName;
213  using Base::maxResidualAllowed;
214 
215  using Base::linsolver_;
216  using Base::residual_;
217  using Base::sd_;
218  using Base::geo_;
219  using Base::ops_;
220  using Base::grid_;
221  using Base::use_threshold_pressure_;
222  using Base::canph_;
223  using Base::active_;
224  using Base::pvdt_;
225  using Base::fluid_;
226  using Base::param_;
228 
229  using Base::isVFPActive;
230  using Base::phaseCondition;
231  using Base::vfp_properties_;
232  using Base::wellsActive;
233 
234  V total_flux_; // HACK, should be part of a revised (transport-specific) SolutionState.
235  V total_wellperf_flux_;
236  DataBlock comp_wellperf_flux_;
237  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind_flags_;
238 
239 
240 
241 
242 
243  SolutionState
244  variableState(const ReservoirState& x,
245  const WellState& xw) const
246  {
247  // As Base::variableState(), except making Pressure, Qs and Bhp constants.
248  std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
249  std::vector<ADB> vars = ADB::variables(vars0);
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);
255  }
256 
257 
258 
259 
260 
261 
262 
263  void assembleMassBalanceEq(const SolutionState& state)
264  {
265  // Compute b_p and the accumulation term b_p*s_p for each phase,
266  // except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
267  // These quantities are stored in sd_.rq[phase].accum[1].
268  // The corresponding accumulation terms from the start of
269  // the timestep (b^0_p*s^0_p etc.) were already computed
270  // on the initial call to assemble() and stored in sd_.rq[phase].accum[0].
271  asImpl().computeAccum(state, 1);
272 
273  // Set up the common parts of the mass balance equations
274  // for each active phase.
275  const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
276  const V trans_nnc = ops_.nnc_trans;
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());
281 
282  // Compute mobilities and heads
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) {
287  // Compute and store mobilities.
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];
292  // Note that the pressure-dependent transmissibility multipliers are considered
293  // part of the mobility here.
294  sd_.rq[ phase_idx ].mob = tr_mult * sd_.rq[phase_idx].kr / sd_.rq[phase_idx].mu;
295 
296  // Compute head differentials. Gravity potential is done using the face average as in eclipse and MRST.
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;
300 
301  if (use_threshold_pressure_) {
302  asImpl().applyThresholdPressures(sd_.rq[ phase_idx ].dh);
303  }
304  }
305 
306  // Extract saturation-dependent part of head differences.
307  const ADB gradp = ops_.grad * state.pressure;
308  std::vector<ADB> dh_sat(numPhases(), ADB::null());
309  for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
310  dh_sat[phase_idx] = gradp - sd_.rq[phase_idx].dh;
311  }
312 
313  // Find upstream directions for each phase.
314  upwind_flags_ = multiPhaseUpwind(dh_sat, trans_all);
315 
316  // Compute (upstream) phase and total mobilities for connections.
317  // Also get upstream b, rs, and rv values to avoid recreating the UpwindSelector.
318  std::vector<ADB> mob(numPhases(), ADB::null());
319  std::vector<ADB> b(numPhases(), ADB::null());
320  ADB rs = ADB::null();
321  ADB rv = ADB::null();
322  ADB tot_mob = ADB::constant(V::Zero(gdz.size()));
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);
330  }
331  if (canph_[phase_idx] == Gas) {
332  rv = upwind.select(state.rv);
333  }
334  }
335 
336  // Compute phase fluxes.
337  for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
338  ADB gflux = ADB::constant(V::Zero(gdz.size()));
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]);
342  }
343  }
344  sd_.rq[phase_idx].mflux = b[phase_idx] * (mob[phase_idx] / tot_mob) * (total_flux_ + trans_all * gflux);
345  }
346 
347 #pragma omp parallel for schedule(static)
348  for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
349  // const int canonical_phase_idx = canph_[ phase_idx ];
350  // const ADB& phase_pressure = state.canonical_phase_pressures[canonical_phase_idx];
351  // asImpl().computeMassFlux(phase_idx, trans_all, kr[canonical_phase_idx], phase_pressure, state);
352 
353  // Material balance equation for this phase.
354  residual_.material_balance_eq[ phase_idx ] =
355  pvdt_ * (sd_.rq[phase_idx].accum[1] - sd_.rq[phase_idx].accum[0])
356  + ops_.div*sd_.rq[phase_idx].mflux;
357  }
358 
359  // -------- Extra (optional) rs and rv contributions to the mass balance equations --------
360 
361  // Add the extra (flux) terms to the mass balance equations
362  // From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
363  // The extra terms in the accumulation part of the equation are already handled.
364  if (active_[ Oil ] && active_[ Gas ]) {
365  const int po = fluid_.phaseUsage().phase_pos[ Oil ];
366  const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
367  residual_.material_balance_eq[ pg ] += ops_.div * (rs * sd_.rq[po].mflux);
368  residual_.material_balance_eq[ po ] += ops_.div * (rv * sd_.rq[pg].mflux);
369  }
370 
371  if (param_.update_equations_scaling_) {
372  asImpl().updateEquationsScaling();
373  }
374 
375  }
376 
377 
378 
379 
380 
381  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> multiPhaseUpwind(const std::vector<ADB>& head_diff,
382  const V& transmissibility)
383  {
384  assert(numPhases() == 3);
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];
390  const int a = ops_.connection_cells(conn, 0); // first cell of connection
391  const int b = ops_.connection_cells(conn, 1); // second cell of connection
392  auto up = connectionMultiPhaseUpwind(
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]}},
396  t,
397  q);
398  for (int ii = 0; ii < numPhases(); ++ii) {
399  upwind(conn, ii) = up[ii];
400  }
401  }
402  return upwind;
403  }
404 
405 
406 
407 
408 
409  void computeWellFlux(const SolutionState& state,
410  const std::vector<ADB>& mob_perfcells,
411  const std::vector<ADB>& b_perfcells,
412  V& /* aliveWells */,
413  std::vector<ADB>& cq_s) const
414  {
415  // Note that use of this function replaces using the well models'
416  // function of the same name.
417  if( ! asImpl().localWellsActive() ) return ;
418 
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();
423 
424  // Compute total mobilities for perforations.
425  ADB totmob_perfcells = ADB::constant(V::Zero(nperf));
426  for (int phase = 0; phase < numPhases(); ++phase) {
427  totmob_perfcells += mob_perfcells[phase];
428  }
429 
430  // Compute fractional flow.
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;
434  }
435 
436  // Identify injecting and producing perforations.
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) {
441  is_inj[c] = 1;
442  } else {
443  is_prod[c] = 1;
444  }
445  }
446 
447  // Compute fluxes for producing perforations.
448  std::vector<ADB> cq_s_prod(3, ADB::null());
449  for (int phase = 0; phase < np; ++phase) {
450  // For producers, we use the total reservoir flux from the pressure solver.
451  cq_s_prod[phase] = b_perfcells[phase] * frac_flow[phase] * total_wellperf_flux_;
452  }
453  if (asImpl().has_disgas_ || asImpl().has_vapoil_) {
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;
460  }
461 
462  // Compute well perforation surface volume fluxes.
463  cq_s.resize(np, ADB::null());
464  for (int phase = 0; phase < np; ++phase) {
465  const int pos = pu.phase_pos[phase];
466  // For injectors, we use the component fluxes computed by the pressure solver.
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;
469  }
470  }
471 
472 
473 
474 
475 
476  bool getConvergence(const SimulatorTimerInterface& timer, const int iteration)
477  {
478  const double dt = timer.currentStepLength();
479  const double tol_mb = param_.tolerance_mb_;
480  const double tol_cnv = param_.tolerance_cnv_;
481 
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);
486 
487  const V& pv = geo_.poreVolume();
488 
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);
496 
497  for ( int idx = 0; idx < nm; ++idx )
498  {
499  const ADB& tempB = sd_.rq[idx].b;
500  B.col(idx) = 1./tempB.value();
501  R.col(idx) = residual_.material_balance_eq[idx].value();
502  tempV.col(idx) = R.col(idx).abs()/pv;
503  }
504 
505  const double pvSum = convergenceReduction(B, tempV, R,
506  R_sum, maxCoeff, B_avg, maxNormWell,
507  nc);
508 
509  std::vector<double> CNV(nm);
510  std::vector<double> mass_balance_residual(nm);
511  std::vector<double> well_flux_residual(np);
512 
513  bool converged_MB = true;
514  bool converged_CNV = true;
515  // Finish computation
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);
521  assert(nm >= np);
522  }
523 
524  const bool converged = converged_MB && converged_CNV;
525 
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));
531  }
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));
536  }
537  }
538 
539  if ( terminal_output_ ) {
540  // Only rank 0 does print to std::cout
541  std::ostringstream os;
542  if (iteration == 0) {
543  os << "\nIter";
544  for (int idx = 1; idx < nm; ++idx) {
545  os << " MB(" << materialName(idx).substr(0, 3) << ") ";
546  }
547  for (int idx = 1; idx < nm; ++idx) {
548  os << " CNV(" << materialName(idx).substr(0, 1) << ") ";
549  }
550  os << '\n';
551  }
552  os.precision(3);
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];
557  }
558  for (int idx = 1; idx < nm; ++idx) {
559  os << std::setw(11) << CNV[idx];
560  }
561  OpmLog::info(os.str());
562  }
563  return converged;
564  }
565  };
566 
567 
569  template <class Grid, class WellModel>
570  struct ModelTraits< BlackoilTransportModel<Grid, WellModel> >
571  {
572  typedef BlackoilState ReservoirState;
576  };
577 
578 } // namespace Opm
579 
580 
581 
582 
583 #endif // OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
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.
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
V nnc_trans
The NNC transmissibilities.
Definition: AutoDiffHelpers.hpp:72
M div
Extract for each cell the sum of its adjacent interior faces&#39; (signed) values.
Definition: AutoDiffHelpers.hpp:60
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...
BlackoilTransportModel(const typename Base::ModelParameters &param, 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
bool localWellsActive() const
return true if wells are available on this process
Definition: BlackoilModelBase.hpp:389
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.
V solveJacobianSystem() const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilTransportModel.hpp:176
virtual double currentStepLength() const =0
Current step length.
int sizeNonLinear() const
The size of the non-linear system.
Definition: LinearisedBlackoilResidual.cpp:22
TwoColInt connection_cells
The set of all connections&#39; cells (face or nnc).
Definition: AutoDiffHelpers.hpp:75
M caver
Extract for each face the average of its adjacent cells&#39; 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
int numPhases() const
The number of active fluid phases in the model.
std::vector< ADB > select(const std::vector< ADB > &xc) const
Apply selector to multiple per-cell quantities.
Definition: AutoDiffHelpers.hpp:230
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
int numMaterials() const
The number of active materials in the model.
Struct for containing iteration variables.
Definition: DefaultBlackoilSolutionState.hpp:29
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
BlackoilTransportModel< Grid, WellModel > & asImpl()
Access the most-derived class used for static polymorphism (CRTP).
Definition: BlackoilModelBase.hpp:370
const std::string & materialName(int material_index) const
The name of an active material in the model.
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&#39; values (second - first).
Definition: AutoDiffHelpers.hpp:56
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
PhaseUsage phaseUsage() const
Definition: BlackoilPropsAdFromDeck.cpp:231
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
bool wellsActive() const
return true if wells are available in the reservoir
Definition: BlackoilModelBase.hpp:386
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
A model implementation for the transport equation in three-phase black oil.
Definition: BlackoilTransportModel.hpp:35