28 #ifndef EWOMS_NCP_MODEL_HH
29 #define EWOMS_NCP_MODEL_HH
31 #include <opm/material/densead/Math.hpp>
50 #include <opm/common/Valgrind.hpp>
51 #include <opm/common/ErrorMacros.hpp>
52 #include <opm/common/Exceptions.hpp>
54 #include <dune/common/fvector.hh>
62 template <
class TypeTag>
67 namespace Properties {
196 template <
class TypeTag>
198 :
public MultiPhaseBaseModel<TypeTag>
200 typedef MultiPhaseBaseModel<TypeTag> ParentType;
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;
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) };
219 typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
221 typedef Opm::MathToolbox<Evaluation> Toolbox;
227 NcpModel(Simulator& simulator)
228 : ParentType(simulator)
238 DiffusionModule::registerParameters();
239 EnergyModule::registerParameters();
258 minActivityCoeff_.resize(this->numGridDof());
259 std::fill(minActivityCoeff_.begin(), minActivityCoeff_.end(), 1.0);
274 if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
277 std::ostringstream oss;
278 if (pvIdx == pressure0Idx)
279 oss <<
"pressure_" << FluidSystem::phaseName(0);
280 else if (saturation0Idx <= pvIdx && pvIdx < saturation0Idx + (numPhases - 1))
281 oss <<
"saturation_" << FluidSystem::phaseName(pvIdx - saturation0Idx);
282 else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents)
283 oss <<
"fugacity^" << FluidSystem::componentName(pvIdx - fugacity0Idx);
296 if (!(s = EnergyModule::eqName(eqIdx)).empty())
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(eqIdx - ncp0EqIdx);
315 ParentType::updateBegin();
320 for (
unsigned dofIdx = 0; dofIdx < this->numGridDof(); ++ dofIdx) {
321 if (this->isLocalDof(dofIdx)) {
323 this->solution(0)[dofIdx][Indices::pressure0Idx];
334 for (
unsigned dofIdx = 0; dofIdx < elemCtx.numDof(0); ++dofIdx) {
335 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
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, 0).fluidState();
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]);
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!");
361 Scalar tmp = EnergyModule::primaryVarWeight(*
this, globalDofIdx, pvIdx);
366 else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
368 unsigned compIdx = pvIdx - fugacity0Idx;
369 assert(0 <= compIdx && compIdx <= numComponents);
371 Opm::Valgrind::CheckDefined(minActivityCoeff_[globalDofIdx][compIdx]);
372 static const Scalar fugacityBaseWeight =
374 result = fugacityBaseWeight / minActivityCoeff_[globalDofIdx][compIdx];
376 else if (Indices::pressure0Idx == pvIdx) {
377 static const Scalar pressureBaseWeight =
GET_PROP_VALUE(TypeTag, NcpPressureBaseWeight);
378 result = pressureBaseWeight / referencePressure_;
382 unsigned phaseIdx = pvIdx - saturation0Idx;
383 assert(0 <= phaseIdx && phaseIdx < numPhases - 1);
387 static const Scalar saturationsBaseWeight =
389 result = saturationsBaseWeight;
392 assert(std::isfinite(result));
401 Scalar
eqWeight(
unsigned globalDofIdx,
unsigned eqIdx)
const
403 Scalar tmp = EnergyModule::eqWeight(*
this, globalDofIdx, eqIdx);
408 else if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
412 unsigned compIdx = eqIdx - Indices::conti0EqIdx;
413 assert(0 <= compIdx && compIdx <= numComponents);
416 return FluidSystem::molarMass(compIdx);
427 {
return minActivityCoeff_[globalDofIdx][compIdx]; }
432 void registerOutputModules_()
434 ParentType::registerOutputModules_();
443 mutable Scalar referencePressure_;
444 mutable std::vector<ComponentVector> minActivityCoeff_;
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition: ncpmodel.hh:313
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition: ncpmodel.hh:293
#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
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
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
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
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
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
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
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
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...
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
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
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition: ncpmodel.hh:359
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
Returns the relative weight of an equation.
Definition: ncpmodel.hh:401
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.
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition: ncpmodel.hh:271