ncplocalresidual.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_LOCAL_RESIDUAL_HH
29#define EWOMS_NCP_LOCAL_RESIDUAL_HH
30
31#include "ncpproperties.hh"
32
35
36#include <opm/material/common/Valgrind.hpp>
37
38namespace Opm {
45template <class TypeTag>
46class NcpLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalResidual>
47{
56
57 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
58 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
59 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
60 enum { ncp0EqIdx = Indices::ncp0EqIdx };
61 enum { conti0EqIdx = Indices::conti0EqIdx };
62
63 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
65
66 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
68
69 using EvalEqVector = Dune::FieldVector<Evaluation, numEq>;
70 using ElemEvalEqVector = Dune::BlockVector<EvalEqVector>;
71 using Toolbox = Opm::MathToolbox<Evaluation>;
72
73public:
77 template <class LhsEval>
78 void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
79 const ElementContext& elemCtx,
80 unsigned dofIdx,
81 unsigned timeIdx,
82 unsigned phaseIdx) const
83 {
84 const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
85 const auto& fluidState = intQuants.fluidState();
86
87 // compute storage term of all components within all phases
88 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
89 unsigned eqIdx = conti0EqIdx + compIdx;
90 storage[eqIdx] +=
91 Toolbox::template decay<LhsEval>(fluidState.molarity(phaseIdx, compIdx))
92 * Toolbox::template decay<LhsEval>(fluidState.saturation(phaseIdx))
93 * Toolbox::template decay<LhsEval>(intQuants.porosity());
94 }
95
96 EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
97 }
98
102 template <class LhsEval>
103 void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
104 const ElementContext& elemCtx,
105 unsigned dofIdx,
106 unsigned timeIdx) const
107 {
108 storage = 0;
109 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
110 addPhaseStorage(storage, elemCtx, dofIdx, timeIdx, phaseIdx);
111
112 EnergyModule::addSolidEnergyStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx));
113 }
114
118 void computeFlux(RateVector& flux,
119 const ElementContext& elemCtx,
120 unsigned scvfIdx,
121 unsigned timeIdx) const
122 {
123 flux = 0.0;
124 addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
125 Opm::Valgrind::CheckDefined(flux);
126
127 addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
128 Opm::Valgrind::CheckDefined(flux);
129 }
130
134 void addAdvectiveFlux(RateVector& flux,
135 const ElementContext& elemCtx,
136 unsigned scvfIdx,
137 unsigned timeIdx) const
138 {
139 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
140
141 unsigned focusDofIdx = elemCtx.focusDofIndex();
142 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
143 // data attached to upstream and the downstream DOFs
144 // of the current phase
145 unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
146 const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
147
148 // this is a bit hacky because it is specific to the element-centered
149 // finite volume scheme. (N.B. that if finite differences are used to
150 // linearize the system of equations, it does not matter.)
151 if (upIdx == focusDofIdx) {
152 Evaluation tmp =
153 up.fluidState().molarDensity(phaseIdx)
154 * extQuants.volumeFlux(phaseIdx);
155
156 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
157 flux[conti0EqIdx + compIdx] +=
158 tmp*up.fluidState().moleFraction(phaseIdx, compIdx);
159 }
160 }
161 else {
162 Evaluation tmp =
163 Toolbox::value(up.fluidState().molarDensity(phaseIdx))
164 * extQuants.volumeFlux(phaseIdx);
165
166 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
167 flux[conti0EqIdx + compIdx] +=
168 tmp*Toolbox::value(up.fluidState().moleFraction(phaseIdx, compIdx));
169 }
170 }
171 }
172
173 EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
174 }
175
179 void addDiffusiveFlux(RateVector& flux,
180 const ElementContext& elemCtx,
181 unsigned scvfIdx,
182 unsigned timeIdx) const
183 {
184 DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
185 EnergyModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
186 }
187
194 void computeSource(RateVector& source,
195 const ElementContext& elemCtx,
196 unsigned dofIdx,
197 unsigned timeIdx) const
198 {
199 Opm::Valgrind::SetUndefined(source);
200 elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
201 Opm::Valgrind::CheckDefined(source);
202
203 // evaluate the NCPs (i.e., the "phase presence" equations)
204 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
205 source[ncp0EqIdx + phaseIdx] =
206 phaseNcp(elemCtx, dofIdx, timeIdx, phaseIdx);
207 }
208 }
209
213 template <class LhsEval = Evaluation>
214 LhsEval phaseNcp(const ElementContext& elemCtx,
215 unsigned dofIdx,
216 unsigned timeIdx,
217 unsigned phaseIdx) const
218 {
219 const auto& fluidState = elemCtx.intensiveQuantities(dofIdx, timeIdx).fluidState();
220 using FluidState = typename std::remove_const<typename std::remove_reference<decltype(fluidState)>::type>::type;
221
222 using LhsToolbox = Opm::MathToolbox<LhsEval>;
223
224 const LhsEval& a = phaseNotPresentIneq_<FluidState, LhsEval>(fluidState, phaseIdx);
225 const LhsEval& b = phasePresentIneq_<FluidState, LhsEval>(fluidState, phaseIdx);
226 return LhsToolbox::min(a, b);
227 }
228
229private:
234 template <class FluidState, class LhsEval>
235 LhsEval phasePresentIneq_(const FluidState& fluidState, unsigned phaseIdx) const
236 {
237 using FsToolbox = Opm::MathToolbox<typename FluidState::Scalar>;
238
239 return FsToolbox::template decay<LhsEval>(fluidState.saturation(phaseIdx));
240 }
241
246 template <class FluidState, class LhsEval>
247 LhsEval phaseNotPresentIneq_(const FluidState& fluidState, unsigned phaseIdx) const
248 {
249 using FsToolbox = Opm::MathToolbox<typename FluidState::Scalar>;
250
251 // difference of sum of mole fractions in the phase from 100%
252 LhsEval a = 1.0;
253 for (unsigned i = 0; i < numComponents; ++i)
254 a -= FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, i));
255 return a;
256 }
257};
258
259} // namespace Opm
260
261#endif
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:48
Provides the auxiliary methods required for consideration of the energy equation.
Definition: energymodule.hh:50
Details needed to calculate the local residual in the compositional multi-phase NCP-model .
Definition: ncplocalresidual.hh:47
LhsEval phaseNcp(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned phaseIdx) const
Returns the value of the NCP-function for a phase.
Definition: ncplocalresidual.hh:214
void addAdvectiveFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx) const
Add the advective mass flux at a given flux integration point.
Definition: ncplocalresidual.hh:134
void addDiffusiveFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx) const
Adds the diffusive flux at a given flux integration point.
Definition: ncplocalresidual.hh:179
void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume.
Definition: ncplocalresidual.hh:118
void addPhaseStorage(Dune::FieldVector< LhsEval, numEq > &storage, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned phaseIdx) const
Adds the amount all conservation quantities (e.g. phase mass) within a single fluid phase.
Definition: ncplocalresidual.hh:78
void computeStorage(Dune::FieldVector< LhsEval, numEq > &storage, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx) const
Evaluate the amount all conservation quantities (e.g. phase mass) within a finite sub-control volume.
Definition: ncplocalresidual.hh:103
void computeSource(RateVector &source, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx) const
Calculate the source term of the equation.
Definition: ncplocalresidual.hh:194
Classes required for molecular diffusion.
Contains the classes required to consider energy as a conservation quantity in a multi-phase module.
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
Declares the properties required for the NCP compositional multi-phase model.