36 #include "libmesh/quadrature.h"
41 template<
class Mu,
class SH,
class TC>
49 template<
class Mu,
class SH,
class TC>
55 template<
class Mu,
class SH,
class TC>
60 #ifdef GRINS_USE_GRVY_TIMERS
61 this->_timer->BeginTimer(
"LowMachNavierStokesVMSStabilization::element_time_derivative");
64 this->assemble_continuity_time_deriv( compute_jacobian, context );
65 this->assemble_momentum_time_deriv( compute_jacobian, context );
66 this->assemble_energy_time_deriv( compute_jacobian, context );
68 #ifdef GRINS_USE_GRVY_TIMERS
69 this->_timer->EndTimer(
"LowMachNavierStokesVMSStabilization::element_time_derivative");
74 template<
class Mu,
class SH,
class TC>
79 #ifdef GRINS_USE_GRVY_TIMERS
80 this->_timer->BeginTimer(
"LowMachNavierStokesVMSStabilization::mass_residual");
83 this->assemble_continuity_mass_residual( compute_jacobian, context );
84 this->assemble_momentum_mass_residual( compute_jacobian, context );
85 this->assemble_energy_mass_residual( compute_jacobian, context );
87 #ifdef GRINS_USE_GRVY_TIMERS
88 this->_timer->EndTimer(
"LowMachNavierStokesVMSStabilization::mass_residual");
93 template<
class Mu,
class SH,
class TC>
98 const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
101 const std::vector<libMesh::Real> &JxW =
102 context.get_element_fe(this->_u_var)->get_JxW();
105 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
106 context.get_element_fe(this->_p_var)->get_dphi();
108 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var);
110 unsigned int n_qpoints = context.get_element_qrule().n_points();
112 for (
unsigned int qp=0; qp != n_qpoints; qp++)
114 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
116 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
117 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
119 libMesh::Real T = context.interior_value( this->_T_var, qp );
121 libMesh::Real mu = this->_mu(T);
123 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
125 libMesh::RealGradient U( context.interior_value( this->_u_var, qp ),
126 context.interior_value( this->_v_var, qp ) );
127 if( this->_dim == 3 )
128 U(2) = context.interior_value( this->_w_var, qp );
130 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
132 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
136 for (
unsigned int i=0; i != n_p_dofs; i++)
138 Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
146 template<
class Mu,
class SH,
class TC>
151 const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
154 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
156 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
159 const std::vector<libMesh::Real> &JxW =
160 context.get_element_fe(this->_u_var)->get_JxW();
163 const std::vector<std::vector<libMesh::Real> >& u_phi =
164 context.get_element_fe(this->_u_var)->get_phi();
167 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
168 context.get_element_fe(this->_u_var)->get_dphi();
170 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var);
171 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var);
172 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var);
174 unsigned int n_qpoints = context.get_element_qrule().n_points();
176 for (
unsigned int qp=0; qp != n_qpoints; qp++)
178 libMesh::Real T = context.interior_value( this->_T_var, qp );
179 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
181 libMesh::RealGradient U( context.interior_value(this->_u_var, qp),
182 context.interior_value(this->_v_var, qp) );
184 libMesh::RealGradient grad_u = context.interior_gradient(this->_u_var, qp);
185 libMesh::RealGradient grad_v = context.interior_gradient(this->_v_var, qp);
186 libMesh::RealGradient grad_w;
188 if( this->_dim == 3 )
190 U(2) = context.interior_value(this->_w_var, qp);
191 grad_w = context.interior_gradient(this->_w_var, qp);
194 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
196 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
197 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
198 libMesh::Real mu = this->_mu(T);
200 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
201 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
203 libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
204 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
206 for (
unsigned int i=0; i != n_u_dofs; i++)
208 Fu(i) += ( -tau_C*RC_s*u_gradphi[i][qp](0)
209 - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp]
210 + rho*tau_M*RM_s*grad_u*u_phi[i][qp]
211 + tau_M*RM_s(0)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
213 Fv(i) += ( -tau_C*RC_s*u_gradphi[i][qp](1)
214 - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp]
215 + rho*tau_M*RM_s*grad_v*u_phi[i][qp]
216 + tau_M*RM_s(1)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
218 if( this->_dim == 3 )
220 Fw(i) += ( -tau_C*RC_s*u_gradphi[i][qp](2)
221 - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp]
222 + rho*tau_M*RM_s*grad_w*u_phi[i][qp]
223 + tau_M*RM_s(2)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
231 template<
class Mu,
class SH,
class TC>
236 const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
239 const std::vector<libMesh::Real> &JxW =
240 context.get_element_fe(this->_T_var)->get_JxW();
243 const std::vector<std::vector<libMesh::Real> >& T_phi =
244 context.get_element_fe(this->_T_var)->get_phi();
247 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
248 context.get_element_fe(this->_T_var)->get_dphi();
250 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var);
252 unsigned int n_qpoints = context.get_element_qrule().n_points();
254 for (
unsigned int qp=0; qp != n_qpoints; qp++)
256 libMesh::Number u, v;
257 u = context.interior_value(this->_u_var, qp);
258 v = context.interior_value(this->_v_var, qp);
260 libMesh::Gradient grad_T = context.interior_gradient(this->_T_var, qp);
262 libMesh::NumberVectorValue U(u,v);
264 U(2) = context.interior_value(this->_w_var, qp);
266 libMesh::Real T = context.interior_value( this->_T_var, qp );
267 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
269 libMesh::Real mu = this->_mu(T);
270 libMesh::Real k = this->_k(T);
271 libMesh::Real cp = this->_cp(T);
273 libMesh::Number rho_cp = rho*this->_cp(T);
275 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
277 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
278 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
280 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
281 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
283 libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
284 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
286 for (
unsigned int i=0; i != n_T_dofs; i++)
288 FT(i) += ( rho_cp*tau_M*RM_s*grad_T*T_phi[i][qp]
289 - rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]
290 + rho_cp*tau_E*RE_s*tau_M*RM_s*T_gradphi[i][qp] )*JxW[qp];
298 template<
class Mu,
class SH,
class TC>
303 const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
306 const std::vector<libMesh::Real> &JxW =
307 context.get_element_fe(this->_u_var)->get_JxW();
310 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
311 context.get_element_fe(this->_p_var)->get_dphi();
313 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var);
315 unsigned int n_qpoints = context.get_element_qrule().n_points();
317 for (
unsigned int qp=0; qp != n_qpoints; qp++)
319 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
321 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
322 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
324 libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
325 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
327 libMesh::Real mu = this->_mu(T);
329 libMesh::RealGradient U( context.fixed_interior_value( this->_u_var, qp ),
330 context.fixed_interior_value( this->_v_var, qp ) );
331 if( this->_dim == 3 )
332 U(2) = context.fixed_interior_value( this->_w_var, qp );
334 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
335 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
339 for (
unsigned int i=0; i != n_p_dofs; i++)
341 Fp(i) += tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
348 template<
class Mu,
class SH,
class TC>
353 const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
356 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
358 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
361 const std::vector<libMesh::Real> &JxW =
362 context.get_element_fe(this->_u_var)->get_JxW();
365 const std::vector<std::vector<libMesh::Real> >& u_phi =
366 context.get_element_fe(this->_u_var)->get_phi();
369 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
370 context.get_element_fe(this->_u_var)->get_dphi();
372 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var);
373 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var);
374 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var);
376 unsigned int n_qpoints = context.get_element_qrule().n_points();
377 for (
unsigned int qp=0; qp != n_qpoints; qp++)
379 libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
380 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
382 libMesh::Real mu = this->_mu(T);
384 libMesh::RealGradient U( context.fixed_interior_value(this->_u_var, qp),
385 context.fixed_interior_value(this->_v_var, qp) );
387 libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_u_var, qp);
388 libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_v_var, qp);
389 libMesh::RealGradient grad_w;
391 if( this->_dim == 3 )
393 U(2) = context.fixed_interior_value(this->_w_var, qp);
394 grad_w = context.fixed_interior_gradient(this->_w_var, qp);
397 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
399 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
400 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
402 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
403 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
405 libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
406 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
407 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
409 for (
unsigned int i=0; i != n_u_dofs; i++)
411 Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
412 + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp]
413 - rho*tau_M*RM_t*grad_u*u_phi[i][qp]
414 - tau_M*(RM_s(0)+RM_t(0))*rho*tau_M*RM_t*u_gradphi[i][qp]
415 - tau_M*RM_t(0)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
417 Fv(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
418 - rho*tau_M*RM_t*grad_v*u_phi[i][qp]
419 + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp]
420 - tau_M*(RM_s(1)+RM_t(1))*rho*tau_M*RM_t*u_gradphi[i][qp]
421 - tau_M*RM_t(1)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
423 if( this->_dim == 3 )
425 Fw(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
426 - rho*tau_M*RM_t*grad_w*u_phi[i][qp]
427 + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp]
428 - tau_M*(RM_s(2)+RM_t(2))*rho*tau_M*RM_t*u_gradphi[i][qp]
429 - tau_M*RM_t(2)*rho*tau_M*RM_s*u_gradphi[i][qp] )*JxW[qp];
437 template<
class Mu,
class SH,
class TC>
442 const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
445 const std::vector<libMesh::Real> &JxW =
446 context.get_element_fe(this->_T_var)->get_JxW();
449 const std::vector<std::vector<libMesh::Real> >& T_phi =
450 context.get_element_fe(this->_T_var)->get_phi();
453 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
454 context.get_element_fe(this->_T_var)->get_dphi();
456 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var);
458 unsigned int n_qpoints = context.get_element_qrule().n_points();
460 for (
unsigned int qp=0; qp != n_qpoints; qp++)
462 libMesh::Number u, v;
463 u = context.fixed_interior_value(this->_u_var, qp);
464 v = context.fixed_interior_value(this->_v_var, qp);
466 libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_T_var, qp);
468 libMesh::NumberVectorValue U(u,v);
470 U(2) = context.fixed_interior_value(this->_w_var, qp);
472 libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
473 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
475 libMesh::Real mu = this->_mu(T);
476 libMesh::Real k = this->_k(T);
477 libMesh::Real cp = this->_cp(T);
479 libMesh::Number rho_cp = rho*this->_cp(T);
481 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
483 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
484 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
486 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
487 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp,
false );
489 libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
490 libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
492 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
493 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
495 for (
unsigned int i=0; i != n_T_dofs; i++)
497 FT(i) -= ( -rho_cp*tau_M*RM_t*grad_T*T_phi[i][qp]
498 +rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]
499 - rho_cp*tau_E*(RE_s+RE_t)*tau_M*RM_t*T_gradphi[i][qp]
500 - rho_cp*tau_E*RE_t*tau_M*RM_s*T_gradphi[i][qp] )*JxW[qp];
Adds VMS-based stabilization to LowMachNavierStokes physics class.
LowMachNavierStokesVMSStabilization()
virtual ~LowMachNavierStokesVMSStabilization()
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...
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context)
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
Adds VMS-based stabilization to LowMachNavierStokes physics class.
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context)
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
void assemble_continuity_time_deriv(bool compute_jacobian, AssemblyContext &context)
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &context)