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