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

Generated on Tue Dec 19 2017 12:47:28 for GRINS-0.8.0 by  doxygen 1.8.9.1