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