3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
multidomain/facet/box/fickslaw.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 *****************************************************************************/
24#ifndef DUMUX_DISCRETIZATION_BOX_FACET_COUPLING_FICKS_LAW_HH
25#define DUMUX_DISCRETIZATION_BOX_FACET_COUPLING_FICKS_LAW_HH
26
27#include <vector>
28#include <cmath>
29
30#include <dune/common/exceptions.hh>
31#include <dune/common/fvector.hh>
32#include <dune/common/float_cmp.hh>
33
36
41
42namespace Dumux {
43
50template<class TypeTag, ReferenceSystemFormulation referenceSystem = ReferenceSystemFormulation::massAveraged>
52: public FicksLawImplementation<TypeTag, DiscretizationMethod::box, referenceSystem>
53{
55
57 using FVElementGeometry = typename GridGeometry::LocalView;
58 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
59 using Extrusion = Extrusion_t<GridGeometry>;
60 using GridView = typename GridGeometry::GridView;
61 using Element = typename GridView::template Codim<0>::Entity;
62
63 static constexpr int dim = GridView::dimension;
64 static constexpr int dimWorld = GridView::dimensionworld;
65
69 static constexpr int numComponents = ModelTraits::numFluidComponents();
70
72 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
73
74public:
75
76 template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
77 static ComponentFluxVector flux(const Problem& problem,
78 const Element& element,
79 const FVElementGeometry& fvGeometry,
80 const ElementVolumeVariables& elemVolVars,
81 const SubControlVolumeFace& scvf,
82 const int phaseIdx,
83 const ElementFluxVarsCache& elemFluxVarCache)
84 {
85 // if this scvf is not on an interior boundary, use the standard law
86 if (!scvf.interiorBoundary())
87 return ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarCache);
88
89 static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0);
90 if ( !Dune::FloatCmp::eq(xi, 1.0, 1e-6) )
91 DUNE_THROW(Dune::NotImplemented, "Xi != 1.0 cannot be used with the Box-Facet-Coupling scheme");
92
93 // get some references for convenience
94 const auto& fluxVarCache = elemFluxVarCache[scvf];
95 const auto& shapeValues = fluxVarCache.shapeValues();
96 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
97 const auto& insideVolVars = elemVolVars[insideScv];
98
99 // interpolate density to scvf integration point
100 Scalar rho = 0.0;
101 for (const auto& scv : scvs(fvGeometry))
102 rho += massOrMolarDensity(elemVolVars[scv], referenceSystem, phaseIdx)*shapeValues[scv.indexInElement()][0];
103
104 // on interior Neumann boundaries, evaluate the flux using the facet effective diffusion coefficient
105 ComponentFluxVector componentFlux(0.0);
106 const auto bcTypes = problem.interiorBoundaryTypes(element, scvf);
107 if (bcTypes.hasOnlyNeumann())
108 {
109 // compute tpfa flux from integration point to facet centerline
110 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
111
112 using std::sqrt;
113 // If this is a surface grid, use the square root of the facet extrusion factor
114 // as an approximate average distance from scvf ip to facet center
115 using std::sqrt;
116 const auto a = facetVolVars.extrusionFactor();
117 auto preGradX = scvf.unitOuterNormal();
118 preGradX *= dim == dimWorld ? 0.5*a : 0.5*sqrt(a);
119 preGradX /= preGradX.two_norm2();
120
121 for (int compIdx = 0; compIdx < numComponents; compIdx++)
122 {
123 if (compIdx == FluidSystem::getMainComponent(phaseIdx))
124 continue;
125
126 // interpolate mole fraction to scvf integration point
127 Scalar x = 0.0;
128 for (const auto& scv : scvs(fvGeometry))
129 x += massOrMoleFraction(elemVolVars[scv], referenceSystem, phaseIdx, compIdx)*shapeValues[scv.indexInElement()][0];
130
131 // compute the diffusive flux by means of a finite difference
132 auto gradX = preGradX;
133 gradX *= (massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx) - x);
134
135 componentFlux[compIdx] = -1.0*rho*Extrusion::area(scvf)
136 *insideVolVars.extrusionFactor()
137 *vtmv(scvf.unitOuterNormal(),
138 facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
139 gradX);
140
141 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
142 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
143 }
144
145 return componentFlux;
146 }
147
148 // on interior Dirichlet boundaries use the facet mass/mole fraction and evaluate flux
149 else if (bcTypes.hasOnlyDirichlet())
150 {
151 for (int compIdx = 0; compIdx < numComponents; compIdx++)
152 {
153 if (compIdx == FluidSystem::getMainComponent(phaseIdx))
154 continue;
155
156 // create vector with nodal mole/mass fractions
157 std::vector<Scalar> xFractions(element.subEntities(dim));
158 for (const auto& scv : scvs(fvGeometry))
159 xFractions[scv.localDofIndex()] = massOrMoleFraction(elemVolVars[scv], referenceSystem, phaseIdx, compIdx);
160
161 // substitute with facet mole/mass fractions for those scvs touching this facet
162 for (const auto& scvfJ : scvfs(fvGeometry))
163 if (scvfJ.interiorBoundary() && scvfJ.facetIndexInElement() == scvf.facetIndexInElement())
164 xFractions[ fvGeometry.scv(scvfJ.insideScvIdx()).localDofIndex() ]
165 = massOrMoleFraction(problem.couplingManager().getLowDimVolVars(element, scvfJ), referenceSystem, phaseIdx, compIdx);
166
167 // evaluate gradX at integration point
168 Dune::FieldVector<Scalar, dimWorld> gradX(0.0);
169 for (const auto& scv : scvs(fvGeometry))
170 gradX.axpy(xFractions[scv.localDofIndex()], fluxVarCache.gradN(scv.indexInElement()));
171
172 // apply matrix diffusion coefficient and return the flux
173 componentFlux[compIdx] = -1.0*rho*Extrusion::area(scvf)
174 *insideVolVars.extrusionFactor()
175 *vtmv(scvf.unitOuterNormal(),
176 insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
177 gradX);
178
179 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
180 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
181 }
182
183 return componentFlux;
184 }
185
186 // mixed boundary types are not supported
187 else
188 DUNE_THROW(Dune::NotImplemented, "Mixed boundary types are not supported");
189 }
190
191 // compute transmissibilities ti for analytical jacobians
192 template<class Problem, class ElementVolumeVariables, class FluxVarCache>
193 static std::vector<Scalar> calculateTransmissibilities(const Problem& problem,
194 const Element& element,
195 const FVElementGeometry& fvGeometry,
196 const ElementVolumeVariables& elemVolVars,
197 const SubControlVolumeFace& scvf,
198 const FluxVarCache& fluxVarCache,
199 unsigned int phaseIdx)
200 {
201 DUNE_THROW(Dune::NotImplemented, "transmissibilty computation for BoxFacetCouplingFicksLaw");
202 }
203};
204
205} // end namespace Dumux
206
207#endif
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
Helper classes to compute the integration elements.
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
VolumeVariables::PrimaryVariables::value_type massOrMolarDensity(const VolumeVariables &volVars, ReferenceSystemFormulation referenceSys, const int phaseIdx)
evaluates the density to be used in Fick's law based on the reference system
Definition: referencesystemformulation.hh:55
Dune::DenseMatrix< MAT >::value_type vtmv(const Dune::DenseVector< V1 > &v1, const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V2 > &v2)
Evaluates the scalar product of a vector v2, projected by a matrix M, with a vector v1.
Definition: math.hh:849
Definition: adapt.hh:29
typename Extrusion< T >::type Extrusion_t
Convenience alias for obtaining the extrusion type.
Definition: extrusion.hh:177
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
forward declaration of the method-specific implemetation
Definition: flux/box/fickslaw.hh:43
Specialization of Fick's Law for the box method.
Definition: flux/box/fickslaw.hh:51
static ComponentFluxVector flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const int phaseIdx, const ElementFluxVariablesCache &elemFluxVarsCache)
Returns the diffusive fluxes of all components within a fluid phase across the given sub-control volu...
Definition: flux/box/fickslaw.hh:92
Ficks's law for the box scheme in the context of coupled models where coupling occurs across the face...
Definition: multidomain/facet/box/fickslaw.hh:53
static std::vector< Scalar > calculateTransmissibilities(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const FluxVarCache &fluxVarCache, unsigned int phaseIdx)
Definition: multidomain/facet/box/fickslaw.hh:193
static ComponentFluxVector flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const int phaseIdx, const ElementFluxVarsCache &elemFluxVarCache)
Definition: multidomain/facet/box/fickslaw.hh:77
Declares all properties used in Dumux.
This file contains the data which is required to calculate diffusive fluxes due to molecular diffusio...