GRINS-0.7.0
inc_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/inc_navier_stokes.h"
28 
29 // GRINS
30 #include "grins/assembly_context.h"
31 #include "grins/generic_ic_handler.h"
32 #include "grins/postprocessed_quantities.h"
33 #include "grins/inc_nav_stokes_macro.h"
34 
35 // libMesh
36 #include "libmesh/quadrature.h"
37 
38 namespace GRINS
39 {
40 
41  template<class Mu>
42  IncompressibleNavierStokes<Mu>::IncompressibleNavierStokes(const std::string& physics_name, const GetPot& input )
43  : IncompressibleNavierStokesBase<Mu>(physics_name,
44  PhysicsNaming::incompressible_navier_stokes(), /* "core" Physics name */
45  input),
46  _p_pinning(input,physics_name),
47  _mu_index(0)
48  {
49  // This is deleted in the base class
50  this->_ic_handler = new GenericICHandler( physics_name, input );
51 
52  this->_pin_pressure = input("Physics/"+PhysicsNaming::incompressible_navier_stokes()+"/pin_pressure", false );
53  }
54 
55  template<class Mu>
56  void IncompressibleNavierStokes<Mu>::auxiliary_init( MultiphysicsSystem& system )
57  {
58  if( _pin_pressure )
59  _p_pinning.check_pin_location(system.get_mesh());
60  }
61 
62  template<class Mu>
63  void IncompressibleNavierStokes<Mu>::register_postprocessing_vars( const GetPot& input,
64  PostProcessedQuantities<libMesh::Real>& postprocessing )
65  {
66  std::string section = "Physics/"+PhysicsNaming::incompressible_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("mu") )
77  {
78  this->_mu_index = postprocessing.register_quantity( name );
79  }
80  else
81  {
82  std::cerr << "Error: Invalid output_vars value for "+PhysicsNaming::incompressible_navier_stokes() << std::endl
83  << " Found " << name << std::endl
84  << " Acceptable values are: mu" << std::endl;
85  libmesh_error();
86  }
87  }
88  }
89 
90  return;
91  }
92 
93  template<class Mu>
94  void IncompressibleNavierStokes<Mu>::element_time_derivative( bool compute_jacobian,
95  AssemblyContext& context,
96  CachedValues& /*cache*/ )
97  {
98 #ifdef GRINS_USE_GRVY_TIMERS
99  this->_timer->BeginTimer("IncompressibleNavierStokes::element_time_derivative");
100 #endif
101 
102  // The number of local degrees of freedom in each variable.
103  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
104  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
105 
106  // Check number of dofs is same for this->_flow_vars.u(), v_var and w_var.
107  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
108  if (this->_dim == 3)
109  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
110 
111  // We get some references to cell-specific data that
112  // will be used to assemble the linear system.
113 
114  // Element Jacobian * quadrature weights for interior integration.
115  const std::vector<libMesh::Real> &JxW =
116  context.get_element_fe(this->_flow_vars.u())->get_JxW();
117 
118  // The velocity shape functions at interior quadrature points.
119  const std::vector<std::vector<libMesh::Real> >& u_phi =
120  context.get_element_fe(this->_flow_vars.u())->get_phi();
121 
122  // The velocity shape function gradients (in global coords.)
123  // at interior quadrature points.
124  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
125  context.get_element_fe(this->_flow_vars.u())->get_dphi();
126 
127  // The pressure shape functions at interior quadrature points.
128  const std::vector<std::vector<libMesh::Real> >& p_phi =
129  context.get_element_fe(this->_press_var.p())->get_phi();
130 
131  const std::vector<libMesh::Point>& u_qpoint =
132  context.get_element_fe(this->_flow_vars.u())->get_xyz();
133 
134  // The subvectors and submatrices we need to fill:
135  //
136  // K_{\alpha \beta} = R_{\alpha},{\beta} = \partial{ R_{\alpha} } / \partial{ {\beta} } (where R denotes residual)
137  // e.g., for \alpha = v and \beta = u we get: K{vu} = R_{v},{u}
138  // Note that Kpu, Kpv, Kpw and Fp comes as constraint.
139 
140  libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u()); // R_{u},{u}
141  libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v()); // R_{u},{v}
142  libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
143 
144  libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u()); // R_{v},{u}
145  libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v()); // R_{v},{v}
146  libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
147 
148  libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
149  libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
150  libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
151 
152  libMesh::DenseSubMatrix<libMesh::Number> &Kup = context.get_elem_jacobian(this->_flow_vars.u(), this->_press_var.p()); // R_{u},{p}
153  libMesh::DenseSubMatrix<libMesh::Number> &Kvp = context.get_elem_jacobian(this->_flow_vars.v(), this->_press_var.p()); // R_{v},{p}
154  libMesh::DenseSubMatrix<libMesh::Number>* Kwp = NULL;
155 
156  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{u}
157  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{v}
158  libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
159 
160  if( this->_dim == 3 )
161  {
162  Kuw = &context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.w()); // R_{u},{w}
163  Kvw = &context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.w()); // R_{v},{w}
164  Kwu = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.u()); // R_{w},{u};
165  Kwv = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.v()); // R_{w},{v};
166  Kww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w()); // R_{w},{w}
167  Kwp = &context.get_elem_jacobian(this->_flow_vars.w(), this->_press_var.p()); // R_{w},{p}
168  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
169  }
170 
171  // Now we will build the element Jacobian and residual.
172  // Constructing the residual requires the solution and its
173  // gradient from the previous timestep. This must be
174  // calculated at each quadrature point by summing the
175  // solution degree-of-freedom values by the appropriate
176  // weight functions.
177  unsigned int n_qpoints = context.get_element_qrule().n_points();
178 
179  for (unsigned int qp=0; qp != n_qpoints; qp++)
180  {
181  // Compute the solution & its gradient at the old Newton iterate.
182  libMesh::Number p, u, v;
183  p = context.interior_value(this->_press_var.p(), qp);
184  u = context.interior_value(this->_flow_vars.u(), qp);
185  v = context.interior_value(this->_flow_vars.v(), qp);
186 
187  libMesh::Gradient grad_u, grad_v, grad_w;
188  grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
189  grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
190  if (this->_dim == 3)
191  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
192 
193  libMesh::NumberVectorValue U(u,v);
194  if (this->_dim == 3)
195  U(2) = context.interior_value(this->_flow_vars.w(), qp); // w
196 
197  const libMesh::Number grad_u_x = grad_u(0);
198  const libMesh::Number grad_u_y = grad_u(1);
199  const libMesh::Number grad_u_z = (this->_dim == 3)?grad_u(2):0;
200  const libMesh::Number grad_v_x = grad_v(0);
201  const libMesh::Number grad_v_y = grad_v(1);
202  const libMesh::Number grad_v_z = (this->_dim == 3)?grad_v(2):0;
203  const libMesh::Number grad_w_x = (this->_dim == 3)?grad_w(0):0;
204  const libMesh::Number grad_w_y = (this->_dim == 3)?grad_w(1):0;
205  const libMesh::Number grad_w_z = (this->_dim == 3)?grad_w(2):0;
206 
207  const libMesh::Number r = u_qpoint[qp](0);
208 
209  libMesh::Real jac = JxW[qp];
210 
211  // Compute the viscosity at this qp
212  libMesh::Real _mu_qp = this->_mu(context, qp);
213 
214  if(Physics::is_axisymmetric())
215  {
216  jac *= r;
217  }
218 
219  // First, an i-loop over the velocity degrees of freedom.
220  // We know that n_u_dofs == n_v_dofs so we can compute contributions
221  // for both at the same time.
222  for (unsigned int i=0; i != n_u_dofs; i++)
223  {
224  Fu(i) += jac *
225  (-this->_rho*u_phi[i][qp]*(U*grad_u) // convection term
226  +p*u_gradphi[i][qp](0) // pressure term
227  -_mu_qp*(u_gradphi[i][qp]*grad_u) ); // diffusion term
228 
230  if(Physics::is_axisymmetric())
231  {
232  Fu(i) += u_phi[i][qp]*( p/r - _mu_qp*U(0)/(r*r) )*jac;
233  }
234 
235  Fv(i) += jac *
236  (-this->_rho*u_phi[i][qp]*(U*grad_v) // convection term
237  +p*u_gradphi[i][qp](1) // pressure term
238  -_mu_qp*(u_gradphi[i][qp]*grad_v) ); // diffusion term
239 
240  if (this->_dim == 3)
241  {
242  (*Fw)(i) += jac *
243  (-this->_rho*u_phi[i][qp]*(U*grad_w) // convection term
244  +p*u_gradphi[i][qp](2) // pressure term
245  -_mu_qp*(u_gradphi[i][qp]*grad_w) ); // diffusion term
246  }
247 
248  if (compute_jacobian)
249  {
250  for (unsigned int j=0; j != n_u_dofs; j++)
251  {
252  // TODO: precompute some terms like:
253  // (Uvec*u_gradphi[j][qp]),
254  // u_phi[i][qp]*u_phi[j][qp],
255  // (u_gradphi[i][qp]*u_gradphi[j][qp])
256 
257  Kuu(i,j) += jac * context.get_elem_solution_derivative() *
258  (-this->_rho*u_phi[i][qp]*(U*u_gradphi[j][qp]) // convection term
259  -this->_rho*u_phi[i][qp]*grad_u_x*u_phi[j][qp] // convection term
260  -_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp])); // diffusion term
261 
262 
263  if(Physics::is_axisymmetric())
264  {
265  Kuu(i,j) -= u_phi[i][qp]*_mu_qp*u_phi[j][qp]/(r*r)*jac * context.get_elem_solution_derivative();
266  }
267 
268  Kuv(i,j) += jac * context.get_elem_solution_derivative() *
269  (-this->_rho*u_phi[i][qp]*grad_u_y*u_phi[j][qp]); // convection term
270 
271  Kvv(i,j) += jac * context.get_elem_solution_derivative() *
272  (-this->_rho*u_phi[i][qp]*(U*u_gradphi[j][qp]) // convection term
273  -this->_rho*u_phi[i][qp]*grad_v_y*u_phi[j][qp] // convection term
274  -_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp])); // diffusion term
275 
276  Kvu(i,j) += jac * context.get_elem_solution_derivative() *
277  (-this->_rho*u_phi[i][qp]*grad_v_x*u_phi[j][qp]); // convection term
278 
279  if (this->_dim == 3)
280  {
281  (*Kuw)(i,j) += jac * context.get_elem_solution_derivative() *
282  (-this->_rho*u_phi[i][qp]*grad_u_z*u_phi[j][qp]); // convection term
283 
284  (*Kvw)(i,j) += jac * context.get_elem_solution_derivative() *
285  (-this->_rho*u_phi[i][qp]*grad_v_z*u_phi[j][qp]); // convection term
286 
287  (*Kww)(i,j) += jac * context.get_elem_solution_derivative() *
288  (-this->_rho*u_phi[i][qp]*(U*u_gradphi[j][qp]) // convection term
289  -this->_rho*u_phi[i][qp]*grad_w_z*u_phi[j][qp] // convection term
290  -_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp])); // diffusion term
291  (*Kwu)(i,j) += jac * context.get_elem_solution_derivative() *
292  (-this->_rho*u_phi[i][qp]*grad_w_x*u_phi[j][qp]); // convection term
293  (*Kwv)(i,j) += jac * context.get_elem_solution_derivative() *
294  (-this->_rho*u_phi[i][qp]*grad_w_y*u_phi[j][qp]); // convection term
295  }
296  } // end of the inner dof (j) loop
297 
298  // Matrix contributions for the up, vp and wp couplings
299  for (unsigned int j=0; j != n_p_dofs; j++)
300  {
301  Kup(i,j) += u_gradphi[i][qp](0)*p_phi[j][qp]*jac * context.get_elem_solution_derivative();
302  Kvp(i,j) += u_gradphi[i][qp](1)*p_phi[j][qp]*jac * context.get_elem_solution_derivative();
303 
304  if (this->_dim == 3)
305  {
306  (*Kwp)(i,j) += u_gradphi[i][qp](2)*p_phi[j][qp]*jac * context.get_elem_solution_derivative();
307  }
308 
309  if(Physics::is_axisymmetric())
310  {
311  Kup(i,j) += u_phi[i][qp]*p_phi[j][qp]/r*jac * context.get_elem_solution_derivative();
312  }
313 
314  } // end of the inner dof (j) loop
315 
316 
317 
318  } // end - if (compute_jacobian)
319 
320  } // end of the outer dof (i) loop
321  } // end of the quadrature point (qp) loop
322 
323 #ifdef GRINS_USE_GRVY_TIMERS
324  this->_timer->EndTimer("IncompressibleNavierStokes::element_time_derivative");
325 #endif
326 
327  return;
328  }
329 
330  template<class Mu>
331  void IncompressibleNavierStokes<Mu>::element_constraint( bool compute_jacobian,
332  AssemblyContext& context,
333  CachedValues& /*cache*/ )
334  {
335 #ifdef GRINS_USE_GRVY_TIMERS
336  this->_timer->BeginTimer("IncompressibleNavierStokes::element_constraint");
337 #endif
338 
339  // The number of local degrees of freedom in each variable.
340  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
341  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
342 
343  // We get some references to cell-specific data that
344  // will be used to assemble the linear system.
345 
346  // Element Jacobian * quadrature weights for interior integration.
347  const std::vector<libMesh::Real> &JxW =
348  context.get_element_fe(this->_flow_vars.u())->get_JxW();
349 
350  // The velocity shape function gradients (in global coords.)
351  // at interior quadrature points.
352  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
353  context.get_element_fe(this->_flow_vars.u())->get_dphi();
354 
355  // The velocity shape function gradients (in global coords.)
356  // at interior quadrature points.
357  const std::vector<std::vector<libMesh::Real> >& u_phi =
358  context.get_element_fe(this->_flow_vars.u())->get_phi();
359 
360  // The pressure shape functions at interior quadrature points.
361  const std::vector<std::vector<libMesh::Real> >& p_phi =
362  context.get_element_fe(this->_press_var.p())->get_phi();
363 
364  const std::vector<libMesh::Point>& u_qpoint =
365  context.get_element_fe(this->_flow_vars.u())->get_xyz();
366 
367  // The subvectors and submatrices we need to fill:
368  //
369  // Kpu, Kpv, Kpw, Fp
370 
371  libMesh::DenseSubMatrix<libMesh::Number> &Kpu = context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u()); // R_{p},{u}
372  libMesh::DenseSubMatrix<libMesh::Number> &Kpv = context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v()); // R_{p},{v}
373  libMesh::DenseSubMatrix<libMesh::Number>* Kpw = NULL;
374 
375  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
376 
377  if( this->_dim == 3 )
378  {
379  Kpw = &context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.w()); // R_{p},{w}
380  }
381 
382  // Add the constraint given by the continuity equation.
383  unsigned int n_qpoints = context.get_element_qrule().n_points();
384  for (unsigned int qp=0; qp != n_qpoints; qp++)
385  {
386  // Compute the velocity gradient at the old Newton iterate.
387  libMesh::Gradient grad_u, grad_v, grad_w;
388  grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
389  grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
390  if (this->_dim == 3)
391  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
392 
393  libMesh::Number divU = grad_u(0) + grad_v(1);
394  if (this->_dim == 3)
395  divU += grad_w(2);
396 
397  const libMesh::Number r = u_qpoint[qp](0);
398 
399  libMesh::Real jac = JxW[qp];
400 
401  if(Physics::is_axisymmetric())
402  {
403  libMesh::Number u = context.interior_value( this->_flow_vars.u(), qp );
404  divU += u/r;
405  jac *= r;
406  }
407 
408  // Now a loop over the pressure degrees of freedom. This
409  // computes the contributions of the continuity equation.
410  for (unsigned int i=0; i != n_p_dofs; i++)
411  {
412  Fp(i) += p_phi[i][qp]*divU*jac;
413 
414  if (compute_jacobian)
415  {
416  libmesh_assert_equal_to (context.get_elem_solution_derivative(), 1.0);
417 
418  for (unsigned int j=0; j != n_u_dofs; j++)
419  {
420  Kpu(i,j) += p_phi[i][qp]*u_gradphi[j][qp](0)*jac * context.get_elem_solution_derivative();
421  Kpv(i,j) += p_phi[i][qp]*u_gradphi[j][qp](1)*jac * context.get_elem_solution_derivative();
422  if (this->_dim == 3)
423  (*Kpw)(i,j) += p_phi[i][qp]*u_gradphi[j][qp](2)*jac * context.get_elem_solution_derivative();
424 
425  if(Physics::is_axisymmetric())
426  {
427  Kpu(i,j) += p_phi[i][qp]*u_phi[j][qp]/r*jac * context.get_elem_solution_derivative();
428  }
429  } // end of the inner dof (j) loop
430 
431  } // end - if (compute_jacobian)
432 
433  } // end of the outer dof (i) loop
434  } // end of the quadrature point (qp) loop
435 
436  // Pin p = p_value at p_point
437  if( _pin_pressure )
438  {
439  _p_pinning.pin_value( context, compute_jacobian, this->_press_var.p() );
440  }
441 
442 #ifdef GRINS_USE_GRVY_TIMERS
443  this->_timer->EndTimer("IncompressibleNavierStokes::element_constraint");
444 #endif
445 
446  return;
447  }
448 
449  template<class Mu>
450  void IncompressibleNavierStokes<Mu>::mass_residual( bool compute_jacobian,
451  AssemblyContext& context,
452  CachedValues& /*cache*/ )
453  {
454  // Element Jacobian * quadrature weights for interior integration
455  // We assume the same for each flow variable
456  const std::vector<libMesh::Real> &JxW =
457  context.get_element_fe(this->_flow_vars.u())->get_JxW();
458 
459  // The shape functions at interior quadrature points.
460  // We assume the same for each flow variable
461  const std::vector<std::vector<libMesh::Real> >& u_phi =
462  context.get_element_fe(this->_flow_vars.u())->get_phi();
463 
464  const std::vector<libMesh::Point>& u_qpoint =
465  context.get_element_fe(this->_flow_vars.u())->get_xyz();
466 
467  // The number of local degrees of freedom in each variable
468  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
469 
470  // The subvectors and submatrices we need to fill:
471  libMesh::DenseSubVector<libMesh::Real> &F_u = context.get_elem_residual(this->_flow_vars.u());
472  libMesh::DenseSubVector<libMesh::Real> &F_v = context.get_elem_residual(this->_flow_vars.v());
473  libMesh::DenseSubVector<libMesh::Real>* F_w = NULL;
474 
475  libMesh::DenseSubMatrix<libMesh::Real> &M_uu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u());
476  libMesh::DenseSubMatrix<libMesh::Real> &M_vv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v());
477  libMesh::DenseSubMatrix<libMesh::Real>* M_ww = NULL;
478 
479  if( this->_dim == 3 )
480  {
481  F_w = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
482  M_ww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w());
483  }
484 
485  unsigned int n_qpoints = context.get_element_qrule().n_points();
486 
487  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
488  {
489  // For the mass residual, we need to be a little careful.
490  // The time integrator is handling the time-discretization
491  // for us so we need to supply M(u_fixed)*u' for the residual.
492  // u_fixed will be given by the fixed_interior_value function
493  // while u' will be given by the interior_rate function.
494  libMesh::Real u_dot, v_dot, w_dot = 0.0;
495  context.interior_rate(this->_flow_vars.u(), qp, u_dot);
496  context.interior_rate(this->_flow_vars.v(), qp, v_dot);
497 
498  if(this->_dim == 3 )
499  context.interior_rate(this->_flow_vars.w(), qp, w_dot);
500 
501  const libMesh::Number r = u_qpoint[qp](0);
502 
503  libMesh::Real jac = JxW[qp];
504 
505  if(Physics::is_axisymmetric())
506  {
507  jac *= r;
508  }
509 
510  for (unsigned int i = 0; i != n_u_dofs; ++i)
511  {
512  F_u(i) -= this->_rho*u_dot*u_phi[i][qp]*jac;
513  F_v(i) -= this->_rho*v_dot*u_phi[i][qp]*jac;
514 
515  if( this->_dim == 3 )
516  (*F_w)(i) -= this->_rho*w_dot*u_phi[i][qp]*jac;
517 
518  if( compute_jacobian )
519  {
520  for (unsigned int j=0; j != n_u_dofs; j++)
521  {
522  // Assuming rho is constant w.r.t. u, v, w
523  // and T (if Boussinesq added).
524  libMesh::Real value = this->_rho*u_phi[i][qp]*u_phi[j][qp]*jac * context.get_elem_solution_rate_derivative();
525 
526  M_uu(i,j) -= value;
527  M_vv(i,j) -= value;
528 
529  if( this->_dim == 3)
530  {
531  (*M_ww)(i,j) -= value;
532  }
533 
534  } // End dof loop
535  } // End Jacobian check
536  } // End dof loop
537  } // End quadrature loop
538 
539  return;
540  }
541 
542  template<class Mu>
543  void IncompressibleNavierStokes<Mu>::compute_postprocessed_quantity( unsigned int quantity_index,
544  const AssemblyContext& context,
545  const libMesh::Point& point,
546  libMesh::Real& value )
547  {
548  if( quantity_index == this->_mu_index )
549  {
550  value = this->_mu(point, context.get_time());
551  }
552 
553  return;
554  }
555 
556 } // namespace GRINS
557 
558 // Instantiate
559 INSTANTIATE_INC_NS_SUBCLASS(IncompressibleNavierStokes);
GRINS::Physics::is_axisymmetric
static bool is_axisymmetric()
Definition: physics.h:133
GRINS::IncompressibleNavierStokes::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: inc_navier_stokes.C:450
GRINS::IncompressibleNavierStokes::compute_postprocessed_quantity
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
Compute value of postprocessed quantities at libMesh::Point.
Definition: inc_navier_stokes.C:543
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::IncompressibleNavierStokes::_pin_pressure
bool _pin_pressure
Enable pressure pinning.
Definition: inc_navier_stokes.h:87
GRINS::IncompressibleNavierStokesBase
Physics class for Incompressible Navier-Stokes.
Definition: inc_navier_stokes_base.h:45
assembly_context.h
INSTANTIATE_INC_NS_SUBCLASS
INSTANTIATE_INC_NS_SUBCLASS(IncompressibleNavierStokes)
GRINS::IncompressibleNavierStokes::register_postprocessing_vars
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for IncompressibleNavierStokes.
Definition: inc_navier_stokes.C:63
GRINS::GenericICHandler
Base class for reading and handling initial conditions for physics classes.
Definition: generic_ic_handler.h:43
postprocessed_quantities.h
GRINS
GRINS namespace.
Definition: arrhenius_catalycity.h:31
inc_navier_stokes.h
GRINS::MultiphysicsSystem
Interface with libMesh for solving Multiphysics problems.
Definition: multiphysics_sys.h:82
GRINS::PhysicsNaming::incompressible_navier_stokes
static PhysicsName incompressible_navier_stokes()
Definition: physics_naming.h:65
GRINS::IncompressibleNavierStokes::element_constraint
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
Definition: inc_navier_stokes.C:331
GRINS::IncompressibleNavierStokes::auxiliary_init
virtual void auxiliary_init(MultiphysicsSystem &system)
Any auxillary initialization a Physics class may need.
Definition: inc_navier_stokes.C:56
generic_ic_handler.h
inc_nav_stokes_macro.h
GRINS::IncompressibleNavierStokes::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: inc_navier_stokes.C:94
Generated on Thu Jun 2 2016 21:52:28 for GRINS-0.7.0 by  doxygen 1.8.10