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

Generated on Thu Jun 2 2016 21:52:28 for GRINS-0.7.0 by  doxygen 1.8.10