All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
BlackoilModelBase.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Statoil ASA.
4  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
5  Copyright 2015 NTNU
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  OPM is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with OPM. If not, see <http://www.gnu.org/licenses/>.
21 */
22 
23 #ifndef OPM_BLACKOILMODELBASE_HEADER_INCLUDED
24 #define OPM_BLACKOILMODELBASE_HEADER_INCLUDED
25 
26 #include <cassert>
27 
28 #include <opm/autodiff/AutoDiffBlock.hpp>
29 #include <opm/autodiff/AutoDiffHelpers.hpp>
30 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
31 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
32 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
33 #include <opm/autodiff/BlackoilModelEnums.hpp>
34 #include <opm/autodiff/VFPProperties.hpp>
36 #include <opm/autodiff/IterationReport.hpp>
37 #include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
38 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
39 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
40 #include <opm/core/simulator/SimulatorReport.hpp>
41 
42 #include <opm/common/data/SimulationDataContainer.hpp>
43 
44 #include <array>
45 
46 struct Wells;
47 
48 namespace Opm {
49 
50  class ParameterGroup;
51  class DerivedGeology;
52  class RockCompressibility;
53  class NewtonIterationBlackoilInterface;
54  class VFPProperties;
55 
59  template <class ConcreteModel>
60  struct ModelTraits;
61 
62 
75  template<class Grid, class WellModel, class Implementation>
77  {
78  public:
79  // --------- Types and enums ---------
80  typedef AutoDiffBlock<double> ADB;
81  typedef ADB::V V;
82  typedef ADB::M M;
83 
86  std::vector<ADB> accum; // Accumulations
87  ADB mflux; // Mass flux (surface conditions)
88  ADB b; // Reciprocal FVF
89  ADB mu; // Viscosities
90  ADB rho; // Densities
91  ADB kr; // Permeabilities
92  ADB dh; // Pressure drop across int. interfaces
93  ADB mob; // Phase mobility (per cell)
94  };
95 
96  struct SimulatorData : public Opm::FIPDataEnums {
97  SimulatorData(int num_phases);
98 
99  using Opm::FIPDataEnums :: FipId ;
100  using Opm::FIPDataEnums :: fipValues ;
101 
102  std::vector<ReservoirResidualQuant> rq;
103  ADB rsSat; // Saturated gas-oil ratio
104  ADB rvSat; // Saturated oil-gas ratio
105 
106  std::vector<double> soMax; // Maximum oil saturation
107 
108  std::vector<double> Pb; // Bubble point pressure
109  std::vector<double> Pd; // Dew point pressure
110 
111  //Hysteresis parameters
112  std::vector<double> krnswdc_ow;
113  std::vector<double> krnswdc_go;
114  std::vector<double> pcswmdc_ow;
115  std::vector<double> pcswmdc_go;
116 
117  std::array<V, fipValues> fip;
118  };
119 
120  typedef Opm::FIPData FIPDataType;
121 
122  typedef typename ModelTraits<Implementation>::ReservoirState ReservoirState;
123  typedef typename ModelTraits<Implementation>::WellState WellState;
124  typedef typename ModelTraits<Implementation>::ModelParameters ModelParameters;
125  typedef typename ModelTraits<Implementation>::SolutionState SolutionState;
126 
127  // For the conversion between the surface volume rate and resrevoir voidage rate
130 
131  // --------- Public methods ---------
132 
148  BlackoilModelBase(const ModelParameters& param,
149  const Grid& grid ,
150  const BlackoilPropsAdFromDeck& fluid,
151  const DerivedGeology& geo ,
152  const RockCompressibility* rock_comp_props,
153  const WellModel& well_model,
154  const NewtonIterationBlackoilInterface& linsolver,
155  std::shared_ptr< const EclipseState > eclState,
156  const bool has_disgas,
157  const bool has_vapoil,
158  const bool terminal_output);
159 
168  void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
169 
174  void prepareStep(const SimulatorTimerInterface& timer,
175  const ReservoirState& reservoir_state,
176  const WellState& well_state);
177 
187  template <class NonlinearSolverType>
188  SimulatorReport nonlinearIteration(const int iteration,
189  const SimulatorTimerInterface& timer,
190  NonlinearSolverType& nonlinear_solver,
191  ReservoirState& reservoir_state,
192  WellState& well_state);
193 
199  void afterStep(const SimulatorTimerInterface& timer,
200  ReservoirState& reservoir_state,
201  WellState& well_state);
202 
207  SimulatorReport
208  assemble(const ReservoirState& reservoir_state,
209  WellState& well_state,
210  const bool initial_assembly);
211 
216  std::vector<double> computeResidualNorms() const;
217 
219  // \return || u^n+1 - u^n || / || u^n+1 ||
220  double relativeChange( const SimulationDataContainer& previous, const SimulationDataContainer& current ) const;
221 
223  int sizeNonLinear() const;
224 
226  int linearIterationsLastSolve() const;
227 
230  V solveJacobianSystem() const;
231 
236  void updateState(const V& dx,
237  ReservoirState& reservoir_state,
238  WellState& well_state);
239 
241  bool isParallel() const;
242 
244  bool terminalOutputEnabled() const;
245 
250  bool getConvergence(const SimulatorTimerInterface& timer, const int iteration);
251 
253  int numPhases() const;
254 
258  int numMaterials() const;
259 
262  const std::string& materialName(int material_index) const;
263 
265  void updateEquationsScaling();
266 
268  WellModel& wellModel() { return well_model_; }
269  const WellModel& wellModel() const { return well_model_; }
270 
272  const SimulatorData& getSimulatorData(const SimulationDataContainer&) const {
273  return sd_;
274  }
275 
278  return FIPDataType( sd_.fip );
279  }
280 
285  std::vector<std::vector<double> >
286  computeFluidInPlace(const ReservoirState& x,
287  const std::vector<int>& fipnum);
288 
292  void computeWellVoidageRates(const ReservoirState& reservoir_state,
293  const WellState& well_state,
294  std::vector<double>& well_voidage_rates,
295  std::vector<double>& voidage_conversion_coeffs);
296 
297 
298  void applyVREPGroupControl(const ReservoirState& reservoir_state,
299  WellState& well_state);
300 
305  const SimulatorReport& failureReport() const
306  { return failureReport_; }
307 
308  protected:
309 
310  // --------- Types and enums ---------
311 
312  typedef Eigen::Array<double,
313  Eigen::Dynamic,
314  Eigen::Dynamic,
315  Eigen::RowMajor> DataBlock;
316 
317 
318  // --------- Data members ---------
319 
320  SimulatorReport failureReport_;
321  const Grid& grid_;
322  const BlackoilPropsAdFromDeck& fluid_;
323  const DerivedGeology& geo_;
324  const RockCompressibility* rock_comp_props_;
325  VFPProperties vfp_properties_;
326  const NewtonIterationBlackoilInterface& linsolver_;
327  // For each canonical phase -> true if active
328  const std::vector<bool> active_;
329  // Size = # active phases. Maps active -> canonical phase indices.
330  const std::vector<int> canph_;
331  const std::vector<int> cells_; // All grid cells
332  HelperOps ops_;
333  const bool has_disgas_;
334  const bool has_vapoil_;
335 
336  ModelParameters param_;
337  bool use_threshold_pressure_;
338  V threshold_pressures_by_connection_;
339 
340  mutable SimulatorData sd_;
341  std::vector<PhasePresence> phaseCondition_;
342 
343  // Well Model
344  WellModel well_model_;
345 
346  V isRs_;
347  V isRv_;
348  V isSg_;
349 
350  LinearisedBlackoilResidual residual_;
351 
356 
357  V pvdt_;
358  std::vector<std::string> material_name_;
359  std::vector<std::vector<double>> residual_norms_history_;
360  double current_relaxation_;
361  V dx_old_;
362 
363  // rate converter between the surface volume rates and reservoir voidage rates
364  RateConverterType rate_converter_;
365 
366  // --------- Protected methods ---------
367 
370  Implementation& asImpl()
371  {
372  return static_cast<Implementation&>(*this);
373  }
374 
377  const Implementation& asImpl() const
378  {
379  return static_cast<const Implementation&>(*this);
380  }
381 
383  const Wells& wells() const { return well_model_.wells(); }
384 
386  bool wellsActive() const { return well_model_.wellsActive(); }
387 
389  bool localWellsActive() const { return well_model_.localWellsActive(); }
390 
391  void
392  makeConstantState(SolutionState& state) const;
393 
394  SolutionState
395  variableState(const ReservoirState& x,
396  const WellState& xw) const;
397 
398  std::vector<V>
399  variableStateInitials(const ReservoirState& x,
400  const WellState& xw) const;
401  void
402  variableReservoirStateInitials(const ReservoirState& x,
403  std::vector<V>& vars0) const;
404 
405  std::vector<int>
406  variableStateIndices() const;
407 
408  SolutionState
409  variableStateExtractVars(const ReservoirState& x,
410  const std::vector<int>& indices,
411  std::vector<ADB>& vars) const;
412 
413  void
414  computeAccum(const SolutionState& state,
415  const int aix );
416 
417  void
418  assembleMassBalanceEq(const SolutionState& state);
419 
420 
421  SimulatorReport
422  solveWellEq(const std::vector<ADB>& mob_perfcells,
423  const std::vector<ADB>& b_perfcells,
424  const ReservoirState& reservoir_state,
425  SolutionState& state,
426  WellState& well_state);
427 
428  void
429  addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
430  const SolutionState& state,
431  const WellState& xw);
432 
433  bool getWellConvergence(const int iteration);
434 
435  bool isVFPActive() const;
436 
437  std::vector<ADB>
438  computePressures(const ADB& po,
439  const ADB& sw,
440  const ADB& so,
441  const ADB& sg) const;
442 
443  V
444  computeGasPressure(const V& po,
445  const V& sw,
446  const V& so,
447  const V& sg) const;
448 
449  std::vector<ADB>
450  computeRelPerm(const SolutionState& state) const;
451 
452  void
453  computeMassFlux(const int actph ,
454  const V& transi,
455  const ADB& kr ,
456  const ADB& mu ,
457  const ADB& rho ,
458  const ADB& p ,
459  const SolutionState& state );
460 
461  void applyThresholdPressures(ADB& dp);
462 
463  ADB
464  fluidViscosity(const int phase,
465  const ADB& p ,
466  const ADB& temp ,
467  const ADB& rs ,
468  const ADB& rv ,
469  const std::vector<PhasePresence>& cond) const;
470 
471  ADB
472  fluidReciprocFVF(const int phase,
473  const ADB& p ,
474  const ADB& temp ,
475  const ADB& rs ,
476  const ADB& rv ,
477  const std::vector<PhasePresence>& cond) const;
478 
479  ADB
480  fluidDensity(const int phase,
481  const ADB& b,
482  const ADB& rs,
483  const ADB& rv) const;
484 
485  V
486  fluidRsSat(const V& p,
487  const V& so,
488  const std::vector<int>& cells) const;
489 
490  ADB
491  fluidRsSat(const ADB& p,
492  const ADB& so,
493  const std::vector<int>& cells) const;
494 
495  V
496  fluidRvSat(const V& p,
497  const V& so,
498  const std::vector<int>& cells) const;
499 
500  ADB
501  fluidRvSat(const ADB& p,
502  const ADB& so,
503  const std::vector<int>& cells) const;
504 
505  ADB
506  poroMult(const ADB& p) const;
507 
508  ADB
509  transMult(const ADB& p) const;
510 
511  const std::vector<PhasePresence>
512  phaseCondition() const {return phaseCondition_;}
513 
514  void
515  classifyCondition(const ReservoirState& state);
516 
517 
520  void
521  updatePrimalVariableFromState(const ReservoirState& state);
522 
525  void
526  updatePhaseCondFromPrimalVariable(const ReservoirState& state);
527 
528  // TODO: added since the interfaces of the function are different
529  // TODO: for StandardWells and MultisegmentWells
530  void
531  computeWellConnectionPressures(const SolutionState& state,
532  const WellState& well_state);
533 
554  double
555  convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
556  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
557  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
558  std::vector<double>& R_sum,
559  std::vector<double>& maxCoeff,
560  std::vector<double>& B_avg,
561  std::vector<double>& maxNormWell,
562  int nc) const;
563 
565  void
566  setupGroupControl(const ReservoirState& reservoir_state,
567  WellState& well_state);
568 
569  double dpMaxRel() const { return param_.dp_max_rel_; }
570  double dbhpMaxRel() const {return param_.dbhp_max_rel_; }
571  double dsMax() const { return param_.ds_max_; }
572  double drMaxRel() const { return param_.dr_max_rel_; }
573  double maxResidualAllowed() const { return param_.max_residual_allowed_; }
574 
575  };
576 } // namespace Opm
577 
578 #include "BlackoilModelBase_impl.hpp"
579 
580 #endif // OPM_BLACKOILMODELBASE_HEADER_INCLUDED
Contains vectors and sparse matrices that represent subsets or operations on (AD or regular) vectors ...
Definition: AutoDiffHelpers.hpp:44
FIPDataType getFIPData() const
Return fluid-in-place data (for output functionality)
Definition: BlackoilModelBase.hpp:277
Definition: BlackoilModelEnums.hpp:66
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
A thin wrapper class that holds one VFPProdProperties and one VFPInjProperties object.
Definition: VFPProperties.hpp:37
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
Definition: BlackoilModelBase.hpp:96
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
const SimulatorData & getSimulatorData(const SimulationDataContainer &) const
Return reservoir simulation data (for output functionality)
Definition: BlackoilModelBase.hpp:272
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
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
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
bool localWellsActive() const
return true if wells are available on this process
Definition: BlackoilModelBase.hpp:389
Facility for converting component rates at surface conditions to phase (voidage) rates at reservoir c...
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
const Wells & wells() const
return the Well struct in the WellModel
Definition: BlackoilModelBase.hpp:383
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilModelBase_impl.hpp:220
Definition: BlackoilModelEnums.hpp:51
void afterStep(const SimulatorTimerInterface &timer, ReservoirState &reservoir_state, WellState &well_state)
Called once after each time step.
Definition: BlackoilModelBase_impl.hpp:333
std::vector< std::vector< double > > computeFluidInPlace(const ReservoirState &x, const std::vector< int > &fipnum)
Compute fluid in place.
Definition: BlackoilModelBase_impl.hpp:2261
int global_nc_
The number of cells of the global grid.
Definition: BlackoilModelBase.hpp:355
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
Definition: BlackoilModelBase.hpp:84
Traits to encapsulate the types used by classes using or extending this model.
Definition: BlackoilModelBase.hpp:60
WellModel & wellModel()
return the WellModel object
Definition: BlackoilModelBase.hpp:268
Residual structure of the fully implicit solver.
Definition: LinearisedBlackoilResidual.hpp:47
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
bool terminal_output_
Whether we print something to std::cout.
Definition: BlackoilModelBase.hpp:353
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
bool isParallel() const
Return true if this is a parallel run.
Definition: BlackoilModelBase_impl.hpp:198
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
Implementation & asImpl()
Access the most-derived class used for static polymorphism (CRTP).
Definition: BlackoilModelBase.hpp:370
const std::string & materialName(int material_index) const
The name of an active material in the model.
Definition: BlackoilModelBase_impl.hpp:406
const SimulatorReport & failureReport() const
return the statistics if the nonlinearIteration() method failed.
Definition: BlackoilModelBase.hpp:305
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< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
bool wellsActive() const
return true if wells are available in the reservoir
Definition: BlackoilModelBase.hpp:386
const Implementation & asImpl() const
Access the most-derived class used for static polymorphism (CRTP).
Definition: BlackoilModelBase.hpp:377
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