ncpmodel.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_NCP_MODEL_HH
29 #define EWOMS_NCP_MODEL_HH
30 
31 #include <opm/material/densead/Math.hpp>
32 
33 #include "ncpproperties.hh"
34 #include "ncplocalresidual.hh"
36 #include "ncpprimaryvariables.hh"
37 #include "ncpboundaryratevector.hh"
38 #include "ncpratevector.hh"
40 #include "ncpnewtonmethod.hh"
41 #include "ncpindices.hh"
42 
49 
50 #include <opm/common/Valgrind.hpp>
51 #include <opm/common/ErrorMacros.hpp>
52 #include <opm/common/Exceptions.hpp>
53 
54 #include <dune/common/fvector.hh>
55 
56 #include <sstream>
57 #include <string>
58 #include <vector>
59 #include <array>
60 
61 namespace Ewoms {
62 template <class TypeTag>
63 class NcpModel;
64 }
65 
66 namespace Ewoms {
67 namespace Properties {
72  VtkComposition,
73  VtkEnergy,
74  VtkDiffusion));
75 
78  LocalResidual,
80 
83 
86 
89 
91 SET_BOOL_PROP(NcpModel, EnableEnergy, false);
92 
94 SET_BOOL_PROP(NcpModel, EnableDiffusion, false);
95 
98 
101 
104 
107 
110 
113 
115 SET_SCALAR_PROP(NcpModel, NcpPressureBaseWeight, 1.0);
117 SET_SCALAR_PROP(NcpModel, NcpSaturationsBaseWeight, 1.0);
119 SET_SCALAR_PROP(NcpModel, NcpFugacitiesBaseWeight, 1.0e-6);
120 
121 } // namespace Properties
122 
196 template <class TypeTag>
197 class NcpModel
198  : public MultiPhaseBaseModel<TypeTag>
199 {
200  typedef MultiPhaseBaseModel<TypeTag> ParentType;
201 
202  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
203  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
204  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
205  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
206  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
207  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
208 
209  enum { numPhases = FluidSystem::numPhases };
210  enum { numComponents = FluidSystem::numComponents };
211  enum { fugacity0Idx = Indices::fugacity0Idx };
212  enum { pressure0Idx = Indices::pressure0Idx };
213  enum { saturation0Idx = Indices::saturation0Idx };
214  enum { conti0EqIdx = Indices::conti0EqIdx };
215  enum { ncp0EqIdx = Indices::ncp0EqIdx };
216  enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
217  enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
218 
219  typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
220 
221  typedef Opm::MathToolbox<Evaluation> Toolbox;
222 
223  typedef Ewoms::EnergyModule<TypeTag, enableEnergy> EnergyModule;
224  typedef Ewoms::DiffusionModule<TypeTag, enableDiffusion> DiffusionModule;
225 
226 public:
227  NcpModel(Simulator& simulator)
228  : ParentType(simulator)
229  {}
230 
234  static void registerParameters()
235  {
237 
238  DiffusionModule::registerParameters();
239  EnergyModule::registerParameters();
240 
241  // register runtime parameters of the VTK output modules
243 
244  if (enableDiffusion)
246 
247  if (enableEnergy)
249  }
250 
254  void finishInit()
255  {
257 
258  minActivityCoeff_.resize(this->numGridDof());
259  std::fill(minActivityCoeff_.begin(), minActivityCoeff_.end(), 1.0);
260  }
261 
265  static std::string name()
266  { return "ncp"; }
267 
271  std::string primaryVarName(unsigned pvIdx) const
272  {
273  std::string s;
274  if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
275  return s;
276 
277  std::ostringstream oss;
278  if (pvIdx == pressure0Idx)
279  oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
280  else if (saturation0Idx <= pvIdx && pvIdx < saturation0Idx + (numPhases - 1))
281  oss << "saturation_" << FluidSystem::phaseName(/*phaseIdx=*/pvIdx - saturation0Idx);
282  else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents)
283  oss << "fugacity^" << FluidSystem::componentName(pvIdx - fugacity0Idx);
284  else
285  assert(false);
286 
287  return oss.str();
288  }
289 
293  std::string eqName(unsigned eqIdx) const
294  {
295  std::string s;
296  if (!(s = EnergyModule::eqName(eqIdx)).empty())
297  return s;
298 
299  std::ostringstream oss;
300  if (conti0EqIdx <= eqIdx && eqIdx < conti0EqIdx + numComponents)
301  oss << "continuity^" << FluidSystem::componentName(eqIdx - conti0EqIdx);
302  else if (ncp0EqIdx <= eqIdx && eqIdx < ncp0EqIdx + numPhases)
303  oss << "ncp_" << FluidSystem::phaseName(/*phaseIdx=*/eqIdx - ncp0EqIdx);
304  else
305  assert(false);
306 
307  return oss.str();
308  }
309 
313  void updateBegin()
314  {
315  ParentType::updateBegin();
316 
317  // find the a reference pressure. The first degree of freedom
318  // might correspond to non-interior entities which would lead
319  // to an undefined value, so we have to iterate...
320  for (unsigned dofIdx = 0; dofIdx < this->numGridDof(); ++ dofIdx) {
321  if (this->isLocalDof(dofIdx)) {
322  referencePressure_ =
323  this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
324  break;
325  }
326  }
327  }
328 
332  void updatePVWeights(const ElementContext& elemCtx) const
333  {
334  for (unsigned dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) {
335  unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
336 
337  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
338  minActivityCoeff_[globalIdx][compIdx] = 1e100;
339  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
340  const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();
341 
342  minActivityCoeff_[globalIdx][compIdx] =
343  std::min(minActivityCoeff_[globalIdx][compIdx],
344  Toolbox::value(fs.fugacityCoefficient(phaseIdx, compIdx))
345  * Toolbox::value(fs.pressure(phaseIdx)));
346  Opm::Valgrind::CheckDefined(minActivityCoeff_[globalIdx][compIdx]);
347  }
348  if (minActivityCoeff_[globalIdx][compIdx] <= 0)
349  OPM_THROW(Opm::NumericalProblem,
350  "The minumum activity coefficient for component " << compIdx
351  << " on DOF " << globalIdx << " is negative or zero!");
352  }
353  }
354  }
355 
359  Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
360  {
361  Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
362  Scalar result;
363  if (tmp > 0)
364  // energy related quantity
365  result = tmp;
366  else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
367  // component fugacity
368  unsigned compIdx = pvIdx - fugacity0Idx;
369  assert(0 <= compIdx && compIdx <= numComponents);
370 
371  Opm::Valgrind::CheckDefined(minActivityCoeff_[globalDofIdx][compIdx]);
372  static const Scalar fugacityBaseWeight =
373  GET_PROP_VALUE(TypeTag, NcpFugacitiesBaseWeight);
374  result = fugacityBaseWeight / minActivityCoeff_[globalDofIdx][compIdx];
375  }
376  else if (Indices::pressure0Idx == pvIdx) {
377  static const Scalar pressureBaseWeight = GET_PROP_VALUE(TypeTag, NcpPressureBaseWeight);
378  result = pressureBaseWeight / referencePressure_;
379  }
380  else {
381 #ifndef NDEBUG
382  unsigned phaseIdx = pvIdx - saturation0Idx;
383  assert(0 <= phaseIdx && phaseIdx < numPhases - 1);
384 #endif
385 
386  // saturation
387  static const Scalar saturationsBaseWeight =
388  GET_PROP_VALUE(TypeTag, NcpSaturationsBaseWeight);
389  result = saturationsBaseWeight;
390  }
391 
392  assert(std::isfinite(result));
393  assert(result > 0);
394 
395  return result;
396  }
397 
401  Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
402  {
403  Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
404  if (tmp > 0)
405  // an energy related equation
406  return tmp;
407  // an NCP
408  else if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
409  return 1.0;
410 
411  // a mass conservation equation
412  unsigned compIdx = eqIdx - Indices::conti0EqIdx;
413  assert(0 <= compIdx && compIdx <= numComponents);
414 
415  // make all kg equal
416  return FluidSystem::molarMass(compIdx);
417  }
418 
426  Scalar minActivityCoeff(unsigned globalDofIdx, unsigned compIdx) const
427  { return minActivityCoeff_[globalDofIdx][compIdx]; }
428 
432  void registerOutputModules_()
433  {
434  ParentType::registerOutputModules_();
435 
436  this->addOutputModule(new Ewoms::VtkCompositionModule<TypeTag>(this->simulator_));
437  if (enableDiffusion)
438  this->addOutputModule(new Ewoms::VtkDiffusionModule<TypeTag>(this->simulator_));
439  if (enableEnergy)
440  this->addOutputModule(new Ewoms::VtkEnergyModule<TypeTag>(this->simulator_));
441  }
442 
443  mutable Scalar referencePressure_;
444  mutable std::vector<ComponentVector> minActivityCoeff_;
445 };
446 
447 } // namespace Ewoms
448 
449 #endif
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition: ncpmodel.hh:313
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Details needed to calculate the local residual in the compositional multi-phase NCP-model ...
Definition: ncplocalresidual.hh:46
Implements a boundary vector for the fully implicit compositional multi-phase NCP model...
void finishInit()
Apply the initial conditions to the model.
Definition: ncpmodel.hh:254
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
Returns the relative weight of an equation.
Definition: ncpmodel.hh:401
Definition: baseauxiliarymodule.hh:37
This template class represents the extensive quantities of the compositional NCP model.
Definition: ncpextensivequantities.hh:47
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkenergymodule.hh:101
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition: ncpmodel.hh:234
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkcompositionmodule.hh:105
Provides the auxiliary methods required for consideration of the energy equation. ...
Definition: energymodule.hh:59
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
A compositional multi-phase model based on non-linear complementarity functions. ...
Definition: ncpmodel.hh:63
Scalar minActivityCoeff(unsigned globalDofIdx, unsigned compIdx) const
Returns the smallest activity coefficient of a component for the most current solution at a vertex...
Definition: ncpmodel.hh:426
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
void updatePVWeights(const ElementContext &elemCtx) const
Update the weights of all primary variables within an element given the complete set of intensive qua...
Definition: ncpmodel.hh:332
VTK output module for quantities which make sense for models which incorperate molecular diffusion...
Implements a boundary vector for the fully implicit compositional multi-phase NCP model...
Definition: ncpboundaryratevector.hh:45
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
Represents the primary variables used by the compositional multi-phase NCP model. ...
Definition: ncpprimaryvariables.hh:54
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition: ncpintensivequantities.hh:53
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition: ncpmodel.hh:293
static std::string name()
Definition: ncpmodel.hh:265
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkdiffusionmodule.hh:101
VTK output module for quantities which make sense for models which assume thermal equilibrium...
void finishInit()
Apply the initial conditions to the model.
Definition: multiphasebasemodel.hh:159
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition: ncpmodel.hh:359
Declares the properties required for the NCP compositional multi-phase model.
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:53
Represents the primary variables used by the compositional multi-phase NCP model. ...
The primary variable and equation indices for the compositional multi-phase NCP model.
Definition: ncpindices.hh:43
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition: ncpmodel.hh:271
VTK output module for the fluid composition.
The multi-dimensional Newton method.
Definition: newtonmethod.hh:56
Implements a vector representing mass, molar or volumetric rates.
Definition: ncpratevector.hh:48
VTK output module for quantities which make sense for models which incorperate molecular diffusion...
Definition: vtkdiffusionmodule.hh:73
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition: multiphasebasemodel.hh:147
A Newton solver specific to the NCP model.
Definition: ncpnewtonmethod.hh:47
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
Definition: multiphasebasemodel.hh:47
Implements a vector representing mass, molar or volumetric rates.
Details needed to calculate the local residual in the compositional multi-phase NCP-model ...
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
VTK output module for quantities which make sense for models which assume thermal equilibrium...
Definition: vtkenergymodule.hh:73
The base class for the problems of ECFV discretizations which deal with a multi-phase flow through a ...
Definition: multiphasebaseproblem.hh:59
A Newton solver specific to the NCP model.
VTK output module for the fluid composition.
Definition: vtkcompositionmodule.hh:77
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
Classes required for molecular diffusion.
This template class represents the extensive quantities of the compositional NCP model.
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
The primary variable and equation indices for the compositional multi-phase NCP model.