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 <dune/common/fvector.hh>
32
33#include <opm/common/Exceptions.hpp>
34
35#include <opm/material/common/Valgrind.hpp>
36#include <opm/material/densead/Math.hpp>
37
47
54
55#include <algorithm>
56#include <cassert>
57#include <cmath>
58#include <memory>
59#include <sstream>
60#include <string>
61#include <tuple>
62#include <vector>
63
64namespace Opm {
65
66template <class TypeTag>
67class NcpModel;
68
69}
70
71namespace Opm::Properties {
72
73namespace TTag {
74
78struct NcpModel { using InheritsFrom = std::tuple<MultiPhaseBaseModel>; };
79
80} // namespace TTag
81
83template<class TypeTag>
84struct LocalResidual<TypeTag, TTag::NcpModel>
86
88template<class TypeTag>
89struct NewtonMethod<TypeTag, TTag::NcpModel>
91
93template<class TypeTag>
94struct Model<TypeTag, TTag::NcpModel>
96
98template<class TypeTag>
99struct BaseProblem<TypeTag, TTag::NcpModel>
101
103template<class TypeTag>
104struct EnableEnergy<TypeTag, TTag::NcpModel>
105{ static constexpr bool value = false; };
106
108template<class TypeTag>
109struct EnableDiffusion<TypeTag, TTag::NcpModel>
110{ static constexpr bool value = false; };
111
113template<class TypeTag>
114struct RateVector<TypeTag, TTag::NcpModel>
116
118template<class TypeTag>
119struct BoundaryRateVector<TypeTag, TTag::NcpModel>
121
123template<class TypeTag>
124struct PrimaryVariables<TypeTag, TTag::NcpModel>
126
128template<class TypeTag>
129struct IntensiveQuantities<TypeTag, TTag::NcpModel>
131
133template<class TypeTag>
134struct ExtensiveQuantities<TypeTag, TTag::NcpModel>
136
138template<class TypeTag>
139struct Indices<TypeTag, TTag::NcpModel>
141
143template<class TypeTag>
144struct NcpPressureBaseWeight<TypeTag, TTag::NcpModel>
145{
147 static constexpr type value = 1.0;
148};
149
151template<class TypeTag>
152struct NcpSaturationsBaseWeight<TypeTag, TTag::NcpModel>
153{
155 static constexpr type value = 1.0;
156};
157
159template<class TypeTag>
160struct NcpFugacitiesBaseWeight<TypeTag, TTag::NcpModel>
161{
163 static constexpr type value = 1.0e-6;
164};
165
166} // namespace Opm::Properties
167
168namespace Opm {
169
244template <class TypeTag>
246 : public MultiPhaseBaseModel<TypeTag>
247{
248 using ParentType = MultiPhaseBaseModel<TypeTag>;
249
256
257 static constexpr int numPhases = FluidSystem::numPhases;
258 static constexpr int numComponents = FluidSystem::numComponents;
259 static constexpr int fugacity0Idx = Indices::fugacity0Idx;
260 enum { pressure0Idx = Indices::pressure0Idx };
261 enum { saturation0Idx = Indices::saturation0Idx };
262 static constexpr int conti0EqIdx = Indices::conti0EqIdx;
263 static constexpr int ncp0EqIdx = Indices::ncp0EqIdx;
264 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
265 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
266
267 using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
268
269 using Toolbox = MathToolbox<Evaluation>;
270
273
274public:
275 explicit NcpModel(Simulator& simulator)
276 : ParentType(simulator)
277 {}
278
282 static void registerParameters()
283 {
285
286 DiffusionModule::registerParameters();
287 EnergyModule::registerParameters();
288
289 // register runtime parameters of the VTK output modules
291
292 if constexpr (enableDiffusion) {
294 }
295
296 if constexpr (enableEnergy) {
298 }
299 }
300
305 {
306 ParentType::finishInit();
307
308 minActivityCoeff_.resize(this->numGridDof());
309 std::fill(minActivityCoeff_.begin(), minActivityCoeff_.end(), 1.0);
310 }
311
313 {
314 ParentType::adaptGrid();
315 minActivityCoeff_.resize(this->numGridDof());
316 }
317
321 static std::string name()
322 { return "ncp"; }
323
327 std::string primaryVarName(unsigned pvIdx) const
328 {
329 const std::string s = EnergyModule::primaryVarName(pvIdx);
330 if (!s.empty()) {
331 return s;
332 }
333
334 std::ostringstream oss;
335 if (pvIdx == pressure0Idx) {
336 oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
337 }
338 else if (saturation0Idx <= pvIdx && pvIdx < saturation0Idx + (numPhases - 1)) {
339 oss << "saturation_" << FluidSystem::phaseName(/*phaseIdx=*/pvIdx - saturation0Idx);
340 }
341 else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
342 oss << "fugacity^" << FluidSystem::componentName(pvIdx - fugacity0Idx);
343 }
344 else {
345 assert(false);
346 }
347
348 return oss.str();
349 }
350
354 std::string eqName(unsigned eqIdx) const
355 {
356 const std::string s = EnergyModule::eqName(eqIdx);
357 if (!s.empty()) {
358 return s;
359 }
360
361 std::ostringstream oss;
362 if (conti0EqIdx <= eqIdx && eqIdx < conti0EqIdx + numComponents) {
363 oss << "continuity^" << FluidSystem::componentName(eqIdx - conti0EqIdx);
364 }
365 else if (ncp0EqIdx <= eqIdx && eqIdx < ncp0EqIdx + numPhases) {
366 oss << "ncp_" << FluidSystem::phaseName(/*phaseIdx=*/eqIdx - ncp0EqIdx);
367 }
368 else {
369 assert(false);
370 }
371
372 return oss.str();
373 }
374
379 {
380 ParentType::updateBegin();
381
382 // find the a reference pressure. The first degree of freedom
383 // might correspond to non-interior entities which would lead
384 // to an undefined value, so we have to iterate...
385 for (unsigned dofIdx = 0; dofIdx < this->numGridDof(); ++dofIdx) {
386 if (this->isLocalDof(dofIdx)) {
388 this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
389 break;
390 }
391 }
392 }
393
397 void updatePVWeights(const ElementContext& elemCtx) const
398 {
399 for (unsigned dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) {
400 const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
401
402 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
403 minActivityCoeff_[globalIdx][compIdx] = 1e100;
404 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
405 const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();
406
407 minActivityCoeff_[globalIdx][compIdx] =
408 std::min(minActivityCoeff_[globalIdx][compIdx],
409 Toolbox::value(fs.fugacityCoefficient(phaseIdx, compIdx)) *
410 Toolbox::value(fs.pressure(phaseIdx)));
411 Valgrind::CheckDefined(minActivityCoeff_[globalIdx][compIdx]);
412 }
413 if (minActivityCoeff_[globalIdx][compIdx] <= 0) {
414 throw NumericalProblem("The minimum activity coefficient for component " +
415 std::to_string(compIdx) + " on DOF " +
416 std::to_string(globalIdx) + " is negative or zero!");
417 }
418 }
419 }
420 }
421
425 Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
426 {
427 const Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
428 Scalar result;
429 if (tmp > 0) {
430 // energy related quantity
431 result = tmp;
432 }
433 else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
434 // component fugacity
435 const unsigned compIdx = pvIdx - fugacity0Idx;
436 assert(compIdx <= numComponents);
437
438 Valgrind::CheckDefined(minActivityCoeff_[globalDofIdx][compIdx]);
439 constexpr Scalar fugacityBaseWeight =
440 getPropValue<TypeTag, Properties::NcpFugacitiesBaseWeight>();
441 result = fugacityBaseWeight / minActivityCoeff_[globalDofIdx][compIdx];
442 }
443 else if (Indices::pressure0Idx == pvIdx) {
444 constexpr Scalar pressureBaseWeight =
445 getPropValue<TypeTag, Properties::NcpPressureBaseWeight>();
446 result = pressureBaseWeight / referencePressure_;
447 }
448 else {
449#ifndef NDEBUG
450 const unsigned phaseIdx = pvIdx - saturation0Idx;
451 assert(phaseIdx < numPhases - 1);
452#endif
453
454 // saturation
455 constexpr Scalar saturationsBaseWeight =
456 getPropValue<TypeTag, Properties::NcpSaturationsBaseWeight>();
457 result = saturationsBaseWeight;
458 }
459
460 assert(std::isfinite(result));
461 assert(result > 0);
462
463 return result;
464 }
465
469 Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
470 {
471 const Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
472 if (tmp > 0) {
473 // an energy related equation
474 return tmp;
475 }
476 // an NCP
477 else if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases) {
478 return 1.0;
479 }
480
481 // a mass conservation equation
482 const unsigned compIdx = eqIdx - Indices::conti0EqIdx;
483 assert(compIdx <= numComponents);
484
485 // make all kg equal
486 return FluidSystem::molarMass(compIdx);
487 }
488
496 Scalar minActivityCoeff(unsigned globalDofIdx, unsigned compIdx) const
497 { return minActivityCoeff_[globalDofIdx][compIdx]; }
498
503 {
505
506 this->addOutputModule(std::make_unique<VtkCompositionModule<TypeTag>>(this->simulator_));
507 if constexpr (enableDiffusion) {
508 this->addOutputModule(std::make_unique<VtkDiffusionModule<TypeTag>>(this->simulator_));
509 }
510 if constexpr (enableEnergy) {
511 this->addOutputModule(std::make_unique<VtkEnergyModule<TypeTag>>(this->simulator_));
512 }
513 }
514
515 mutable Scalar referencePressure_;
516 mutable std::vector<ComponentVector> minActivityCoeff_;
517};
518
519} // namespace Opm
520
521#endif
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:51
Provides the auxiliary methods required for consideration of the energy equation.
Definition: energymodule.hh:54
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
Definition: multiphasebasemodel.hh:168
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition: multiphasebasemodel.hh:194
void registerOutputModules_()
Definition: multiphasebasemodel.hh:270
The base class for the problems of ECFV discretizations which deal with a multi-phase flow through a ...
Definition: multiphasebaseproblem.hh:65
Implements a boundary vector for the fully implicit compositional multi-phase NCP model.
Definition: ncpboundaryratevector.hh:50
This template class represents the extensive quantities of the compositional NCP model.
Definition: ncpextensivequantities.hh:49
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition: ncpintensivequantities.hh:61
Details needed to calculate the local residual in the compositional multi-phase NCP-model .
Definition: ncplocalresidual.hh:53
A compositional multi-phase model based on non-linear complementarity functions.
Definition: ncpmodel.hh:247
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:397
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition: ncpmodel.hh:354
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition: ncpmodel.hh:327
NcpModel(Simulator &simulator)
Definition: ncpmodel.hh:275
static std::string name()
Definition: ncpmodel.hh:321
Scalar referencePressure_
Definition: ncpmodel.hh:515
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition: ncpmodel.hh:425
void adaptGrid()
Definition: ncpmodel.hh:312
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:496
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition: ncpmodel.hh:282
void finishInit()
Apply the initial conditions to the model.
Definition: ncpmodel.hh:304
void updateBegin()
Called by the update() method before it tries to apply the newton method. This is primary a hook whic...
Definition: ncpmodel.hh:378
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
Returns the relative weight of an equation.
Definition: ncpmodel.hh:469
void registerOutputModules_()
Definition: ncpmodel.hh:502
std::vector< ComponentVector > minActivityCoeff_
Definition: ncpmodel.hh:516
A Newton solver specific to the NCP model.
Definition: ncpnewtonmethod.hh:56
Represents the primary variables used by the compositional multi-phase NCP model.
Definition: ncpprimaryvariables.hh:56
Implements a vector representing mass, molar or volumetric rates.
Definition: ncpratevector.hh:54
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:87
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:88
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:87
Classes required for molecular diffusion.
Contains the classes required to consider energy as a conservation quantity in a multi-phase module.
Definition: blackoilmodel.hh:79
Definition: blackoilboundaryratevector.hh:39
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:233
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Declares the properties required for the NCP compositional multi-phase model.
The primary variable and equation indices for the compositional multi-phase NCP model.
Definition: ncpindices.hh:47
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:91
Specify whether energy should be considered as a conservation quantity or not.
Definition: multiphasebaseproperties.hh:87
Data required to calculate a flux over a face.
Definition: fvbaseproperties.hh:149
Enumerations used by the model.
Definition: multiphasebaseproperties.hh:51
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
GetPropType< TypeTag, Scalar > type
Definition: ncpmodel.hh:162
The unmodified weight for the fugacity primary variables.
Definition: ncpproperties.hh:47
GetPropType< TypeTag, Scalar > type
Definition: ncpmodel.hh:146
The unmodified weight for the pressure primary variable.
Definition: ncpproperties.hh:39
GetPropType< TypeTag, Scalar > type
Definition: ncpmodel.hh:154
The weight for the saturation primary variables.
Definition: ncpproperties.hh:43
Specifies the type of the actual Newton method.
Definition: newtonmethodproperties.hh:32
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:78
std::tuple< MultiPhaseBaseModel > InheritsFrom
Definition: ncpmodel.hh:78