GRINS-0.7.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-2016 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
27 #include "grins/low_mach_navier_stokes_vms_stab.h"
28 
29 // GRINS
30 #include "grins/assembly_context.h"
31 #include "grins/constant_viscosity.h"
32 #include "grins/constant_specific_heat.h"
33 #include "grins/constant_conductivity.h"
34 
35 // libMesh
36 #include "libmesh/quadrature.h"
37 
38 namespace GRINS
39 {
40 
41  template<class Mu, class SH, class TC>
42  LowMachNavierStokesVMSStabilization<Mu,SH,TC>::LowMachNavierStokesVMSStabilization( const std::string& physics_name,
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>
50  LowMachNavierStokesVMSStabilization<Mu,SH,TC>::~LowMachNavierStokesVMSStabilization()
51  {
52  return;
53  }
54 
55  template<class Mu, class SH, class TC>
56  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::element_time_derivative( bool compute_jacobian,
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>
75  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::mass_residual( bool compute_jacobian,
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>
94  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::assemble_continuity_time_deriv( bool /*compute_jacobian*/,
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->_press_var.p()).size();
99 
100  // Element Jacobian * quadrature weights for interior integration.
101  const std::vector<libMesh::Real> &JxW =
102  context.get_element_fe(this->_flow_vars.u())->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->_press_var.p())->get_dphi();
107 
108  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // 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->_flow_vars.u());
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->_temp_vars.T(), 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->_flow_vars.u(), qp ),
126  context.interior_value( this->_flow_vars.v(), qp ) );
127  if( this->_dim == 3 )
128  U(2) = context.interior_value( this->_flow_vars.w(), 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>
147  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::assemble_momentum_time_deriv( bool /*compute_jacobian*/,
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->_flow_vars.u()).size();
152 
153  // Check number of dofs is same for _flow_vars.u(), v_var and w_var.
154  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
155  if (this->_dim == 3)
156  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
157 
158  // Element Jacobian * quadrature weights for interior integration.
159  const std::vector<libMesh::Real> &JxW =
160  context.get_element_fe(this->_flow_vars.u())->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->_flow_vars.u())->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->_flow_vars.u())->get_dphi();
169 
170  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{u}
171  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{v}
172  libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
173 
174  if( this->_dim == 3 )
175  {
176  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
177  }
178 
179  unsigned int n_qpoints = context.get_element_qrule().n_points();
180 
181  for (unsigned int qp=0; qp != n_qpoints; qp++)
182  {
183  libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
184  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
185 
186  libMesh::RealGradient U( context.interior_value(this->_flow_vars.u(), qp),
187  context.interior_value(this->_flow_vars.v(), qp) );
188 
189  libMesh::RealGradient grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
190  libMesh::RealGradient grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
191  libMesh::RealGradient grad_w;
192 
193  if( this->_dim == 3 )
194  {
195  U(2) = context.interior_value(this->_flow_vars.w(), qp);
196  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
197  }
198 
199  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
200 
201  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
202  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
203  libMesh::Real mu = this->_mu(T);
204 
205  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
206  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
207 
208  libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
209  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
210 
211  for (unsigned int i=0; i != n_u_dofs; i++)
212  {
213  Fu(i) += ( -tau_C*RC_s*u_gradphi[i][qp](0)
214  - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp]
215  + rho*tau_M*RM_s*grad_u*u_phi[i][qp]
216  + tau_M*RM_s(0)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
217 
218  Fv(i) += ( -tau_C*RC_s*u_gradphi[i][qp](1)
219  - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp]
220  + rho*tau_M*RM_s*grad_v*u_phi[i][qp]
221  + tau_M*RM_s(1)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
222 
223  if( this->_dim == 3 )
224  {
225  (*Fw)(i) += ( -tau_C*RC_s*u_gradphi[i][qp](2)
226  - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp]
227  + rho*tau_M*RM_s*grad_w*u_phi[i][qp]
228  + tau_M*RM_s(2)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
229  }
230  }
231 
232  }
233  return;
234  }
235 
236  template<class Mu, class SH, class TC>
237  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::assemble_energy_time_deriv( bool /*compute_jacobian*/,
238  AssemblyContext& context )
239  {
240  // The number of local degrees of freedom in each variable.
241  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
242 
243  // Element Jacobian * quadrature weights for interior integration.
244  const std::vector<libMesh::Real> &JxW =
245  context.get_element_fe(this->_temp_vars.T())->get_JxW();
246 
247  // The temperature shape functions at interior quadrature points.
248  const std::vector<std::vector<libMesh::Real> >& T_phi =
249  context.get_element_fe(this->_temp_vars.T())->get_phi();
250 
251  // The temperature shape functions gradients at interior quadrature points.
252  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
253  context.get_element_fe(this->_temp_vars.T())->get_dphi();
254 
255  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T()); // R_{T}
256 
257  unsigned int n_qpoints = context.get_element_qrule().n_points();
258 
259  for (unsigned int qp=0; qp != n_qpoints; qp++)
260  {
261  libMesh::Number u, v;
262  u = context.interior_value(this->_flow_vars.u(), qp);
263  v = context.interior_value(this->_flow_vars.v(), qp);
264 
265  libMesh::Gradient grad_T = context.interior_gradient(this->_temp_vars.T(), qp);
266 
267  libMesh::NumberVectorValue U(u,v);
268  if (this->_dim == 3)
269  U(2) = context.interior_value(this->_flow_vars.w(), qp); // w
270 
271  libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
272  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
273 
274  libMesh::Real mu = this->_mu(T);
275  libMesh::Real k = this->_k(T);
276  libMesh::Real cp = this->_cp(T);
277 
278  libMesh::Number rho_cp = rho*this->_cp(T);
279 
280  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
281 
282  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
283  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
284 
285  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
286  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
287 
288  libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
289  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
290 
291  for (unsigned int i=0; i != n_T_dofs; i++)
292  {
293  FT(i) += ( rho_cp*tau_M*RM_s*grad_T*T_phi[i][qp]
294  - rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]
295  + rho_cp*tau_E*RE_s*tau_M*RM_s*T_gradphi[i][qp] )*JxW[qp];
296  }
297 
298  }
299 
300  return;
301  }
302 
303  template<class Mu, class SH, class TC>
304  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::assemble_continuity_mass_residual( bool /*compute_jacobian*/,
305  AssemblyContext& context )
306  {
307  // The number of local degrees of freedom in each variable.
308  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
309 
310  // Element Jacobian * quadrature weights for interior integration.
311  const std::vector<libMesh::Real> &JxW =
312  context.get_element_fe(this->_flow_vars.u())->get_JxW();
313 
314  // The pressure shape functions at interior quadrature points.
315  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
316  context.get_element_fe(this->_press_var.p())->get_dphi();
317 
318  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
319 
320  unsigned int n_qpoints = context.get_element_qrule().n_points();
321 
322  for (unsigned int qp=0; qp != n_qpoints; qp++)
323  {
324  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
325 
326  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
327  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
328 
329  libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
330  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
331 
332  libMesh::Real mu = this->_mu(T);
333 
334  libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
335  context.fixed_interior_value( this->_flow_vars.v(), qp ) );
336  if( this->_dim == 3 )
337  U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
338 
339  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
340  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
341 
342  // Now a loop over the pressure degrees of freedom. This
343  // computes the contributions of the continuity equation.
344  for (unsigned int i=0; i != n_p_dofs; i++)
345  {
346  Fp(i) += tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
347  }
348  }
349 
350  return;
351  }
352 
353  template<class Mu, class SH, class TC>
354  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::assemble_momentum_mass_residual( bool /*compute_jacobian*/,
355  AssemblyContext& context )
356  {
357  // The number of local degrees of freedom in each variable.
358  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
359 
360  // Check number of dofs is same for _flow_vars.u(), v_var and w_var.
361  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
362  if (this->_dim == 3)
363  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
364 
365  // Element Jacobian * quadrature weights for interior integration.
366  const std::vector<libMesh::Real> &JxW =
367  context.get_element_fe(this->_flow_vars.u())->get_JxW();
368 
369  // The pressure shape functions at interior quadrature points.
370  const std::vector<std::vector<libMesh::Real> >& u_phi =
371  context.get_element_fe(this->_flow_vars.u())->get_phi();
372 
373  // The velocity shape function gradients at interior quadrature points.
374  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
375  context.get_element_fe(this->_flow_vars.u())->get_dphi();
376 
377  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{u}
378  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{v}
379  libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
380 
381  if( this->_dim == 3 )
382  {
383  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
384  }
385 
386  unsigned int n_qpoints = context.get_element_qrule().n_points();
387  for (unsigned int qp=0; qp != n_qpoints; qp++)
388  {
389  libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
390  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
391 
392  libMesh::Real mu = this->_mu(T);
393 
394  libMesh::RealGradient U( context.fixed_interior_value(this->_flow_vars.u(), qp),
395  context.fixed_interior_value(this->_flow_vars.v(), qp) );
396 
397  libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_flow_vars.u(), qp);
398  libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_flow_vars.v(), qp);
399  libMesh::RealGradient grad_w;
400 
401  if( this->_dim == 3 )
402  {
403  U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
404  grad_w = context.fixed_interior_gradient(this->_flow_vars.w(), qp);
405  }
406 
407  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
408 
409  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
410  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
411 
412  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
413  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
414 
415  libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
416  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
417  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
418 
419  for (unsigned int i=0; i != n_u_dofs; i++)
420  {
421  Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
422  + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp]
423  - rho*tau_M*RM_t*grad_u*u_phi[i][qp]
424  - tau_M*(RM_s(0)+RM_t(0))*rho*tau_M*RM_t*u_gradphi[i][qp]
425  - tau_M*RM_t(0)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
426 
427  Fv(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
428  - rho*tau_M*RM_t*grad_v*u_phi[i][qp]
429  + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp]
430  - tau_M*(RM_s(1)+RM_t(1))*rho*tau_M*RM_t*u_gradphi[i][qp]
431  - tau_M*RM_t(1)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
432 
433  if( this->_dim == 3 )
434  {
435  (*Fw)(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
436  - rho*tau_M*RM_t*grad_w*u_phi[i][qp]
437  + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp]
438  - tau_M*(RM_s(2)+RM_t(2))*rho*tau_M*RM_t*u_gradphi[i][qp]
439  - tau_M*RM_t(2)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
440  }
441  }
442 
443  }
444  return;
445  }
446 
447  template<class Mu, class SH, class TC>
448  void LowMachNavierStokesVMSStabilization<Mu,SH,TC>::assemble_energy_mass_residual( bool /*compute_jacobian*/,
449  AssemblyContext& context )
450  {
451  // The number of local degrees of freedom in each variable.
452  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
453 
454  // Element Jacobian * quadrature weights for interior integration.
455  const std::vector<libMesh::Real> &JxW =
456  context.get_element_fe(this->_temp_vars.T())->get_JxW();
457 
458  // The temperature shape functions at interior quadrature points.
459  const std::vector<std::vector<libMesh::Real> >& T_phi =
460  context.get_element_fe(this->_temp_vars.T())->get_phi();
461 
462  // The temperature shape functions gradients at interior quadrature points.
463  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
464  context.get_element_fe(this->_temp_vars.T())->get_dphi();
465 
466  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T()); // R_{T}
467 
468  unsigned int n_qpoints = context.get_element_qrule().n_points();
469 
470  for (unsigned int qp=0; qp != n_qpoints; qp++)
471  {
472  libMesh::Number u, v;
473  u = context.fixed_interior_value(this->_flow_vars.u(), qp);
474  v = context.fixed_interior_value(this->_flow_vars.v(), qp);
475 
476  libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_temp_vars.T(), qp);
477 
478  libMesh::NumberVectorValue U(u,v);
479  if (this->_dim == 3)
480  U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp); // w
481 
482  libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
483  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
484 
485  libMesh::Real mu = this->_mu(T);
486  libMesh::Real k = this->_k(T);
487  libMesh::Real cp = this->_cp(T);
488 
489  libMesh::Number rho_cp = rho*this->_cp(T);
490 
491  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
492 
493  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
494  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
495 
496  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
497  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, false );
498 
499  libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
500  libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
501 
502  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
503  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
504 
505  for (unsigned int i=0; i != n_T_dofs; i++)
506  {
507  FT(i) -= ( -rho_cp*tau_M*RM_t*grad_T*T_phi[i][qp]
508  +rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]
509  - rho_cp*tau_E*(RE_s+RE_t)*tau_M*RM_t*T_gradphi[i][qp]
510  - rho_cp*tau_E*RE_t*tau_M*RM_s*T_gradphi[i][qp] )*JxW[qp];
511  }
512 
513  }
514 
515  return;
516  }
517 
518 } // namespace GRINS
519 
520 // Instantiate
521 template class GRINS::LowMachNavierStokesVMSStabilization<GRINS::ConstantViscosity,GRINS::ConstantSpecificHeat,GRINS::ConstantConductivity>;
constant_conductivity.h
assembly_context.h
GRINS::LowMachNavierStokesVMSStabilization
Adds VMS-based stabilization to LowMachNavierStokes physics class.
Definition: low_mach_navier_stokes_vms_stab.h:36
constant_specific_heat.h
GRINS::LowMachNavierStokesVMSStabilization::mass_residual
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...
Definition: low_mach_navier_stokes_vms_stab.C:75
GRINS
GRINS namespace.
Definition: arrhenius_catalycity.h:31
GRINS::LowMachNavierStokesVMSStabilization::element_time_derivative
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
Definition: low_mach_navier_stokes_vms_stab.C:56
GRINS::LowMachNavierStokesVMSStabilization::assemble_energy_time_deriv
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_vms_stab.C:237
GRINS::LowMachNavierStokesVMSStabilization::assemble_momentum_mass_residual
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_vms_stab.C:354
GRINS::LowMachNavierStokesStabilizationBase
Adds VMS-based stabilization to LowMachNavierStokes physics class.
Definition: low_mach_navier_stokes_stab_base.h:37
GRINS::LowMachNavierStokesVMSStabilization::assemble_momentum_time_deriv
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_vms_stab.C:147
GRINS::LowMachNavierStokesVMSStabilization::assemble_continuity_mass_residual
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_vms_stab.C:304
GRINS::LowMachNavierStokesVMSStabilization::assemble_continuity_time_deriv
void assemble_continuity_time_deriv(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_vms_stab.C:94
low_mach_navier_stokes_vms_stab.h
constant_viscosity.h
GRINS::LowMachNavierStokesVMSStabilization::assemble_energy_mass_residual
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_vms_stab.C:448
Generated on Thu Jun 2 2016 21:52:28 for GRINS-0.7.0 by  doxygen 1.8.10