All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
BlackoilModelBase_impl.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2014, 2015 Statoil ASA.
5  Copyright 2015 NTNU
6  Copyright 2015 IRIS AS
7 
8  This file is part of the Open Porous Media project (OPM).
9 
10  OPM is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  OPM is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with OPM. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 #ifndef OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
25 #define OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
26 
27 #include <opm/autodiff/BlackoilModelBase.hpp>
28 #include <opm/autodiff/BlackoilDetails.hpp>
29 #include <opm/autodiff/BlackoilLegacyDetails.hpp>
30 
31 #include <opm/autodiff/AutoDiffBlock.hpp>
32 #include <opm/autodiff/AutoDiffHelpers.hpp>
33 #include <opm/autodiff/GridHelpers.hpp>
34 #include <opm/autodiff/WellHelpers.hpp>
35 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
36 #include <opm/autodiff/GeoProps.hpp>
37 #include <opm/autodiff/WellDensitySegmented.hpp>
38 #include <opm/autodiff/VFPProperties.hpp>
39 #include <opm/autodiff/VFPProdProperties.hpp>
40 #include <opm/autodiff/VFPInjProperties.hpp>
41 
42 #include <opm/core/grid.h>
43 #include <opm/core/linalg/LinearSolverInterface.hpp>
44 #include <opm/core/linalg/ParallelIstlInformation.hpp>
45 #include <opm/core/props/rock/RockCompressibility.hpp>
46 #include <opm/common/ErrorMacros.hpp>
47 #include <opm/common/Exceptions.hpp>
48 #include <opm/common/OpmLog/OpmLog.hpp>
49 #include <opm/parser/eclipse/Units/Units.hpp>
50 #include <opm/core/well_controls.h>
51 #include <opm/core/utility/parameters/ParameterGroup.hpp>
52 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
53 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
54 
55 #include <dune/common/timer.hh>
56 
57 #include <cassert>
58 #include <cmath>
59 #include <iostream>
60 #include <iomanip>
61 #include <limits>
62 #include <vector>
63 #include <algorithm>
64 //#include <fstream>
65 
66 // A debugging utility.
67 #define OPM_AD_DUMP(foo) \
68  do { \
69  std::cout << "==========================================\n" \
70  << #foo ":\n" \
71  << collapseJacs(foo) << std::endl; \
72  } while (0)
73 
74 #define OPM_AD_DUMPVAL(foo) \
75  do { \
76  std::cout << "==========================================\n" \
77  << #foo ":\n" \
78  << foo.value() << std::endl; \
79  } while (0)
80 
81 #define OPM_AD_DISKVAL(foo) \
82  do { \
83  std::ofstream os(#foo); \
84  os.precision(16); \
85  os << foo.value() << std::endl; \
86  } while (0)
87 
88 
89 namespace Opm {
90 
91 typedef AutoDiffBlock<double> ADB;
92 typedef ADB::V V;
93 typedef ADB::M M;
94 typedef Eigen::Array<double,
95  Eigen::Dynamic,
96  Eigen::Dynamic,
97  Eigen::RowMajor> DataBlock;
98 
99  template <class Grid, class WellModel, class Implementation>
101  BlackoilModelBase(const ModelParameters& param,
102  const Grid& grid ,
103  const BlackoilPropsAdFromDeck& fluid,
104  const DerivedGeology& geo ,
105  const RockCompressibility* rock_comp_props,
106  const WellModel& well_model,
107  const NewtonIterationBlackoilInterface& linsolver,
108  std::shared_ptr< const Opm::EclipseState > eclState,
109  const bool has_disgas,
110  const bool has_vapoil,
111  const bool terminal_output)
112  : grid_ (grid)
113  , fluid_ (fluid)
114  , geo_ (geo)
115  , rock_comp_props_(rock_comp_props)
116  , vfp_properties_(
117  eclState->getTableManager().getVFPInjTables(),
118  eclState->getTableManager().getVFPProdTables())
119  , linsolver_ (linsolver)
120  , active_(detail::activePhases(fluid.phaseUsage()))
121  , canph_ (detail::active2Canonical(fluid.phaseUsage()))
122  , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
123  , ops_ (grid, geo.nnc())
124  , has_disgas_(has_disgas)
125  , has_vapoil_(has_vapoil)
126  , param_( param )
127  , use_threshold_pressure_(false)
128  , sd_ (fluid.numPhases())
129  , phaseCondition_(AutoDiffGrid::numCells(grid))
130  , well_model_ (well_model)
131  , isRs_(V::Zero(AutoDiffGrid::numCells(grid)))
132  , isRv_(V::Zero(AutoDiffGrid::numCells(grid)))
133  , isSg_(V::Zero(AutoDiffGrid::numCells(grid)))
134  , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
135  ADB::null(),
136  ADB::null(),
137  { 1.1169, 1.0031, 0.0031 }, // the default magic numbers
138  false } )
139  , terminal_output_ (terminal_output)
140  , material_name_(0)
141  , current_relaxation_(1.0)
142  // only one region 0 used, which means average reservoir hydrocarbon conditions in
143  // the field will be calculated.
144  // TODO: more delicate implementation will be required if we want to handle different
145  // FIP regions specified from the well specifications.
146  , rate_converter_(fluid_.phaseUsage(), std::vector<int>(AutoDiffGrid::numCells(grid_),0))
147  {
148  if (active_[Water]) {
149  material_name_.push_back("Water");
150  }
151  if (active_[Oil]) {
152  material_name_.push_back("Oil");
153  }
154  if (active_[Gas]) {
155  material_name_.push_back("Gas");
156  }
157 
158  assert(numMaterials() == std::accumulate(active_.begin(), active_.end(), 0)); // Due to the material_name_ init above.
159 
160  const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
161  const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_);
162 
163  well_model_.init(&fluid_, &active_, &phaseCondition_, &vfp_properties_, gravity, depth);
164 
165  // TODO: put this for now to avoid modify the following code.
166  // TODO: this code can be fragile.
167 #if HAVE_MPI
168  const Wells* wells_arg = asImpl().well_model_.wellsPointer();
169 
170  if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
171  {
172  const ParallelISTLInformation& info =
173  boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
174  if ( terminal_output_ ) {
175  // Only rank 0 does print to std::cout if terminal_output is enabled
176  terminal_output_ = (info.communicator().rank()==0);
177  }
178  int local_number_of_wells = localWellsActive() ? wells().number_of_wells : 0;
179  int global_number_of_wells = info.communicator().sum(local_number_of_wells);
180  const bool wells_active = ( wells_arg && global_number_of_wells > 0 );
181  wellModel().setWellsActive(wells_active);
182  // Compute the global number of cells
183  std::vector<int> v( Opm::AutoDiffGrid::numCells(grid_), 1);
184  global_nc_ = 0;
185  info.computeReduction(v, Opm::Reduction::makeGlobalSumFunctor<int>(), global_nc_);
186  }else
187 #endif
188  {
189  wellModel().setWellsActive( localWellsActive() );
190  global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
191  }
192  }
193 
194 
195  template <class Grid, class WellModel, class Implementation>
196  bool
198  isParallel() const
199  {
200 #if HAVE_MPI
201  if ( linsolver_.parallelInformation().type() !=
202  typeid(ParallelISTLInformation) )
203  {
204  return false;
205  }
206  else
207  {
208  const auto& comm =boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation()).communicator();
209  return comm.size() > 1;
210  }
211 #else
212  return false;
213 #endif
214  }
215 
216 
217  template <class Grid, class WellModel, class Implementation>
218  void
221  const ReservoirState& reservoir_state,
222  const WellState& /* well_state */)
223  {
224  const double dt = timer.currentStepLength();
225 
226  pvdt_ = geo_.poreVolume() / dt;
227  if (active_[Gas]) {
228  updatePrimalVariableFromState(reservoir_state);
229  }
230  }
231 
232 
233 
234 
235 
236  template <class Grid, class WellModel, class Implementation>
237  template <class NonlinearSolverType>
238  SimulatorReport
240  nonlinearIteration(const int iteration,
241  const SimulatorTimerInterface& timer,
242  NonlinearSolverType& nonlinear_solver,
243  ReservoirState& reservoir_state,
244  WellState& well_state)
245  {
246  SimulatorReport report;
247  Dune::Timer perfTimer;
248 
249  perfTimer.start();
250  const double dt = timer.currentStepLength();
251 
252  if (iteration == 0) {
253  // For each iteration we store in a vector the norms of the residual of
254  // the mass balance for each active phase, the well flux and the well equations.
255  residual_norms_history_.clear();
256  current_relaxation_ = 1.0;
257  dx_old_ = V::Zero(sizeNonLinear());
258  }
259  try {
260  report += asImpl().assemble(reservoir_state, well_state, iteration == 0);
261  report.assemble_time += perfTimer.stop();
262  }
263  catch (...) {
264  report.assemble_time += perfTimer.stop();
265  throw;
266  }
267 
268  report.total_linearizations = 1;
269  perfTimer.reset();
270  perfTimer.start();
271  report.converged = asImpl().getConvergence(timer, iteration);
272  residual_norms_history_.push_back(asImpl().computeResidualNorms());
273  report.update_time += perfTimer.stop();
274 
275  const bool must_solve = (iteration < nonlinear_solver.minIter()) || (!report.converged);
276  if (must_solve) {
277  perfTimer.reset();
278  perfTimer.start();
279  report.total_newton_iterations = 1;
280 
281  // enable single precision for solvers when dt is smaller then maximal time step for single precision
282  residual_.singlePrecision = ( dt < param_.maxSinglePrecisionTimeStep_ );
283 
284  // Compute the nonlinear update.
285  V dx;
286  try {
287  dx = asImpl().solveJacobianSystem();
288  report.linear_solve_time += perfTimer.stop();
289  report.total_linear_iterations += linearIterationsLastSolve();
290  }
291  catch (...) {
292  report.linear_solve_time += perfTimer.stop();
293  report.total_linear_iterations += linearIterationsLastSolve();
294  throw;
295  }
296 
297  perfTimer.reset();
298  perfTimer.start();
299 
300  if (param_.use_update_stabilization_) {
301  // Stabilize the nonlinear update.
302  bool isOscillate = false;
303  bool isStagnate = false;
304  nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
305  if (isOscillate) {
306  current_relaxation_ -= nonlinear_solver.relaxIncrement();
307  current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
308  if (terminalOutputEnabled()) {
309  std::string msg = " Oscillating behavior detected: Relaxation set to "
310  + std::to_string(current_relaxation_);
311  OpmLog::info(msg);
312  }
313  }
314  nonlinear_solver.stabilizeNonlinearUpdate(dx, dx_old_, current_relaxation_);
315  }
316 
317  // Apply the update, applying model-dependent
318  // limitations and chopping of the update.
319  asImpl().updateState(dx, reservoir_state, well_state);
320  report.update_time += perfTimer.stop();
321  }
322 
323  return report;
324  }
325 
326 
327 
328 
329 
330  template <class Grid, class WellModel, class Implementation>
331  void
334  ReservoirState& /* reservoir_state */,
335  WellState& /* well_state */)
336  {
337  // Does nothing in this model.
338  }
339 
340 
341 
342 
343 
344  template <class Grid, class WellModel, class Implementation>
345  int
348  {
349  return residual_.sizeNonLinear();
350  }
351 
352 
353 
354 
355 
356  template <class Grid, class WellModel, class Implementation>
357  int
360  {
361  return linsolver_.iterations();
362  }
363 
364 
365 
366 
367 
368  template <class Grid, class WellModel, class Implementation>
369  bool
372  {
373  return terminal_output_;
374  }
375 
376 
377 
378 
379 
380  template <class Grid, class WellModel, class Implementation>
381  int
383  numPhases() const
384  {
385  return fluid_.numPhases();
386  }
387 
388 
389 
390 
391  template <class Grid, class WellModel, class Implementation>
392  int
395  {
396  return material_name_.size();
397  }
398 
399 
400 
401 
402 
403  template <class Grid, class WellModel, class Implementation>
404  const std::string&
406  materialName(int material_index) const
407  {
408  assert(material_index < numMaterials());
409  return material_name_[material_index];
410  }
411 
412 
413 
414 
415 
416  template <class Grid, class WellModel, class Implementation>
417  void
419  setThresholdPressures(const std::vector<double>& threshold_pressures)
420  {
421  const int num_faces = AutoDiffGrid::numFaces(grid_);
422  const int num_nnc = geo_.nnc().numNNC();
423  const int num_connections = num_faces + num_nnc;
424  if (int(threshold_pressures.size()) != num_connections) {
425  OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures input ( " << threshold_pressures.size()
426  << " ), must be equal to number of faces + nncs ( " << num_faces << " + " << num_nnc << " ).");
427  }
428  use_threshold_pressure_ = true;
429  // Map to interior faces.
430  const int num_ifaces = ops_.internal_faces.size();
431  threshold_pressures_by_connection_.resize(num_ifaces + num_nnc);
432  for (int ii = 0; ii < num_ifaces; ++ii) {
433  threshold_pressures_by_connection_[ii] = threshold_pressures[ops_.internal_faces[ii]];
434  }
435  // Handle NNCs
436  // Note: the nnc threshold pressures is appended after the face threshold pressures
437  for (int ii = 0; ii < num_nnc; ++ii) {
438  threshold_pressures_by_connection_[ii + num_ifaces] = threshold_pressures[ii + num_faces];
439  }
440 
441  }
442 
443 
444 
445 
446  template <class Grid, class WellModel, class Implementation>
449  : accum(2, ADB::null())
450  , mflux( ADB::null())
451  , b ( ADB::null())
452  , mu ( ADB::null())
453  , rho ( ADB::null())
454  , kr ( ADB::null())
455  , dh ( ADB::null())
456  , mob ( ADB::null())
457  {
458  }
459 
460  template <class Grid, class WellModel, class Implementation>
461  BlackoilModelBase<Grid, WellModel, Implementation>::
462  SimulatorData::SimulatorData(int num_phases)
463  : rq(num_phases)
464  , rsSat(ADB::null())
465  , rvSat(ADB::null())
466  , soMax()
467  , Pb()
468  , Pd()
469  , krnswdc_ow()
470  , krnswdc_go()
471  , pcswmdc_ow()
472  , pcswmdc_go()
473  , fip()
474  {
475  }
476 
477 
478 
479 
480 
481  template <class Grid, class WellModel, class Implementation>
482  void
483  BlackoilModelBase<Grid, WellModel, Implementation>::
484  makeConstantState(SolutionState& state) const
485  {
486  // HACK: throw away the derivatives. this may not be the most
487  // performant way to do things, but it will make the state
488  // automatically consistent with variableState() (and doing
489  // things automatically is all the rage in this module ;)
490  state.pressure = ADB::constant(state.pressure.value());
491  state.temperature = ADB::constant(state.temperature.value());
492  state.rs = ADB::constant(state.rs.value());
493  state.rv = ADB::constant(state.rv.value());
494  const int num_phases = state.saturation.size();
495  for (int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) {
496  state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value());
497  }
498  state.qs = ADB::constant(state.qs.value());
499  state.bhp = ADB::constant(state.bhp.value());
500  assert(state.canonical_phase_pressures.size() == static_cast<std::size_t>(Opm::BlackoilPhases::MaxNumPhases));
501  for (int canphase = 0; canphase < Opm::BlackoilPhases::MaxNumPhases; ++canphase) {
502  ADB& pp = state.canonical_phase_pressures[canphase];
503  pp = ADB::constant(pp.value());
504  }
505  }
506 
507 
508 
509 
510 
511  template <class Grid, class WellModel, class Implementation>
512  typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
513  BlackoilModelBase<Grid, WellModel, Implementation>::
514  variableState(const ReservoirState& x,
515  const WellState& xw) const
516  {
517  std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
518  std::vector<ADB> vars = ADB::variables(vars0);
519  return asImpl().variableStateExtractVars(x, asImpl().variableStateIndices(), vars);
520  }
521 
522 
523 
524 
525 
526  template <class Grid, class WellModel, class Implementation>
527  std::vector<V>
528  BlackoilModelBase<Grid, WellModel, Implementation>::
529  variableStateInitials(const ReservoirState& x,
530  const WellState& xw) const
531  {
532  assert(active_[ Oil ]);
533 
534  const int np = x.numPhases();
535 
536  std::vector<V> vars0;
537  // p, Sw and Rs, Rv or Sg is used as primary depending on solution conditions
538  // and bhp and Q for the wells
539  vars0.reserve(np + 1);
540  variableReservoirStateInitials(x, vars0);
541  asImpl().wellModel().variableWellStateInitials(xw, vars0);
542  return vars0;
543  }
544 
545 
546 
547 
548 
549  template <class Grid, class WellModel, class Implementation>
550  void
551  BlackoilModelBase<Grid, WellModel, Implementation>::
552  variableReservoirStateInitials(const ReservoirState& x, std::vector<V>& vars0) const
553  {
554  using namespace Opm::AutoDiffGrid;
555  const int nc = numCells(grid_);
556  const int np = x.numPhases();
557  // Initial pressure.
558  assert (not x.pressure().empty());
559  const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
560  vars0.push_back(p);
561 
562  // Initial saturation.
563  assert (not x.saturation().empty());
564  const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
565  const Opm::PhaseUsage pu = fluid_.phaseUsage();
566  // We do not handle a Water/Gas situation correctly, guard against it.
567  assert (active_[ Oil]);
568  if (active_[ Water ]) {
569  const V sw = s.col(pu.phase_pos[ Water ]);
570  vars0.push_back(sw);
571  }
572 
573  if (active_[ Gas ]) {
574  // define new primary variable xvar depending on solution condition
575  V xvar(nc);
576  const V sg = s.col(pu.phase_pos[ Gas ]);
577  const V rs = Eigen::Map<const V>(& x.gasoilratio()[0], x.gasoilratio().size());
578  const V rv = Eigen::Map<const V>(& x.rv()[0], x.rv().size());
579  xvar = isRs_*rs + isRv_*rv + isSg_*sg;
580  vars0.push_back(xvar);
581  }
582  }
583 
584 
585 
586 
587 
588  template <class Grid, class WellModel, class Implementation>
589  std::vector<int>
590  BlackoilModelBase<Grid, WellModel, Implementation>::
591  variableStateIndices() const
592  {
593  assert(active_[Oil]);
594  std::vector<int> indices(5, -1);
595  int next = 0;
596  indices[Pressure] = next++;
597  if (active_[Water]) {
598  indices[Sw] = next++;
599  }
600  if (active_[Gas]) {
601  indices[Xvar] = next++;
602  }
603  asImpl().wellModel().variableStateWellIndices(indices, next);
604  assert(next == fluid_.numPhases() + 2);
605  return indices;
606  }
607 
608 
609 
610 
611 
612  template <class Grid, class WellModel, class Implementation>
613  typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
614  BlackoilModelBase<Grid, WellModel, Implementation>::
615  variableStateExtractVars(const ReservoirState& x,
616  const std::vector<int>& indices,
617  std::vector<ADB>& vars) const
618  {
619  //using namespace Opm::AutoDiffGrid;
620  const int nc = Opm::AutoDiffGrid::numCells(grid_);
621  const Opm::PhaseUsage pu = fluid_.phaseUsage();
622 
623  SolutionState state(fluid_.numPhases());
624 
625  // Pressure.
626  state.pressure = std::move(vars[indices[Pressure]]);
627 
628  // Temperature cannot be a variable at this time (only constant).
629  const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
630  state.temperature = ADB::constant(temp);
631 
632  // Saturations
633  {
634  ADB so = ADB::constant(V::Ones(nc, 1));
635 
636  if (active_[ Water ]) {
637  state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]);
638  const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
639  so -= sw;
640  }
641 
642  if (active_[ Gas ]) {
643  // Define Sg Rs and Rv in terms of xvar.
644  // Xvar is only defined if gas phase is active
645  const ADB& xvar = vars[indices[Xvar]];
646  ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
647  sg = isSg_*xvar + isRv_*so;
648  so -= sg;
649 
650  //Compute the phase pressures before computing RS/RV
651  {
652  const ADB& sw = (active_[ Water ]
653  ? state.saturation[ pu.phase_pos[ Water ] ]
654  : ADB::null());
655  state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
656  }
657 
658  if (active_[ Oil ]) {
659  // RS and RV is only defined if both oil and gas phase are active.
660  sd_.rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
661  if (has_disgas_) {
662  state.rs = (1-isRs_)*sd_.rsSat + isRs_*xvar;
663  } else {
664  state.rs = sd_.rsSat;
665  }
666  sd_.rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
667  if (has_vapoil_) {
668  state.rv = (1-isRv_)*sd_.rvSat + isRv_*xvar;
669  } else {
670  state.rv = sd_.rvSat;
671  }
672  sd_.soMax = fluid_.satOilMax();
673  fluid_.getGasOilHystParams(sd_.krnswdc_go, sd_.pcswmdc_go, cells_);
674  fluid_.getOilWaterHystParams(sd_.krnswdc_ow, sd_.pcswmdc_ow, cells_);
675 
676  sd_.Pb = fluid_.bubblePointPressure(cells_,
677  state.temperature.value(),
678  state.rs.value());
679  sd_.Pd = fluid_.dewPointPressure(cells_,
680  state.temperature.value(),
681  state.rv.value());
682  }
683  }
684  else {
685  // Compute phase pressures also if gas phase is not active
686  const ADB& sw = (active_[ Water ]
687  ? state.saturation[ pu.phase_pos[ Water ] ]
688  : ADB::null());
689  const ADB& sg = ADB::null();
690  state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
691  }
692 
693  if (active_[ Oil ]) {
694  // Note that so is never a primary variable.
695  state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
696  }
697  }
698  // wells
699  asImpl().wellModel().variableStateExtractWellsVars(indices, vars, state);
700  return state;
701  }
702 
703 
704 
705 
706 
707  template <class Grid, class WellModel, class Implementation>
708  void
709  BlackoilModelBase<Grid, WellModel, Implementation>::
710  computeAccum(const SolutionState& state,
711  const int aix )
712  {
713  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
714 
715  const ADB& press = state.pressure;
716  const ADB& temp = state.temperature;
717  const std::vector<ADB>& sat = state.saturation;
718  const ADB& rs = state.rs;
719  const ADB& rv = state.rv;
720 
721  const std::vector<PhasePresence> cond = phaseCondition();
722 
723  const ADB pv_mult = poroMult(press);
724 
725  const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
726  for (int phase = 0; phase < maxnp; ++phase) {
727  if (active_[ phase ]) {
728  const int pos = pu.phase_pos[ phase ];
729  sd_.rq[pos].b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
730  sd_.rq[pos].accum[aix] = pv_mult * sd_.rq[pos].b * sat[pos];
731  // OPM_AD_DUMP(sd_.rq[pos].b);
732  // OPM_AD_DUMP(sd_.rq[pos].accum[aix]);
733  }
734  }
735 
736  if (active_[ Oil ] && active_[ Gas ]) {
737  // Account for gas dissolved in oil and vaporized oil
738  const int po = pu.phase_pos[ Oil ];
739  const int pg = pu.phase_pos[ Gas ];
740 
741  // Temporary copy to avoid contribution of dissolved gas in the vaporized oil
742  // when both dissolved gas and vaporized oil are present.
743  const ADB accum_gas_copy =sd_.rq[pg].accum[aix];
744 
745  sd_.rq[pg].accum[aix] += state.rs * sd_.rq[po].accum[aix];
746  sd_.rq[po].accum[aix] += state.rv * accum_gas_copy;
747  // OPM_AD_DUMP(sd_.rq[pg].accum[aix]);
748  }
749  }
750 
751 
752 
753 
754 
755  template <class Grid, class WellModel, class Implementation>
756  SimulatorReport
758  assemble(const ReservoirState& reservoir_state,
759  WellState& well_state,
760  const bool initial_assembly)
761  {
762  using namespace Opm::AutoDiffGrid;
763 
764  SimulatorReport report;
765 
766  // If we have VFP tables, we need the well connection
767  // pressures for the "simple" hydrostatic correction
768  // between well depth and vfp table depth.
769  if (isVFPActive()) {
770  SolutionState state = asImpl().variableState(reservoir_state, well_state);
771  SolutionState state0 = state;
772  asImpl().makeConstantState(state0);
773  asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
774  }
775 
776  // set up the guide rate and group control
777  if (asImpl().wellModel().wellCollection()->groupControlActive() && initial_assembly) {
778  setupGroupControl(reservoir_state, well_state);
779  }
780 
781  // Possibly switch well controls and updating well state to
782  // get reasonable initial conditions for the wells
783  asImpl().wellModel().updateWellControls(well_state);
784 
785  if (asImpl().wellModel().wellCollection()->groupControlActive()) {
786  // enforce VREP control when necessary.
787  applyVREPGroupControl(reservoir_state, well_state);
788 
789  asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
790  }
791 
792  // Create the primary variables.
793  SolutionState state = asImpl().variableState(reservoir_state, well_state);
794 
795  if (initial_assembly) {
796  // Create the (constant, derivativeless) initial state.
797  SolutionState state0 = state;
798  asImpl().makeConstantState(state0);
799  // Compute initial accumulation contributions
800  // and well connection pressures.
801  asImpl().computeAccum(state0, 0);
802  asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
803  }
804 
805  // OPM_AD_DISKVAL(state.pressure);
806  // OPM_AD_DISKVAL(state.saturation[0]);
807  // OPM_AD_DISKVAL(state.saturation[1]);
808  // OPM_AD_DISKVAL(state.saturation[2]);
809  // OPM_AD_DISKVAL(state.rs);
810  // OPM_AD_DISKVAL(state.rv);
811  // OPM_AD_DISKVAL(state.qs);
812  // OPM_AD_DISKVAL(state.bhp);
813 
814  // -------- Mass balance equations --------
815  asImpl().assembleMassBalanceEq(state);
816 
817  // -------- Well equations ----------
818  if ( ! wellsActive() ) {
819  return report;
820  }
821 
822  std::vector<ADB> mob_perfcells;
823  std::vector<ADB> b_perfcells;
824  asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
825  if (param_.solve_welleq_initially_ && initial_assembly) {
826  // solve the well equations as a pre-processing step
827  report += asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
828  }
829  V aliveWells;
830  std::vector<ADB> cq_s;
831  asImpl().wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
832  asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
833  asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
834  asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
835  asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
836 
837  return report;
838  }
839 
840 
841 
842 
843  template <class Grid, class WellModel, class Implementation>
844  void
846  assembleMassBalanceEq(const SolutionState& state)
847  {
848  // Compute b_p and the accumulation term b_p*s_p for each phase,
849  // except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
850  // These quantities are stored in sd_.rq[phase].accum[1].
851  // The corresponding accumulation terms from the start of
852  // the timestep (b^0_p*s^0_p etc.) were already computed
853  // on the initial call to assemble() and stored in sd_.rq[phase].accum[0].
854  asImpl().computeAccum(state, 1);
855 
856  // Set up the common parts of the mass balance equations
857  // for each active phase.
858  const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
859  const V trans_nnc = ops_.nnc_trans;
860  V trans_all(transi.size() + trans_nnc.size());
861  trans_all << transi, trans_nnc;
862 
863 
864  {
865  const std::vector<ADB> kr = asImpl().computeRelPerm(state);
866  for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
867  sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
868  }
869  }
870 #pragma omp parallel for schedule(static)
871  for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
872  const std::vector<PhasePresence>& cond = phaseCondition();
873  sd_.rq[phaseIdx].mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
874  sd_.rq[phaseIdx].rho = asImpl().fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
875  asImpl().computeMassFlux(phaseIdx, trans_all, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
876 
877  residual_.material_balance_eq[ phaseIdx ] =
878  pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
879  + ops_.div*sd_.rq[phaseIdx].mflux;
880  }
881 
882  // -------- Extra (optional) rs and rv contributions to the mass balance equations --------
883 
884  // Add the extra (flux) terms to the mass balance equations
885  // From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
886  // The extra terms in the accumulation part of the equation are already handled.
887  if (active_[ Oil ] && active_[ Gas ]) {
888  const int po = fluid_.phaseUsage().phase_pos[ Oil ];
889  const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
890 
891  const UpwindSelector<double> upwindOil(grid_, ops_,
892  sd_.rq[po].dh.value());
893  const ADB rs_face = upwindOil.select(state.rs);
894 
895  const UpwindSelector<double> upwindGas(grid_, ops_,
896  sd_.rq[pg].dh.value());
897  const ADB rv_face = upwindGas.select(state.rv);
898 
899  residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
900  residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
901 
902  // OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
903 
904  }
905 
906 
907  if (param_.update_equations_scaling_) {
908  asImpl().updateEquationsScaling();
909  }
910 
911  }
912 
913 
914 
915 
916 
917  template <class Grid, class WellModel, class Implementation>
918  void
921  ADB::V B;
922  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
923  for ( int idx=0; idx<MaxNumPhases; ++idx )
924  {
925  if (active_[idx]) {
926  const int pos = pu.phase_pos[idx];
927  const ADB& temp_b = sd_.rq[pos].b;
928  B = 1. / temp_b.value();
929 #if HAVE_MPI
930  if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
931  {
932  const ParallelISTLInformation& real_info =
933  boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
934  double B_global_sum = 0;
935  real_info.computeReduction(B, Reduction::makeGlobalSumFunctor<double>(), B_global_sum);
936  residual_.matbalscale[idx] = B_global_sum / global_nc_;
937  }
938  else
939 #endif
940  {
941  residual_.matbalscale[idx] = B.mean();
942  }
943  }
944  }
945  }
946 
947 
948 
949 
950 
951  template <class Grid, class WellModel, class Implementation>
952  void
954  addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
955  const SolutionState&,
956  const WellState&)
957  {
958  if ( !asImpl().localWellsActive() )
959  {
960  // If there are no wells in the subdomain of the proces then
961  // cq_s has zero size and will cause a segmentation fault below.
962  return;
963  }
964 
965  // Add well contributions to mass balance equations
966  const int nc = Opm::AutoDiffGrid::numCells(grid_);
967  const int np = asImpl().numPhases();
968  const V& efficiency_factors = wellModel().wellPerfEfficiencyFactors();
969  for (int phase = 0; phase < np; ++phase) {
970  residual_.material_balance_eq[phase] -= superset(efficiency_factors * cq_s[phase],
971  wellModel().wellOps().well_cells, nc);
972  }
973  }
974 
975 
976 
977 
978 
979  template <class Grid, class WellModel, class Implementation>
980  bool
981  BlackoilModelBase<Grid, WellModel, Implementation>::
982  isVFPActive() const
983  {
984  if( ! localWellsActive() ) {
985  return false;
986  }
987 
988  if ( vfp_properties_.getProd()->empty() && vfp_properties_.getInj()->empty() ) {
989  return false;
990  }
991 
992  const int nw = wells().number_of_wells;
993  //Loop over all wells
994  for (int w = 0; w < nw; ++w) {
995  const WellControls* wc = wells().ctrls[w];
996 
997  const int nwc = well_controls_get_num(wc);
998 
999  //Loop over all controls
1000  for (int c=0; c < nwc; ++c) {
1001  const WellControlType ctrl_type = well_controls_iget_type(wc, c);
1002 
1003  if (ctrl_type == THP) {
1004  return true;
1005  }
1006  }
1007  }
1008 
1009  return false;
1010  }
1011 
1012 
1013 
1014 
1015  template <class Grid, class WellModel, class Implementation>
1016  SimulatorReport
1017  BlackoilModelBase<Grid, WellModel, Implementation>::
1018  solveWellEq(const std::vector<ADB>& mob_perfcells,
1019  const std::vector<ADB>& b_perfcells,
1020  const ReservoirState& reservoir_state,
1021  SolutionState& state,
1022  WellState& well_state)
1023  {
1024  V aliveWells;
1025  const int np = wells().number_of_phases;
1026  std::vector<ADB> cq_s(np, ADB::null());
1027  std::vector<int> indices = asImpl().wellModel().variableWellStateIndices();
1028  SolutionState state0 = state;
1029  WellState well_state0 = well_state;
1030  asImpl().makeConstantState(state0);
1031 
1032  std::vector<ADB> mob_perfcells_const(np, ADB::null());
1033  std::vector<ADB> b_perfcells_const(np, ADB::null());
1034 
1035  if (asImpl().localWellsActive() ){
1036  // If there are non well in the sudomain of the process
1037  // thene mob_perfcells_const and b_perfcells_const would be empty
1038  for (int phase = 0; phase < np; ++phase) {
1039  mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value());
1040  b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value());
1041  }
1042  }
1043 
1044  int it = 0;
1045  bool converged;
1046  do {
1047  // bhp and Q for the wells
1048  std::vector<V> vars0;
1049  vars0.reserve(2);
1050  asImpl().wellModel().variableWellStateInitials(well_state, vars0);
1051  std::vector<ADB> vars = ADB::variables(vars0);
1052 
1053  SolutionState wellSolutionState = state0;
1054  asImpl().wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState);
1055  asImpl().wellModel().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
1056  asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
1057  asImpl().wellModel().addWellFluxEq(cq_s, wellSolutionState, residual_);
1058  asImpl().wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_);
1059  converged = getWellConvergence(it);
1060 
1061  if (converged) {
1062  break;
1063  }
1064 
1065  ++it;
1066  if( localWellsActive() )
1067  {
1068  std::vector<ADB> eqs;
1069  eqs.reserve(2);
1070  eqs.push_back(residual_.well_flux_eq);
1071  eqs.push_back(residual_.well_eq);
1072  ADB total_residual = vertcatCollapseJacs(eqs);
1073  const std::vector<M>& Jn = total_residual.derivative();
1074  typedef Eigen::SparseMatrix<double> Sp;
1075  Sp Jn0;
1076  Jn[0].toSparse(Jn0);
1077  const Eigen::SparseLU< Sp > solver(Jn0);
1078  ADB::V total_residual_v = total_residual.value();
1079  const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
1080  assert(dx.size() == total_residual_v.size());
1081  asImpl().wellModel().updateWellState(dx.array(), dbhpMaxRel(), well_state);
1082  }
1083  // We have to update the well controls regardless whether there are local
1084  // wells active or not as parallel logging will take place that needs to
1085  // communicate with all processes.
1086  asImpl().wellModel().updateWellControls(well_state);
1087 
1088  if (asImpl().wellModel().wellCollection()->groupControlActive()) {
1089  // Enforce the VREP control
1090  applyVREPGroupControl(reservoir_state, well_state);
1091  asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
1092  }
1093  } while (it < 15);
1094 
1095  if (converged) {
1096  if (terminalOutputEnabled())
1097  {
1098  OpmLog::note("well converged iter: " + std::to_string(it));
1099  }
1100 
1101  const int nw = wells().number_of_wells;
1102  {
1103  // We will set the bhp primary variable to the new ones,
1104  // but we do not change the derivatives here.
1105  ADB::V new_bhp = Eigen::Map<ADB::V>(well_state.bhp().data(), nw);
1106  // Avoiding the copy below would require a value setter method
1107  // in AutoDiffBlock.
1108  std::vector<ADB::M> old_derivs = state.bhp.derivative();
1109  state.bhp = ADB::function(std::move(new_bhp), std::move(old_derivs));
1110  }
1111  {
1112  // Need to reshuffle well rates, from phase running fastest
1113  // to wells running fastest.
1114  // The transpose() below switches the ordering.
1115  const DataBlock wrates = Eigen::Map<const DataBlock>(well_state.wellRates().data(), nw, np).transpose();
1116  ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np);
1117  std::vector<ADB::M> old_derivs = state.qs.derivative();
1118  state.qs = ADB::function(std::move(new_qs), std::move(old_derivs));
1119  }
1120  asImpl().computeWellConnectionPressures(state, well_state);
1121  }
1122 
1123  if (!converged) {
1124  well_state = well_state0;
1125  }
1126 
1127  SimulatorReport report;
1128  report.total_well_iterations = it;
1129  report.converged = converged;
1130  return report;
1131  }
1132 
1133 
1134 
1135 
1136 
1137  template <class Grid, class WellModel, class Implementation>
1138  V
1141  {
1142  return linsolver_.computeNewtonIncrement(residual_);
1143  }
1144 
1145  template <class Grid, class WellModel, class Implementation>
1146  void
1148  updateState(const V& dx,
1149  ReservoirState& reservoir_state,
1150  WellState& well_state)
1151  {
1152  using namespace Opm::AutoDiffGrid;
1153  const int np = fluid_.numPhases();
1154  const int nc = numCells(grid_);
1155  const V null;
1156  assert(null.size() == 0);
1157  const V zero = V::Zero(nc);
1158  const V ones = V::Constant(nc,1.0);
1159 
1160  // Extract parts of dx corresponding to each part.
1161  const V dp = subset(dx, Span(nc));
1162  int varstart = nc;
1163  const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
1164  varstart += dsw.size();
1165 
1166  const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
1167  varstart += dxvar.size();
1168 
1169  // Extract well parts np phase rates + bhp
1170  const V dwells = subset(dx, Span(asImpl().wellModel().numWellVars(), 1, varstart));
1171  varstart += dwells.size();
1172 
1173  assert(varstart == dx.size());
1174 
1175  // Pressure update.
1176  const double dpmaxrel = dpMaxRel();
1177  const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
1178  const V absdpmax = dpmaxrel*p_old.abs();
1179  const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
1180  const V p = (p_old - dp_limited).max(zero);
1181  std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
1182 
1183  // Saturation updates.
1184  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
1185  const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
1186  const double dsmax = dsMax();
1187 
1188  V so;
1189  V sw;
1190  V sg;
1191 
1192  {
1193  V maxVal = zero;
1194  V dso = zero;
1195  if (active_[Water]){
1196  maxVal = dsw.abs().max(maxVal);
1197  dso = dso - dsw;
1198  }
1199 
1200  V dsg;
1201  if (active_[Gas]){
1202  dsg = isSg_ * dxvar - isRv_ * dsw;
1203  maxVal = dsg.abs().max(maxVal);
1204  dso = dso - dsg;
1205  }
1206 
1207  maxVal = dso.abs().max(maxVal);
1208 
1209  V step = dsmax/maxVal;
1210  step = step.min(1.);
1211 
1212  if (active_[Water]) {
1213  const int pos = pu.phase_pos[ Water ];
1214  const V sw_old = s_old.col(pos);
1215  sw = sw_old - step * dsw;
1216  }
1217 
1218  if (active_[Gas]) {
1219  const int pos = pu.phase_pos[ Gas ];
1220  const V sg_old = s_old.col(pos);
1221  sg = sg_old - step * dsg;
1222  }
1223 
1224  assert(active_[Oil]);
1225  const int pos = pu.phase_pos[ Oil ];
1226  const V so_old = s_old.col(pos);
1227  so = so_old - step * dso;
1228  }
1229 
1230  if (active_[Gas]) {
1231  auto ixg = sg < 0;
1232  for (int c = 0; c < nc; ++c) {
1233  if (ixg[c]) {
1234  if (active_[Water]) {
1235  sw[c] = sw[c] / (1-sg[c]);
1236  }
1237  so[c] = so[c] / (1-sg[c]);
1238  sg[c] = 0;
1239  }
1240  }
1241  }
1242 
1243  if (active_[Oil]) {
1244  auto ixo = so < 0;
1245  for (int c = 0; c < nc; ++c) {
1246  if (ixo[c]) {
1247  if (active_[Water]) {
1248  sw[c] = sw[c] / (1-so[c]);
1249  }
1250  if (active_[Gas]) {
1251  sg[c] = sg[c] / (1-so[c]);
1252  }
1253  so[c] = 0;
1254  }
1255  }
1256  }
1257 
1258  if (active_[Water]) {
1259  auto ixw = sw < 0;
1260  for (int c = 0; c < nc; ++c) {
1261  if (ixw[c]) {
1262  so[c] = so[c] / (1-sw[c]);
1263  if (active_[Gas]) {
1264  sg[c] = sg[c] / (1-sw[c]);
1265  }
1266  sw[c] = 0;
1267  }
1268  }
1269  }
1270 
1271  // Update rs and rv
1272  const double drmaxrel = drMaxRel();
1273  V rs;
1274  if (has_disgas_) {
1275  const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
1276  const V drs = isRs_ * dxvar;
1277  const V drs_limited = sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
1278  rs = rs_old - drs_limited;
1279  rs = rs.max(zero);
1280  }
1281  V rv;
1282  if (has_vapoil_) {
1283  const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
1284  const V drv = isRv_ * dxvar;
1285  const V drv_limited = sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
1286  rv = rv_old - drv_limited;
1287  rv = rv.max(zero);
1288  }
1289 
1290  // Sg is used as primal variable for water only cells.
1291  const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
1292  auto watOnly = sw > (1 - epsilon);
1293 
1294  // phase translation sg <-> rs
1295  std::vector<HydroCarbonState>& hydroCarbonState = reservoir_state.hydroCarbonState();
1296  std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
1297 
1298  if (has_disgas_) {
1299  const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
1300  const V rsSat = fluidRsSat(p, so, cells_);
1301  sd_.rsSat = ADB::constant(rsSat);
1302  // The obvious case
1303  auto hasGas = (sg > 0 && isRs_ == 0);
1304 
1305  // Set oil saturated if previous rs is sufficiently large
1306  const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
1307  auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
1308  auto useSg = watOnly || hasGas || gasVaporized;
1309  for (int c = 0; c < nc; ++c) {
1310  if (useSg[c]) {
1311  rs[c] = rsSat[c];
1312  if (watOnly[c]) {
1313  so[c] = 0;
1314  sg[c] = 0;
1315  rs[c] = 0;
1316  }
1317  } else {
1318  hydroCarbonState[c] = HydroCarbonState::OilOnly;
1319  }
1320  }
1321  //rs = rs.min(rsSat);
1322  }
1323 
1324  // phase transitions so <-> rv
1325  if (has_vapoil_) {
1326 
1327  // The gas pressure is needed for the rvSat calculations
1328  const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
1329  const V gaspress = computeGasPressure(p, sw, so, sg);
1330  const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
1331  const V rvSat = fluidRvSat(gaspress, so, cells_);
1332  sd_.rvSat = ADB::constant(rvSat);
1333 
1334  // The obvious case
1335  auto hasOil = (so > 0 && isRv_ == 0);
1336 
1337  // Set oil saturated if previous rv is sufficiently large
1338  const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
1339  auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
1340  auto useSg = watOnly || hasOil || oilCondensed;
1341  for (int c = 0; c < nc; ++c) {
1342  if (useSg[c]) {
1343  rv[c] = rvSat[c];
1344  if (watOnly[c]) {
1345  so[c] = 0;
1346  sg[c] = 0;
1347  rv[c] = 0;
1348  }
1349  } else {
1350  hydroCarbonState[c] = HydroCarbonState::GasOnly;
1351  }
1352  }
1353  //rv = rv.min(rvSat);
1354  }
1355 
1356  // Update the reservoir_state
1357  if (active_[Water]) {
1358  for (int c = 0; c < nc; ++c) {
1359  reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
1360  }
1361  }
1362 
1363  if (active_[Gas]) {
1364  for (int c = 0; c < nc; ++c) {
1365  reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
1366  }
1367  }
1368 
1369  if (active_[ Oil ]) {
1370  for (int c = 0; c < nc; ++c) {
1371  reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
1372  }
1373  }
1374 
1375  if (has_disgas_) {
1376  std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
1377  }
1378 
1379  if (has_vapoil_) {
1380  std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
1381  }
1382 
1383  asImpl().wellModel().updateWellState(dwells, dbhpMaxRel(), well_state);
1384 
1385  // Update phase conditions used for property calculations.
1386  updatePhaseCondFromPrimalVariable(reservoir_state);
1387  }
1388 
1389 
1390 
1391 
1392 
1393  template <class Grid, class WellModel, class Implementation>
1394  std::vector<ADB>
1396  computeRelPerm(const SolutionState& state) const
1397  {
1398  using namespace Opm::AutoDiffGrid;
1399  const int nc = numCells(grid_);
1400 
1401  const ADB zero = ADB::constant(V::Zero(nc));
1402 
1403  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
1404  const ADB& sw = (active_[ Water ]
1405  ? state.saturation[ pu.phase_pos[ Water ] ]
1406  : zero);
1407 
1408  const ADB& so = (active_[ Oil ]
1409  ? state.saturation[ pu.phase_pos[ Oil ] ]
1410  : zero);
1411 
1412  const ADB& sg = (active_[ Gas ]
1413  ? state.saturation[ pu.phase_pos[ Gas ] ]
1414  : zero);
1415 
1416  return fluid_.relperm(sw, so, sg, cells_);
1417  }
1418 
1419 
1420 
1421 
1422 
1423  template <class Grid, class WellModel, class Implementation>
1424  std::vector<ADB>
1425  BlackoilModelBase<Grid, WellModel, Implementation>::
1426  computePressures(const ADB& po,
1427  const ADB& sw,
1428  const ADB& so,
1429  const ADB& sg) const
1430  {
1431  // convert the pressure offsets to the capillary pressures
1432  std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
1433  for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
1434  // The reference pressure is always the liquid phase (oil) pressure.
1435  if (phaseIdx == BlackoilPhases::Liquid)
1436  continue;
1437  if (active_[phaseIdx]) {
1438  pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
1439  }
1440  }
1441 
1442  // Since pcow = po - pw, but pcog = pg - po,
1443  // we have
1444  // pw = po - pcow
1445  // pg = po + pcgo
1446  // This is an unfortunate inconsistency, but a convention we must handle.
1447  for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
1448  if (active_[phaseIdx]) {
1449  if (phaseIdx == BlackoilPhases::Aqua) {
1450  pressure[phaseIdx] = po - pressure[phaseIdx];
1451  } else {
1452  pressure[phaseIdx] += po;
1453  }
1454  }
1455  }
1456 
1457  return pressure;
1458  }
1459 
1460 
1461 
1462 
1463 
1464  template <class Grid, class WellModel, class Implementation>
1465  V
1466  BlackoilModelBase<Grid, WellModel, Implementation>::
1467  computeGasPressure(const V& po,
1468  const V& sw,
1469  const V& so,
1470  const V& sg) const
1471  {
1472  assert (active_[Gas]);
1473  std::vector<ADB> cp = fluid_.capPress(ADB::constant(sw),
1474  ADB::constant(so),
1475  ADB::constant(sg),
1476  cells_);
1477  return cp[Gas].value() + po;
1478  }
1479 
1480 
1481 
1482  template <class Grid, class WellModel, class Implementation>
1483  void
1484  BlackoilModelBase<Grid, WellModel, Implementation>::
1485  computeMassFlux(const int actph ,
1486  const V& transi,
1487  const ADB& kr ,
1488  const ADB& mu ,
1489  const ADB& rho ,
1490  const ADB& phasePressure,
1491  const SolutionState& state)
1492  {
1493  // Compute and store mobilities.
1494  const ADB tr_mult = transMult(state.pressure);
1495  sd_.rq[ actph ].mob = tr_mult * kr / mu;
1496 
1497  // Compute head differentials. Gravity potential is done using the face average as in eclipse and MRST.
1498  const ADB rhoavg = ops_.caver * rho;
1499  sd_.rq[ actph ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
1500  if (use_threshold_pressure_) {
1501  applyThresholdPressures(sd_.rq[ actph ].dh);
1502  }
1503 
1504  // Compute phase fluxes with upwinding of formation value factor and mobility.
1505  const ADB& b = sd_.rq[ actph ].b;
1506  const ADB& mob = sd_.rq[ actph ].mob;
1507  const ADB& dh = sd_.rq[ actph ].dh;
1508  UpwindSelector<double> upwind(grid_, ops_, dh.value());
1509  sd_.rq[ actph ].mflux = upwind.select(b * mob) * (transi * dh);
1510  }
1511 
1512 
1513 
1514 
1515 
1516  template <class Grid, class WellModel, class Implementation>
1517  void
1518  BlackoilModelBase<Grid, WellModel, Implementation>::
1519  applyThresholdPressures(ADB& dp)
1520  {
1521  // We support reversible threshold pressures only.
1522  // Method: if the potential difference is lower (in absolute
1523  // value) than the threshold for any face, then the potential
1524  // (and derivatives) is set to zero. If it is above the
1525  // threshold, the threshold pressure is subtracted from the
1526  // absolute potential (the potential is moved towards zero).
1527 
1528  // Identify the set of faces where the potential is under the
1529  // threshold, that shall have zero flow. Storing the bool
1530  // Array as a V (a double Array) with 1 and 0 elements, a
1531  // 1 where flow is allowed, a 0 where it is not.
1532  const V high_potential = (dp.value().abs() >= threshold_pressures_by_connection_).template cast<double>();
1533 
1534  // Create a sparse vector that nullifies the low potential elements.
1535  const M keep_high_potential(high_potential.matrix().asDiagonal());
1536 
1537  // Find the current sign for the threshold modification
1538  const V sign_dp = sign(dp.value());
1539  const V threshold_modification = sign_dp * threshold_pressures_by_connection_;
1540 
1541  // Modify potential and nullify where appropriate.
1542  dp = keep_high_potential * (dp - threshold_modification);
1543  }
1544 
1545 
1546 
1547 
1548 
1549  template <class Grid, class WellModel, class Implementation>
1550  std::vector<double>
1553  {
1554  std::vector<double> residualNorms;
1555 
1556  std::vector<ADB>::const_iterator massBalanceIt = residual_.material_balance_eq.begin();
1557  const std::vector<ADB>::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
1558 
1559  for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
1560  const double massBalanceResid = detail::infinityNorm( (*massBalanceIt),
1561  linsolver_.parallelInformation() );
1562  if (!std::isfinite(massBalanceResid)) {
1563  OPM_THROW(Opm::NumericalProblem,
1564  "Encountered a non-finite residual");
1565  }
1566  residualNorms.push_back(massBalanceResid);
1567  }
1568 
1569  // the following residuals are not used in the oscillation detection now
1570  const double wellFluxResid = detail::infinityNormWell( residual_.well_flux_eq,
1571  linsolver_.parallelInformation() );
1572  if (!std::isfinite(wellFluxResid)) {
1573  OPM_THROW(Opm::NumericalProblem,
1574  "Encountered a non-finite residual");
1575  }
1576  residualNorms.push_back(wellFluxResid);
1577 
1578  const double wellResid = detail::infinityNormWell( residual_.well_eq,
1579  linsolver_.parallelInformation() );
1580  if (!std::isfinite(wellResid)) {
1581  OPM_THROW(Opm::NumericalProblem,
1582  "Encountered a non-finite residual");
1583  }
1584  residualNorms.push_back(wellResid);
1585 
1586  return residualNorms;
1587  }
1588 
1589 
1590  template <class Grid, class WellModel, class Implementation>
1591  double
1593  relativeChange(const SimulationDataContainer& previous,
1594  const SimulationDataContainer& current ) const
1595  {
1596  std::vector< double > p0 ( previous.pressure() );
1597  std::vector< double > sat0( previous.saturation() );
1598 
1599  const std::size_t pSize = p0.size();
1600  const std::size_t satSize = sat0.size();
1601 
1602  // compute u^n - u^n+1
1603  for( std::size_t i=0; i<pSize; ++i ) {
1604  p0[ i ] -= current.pressure()[ i ];
1605  }
1606 
1607  for( std::size_t i=0; i<satSize; ++i ) {
1608  sat0[ i ] -= current.saturation()[ i ];
1609  }
1610 
1611  // compute || u^n - u^n+1 ||
1612  const double stateOld = detail::euclidianNormSquared( p0.begin(), p0.end(), 1, linsolver_.parallelInformation() ) +
1613  detail::euclidianNormSquared( sat0.begin(), sat0.end(),
1614  current.numPhases(),
1615  linsolver_.parallelInformation() );
1616 
1617  // compute || u^n+1 ||
1618  const double stateNew = detail::euclidianNormSquared( current.pressure().begin(), current.pressure().end(), 1, linsolver_.parallelInformation() ) +
1619  detail::euclidianNormSquared( current.saturation().begin(), current.saturation().end(),
1620  current.numPhases(),
1621  linsolver_.parallelInformation() );
1622 
1623  if( stateNew > 0.0 ) {
1624  return stateOld / stateNew ;
1625  }
1626  else {
1627  return 0.0;
1628  }
1629  }
1630 
1631  template <class Grid, class WellModel, class Implementation>
1632  double
1634  convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
1635  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
1636  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
1637  std::vector<double>& R_sum,
1638  std::vector<double>& maxCoeff,
1639  std::vector<double>& B_avg,
1640  std::vector<double>& maxNormWell,
1641  int nc) const
1642  {
1643  const int np = asImpl().numPhases();
1644  const int nm = asImpl().numMaterials();
1645  const int nw = residual_.well_flux_eq.size() / np;
1646  assert(nw * np == int(residual_.well_flux_eq.size()));
1647 
1648  // Do the global reductions
1649 #if HAVE_MPI
1650  if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
1651  {
1652  const ParallelISTLInformation& info =
1653  boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
1654 
1655  // Compute the global number of cells and porevolume
1656  std::vector<int> v(nc, 1);
1657  auto nc_and_pv = std::tuple<int, double>(0, 0.0);
1658  auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
1659  Opm::Reduction::makeGlobalSumFunctor<double>());
1660  auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume());
1661  info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
1662 
1663  for ( int idx = 0; idx < nm; ++idx )
1664  {
1665  auto values = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
1666  auto containers = std::make_tuple(B.col(idx),
1667  tempV.col(idx),
1668  R.col(idx));
1669  auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
1670  Opm::Reduction::makeGlobalMaxFunctor<double>(),
1671  Opm::Reduction::makeGlobalSumFunctor<double>());
1672  info.computeReduction(containers, operators, values);
1673  B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
1674  maxCoeff[idx] = std::get<1>(values);
1675  R_sum[idx] = std::get<2>(values);
1676  assert(nm >= np);
1677  if (idx < np) {
1678  maxNormWell[idx] = 0.0;
1679  for ( int w = 0; w < nw; ++w ) {
1680  maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
1681  }
1682  }
1683  }
1684  info.communicator().max(maxNormWell.data(), np);
1685  // Compute pore volume
1686  return std::get<1>(nc_and_pv);
1687  }
1688  else
1689 #endif
1690  {
1691  B_avg.resize(nm);
1692  maxCoeff.resize(nm);
1693  R_sum.resize(nm);
1694  maxNormWell.resize(np);
1695  for ( int idx = 0; idx < nm; ++idx )
1696  {
1697  B_avg[idx] = B.col(idx).sum()/nc;
1698  maxCoeff[idx] = tempV.col(idx).maxCoeff();
1699  R_sum[idx] = R.col(idx).sum();
1700 
1701  assert(nm >= np);
1702  if (idx < np) {
1703  maxNormWell[idx] = 0.0;
1704  for ( int w = 0; w < nw; ++w ) {
1705  maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
1706  }
1707  }
1708  }
1709  // Compute total pore volume
1710  return geo_.poreVolume().sum();
1711  }
1712  }
1713 
1714 
1715 
1716 
1717 
1718  template <class Grid, class WellModel, class Implementation>
1719  bool
1721  getConvergence(const SimulatorTimerInterface& timer, const int iteration)
1722  {
1723  const double dt = timer.currentStepLength();
1724  const double tol_mb = param_.tolerance_mb_;
1725  const double tol_cnv = param_.tolerance_cnv_;
1726  const double tol_wells = param_.tolerance_wells_;
1727  const double tol_well_control = param_.tolerance_well_control_;
1728 
1729  const int nc = Opm::AutoDiffGrid::numCells(grid_);
1730  const int np = asImpl().numPhases();
1731  const int nm = asImpl().numMaterials();
1732  assert(int(sd_.rq.size()) == nm);
1733 
1734  const V& pv = geo_.poreVolume();
1735 
1736  std::vector<double> R_sum(nm);
1737  std::vector<double> B_avg(nm);
1738  std::vector<double> maxCoeff(nm);
1739  std::vector<double> maxNormWell(np);
1740  Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
1741  Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
1742  Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
1743 
1744  for ( int idx = 0; idx < nm; ++idx )
1745  {
1746  const ADB& tempB = sd_.rq[idx].b;
1747  B.col(idx) = 1./tempB.value();
1748  R.col(idx) = residual_.material_balance_eq[idx].value();
1749  tempV.col(idx) = R.col(idx).abs()/pv;
1750  }
1751 
1752  const double pvSum = convergenceReduction(B, tempV, R,
1753  R_sum, maxCoeff, B_avg, maxNormWell,
1754  nc);
1755 
1756  std::vector<double> CNV(nm);
1757  std::vector<double> mass_balance_residual(nm);
1758  std::vector<double> well_flux_residual(np);
1759 
1760  bool converged_MB = true;
1761  bool converged_CNV = true;
1762  bool converged_Well = true;
1763  // Finish computation
1764  for ( int idx = 0; idx < nm; ++idx )
1765  {
1766  CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
1767  mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
1768  converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
1769  converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
1770  // Well flux convergence is only for fluid phases, not other materials
1771  // in our current implementation.
1772  assert(nm >= np);
1773  if (idx < np) {
1774  well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
1775  converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
1776  }
1777  }
1778 
1779  const double residualWell = detail::infinityNormWell(residual_.well_eq,
1780  linsolver_.parallelInformation());
1781  converged_Well = converged_Well && (residualWell < tol_well_control);
1782 
1783  const bool converged = converged_MB && converged_CNV && converged_Well;
1784 
1785  // Residual in Pascal can have high values and still be ok.
1786  const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed();
1787 
1788  if ( terminal_output_ )
1789  {
1790  // Only rank 0 does print to std::cout
1791  if (iteration == 0) {
1792  std::string msg = "Iter";
1793  for (int idx = 0; idx < nm; ++idx) {
1794  msg += " MB(" + materialName(idx).substr(0, 3) + ") ";
1795  }
1796  for (int idx = 0; idx < nm; ++idx) {
1797  msg += " CNV(" + materialName(idx).substr(0, 1) + ") ";
1798  }
1799  for (int idx = 0; idx < np; ++idx) {
1800  msg += " W-FLUX(" + materialName(idx).substr(0, 1) + ")";
1801  }
1802  msg += " WELL-CONT";
1803  // std::cout << " WELL-CONT ";
1804  OpmLog::note(msg);
1805  }
1806  std::ostringstream ss;
1807  const std::streamsize oprec = ss.precision(3);
1808  const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1809  ss << std::setw(4) << iteration;
1810  for (int idx = 0; idx < nm; ++idx) {
1811  ss << std::setw(11) << mass_balance_residual[idx];
1812  }
1813  for (int idx = 0; idx < nm; ++idx) {
1814  ss << std::setw(11) << CNV[idx];
1815  }
1816  for (int idx = 0; idx < np; ++idx) {
1817  ss << std::setw(11) << well_flux_residual[idx];
1818  }
1819  ss << std::setw(11) << residualWell;
1820  // std::cout << std::setw(11) << residualWell;
1821  ss.precision(oprec);
1822  ss.flags(oflags);
1823  OpmLog::note(ss.str());
1824  }
1825 
1826  for (int idx = 0; idx < nm; ++idx) {
1827  if (std::isnan(mass_balance_residual[idx])
1828  || std::isnan(CNV[idx])
1829  || (idx < np && std::isnan(well_flux_residual[idx]))) {
1830  const auto msg = std::string("NaN residual for phase ") + materialName(idx);
1831  if (terminal_output_) {
1832  OpmLog::bug(msg);
1833  }
1834  OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1835  }
1836  if (mass_balance_residual[idx] > maxResidualAllowed()
1837  || CNV[idx] > maxResidualAllowed()
1838  || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
1839  const auto msg = std::string("Too large residual for phase ") + materialName(idx);
1840  if (terminal_output_) {
1841  OpmLog::problem(msg);
1842  }
1843  OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1844  }
1845  }
1846  if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) {
1847  const auto msg = std::string("NaN or too large residual for well control equation");
1848  if (terminal_output_) {
1849  OpmLog::problem(msg);
1850  }
1851  OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1852  }
1853 
1854  return converged;
1855  }
1856 
1857 
1858 
1859 
1860 
1861  template <class Grid, class WellModel, class Implementation>
1862  bool
1864  getWellConvergence(const int iteration)
1865  {
1866  const double tol_wells = param_.tolerance_wells_;
1867  const double tol_well_control = param_.tolerance_well_control_;
1868 
1869  const int nc = Opm::AutoDiffGrid::numCells(grid_);
1870  const int np = asImpl().numPhases();
1871  const int nm = asImpl().numMaterials();
1872 
1873  const V& pv = geo_.poreVolume();
1874  std::vector<double> R_sum(nm);
1875  std::vector<double> B_avg(nm);
1876  std::vector<double> maxCoeff(nm);
1877  std::vector<double> maxNormWell(np);
1878  Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
1879  Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
1880  Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
1881  for ( int idx = 0; idx < nm; ++idx )
1882  {
1883  const ADB& tempB = sd_.rq[idx].b;
1884  B.col(idx) = 1./tempB.value();
1885  R.col(idx) = residual_.material_balance_eq[idx].value();
1886  tempV.col(idx) = R.col(idx).abs()/pv;
1887  }
1888 
1889  convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc);
1890 
1891  std::vector<double> well_flux_residual(np);
1892  bool converged_Well = true;
1893  // Finish computation
1894  for ( int idx = 0; idx < np; ++idx )
1895  {
1896  well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
1897  converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
1898  }
1899 
1900  const double residualWell = detail::infinityNormWell(residual_.well_eq,
1901  linsolver_.parallelInformation());
1902  converged_Well = converged_Well && (residualWell < tol_well_control);
1903 
1904  const bool converged = converged_Well;
1905 
1906  // if one of the residuals is NaN, throw exception, so that the solver can be restarted
1907  for (int idx = 0; idx < np; ++idx) {
1908  if (std::isnan(well_flux_residual[idx])) {
1909  const auto msg = std::string("NaN residual for phase ") + materialName(idx);
1910  if (terminal_output_) {
1911  OpmLog::bug(msg);
1912  }
1913  OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1914  }
1915  if (well_flux_residual[idx] > maxResidualAllowed()) {
1916  const auto msg = std::string("Too large residual for phase ") + materialName(idx);
1917  if (terminal_output_) {
1918  OpmLog::problem(msg);
1919  }
1920  OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1921  }
1922  }
1923 
1924  if ( terminal_output_ )
1925  {
1926  // Only rank 0 does print to std::cout
1927  if (iteration == 0) {
1928  std::string msg;
1929  msg = "Iter";
1930  for (int idx = 0; idx < np; ++idx) {
1931  msg += " W-FLUX(" + materialName(idx).substr(0, 1) + ")";
1932  }
1933  msg += " WELL-CONT";
1934  OpmLog::note(msg);
1935  }
1936  std::ostringstream ss;
1937  const std::streamsize oprec = ss.precision(3);
1938  const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1939  ss << std::setw(4) << iteration;
1940  for (int idx = 0; idx < np; ++idx) {
1941  ss << std::setw(11) << well_flux_residual[idx];
1942  }
1943  ss << std::setw(11) << residualWell;
1944  ss.precision(oprec);
1945  ss.flags(oflags);
1946  OpmLog::note(ss.str());
1947  }
1948  return converged;
1949  }
1950 
1951 
1952 
1953 
1954 
1955  template <class Grid, class WellModel, class Implementation>
1956  ADB
1957  BlackoilModelBase<Grid, WellModel, Implementation>::
1958  fluidViscosity(const int phase,
1959  const ADB& p ,
1960  const ADB& temp ,
1961  const ADB& rs ,
1962  const ADB& rv ,
1963  const std::vector<PhasePresence>& cond) const
1964  {
1965  switch (phase) {
1966  case Water:
1967  return fluid_.muWat(p, temp, cells_);
1968  case Oil:
1969  return fluid_.muOil(p, temp, rs, cond, cells_);
1970  case Gas:
1971  return fluid_.muGas(p, temp, rv, cond, cells_);
1972  default:
1973  OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
1974  }
1975  }
1976 
1977 
1978 
1979 
1980 
1981  template <class Grid, class WellModel, class Implementation>
1982  ADB
1983  BlackoilModelBase<Grid, WellModel, Implementation>::
1984  fluidReciprocFVF(const int phase,
1985  const ADB& p ,
1986  const ADB& temp ,
1987  const ADB& rs ,
1988  const ADB& rv ,
1989  const std::vector<PhasePresence>& cond) const
1990  {
1991  switch (phase) {
1992  case Water:
1993  return fluid_.bWat(p, temp, cells_);
1994  case Oil:
1995  return fluid_.bOil(p, temp, rs, cond, cells_);
1996  case Gas:
1997  return fluid_.bGas(p, temp, rv, cond, cells_);
1998  default:
1999  OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
2000  }
2001  }
2002 
2003 
2004 
2005 
2006 
2007  template <class Grid, class WellModel, class Implementation>
2008  ADB
2009  BlackoilModelBase<Grid, WellModel, Implementation>::
2010  fluidDensity(const int phase,
2011  const ADB& b,
2012  const ADB& rs,
2013  const ADB& rv) const
2014  {
2015  const V& rhos = fluid_.surfaceDensity(phase, cells_);
2016  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
2017  ADB rho = rhos * b;
2018  if (phase == Oil && active_[Gas]) {
2019  rho += fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * rs * b;
2020  }
2021  if (phase == Gas && active_[Oil]) {
2022  rho += fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * rv * b;
2023  }
2024  return rho;
2025  }
2026 
2027 
2028 
2029 
2030 
2031  template <class Grid, class WellModel, class Implementation>
2032  V
2033  BlackoilModelBase<Grid, WellModel, Implementation>::
2034  fluidRsSat(const V& p,
2035  const V& satOil,
2036  const std::vector<int>& cells) const
2037  {
2038  return fluid_.rsSat(ADB::constant(p), ADB::constant(satOil), cells).value();
2039  }
2040 
2041 
2042 
2043 
2044 
2045  template <class Grid, class WellModel, class Implementation>
2046  ADB
2047  BlackoilModelBase<Grid, WellModel, Implementation>::
2048  fluidRsSat(const ADB& p,
2049  const ADB& satOil,
2050  const std::vector<int>& cells) const
2051  {
2052  return fluid_.rsSat(p, satOil, cells);
2053  }
2054 
2055 
2056 
2057 
2058 
2059  template <class Grid, class WellModel, class Implementation>
2060  V
2061  BlackoilModelBase<Grid, WellModel, Implementation>::
2062  fluidRvSat(const V& p,
2063  const V& satOil,
2064  const std::vector<int>& cells) const
2065  {
2066  return fluid_.rvSat(ADB::constant(p), ADB::constant(satOil), cells).value();
2067  }
2068 
2069 
2070 
2071 
2072 
2073  template <class Grid, class WellModel, class Implementation>
2074  ADB
2075  BlackoilModelBase<Grid, WellModel, Implementation>::
2076  fluidRvSat(const ADB& p,
2077  const ADB& satOil,
2078  const std::vector<int>& cells) const
2079  {
2080  return fluid_.rvSat(p, satOil, cells);
2081  }
2082 
2083 
2084 
2085 
2086 
2087  template <class Grid, class WellModel, class Implementation>
2088  ADB
2089  BlackoilModelBase<Grid, WellModel, Implementation>::
2090  poroMult(const ADB& p) const
2091  {
2092  const int n = p.size();
2093  if (rock_comp_props_ && rock_comp_props_->isActive()) {
2094  V pm(n);
2095  V dpm(n);
2096 #pragma omp parallel for schedule(static)
2097  for (int i = 0; i < n; ++i) {
2098  pm[i] = rock_comp_props_->poroMult(p.value()[i]);
2099  dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
2100  }
2101  ADB::M dpm_diag(dpm.matrix().asDiagonal());
2102  const int num_blocks = p.numBlocks();
2103  std::vector<ADB::M> jacs(num_blocks);
2104 #pragma omp parallel for schedule(dynamic)
2105  for (int block = 0; block < num_blocks; ++block) {
2106  fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]);
2107  }
2108  return ADB::function(std::move(pm), std::move(jacs));
2109  } else {
2110  return ADB::constant(V::Constant(n, 1.0));
2111  }
2112  }
2113 
2114 
2115 
2116 
2117 
2118  template <class Grid, class WellModel, class Implementation>
2119  ADB
2120  BlackoilModelBase<Grid, WellModel, Implementation>::
2121  transMult(const ADB& p) const
2122  {
2123  const int n = p.size();
2124  if (rock_comp_props_ && rock_comp_props_->isActive()) {
2125  V tm(n);
2126  V dtm(n);
2127 #pragma omp parallel for schedule(static)
2128  for (int i = 0; i < n; ++i) {
2129  tm[i] = rock_comp_props_->transMult(p.value()[i]);
2130  dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
2131  }
2132  ADB::M dtm_diag(dtm.matrix().asDiagonal());
2133  const int num_blocks = p.numBlocks();
2134  std::vector<ADB::M> jacs(num_blocks);
2135 #pragma omp parallel for schedule(dynamic)
2136  for (int block = 0; block < num_blocks; ++block) {
2137  fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]);
2138  }
2139  return ADB::function(std::move(tm), std::move(jacs));
2140  } else {
2141  return ADB::constant(V::Constant(n, 1.0));
2142  }
2143  }
2144 
2145 
2146 
2147 
2148 
2149  template <class Grid, class WellModel, class Implementation>
2150  void
2151  BlackoilModelBase<Grid, WellModel, Implementation>::
2152  classifyCondition(const ReservoirState& state)
2153  {
2154  using namespace Opm::AutoDiffGrid;
2155  const int nc = numCells(grid_);
2156  const int np = state.numPhases();
2157 
2158  const PhaseUsage& pu = fluid_.phaseUsage();
2159  const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
2160  if (active_[ Gas ]) {
2161  // Oil/Gas or Water/Oil/Gas system
2162  const V so = s.col(pu.phase_pos[ Oil ]);
2163  const V sg = s.col(pu.phase_pos[ Gas ]);
2164 
2165  for (V::Index c = 0, e = sg.size(); c != e; ++c) {
2166  if (so[c] > 0) { phaseCondition_[c].setFreeOil (); }
2167  if (sg[c] > 0) { phaseCondition_[c].setFreeGas (); }
2168  if (active_[ Water ]) { phaseCondition_[c].setFreeWater(); }
2169  }
2170  }
2171  else {
2172  // Water/Oil system
2173  assert (active_[ Water ]);
2174 
2175  const V so = s.col(pu.phase_pos[ Oil ]);
2176 
2177 
2178  for (V::Index c = 0, e = so.size(); c != e; ++c) {
2179  phaseCondition_[c].setFreeWater();
2180 
2181  if (so[c] > 0) { phaseCondition_[c].setFreeOil(); }
2182  }
2183  }
2184  }
2185 
2186 
2187 
2188 
2189 
2190  template <class Grid, class WellModel, class Implementation>
2191  void
2193  updatePrimalVariableFromState(const ReservoirState& state)
2194  {
2195  updatePhaseCondFromPrimalVariable(state);
2196  }
2197 
2198 
2199 
2200 
2201 
2203  template <class Grid, class WellModel, class Implementation>
2204  void
2206  updatePhaseCondFromPrimalVariable(const ReservoirState& state)
2207  {
2208  const int nc = Opm::AutoDiffGrid::numCells(grid_);
2209  isRs_ = V::Zero(nc);
2210  isRv_ = V::Zero(nc);
2211  isSg_ = V::Zero(nc);
2212 
2213  if (! (active_[Gas] && active_[Oil])) {
2214  // updatePhaseCondFromPrimarVariable() logic requires active gas and oil phase.
2215  phaseCondition_.assign(nc, PhasePresence());
2216  return;
2217  }
2218  for (int c = 0; c < nc; ++c) {
2219  phaseCondition_[c] = PhasePresence(); // No free phases.
2220  phaseCondition_[c].setFreeWater(); // Not necessary for property calculation usage.
2221  switch (state.hydroCarbonState()[c]) {
2222  case HydroCarbonState::GasAndOil:
2223  phaseCondition_[c].setFreeOil();
2224  phaseCondition_[c].setFreeGas();
2225  isSg_[c] = 1;
2226  break;
2227  case HydroCarbonState::OilOnly:
2228  phaseCondition_[c].setFreeOil();
2229  isRs_[c] = 1;
2230  break;
2231  case HydroCarbonState::GasOnly:
2232  phaseCondition_[c].setFreeGas();
2233  isRv_[c] = 1;
2234  break;
2235  default:
2236  OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << c << ": " << state.hydroCarbonState()[c]);
2237  }
2238  }
2239  }
2240 
2241 
2242 
2243 
2244 
2245  template <class Grid, class WellModel, class Implementation>
2246  void
2248  computeWellConnectionPressures(const SolutionState& state,
2249  const WellState& well_state)
2250  {
2251  asImpl().wellModel().computeWellConnectionPressures(state, well_state);
2252  }
2253 
2254 
2255 
2256 
2257 
2258  template <class Grid, class WellModel, class Implementation>
2259  std::vector<std::vector<double> >
2261  computeFluidInPlace(const ReservoirState& x,
2262  const std::vector<int>& fipnum)
2263  {
2264  using namespace Opm::AutoDiffGrid;
2265  const int nc = numCells(grid_);
2266  std::vector<ADB> saturation(3, ADB::null());
2267  const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, x.numPhases());
2268  const ADB pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
2269  const ADB temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
2270  saturation[Water] = active_[Water] ? ADB::constant(s.col(Water)) : ADB::null();
2271  saturation[Oil] = active_[Oil] ? ADB::constant(s.col(Oil)) : ADB::constant(V::Zero(nc));
2272  saturation[Gas] = active_[Gas] ? ADB::constant(s.col(Gas)) : ADB::constant(V::Zero(nc));
2273  const ADB rs = ADB::constant(Eigen::Map<const V>(& x.gasoilratio()[0], nc, 1));
2274  const ADB rv = ADB::constant(Eigen::Map<const V>(& x.rv()[0], nc, 1));
2275  const auto canonical_phase_pressures = computePressures(pressure, saturation[Water], saturation[Oil], saturation[Gas]);
2276  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
2277  const std::vector<PhasePresence> cond = phaseCondition();
2278 
2279  const ADB pv_mult = poroMult(pressure);
2280  const V& pv = geo_.poreVolume();
2281  const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
2282  for (int phase = 0; phase < maxnp; ++phase) {
2283  if (active_[ phase ]) {
2284  const int pos = pu.phase_pos[ phase ];
2285  const auto& b = asImpl().fluidReciprocFVF(phase, canonical_phase_pressures[phase], temperature, rs, rv, cond);
2286  sd_.fip[phase] = ((pv_mult * b * saturation[pos] * pv).value());
2287  }
2288  }
2289 
2290  if (active_[ Oil ] && active_[ Gas ]) {
2291  // Account for gas dissolved in oil and vaporized oil
2292  sd_.fip[SimulatorData::FIP_DISSOLVED_GAS] = rs.value() * sd_.fip[SimulatorData::FIP_LIQUID];
2293  sd_.fip[SimulatorData::FIP_VAPORIZED_OIL] = rv.value() * sd_.fip[SimulatorData::FIP_VAPOUR];
2294  }
2295 
2296  // For a parallel run this is just a local maximum and needs to be updated later
2297  int dims = *std::max_element(fipnum.begin(), fipnum.end());
2298  std::vector<std::vector<double> > values(dims);
2299  for (int i=0; i < dims; ++i) {
2300  values[i].resize(7, 0.0);
2301  }
2302 
2303  const V hydrocarbon = saturation[Oil].value() + saturation[Gas].value();
2304  V hcpv;
2305  V pres;
2306 
2307  if ( !isParallel() )
2308  {
2309  //Accumulate phases for each region
2310  for (int phase = 0; phase < maxnp; ++phase) {
2311  if (active_[ phase ]) {
2312  for (int c = 0; c < nc; ++c) {
2313  const int region = fipnum[c] - 1;
2314  if (region != -1) {
2315  values[region][phase] += sd_.fip[phase][c];
2316  }
2317  }
2318  }
2319  }
2320 
2321  //Accumulate RS and RV-volumes for each region
2322  if (active_[ Oil ] && active_[ Gas ]) {
2323  for (int c = 0; c < nc; ++c) {
2324  const int region = fipnum[c] - 1;
2325  if (region != -1) {
2326  values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
2327  values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
2328  }
2329  }
2330  }
2331 
2332 
2333  hcpv = V::Zero(dims);
2334  pres = V::Zero(dims);
2335 
2336  for (int c = 0; c < nc; ++c) {
2337  const int region = fipnum[c] - 1;
2338  if (region != -1) {
2339  hcpv[region] += pv[c] * hydrocarbon[c];
2340  pres[region] += pv[c] * pressure.value()[c];
2341  }
2342  }
2343 
2344  sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
2345  sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
2346 
2347  for (int c = 0; c < nc; ++c) {
2348  const int region = fipnum[c] - 1;
2349  if (region != -1) {
2350  sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
2351 
2352  //Compute hydrocarbon pore volume weighted average pressure.
2353  //If we have no hydrocarbon in region, use pore volume weighted average pressure instead
2354  if (hcpv[region] != 0) {
2355  sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
2356  } else {
2357  sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
2358  }
2359 
2360  values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
2361  values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
2362  }
2363  }
2364  }
2365  else
2366  {
2367 #if HAVE_MPI
2368  // mask[c] is 1 if we need to compute something in parallel
2369  const auto & pinfo =
2370  boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
2371  const auto& mask = pinfo.getOwnerMask();
2372  auto comm = pinfo.communicator();
2373  // Compute the global dims value and resize values accordingly.
2374  dims = comm.max(dims);
2375  values.resize(dims);
2376  for (int i=0; i < dims; ++i) {
2377  values[i].resize(7);
2378  std::fill(values[i].begin(), values[i].end(), 0.0);
2379  }
2380 
2381  //Accumulate phases for each region
2382  for (int phase = 0; phase < maxnp; ++phase) {
2383  for (int c = 0; c < nc; ++c) {
2384  const int region = fipnum[c] - 1;
2385  if (region != -1 && mask[c]) {
2386  values[region][phase] += sd_.fip[phase][c];
2387  }
2388  }
2389  }
2390 
2391  //Accumulate RS and RV-volumes for each region
2392  if (active_[ Oil ] && active_[ Gas ]) {
2393  for (int c = 0; c < nc; ++c) {
2394  const int region = fipnum[c] - 1;
2395  if (region != -1 && mask[c]) {
2396  values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
2397  values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
2398  }
2399  }
2400  }
2401 
2402  hcpv = V::Zero(dims);
2403  pres = V::Zero(dims);
2404 
2405  for (int c = 0; c < nc; ++c) {
2406  const int region = fipnum[c] - 1;
2407  if (region != -1 && mask[c]) {
2408  hcpv[region] += pv[c] * hydrocarbon[c];
2409  pres[region] += pv[c] * pressure.value()[c];
2410  }
2411  }
2412 
2413  comm.sum(hcpv.data(), hcpv.size());
2414  comm.sum(pres.data(), pres.size());
2415 
2416  sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
2417  sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
2418 
2419  for (int c = 0; c < nc; ++c) {
2420  const int region = fipnum[c] - 1;
2421  if (region != -1 && mask[c]) {
2422  sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
2423 
2424  if (hcpv[region] != 0) {
2425  sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
2426  } else {
2427  sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
2428  }
2429 
2430  values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
2431  values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
2432  }
2433  }
2434 
2435  // For the frankenstein branch we hopefully can turn values into a vanilla
2436  // std::vector<double>, use some index magic above, use one communication
2437  // to sum up the vector entries instead of looping over the regions.
2438  for(int reg=0; reg < dims; ++reg)
2439  {
2440  comm.sum(values[reg].data(), values[reg].size());
2441  }
2442 #else
2443  // This should never happen!
2444  OPM_THROW(std::logic_error, "HAVE_MPI should be defined if we are running in parallel");
2445 #endif
2446  }
2447 
2448  return values;
2449  }
2450 
2451 
2452 
2453 
2454 
2455  template <class Grid, class WellModel, class Implementation>
2456  void
2458  computeWellVoidageRates(const ReservoirState& reservoir_state,
2459  const WellState& well_state,
2460  std::vector<double>& well_voidage_rates,
2461  std::vector<double>& voidage_conversion_coeffs)
2462  {
2463  // TODO: for now, we store the voidage rates for all the production wells.
2464  // For injection wells, the rates are stored as zero.
2465  // We only store the conversion coefficients for all the injection wells.
2466  // Later, more delicate model will be implemented here.
2467  // And for the moment, group control can only work for serial running.
2468  const int nw = well_state.numWells();
2469  const int np = numPhases();
2470 
2471  const Wells* wells = asImpl().wellModel().wellsPointer();
2472 
2473  // we calculate the voidage rate for each well, that means the sum of all the phases.
2474  well_voidage_rates.resize(nw, 0);
2475  // store the conversion coefficients, while only for the use of injection wells.
2476  voidage_conversion_coeffs.resize(nw * np, 1.0);
2477 
2478  int global_number_wells = nw;
2479 
2480 #if HAVE_MPI
2481  if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
2482  {
2483  const auto& info =
2484  boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
2485  global_number_wells = info.communicator().sum(global_number_wells);
2486  if ( global_number_wells )
2487  {
2488  rate_converter_.defineState(reservoir_state, boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation()));
2489  }
2490  }
2491  else
2492 #endif
2493  {
2494  if ( global_number_wells )
2495  {
2496  rate_converter_.defineState(reservoir_state);
2497  }
2498  }
2499 
2500  std::vector<double> well_rates(np, 0.0);
2501  std::vector<double> convert_coeff(np, 1.0);
2502 
2503 
2504  if ( !well_voidage_rates.empty() ) {
2505  for (int w = 0; w < nw; ++w) {
2506  const bool is_producer = wells->type[w] == PRODUCER;
2507 
2508  // not sure necessary to change all the value to be positive
2509  if (is_producer) {
2510  std::transform(well_state.wellRates().begin() + np * w,
2511  well_state.wellRates().begin() + np * (w + 1),
2512  well_rates.begin(), std::negate<double>());
2513 
2514  // the average hydrocarbon conditions of the whole field will be used
2515  const int fipreg = 0; // Not considering FIP for the moment.
2516  const int well_cell_top = wells->well_cells[wells->well_connpos[w]];
2517  const int pvtreg = fluid_.cellPvtRegionIndex()[well_cell_top];
2518 
2519  rate_converter_.calcCoeff(fipreg, pvtreg, convert_coeff);
2520  well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(),
2521  convert_coeff.begin(), 0.0);
2522  } else {
2523  // TODO: Not sure whether will encounter situation with all zero rates
2524  // and whether it will cause problem here.
2525  std::copy(well_state.wellRates().begin() + np * w,
2526  well_state.wellRates().begin() + np * (w + 1),
2527  well_rates.begin());
2528  // the average hydrocarbon conditions of the whole field will be used
2529  const int fipreg = 0; // Not considering FIP for the moment.
2530  const int well_cell_top = wells->well_cells[wells->well_connpos[w]];
2531  const int pvtreg = fluid_.cellPvtRegionIndex()[well_cell_top];
2532  rate_converter_.calcCoeff(fipreg, pvtreg, convert_coeff);
2533  std::copy(convert_coeff.begin(), convert_coeff.end(),
2534  voidage_conversion_coeffs.begin() + np * w);
2535  }
2536  }
2537  }
2538  }
2539 
2540 
2541 
2542 
2543 
2544  template <class Grid, class WellModel, class Implementation>
2545  void
2547  applyVREPGroupControl(const ReservoirState& reservoir_state,
2548  WellState& well_state)
2549  {
2550  if (asImpl().wellModel().wellCollection()->havingVREPGroups()) {
2551  std::vector<double> well_voidage_rates;
2552  std::vector<double> voidage_conversion_coeffs;
2553  computeWellVoidageRates(reservoir_state, well_state, well_voidage_rates, voidage_conversion_coeffs);
2554  asImpl().wellModel().wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
2555 
2556  // for the wells under group control, update the currentControls for the well_state
2557  for (const WellNode* well_node : asImpl().wellModel().wellCollection()->getLeafNodes()) {
2558  if (well_node->isInjector() && !well_node->individualControl()) {
2559  const int well_index = well_node->selfIndex();
2560  well_state.currentControls()[well_index] = well_node->groupControlIndex();
2561  }
2562  }
2563  }
2564  }
2565 
2566 
2567 
2568 
2569 
2570  template <class Grid, class WellModel, class Implementation>
2571  void
2573  setupGroupControl(const ReservoirState& reservoir_state,
2574  WellState& well_state)
2575  {
2576  if (asImpl().wellModel().wellCollection()->requireWellPotentials()) {
2577  SolutionState state = asImpl().variableState(reservoir_state, well_state);
2578  asImpl().makeConstantState(state);
2579  asImpl().wellModel().computeWellConnectionPressures(state, well_state);
2580 
2581  const int np = numPhases();
2582  std::vector<ADB> b(np, ADB::null());
2583  std::vector<ADB> mob(np, ADB::null());
2584 
2585  const ADB& press = state.pressure;
2586  const ADB& temp = state.temperature;
2587  const ADB& rs = state.rs;
2588  const ADB& rv = state.rv;
2589 
2590  const std::vector<PhasePresence> cond = phaseCondition();
2591 
2592  const ADB pv_mult = poroMult(press);
2593  const ADB tr_mult = transMult(press);
2594  const std::vector<ADB> kr = asImpl().computeRelPerm(state);
2595 
2596  std::vector<ADB> mob_perfcells(np, ADB::null());
2597  std::vector<ADB> b_perfcells(np, ADB::null());
2598 
2599  for (int phase = 0; phase < np; ++phase) {
2600  if (active_[phase]) {
2601  const std::vector<int>& well_cells = asImpl().wellModel().wellOps().well_cells;
2602  const ADB mu = asImpl().fluidViscosity(canph_[phase], state.canonical_phase_pressures[canph_[phase]],
2603  temp, rs, rv, cond);
2604  mob[phase] = tr_mult * kr[canph_[phase]] / mu;
2605  mob_perfcells[phase] = subset(mob[phase], well_cells);
2606 
2607  b[phase] = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
2608  b_perfcells[phase] = subset(b[phase], well_cells);
2609  }
2610  }
2611 
2612  // well potentials for each well
2613  std::vector<double> well_potentials;
2614  asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, well_state, state, well_potentials);
2615  asImpl().wellModel().wellCollection()->setGuideRatesWithPotentials(asImpl().wellModel().wellsPointer(),
2616  fluid_.phaseUsage(), well_potentials);
2617  } // end of if
2618 
2619  applyVREPGroupControl(reservoir_state, well_state);
2620 
2621  if (asImpl().wellModel().wellCollection()->groupControlApplied()) {
2622  asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
2623  } else {
2624  asImpl().wellModel().wellCollection()->applyGroupControls();
2625 
2626  // the well collections do not have access to Well State, so the currentControls() of Well State need to
2627  // be updated based on the group control setup
2628  const int nw = wells().number_of_wells;
2629  for (int w = 0; w < nw; ++w) {
2630  const WellNode& well_node = asImpl().wellModel().wellCollection()->findWellNode(wells().name[w]);
2631  if (!well_node.individualControl()) {
2632  well_state.currentControls()[w] = well_node.groupControlIndex();
2633  }
2634  }
2635  } // end of else
2636  }
2637 
2638 } // namespace Opm
2639 
2640 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
int sizeNonLinear() const
The size (number of unknowns) of the nonlinear system of equations.
Definition: BlackoilModelBase_impl.hpp:347
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
void computeWellVoidageRates(const ReservoirState &reservoir_state, const WellState &well_state, std::vector< double > &well_voidage_rates, std::vector< double > &voidage_conversion_coeffs)
Function to compute the resevoir voidage for the production wells.
Definition: BlackoilModelBase_impl.hpp:2458
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
bool getConvergence(const SimulatorTimerInterface &timer, const int iteration)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv)...
Definition: BlackoilModelBase_impl.hpp:1721
SimulatorReport assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilModelBase_impl.hpp:758
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
AutoDiffMatrix M
Underlying type for jacobians.
Definition: AutoDiffBlock.hpp:101
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModelBase_impl.hpp:371
double relativeChange(const SimulationDataContainer &previous, const SimulationDataContainer &current) const
compute the relative change between to simulation states
Definition: BlackoilModelBase_impl.hpp:1593
int numPhases() const
Definition: BlackoilPropsAdFromDeck.cpp:225
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilModelBase_impl.hpp:1148
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilModelBase_impl.hpp:220
void afterStep(const SimulatorTimerInterface &timer, ReservoirState &reservoir_state, WellState &well_state)
Called once after each time step.
Definition: BlackoilModelBase_impl.hpp:333
virtual double currentStepLength() const =0
Current step length.
std::vector< std::vector< double > > computeFluidInPlace(const ReservoirState &x, const std::vector< int > &fipnum)
Compute fluid in place.
Definition: BlackoilModelBase_impl.hpp:2261
void updatePrimalVariableFromState(const ReservoirState &state)
update the primal variable for Sg, Rv or Rs.
Definition: BlackoilModelBase_impl.hpp:2193
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelBase_impl.hpp:383
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:319
WellModel & wellModel()
return the WellModel object
Definition: BlackoilModelBase.hpp:268
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
BlackoilModelBase(const ModelParameters &param, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const WellModel &well_model, 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: BlackoilModelBase_impl.hpp:101
void setThresholdPressures(const std::vector< double > &threshold_pressures_by_face)
Set threshold pressures that prevent or reduce flow.
Definition: BlackoilModelBase_impl.hpp:419
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:1552
void updateEquationsScaling()
Update the scaling factors for mass balance equations.
Definition: BlackoilModelBase_impl.hpp:920
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
int numMaterials() const
The number of active materials in the model.
Definition: BlackoilModelBase_impl.hpp:394
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition: BlackoilModelBase_impl.hpp:359
Definition: AutoDiffHelpers.hpp:623
bool isParallel() const
Return true if this is a parallel run.
Definition: BlackoilModelBase_impl.hpp:198
AutoDiffBlock< double > vertcatCollapseJacs(const std::vector< AutoDiffBlock< double > > &x)
Returns the vertical concatenation [ x[0]; x[1]; ...; x[n-1] ] of the inputs.
Definition: AutoDiffHelpers.hpp:530
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
const std::string & materialName(int material_index) const
The name of an active material in the model.
Definition: BlackoilModelBase_impl.hpp:406
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
void setupGroupControl(const ReservoirState &reservoir_state, WellState &well_state)
Set up the group control related at the beginning of each time step.
Definition: BlackoilModelBase_impl.hpp:2573
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
V solveJacobianSystem() const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilModelBase_impl.hpp:1140
void updatePhaseCondFromPrimalVariable(const ReservoirState &state)
Update the phaseCondition_ member based on the primalVariable_ member.
Definition: BlackoilModelBase_impl.hpp:2206
Eigen::Array< double, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Utility function to lessen code changes required elsewhere.
Definition: AutoDiffMatrix.hpp:708
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
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:719
SimulatorReport nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver, ReservoirState &reservoir_state, WellState &well_state)
Called once per nonlinear iteration.
Definition: BlackoilModelBase_impl.hpp:240