ug4
darcy_velocity_linker.h
Go to the documentation of this file.
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
39 #endif
40 
41 namespace ug{
42 
43 
45 // Darcy Velocity linker
47 
49 
62 template <int dim>
64  : public StdDataLinker< DarcyVelocityLinker<dim>, MathVector<dim>, dim>
65 {
68 
69  public:
72  m_spViscosity(NULL), m_spDViscosity(NULL),
73  m_spDensity(NULL), m_spDDensity(NULL),
74  m_spGravity(NULL), m_spDGravity(NULL),
76  {
77  // this linker needs exactly five input
78  this->set_num_input(5);
79  }
80 
81 
82  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  (*m_spDensity)(density, globIP, time, si);
93  (*m_spViscosity)(viscosity, globIP, time, si);
94  (*m_spGravity)(gravity, globIP, time, si);
95  (*m_spPressureGrad)(pressureGrad, globIP, time, si);
96  (*m_spPermeability)(permeability, globIP, time, si);
97 
98  // Variables
99  MathVector<dim> Vel;
100 
101  // compute rho*g
102  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  VecScale(value, value, 1./viscosity);
110  }
111 
112  template <int refDim>
113  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  std::vector<number> vDensity(nip);
124  std::vector<number> vViscosity(nip);
125  std::vector<MathVector<dim> > vGravity(nip);
126  std::vector<MathVector<dim> > vPressureGrad(nip);
127  std::vector<MathMatrix<dim,dim> > vPermeability(nip);
128 
129  (*m_spDensity)(&vDensity[0], vGlobIP, time, si,
130  elem, vCornerCoords, vLocIP, nip, u, vJT);
131  (*m_spViscosity)(&vViscosity[0], vGlobIP, time, si,
132  elem, vCornerCoords, vLocIP, nip, u, vJT);
133  (*m_spGravity)(&vGravity[0], vGlobIP, time, si,
134  elem, vCornerCoords, vLocIP, nip, u, vJT);
135  (*m_spPressureGrad)(&vPressureGrad[0], vGlobIP, time, si,
136  elem, vCornerCoords, vLocIP, nip, u, vJT);
137  (*m_spPermeability)(&vPermeability[0], vGlobIP, time, si,
138  elem, vCornerCoords, vLocIP, nip, u, vJT);
139 
140  for(size_t ip = 0; ip < nip; ++ip)
141  {
142  // Variables
143  MathVector<dim> Vel;
144 
145  // compute rho*g
146  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  MatVecMult(vValue[ip], vPermeability[ip], Vel);
153  VecScale(vValue[ip], vValue[ip], 1./vViscosity[ip]);
154  }
155  }
156 
157  template <int refDim>
158  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  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  for(size_t ip = 0; ip < nip; ++ip)
179  {
180  // Variables
181  MathVector<dim> Vel;
182 
183  // compute rho*g
184  VecScale(Vel, vGravity[ip], vDensity[ip]);
185 
186  // compute rho*g - \nabla p
187  VecSubtract(Vel, Vel, vPressureGrad[ip]);
188 
189  // compute Darcy velocity q := K / mu * (rho*g - \nabla p)
190  MatVecMult(vDarcyVel[ip], vPermeability[ip], Vel);
191  VecScale(vDarcyVel[ip], vDarcyVel[ip], 1./vViscosity[ip]);
192  }
193 
194  // Compute the derivatives at all ips //
196 
197  // check if something to do
198  if(!bDeriv || this->zero_derivative()) return;
199 
200  // clear all derivative values
201  this->set_zero(vvvDeriv, nip);
202 
203  // Derivatives of Viscosity
204  if(m_partialDerivMask == 0 && m_spDViscosity.valid() && !m_spDViscosity->zero_derivative())
205  for(size_t ip = 0; ip < nip; ++ip)
206  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  for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
216  {
217  // DarcyVel_fct[sh] -= mu_fct_sh / mu * q
218  VecScaleAppend(vvvDeriv[ip][commonFct][sh], -vDViscosity[sh] / vViscosity[ip], vDarcyVel[ip]);
219  }
220  }
221 
222  // Derivatives of Density
223  if( m_partialDerivMask == 0 && m_spDDensity.valid() && !m_spDDensity->zero_derivative())
224  for(size_t ip = 0; ip < nip; ++ip)
225  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  MatVecMult(Kmug, vPermeability[ip], vGravity[ip]);
238 
239  // b) compute K* g / mu
240  VecScale(Kmug, Kmug, 1./vViscosity[ip]);
241 
242  // loop all shapes and set the derivative
243  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  VecScaleAppend(vvvDeriv[ip][commonFct][sh],
249  vDDensity[sh], Kmug);
250  }
251  }
252 
253  // Derivatives of Gravity
254  if(m_spDGravity.valid() && !m_spDGravity->zero_derivative())
255  for(size_t ip = 0; ip < nip; ++ip)
256  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  MatScale(Kmurho, vDensity[ip]/vViscosity[ip],vPermeability[ip]);
269 
270  // loop all shapes and set the derivative
271  for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
272  {
273  MathVector<dim> tmp;
274  MatVecMult(tmp, Kmurho, vDGravity[sh]);
275 
276  vvvDeriv[ip][commonFct][sh] += tmp;
277  }
278  }
279 
280  // Derivatives of Pressure
281  if(m_partialDerivMask ==0 && m_spDPressureGrad.valid() && !m_spDPressureGrad->zero_derivative() )
282  for(size_t ip = 0; ip < nip; ++ip)
283  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
293 
294  // a) compute -K/mu
295  MatScale(Kmu, -1.0/vViscosity[ip],vPermeability[ip]);
296 
297  // loop all shapes and set the derivative
298  for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
299  {
300  MathVector<dim> tmp;
301  MatVecMult(tmp, Kmu, vDPressureGrad[sh]);
302 
303  vvvDeriv[ip][commonFct][sh] += tmp;
304  }
305  }
306 
307  // Derivatives of Permeability
308  if(m_spDPermeability.valid() && !m_spDPermeability->zero_derivative())
309  for(size_t ip = 0; ip < nip; ++ip)
310  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  VecScale(Vel, vGravity[ip], vDensity[ip]);
323 
324  // compute rho*g - \nabla p
325  VecSubtract(Vel, Vel, vPressureGrad[ip]);
326 
327  // compute Darcy velocity q := K / mu * (rho*g - \nabla p)
328  VecScale(Vel, Vel, 1./vViscosity[ip]);
329 
330  // loop all shapes and set the derivative
331  for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
332  {
333  MathVector<dim> tmp;
334  MatVecMult(tmp, vDPermeability[sh], Vel);
335 
336  vvvDeriv[ip][commonFct][sh] += tmp;
337  }
338  }
339  }
340 
341  public:
344  {
345  m_spPermeability = data;
346  m_spDPermeability = data.template cast_dynamic<DependentUserData<MathMatrix<dim,dim>, dim> >();
347  base_type::set_input(_K_, data, data);
348  }
349 
351  {
353  }
354 
355  #ifdef UG_FOR_LUA
356  void set_permeability(const char* fctName)
357  {
359  }
360 
361  void set_permeability(LuaFunctionHandle fct)
362  {
363  set_permeability(make_sp(new LuaUserData<MathMatrix<dim,dim>, dim>(fct)));
364  }
365  #endif
366 
369  {
370  m_spViscosity = data;
371  m_spDViscosity = data.template cast_dynamic<DependentUserData<number, dim> >();
372  base_type::set_input(_MU_, data, data);
373  }
374 
376  {
378  }
379 
382  {
383  m_spDensity = data;
384  m_spDDensity = data.template cast_dynamic<DependentUserData<number, dim> >();
385  base_type::set_input(_RHO_, data, data);
386  }
387 
388  void set_density(number val)
389  {
391  }
392 
395  {
396  m_spGravity = data;
397  m_spDGravity = data.template cast_dynamic<DependentUserData<MathVector<dim>, dim> >();
398  base_type::set_input(_G_, data, data);
399  }
400 
403  {
404  m_spPressureGrad = data;
405  m_spDPressureGrad = data.template cast_dynamic<DependentUserData<MathVector<dim>, dim> >();
406  base_type::set_input(_DP_, data, data);
407  }
408 
409  protected:
411  static const size_t _K_ = 0;
414 
416  static const size_t _MU_ = 1;
419 
421  static const size_t _RHO_ = 2;
424 
426  static const size_t _G_ = 3;
429 
431  static const size_t _DP_ = 4;
434 
435 
436  public:
437  void set_derivative_mask(int mask) {
438  m_partialDerivMask = mask;
439  std::cerr << "Setting some derivatives: "<< m_partialDerivMask << "(" << this <<")" << std::endl;
440  }
441 
442  protected:
443  // disable certain derivatives
445 };
446 
447 } // end namespace ug
448 
449 #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DARCY_VELOCITY_LINKER__ */
parameterString s
Definition: smart_pointer.h:108
constant matrix user data
Definition: const_user_data.h:232
constant scalar user data
Definition: const_user_data.h:153
Type based UserData.
Definition: user_data.h:501
const TData & value(size_t s, size_t ip) const
returns the value at ip
Definition: user_data.h:512
Hard Coded Linker for the Darcy velocity.
Definition: darcy_velocity_linker.h:65
SmartPtr< CplUserData< MathVector< dim >, dim > > m_spPressureGrad
Definition: darcy_velocity_linker.h:432
void set_density(number val)
Definition: darcy_velocity_linker.h:388
void set_pressure_gradient(SmartPtr< CplUserData< MathVector< dim >, dim > > data)
set pressure gradient import
Definition: darcy_velocity_linker.h:402
SmartPtr< CplUserData< MathMatrix< dim, dim >, dim > > m_spPermeability
Definition: darcy_velocity_linker.h:412
SmartPtr< CplUserData< MathVector< dim >, dim > > m_spGravity
Definition: darcy_velocity_linker.h:427
void set_viscosity(number val)
Definition: darcy_velocity_linker.h:375
SmartPtr< DependentUserData< MathVector< dim >, dim > > m_spDPressureGrad
Definition: darcy_velocity_linker.h:433
void evaluate(MathVector< dim > vValue[], const MathVector< dim > vGlobIP[], number time, int si, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< refDim > vLocIP[], const size_t nip, LocalVector *u, const MathMatrix< refDim, dim > *vJT=NULL) const
Definition: darcy_velocity_linker.h:113
void set_permeability(SmartPtr< CplUserData< MathMatrix< dim, dim >, dim > > data)
set permeability import
Definition: darcy_velocity_linker.h:343
static const size_t _G_
import for gravity
Definition: darcy_velocity_linker.h:426
void set_viscosity(SmartPtr< CplUserData< number, dim > > data)
set permeability import
Definition: darcy_velocity_linker.h:368
SmartPtr< CplUserData< number, dim > > m_spViscosity
Definition: darcy_velocity_linker.h:417
SmartPtr< CplUserData< number, dim > > m_spDensity
Definition: darcy_velocity_linker.h:422
void evaluate(MathVector< dim > &value, const MathVector< dim > &globIP, number time, int si) const
Definition: darcy_velocity_linker.h:82
void set_derivative_mask(int mask)
Definition: darcy_velocity_linker.h:437
static const size_t _DP_
import for pressure gradient
Definition: darcy_velocity_linker.h:431
void set_permeability(number val)
Definition: darcy_velocity_linker.h:350
static const size_t _RHO_
import for density
Definition: darcy_velocity_linker.h:421
static const size_t _K_
import for permeability
Definition: darcy_velocity_linker.h:411
SmartPtr< DependentUserData< number, dim > > m_spDDensity
Definition: darcy_velocity_linker.h:423
StdDataLinker< DarcyVelocityLinker< dim >, MathVector< dim >, dim > base_type
Base class type.
Definition: darcy_velocity_linker.h:67
void set_gravity(SmartPtr< CplUserData< MathVector< dim >, dim > > data)
set gravity import
Definition: darcy_velocity_linker.h:394
static const size_t _MU_
import for viscosity
Definition: darcy_velocity_linker.h:416
SmartPtr< DependentUserData< MathMatrix< dim, dim >, dim > > m_spDPermeability
Definition: darcy_velocity_linker.h:413
void eval_and_deriv(MathVector< dim > vDarcyVel[], const MathVector< dim > vGlobIP[], number time, int si, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< refDim > vLocIP[], const size_t nip, LocalVector *u, bool bDeriv, int s, std::vector< std::vector< MathVector< dim > > > vvvDeriv[], const MathMatrix< refDim, dim > *vJT=NULL) const
Definition: darcy_velocity_linker.h:158
void set_density(SmartPtr< CplUserData< number, dim > > data)
set density import
Definition: darcy_velocity_linker.h:381
int m_partialDerivMask
Definition: darcy_velocity_linker.h:444
DarcyVelocityLinker()
Definition: darcy_velocity_linker.h:70
SmartPtr< DependentUserData< MathVector< dim >, dim > > m_spDGravity
Definition: darcy_velocity_linker.h:428
SmartPtr< DependentUserData< number, dim > > m_spDViscosity
Definition: darcy_velocity_linker.h:418
The base class for all geometric objects, such as vertices, edges, faces, volumes,...
Definition: grid_base_objects.h:157
const MathVector< dim > & ip(size_t s, size_t ip) const
returns global ip
Definition: user_data.h:401
number time() const
get the current evaluation time
Definition: user_data.h:285
Definition: local_algebra.h:198
Factory providing LuaUserData.
Definition: lua_user_data.h:180
combines several UserDatas to a new UserData of a specified type
Definition: linker.h:54
void set_num_input(size_t num)
sets the number of inputs
Definition: linker.h:107
size_t input_common_fct(size_t i, size_t fct) const
returns the number in the common FctGrp for a fct of an input
Definition: linker.h:153
virtual void set_input(size_t i, SmartPtr< ICplUserData< dim > > input, SmartPtr< UserDataInfo > info)
sets an input
Definition: linker.h:114
virtual bool zero_derivative() const
returns if derivative is zero
Definition: linker_impl.h:179
void MatScale(matrix_t &mOut, typename matrix_t::value_type s, const matrix_t &m)
scales a matrix_t
Definition: math_matrix_functions_common_impl.hpp:317
static const int dim
#define UG_ASSERT(expr, msg)
Definition: assert.h:70
double number
Definition: types.h:124
void MatVecMult(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Matrix - Vector Multiplication.
Definition: math_matrix_vector_functions_common_impl.hpp:49
void VecScaleAppend(vector_t &vOut, typename vector_t::value_type s1, const vector_t &v1)
Scales a Vector and adds it to a second vector.
Definition: math_vector_functions_common_impl.hpp:126
void VecSubtract(vector_t &vOut, const vector_t &v1, const vector_t &v2)
subtracts v2 from v1 and stores the result in a vOut
Definition: math_vector_functions_common_impl.hpp:226
void VecScale(vector_t &vOut, const vector_t &v, typename vector_t::value_type s)
scales a MathVector<N>
Definition: math_vector_functions_common_impl.hpp:252
function util LuaCallbackHelper create(func)
the ug namespace
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition: smart_pointer.h:836