Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Jonas Simon
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__BINGHAM_VISCOSITY_LINKER__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__BINGHAM_VISCOSITY_LINKER__
35 :
36 : #include "linker.h"
37 :
38 : namespace ug{
39 :
40 :
41 : ////////////////////////////////////////////////////////////////////////////////
42 : // Bingham Viscosity linker
43 : ////////////////////////////////////////////////////////////////////////////////
44 :
45 : /// Linker for the Bingham viscosity
46 : /**
47 : * This linker computes the Bingham viscosity \f$ \mathbf{q} = - \frac{\mathbf{K}}{\mu} ( \nabla p - \rho \mathbf{g} ) \f$,
48 : * where
49 : * <ul>
50 : * <li> \f$ \mathbf{K} \f$ permeability
51 : * <li> \f$ \mu \f$ viscosity
52 : * <li> \f$ \rho \f$ density
53 : * <li> \f$ \mathbf{g} \f$ gravity
54 : * <li> \f$ \nabla p \f$ pressure gradient
55 : * </ul>
56 : * are input parameters.
57 : */
58 :
59 : template <int dim>
60 : class BinghamViscosityLinker
61 : : public StdDataLinker< BinghamViscosityLinker<dim>, number, dim>
62 : {
63 : // Base class type
64 : typedef StdDataLinker< BinghamViscosityLinker<dim>, number, dim> base_type;
65 :
66 : // Constructor
67 : public:
68 0 : BinghamViscosityLinker() :
69 : m_spDensity(NULL), m_spDDensity(NULL),
70 : m_spViscosity(NULL), m_spDViscosity(NULL),
71 : m_spYieldStress(NULL), m_spDYieldStress(NULL),
72 0 : m_spVelocityGrad(NULL), m_spDVelocityGrad(NULL)
73 : {
74 : // this linker needs exactly four input
75 : this->set_num_input(4);
76 0 : }
77 :
78 : // function for evaluation at single ip?
79 0 : inline void evaluate (number& value,
80 : const MathVector<dim>& globIP,
81 : number time, int si) const
82 : {
83 : UG_LOG("BinghamViscosityLinker::evaluate single called");
84 : number density;
85 : number viscosity;
86 : number yieldStress;
87 : MathMatrix<dim,dim> velocityGrad;
88 :
89 0 : (*m_spDensity)(density, globIP, time, si);
90 0 : (*m_spViscosity)(viscosity, globIP, time, si);
91 0 : (*m_spYieldStress)(yieldStress, globIP, time, si);
92 0 : (*m_spVelocityGrad)(velocityGrad, globIP, time, si);
93 :
94 : number innerSum = 0.0;
95 :
96 : // compute inner sum
97 0 : for(int d1 = 0; d1 < dim; ++d1)
98 : {
99 0 : for(int d2 = 0; d2 < dim; ++d2)
100 : {
101 0 : innerSum += pow(velocityGrad(d1,d2) + velocityGrad(d2,d1),2);
102 : }
103 : }
104 :
105 : // compute mu = eta + sigma / \sqrt( delta + 1 / I)
106 0 : value = viscosity + yieldStress/sqrt(0.5+(1.0/(pow(2, dim))*innerSum));
107 0 : value *= 1./density;
108 0 : }
109 :
110 : // function for evaluation at multiple ips??
111 : template <int refDim>
112 0 : inline void evaluate(number vValue[],
113 : const MathVector<dim> vGlobIP[],
114 : number time, int si,
115 : GridObject* elem,
116 : const MathVector<dim> vCornerCoords[],
117 : const MathVector<refDim> vLocIP[],
118 : const size_t nip,
119 : LocalVector* u,
120 : const MathMatrix<refDim, dim>* vJT = NULL) const
121 : {
122 : UG_LOG("BinghamViscosityLinker::evaluate called");
123 0 : std::vector<number> vDensity(nip);
124 0 : std::vector<number> vViscosity(nip);
125 0 : std::vector<number> vYieldStress(nip);
126 0 : std::vector<MathMatrix<dim,dim> > vVelocityGrad(nip);
127 :
128 0 : (*m_spDensity)(&vDensity[0], vGlobIP, time, si,
129 : elem, vCornerCoords, vLocIP, nip, u, vJT);
130 0 : (*m_spViscosity)(&vViscosity[0], vGlobIP, time, si,
131 : elem, vCornerCoords, vLocIP, nip, u, vJT);
132 0 : (*m_spYieldStress)(&vYieldStress[0], vGlobIP, time, si,
133 : elem, vCornerCoords, vLocIP, nip, u, vJT);
134 0 : (*m_spVelocityGrad)(&vVelocityGrad[0], vGlobIP, time, si,
135 : elem, vCornerCoords, vLocIP, nip, u, vJT);
136 :
137 0 : for(size_t ip = 0; ip < nip; ++ip)
138 : {
139 : number mu = 0.0;
140 :
141 : // compute inner sum
142 0 : for(int d1 = 0; d1 < dim; ++d1)
143 : {
144 0 : for(int d2 = 0; d2 < dim; ++d2)
145 : {
146 0 : mu += pow(vVelocityGrad[ip](d1,d2) + vVelocityGrad[ip](d2,d1),2);
147 : }
148 : }
149 :
150 : // compute mu = eta + sigma
151 0 : mu = vViscosity[ip] + vYieldStress[ip]/sqrt(0.5+(1.0/(pow(2, dim))*mu));
152 0 : vValue[ip] = mu * 1./vDensity[ip];
153 : }
154 0 : }
155 :
156 : template <int refDim>
157 0 : void eval_and_deriv(number vValue[],
158 : const MathVector<dim> vGlobIP[],
159 : number time, int si,
160 : GridObject* elem,
161 : const MathVector<dim> vCornerCoords[],
162 : const MathVector<refDim> vLocIP[],
163 : const size_t nip,
164 : LocalVector* u,
165 : bool bDeriv,
166 : int s,
167 : std::vector<std::vector<number> > vvvDeriv[],
168 : const MathMatrix<refDim, dim>* vJT = NULL) const
169 : {
170 :
171 : UG_LOG("BinghamViscosityLinker::eval_and_deriv called");
172 : // get the data of the ip series
173 0 : const number* vDensity = m_spDensity->values(s);
174 : const number* vViscosity = m_spViscosity->values(s);
175 : const number* vYieldStress = m_spYieldStress->values(s);
176 : const MathMatrix<dim,dim>* vVelocityGrad = m_spVelocityGrad->values(s);
177 :
178 0 : for(size_t ip = 0; ip < nip; ++ip)
179 : {
180 : number mu = 0.0;
181 :
182 : // compute mu =
183 0 : for(int d1 = 0; d1 < dim; ++d1)
184 : {
185 0 : for(int d2 = 0; d2 < dim; ++d2)
186 : {
187 0 : mu += pow(vVelocityGrad[ip](d1,d2) + vVelocityGrad[ip](d2,d1),2);
188 : }
189 : }
190 :
191 : // compute mu = (eta + tau_F / \sqrt(delta + TrD^2)) / rho
192 0 : mu = vViscosity[ip] + vYieldStress[ip]/sqrt(0.5+(1.0/(pow(2, dim))*mu));
193 0 : vValue[ip] = mu * 1./vDensity[ip];
194 : }
195 :
196 : // Compute the derivatives at all ips //
197 : /////////////////////////////////////////////
198 :
199 : // check if something is left to do
200 0 : if(!bDeriv || this->zero_derivative()) return;
201 :
202 : // clear all derivative values
203 0 : this->set_zero(vvvDeriv, nip);
204 :
205 : // Derivatives of Density
206 0 : if(m_spDDensity.valid() && !m_spDDensity->zero_derivative())
207 : {
208 0 : for(size_t ip = 0; ip < nip; ++ip)
209 : {
210 0 : for(size_t fct = 0; fct < m_spDDensity->num_fct(); ++fct)
211 : {
212 : // get derivative of density w.r.t. to all functions
213 : const number* vDDensity = m_spDDensity->deriv(s, ip, fct);
214 :
215 : // get common fct id for this function
216 : const size_t commonFct = this->input_common_fct(_RHO_, fct);
217 :
218 : // loop all shapes and set the derivative
219 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
220 : {
221 0 : vvvDeriv[ip][commonFct][sh] -= vDDensity[sh] / vDensity[ip] * vValue[ip];
222 : }
223 : }
224 : }
225 : }
226 :
227 : // Derivatives of Viscosity
228 0 : if(m_spDViscosity.valid() && !m_spDViscosity->zero_derivative())
229 : {
230 0 : for(size_t ip = 0; ip < nip; ++ip)
231 : {
232 0 : for(size_t fct = 0; fct < m_spDViscosity->num_fct(); ++fct)
233 : {
234 : // get derivative of viscosity w.r.t. to all functions
235 : const number* vDViscosity = m_spDViscosity->deriv(s, ip, fct);
236 :
237 : // get common fct id for this function
238 : const size_t commonFct = this->input_common_fct(_ETA_, fct);
239 :
240 : // loop all shapes and set the derivative
241 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
242 : {
243 0 : vvvDeriv[ip][commonFct][sh] += vDViscosity[sh] / vDensity[ip];
244 : }
245 : }
246 : }
247 : }
248 :
249 : // Derivatives of yield stress
250 0 : if(m_spDYieldStress.valid() && !m_spDYieldStress->zero_derivative())
251 : {
252 0 : for(size_t ip = 0; ip < nip; ++ip)
253 : {
254 : // Precompute 1 / (rho * \sqrt(delta + TrD^2))
255 0 : number rInvariant = vValue[ip] - (vViscosity[ip] / vDensity[ip]);
256 :
257 0 : for(size_t fct = 0; fct < m_spDYieldStress->num_fct(); ++fct)
258 : {
259 : // get derivative of yield stress w.r.t. to all functions
260 : const number* vDYieldStress = m_spDYieldStress->deriv(s, ip, fct);
261 :
262 : // get common fct id for this function
263 : const size_t commonFct = this->input_common_fct(_TAU_, fct);
264 :
265 : // loop all shapes and set the derivative
266 0 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
267 : {
268 0 : vvvDeriv[ip][commonFct][sh] += vDYieldStress[sh] * rInvariant;
269 : }
270 : }
271 : }
272 : }
273 :
274 : // Derivatives of velocity gradient
275 : UG_LOG("Derivatives of velocity gradient missing");
276 : /*if(m_spDVelocityGrad.valid() && !m_spDVelocityGrad->zero_derivative())
277 : {
278 : for(size_t ip = 0; ip < nip; ++ip)
279 : {
280 : for(size_t fct = 0; fct < m_spDYieldStress->num_fct(); ++fct)
281 : {
282 : // get derivative of velocity gradient w.r.t. to all functions
283 : const number* vDYieldStress = m_spDYieldStress->deriv(s, ip, fct);
284 :
285 : // get common fct id for this function
286 : const size_t commonFct = this->input_common_fct(_DV_, fct);
287 :
288 : // loop all shapes and set the derivative
289 : for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
290 : {
291 : VecScaleAppend(vvvDeriv[ip][commonFct][sh], vDYieldStress[sh], rInvariant);
292 : }
293 : }
294 : }
295 : }*/
296 : }
297 :
298 :
299 :
300 : // Setter functions for imports
301 : /// set density import
302 0 : void set_density(SmartPtr<CplUserData<number, dim> > data)
303 : {
304 0 : m_spDensity = data;
305 0 : m_spDDensity = data.template cast_dynamic<DependentUserData<number, dim> >();
306 0 : base_type::set_input(_RHO_, data, data);
307 0 : }
308 :
309 0 : void set_density(number val)
310 : {
311 0 : set_density(make_sp(new ConstUserNumber<dim>(val)));
312 0 : }
313 :
314 : /// set viscosity import
315 0 : void set_viscosity(SmartPtr<CplUserData<number, dim> > data)
316 : {
317 0 : m_spViscosity = data;
318 0 : m_spDViscosity = data.template cast_dynamic<DependentUserData<number, dim> >();
319 0 : base_type::set_input(_ETA_, data, data);
320 0 : }
321 :
322 0 : void set_viscosity(number val)
323 : {
324 0 : set_viscosity(make_sp(new ConstUserNumber<dim>(val)));
325 0 : }
326 :
327 : /// set density import
328 0 : void set_yield_stress(SmartPtr<CplUserData<number, dim> > data)
329 : {
330 0 : m_spYieldStress = data;
331 0 : m_spDYieldStress = data.template cast_dynamic<DependentUserData<number, dim> >();
332 0 : base_type::set_input(_TAU_, data, data);
333 0 : }
334 :
335 0 : void set_yield_stress(number val)
336 : {
337 0 : set_yield_stress(make_sp(new ConstUserNumber<dim>(val)));
338 0 : }
339 :
340 : /// set velocity gradient import
341 0 : void set_velocity_gradient(SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > data)
342 : {
343 : UG_LOG("Debug mark 1\n");
344 0 : m_spVelocityGrad = data;
345 : UG_LOG("Debug mark 2\n");
346 : try{
347 0 : UG_LOG(typeid(data).name() << "\n");
348 0 : m_spDVelocityGrad = data.template cast_dynamic<DependentUserData<MathMatrix<dim,dim>, dim> >();
349 0 : } catch (std::exception& e) {
350 0 : UG_LOG("Error: " << e.what());
351 : }
352 : UG_LOG("Debug mark 3\n");
353 0 : base_type::set_input(_DV_, data, data);
354 0 : }
355 :
356 : protected:
357 : // variables for storing imports
358 : /// import for density
359 : static const size_t _RHO_ = 0;
360 : SmartPtr<CplUserData<number, dim> > m_spDensity;
361 : SmartPtr<DependentUserData<number, dim> > m_spDDensity;
362 : /// import for viscosity
363 : static const size_t _ETA_ = 1;
364 : SmartPtr<CplUserData<number, dim> > m_spViscosity;
365 : SmartPtr<DependentUserData<number, dim> > m_spDViscosity;
366 : /// import for yield stress
367 : static const size_t _TAU_ = 2;
368 : SmartPtr<CplUserData<number, dim> > m_spYieldStress;
369 : SmartPtr<DependentUserData<number, dim> > m_spDYieldStress;
370 : /// import for velocity gradient
371 : static const size_t _DV_ = 3;
372 : SmartPtr<CplUserData<MathMatrix<dim,dim>, dim> > m_spVelocityGrad;
373 : SmartPtr<DependentUserData<MathMatrix<dim,dim>, dim> > m_spDVelocityGrad;
374 : };
375 :
376 : } // end of namespace ug
377 :
378 : #endif // __H__UG__LIB_DISC__SPATIAL_DISC__BINGHAM_VISCOSITY_LINKER__
|
