Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONV_SHAPE_INTERFACE__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONV_SHAPE_INTERFACE__
35 :
36 : #include "common/common.h"
37 : #include "common/util/provider.h"
38 : #include "common/math/ugmath_types.h"
39 : #include "lib_disc/spatial_disc/disc_util/fv_geom_base.h"
40 :
41 : namespace ug{
42 :
43 : /////////////////////////////////////////////////////////////////////////////
44 : // Interface for Convection shapes
45 : /////////////////////////////////////////////////////////////////////////////
46 :
47 : /// Interface class for upwind methods
48 : /**
49 : * Upwind is a way of stabilization of discretizations of convection(-diffusion)
50 : * PDEs. For a detailed description of the idea, cf.
51 : * <ul>
52 : * <li> P. Frolkovic, H. De Schepper, Numerical Modelling of Convection Dominated
53 : * Transport Coupled with Density Driven Flow in Porous Media.
54 : * Advances in Water Resources 24 (1): 63–72, 2000, DOI: 10.1016/S0309-1708(00)00025-7 </li>
55 : * <li> K. W. Morton, D. Mayers, Numerical Solution of Partial Differential Equations:
56 : * An Introduction (2nd ed.), Cambridge: Cambridge University Press, 2005, DOI: 10.1017/CBO9780511812248 </li>
57 : * <li> R. J. LeVeque, Numerical Methods for Conservation Laws (2nd ed.), Springer Basel AG, 1992, DOI: 10.1007/978-3-0348-8629-1 </li>
58 : * </ul>
59 : * Convection shapes is one of possible implementations of upwind methods,
60 : * in which the upwind is considered as a special kind of interpolation of degrees
61 : * of freedom of the unknown function (for ex., concentration at the corners of
62 : * the element) into integration points inside the grid element. E.g. the so-called
63 : * 'full upwind' method for the vertex-centered FV discretizations sets the
64 : * concentration in the whole grid element be equal to the concentration at a
65 : * certain specially selected 'upwind corner' chosen according to the direction
66 : * of the convection velocity. Like every interpolation, this one can be introduced
67 : * by a set of shape functions associated with the degrees of freedom of the
68 : * grid function. These shape functions for an upwind method are said to be
69 : * its convection shapes.
70 : *
71 : * But in terms of this interface, the convection shapes are the values of the
72 : * convection shape functions at the integration points premultiplied by the norm
73 : * of the projection of the velocity to the normal to the corresponding
74 : * subcontrol volume faces and by the area of these faces. Thus, for a convective
75 : * term of the form
76 : * \f[ \nabla \cdot (\mathbf{v} c) \f]
77 : * (with \f$\mathbf{v}\f$ being the convective velocity) the convective flux
78 : * through a subcontrol volume face with the integration point \f$ip\f$ is
79 : * \f[ \sum_{s=1}^{n_s} \phi_s (ip) \, c_s, \f]
80 : * where \f$\phi_s (ip)\f$ are the values of the \f$n_s\f$ ``convection shapes''
81 : * at \f$ip\f$, whereas \f$c_s\f$ are the degrees of freedom of \f$c\f$ corresponding
82 : * to the shapes.
83 : *
84 : * Although the convection shapes for every upwind method are based on the
85 : * same idea, they depend on the geometry of the secondary grid (i.e., on the
86 : * finite volume geometry). For this, for every FV geometry, an own object
87 : * of the convection shapes class is registered.
88 : *
89 : * Note furthermore that the convection shapes may depend not only on the
90 : * velocity but also on the diffusion. Such methods raise up the asymptotic
91 : * convergence of the discretization of convection-diffusion problems. As
92 : * the diffusion tensor may also depend on the unknowns (via its dispersion
93 : * part), the derivatives of the shapes w.r.t. the components of the diffusion
94 : * tensor are taken into account.
95 : *
96 : * \tparam dim the world dimension (in which the velocity is given)
97 : */
98 : template <int dim>
99 : class IConvectionShapes
100 : {
101 : public:
102 : /// Abbreviation for own type
103 : typedef IConvectionShapes<dim> this_type;
104 :
105 : /// type of update function
106 : typedef bool (this_type::*UpdateFunc)
107 : (const FVGeometryBase* geo,
108 : const MathVector<dim>* Velocity,
109 : const MathMatrix<dim, dim>* Diffusion,
110 : bool computeDeriv);
111 :
112 : public:
113 : /// constructor
114 0 : IConvectionShapes()
115 0 : : m_numScvf(0), m_numSh(0), m_bNonZeroDerivDiffusion(true)
116 : {
117 : m_vUpdateFunc.clear();
118 : m_vConvShape.clear();
119 : m_vConvShapeVel.clear();
120 : m_vConvShapeDiffusion.clear();
121 : }
122 :
123 : /// destructor
124 0 : virtual ~IConvectionShapes() {};
125 :
126 : /// returns number of shapes
127 0 : size_t num_sh() const {return m_numSh;}
128 :
129 : /// returns number of sub control volume faces
130 : size_t num_scvf() const {return m_numScvf;}
131 :
132 : /// shape value
133 : number operator()(size_t scvf, size_t sh) const
134 : {
135 : UG_ASSERT(scvf < m_vConvShape.size(), "Invalid index");
136 : UG_ASSERT(sh < m_vConvShape[scvf].size(), "Invalid index");
137 0 : return m_vConvShape[scvf][sh];
138 : }
139 :
140 : /// upwind shape for corner vel
141 : const MathVector<dim>& D_vel(size_t scvf, size_t sh) const
142 : {
143 : UG_ASSERT(scvf < m_vConvShapeVel.size(), "Invalid index");
144 : UG_ASSERT(sh < m_vConvShapeVel[scvf].size(), "Invalid index");
145 : return m_vConvShapeVel[scvf][sh];
146 : }
147 :
148 : /// returns if upwind shape w.r.t. ip vel is non-zero
149 0 : bool non_zero_deriv_diffusion() const {return m_bNonZeroDerivDiffusion;}
150 :
151 : /// upwind shapes for ip vel
152 : const MathMatrix<dim,dim>& D_diffusion(size_t scvf, size_t sh) const
153 : {
154 : UG_ASSERT(scvf < m_vConvShapeDiffusion.size(), "Invalid index");
155 : UG_ASSERT(sh < m_vConvShapeDiffusion[scvf].size(), "Invalid index");
156 : return m_vConvShapeDiffusion[scvf][sh];
157 : }
158 :
159 0 : bool update(const FVGeometryBase* geo,
160 : const MathVector<dim>* Velocity,
161 : const MathMatrix<dim, dim>* DiffDisp,
162 : bool computeDeriv)
163 0 : {return (this->*(m_vUpdateFunc[m_id])) (geo, Velocity, DiffDisp, computeDeriv);}
164 :
165 : //////////////////////////
166 : // internal handling
167 : //////////////////////////
168 :
169 : protected:
170 : /// resize the data arrays
171 : void set_sizes(size_t numScvf, size_t numSh);
172 :
173 : /// sets the shape ip flag
174 0 : void set_non_zero_deriv_diffusion_flag(bool flag) {m_bNonZeroDerivDiffusion = flag;}
175 :
176 : /// non-const access to ip velocity (i.e. interpolated velocity at ip)
177 : number& conv_shape(size_t scvf, size_t sh)
178 : {
179 : UG_ASSERT(scvf < m_vConvShape.size(), "Invalid index");
180 : UG_ASSERT(sh < m_vConvShape[scvf].size(), "Invalid index");
181 : return m_vConvShape[scvf][sh];
182 : }
183 :
184 : /// non-const access to upwind shapes for corner vel
185 : MathVector<dim>& D_vel(size_t scvf, size_t sh)
186 : {
187 : UG_ASSERT(scvf < m_vConvShapeVel.size(), "Invalid index");
188 : UG_ASSERT(sh < m_vConvShapeVel[scvf].size(), "Invalid index");
189 : return m_vConvShapeVel[scvf][sh];
190 : }
191 :
192 : /// non-const access to upwind shapes for ip vel
193 : MathMatrix<dim,dim>& conv_shape_diffusion(size_t scvf, size_t sh)
194 : {
195 : UG_ASSERT(scvf < m_vConvShapeDiffusion.size(), "Invalid index");
196 : UG_ASSERT(sh < m_vConvShapeDiffusion[scvf].size(), "Invalid index");
197 : return m_vConvShapeDiffusion[scvf][sh];
198 : }
199 :
200 :
201 : /// number of current scvf
202 : size_t m_numScvf;
203 :
204 : /// number of current shape functions (usually in corners)
205 : size_t m_numSh;
206 :
207 : /// upwind shapes for corners shape functions
208 : std::vector<std::vector<number> > m_vConvShape;
209 :
210 : /// upwind shapes for corners shape functions
211 : std::vector<std::vector<MathVector<dim> > > m_vConvShapeVel;
212 :
213 : /// upwind shapes for corners shape functions
214 : std::vector<std::vector<MathMatrix<dim,dim> > > m_vConvShapeDiffusion;
215 :
216 : /// flag if ip shapes are non-zero
217 : bool m_bNonZeroDerivDiffusion;
218 :
219 : //////////////////////////
220 : // registering process
221 : //////////////////////////
222 :
223 : public:
224 : /// register a update function for a Geometry
225 : template <typename TFVGeom, typename TAssFunc>
226 : void register_update_func(TAssFunc func);
227 :
228 : /// set the Geometry type to use for next updates
229 : template <typename TFVGeom>
230 : bool set_geometry_type(const TFVGeom& geo);
231 :
232 : protected:
233 : /// Vector holding all update functions
234 : std::vector<UpdateFunc> m_vUpdateFunc;
235 :
236 : /// id of current geometry type
237 : int m_id;
238 : };
239 :
240 : /////////////////////////////////////////////////////////////////////////////
241 : // Interface for Stabilization
242 : /////////////////////////////////////////////////////////////////////////////
243 :
244 : // register a update function for a Geometry
245 : template <int dim>
246 : template <typename TFVGeom, typename TAssFunc>
247 : void
248 0 : IConvectionShapes<dim>::
249 : register_update_func(TAssFunc func)
250 : {
251 : // get unique geometry id
252 0 : size_t id = GetUniqueFVGeomID<TFVGeom>();
253 :
254 : // make sure that there is enough space
255 0 : if((size_t)id >= m_vUpdateFunc.size())
256 0 : m_vUpdateFunc.resize(id+1, NULL);
257 :
258 : // set pointer
259 0 : m_vUpdateFunc[id] = (UpdateFunc)func;
260 0 : }
261 :
262 : // set the Geometry type to use for next updates
263 : template <int dim>
264 : template <typename TFVGeom>
265 : bool
266 0 : IConvectionShapes<dim>::
267 : set_geometry_type(const TFVGeom& geo)
268 : {
269 : // get unique geometry id
270 0 : size_t id = GetUniqueFVGeomID<TFVGeom>();
271 :
272 : // check that function exists
273 0 : if(id >= m_vUpdateFunc.size() || m_vUpdateFunc[id] == NULL)
274 : {
275 : UG_LOG("ERROR in 'IConvectionShapes::set_geometry_type':"
276 : " No update function registered for this Geometry.\n");
277 0 : return false;
278 : }
279 :
280 : // set current geometry
281 0 : m_id = id;
282 :
283 : // set sizes
284 0 : set_sizes(geo.num_scvf(), geo.num_sh());
285 :
286 : // we're done
287 0 : return true;
288 : }
289 :
290 :
291 : // resize the data arrays
292 : template <int dim>
293 : void
294 0 : IConvectionShapes<dim>::
295 : set_sizes(size_t numScvf, size_t numSh)
296 : {
297 : // remember sizes
298 0 : m_numScvf = numScvf;
299 0 : m_numSh = numSh;
300 :
301 : // adjust arrays
302 0 : m_vConvShape.resize(m_numScvf);
303 0 : m_vConvShapeVel.resize(m_numScvf);
304 0 : m_vConvShapeDiffusion.resize(m_numScvf);
305 0 : for(size_t scvf = 0; scvf < m_numScvf; ++scvf)
306 : {
307 0 : m_vConvShape[scvf].resize(m_numSh);
308 0 : m_vConvShapeVel[scvf].resize(m_numSh);
309 0 : m_vConvShapeDiffusion[scvf].resize(m_numSh);
310 : }
311 0 : }
312 :
313 :
314 :
315 : } // end namespace ug
316 :
317 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__DENSITY_DRIVEN_FLOW__FV1__CONV_SHAPE_INTERFACE__ */
|
