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