GRINS-0.6.0
low_mach_navier_stokes_spgsm_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
27 #include "grins/low_mach_navier_stokes_spgsm_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  template<class Mu, class SH, class TC>
41  LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::LowMachNavierStokesSPGSMStabilization( const std::string& physics_name,
42  const GetPot& input )
43  : LowMachNavierStokesStabilizationBase<Mu,SH,TC>(physics_name,input)
44  {
45  }
46 
47  template<class Mu, class SH, class TC>
48  LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::~LowMachNavierStokesSPGSMStabilization()
49  {
50  return;
51  }
52 
53  template<class Mu, class SH, class TC>
54  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::element_time_derivative( bool compute_jacobian,
55  AssemblyContext& context,
56  CachedValues& /*cache*/ )
57  {
58 #ifdef GRINS_USE_GRVY_TIMERS
59  this->_timer->BeginTimer("LowMachNavierStokesSPGSMStabilization::element_time_derivative");
60 #endif
61 
62  this->assemble_continuity_time_deriv( compute_jacobian, context );
63  this->assemble_momentum_time_deriv( compute_jacobian, context );
64  this->assemble_energy_time_deriv( compute_jacobian, context );
65 
66 #ifdef GRINS_USE_GRVY_TIMERS
67  this->_timer->EndTimer("LowMachNavierStokesSPGSMStabilization::element_time_derivative");
68 #endif
69  return;
70  }
71 
72  template<class Mu, class SH, class TC>
73  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::mass_residual( bool compute_jacobian,
74  AssemblyContext& context,
75  CachedValues& /*cache*/ )
76  {
77 #ifdef GRINS_USE_GRVY_TIMERS
78  this->_timer->BeginTimer("LowMachNavierStokesSPGSMStabilization::mass_residual");
79 #endif
80 
81  this->assemble_continuity_mass_residual( compute_jacobian, context );
82  this->assemble_momentum_mass_residual( compute_jacobian, context );
83  this->assemble_energy_mass_residual( compute_jacobian, context );
84 
85 #ifdef GRINS_USE_GRVY_TIMERS
86  this->_timer->EndTimer("LowMachNavierStokesSPGSMStabilization::mass_residual");
87 #endif
88  return;
89  }
90 
91  template<class Mu, class SH, class TC>
92  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::assemble_continuity_time_deriv( bool /*compute_jacobian*/,
93  AssemblyContext& context )
94  {
95  // The number of local degrees of freedom in each variable.
96  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
97 
98  // Element Jacobian * quadrature weights for interior integration.
99  const std::vector<libMesh::Real> &JxW =
100  context.get_element_fe(this->_u_var)->get_JxW();
101 
102  // The pressure shape functions at interior quadrature points.
103  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
104  context.get_element_fe(this->_p_var)->get_dphi();
105 
106  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var); // R_{p}
107 
108  unsigned int n_qpoints = context.get_element_qrule().n_points();
109 
110  for (unsigned int qp=0; qp != n_qpoints; qp++)
111  {
112  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
113 
114  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
115  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
116 
117  libMesh::Real T = context.interior_value( this->_T_var, qp );
118 
119  libMesh::Real mu = this->_mu(T);
120 
121  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
122 
123  libMesh::RealGradient U( context.interior_value( this->_u_var, qp ),
124  context.interior_value( this->_v_var, qp ) );
125  if( this->_dim == 3 )
126  U(2) = context.interior_value( this->_w_var, qp );
127 
128  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
129 
130  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
131 
132  // Now a loop over the pressure degrees of freedom. This
133  // computes the contributions of the continuity equation.
134  for (unsigned int i=0; i != n_p_dofs; i++)
135  {
136  Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
137  }
138 
139  }
140 
141  return;
142  }
143 
144  template<class Mu, class SH, class TC>
145  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::assemble_momentum_time_deriv( bool /*compute_jacobian*/,
146  AssemblyContext& context )
147  {
148  // The number of local degrees of freedom in each variable.
149  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
150 
151  // Check number of dofs is same for _u_var, v_var and w_var.
152  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
153  if (this->_dim == 3)
154  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
155 
156  // Element Jacobian * quadrature weights for interior integration.
157  const std::vector<libMesh::Real> &JxW =
158  context.get_element_fe(this->_u_var)->get_JxW();
159 
160  // The velocity shape function gradients at interior quadrature points.
161  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
162  context.get_element_fe(this->_u_var)->get_dphi();
163 
164  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var); // R_{u}
165  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var); // R_{v}
166  libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var); // R_{w}
167 
168  unsigned int n_qpoints = context.get_element_qrule().n_points();
169 
170  for (unsigned int qp=0; qp != n_qpoints; qp++)
171  {
172  libMesh::Real T = context.interior_value( this->_T_var, qp );
173  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
174 
175  libMesh::RealGradient U( context.interior_value(this->_u_var, qp),
176  context.interior_value(this->_v_var, qp) );
177 
178  libMesh::RealGradient grad_u = context.interior_gradient(this->_u_var, qp);
179  libMesh::RealGradient grad_v = context.interior_gradient(this->_v_var, qp);
180  libMesh::RealGradient grad_w;
181 
182  if( this->_dim == 3 )
183  {
184  U(2) = context.interior_value(this->_w_var, qp);
185  grad_w = context.interior_gradient(this->_w_var, qp);
186  }
187 
188  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
189 
190  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
191  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
192  libMesh::Real mu = this->_mu(T);
193 
194  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
195  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
196 
197  libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
198  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
199 
200  for (unsigned int i=0; i != n_u_dofs; i++)
201  {
202  Fu(i) += ( - tau_C*RC_s*u_gradphi[i][qp](0)
203  - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp] )*JxW[qp];
204 
205  Fv(i) += ( - tau_C*RC_s*u_gradphi[i][qp](1)
206  - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
207 
208  if( this->_dim == 3 )
209  {
210  Fw(i) += ( - tau_C*RC_s*u_gradphi[i][qp](2)
211  - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
212  }
213  }
214 
215  }
216  return;
217  }
218 
219  template<class Mu, class SH, class TC>
220  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::assemble_energy_time_deriv( bool /*compute_jacobian*/,
221  AssemblyContext& context )
222  {
223  // The number of local degrees of freedom in each variable.
224  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
225 
226  // Element Jacobian * quadrature weights for interior integration.
227  const std::vector<libMesh::Real> &JxW =
228  context.get_element_fe(this->_T_var)->get_JxW();
229 
230  // The temperature shape functions gradients at interior quadrature points.
231  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
232  context.get_element_fe(this->_T_var)->get_dphi();
233 
234  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var); // R_{T}
235 
236  unsigned int n_qpoints = context.get_element_qrule().n_points();
237 
238  for (unsigned int qp=0; qp != n_qpoints; qp++)
239  {
240  libMesh::Number u, v;
241  u = context.interior_value(this->_u_var, qp);
242  v = context.interior_value(this->_v_var, qp);
243 
244  libMesh::Gradient grad_T = context.interior_gradient(this->_T_var, qp);
245 
246  libMesh::NumberVectorValue U(u,v);
247  if (this->_dim == 3)
248  U(2) = context.interior_value(this->_w_var, qp); // w
249 
250  libMesh::Real T = context.interior_value( this->_T_var, qp );
251  libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
252 
253  libMesh::Real k = this->_k(T);
254  libMesh::Real cp = this->_cp(T);
255 
256  libMesh::Number rho_cp = rho*this->_cp(T);
257 
258  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
259 
260  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
261  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
262 
263  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
264 
265  libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
266 
267  for (unsigned int i=0; i != n_T_dofs; i++)
268  {
269  FT(i) -= rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]*JxW[qp];
270  }
271 
272  }
273 
274  return;
275  }
276 
277  template<class Mu, class SH, class TC>
278  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::assemble_continuity_mass_residual( bool /*compute_jacobian*/,
279  AssemblyContext& context)
280  {
281  // The number of local degrees of freedom in each variable.
282  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
283 
284  // Element Jacobian * quadrature weights for interior integration.
285  const std::vector<libMesh::Real> &JxW =
286  context.get_element_fe(this->_u_var)->get_JxW();
287 
288  // The pressure shape functions at interior quadrature points.
289  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
290  context.get_element_fe(this->_p_var)->get_dphi();
291 
292  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var); // R_{p}
293 
294  unsigned int n_qpoints = context.get_element_qrule().n_points();
295 
296  for (unsigned int qp=0; qp != n_qpoints; qp++)
297  {
298  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
299 
300  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
301  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
302 
303  libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
304  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
305 
306  libMesh::Real mu = this->_mu(T);
307 
308  libMesh::RealGradient U( context.fixed_interior_value( this->_u_var, qp ),
309  context.fixed_interior_value( this->_v_var, qp ) );
310  if( this->_dim == 3 )
311  U(2) = context.fixed_interior_value( this->_w_var, qp );
312 
313  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
314  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
315 
316  // Now a loop over the pressure degrees of freedom. This
317  // computes the contributions of the continuity equation.
318  for (unsigned int i=0; i != n_p_dofs; i++)
319  {
320  Fp(i) += tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
321  }
322  }
323 
324  return;
325  }
326 
327  template<class Mu, class SH, class TC>
328  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::assemble_momentum_mass_residual( bool /*compute_jacobian*/,
329  AssemblyContext& context )
330  {
331  // The number of local degrees of freedom in each variable.
332  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
333 
334  // Check number of dofs is same for _u_var, v_var and w_var.
335  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
336  if (this->_dim == 3)
337  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
338 
339  // Element Jacobian * quadrature weights for interior integration.
340  const std::vector<libMesh::Real> &JxW =
341  context.get_element_fe(this->_u_var)->get_JxW();
342 
343  // The velocity shape function gradients at interior quadrature points.
344  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
345  context.get_element_fe(this->_u_var)->get_dphi();
346 
347  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var); // R_{u}
348  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var); // R_{v}
349  libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var); // R_{w}
350 
351  unsigned int n_qpoints = context.get_element_qrule().n_points();
352  for (unsigned int qp=0; qp != n_qpoints; qp++)
353  {
354  libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
355  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
356 
357  libMesh::Real mu = this->_mu(T);
358 
359  libMesh::RealGradient U( context.fixed_interior_value(this->_u_var, qp),
360  context.fixed_interior_value(this->_v_var, qp) );
361 
362  libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_u_var, qp);
363  libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_v_var, qp);
364  libMesh::RealGradient grad_w;
365 
366  if( this->_dim == 3 )
367  {
368  U(2) = context.fixed_interior_value(this->_w_var, qp);
369  grad_w = context.fixed_interior_gradient(this->_w_var, qp);
370  }
371 
372  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
373 
374  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
375  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
376 
377  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
378  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
379 
380  libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
381  libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
382  libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
383 
384  for (unsigned int i=0; i != n_u_dofs; i++)
385  {
386  Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
387  + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp] )*JxW[qp];
388 
389  Fv(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
390  + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
391 
392  if( this->_dim == 3 )
393  {
394  Fw(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
395  + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
396  }
397  }
398 
399  }
400  return;
401  }
402 
403  template<class Mu, class SH, class TC>
404  void LowMachNavierStokesSPGSMStabilization<Mu,SH,TC>::assemble_energy_mass_residual( bool /*compute_jacobian*/,
405  AssemblyContext& context )
406  {
407  // The number of local degrees of freedom in each variable.
408  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
409 
410  // Element Jacobian * quadrature weights for interior integration.
411  const std::vector<libMesh::Real> &JxW =
412  context.get_element_fe(this->_T_var)->get_JxW();
413 
414  // The temperature shape functions gradients at interior quadrature points.
415  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
416  context.get_element_fe(this->_T_var)->get_dphi();
417 
418  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var); // R_{T}
419 
420  unsigned int n_qpoints = context.get_element_qrule().n_points();
421 
422  for (unsigned int qp=0; qp != n_qpoints; qp++)
423  {
424  libMesh::Number u, v;
425  u = context.fixed_interior_value(this->_u_var, qp);
426  v = context.fixed_interior_value(this->_v_var, qp);
427 
428  libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_T_var, qp);
429 
430  libMesh::NumberVectorValue U(u,v);
431  if (this->_dim == 3)
432  U(2) = context.fixed_interior_value(this->_w_var, qp); // w
433 
434  libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
435  libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
436 
437  libMesh::Real k = this->_k(T);
438  libMesh::Real cp = this->_cp(T);
439 
440  libMesh::Number rho_cp = rho*this->_cp(T);
441 
442  libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
443 
444  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
445  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
446 
447  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, false );
448 
449  libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
450 
451  for (unsigned int i=0; i != n_T_dofs; i++)
452  {
453  FT(i) -= rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]*JxW[qp];
454  }
455 
456  }
457 
458  return;
459  }
460 
461 } // namespace GRINS
462 
463 // Instantiate
464 template class GRINS::LowMachNavierStokesSPGSMStabilization<GRINS::ConstantViscosity,GRINS::ConstantSpecificHeat,GRINS::ConstantConductivity>;
constant_conductivity.h
GRINS::LowMachNavierStokesSPGSMStabilization
Adds SPGSM-based stabilization to LowMachNavierStokes physics class.
Definition: low_mach_navier_stokes_spgsm_stab.h:35
assembly_context.h
constant_specific_heat.h
GRINS::LowMachNavierStokesSPGSMStabilization::assemble_continuity_mass_residual
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_spgsm_stab.C:278
low_mach_navier_stokes_spgsm_stab.h
GRINS::LowMachNavierStokesSPGSMStabilization::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_spgsm_stab.C:54
GRINS::LowMachNavierStokesSPGSMStabilization::assemble_momentum_mass_residual
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_spgsm_stab.C:328
GRINS::LowMachNavierStokesSPGSMStabilization::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_spgsm_stab.C:73
GRINS::LowMachNavierStokesSPGSMStabilization::assemble_momentum_time_deriv
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_spgsm_stab.C:145
GRINS
GRINS namespace.
Definition: axisym_heat_transfer_bc_handling.h:31
GRINS::LowMachNavierStokesSPGSMStabilization::assemble_energy_mass_residual
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_spgsm_stab.C:404
GRINS::LowMachNavierStokesStabilizationBase
Adds VMS-based stabilization to LowMachNavierStokes physics class.
Definition: low_mach_navier_stokes_stab_base.h:37
GRINS::LowMachNavierStokesSPGSMStabilization::assemble_continuity_time_deriv
void assemble_continuity_time_deriv(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_spgsm_stab.C:92
GRINS::LowMachNavierStokesSPGSMStabilization::assemble_energy_time_deriv
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context)
Definition: low_mach_navier_stokes_spgsm_stab.C:220
constant_viscosity.h
Generated on Mon Jun 22 2015 21:32:20 for GRINS-0.6.0 by  doxygen 1.8.9.1