Loading...
Searching...
No Matches
ncpmodel.hh
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef EWOMS_NCP_MODEL_HH
#define EWOMS_NCP_MODEL_HH
#include <dune/common/fvector.hh>
#include <opm/common/Exceptions.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/models/ncp/ncpproperties.hh>
#include <opm/models/ncp/ncplocalresidual.hh>
#include <opm/models/ncp/ncpextensivequantities.hh>
#include <opm/models/ncp/ncpprimaryvariables.hh>
#include <opm/models/ncp/ncpboundaryratevector.hh>
#include <opm/models/ncp/ncpratevector.hh>
#include <opm/models/ncp/ncpintensivequantities.hh>
#include <opm/models/ncp/ncpnewtonmethod.hh>
#include <opm/models/ncp/ncpindices.hh>
#include <opm/models/common/diffusionmodule.hh>
#include <opm/models/common/energymodule.hh>
#include <opm/models/common/multiphasebasemodel.hh>
#include <opm/models/io/vtkcompositionmodule.hpp>
#include <opm/models/io/vtkdiffusionmodule.hpp>
#include <opm/models/io/vtkenergymodule.hpp>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <memory>
#include <sstream>
#include <string>
#include <tuple>
#include <vector>
namespace Opm {
template <class TypeTag>
class NcpModel;
}
namespace Opm::Properties {
namespace TTag {
struct NcpModel { using InheritsFrom = std::tuple<MultiPhaseBaseModel>; };
} // namespace TTag
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::NcpModel>
{ using type = NcpLocalResidual<TypeTag>; };
template<class TypeTag>
struct NewtonMethod<TypeTag, TTag::NcpModel>
{ using type = NcpNewtonMethod<TypeTag>; };
template<class TypeTag>
struct Model<TypeTag, TTag::NcpModel>
{ using type = NcpModel<TypeTag>; };
template<class TypeTag>
struct BaseProblem<TypeTag, TTag::NcpModel>
{ using type = MultiPhaseBaseProblem<TypeTag>; };
template<class TypeTag>
struct EnableEnergy<TypeTag, TTag::NcpModel>
{ static constexpr bool value = false; };
template<class TypeTag>
struct EnableDiffusion<TypeTag, TTag::NcpModel>
{ static constexpr bool value = false; };
template<class TypeTag>
struct RateVector<TypeTag, TTag::NcpModel>
{ using type = NcpRateVector<TypeTag>; };
template<class TypeTag>
struct BoundaryRateVector<TypeTag, TTag::NcpModel>
{ using type = NcpBoundaryRateVector<TypeTag>; };
template<class TypeTag>
struct PrimaryVariables<TypeTag, TTag::NcpModel>
{ using type = NcpPrimaryVariables<TypeTag>; };
template<class TypeTag>
struct IntensiveQuantities<TypeTag, TTag::NcpModel>
{ using type = NcpIntensiveQuantities<TypeTag>; };
template<class TypeTag>
struct ExtensiveQuantities<TypeTag, TTag::NcpModel>
{ using type = NcpExtensiveQuantities<TypeTag>; };
template<class TypeTag>
struct Indices<TypeTag, TTag::NcpModel>
{ using type = NcpIndices<TypeTag, 0>; };
template<class TypeTag>
struct NcpPressureBaseWeight<TypeTag, TTag::NcpModel>
{
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1.0;
};
template<class TypeTag>
struct NcpSaturationsBaseWeight<TypeTag, TTag::NcpModel>
{
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1.0;
};
template<class TypeTag>
struct NcpFugacitiesBaseWeight<TypeTag, TTag::NcpModel>
{
using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1.0e-6;
};
} // namespace Opm::Properties
namespace Opm {
template <class TypeTag>
class NcpModel
: public MultiPhaseBaseModel<TypeTag>
{
using ParentType = MultiPhaseBaseModel<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
static constexpr int numPhases = FluidSystem::numPhases;
static constexpr int numComponents = FluidSystem::numComponents;
static constexpr int fugacity0Idx = Indices::fugacity0Idx;
enum { pressure0Idx = Indices::pressure0Idx };
enum { saturation0Idx = Indices::saturation0Idx };
static constexpr int conti0EqIdx = Indices::conti0EqIdx;
static constexpr int ncp0EqIdx = Indices::ncp0EqIdx;
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
using Toolbox = MathToolbox<Evaluation>;
using DiffusionModule = ::Opm::DiffusionModule<TypeTag, enableDiffusion>;
using EnergyModule = ::Opm::EnergyModule<TypeTag, enableEnergy>;
public:
explicit NcpModel(Simulator& simulator)
: ParentType(simulator)
{}
static void registerParameters()
{
ParentType::registerParameters();
DiffusionModule::registerParameters();
EnergyModule::registerParameters();
// register runtime parameters of the VTK output modules
VtkCompositionModule<TypeTag>::registerParameters();
if constexpr (enableDiffusion) {
VtkDiffusionModule<TypeTag>::registerParameters();
}
if constexpr (enableEnergy) {
VtkEnergyModule<TypeTag>::registerParameters();
}
}
void finishInit()
{
ParentType::finishInit();
minActivityCoeff_.resize(this->numGridDof());
std::fill(minActivityCoeff_.begin(), minActivityCoeff_.end(), 1.0);
}
void adaptGrid()
{
ParentType::adaptGrid();
minActivityCoeff_.resize(this->numGridDof());
}
static std::string name()
{ return "ncp"; }
std::string primaryVarName(unsigned pvIdx) const
{
const std::string s = EnergyModule::primaryVarName(pvIdx);
if (!s.empty()) {
return s;
}
std::ostringstream oss;
if (pvIdx == pressure0Idx) {
oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
}
else if (saturation0Idx <= pvIdx && pvIdx < saturation0Idx + (numPhases - 1)) {
oss << "saturation_" << FluidSystem::phaseName(/*phaseIdx=*/pvIdx - saturation0Idx);
}
else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
oss << "fugacity^" << FluidSystem::componentName(pvIdx - fugacity0Idx);
}
else {
assert(false);
}
return oss.str();
}
std::string eqName(unsigned eqIdx) const
{
const std::string s = EnergyModule::eqName(eqIdx);
if (!s.empty()) {
return s;
}
std::ostringstream oss;
if (conti0EqIdx <= eqIdx && eqIdx < conti0EqIdx + numComponents) {
oss << "continuity^" << FluidSystem::componentName(eqIdx - conti0EqIdx);
}
else if (ncp0EqIdx <= eqIdx && eqIdx < ncp0EqIdx + numPhases) {
oss << "ncp_" << FluidSystem::phaseName(/*phaseIdx=*/eqIdx - ncp0EqIdx);
}
else {
assert(false);
}
return oss.str();
}
void updateBegin()
{
ParentType::updateBegin();
// find the a reference pressure. The first degree of freedom
// might correspond to non-interior entities which would lead
// to an undefined value, so we have to iterate...
for (unsigned dofIdx = 0; dofIdx < this->numGridDof(); ++dofIdx) {
if (this->isLocalDof(dofIdx)) {
referencePressure_ =
this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
break;
}
}
}
void updatePVWeights(const ElementContext& elemCtx) const
{
for (unsigned dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) {
const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
minActivityCoeff_[globalIdx][compIdx] = 1e100;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();
minActivityCoeff_[globalIdx][compIdx] =
std::min(minActivityCoeff_[globalIdx][compIdx],
Toolbox::value(fs.fugacityCoefficient(phaseIdx, compIdx)) *
Toolbox::value(fs.pressure(phaseIdx)));
Valgrind::CheckDefined(minActivityCoeff_[globalIdx][compIdx]);
}
if (minActivityCoeff_[globalIdx][compIdx] <= 0) {
throw NumericalProblem("The minimum activity coefficient for component " +
std::to_string(compIdx) + " on DOF " +
std::to_string(globalIdx) + " is negative or zero!");
}
}
}
}
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
{
const Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
Scalar result;
if (tmp > 0) {
// energy related quantity
result = tmp;
}
else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
// component fugacity
const unsigned compIdx = pvIdx - fugacity0Idx;
assert(compIdx <= numComponents);
Valgrind::CheckDefined(minActivityCoeff_[globalDofIdx][compIdx]);
constexpr Scalar fugacityBaseWeight =
getPropValue<TypeTag, Properties::NcpFugacitiesBaseWeight>();
result = fugacityBaseWeight / minActivityCoeff_[globalDofIdx][compIdx];
}
else if (Indices::pressure0Idx == pvIdx) {
constexpr Scalar pressureBaseWeight =
getPropValue<TypeTag, Properties::NcpPressureBaseWeight>();
result = pressureBaseWeight / referencePressure_;
}
else {
#ifndef NDEBUG
const unsigned phaseIdx = pvIdx - saturation0Idx;
assert(phaseIdx < numPhases - 1);
#endif
// saturation
constexpr Scalar saturationsBaseWeight =
getPropValue<TypeTag, Properties::NcpSaturationsBaseWeight>();
result = saturationsBaseWeight;
}
assert(std::isfinite(result));
assert(result > 0);
return result;
}
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
{
const Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
if (tmp > 0) {
// an energy related equation
return tmp;
}
// an NCP
else if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases) {
return 1.0;
}
// a mass conservation equation
const unsigned compIdx = eqIdx - Indices::conti0EqIdx;
assert(compIdx <= numComponents);
// make all kg equal
return FluidSystem::molarMass(compIdx);
}
Scalar minActivityCoeff(unsigned globalDofIdx, unsigned compIdx) const
{ return minActivityCoeff_[globalDofIdx][compIdx]; }
void registerOutputModules_()
{
ParentType::registerOutputModules_();
this->addOutputModule(std::make_unique<VtkCompositionModule<TypeTag>>(this->simulator_));
if constexpr (enableDiffusion) {
this->addOutputModule(std::make_unique<VtkDiffusionModule<TypeTag>>(this->simulator_));
}
if constexpr (enableEnergy) {
this->addOutputModule(std::make_unique<VtkEnergyModule<TypeTag>>(this->simulator_));
}
}
mutable Scalar referencePressure_;
mutable std::vector<ComponentVector> minActivityCoeff_;
};
} // namespace Opm
#endif
Opm::DiffusionModule
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition diffusionmodule.hh:51
Opm::EnergyModule
Provides the auxiliary methods required for consideration of the energy equation.
Definition energymodule.hh:54
Opm::MultiPhaseBaseModel::registerParameters
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition multiphasebasemodel.hh:190
Opm::MultiPhaseBaseProblem
The base class for the problems of ECFV discretizations which deal with a multi-phase flow through a ...
Definition multiphasebaseproblem.hh:65
Opm::NcpBoundaryRateVector
Implements a boundary vector for the fully implicit compositional multi-phase NCP model.
Definition ncpboundaryratevector.hh:50
Opm::NcpExtensiveQuantities
This template class represents the extensive quantities of the compositional NCP model.
Definition ncpextensivequantities.hh:49
Opm::NcpIntensiveQuantities
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition ncpintensivequantities.hh:61
Opm::NcpLocalResidual
Details needed to calculate the local residual in the compositional multi-phase NCP-model .
Definition ncplocalresidual.hh:53
Opm::NcpModel
A compositional multi-phase model based on non-linear complementarity functions.
Definition ncpmodel.hh:247
Opm::NcpModel::updatePVWeights
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
Opm::NcpModel::eqName
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition ncpmodel.hh:354
Opm::NcpModel::primaryVarName
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition ncpmodel.hh:327
Opm::NcpModel::name
static std::string name()
Definition ncpmodel.hh:321
Opm::NcpModel::primaryVarWeight
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition ncpmodel.hh:425
Opm::NcpModel::minActivityCoeff
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:499
Opm::NcpModel::registerParameters
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition ncpmodel.hh:282
Opm::NcpModel::finishInit
void finishInit()
Apply the initial conditions to the model.
Definition ncpmodel.hh:304
Opm::NcpModel::updateBegin
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition ncpmodel.hh:378
Opm::NcpModel::eqWeight
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
Returns the relative weight of an equation.
Definition ncpmodel.hh:472
Opm::NcpNewtonMethod
A Newton solver specific to the NCP model.
Definition ncpnewtonmethod.hh:56
Opm::NcpPrimaryVariables
Represents the primary variables used by the compositional multi-phase NCP model.
Definition ncpprimaryvariables.hh:56
Opm::NcpRateVector
Implements a vector representing mass, molar or volumetric rates.
Definition ncpratevector.hh:54
Opm::VtkCompositionModule
VTK output module for the fluid composition.
Definition vtkcompositionmodule.hpp:57
Opm::VtkCompositionModule::registerParameters
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkcompositionmodule.hpp:87
Opm::VtkDiffusionModule
VTK output module for quantities which make sense for models which incorperate molecular diffusion.
Definition vtkdiffusionmodule.hpp:58
Opm::VtkDiffusionModule::registerParameters
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkdiffusionmodule.hpp:88
Opm::VtkEnergyModule::registerParameters
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkenergymodule.hpp:87
diffusionmodule.hh
Classes required for molecular diffusion.
energymodule.hh
Contains the classes required to consider energy as a conservation quantity in a multi-phase module.
multiphasebasemodel.hh
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
Opm::Properties::TTag
The generic type tag for problems using the immiscible multi-phase model.
Definition blackoilmodel.hh:81
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
Opm::GetPropType
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
Opm::getPropValue
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:240
ncpboundaryratevector.hh
Implements a boundary vector for the fully implicit compositional multi-phase NCP model.
ncpextensivequantities.hh
This template class represents the extensive quantities of the compositional NCP model.
ncpindices.hh
The primary variable and equation indices for the compositional multi-phase NCP model.
ncpintensivequantities.hh
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
ncplocalresidual.hh
Details needed to calculate the local residual in the compositional multi-phase NCP-model .
ncpnewtonmethod.hh
A Newton solver specific to the NCP model.
ncpprimaryvariables.hh
Represents the primary variables used by the compositional multi-phase NCP model.
ncpproperties.hh
Declares the properties required for the NCP compositional multi-phase model.
ncpratevector.hh
Implements a vector representing mass, molar or volumetric rates.
Opm::NcpIndices
The primary variable and equation indices for the compositional multi-phase NCP model.
Definition ncpindices.hh:47
Opm::Properties::BaseProblem
The type of the base class for all problems which use this model.
Definition fvbaseproperties.hh:84
Opm::Properties::BoundaryRateVector
Type of object for specifying boundary conditions.
Definition fvbaseproperties.hh:119
Opm::Properties::EnableDiffusion
Enable diffusive fluxes?
Definition multiphasebaseproperties.hh:91
Opm::Properties::EnableEnergy
Specify whether energy should be considered as a conservation quantity or not.
Definition multiphasebaseproperties.hh:87
Opm::Properties::ExtensiveQuantities
Data required to calculate a flux over a face.
Definition fvbaseproperties.hh:149
Opm::Properties::Indices
Enumerations used by the model.
Definition multiphasebaseproperties.hh:51
Opm::Properties::IntensiveQuantities
The secondary variables within a sub-control volume.
Definition fvbaseproperties.hh:133
Opm::Properties::LocalResidual
The type of the local residual function.
Definition fvbaseproperties.hh:94
Opm::Properties::Model
The type of the model.
Definition basicproperties.hh:88
Opm::Properties::NcpFugacitiesBaseWeight
The unmodified weight for the fugacity primary variables.
Definition ncpproperties.hh:47
Opm::Properties::NcpPressureBaseWeight
The unmodified weight for the pressure primary variable.
Definition ncpproperties.hh:39
Opm::Properties::NcpSaturationsBaseWeight
The weight for the saturation primary variables.
Definition ncpproperties.hh:43
Opm::Properties::NewtonMethod
Specifies the type of the actual Newton method.
Definition newtonmethodproperties.hh:32
Opm::Properties::PrimaryVariables
A vector of primary variables within a sub-control volume.
Definition fvbaseproperties.hh:130
Opm::Properties::RateVector
Vector containing volumetric or areal rates of quantities.
Definition fvbaseproperties.hh:116
Opm::Properties::TTag::NcpModel
Define the type tag for the compositional NCP model.
Definition ncpmodel.hh:78
vtkcompositionmodule.hpp
VTK output module for the fluid composition.
vtkdiffusionmodule.hpp
VTK output module for quantities which make sense for models which incorperate molecular diffusion.
vtkenergymodule.hpp
VTK output module for quantities which make sense for models which assume thermal equilibrium.