GRINS-0.6.0
low_mach_navier_stokes.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.h"
28 
29 // GRINS
30 #include "grins_config.h"
31 #include "grins/assembly_context.h"
32 #include "grins/low_mach_navier_stokes_bc_handling.h"
33 #include "grins/constant_viscosity.h"
34 #include "grins/constant_specific_heat.h"
35 #include "grins/constant_conductivity.h"
36 #include "grins/generic_ic_handler.h"
37 #include "grins/postprocessed_quantities.h"
38 
39 // libMesh
40 #include "libmesh/quadrature.h"
41 
42 namespace GRINS
43 {
44 
45  template<class Mu, class SH, class TC>
46  LowMachNavierStokes<Mu,SH,TC>::LowMachNavierStokes(const std::string& physics_name, const GetPot& input)
47  : LowMachNavierStokesBase<Mu,SH,TC>(physics_name,input),
48  _p_pinning(input,physics_name),
49  _rho_index(0) // Initialize to zero
50  {
51  this->read_input_options(input);
52 
53  // This is deleted in the base class
54  this->_bc_handler = new LowMachNavierStokesBCHandling( physics_name, input );
55  this->_ic_handler = new GenericICHandler( physics_name, input );
56 
57  return;
58  }
59 
60  template<class Mu, class SH, class TC>
61  LowMachNavierStokes<Mu,SH,TC>::~LowMachNavierStokes()
62  {
63  return;
64  }
65 
66  template<class Mu, class SH, class TC>
67  void LowMachNavierStokes<Mu,SH,TC>::read_input_options( const GetPot& input )
68  {
69  // Other quantities read in base class
70 
71  // Read pressure pinning information
72  this->_pin_pressure = input("Physics/"+low_mach_navier_stokes+"/pin_pressure", false );
73 
74  return;
75  }
76 
77  template<class Mu, class SH, class TC>
78  void LowMachNavierStokes<Mu,SH,TC>::register_postprocessing_vars( const GetPot& input,
79  PostProcessedQuantities<libMesh::Real>& postprocessing )
80  {
81  std::string section = "Physics/"+low_mach_navier_stokes+"/output_vars";
82 
83  if( input.have_variable(section) )
84  {
85  unsigned int n_vars = input.vector_variable_size(section);
86 
87  for( unsigned int v = 0; v < n_vars; v++ )
88  {
89  std::string name = input(section,"DIE!",v);
90 
91  if( name == std::string("rho") )
92  {
93  this->_rho_index = postprocessing.register_quantity( name );
94  }
95  else
96  {
97  std::cerr << "Error: Invalue output_vars value for "+low_mach_navier_stokes << std::endl
98  << " Found " << name << std::endl
99  << " Acceptable values are: rho" << std::endl;
100  libmesh_error();
101  }
102  }
103  }
104 
105  return;
106  }
107 
108  template<class Mu, class SH, class TC>
109  void LowMachNavierStokes<Mu,SH,TC>::init_context( AssemblyContext& context )
110  {
111  // First call base class
112  LowMachNavierStokesBase<Mu,SH,TC>::init_context(context);
113 
114  // We also need the side shape functions, etc.
115  context.get_side_fe(this->_u_var)->get_JxW();
116  context.get_side_fe(this->_u_var)->get_phi();
117  context.get_side_fe(this->_u_var)->get_dphi();
118  context.get_side_fe(this->_u_var)->get_xyz();
119 
120  context.get_side_fe(this->_T_var)->get_JxW();
121  context.get_side_fe(this->_T_var)->get_phi();
122  context.get_side_fe(this->_T_var)->get_dphi();
123  context.get_side_fe(this->_T_var)->get_xyz();
124 
125  return;
126  }
127 
128 
129  template<class Mu, class SH, class TC>
130  void LowMachNavierStokes<Mu,SH,TC>::element_time_derivative( bool compute_jacobian,
131  AssemblyContext& context,
132  CachedValues& cache )
133  {
134 #ifdef GRINS_USE_GRVY_TIMERS
135  this->_timer->BeginTimer("LowMachNavierStokes::element_time_derivative");
136 #endif
137 
138  this->assemble_mass_time_deriv( compute_jacobian, context, cache );
139  this->assemble_momentum_time_deriv( compute_jacobian, context, cache );
140  this->assemble_energy_time_deriv( compute_jacobian, context, cache );
141 
142  if( this->_enable_thermo_press_calc )
143  this->assemble_thermo_press_elem_time_deriv( compute_jacobian, context );
144 
145 #ifdef GRINS_USE_GRVY_TIMERS
146  this->_timer->EndTimer("LowMachNavierStokes::element_time_derivative");
147 #endif
148 
149  return;
150  }
151 
152  template<class Mu, class SH, class TC>
153  void LowMachNavierStokes<Mu,SH,TC>::side_time_derivative( bool compute_jacobian,
154  AssemblyContext& context,
155  CachedValues& /*cache*/ )
156  {
157  if( this->_enable_thermo_press_calc )
158  {
159 #ifdef GRINS_USE_GRVY_TIMERS
160  this->_timer->BeginTimer("LowMachNavierStokes::side_time_derivative");
161 #endif
162 
163  this->assemble_thermo_press_side_time_deriv( compute_jacobian, context );
164 
165 #ifdef GRINS_USE_GRVY_TIMERS
166  this->_timer->EndTimer("LowMachNavierStokes::side_time_derivative");
167 #endif
168  }
169 
170  return;
171  }
172 
173  template<class Mu, class SH, class TC>
174  void LowMachNavierStokes<Mu,SH,TC>::side_constraint( bool compute_jacobian,
175  AssemblyContext& context,
176  CachedValues& /* cache */ )
177  {
178  // Pin p = p_value at p_point
179  if( this->_pin_pressure )
180  {
181  this->_p_pinning.pin_value( context, compute_jacobian, this->_p_var);
182  }
183 
184  return;
185  }
186 
187  template<class Mu, class SH, class TC>
188  void LowMachNavierStokes<Mu,SH,TC>::mass_residual( bool compute_jacobian,
189  AssemblyContext& context,
190  CachedValues& /*cache*/ )
191  {
192  this->assemble_continuity_mass_residual( compute_jacobian, context );
193 
194  this->assemble_momentum_mass_residual( compute_jacobian, context );
195 
196  this->assemble_energy_mass_residual( compute_jacobian, context );
197 
198  if( this->_enable_thermo_press_calc )
199  this->assemble_thermo_press_mass_residual( compute_jacobian, context );
200 
201  return;
202  }
203 
204  template<class Mu, class SH, class TC>
205  void LowMachNavierStokes<Mu,SH,TC>::assemble_mass_time_deriv( bool compute_jacobian,
206  AssemblyContext& context,
207  CachedValues& cache )
208  {
209  // The number of local degrees of freedom in each variable.
210  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
211 
212  // Element Jacobian * quadrature weights for interior integration.
213  const std::vector<libMesh::Real> &JxW =
214  context.get_element_fe(this->_u_var)->get_JxW();
215 
216  // The pressure shape functions at interior quadrature points.
217  const std::vector<std::vector<libMesh::Real> >& p_phi =
218  context.get_element_fe(this->_p_var)->get_phi();
219 
220 
221 
222  // The velocity shape functions at interior quadrature points.
223  const std::vector<std::vector<libMesh::Real> >& u_phi =
224  context.get_element_fe(this->_u_var)->get_phi();
225 
226  // The velocity shape function gradients (in global coords.)
227  // at interior quadrature points.
228  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
229  context.get_element_fe(this->_u_var)->get_dphi();
230 
231  // The temperature shape functions at interior quadrature points.
232  const std::vector<std::vector<libMesh::Real> >& T_phi =
233  context.get_element_fe(this->_T_var)->get_phi();
234 
235  // The temperature shape function gradients (in global coords.)
236  // at interior quadrature points.
237  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
238  context.get_element_fe(this->_T_var)->get_dphi();
239 
240 
241 
242  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var); // R_{p}
243 
244  unsigned int n_qpoints = context.get_element_qrule().n_points();
245 
246  // The number of local degrees of freedom in each variable.
247  const unsigned int n_t_dofs = context.get_dof_indices(this->_T_var).size();
248 
249  // The number of local degrees of freedom in each variable.
250  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
251 
252  // Check number of dofs is same for _u_var, v_var and w_var.
253  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
254  if (this->_dim == 3)
255  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
256 
257  for (unsigned int qp=0; qp != n_qpoints; qp++)
258  {
259  libMesh::Number u, v, T;
260  u = cache.get_cached_values(Cache::X_VELOCITY)[qp];
261  v = cache.get_cached_values(Cache::Y_VELOCITY)[qp];
262 
263  T = cache.get_cached_values(Cache::TEMPERATURE)[qp];
264 
265  libMesh::Gradient grad_u = cache.get_cached_gradient_values(Cache::X_VELOCITY_GRAD)[qp];
266  libMesh::Gradient grad_v = cache.get_cached_gradient_values(Cache::Y_VELOCITY_GRAD)[qp];
267 
268  libMesh::Gradient grad_T = cache.get_cached_gradient_values(Cache::TEMPERATURE_GRAD)[qp];
269 
270  libMesh::NumberVectorValue U(u,v);
271  if (this->_dim == 3)
272  U(2) = cache.get_cached_values(Cache::Z_VELOCITY)[qp]; // w
273 
274  libMesh::Number divU = grad_u(0) + grad_v(1);
275  if (this->_dim == 3)
276  {
277  libMesh::Gradient grad_w = cache.get_cached_gradient_values(Cache::Z_VELOCITY_GRAD)[qp];
278  divU += grad_w(2);
279  }
280 
281 
282  libMesh::DenseSubMatrix<libMesh::Number> &KPu = context.get_elem_jacobian(this->_p_var, this->_u_var);
283  libMesh::DenseSubMatrix<libMesh::Number> &KPv = context.get_elem_jacobian(this->_p_var, this->_v_var);
284  libMesh::DenseSubMatrix<libMesh::Number> &KPT = context.get_elem_jacobian(this->_p_var, this->_T_var);
285 
286  libMesh::DenseSubMatrix<libMesh::Number>* KPw = NULL;
287 
288  if( this->_dim == 3 )
289  {
290  KPw = &context.get_elem_jacobian(this->_p_var, this->_w_var);
291  }
292 
293  // Now a loop over the pressure degrees of freedom. This
294  // computes the contributions of the continuity equation.
295  for (unsigned int i=0; i != n_p_dofs; i++)
296  {
297  Fp(i) += (-U*grad_T/T + divU)*p_phi[i][qp]*JxW[qp];
298 
299 
300  if (compute_jacobian)
301  {
302 
303  for (unsigned int j=0; j!=n_u_dofs; j++)
304  {
305  KPu(i,j) += JxW[qp]*(
306  +u_gradphi[j][qp](0)*p_phi[i][qp]
307  -u_phi[j][qp]*p_phi[i][qp]*grad_T(0)/T
308  );
309 
310  KPv(i,j) += JxW[qp]*(
311  +u_gradphi[j][qp](1)*p_phi[i][qp]
312  -u_phi[j][qp]*p_phi[i][qp]*grad_T(1)/T
313  );
314 
315  if (this->_dim == 3)
316  {
317  (*KPw)(i,j) += JxW[qp]*(
318  +u_gradphi[j][qp](2)*p_phi[i][qp]
319  -u_phi[j][qp]*p_phi[i][qp]*grad_T(2)/T
320  );
321  }
322 
323  }
324 
325  for (unsigned int j=0; j!=n_t_dofs; j++)
326  {
327  KPT(i,j) += JxW[qp]*(
328  -T_gradphi[j][qp]*U*p_phi[i][qp]/T
329  +U*p_phi[i][qp]*grad_T*T_phi[j][qp])/(T*T);
330  }
331  } // end if compute_jacobian
332  } // end p_dofs loop
333  } // end qp loop
334 
335  return;
336  }
337 
338  template<class Mu, class SH, class TC>
339  void LowMachNavierStokes<Mu,SH,TC>::assemble_momentum_time_deriv( bool compute_jacobian,
340  AssemblyContext& context,
341  CachedValues& cache )
342  {
343  // The number of local degrees of freedom in each variable.
344  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
345  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
346  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
347 
348  // Check number of dofs is same for _u_var, v_var and w_var.
349  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
350  if (this->_dim == 3)
351  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
352 
353  // Element Jacobian * quadrature weights for interior integration.
354  const std::vector<libMesh::Real> &JxW =
355  context.get_element_fe(this->_u_var)->get_JxW();
356 
357  // The pressure shape functions at interior quadrature points.
358  const std::vector<std::vector<libMesh::Real> >& u_phi =
359  context.get_element_fe(this->_u_var)->get_phi();
360  const std::vector<std::vector<libMesh::Real> >& p_phi =
361  context.get_element_fe(this->_p_var)->get_phi();
362  const std::vector<std::vector<libMesh::Real> >& T_phi =
363  context.get_element_fe(this->_T_var)->get_phi();
364 
365  // The velocity shape function gradients at interior quadrature points.
366  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
367  context.get_element_fe(this->_u_var)->get_dphi();
368 
369  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var); // R_{u}
370  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var); // R_{v}
371  libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var); // R_{w}
372 
373  unsigned int n_qpoints = context.get_element_qrule().n_points();
374  for (unsigned int qp=0; qp != n_qpoints; qp++)
375  {
376  libMesh::Number u, v, p, p0, T;
377  u = cache.get_cached_values(Cache::X_VELOCITY)[qp];
378  v = cache.get_cached_values(Cache::Y_VELOCITY)[qp];
379 
380  T = cache.get_cached_values(Cache::TEMPERATURE)[qp];
381  p = cache.get_cached_values(Cache::PRESSURE)[qp];
382  p0 = cache.get_cached_values(Cache::THERMO_PRESSURE)[qp];
383 
384  libMesh::Gradient grad_u = cache.get_cached_gradient_values(Cache::X_VELOCITY_GRAD)[qp];
385  libMesh::Gradient grad_v = cache.get_cached_gradient_values(Cache::Y_VELOCITY_GRAD)[qp];
386 
387  libMesh::Gradient grad_w;
388  if (this->_dim == 3)
389  grad_w = cache.get_cached_gradient_values(Cache::Z_VELOCITY_GRAD)[qp];
390 
391  libMesh::NumberVectorValue grad_uT( grad_u(0), grad_v(0) );
392  libMesh::NumberVectorValue grad_vT( grad_u(1), grad_v(1) );
393  libMesh::NumberVectorValue grad_wT;
394  if( this->_dim == 3 )
395  {
396  grad_uT(2) = grad_w(0);
397  grad_vT(2) = grad_w(1);
398  grad_wT = libMesh::NumberVectorValue( grad_u(2), grad_v(2), grad_w(2) );
399  }
400 
401  libMesh::NumberVectorValue U(u,v);
402  if (this->_dim == 3)
403  U(2) = cache.get_cached_values(Cache::Z_VELOCITY)[qp]; // w
404 
405  libMesh::Number divU = grad_u(0) + grad_v(1);
406  if (this->_dim == 3)
407  divU += grad_w(2);
408 
409  libMesh::Number rho = this->rho( T, p0 );
410  libMesh::Number d_rho = this->d_rho_dT( T, p0 );
411  libMesh::Number d_mu = this->_mu.deriv(T);
412 
413  libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_u_var, this->_u_var); // R_{u},{u}
414  libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_u_var, this->_v_var); // R_{u},{v}
415  libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
416 
417  libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_v_var, this->_u_var); // R_{v},{u}
418  libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_v_var, this->_v_var); // R_{v},{v}
419  libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
420 
421  libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
422  libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
423  libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
424 
425  libMesh::DenseSubMatrix<libMesh::Number> &Kup = context.get_elem_jacobian(this->_u_var, this->_p_var); // R_{u},{p}
426  libMesh::DenseSubMatrix<libMesh::Number> &Kvp = context.get_elem_jacobian(this->_v_var, this->_p_var); // R_{v},{p}
427  libMesh::DenseSubMatrix<libMesh::Number>* Kwp = NULL;
428 
429  libMesh::DenseSubMatrix<libMesh::Number> &KuT = context.get_elem_jacobian(this->_u_var, this->_T_var); // R_{u},{p}
430  libMesh::DenseSubMatrix<libMesh::Number> &KvT = context.get_elem_jacobian(this->_v_var, this->_T_var); // R_{v},{p}
431  libMesh::DenseSubMatrix<libMesh::Number>* KwT = NULL;
432 
433  if( this->_dim == 3 )
434  {
435  Kuw = &context.get_elem_jacobian(this->_u_var, this->_w_var); // R_{u},{w}
436  Kvw = &context.get_elem_jacobian(this->_v_var, this->_w_var); // R_{v},{w}
437  Kwu = &context.get_elem_jacobian(this->_w_var, this->_u_var); // R_{w},{u};
438  Kwv = &context.get_elem_jacobian(this->_w_var, this->_v_var); // R_{w},{v};
439  Kww = &context.get_elem_jacobian(this->_w_var, this->_w_var); // R_{w},{w}
440  Kwp = &context.get_elem_jacobian(this->_w_var, this->_p_var); // R_{w},{p}
441  KwT = &context.get_elem_jacobian(this->_w_var, this->_T_var); // R_{w},{T}
442  }
443 
444  // Now a loop over the pressure degrees of freedom. This
445  // computes the contributions of the continuity equation.
446  for (unsigned int i=0; i != n_u_dofs; i++)
447  {
448  Fu(i) += ( -rho*U*grad_u*u_phi[i][qp] // convection term
449  + p*u_gradphi[i][qp](0) // pressure term
450  - this->_mu(T)*(u_gradphi[i][qp]*grad_u + u_gradphi[i][qp]*grad_uT
451  - 2.0/3.0*divU*u_gradphi[i][qp](0) ) // diffusion term
452  + rho*this->_g(0)*u_phi[i][qp] // hydrostatic term
453  )*JxW[qp];
454 
455  Fv(i) += ( -rho*U*grad_v*u_phi[i][qp] // convection term
456  + p*u_gradphi[i][qp](1) // pressure term
457  - this->_mu(T)*(u_gradphi[i][qp]*grad_v + u_gradphi[i][qp]*grad_vT
458  - 2.0/3.0*divU*u_gradphi[i][qp](1) ) // diffusion term
459  + rho*this->_g(1)*u_phi[i][qp] // hydrostatic term
460  )*JxW[qp];
461  if (this->_dim == 3)
462  {
463  Fw(i) += ( -rho*U*grad_w*u_phi[i][qp] // convection term
464  + p*u_gradphi[i][qp](2) // pressure term
465  - this->_mu(T)*(u_gradphi[i][qp]*grad_w + u_gradphi[i][qp]*grad_wT
466  - 2.0/3.0*divU*u_gradphi[i][qp](2) ) // diffusion term
467  + rho*this->_g(2)*u_phi[i][qp] // hydrostatic term
468  )*JxW[qp];
469  }
470 
471  if (compute_jacobian && context.get_elem_solution_derivative())
472  {
473  libmesh_assert (context.get_elem_solution_derivative() == 1.0);
474 
475  for (unsigned int j=0; j != n_u_dofs; j++)
476  {
477  //precompute repeated terms
478  libMesh::Number r0 = rho*U*u_phi[i][qp]*u_gradphi[j][qp];
479  libMesh::Number r1 = u_gradphi[i][qp]*u_gradphi[j][qp];
480  libMesh::Number r2 = rho*u_phi[i][qp]*u_phi[j][qp];
481 
482 
483  Kuu(i,j) += JxW[qp]*(
484  -r0
485  //-rho*U*u_gradphi[j][qp]*u_phi[i][qp]
486  -r2*grad_u(0)
487  //-rho*u_phi[i][qp]*grad_u(0)*u_phi[j][qp]
488  -this->_mu(T)*(
489  r1
490  //u_gradphi[i][qp]*u_gradphi[j][qp]
491  + u_gradphi[i][qp](0)*u_gradphi[j][qp](0) // transpose
492  - 2.0/3.0*u_gradphi[i][qp](0)*u_gradphi[j][qp](0)
493  ));
494  Kvv(i,j) += JxW[qp]*(
495  -r0
496  //-rho*U*u_gradphi[j][qp]*u_phi[i][qp]
497  -r2*grad_v(1)
498  //-rho*u_phi[i][qp]*grad_v(1)*u_phi[j][qp]
499  -this->_mu(T)*(
500  r1
501  //u_gradphi[i][qp]*u_gradphi[j][qp]
502  + u_gradphi[i][qp](1)*u_gradphi[j][qp](1) // transpose
503  - 2.0/3.0*u_gradphi[i][qp](1)*u_gradphi[j][qp](1)
504  ));
505 
506  Kuv(i,j) += JxW[qp]*(
507  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](1)
508  -this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](0)
509  -r2*grad_u(1)
510  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_u(1));
511  );
512 
513  Kvu(i,j) += JxW[qp]*(
514  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](0)
515  -this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](1)
516  -r2*grad_v(0)
517  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_v(0));
518  );
519 
520 
521 
522  if (this->_dim == 3)
523  {
524  (*Kuw)(i,j) += JxW[qp]*(
525  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](2)
526  -this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](0)
527  -r2*grad_u(2)
528  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_u(2)
529  );
530 
531  (*Kvw)(i,j) += JxW[qp]*(
532  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](2)
533  -this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](1)
534  -r2*grad_v(2)
535  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_v(2)
536  );
537 
538  (*Kwu)(i,j) += JxW[qp]*(
539  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](0)
540  -this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](2)
541  -r2*grad_w(0)
542  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_w(0)
543  );
544 
545  (*Kwv)(i,j) += JxW[qp]*(
546  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](1)
547  -this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](2)
548  -r2*grad_w(1)
549  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_w(1)
550  );
551 
552  (*Kww)(i,j) += JxW[qp]*(
553  -r0
554  //-rho*U*u_gradphi[j][qp]*u_phi[i][qp]
555  -r2*grad_w(2)
556  //-rho*u_phi[i][qp]*grad_w(2)*u_phi[j][qp]
557  -this->_mu(T)*(
558  r1
559  //u_gradphi[i][qp]*u_gradphi[j][qp]
560  + u_gradphi[i][qp](2)*u_gradphi[j][qp](2) // transpose
561  - 2.0/3.0*u_gradphi[i][qp](2)*u_gradphi[j][qp](2)
562  ));
563 
564  } // end if _dim==3
565  } // end of the inner dof (j) loop
566 
567  for (unsigned int j=0; j!=n_T_dofs; j++)
568  {
569 
570  //precompute repeated term
571  libMesh:: Number r3 = d_rho*u_phi[i][qp]*T_phi[j][qp];
572 
573  // Analytical Jacobains
574  KuT(i,j) += JxW[qp]*(
575  -r3*U*grad_u
576  //-d_rho*T_phi[j][qp]*U*grad_u*u_phi[i][qp]
577  -d_mu*T_phi[j][qp]*grad_u*u_gradphi[i][qp]
578  -d_mu*T_phi[j][qp]*grad_u(0)*u_gradphi[i][qp](0) // transpose
579  +2.0/3.0*d_mu*T_phi[j][qp]*divU*u_gradphi[i][qp](0)
580  +r3*this->_g(0)
581  //+d_rho*T_phi[j][qp]*this->_g(0)*u_phi[i][qp]
582  );
583 
584  KvT(i,j) += JxW[qp]*(
585  -r3*U*grad_v
586  //-d_rho*T_phi[j][qp]*U*grad_v*u_phi[i][qp]
587  -d_mu*T_phi[j][qp]*grad_v*u_gradphi[i][qp]
588  -d_mu*T_phi[j][qp]*grad_v(1)*u_gradphi[i][qp](1) // transpose
589  +2.0/3.0*d_mu*T_phi[j][qp]*divU*u_gradphi[i][qp](1)
590  +r3*this->_g(1)
591  //+d_rho*T_phi[j][qp]*this->_g(1)*u_phi[i][qp]
592  );
593 
594  if (this->_dim == 3)
595  {
596  (*KwT)(i,j) += JxW[qp]*(
597  -r3*U*grad_w
598  //-d_rho*T_phi[j][qp]*U*grad_v*u_phi[i][qp]
599  -d_mu*T_phi[j][qp]*grad_w*u_gradphi[i][qp]
600  -d_mu*T_phi[j][qp]*grad_w(2)*u_gradphi[i][qp](2) // transpose
601  +2.0/3.0*d_mu*T_phi[j][qp]*divU*u_gradphi[i][qp](2)
602  +r3*this->_g(2)
603  //+d_rho*T_phi[j][qp]*this->_g(2)*u_phi[i][qp]
604  );
605 
606  } // end if _dim==3
607  } // end T_dofs loop
608 
609  // Matrix contributions for the up, vp and wp couplings
610  for (unsigned int j=0; j != n_p_dofs; j++)
611  {
612  Kup(i,j) += JxW[qp]*p_phi[j][qp]*u_gradphi[i][qp](0);
613  Kvp(i,j) += JxW[qp]*p_phi[j][qp]*u_gradphi[i][qp](1);
614  if (this->_dim == 3)
615  (*Kwp)(i,j) += JxW[qp]*p_phi[j][qp]*u_gradphi[i][qp](2);
616  } // end of the inner dof (j) loop
617 
618  } // end - if (compute_jacobian && context.get_elem_solution_derivative())
619 
620  } // end of the outer dof (i) loop
621  } // end of the quadrature point (qp) loop
622 
623  return;
624  }
625 
626  template<class Mu, class SH, class TC>
627  void LowMachNavierStokes<Mu,SH,TC>::assemble_energy_time_deriv( bool compute_jacobian,
628  AssemblyContext& context,
629  CachedValues& cache )
630  {
631  // The number of local degrees of freedom in each variable.
632  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
633  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
634 
635  // Element Jacobian * quadrature weights for interior integration.
636  const std::vector<libMesh::Real> &JxW =
637  context.get_element_fe(this->_T_var)->get_JxW();
638 
639  // The temperature shape functions at interior quadrature points.
640  const std::vector<std::vector<libMesh::Real> >& T_phi =
641  context.get_element_fe(this->_T_var)->get_phi();
642  const std::vector<std::vector<libMesh::Real> >& u_phi =
643  context.get_element_fe(this->_u_var)->get_phi();
644 
645  // The temperature shape functions gradients at interior quadrature points.
646  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
647  context.get_element_fe(this->_T_var)->get_dphi();
648 
649  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var); // R_{T}
650 
651  unsigned int n_qpoints = context.get_element_qrule().n_points();
652  for (unsigned int qp=0; qp != n_qpoints; qp++)
653  {
654  libMesh::Number u, v, T, p0;
655  u = cache.get_cached_values(Cache::X_VELOCITY)[qp];
656  v = cache.get_cached_values(Cache::Y_VELOCITY)[qp];
657  T = cache.get_cached_values(Cache::TEMPERATURE)[qp];
658  p0 = cache.get_cached_values(Cache::THERMO_PRESSURE)[qp];
659 
660  libMesh::Gradient grad_T = cache.get_cached_gradient_values(Cache::TEMPERATURE_GRAD)[qp];
661 
662  libMesh::NumberVectorValue U(u,v);
663  if (this->_dim == 3)
664  U(2) = cache.get_cached_values(Cache::Z_VELOCITY)[qp]; // w
665 
666  libMesh::DenseSubMatrix<libMesh::Number> &KTu = context.get_elem_jacobian(this->_T_var, this->_u_var); // R_{u},{u}
667  libMesh::DenseSubMatrix<libMesh::Number> &KTv = context.get_elem_jacobian(this->_T_var, this->_v_var); // R_{u},{u}
668  libMesh::DenseSubMatrix<libMesh::Number>* KTw = NULL;
669 
670  libMesh::DenseSubMatrix<libMesh::Number> &KTT = context.get_elem_jacobian(this->_T_var, this->_T_var); // R_{u},{u}
671 
672  if( this->_dim == 3 )
673  {
674  KTw = &context.get_elem_jacobian(this->_T_var, this->_w_var); // R_{u},{w}
675  }
676 
677  libMesh::Number k = this->_k(T);
678  libMesh::Number dk_dT = this->_k.deriv(T);
679  libMesh::Number cp = this->_cp(T);
680  libMesh::Number d_cp = this->_cp.deriv(T);
681  libMesh::Number rho = this->rho( T, p0 );
682  libMesh::Number d_rho = this->d_rho_dT( T, p0 );
683 
684  // Now a loop over the pressure degrees of freedom. This
685  // computes the contributions of the continuity equation.
686  for (unsigned int i=0; i != n_T_dofs; i++)
687  {
688  FT(i) += ( -rho*cp*U*grad_T*T_phi[i][qp] // convection term
689  - k*grad_T*T_gradphi[i][qp] // diffusion term
690  )*JxW[qp];
691 
692 
693  if(compute_jacobian)
694  {
695  for (unsigned int j=0; j!=n_u_dofs; j++)
696  {
697  //pre-compute repeated term
698  libMesh::Number r0 = rho*cp*T_phi[i][qp]*u_phi[j][qp];
699 
700  KTu(i,j) += JxW[qp]*
701  -r0*grad_T(0);
702  //-rho*cp*u_phi[j][qp]*grad_T(0)*T_phi[i][qp];
703 
704  KTv(i,j) += JxW[qp]*
705  -r0*grad_T(1);
706  //-rho*cp*u_phi[j][qp]*grad_T(1)*T_phi[i][qp];
707 
708  if (this->_dim == 3)
709  {
710  (*KTw)(i,j) += JxW[qp]*
711  -r0*grad_T(2);
712  //-rho*cp*u_phi[j][qp]*grad_T(2)*T_phi[i][qp];
713  }
714 
715  } // end u_dofs loop (j)
716 
717 
718  for (unsigned int j=0; j!=n_T_dofs; j++)
719  {
720  KTT(i,j) += JxW[qp]* (
721  -rho*(
722  cp*U*T_phi[i][qp]*T_gradphi[j][qp]
723  + U*grad_T*T_phi[i][qp]*d_cp*T_phi[j][qp]
724  )
725  -cp*U*grad_T*T_phi[i][qp]*d_rho*T_phi[j][qp]
726  -k*T_gradphi[i][qp]*T_gradphi[j][qp]
727  -grad_T*T_gradphi[i][qp]*dk_dT*T_phi[j][qp]
728  );
729  } // end T_dofs loop (j)
730 
731  } // end if compute_jacobian
732  } // end outer T_dofs loop (i)
733  } //end qp loop
734 
735  return;
736  }
737 
738  template<class Mu, class SH, class TC>
739  void LowMachNavierStokes<Mu,SH,TC>::assemble_continuity_mass_residual( bool /*compute_jacobian*/,
740  AssemblyContext& context )
741  {
742  // Element Jacobian * quadrature weights for interior integration
743  const std::vector<libMesh::Real> &JxW =
744  context.get_element_fe(this->_u_var)->get_JxW();
745 
746  // The shape functions at interior quadrature points.
747  const std::vector<std::vector<libMesh::Real> >& p_phi =
748  context.get_element_fe(this->_p_var)->get_phi();
749 
750  // The number of local degrees of freedom in each variable
751  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
752 
753  // The subvectors and submatrices we need to fill:
754  libMesh::DenseSubVector<libMesh::Real> &F_p = context.get_elem_residual(this->_p_var);
755 
756  unsigned int n_qpoints = context.get_element_qrule().n_points();
757 
758  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
759  {
760  // For the mass residual, we need to be a little careful.
761  // The time integrator is handling the time-discretization
762  // for us so we need to supply M(u_fixed)*u' for the residual.
763  // u_fixed will be given by the fixed_interior_value function
764  // while u' will be given by the interior_rate function.
765  libMesh::Real T_dot;
766  context.interior_rate(this->_T_var, qp, T_dot);
767 
768  libMesh::Real T = context.fixed_interior_value(this->_T_var, qp);
769 
770  for (unsigned int i = 0; i != n_p_dofs; ++i)
771  {
772  F_p(i) -= T_dot/T*p_phi[i][qp]*JxW[qp];
773  } // End DoF loop i
774 
775  } // End quadrature loop qp
776 
777  return;
778  }
779 
780  template<class Mu, class SH, class TC>
781  void LowMachNavierStokes<Mu,SH,TC>::assemble_momentum_mass_residual( bool /*compute_jacobian*/,
782  AssemblyContext& context )
783  {
784  // Element Jacobian * quadrature weights for interior integration
785  const std::vector<libMesh::Real> &JxW =
786  context.get_element_fe(this->_u_var)->get_JxW();
787 
788  // The shape functions at interior quadrature points.
789  const std::vector<std::vector<libMesh::Real> >& u_phi =
790  context.get_element_fe(this->_u_var)->get_phi();
791 
792  // The number of local degrees of freedom in each variable
793  const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
794 
795  // for convenience
796  if (this->_dim != 3)
797  this->_w_var = this->_u_var;
798 
799  // The subvectors and submatrices we need to fill:
800  libMesh::DenseSubVector<libMesh::Real> &F_u = context.get_elem_residual(this->_u_var);
801  libMesh::DenseSubVector<libMesh::Real> &F_v = context.get_elem_residual(this->_v_var);
802  libMesh::DenseSubVector<libMesh::Real> &F_w = context.get_elem_residual(this->_w_var);
803 
804  unsigned int n_qpoints = context.get_element_qrule().n_points();
805 
806  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
807  {
808  // For the mass residual, we need to be a little careful.
809  // The time integrator is handling the time-discretization
810  // for us so we need to supply M(u_fixed)*u' for the residual.
811  // u_fixed will be given by the fixed_interior_value function
812  // while u' will be given by the interior_rate function.
813  libMesh::Real u_dot, v_dot, w_dot = 0.0;
814  context.interior_rate(this->_u_var, qp, u_dot);
815  context.interior_rate(this->_v_var, qp, v_dot);
816 
817  if( this->_dim == 3 )
818  context.interior_rate(this->_w_var, qp, w_dot);
819 
820  libMesh::Real T = context.fixed_interior_value(this->_T_var, qp);
821 
822  libMesh::Number rho = this->rho(T, this->get_p0_transient(context, qp));
823 
824  for (unsigned int i = 0; i != n_u_dofs; ++i)
825  {
826  F_u(i) -= rho*u_dot*u_phi[i][qp]*JxW[qp];
827  F_v(i) -= rho*v_dot*u_phi[i][qp]*JxW[qp];
828 
829  if( this->_dim == 3 )
830  F_w(i) -= rho*w_dot*u_phi[i][qp]*JxW[qp];
831 
832  /*
833  if( compute_jacobian )
834  {
835  for (unsigned int j=0; j != n_u_dofs; j++)
836  {
837  // Assuming rho is constant w.r.t. u, v, w
838  // and T (if Boussinesq added).
839  libMesh::Real value = JxW[qp]*_rho*u_phi[i][qp]*u_phi[j][qp];
840 
841  M_uu(i,j) += value;
842  M_vv(i,j) += value;
843 
844  if( _dim == 3)
845  {
846  M_ww(i,j) += value;
847  }
848 
849  } // End DoF loop j
850  } // End Jacobian check
851  */
852 
853  } // End DoF loop i
854  } // End quadrature loop qp
855 
856  return;
857  }
858 
859  template<class Mu, class SH, class TC>
860  void LowMachNavierStokes<Mu,SH,TC>::assemble_energy_mass_residual( bool /*compute_jacobian*/,
861  AssemblyContext& context )
862  {
863  // Element Jacobian * quadrature weights for interior integration
864  const std::vector<libMesh::Real> &JxW =
865  context.get_element_fe(this->_T_var)->get_JxW();
866 
867  // The shape functions at interior quadrature points.
868  const std::vector<std::vector<libMesh::Real> >& T_phi =
869  context.get_element_fe(this->_T_var)->get_phi();
870 
871  // The number of local degrees of freedom in each variable
872  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
873 
874  // The subvectors and submatrices we need to fill:
875  libMesh::DenseSubVector<libMesh::Real> &F_T = context.get_elem_residual(this->_T_var);
876 
877  unsigned int n_qpoints = context.get_element_qrule().n_points();
878 
879  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
880  {
881  // For the mass residual, we need to be a little careful.
882  // The time integrator is handling the time-discretization
883  // for us so we need to supply M(u_fixed)*u' for the residual.
884  // u_fixed will be given by the fixed_interior_value function
885  // while u will be given by the interior_rate function.
886  libMesh::Real T_dot;
887  context.interior_rate(this->_T_var, qp, T_dot);
888 
889  libMesh::Real T = context.fixed_interior_value(this->_T_var, qp);
890 
891  libMesh::Real cp = this->_cp(T);
892 
893  libMesh::Number rho = this->rho(T, this->get_p0_transient(context, qp));
894 
895  for (unsigned int i = 0; i != n_T_dofs; ++i)
896  {
897  F_T(i) -= rho*cp*T_dot*T_phi[i][qp]*JxW[qp];
898  } // End DoF loop i
899 
900  } // End quadrature loop qp
901 
902  return;
903  }
904 
905  template<class Mu, class SH, class TC>
906  void LowMachNavierStokes<Mu,SH,TC>::assemble_thermo_press_elem_time_deriv( bool /*compute_jacobian*/,
907  AssemblyContext& context )
908  {
909  // Element Jacobian * quadrature weights for interior integration
910  const std::vector<libMesh::Real> &JxW =
911  context.get_element_fe(this->_T_var)->get_JxW();
912 
913  // The number of local degrees of freedom in each variable
914  const unsigned int n_p0_dofs = context.get_dof_indices(this->_p0_var).size();
915 
916  // The subvectors and submatrices we need to fill:
917  libMesh::DenseSubVector<libMesh::Real> &F_p0 = context.get_elem_residual(this->_p0_var);
918 
919  unsigned int n_qpoints = context.get_element_qrule().n_points();
920 
921  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
922  {
923  libMesh::Number T;
924  T = context.interior_value(this->_T_var, qp);
925 
926  libMesh::Gradient grad_u, grad_v, grad_w;
927  grad_u = context.interior_gradient(this->_u_var, qp);
928  grad_v = context.interior_gradient(this->_v_var, qp);
929  if (this->_dim == 3)
930  grad_w = context.interior_gradient(this->_w_var, qp);
931 
932  libMesh::Number divU = grad_u(0) + grad_v(1);
933  if(this->_dim==3)
934  divU += grad_w(2);
935 
936  //libMesh::Number cp = this->_cp(T);
937  //libMesh::Number cv = cp + this->_R;
938  //libMesh::Number gamma = cp/cv;
939  //libMesh::Number gamma_ratio = gamma/(gamma-1.0);
940 
941  libMesh::Number p0 = context.interior_value( this->_p0_var, qp );
942 
943  for (unsigned int i = 0; i != n_p0_dofs; ++i)
944  {
945  F_p0(i) += (p0/T - this->_p0/this->_T0)*JxW[qp];
946  //F_p0(i) -= p0*gamma_ratio*divU*JxW[qp];
947  } // End DoF loop i
948  }
949 
950  return;
951  }
952 
953  template<class Mu, class SH, class TC>
954  void LowMachNavierStokes<Mu,SH,TC>::assemble_thermo_press_side_time_deriv( bool /*compute_jacobian*/,
955  AssemblyContext& context )
956  {
957  // The number of local degrees of freedom in each variable.
958  const unsigned int n_p0_dofs = context.get_dof_indices(this->_p0_var).size();
959 
960  // Element Jacobian * quadrature weight for side integration.
961  //const std::vector<libMesh::Real> &JxW_side = context.get_side_fe(this->_T_var)->get_JxW();
962 
963  //const std::vector<Point> &normals = context.get_side_fe(this->_T_var)->get_normals();
964 
965  //libMesh::DenseSubVector<libMesh::Number> &F_p0 = context.get_elem_residual(this->_p0_var); // residual
966 
967  // Physical location of the quadrature points
968  //const std::vector<libMesh::Point>& qpoint = context.get_side_fe(this->_T_var)->get_xyz();
969 
970  unsigned int n_qpoints = context.get_side_qrule().n_points();
971  for (unsigned int qp=0; qp != n_qpoints; qp++)
972  {
973  /*
974  libMesh::Number T = context.side_value( this->_T_var, qp );
975  libMesh::Gradient U = ( context.side_value( this->_u_var, qp ),
976  context.side_value( this->_v_var, qp ) );
977  libMesh::Gradient grad_T = context.side_gradient( this->_T_var, qp );
978 
979 
980  libMesh::Number p0 = context.side_value( this->_p0_var, qp );
981 
982  libMesh::Number k = this->_k(T);
983  libMesh::Number cp = this->_cp(T);
984 
985  libMesh::Number cv = cp + this->_R;
986  libMesh::Number gamma = cp/cv;
987  libMesh::Number gamma_ratio = gamma/(gamma-1.0);
988  */
989 
990  //std::cout << "U = " << U << std::endl;
991 
992  //std::cout << "x = " << qpoint[qp] << ", grad_T = " << grad_T << std::endl;
993 
994  for (unsigned int i=0; i != n_p0_dofs; i++)
995  {
996  //F_p0(i) += (k*grad_T*normals[qp] - p0*gamma_ratio*U*normals[qp] )*JxW_side[qp];
997  }
998  }
999 
1000  return;
1001  }
1002 
1003  template<class Mu, class SH, class TC>
1004  void LowMachNavierStokes<Mu,SH,TC>::assemble_thermo_press_mass_residual( bool /*compute_jacobian*/,
1005  AssemblyContext& context )
1006  {
1007  // The number of local degrees of freedom in each variable.
1008  const unsigned int n_p0_dofs = context.get_dof_indices(this->_p0_var).size();
1009  const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
1010  const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
1011 
1012  // Element Jacobian * quadrature weights for interior integration
1013  const std::vector<libMesh::Real> &JxW =
1014  context.get_element_fe(this->_T_var)->get_JxW();
1015 
1016  // The temperature shape functions at interior quadrature points.
1017  const std::vector<std::vector<libMesh::Real> >& T_phi =
1018  context.get_element_fe(this->_T_var)->get_phi();
1019 
1020  // The temperature shape functions at interior quadrature points.
1021  const std::vector<std::vector<libMesh::Real> >& p_phi =
1022  context.get_element_fe(this->_p_var)->get_phi();
1023 
1024  // The subvectors and submatrices we need to fill:
1025  libMesh::DenseSubVector<libMesh::Real> &F_p0 = context.get_elem_residual(this->_p0_var);
1026  libMesh::DenseSubVector<libMesh::Real> &F_T = context.get_elem_residual(this->_T_var);
1027  libMesh::DenseSubVector<libMesh::Real> &F_p = context.get_elem_residual(this->_p_var);
1028 
1029  unsigned int n_qpoints = context.get_element_qrule().n_points();
1030 
1031  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
1032  {
1033  libMesh::Number T;
1034  T = context.fixed_interior_value(this->_T_var, qp);
1035 
1036  libMesh::Number cp = this->_cp(T);
1037  libMesh::Number cv = cp + this->_R;
1038  libMesh::Number gamma = cp/cv;
1039  libMesh::Number one_over_gamma = 1.0/(gamma-1.0);
1040 
1041  libMesh::Number p0_dot;
1042  context.interior_rate(this->_p0_var, qp, p0_dot);
1043 
1044  libMesh::Number p0 = context.fixed_interior_value(this->_p0_var, qp );
1045 
1046  for (unsigned int i=0; i != n_p0_dofs; i++)
1047  {
1048  F_p0(i) -= p0_dot*one_over_gamma*JxW[qp];
1049  }
1050 
1051  for (unsigned int i=0; i != n_T_dofs; i++)
1052  {
1053  F_T(i) += p0_dot*T_phi[i][qp]*JxW[qp];
1054  }
1055 
1056  for (unsigned int i=0; i != n_p_dofs; i++)
1057  {
1058  F_p(i) += p0_dot/p0*p_phi[i][qp]*JxW[qp];
1059  }
1060 
1061  }
1062  return;
1063  }
1064 
1065  template<class Mu, class SH, class TC>
1066  void LowMachNavierStokes<Mu,SH,TC>::compute_element_time_derivative_cache( const AssemblyContext& context,
1067  CachedValues& cache )
1068  {
1069  const unsigned int n_qpoints = context.get_element_qrule().n_points();
1070 
1071  std::vector<libMesh::Real> u, v, w, T, p, p0;
1072  u.resize(n_qpoints);
1073  v.resize(n_qpoints);
1074  if( this->_dim > 2 )
1075  w.resize(n_qpoints);
1076 
1077  T.resize(n_qpoints);
1078  p.resize(n_qpoints);
1079  p0.resize(n_qpoints);
1080 
1081  std::vector<libMesh::Gradient> grad_u, grad_v, grad_w, grad_T;
1082  grad_u.resize(n_qpoints);
1083  grad_v.resize(n_qpoints);
1084  if( this->_dim > 2 )
1085  grad_w.resize(n_qpoints);
1086 
1087  grad_T.resize(n_qpoints);
1088 
1089  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
1090  {
1091  u[qp] = context.interior_value(this->_u_var, qp);
1092  v[qp] = context.interior_value(this->_v_var, qp);
1093 
1094  grad_u[qp] = context.interior_gradient(this->_u_var, qp);
1095  grad_v[qp] = context.interior_gradient(this->_v_var, qp);
1096  if( this->_dim > 2 )
1097  {
1098  w[qp] = context.interior_value(this->_w_var, qp);
1099  grad_w[qp] = context.interior_gradient(this->_w_var, qp);
1100  }
1101  T[qp] = context.interior_value(this->_T_var, qp);
1102  grad_T[qp] = context.interior_gradient(this->_T_var, qp);
1103 
1104  p[qp] = context.interior_value(this->_p_var, qp);
1105  p0[qp] = this->get_p0_steady(context, qp);
1106  }
1107 
1108  cache.set_values(Cache::X_VELOCITY, u);
1109  cache.set_values(Cache::Y_VELOCITY, v);
1110 
1111  cache.set_gradient_values(Cache::X_VELOCITY_GRAD, grad_u);
1112  cache.set_gradient_values(Cache::Y_VELOCITY_GRAD, grad_v);
1113 
1114  if(this->_dim > 2)
1115  {
1116  cache.set_values(Cache::Z_VELOCITY, w);
1117  cache.set_gradient_values(Cache::Z_VELOCITY_GRAD, grad_w);
1118  }
1119 
1120  cache.set_values(Cache::TEMPERATURE, T);
1121  cache.set_gradient_values(Cache::TEMPERATURE_GRAD, grad_T);
1122 
1123  cache.set_values(Cache::PRESSURE, p);
1124  cache.set_values(Cache::THERMO_PRESSURE, p0);
1125 
1126  return;
1127  }
1128 
1129  template<class Mu, class SH, class TC>
1130  void LowMachNavierStokes<Mu,SH,TC>::compute_postprocessed_quantity( unsigned int quantity_index,
1131  const AssemblyContext& context,
1132  const libMesh::Point& point,
1133  libMesh::Real& value )
1134  {
1135  if( quantity_index == this->_rho_index )
1136  {
1137  libMesh::Real T = this->T(point,context);
1138  libMesh::Real p0 = this->get_p0_steady(context,point);
1139 
1140  value = this->rho( T, p0 );
1141  }
1142 
1143  return;
1144  }
1145 
1146 } // namespace GRINS
1147 
1148 // Instantiate
1149 template class GRINS::LowMachNavierStokes<GRINS::ConstantViscosity,GRINS::ConstantSpecificHeat,GRINS::ConstantConductivity>;
GRINS::LowMachNavierStokes::compute_element_time_derivative_cache
virtual void compute_element_time_derivative_cache(const AssemblyContext &context, CachedValues &cache)
Definition: low_mach_navier_stokes.C:1066
constant_conductivity.h
GRINS::LowMachNavierStokes::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.C:130
GRINS::PostProcessedQuantities::register_quantity
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
Definition: postprocessed_quantities.C:58
GRINS::Physics::_ic_handler
GRINS::ICHandlingBase * _ic_handler
Definition: physics.h:258
GRINS::LowMachNavierStokesBCHandling
Base class for reading and handling boundary conditions for physics classes.
Definition: low_mach_navier_stokes_bc_handling.h:35
GRINS::CachedValues::set_values
void set_values(unsigned int quantity, std::vector< libMesh::Number > &values)
Definition: cached_values.C:72
GRINS::Physics::_bc_handler
GRINS::BCHandlingBase * _bc_handler
Definition: physics.h:256
GRINS::LowMachNavierStokes::assemble_energy_time_deriv
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Helper function.
Definition: low_mach_navier_stokes.C:627
assembly_context.h
GRINS::LowMachNavierStokes::side_time_derivative
virtual void side_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for boundaries of elements on the domain boundary.
Definition: low_mach_navier_stokes.C:153
constant_specific_heat.h
GRINS::LowMachNavierStokes::read_input_options
virtual void read_input_options(const GetPot &input)
Read options from GetPot input file.
Definition: low_mach_navier_stokes.C:67
low_mach_navier_stokes.h
GRINS::GenericICHandler
Base class for reading and handling initial conditions for physics classes.
Definition: generic_ic_handler.h:44
GRINS::LowMachNavierStokes::compute_postprocessed_quantity
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
Definition: low_mach_navier_stokes.C:1130
postprocessed_quantities.h
GRINS::LowMachNavierStokes::assemble_thermo_press_mass_residual
void assemble_thermo_press_mass_residual(bool compute_jacobian, AssemblyContext &c)
Definition: low_mach_navier_stokes.C:1004
GRINS
GRINS namespace.
Definition: axisym_heat_transfer_bc_handling.h:31
GRINS::LowMachNavierStokes::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.C:188
GRINS::LowMachNavierStokes
Physics class for Incompressible Navier-Stokes.
Definition: low_mach_navier_stokes.h:41
GRINS::LowMachNavierStokes::assemble_thermo_press_elem_time_deriv
void assemble_thermo_press_elem_time_deriv(bool compute_jacobian, AssemblyContext &c)
Definition: low_mach_navier_stokes.C:906
GRINS::LowMachNavierStokes::assemble_momentum_time_deriv
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Helper function.
Definition: low_mach_navier_stokes.C:339
GRINS::LowMachNavierStokes::assemble_momentum_mass_residual
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &c)
Helper function.
Definition: low_mach_navier_stokes.C:781
GRINS::CachedValues::get_cached_gradient_values
const std::vector< libMesh::Gradient > & get_cached_gradient_values(unsigned int quantity) const
Definition: cached_values.C:105
GRINS::low_mach_navier_stokes
const PhysicsName low_mach_navier_stokes
Definition: grins_physics_names.h:67
GRINS::LowMachNavierStokes::side_constraint
virtual void side_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for boundaries of elements on the domain boundary.
Definition: low_mach_navier_stokes.C:174
GRINS::LowMachNavierStokes::assemble_mass_time_deriv
void assemble_mass_time_deriv(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Helper function.
Definition: low_mach_navier_stokes.C:205
GRINS::LowMachNavierStokes::register_postprocessing_vars
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for LowMachNavierStokes.
Definition: low_mach_navier_stokes.C:78
low_mach_navier_stokes_bc_handling.h
generic_ic_handler.h
GRINS::LowMachNavierStokesBase
Physics class for Incompressible Navier-Stokes.
Definition: low_mach_navier_stokes_base.h:47
GRINS::LowMachNavierStokes::init_context
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
Definition: low_mach_navier_stokes.C:109
GRINS::LowMachNavierStokesBase::init_context
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
Definition: low_mach_navier_stokes_base.C:175
GRINS::LowMachNavierStokes::assemble_energy_mass_residual
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &c)
Helper function.
Definition: low_mach_navier_stokes.C:860
GRINS::CachedValues::set_gradient_values
void set_gradient_values(unsigned int quantity, std::vector< libMesh::Gradient > &values)
Definition: cached_values.C:78
GRINS::LowMachNavierStokes::assemble_continuity_mass_residual
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &c)
Helper function.
Definition: low_mach_navier_stokes.C:739
GRINS::CachedValues::get_cached_values
const std::vector< libMesh::Number > & get_cached_values(unsigned int quantity) const
Definition: cached_values.C:99
constant_viscosity.h
GRINS::LowMachNavierStokes::assemble_thermo_press_side_time_deriv
void assemble_thermo_press_side_time_deriv(bool compute_jacobian, AssemblyContext &c)
Definition: low_mach_navier_stokes.C:954
Generated on Mon Jun 22 2015 21:32:20 for GRINS-0.6.0 by  doxygen 1.8.9.1