My Project
Loading...
Searching...
No Matches
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"
38#include "ncpratevector.hh"
40#include "ncpnewtonmethod.hh"
41#include "ncpindices.hh"
42
43#include <opm/common/Exceptions.hpp>
44
51
52#include <opm/material/common/Valgrind.hpp>
53
54#include <dune/common/fvector.hh>
55
56#include <sstream>
57#include <string>
58#include <vector>
59
60namespace Opm {
61template <class TypeTag>
62class NcpModel;
63}
64
65namespace Opm::Properties {
66
67namespace TTag {
71struct NcpModel { using InheritsFrom = std::tuple<MultiPhaseBaseModel>; };
72} // namespace TTag
73
75template<class TypeTag>
76struct LocalResidual<TypeTag, TTag::NcpModel> { using type = NcpLocalResidual<TypeTag>; };
77
79template<class TypeTag>
80struct NewtonMethod<TypeTag, TTag::NcpModel> { using type = NcpNewtonMethod<TypeTag>; };
81
83template<class TypeTag>
84struct Model<TypeTag, TTag::NcpModel> { using type = NcpModel<TypeTag>; };
85
87template<class TypeTag>
88struct BaseProblem<TypeTag, TTag::NcpModel> { using type = MultiPhaseBaseProblem<TypeTag>; };
89
91template<class TypeTag>
92struct EnableEnergy<TypeTag, TTag::NcpModel> { static constexpr bool value = false; };
93
95template<class TypeTag>
96struct EnableDiffusion<TypeTag, TTag::NcpModel> { static constexpr bool value = false; };
97
99template<class TypeTag>
100struct RateVector<TypeTag, TTag::NcpModel> { using type = NcpRateVector<TypeTag>; };
101
103template<class TypeTag>
105
107template<class TypeTag>
108struct PrimaryVariables<TypeTag, TTag::NcpModel> { using type = NcpPrimaryVariables<TypeTag>; };
109
111template<class TypeTag>
113
115template<class TypeTag>
117
119template<class TypeTag>
120struct Indices<TypeTag, TTag::NcpModel> { using type = NcpIndices<TypeTag, 0>; };
121
123template<class TypeTag>
124struct NcpPressureBaseWeight<TypeTag, TTag::NcpModel>
125{
126 using type = GetPropType<TypeTag, Scalar>;
127 static constexpr type value = 1.0;
128};
130template<class TypeTag>
131struct NcpSaturationsBaseWeight<TypeTag, TTag::NcpModel>
132{
133 using type = GetPropType<TypeTag, Scalar>;
134 static constexpr type value = 1.0;
135};
137template<class TypeTag>
138struct NcpFugacitiesBaseWeight<TypeTag, TTag::NcpModel>
139{
140 using type = GetPropType<TypeTag, Scalar>;
141 static constexpr type value = 1.0e-6;
142};
143
144} // namespace Opm::Properties
145
146namespace Opm {
147
222template <class TypeTag>
224 : public MultiPhaseBaseModel<TypeTag>
225{
227
234
235 enum { numPhases = FluidSystem::numPhases };
236 enum { numComponents = FluidSystem::numComponents };
237 enum { fugacity0Idx = Indices::fugacity0Idx };
238 enum { pressure0Idx = Indices::pressure0Idx };
239 enum { saturation0Idx = Indices::saturation0Idx };
240 enum { conti0EqIdx = Indices::conti0EqIdx };
241 enum { ncp0EqIdx = Indices::ncp0EqIdx };
242 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
243 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
244
245 using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
246
247 using Toolbox = MathToolbox<Evaluation>;
248
251
252public:
253 NcpModel(Simulator& simulator)
254 : ParentType(simulator)
255 {}
256
260 static void registerParameters()
261 {
263
264 DiffusionModule::registerParameters();
265 EnergyModule::registerParameters();
266
267 // register runtime parameters of the VTK output modules
269
270 if (enableDiffusion)
272
273 if (enableEnergy)
275 }
276
281 {
282 ParentType::finishInit();
283
284 minActivityCoeff_.resize(this->numGridDof());
285 std::fill(minActivityCoeff_.begin(), minActivityCoeff_.end(), 1.0);
286 }
287
288 void adaptGrid()
289 {
290 ParentType::adaptGrid();
291 minActivityCoeff_.resize(this->numGridDof());
292 }
293
297 static std::string name()
298 { return "ncp"; }
299
303 std::string primaryVarName(unsigned pvIdx) const
304 {
305 std::string s;
306 if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
307 return s;
308
309 std::ostringstream oss;
310 if (pvIdx == pressure0Idx)
311 oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
312 else if (saturation0Idx <= pvIdx && pvIdx < saturation0Idx + (numPhases - 1))
313 oss << "saturation_" << FluidSystem::phaseName(/*phaseIdx=*/pvIdx - saturation0Idx);
314 else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents)
315 oss << "fugacity^" << FluidSystem::componentName(pvIdx - fugacity0Idx);
316 else
317 assert(false);
318
319 return oss.str();
320 }
321
325 std::string eqName(unsigned eqIdx) const
326 {
327 std::string s;
328 if (!(s = EnergyModule::eqName(eqIdx)).empty())
329 return s;
330
331 std::ostringstream oss;
332 if (conti0EqIdx <= eqIdx && eqIdx < conti0EqIdx + numComponents)
333 oss << "continuity^" << FluidSystem::componentName(eqIdx - conti0EqIdx);
334 else if (ncp0EqIdx <= eqIdx && eqIdx < ncp0EqIdx + numPhases)
335 oss << "ncp_" << FluidSystem::phaseName(/*phaseIdx=*/eqIdx - ncp0EqIdx);
336 else
337 assert(false);
338
339 return oss.str();
340 }
341
346 {
347 ParentType::updateBegin();
348
349 // find the a reference pressure. The first degree of freedom
350 // might correspond to non-interior entities which would lead
351 // to an undefined value, so we have to iterate...
352 for (unsigned dofIdx = 0; dofIdx < this->numGridDof(); ++ dofIdx) {
353 if (this->isLocalDof(dofIdx)) {
354 referencePressure_ =
355 this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
356 break;
357 }
358 }
359 }
360
364 void updatePVWeights(const ElementContext& elemCtx) const
365 {
366 for (unsigned dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) {
367 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
368
369 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
370 minActivityCoeff_[globalIdx][compIdx] = 1e100;
371 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
372 const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();
373
374 minActivityCoeff_[globalIdx][compIdx] =
375 std::min(minActivityCoeff_[globalIdx][compIdx],
376 Toolbox::value(fs.fugacityCoefficient(phaseIdx, compIdx))
377 * Toolbox::value(fs.pressure(phaseIdx)));
378 Valgrind::CheckDefined(minActivityCoeff_[globalIdx][compIdx]);
379 }
380 if (minActivityCoeff_[globalIdx][compIdx] <= 0)
381 throw NumericalProblem("The minimum activity coefficient for component "+std::to_string(compIdx)
382 +" on DOF "+std::to_string(globalIdx)+" is negative or zero!");
383 }
384 }
385 }
386
390 Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
391 {
392 Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
393 Scalar result;
394 if (tmp > 0)
395 // energy related quantity
396 result = tmp;
397 else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
398 // component fugacity
399 unsigned compIdx = pvIdx - fugacity0Idx;
400 assert(compIdx <= numComponents);
401
402 Valgrind::CheckDefined(minActivityCoeff_[globalDofIdx][compIdx]);
403 static const Scalar fugacityBaseWeight =
405 result = fugacityBaseWeight / minActivityCoeff_[globalDofIdx][compIdx];
406 }
407 else if (Indices::pressure0Idx == pvIdx) {
408 static const Scalar pressureBaseWeight = getPropValue<TypeTag, Properties::NcpPressureBaseWeight>();
409 result = pressureBaseWeight / referencePressure_;
410 }
411 else {
412#ifndef NDEBUG
413 unsigned phaseIdx = pvIdx - saturation0Idx;
414 assert(phaseIdx < numPhases - 1);
415#endif
416
417 // saturation
418 static const Scalar saturationsBaseWeight =
420 result = saturationsBaseWeight;
421 }
422
423 assert(std::isfinite(result));
424 assert(result > 0);
425
426 return result;
427 }
428
432 Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
433 {
434 Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
435 if (tmp > 0)
436 // an energy related equation
437 return tmp;
438 // an NCP
439 else if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
440 return 1.0;
441
442 // a mass conservation equation
443 unsigned compIdx = eqIdx - Indices::conti0EqIdx;
444 assert(compIdx <= numComponents);
445
446 // make all kg equal
447 return FluidSystem::molarMass(compIdx);
448 }
449
457 Scalar minActivityCoeff(unsigned globalDofIdx, unsigned compIdx) const
458 { return minActivityCoeff_[globalDofIdx][compIdx]; }
459
463 void registerOutputModules_()
464 {
465 ParentType::registerOutputModules_();
466
467 this->addOutputModule(new VtkCompositionModule<TypeTag>(this->simulator_));
468 if (enableDiffusion)
469 this->addOutputModule(new VtkDiffusionModule<TypeTag>(this->simulator_));
470 if (enableEnergy)
471 this->addOutputModule(new VtkEnergyModule<TypeTag>(this->simulator_));
472 }
473
474 mutable Scalar referencePressure_;
475 mutable std::vector<ComponentVector> minActivityCoeff_;
476};
477
478} // namespace Opm
479
480#endif
Provides the auxiliary methods required for consideration of the diffusion equation.
Provides the auxiliary methods required for consideration of the energy equation.
Definition energymodule.hh:50
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
Definition multiphasebasemodel.hh:153
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition multiphasebasemodel.hh:179
The base class for the problems of ECFV discretizations which deal with a multi-phase flow through a ...
Definition multiphasebaseproblem.hh:60
Implements a boundary vector for the fully implicit compositional multi-phase NCP model.
Definition ncpboundaryratevector.hh:46
This template class represents the extensive quantities of the compositional NCP model.
Definition ncpextensivequantities.hh:51
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition ncpintensivequantities.hh:58
Details needed to calculate the local residual in the compositional multi-phase NCP-model .
Definition ncplocalresidual.hh:47
A compositional multi-phase model based on non-linear complementarity functions.
Definition ncpmodel.hh:225
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:364
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition ncpmodel.hh:325
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition ncpmodel.hh:303
static std::string name()
Definition ncpmodel.hh:297
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition ncpmodel.hh:390
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:457
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition ncpmodel.hh:260
void finishInit()
Apply the initial conditions to the model.
Definition ncpmodel.hh:280
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition ncpmodel.hh:345
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
Returns the relative weight of an equation.
Definition ncpmodel.hh:432
A Newton solver specific to the NCP model.
Definition ncpnewtonmethod.hh:55
Represents the primary variables used by the compositional multi-phase NCP model.
Definition ncpprimaryvariables.hh:55
Implements a vector representing mass, molar or volumetric rates.
Definition ncpratevector.hh:51
VTK output module for the fluid composition.
Definition vtkcompositionmodule.hpp:57
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkcompositionmodule.hpp:86
VTK output module for quantities which make sense for models which incorperate molecular diffusion.
Definition vtkdiffusionmodule.hpp:58
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkdiffusionmodule.hpp:87
VTK output module for quantities which make sense for models which assume thermal equilibrium.
Definition vtkenergymodule.hpp:58
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkenergymodule.hpp:86
Classes required for molecular diffusion.
Contains the classes required to consider energy as a conservation quantity in a multi-phase module.
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Implements a boundary vector for the fully implicit compositional multi-phase NCP model.
This template class represents the extensive quantities of the compositional NCP model.
The primary variable and equation indices for the compositional multi-phase NCP model.
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Details needed to calculate the local residual in the compositional multi-phase NCP-model .
A Newton solver specific to the NCP model.
Represents the primary variables used by the compositional multi-phase NCP model.
Declares the properties required for the NCP compositional multi-phase model.
Implements a vector representing mass, molar or volumetric rates.
The primary variable and equation indices for the compositional multi-phase NCP model.
Definition ncpindices.hh:48
The type of the base class for all problems which use this model.
Definition fvbaseproperties.hh:84
Type of object for specifying boundary conditions.
Definition fvbaseproperties.hh:119
Enable diffusive fluxes?
Definition multiphasebaseproperties.hh:79
Specify whether energy should be considered as a conservation quantity or not.
Definition multiphasebaseproperties.hh:76
Data required to calculate a flux over a face.
Definition fvbaseproperties.hh:149
Enumerations used by the model.
Definition multiphasebaseproperties.hh:48
The secondary variables within a sub-control volume.
Definition fvbaseproperties.hh:133
The type of the local residual function.
Definition fvbaseproperties.hh:94
The type of the model.
Definition basicproperties.hh:88
The unmodified weight for the fugacity primary variables.
Definition ncpproperties.hh:48
The unmodified weight for the pressure primary variable.
Definition ncpproperties.hh:42
The weight for the saturation primary variables.
Definition ncpproperties.hh:45
Specifies the type of the actual Newton method.
Definition nullconvergencewriter.hh:39
A vector of primary variables within a sub-control volume.
Definition fvbaseproperties.hh:130
Vector containing volumetric or areal rates of quantities.
Definition fvbaseproperties.hh:116
Define the type tag for the compositional NCP model.
Definition ncpmodel.hh:71
VTK output module for the fluid composition.
VTK output module for quantities which make sense for models which incorperate molecular diffusion.
VTK output module for quantities which make sense for models which assume thermal equilibrium.