Line data Source code
1 : /*
2 : * Copyright (c) 2013-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__DARCY_VELOCITY_LINKER__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DARCY_VELOCITY_LINKER__
35 :
36 : #include "linker.h"
37 : #ifdef UG_FOR_LUA
38 : #include "bindings/lua/lua_user_data.h"
39 : #endif
40 :
41 : namespace ug{
42 :
43 :
44 : ////////////////////////////////////////////////////////////////////////////////
45 : // Darcy Velocity linker
46 : ////////////////////////////////////////////////////////////////////////////////
47 :
48 : /// Hard Coded Linker for the Darcy velocity
49 : /**
50 : * This linker computes the Darcy velocity \f$ \mathbf{q} = - \frac{\mathbf{K}}{\mu} ( \nabla p - \rho \mathbf{g} ) \f$,
51 : * where
52 : * <ul>
53 : * <li> \f$ \mathbf{K} \f$ permeability
54 : * <li> \f$ \mu \f$ viscosity
55 : * <li> \f$ \rho \f$ density
56 : * <li> \f$ \mathbf{g} \f$ gravity
57 : * <li> \f$ \nabla p \f$ pressure gradient
58 : * </ul>
59 : * are input parameters. This linker can be composed as a tree of other linkers,
60 : * but is computationally cheaper.
61 : */
62 : template <int dim>
63 : class DarcyVelocityLinker
64 : : public StdDataLinker< DarcyVelocityLinker<dim>, MathVector<dim>, dim>
65 : {
66 : /// Base class type
67 : typedef StdDataLinker< DarcyVelocityLinker<dim>, MathVector<dim>, dim> base_type;
68 :
69 : public:
70 0 : DarcyVelocityLinker() :
71 : m_spPermeability(NULL), m_spDPermeability(NULL),
72 : m_spViscosity(NULL), m_spDViscosity(NULL),
73 : m_spDensity(NULL), m_spDDensity(NULL),
74 : m_spGravity(NULL), m_spDGravity(NULL),
75 0 : m_spPressureGrad(NULL), m_spDPressureGrad(NULL), m_partialDerivMask(0)
76 : {
77 : // this linker needs exactly five input
78 : this->set_num_input(5);
79 0 : }
80 :
81 :
82 0 : inline void evaluate (MathVector<dim>& value,
83 : const MathVector<dim>& globIP,
84 : number time, int si) const
85 : {
86 : number density;
87 : number viscosity;
88 : MathVector<dim> gravity;
89 : MathVector<dim> pressureGrad;
90 : MathMatrix<dim,dim> permeability;
91 :
92 0 : (*m_spDensity)(density, globIP, time, si);
93 0 : (*m_spViscosity)(viscosity, globIP, time, si);
94 0 : (*m_spGravity)(gravity, globIP, time, si);
95 0 : (*m_spPressureGrad)(pressureGrad, globIP, time, si);
96 0 : (*m_spPermeability)(permeability, globIP, time, si);
97 :
98 : // Variables
99 : MathVector<dim> Vel;
100 :
101 : // compute rho*g
102 0 : VecScale(Vel, gravity, density);
103 :
104 : // compute rho*g - \nabla p
105 : VecSubtract(Vel, Vel, pressureGrad);
106 :
107 : // compute Darcy velocity q := K / mu * (rho*g - \nabla p)
108 : MatVecMult(value, permeability, Vel);
109 0 : VecScale(value, value, 1./viscosity);
110 0 : }
111 :
112 : template <int refDim>
113 0 : inline void evaluate(MathVector<dim> vValue[],
114 : const MathVector<dim> vGlobIP[],
115 : number time, int si,
116 : GridObject* elem,
117 : const MathVector<dim> vCornerCoords[],
118 : const MathVector<refDim> vLocIP[],
119 : const size_t nip,
120 : LocalVector* u,
121 : const MathMatrix<refDim, dim>* vJT = NULL) const
122 : {
123 0 : std::vector<number> vDensity(nip);
124 0 : std::vector<number> vViscosity(nip);
125 0 : std::vector<MathVector<dim> > vGravity(nip);
126 0 : std::vector<MathVector<dim> > vPressureGrad(nip);
127 0 : std::vector<MathMatrix<dim,dim> > vPermeability(nip);
128 :
129 0 : (*m_spDensity)(&vDensity[0], vGlobIP, time, si,
130 : elem, vCornerCoords, vLocIP, nip, u, vJT);
131 0 : (*m_spViscosity)(&vViscosity[0], vGlobIP, time, si,
132 : elem, vCornerCoords, vLocIP, nip, u, vJT);
133 0 : (*m_spGravity)(&vGravity[0], vGlobIP, time, si,
134 : elem, vCornerCoords, vLocIP, nip, u, vJT);
135 0 : (*m_spPressureGrad)(&vPressureGrad[0], vGlobIP, time, si,
136 : elem, vCornerCoords, vLocIP, nip, u, vJT);
137 0 : (*m_spPermeability)(&vPermeability[0], vGlobIP, time, si,
138 : elem, vCornerCoords, vLocIP, nip, u, vJT);
139 :
140 0 : for(size_t ip = 0; ip < nip; ++ip)
141 : {
142 : // Variables
143 : MathVector<dim> Vel;
144 :
145 : // compute rho*g
146 0 : VecScale(Vel, vGravity[ip], vDensity[ip]);
147 :
148 : // compute rho*g - \nabla p
149 : VecSubtract(Vel, Vel, vPressureGrad[ip]);
150 :
151 : // compute Darcy velocity q := K / mu * (rho*g - \nabla p)
152 0 : MatVecMult(vValue[ip], vPermeability[ip], Vel);
153 0 : VecScale(vValue[ip], vValue[ip], 1./vViscosity[ip]);
154 : }
155 0 : }
156 :
157 : template <int refDim>
158 0 : void eval_and_deriv(MathVector<dim> vDarcyVel[],
159 : const MathVector<dim> vGlobIP[],
160 : number time, int si,
161 : GridObject* elem,
162 : const MathVector<dim> vCornerCoords[],
163 : const MathVector<refDim> vLocIP[],
164 : const size_t nip,
165 : LocalVector* u,
166 : bool bDeriv,
167 : int s,
168 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[],
169 : const MathMatrix<refDim, dim>* vJT = NULL) const
170 : {
171 : // get the data of the ip series
172 0 : const number* vDensity = m_spDensity->values(s);
173 : const number* vViscosity = m_spViscosity->values(s);
174 : const MathVector<dim>* vGravity = m_spGravity->values(s);
175 : const MathVector<dim>* vPressureGrad = m_spPressureGrad->values(s);
176 : const MathMatrix<dim,dim>* vPermeability = m_spPermeability->values(s);
177 :
178 0 : for(size_t ip = 0; ip < nip; ++ip)
179 : {
180 : // Variables
181 : MathVector<dim> Vel;
182 :
183 : // compute rho*g
184 0 : VecScale(Vel, vGravity[ip], vDensity[ip]);
185 :
186 : // compute rho*g - \nabla p
187 0 : VecSubtract(Vel, Vel, vPressureGrad[ip]);
188 :
189 : // compute Darcy velocity q := K / mu * (rho*g - \nabla p)
190 0 : MatVecMult(vDarcyVel[ip], vPermeability[ip], Vel);
191 0 : VecScale(vDarcyVel[ip], vDarcyVel[ip], 1./vViscosity[ip]);
192 : }
193 :
194 : // Compute the derivatives at all ips //
195 : /////////////////////////////////////////////
196 :
197 : // check if something to do
198 0 : if(!bDeriv || this->zero_derivative()) return;
199 :
200 : // clear all derivative values
201 0 : this->set_zero(vvvDeriv, nip);
202 :
203 : // Derivatives of Viscosity
204 0 : if(m_partialDerivMask == 0 && m_spDViscosity.valid() && !m_spDViscosity->zero_derivative())
205 0 : for(size_t ip = 0; ip < nip; ++ip)
206 0 : for(size_t fct = 0; fct < m_spDViscosity->num_fct(); ++fct)
207 : {
208 : // get derivative of viscosity w.r.t. to all functions
209 : const number* vDViscosity = m_spDViscosity->deriv(s, ip, fct);
210 :
211 : // get common fct id for this function
212 : const size_t commonFct = this->input_common_fct(_MU_, fct);
213 :
214 : // loop all shapes and set the derivative
215 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
216 : {
217 : // DarcyVel_fct[sh] -= mu_fct_sh / mu * q
218 0 : VecScaleAppend(vvvDeriv[ip][commonFct][sh], -vDViscosity[sh] / vViscosity[ip], vDarcyVel[ip]);
219 : }
220 : }
221 :
222 : // Derivatives of Density
223 0 : if( m_partialDerivMask == 0 && m_spDDensity.valid() && !m_spDDensity->zero_derivative())
224 0 : for(size_t ip = 0; ip < nip; ++ip)
225 0 : for(size_t fct = 0; fct < m_spDDensity->num_fct(); ++fct)
226 : {
227 : // get derivative of viscosity w.r.t. to all functions
228 : const number* vDDensity = m_spDDensity->deriv(s, ip, fct);
229 :
230 : // get common fct id for this function
231 : const size_t commonFct = this->input_common_fct(_RHO_, fct);
232 :
233 : // Precompute K/mu * g
234 : MathVector<dim> Kmug;
235 :
236 : // a) compute K * g
237 0 : MatVecMult(Kmug, vPermeability[ip], vGravity[ip]);
238 :
239 : // b) compute K* g / mu
240 0 : VecScale(Kmug, Kmug, 1./vViscosity[ip]);
241 :
242 : // loop all shapes and set the derivative
243 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
244 : {
245 : UG_ASSERT(commonFct < vvvDeriv[ip].size(), commonFct<<", "<<vvvDeriv[ip].size());
246 : UG_ASSERT(sh < vvvDeriv[ip][commonFct].size(), sh<<", "<<vvvDeriv[ip][commonFct].size());
247 : // DarcyVel_fct[sh] += K/mu * (rho_fct_sh * g)
248 0 : VecScaleAppend(vvvDeriv[ip][commonFct][sh],
249 0 : vDDensity[sh], Kmug);
250 : }
251 : }
252 :
253 : // Derivatives of Gravity
254 0 : if(m_spDGravity.valid() && !m_spDGravity->zero_derivative())
255 0 : for(size_t ip = 0; ip < nip; ++ip)
256 0 : for(size_t fct = 0; fct < m_spDGravity->num_fct(); ++fct)
257 : {
258 : // get derivative of viscosity w.r.t. to all functions
259 : const MathVector<dim>* vDGravity = m_spDGravity->deriv(s, ip, fct);
260 :
261 : // get common fct id for this function
262 : const size_t commonFct = this->input_common_fct(_G_, fct);
263 :
264 : // Precompute K/mu * rho
265 : MathMatrix<dim,dim> Kmurho;
266 :
267 : // a) compute K/mu * rho
268 0 : MatScale(Kmurho, vDensity[ip]/vViscosity[ip],vPermeability[ip]);
269 :
270 : // loop all shapes and set the derivative
271 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
272 : {
273 : MathVector<dim> tmp;
274 0 : MatVecMult(tmp, Kmurho, vDGravity[sh]);
275 :
276 0 : vvvDeriv[ip][commonFct][sh] += tmp;
277 : }
278 : }
279 :
280 : // Derivatives of Pressure
281 0 : if(m_partialDerivMask ==0 && m_spDPressureGrad.valid() && !m_spDPressureGrad->zero_derivative() )
282 0 : for(size_t ip = 0; ip < nip; ++ip)
283 0 : for(size_t fct = 0; fct < m_spDPressureGrad->num_fct(); ++fct)
284 : {
285 : // get derivative of viscosity w.r.t. to all functions
286 : const MathVector<dim>* vDPressureGrad = m_spDPressureGrad->deriv(s, ip, fct);
287 :
288 : // get common fct id for this function
289 : const size_t commonFct = this->input_common_fct(_DP_, fct);
290 :
291 : // Precompute -K/mu
292 : MathMatrix<dim,dim> Kmu;
293 :
294 : // a) compute -K/mu
295 0 : MatScale(Kmu, -1.0/vViscosity[ip],vPermeability[ip]);
296 :
297 : // loop all shapes and set the derivative
298 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
299 : {
300 : MathVector<dim> tmp;
301 0 : MatVecMult(tmp, Kmu, vDPressureGrad[sh]);
302 :
303 0 : vvvDeriv[ip][commonFct][sh] += tmp;
304 : }
305 : }
306 :
307 : // Derivatives of Permeability
308 0 : if(m_spDPermeability.valid() && !m_spDPermeability->zero_derivative())
309 0 : for(size_t ip = 0; ip < nip; ++ip)
310 0 : for(size_t fct = 0; fct < m_spDPermeability->num_fct(); ++fct)
311 : {
312 : // get derivative of viscosity w.r.t. to all functions
313 : const MathMatrix<dim,dim>* vDPermeability = m_spDPermeability->deriv(s, ip, fct);
314 :
315 : // get common fct id for this function
316 : const size_t commonFct = this->input_common_fct(_K_, fct);
317 :
318 : // Variables
319 : MathVector<dim> Vel;
320 :
321 : // compute rho*g
322 0 : VecScale(Vel, vGravity[ip], vDensity[ip]);
323 :
324 : // compute rho*g - \nabla p
325 0 : VecSubtract(Vel, Vel, vPressureGrad[ip]);
326 :
327 : // compute Darcy velocity q := K / mu * (rho*g - \nabla p)
328 0 : VecScale(Vel, Vel, 1./vViscosity[ip]);
329 :
330 : // loop all shapes and set the derivative
331 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
332 : {
333 : MathVector<dim> tmp;
334 0 : MatVecMult(tmp, vDPermeability[sh], Vel);
335 :
336 0 : vvvDeriv[ip][commonFct][sh] += tmp;
337 : }
338 : }
339 : }
340 :
341 : public:
342 : /// set permeability import
343 0 : void set_permeability(SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > data)
344 : {
345 0 : m_spPermeability = data;
346 0 : m_spDPermeability = data.template cast_dynamic<DependentUserData<MathMatrix<dim,dim>, dim> >();
347 0 : base_type::set_input(_K_, data, data);
348 0 : }
349 :
350 0 : void set_permeability(number val)
351 : {
352 0 : set_permeability(make_sp(new ConstUserMatrix<dim>(val)));
353 0 : }
354 :
355 : #ifdef UG_FOR_LUA
356 0 : void set_permeability(const char* fctName)
357 : {
358 0 : set_permeability(LuaUserDataFactory<MathMatrix<dim,dim>, dim>::create(fctName));
359 0 : }
360 :
361 0 : void set_permeability(LuaFunctionHandle fct)
362 : {
363 0 : set_permeability(make_sp(new LuaUserData<MathMatrix<dim,dim>, dim>(fct)));
364 0 : }
365 : #endif
366 :
367 : /// set permeability import
368 0 : void set_viscosity(SmartPtr<CplUserData<number, dim> > data)
369 : {
370 0 : m_spViscosity = data;
371 0 : m_spDViscosity = data.template cast_dynamic<DependentUserData<number, dim> >();
372 0 : base_type::set_input(_MU_, data, data);
373 0 : }
374 :
375 0 : void set_viscosity(number val)
376 : {
377 0 : set_viscosity(make_sp(new ConstUserNumber<dim>(val)));
378 0 : }
379 :
380 : /// set density import
381 0 : void set_density(SmartPtr<CplUserData<number, dim> > data)
382 : {
383 0 : m_spDensity = data;
384 0 : m_spDDensity = data.template cast_dynamic<DependentUserData<number, dim> >();
385 0 : base_type::set_input(_RHO_, data, data);
386 0 : }
387 :
388 0 : void set_density(number val)
389 : {
390 0 : set_density(make_sp(new ConstUserNumber<dim>(val)));
391 0 : }
392 :
393 : /// set gravity import
394 0 : void set_gravity(SmartPtr<CplUserData<MathVector<dim>, dim> > data)
395 : {
396 0 : m_spGravity = data;
397 0 : m_spDGravity = data.template cast_dynamic<DependentUserData<MathVector<dim>, dim> >();
398 0 : base_type::set_input(_G_, data, data);
399 0 : }
400 :
401 : /// set pressure gradient import
402 0 : void set_pressure_gradient(SmartPtr<CplUserData<MathVector<dim>, dim> > data)
403 : {
404 0 : m_spPressureGrad = data;
405 0 : m_spDPressureGrad = data.template cast_dynamic<DependentUserData<MathVector<dim>, dim> >();
406 0 : base_type::set_input(_DP_, data, data);
407 0 : }
408 :
409 : protected:
410 : /// import for permeability
411 : static const size_t _K_ = 0;
412 : SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > m_spPermeability;
413 : SmartPtr<DependentUserData<MathMatrix<dim,dim>, dim> > m_spDPermeability;
414 :
415 : /// import for viscosity
416 : static const size_t _MU_ = 1;
417 : SmartPtr<CplUserData<number, dim> > m_spViscosity;
418 : SmartPtr<DependentUserData<number, dim> > m_spDViscosity;
419 :
420 : /// import for density
421 : static const size_t _RHO_ = 2;
422 : SmartPtr<CplUserData<number, dim> > m_spDensity;
423 : SmartPtr<DependentUserData<number, dim> > m_spDDensity;
424 :
425 : /// import for gravity
426 : static const size_t _G_ = 3;
427 : SmartPtr<CplUserData<MathVector<dim>, dim> > m_spGravity;
428 : SmartPtr<DependentUserData<MathVector<dim>, dim> > m_spDGravity;
429 :
430 : /// import for pressure gradient
431 : static const size_t _DP_ = 4;
432 : SmartPtr<CplUserData<MathVector<dim>, dim> > m_spPressureGrad;
433 : SmartPtr<DependentUserData<MathVector<dim>, dim> > m_spDPressureGrad;
434 :
435 :
436 : public:
437 0 : void set_derivative_mask(int mask) {
438 0 : m_partialDerivMask = mask;
439 0 : std::cerr << "Setting some derivatives: "<< m_partialDerivMask << "(" << this <<")" << std::endl;
440 0 : }
441 :
442 : protected:
443 : // disable certain derivatives
444 : int m_partialDerivMask;
445 : };
446 :
447 : } // end namespace ug
448 :
449 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DARCY_VELOCITY_LINKER__ */
|
