3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
porousmediumflow/compositional/localresidual.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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program 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 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
26#ifndef DUMUX_COMPOSITIONAL_LOCAL_RESIDUAL_HH
27#define DUMUX_COMPOSITIONAL_LOCAL_RESIDUAL_HH
28
29#include <vector>
30#include <dune/common/exceptions.hh>
33
34namespace Dumux {
35
41template<class TypeTag>
42class CompositionalLocalResidual: public GetPropType<TypeTag, Properties::BaseLocalResidual>
43{
48 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
49 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
50 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
53 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
55 using Element = typename GridView::template Codim<0>::Entity;
56 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
61 using Indices = typename ModelTraits::Indices;
62
63 static constexpr int numPhases = ModelTraits::numFluidPhases();
64 static constexpr int numComponents = ModelTraits::numFluidComponents();
65 static constexpr bool useMoles = ModelTraits::useMoles();
66
67 enum { conti0EqIdx = Indices::conti0EqIdx };
68
70 static constexpr int replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
71 static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < numComponents;
72
73public:
74 using ParentType::ParentType;
75
87 NumEqVector computeStorage(const Problem& problem,
88 const SubControlVolume& scv,
89 const VolumeVariables& volVars) const
90 {
91 NumEqVector storage(0.0);
92
93 const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
94 { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
95
96 const auto massOrMoleFraction= [](const auto& volVars, const int phaseIdx, const int compIdx)
97 { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
98
99 // compute storage term of all components within all phases
100 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
101 {
102 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
103 {
104 auto eqIdx = conti0EqIdx + compIdx;
105 if (eqIdx != replaceCompEqIdx)
106 storage[eqIdx] += volVars.porosity()
107 * volVars.saturation(phaseIdx)
108 * massOrMoleDensity(volVars, phaseIdx)
109 * massOrMoleFraction(volVars, phaseIdx, compIdx);
110 }
111
112 // in case one balance is substituted by the total mole balance
113 if (useTotalMoleOrMassBalance)
114 storage[replaceCompEqIdx] += massOrMoleDensity(volVars, phaseIdx)
115 * volVars.porosity()
116 * volVars.saturation(phaseIdx);
117
119 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
120 }
121
123 EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
124
125 return storage;
126 }
127
139 NumEqVector computeFlux(const Problem& problem,
140 const Element& element,
141 const FVElementGeometry& fvGeometry,
142 const ElementVolumeVariables& elemVolVars,
143 const SubControlVolumeFace& scvf,
144 const ElementFluxVariablesCache& elemFluxVarsCache) const
145 {
146 FluxVariables fluxVars;
147 fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
148 static constexpr auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
149 // get upwind weights into local scope
150 NumEqVector flux(0.0);
151
152 const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
153 { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
154
155 const auto massOrMoleFraction= [](const auto& volVars, const int phaseIdx, const int compIdx)
156 { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
157
158 // advective fluxes
159 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
160 {
161 const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
162 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
163 {
164 // get equation index
165 const auto eqIdx = conti0EqIdx + compIdx;
166
167 // the physical quantities for which we perform upwinding
168 const auto upwindTerm = [&massOrMoleDensity, &massOrMoleFraction, phaseIdx, compIdx] (const auto& volVars)
169 { return massOrMoleDensity(volVars, phaseIdx)*massOrMoleFraction(volVars, phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
170
171 if (eqIdx != replaceCompEqIdx)
172 flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
173
174 // diffusive fluxes (only for the component balances)
175 if(eqIdx != replaceCompEqIdx)
176 {
177 //check for the reference system and adapt units of the diffusive flux accordingly.
178 if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
179 flux[eqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx)
180 : diffusiveFluxes[compIdx];
181 else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
182 flux[eqIdx] += useMoles ? diffusiveFluxes[compIdx]
183 : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
184 else
185 DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
186 }
187 }
188
189 // in case one balance is substituted by the total mole balance
190 if (useTotalMoleOrMassBalance)
191 {
192 // the physical quantities for which we perform upwinding
193 const auto upwindTerm = [&massOrMoleDensity, phaseIdx] (const auto& volVars)
194 { return massOrMoleDensity(volVars, phaseIdx)*volVars.mobility(phaseIdx); };
195
196 flux[replaceCompEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
197
198 for(int compIdx = 0; compIdx < numComponents; ++compIdx)
199 {
200 //check for the reference system and adapt units of the diffusive flux accordingly.
201 if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
202 flux[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
203 else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
204 flux[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]
205 : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
206 else
207 DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
208 }
209 }
210
212 EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
213 }
214
216 EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
217
218 return flux;
219 }
220
221protected:
222 Implementation *asImp_()
223 { return static_cast<Implementation *> (this); }
224
225 const Implementation *asImp_() const
226 { return static_cast<const Implementation *> (this); }
227};
228
229} // end namespace Dumux
230
231#endif
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
VolumeVariables::PrimaryVariables::value_type massOrMoleFraction(const VolumeVariables &volVars, ReferenceSystemFormulation referenceSys, const int phaseIdx, const int compIdx)
returns the mass or mole fraction to be used in Fick's law based on the reference system
Definition: referencesystemformulation.hh:66
Definition: adapt.hh:29
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:149
Element-wise calculation of the local residual for problems using compositional fully implicit model.
Definition: porousmediumflow/compositional/localresidual.hh:43
NumEqVector computeFlux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const ElementFluxVariablesCache &elemFluxVarsCache) const
Evaluates the total flux of all conservation quantities over a face of a sub-control volume.
Definition: porousmediumflow/compositional/localresidual.hh:139
NumEqVector computeStorage(const Problem &problem, const SubControlVolume &scv, const VolumeVariables &volVars) const
Evaluates the amount of all conservation quantities (e.g. phase mass) within a sub-control volume.
Definition: porousmediumflow/compositional/localresidual.hh:87
Implementation * asImp_()
Definition: porousmediumflow/compositional/localresidual.hh:222
const Implementation * asImp_() const
Definition: porousmediumflow/compositional/localresidual.hh:225
Declares all properties used in Dumux.