36 #include "libmesh/quadrature.h"
40 template<
class Mu,
class SH,
class TC>
47 template<
class Mu,
class SH,
class TC>
53 template<
class Mu,
class SH,
class TC>
58 #ifdef GRINS_USE_GRVY_TIMERS
59 this->_timer->BeginTimer(
"LowMachNavierStokesSPGSMStabilization::element_time_derivative");
62 this->assemble_continuity_time_deriv( compute_jacobian, context );
63 this->assemble_momentum_time_deriv( compute_jacobian, context );
64 this->assemble_energy_time_deriv( compute_jacobian, context );
66 #ifdef GRINS_USE_GRVY_TIMERS
67 this->_timer->EndTimer(
"LowMachNavierStokesSPGSMStabilization::element_time_derivative");
72 template<
class Mu,
class SH,
class TC>
77 #ifdef GRINS_USE_GRVY_TIMERS
78 this->_timer->BeginTimer(
"LowMachNavierStokesSPGSMStabilization::mass_residual");
81 this->assemble_continuity_mass_residual( compute_jacobian, context );
82 this->assemble_momentum_mass_residual( compute_jacobian, context );
83 this->assemble_energy_mass_residual( compute_jacobian, context );
85 #ifdef GRINS_USE_GRVY_TIMERS
86 this->_timer->EndTimer(
"LowMachNavierStokesSPGSMStabilization::mass_residual");
91 template<
class Mu,
class SH,
class TC>
96 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
99 const std::vector<libMesh::Real> &JxW =
100 context.get_element_fe(this->_flow_vars.u())->get_JxW();
103 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
104 context.get_element_fe(this->_press_var.p())->get_dphi();
106 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
108 unsigned int n_qpoints = context.get_element_qrule().n_points();
110 for (
unsigned int qp=0; qp != n_qpoints; qp++)
112 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
114 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
115 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
117 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
119 libMesh::Real mu = this->_mu(T);
121 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
123 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
124 context.interior_value( this->_flow_vars.v(), qp ) );
125 if( this->_dim == 3 )
126 U(2) = context.interior_value( this->_flow_vars.w(), qp );
128 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
130 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
134 for (
unsigned int i=0; i != n_p_dofs; i++)
136 Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
144 template<
class Mu,
class SH,
class TC>
149 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
152 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
154 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
157 const std::vector<libMesh::Real> &JxW =
158 context.get_element_fe(this->_flow_vars.u())->get_JxW();
161 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
162 context.get_element_fe(this->_flow_vars.u())->get_dphi();
164 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
165 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
166 libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
168 if( this->_dim == 3 )
170 Fw = &context.get_elem_residual(this->_flow_vars.w());
173 unsigned int n_qpoints = context.get_element_qrule().n_points();
175 for (
unsigned int qp=0; qp != n_qpoints; qp++)
177 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
178 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
180 libMesh::RealGradient U( context.interior_value(this->_flow_vars.u(), qp),
181 context.interior_value(this->_flow_vars.v(), qp) );
183 libMesh::RealGradient grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
184 libMesh::RealGradient grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
185 libMesh::RealGradient grad_w;
187 if( this->_dim == 3 )
189 U(2) = context.interior_value(this->_flow_vars.w(), qp);
190 grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
193 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
195 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
196 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
197 libMesh::Real mu = this->_mu(T);
199 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
200 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
202 libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
203 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
205 for (
unsigned int i=0; i != n_u_dofs; i++)
207 Fu(i) += ( - tau_C*RC_s*u_gradphi[i][qp](0)
208 - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp] )*JxW[qp];
210 Fv(i) += ( - tau_C*RC_s*u_gradphi[i][qp](1)
211 - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
213 if( this->_dim == 3 )
215 (*Fw)(i) += ( - tau_C*RC_s*u_gradphi[i][qp](2)
216 - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
224 template<
class Mu,
class SH,
class TC>
229 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
232 const std::vector<libMesh::Real> &JxW =
233 context.get_element_fe(this->_temp_vars.T())->get_JxW();
236 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
237 context.get_element_fe(this->_temp_vars.T())->get_dphi();
239 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
241 unsigned int n_qpoints = context.get_element_qrule().n_points();
243 for (
unsigned int qp=0; qp != n_qpoints; qp++)
245 libMesh::Number u, v;
246 u = context.interior_value(this->_flow_vars.u(), qp);
247 v = context.interior_value(this->_flow_vars.v(), qp);
249 libMesh::Gradient grad_T = context.interior_gradient(this->_temp_vars.T(), qp);
251 libMesh::NumberVectorValue U(u,v);
253 U(2) = context.interior_value(this->_flow_vars.w(), qp);
255 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
256 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
258 libMesh::Real k = this->_k(T);
259 libMesh::Real cp = this->_cp(T);
261 libMesh::Number rho_cp = rho*this->_cp(T);
263 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
265 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
266 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
268 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
270 libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
272 for (
unsigned int i=0; i != n_T_dofs; i++)
274 FT(i) -= rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]*JxW[qp];
282 template<
class Mu,
class SH,
class TC>
287 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
290 const std::vector<libMesh::Real> &JxW =
291 context.get_element_fe(this->_flow_vars.u())->get_JxW();
294 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
295 context.get_element_fe(this->_press_var.p())->get_dphi();
297 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
299 unsigned int n_qpoints = context.get_element_qrule().n_points();
301 for (
unsigned int qp=0; qp != n_qpoints; qp++)
303 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
305 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
306 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
308 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
309 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
311 libMesh::Real mu = this->_mu(T);
313 libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
314 context.fixed_interior_value( this->_flow_vars.v(), qp ) );
315 if( this->_dim == 3 )
316 U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
318 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
319 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
323 for (
unsigned int i=0; i != n_p_dofs; i++)
325 Fp(i) += tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
332 template<
class Mu,
class SH,
class TC>
337 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
340 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
342 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
345 const std::vector<libMesh::Real> &JxW =
346 context.get_element_fe(this->_flow_vars.u())->get_JxW();
349 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
350 context.get_element_fe(this->_flow_vars.u())->get_dphi();
352 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
353 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
354 libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
356 if( this->_dim == 3 )
358 Fw = &context.get_elem_residual(this->_flow_vars.w());
361 unsigned int n_qpoints = context.get_element_qrule().n_points();
362 for (
unsigned int qp=0; qp != n_qpoints; qp++)
364 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
365 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
367 libMesh::Real mu = this->_mu(T);
369 libMesh::RealGradient U( context.fixed_interior_value(this->_flow_vars.u(), qp),
370 context.fixed_interior_value(this->_flow_vars.v(), qp) );
372 libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_flow_vars.u(), qp);
373 libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_flow_vars.v(), qp);
374 libMesh::RealGradient grad_w;
376 if( this->_dim == 3 )
378 U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
379 grad_w = context.fixed_interior_gradient(this->_flow_vars.w(), qp);
382 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
384 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
385 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
387 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
388 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
390 libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
391 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
392 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
394 for (
unsigned int i=0; i != n_u_dofs; i++)
396 Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
397 + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp] )*JxW[qp];
399 Fv(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
400 + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
402 if( this->_dim == 3 )
404 (*Fw)(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
405 + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
413 template<
class Mu,
class SH,
class TC>
418 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
421 const std::vector<libMesh::Real> &JxW =
422 context.get_element_fe(this->_temp_vars.T())->get_JxW();
425 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
426 context.get_element_fe(this->_temp_vars.T())->get_dphi();
428 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
430 unsigned int n_qpoints = context.get_element_qrule().n_points();
432 for (
unsigned int qp=0; qp != n_qpoints; qp++)
434 libMesh::Number u, v;
435 u = context.fixed_interior_value(this->_flow_vars.u(), qp);
436 v = context.fixed_interior_value(this->_flow_vars.v(), qp);
438 libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_temp_vars.T(), qp);
440 libMesh::NumberVectorValue U(u,v);
442 U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
444 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
445 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
447 libMesh::Real k = this->_k(T);
448 libMesh::Real cp = this->_cp(T);
450 libMesh::Number rho_cp = rho*this->_cp(T);
452 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
454 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
455 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
457 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp,
false );
459 libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
461 for (
unsigned int i=0; i != n_T_dofs; i++)
463 FT(i) -= rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]*JxW[qp];
virtual ~LowMachNavierStokesSPGSMStabilization()
Adds SPGSM-based stabilization to LowMachNavierStokes physics class.
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
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...
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context)
LowMachNavierStokesSPGSMStabilization()
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &context)
Adds VMS-based stabilization to LowMachNavierStokes physics class.
void assemble_continuity_time_deriv(bool compute_jacobian, AssemblyContext &context)
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context)