GRINS-0.6.0
low_mach_navier_stokes_vms_stab.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // GRINS - General Reacting Incompressible Navier-Stokes
5 //
6 // Copyright (C) 2014-2015 Paul T. Bauman, Roy H. Stogner
7 // Copyright (C) 2010-2013 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 
26 // This class
28 
29 // GRINS
30 #include "grins/assembly_context.h"
34 
35 // libMesh
36 #include "libmesh/quadrature.h"
37 
38 namespace GRINS
39 {
40 
41  template<class Mu, class SH, class TC>
43  const GetPot& input )
44  : LowMachNavierStokesStabilizationBase<Mu,SH,TC>(physics_name,input)
45  {
46  return;
47  }
48 
49  template<class Mu, class SH, class TC>
51  {
52  return;
53  }
54 
55  template<class Mu, class SH, class TC>
57  AssemblyContext& context,
58  CachedValues& /*cache*/ )
59  {
60 #ifdef GRINS_USE_GRVY_TIMERS
61  this->_timer->BeginTimer("LowMachNavierStokesVMSStabilization::element_time_derivative");
62 #endif
63 
64  this->assemble_continuity_time_deriv( compute_jacobian, context );
65  this->assemble_momentum_time_deriv( compute_jacobian, context );
66  this->assemble_energy_time_deriv( compute_jacobian, context );
67 
68 #ifdef GRINS_USE_GRVY_TIMERS
69  this->_timer->EndTimer("LowMachNavierStokesVMSStabilization::element_time_derivative");
70 #endif
71  return;
72  }
73 
74  template<class Mu, class SH, class TC>
76  AssemblyContext& context,
77  CachedValues& /*cache*/ )
78  {
79 #ifdef GRINS_USE_GRVY_TIMERS
80  this->_timer->BeginTimer("LowMachNavierStokesVMSStabilization::mass_residual");
81 #endif
82 
83  this->assemble_continuity_mass_residual( compute_jacobian, context );
84  this->assemble_momentum_mass_residual( compute_jacobian, context );
85  this->assemble_energy_mass_residual( compute_jacobian, context );
86 
87 #ifdef GRINS_USE_GRVY_TIMERS
88  this->_timer->EndTimer("LowMachNavierStokesVMSStabilization::mass_residual");
89 #endif
90  return;
91  }
92 
93  template<class Mu, class SH, class TC>
95  AssemblyContext& context )
96  {
97  // The number of local degrees of freedom in each variable.
98  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
99 
100  // Element Jacobian * quadrature weights for interior integration.
101  const std::vector<libMesh::Real> &JxW =
102  context.get_element_fe(this->_u_var)->get_JxW();
103 
104  // The pressure shape functions at interior quadrature points.
105  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
106  context.get_element_fe(this->_p_var)->get_dphi();
107 
108  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var); // R_{p}
109 
110  unsigned int n_qpoints = context.get_element_qrule().n_points();
111 
112  for (unsigned int qp=0; qp != n_qpoints; qp++)
113  {
114  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
115 
116  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
117  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
118 
119  libMesh::Real T = context.interior_value( this->_T_var, qp );
120 
121  libMesh::Real mu = this->_mu(T);
122 
123  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
124 
125  libMesh::RealGradient U( context.interior_value( this->_u_var, qp ),
126  context.interior_value( this->_v_var, qp ) );
127  if( this->_dim == 3 )
128  U(2) = context.interior_value( this->_w_var, qp );
129 
130  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
131 
132  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
133 
134  // Now a loop over the pressure degrees of freedom. This
135  // computes the contributions of the continuity equation.
136  for (unsigned int i=0; i != n_p_dofs; i++)
137  {
138  Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
139  }
140 
141  }
142 
143  return;
144  }
145 
146  template<class Mu, class SH, class TC>
148  AssemblyContext& context )
149  {
150  // The number of local degrees of freedom in each variable.
151  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
152 
153  // Check number of dofs is same for _u_var, v_var and w_var.
154  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
155  if (this->_dim == 3)
156  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
157 
158  // Element Jacobian * quadrature weights for interior integration.
159  const std::vector<libMesh::Real> &JxW =
160  context.get_element_fe(this->_u_var)->get_JxW();
161 
162  // The pressure shape functions at interior quadrature points.
163  const std::vector<std::vector<libMesh::Real> >& u_phi =
164  context.get_element_fe(this->_u_var)->get_phi();
165 
166  // The velocity shape function gradients at interior quadrature points.
167  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
168  context.get_element_fe(this->_u_var)->get_dphi();
169 
170  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var); // R_{u}
171  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var); // R_{v}
172  libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var); // R_{w}
173 
174  unsigned int n_qpoints = context.get_element_qrule().n_points();
175 
176  for (unsigned int qp=0; qp != n_qpoints; qp++)
177  {
178  libMesh::Real T = context.interior_value( this->_T_var, qp );
179  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
180 
181  libMesh::RealGradient U( context.interior_value(this->_u_var, qp),
182  context.interior_value(this->_v_var, qp) );
183 
184  libMesh::RealGradient grad_u = context.interior_gradient(this->_u_var, qp);
185  libMesh::RealGradient grad_v = context.interior_gradient(this->_v_var, qp);
186  libMesh::RealGradient grad_w;
187 
188  if( this->_dim == 3 )
189  {
190  U(2) = context.interior_value(this->_w_var, qp);
191  grad_w = context.interior_gradient(this->_w_var, qp);
192  }
193 
194  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
195 
196  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
197  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
198  libMesh::Real mu = this->_mu(T);
199 
200  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
201  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
202 
203  libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
204  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
205 
206  for (unsigned int i=0; i != n_u_dofs; i++)
207  {
208  Fu(i) += ( -tau_C*RC_s*u_gradphi[i][qp](0)
209  - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp]
210  + rho*tau_M*RM_s*grad_u*u_phi[i][qp]
211  + tau_M*RM_s(0)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
212 
213  Fv(i) += ( -tau_C*RC_s*u_gradphi[i][qp](1)
214  - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp]
215  + rho*tau_M*RM_s*grad_v*u_phi[i][qp]
216  + tau_M*RM_s(1)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
217 
218  if( this->_dim == 3 )
219  {
220  Fw(i) += ( -tau_C*RC_s*u_gradphi[i][qp](2)
221  - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp]
222  + rho*tau_M*RM_s*grad_w*u_phi[i][qp]
223  + tau_M*RM_s(2)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
224  }
225  }
226 
227  }
228  return;
229  }
230 
231  template<class Mu, class SH, class TC>
233  AssemblyContext& context )
234  {
235  // The number of local degrees of freedom in each variable.
236  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
237 
238  // Element Jacobian * quadrature weights for interior integration.
239  const std::vector<libMesh::Real> &JxW =
240  context.get_element_fe(this->_T_var)->get_JxW();
241 
242  // The temperature shape functions at interior quadrature points.
243  const std::vector<std::vector<libMesh::Real> >& T_phi =
244  context.get_element_fe(this->_T_var)->get_phi();
245 
246  // The temperature shape functions gradients at interior quadrature points.
247  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
248  context.get_element_fe(this->_T_var)->get_dphi();
249 
250  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var); // R_{T}
251 
252  unsigned int n_qpoints = context.get_element_qrule().n_points();
253 
254  for (unsigned int qp=0; qp != n_qpoints; qp++)
255  {
256  libMesh::Number u, v;
257  u = context.interior_value(this->_u_var, qp);
258  v = context.interior_value(this->_v_var, qp);
259 
260  libMesh::Gradient grad_T = context.interior_gradient(this->_T_var, qp);
261 
262  libMesh::NumberVectorValue U(u,v);
263  if (this->_dim == 3)
264  U(2) = context.interior_value(this->_w_var, qp); // w
265 
266  libMesh::Real T = context.interior_value( this->_T_var, qp );
267  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
268 
269  libMesh::Real mu = this->_mu(T);
270  libMesh::Real k = this->_k(T);
271  libMesh::Real cp = this->_cp(T);
272 
273  libMesh::Number rho_cp = rho*this->_cp(T);
274 
275  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
276 
277  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
278  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
279 
280  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
281  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
282 
283  libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
284  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
285 
286  for (unsigned int i=0; i != n_T_dofs; i++)
287  {
288  FT(i) += ( rho_cp*tau_M*RM_s*grad_T*T_phi[i][qp]
289  - rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]
290  + rho_cp*tau_E*RE_s*tau_M*RM_s*T_gradphi[i][qp] )*JxW[qp];
291  }
292 
293  }
294 
295  return;
296  }
297 
298  template<class Mu, class SH, class TC>
300  AssemblyContext& context )
301  {
302  // The number of local degrees of freedom in each variable.
303  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
304 
305  // Element Jacobian * quadrature weights for interior integration.
306  const std::vector<libMesh::Real> &JxW =
307  context.get_element_fe(this->_u_var)->get_JxW();
308 
309  // The pressure shape functions at interior quadrature points.
310  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
311  context.get_element_fe(this->_p_var)->get_dphi();
312 
313  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var); // R_{p}
314 
315  unsigned int n_qpoints = context.get_element_qrule().n_points();
316 
317  for (unsigned int qp=0; qp != n_qpoints; qp++)
318  {
319  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
320 
321  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
322  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
323 
324  libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
325  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
326 
327  libMesh::Real mu = this->_mu(T);
328 
329  libMesh::RealGradient U( context.fixed_interior_value( this->_u_var, qp ),
330  context.fixed_interior_value( this->_v_var, qp ) );
331  if( this->_dim == 3 )
332  U(2) = context.fixed_interior_value( this->_w_var, qp );
333 
334  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
335  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
336 
337  // Now a loop over the pressure degrees of freedom. This
338  // computes the contributions of the continuity equation.
339  for (unsigned int i=0; i != n_p_dofs; i++)
340  {
341  Fp(i) += tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
342  }
343  }
344 
345  return;
346  }
347 
348  template<class Mu, class SH, class TC>
350  AssemblyContext& context )
351  {
352  // The number of local degrees of freedom in each variable.
353  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
354 
355  // Check number of dofs is same for _u_var, v_var and w_var.
356  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
357  if (this->_dim == 3)
358  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
359 
360  // Element Jacobian * quadrature weights for interior integration.
361  const std::vector<libMesh::Real> &JxW =
362  context.get_element_fe(this->_u_var)->get_JxW();
363 
364  // The pressure shape functions at interior quadrature points.
365  const std::vector<std::vector<libMesh::Real> >& u_phi =
366  context.get_element_fe(this->_u_var)->get_phi();
367 
368  // The velocity shape function gradients at interior quadrature points.
369  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
370  context.get_element_fe(this->_u_var)->get_dphi();
371 
372  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var); // R_{u}
373  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var); // R_{v}
374  libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var); // R_{w}
375 
376  unsigned int n_qpoints = context.get_element_qrule().n_points();
377  for (unsigned int qp=0; qp != n_qpoints; qp++)
378  {
379  libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
380  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
381 
382  libMesh::Real mu = this->_mu(T);
383 
384  libMesh::RealGradient U( context.fixed_interior_value(this->_u_var, qp),
385  context.fixed_interior_value(this->_v_var, qp) );
386 
387  libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_u_var, qp);
388  libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_v_var, qp);
389  libMesh::RealGradient grad_w;
390 
391  if( this->_dim == 3 )
392  {
393  U(2) = context.fixed_interior_value(this->_w_var, qp);
394  grad_w = context.fixed_interior_gradient(this->_w_var, qp);
395  }
396 
397  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
398 
399  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
400  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
401 
402  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
403  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
404 
405  libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
406  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
407  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
408 
409  for (unsigned int i=0; i != n_u_dofs; i++)
410  {
411  Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
412  + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp]
413  - rho*tau_M*RM_t*grad_u*u_phi[i][qp]
414  - tau_M*(RM_s(0)+RM_t(0))*rho*tau_M*RM_t*u_gradphi[i][qp]
415  - tau_M*RM_t(0)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
416 
417  Fv(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
418  - rho*tau_M*RM_t*grad_v*u_phi[i][qp]
419  + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp]
420  - tau_M*(RM_s(1)+RM_t(1))*rho*tau_M*RM_t*u_gradphi[i][qp]
421  - tau_M*RM_t(1)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
422 
423  if( this->_dim == 3 )
424  {
425  Fw(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
426  - rho*tau_M*RM_t*grad_w*u_phi[i][qp]
427  + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp]
428  - tau_M*(RM_s(2)+RM_t(2))*rho*tau_M*RM_t*u_gradphi[i][qp]
429  - tau_M*RM_t(2)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
430  }
431  }
432 
433  }
434  return;
435  }
436 
437  template<class Mu, class SH, class TC>
439  AssemblyContext& context )
440  {
441  // The number of local degrees of freedom in each variable.
442  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
443 
444  // Element Jacobian * quadrature weights for interior integration.
445  const std::vector<libMesh::Real> &JxW =
446  context.get_element_fe(this->_T_var)->get_JxW();
447 
448  // The temperature shape functions at interior quadrature points.
449  const std::vector<std::vector<libMesh::Real> >& T_phi =
450  context.get_element_fe(this->_T_var)->get_phi();
451 
452  // The temperature shape functions gradients at interior quadrature points.
453  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
454  context.get_element_fe(this->_T_var)->get_dphi();
455 
456  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var); // R_{T}
457 
458  unsigned int n_qpoints = context.get_element_qrule().n_points();
459 
460  for (unsigned int qp=0; qp != n_qpoints; qp++)
461  {
462  libMesh::Number u, v;
463  u = context.fixed_interior_value(this->_u_var, qp);
464  v = context.fixed_interior_value(this->_v_var, qp);
465 
466  libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_T_var, qp);
467 
468  libMesh::NumberVectorValue U(u,v);
469  if (this->_dim == 3)
470  U(2) = context.fixed_interior_value(this->_w_var, qp); // w
471 
472  libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
473  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
474 
475  libMesh::Real mu = this->_mu(T);
476  libMesh::Real k = this->_k(T);
477  libMesh::Real cp = this->_cp(T);
478 
479  libMesh::Number rho_cp = rho*this->_cp(T);
480 
481  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
482 
483  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
484  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
485 
486  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
487  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, false );
488 
489  libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
490  libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
491 
492  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
493  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
494 
495  for (unsigned int i=0; i != n_T_dofs; i++)
496  {
497  FT(i) -= ( -rho_cp*tau_M*RM_t*grad_T*T_phi[i][qp]
498  +rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]
499  - rho_cp*tau_E*(RE_s+RE_t)*tau_M*RM_t*T_gradphi[i][qp]
500  - rho_cp*tau_E*RE_t*tau_M*RM_s*T_gradphi[i][qp] )*JxW[qp];
501  }
502 
503  }
504 
505  return;
506  }
507 
508 } // namespace GRINS
509 
510 // Instantiate
Adds VMS-based stabilization to LowMachNavierStokes physics class.
virtual void mass_residual(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part...
GRINS namespace.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context)
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
Adds VMS-based stabilization to LowMachNavierStokes physics class.
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context)
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
void assemble_continuity_time_deriv(bool compute_jacobian, AssemblyContext &context)
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &context)

Generated on Mon Jun 22 2015 21:32:20 for GRINS-0.6.0 by  doxygen 1.8.9.1