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>
55 (
bool compute_jacobian,
58 this->assemble_continuity_time_deriv( compute_jacobian, context );
59 this->assemble_momentum_time_deriv( compute_jacobian, context );
60 this->assemble_energy_time_deriv( compute_jacobian, context );
63 template<
class Mu,
class SH,
class TC>
67 this->assemble_continuity_mass_residual( compute_jacobian, context );
68 this->assemble_momentum_mass_residual( compute_jacobian, context );
69 this->assemble_energy_mass_residual( compute_jacobian, context );
72 template<
class Mu,
class SH,
class TC>
77 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
80 const std::vector<libMesh::Real> &JxW =
81 context.get_element_fe(this->_flow_vars.u())->get_JxW();
84 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
85 context.get_element_fe(this->_press_var.p())->get_dphi();
87 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
89 unsigned int n_qpoints = context.get_element_qrule().n_points();
91 for (
unsigned int qp=0; qp != n_qpoints; qp++)
93 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
95 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
96 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
98 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
100 libMesh::Real mu = this->_mu(T);
102 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
104 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
105 context.interior_value( this->_flow_vars.v(), qp ) );
106 if( this->_flow_vars.dim() == 3 )
107 U(2) = context.interior_value( this->_flow_vars.w(), qp );
109 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
111 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
115 for (
unsigned int i=0; i != n_p_dofs; i++)
117 Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
125 template<
class Mu,
class SH,
class TC>
130 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
133 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
134 if (this->_flow_vars.dim() == 3)
135 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
138 const std::vector<libMesh::Real> &JxW =
139 context.get_element_fe(this->_flow_vars.u())->get_JxW();
142 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
143 context.get_element_fe(this->_flow_vars.u())->get_dphi();
145 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
146 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
147 libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
149 if( this->_flow_vars.dim() == 3 )
151 Fw = &context.get_elem_residual(this->_flow_vars.w());
154 unsigned int n_qpoints = context.get_element_qrule().n_points();
156 for (
unsigned int qp=0; qp != n_qpoints; qp++)
158 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
159 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
161 libMesh::RealGradient U( context.interior_value(this->_flow_vars.u(), qp),
162 context.interior_value(this->_flow_vars.v(), qp) );
164 libMesh::RealGradient grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
165 libMesh::RealGradient grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
166 libMesh::RealGradient grad_w;
168 if( this->_flow_vars.dim() == 3 )
170 U(2) = context.interior_value(this->_flow_vars.w(), qp);
171 grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
174 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
176 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
177 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
178 libMesh::Real mu = this->_mu(T);
180 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
181 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
183 libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
184 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
186 for (
unsigned int i=0; i != n_u_dofs; i++)
188 Fu(i) += ( - tau_C*RC_s*u_gradphi[i][qp](0)
189 - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp] )*JxW[qp];
191 Fv(i) += ( - tau_C*RC_s*u_gradphi[i][qp](1)
192 - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
194 if( this->_flow_vars.dim() == 3 )
196 (*Fw)(i) += ( - tau_C*RC_s*u_gradphi[i][qp](2)
197 - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
205 template<
class Mu,
class SH,
class TC>
210 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
213 const std::vector<libMesh::Real> &JxW =
214 context.get_element_fe(this->_temp_vars.T())->get_JxW();
217 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
218 context.get_element_fe(this->_temp_vars.T())->get_dphi();
220 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
222 unsigned int n_qpoints = context.get_element_qrule().n_points();
224 for (
unsigned int qp=0; qp != n_qpoints; qp++)
226 libMesh::Number u, v;
227 u = context.interior_value(this->_flow_vars.u(), qp);
228 v = context.interior_value(this->_flow_vars.v(), qp);
230 libMesh::Gradient grad_T = context.interior_gradient(this->_temp_vars.T(), qp);
232 libMesh::NumberVectorValue U(u,v);
233 if (this->_flow_vars.dim() == 3)
234 U(2) = context.interior_value(this->_flow_vars.w(), qp);
236 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
237 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
239 libMesh::Real k = this->_k(T);
240 libMesh::Real cp = this->_cp(T);
242 libMesh::Number rho_cp = rho*this->_cp(T);
244 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
246 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
247 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
249 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
251 libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
253 for (
unsigned int i=0; i != n_T_dofs; i++)
255 FT(i) -= rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]*JxW[qp];
263 template<
class Mu,
class SH,
class TC>
268 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
271 const std::vector<libMesh::Real> &JxW =
272 context.get_element_fe(this->_flow_vars.u())->get_JxW();
275 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
276 context.get_element_fe(this->_press_var.p())->get_dphi();
278 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
280 unsigned int n_qpoints = context.get_element_qrule().n_points();
282 for (
unsigned int qp=0; qp != n_qpoints; qp++)
284 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
286 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
287 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
289 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
290 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
292 libMesh::Real mu = this->_mu(T);
294 libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
295 context.fixed_interior_value( this->_flow_vars.v(), qp ) );
296 if( this->_flow_vars.dim() == 3 )
297 U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
299 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
300 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
304 for (
unsigned int i=0; i != n_p_dofs; i++)
306 Fp(i) += tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
313 template<
class Mu,
class SH,
class TC>
318 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
321 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
322 if (this->_flow_vars.dim() == 3)
323 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
326 const std::vector<libMesh::Real> &JxW =
327 context.get_element_fe(this->_flow_vars.u())->get_JxW();
330 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
331 context.get_element_fe(this->_flow_vars.u())->get_dphi();
333 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
334 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
335 libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
337 if( this->_flow_vars.dim() == 3 )
339 Fw = &context.get_elem_residual(this->_flow_vars.w());
342 unsigned int n_qpoints = context.get_element_qrule().n_points();
343 for (
unsigned int qp=0; qp != n_qpoints; qp++)
345 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
346 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
348 libMesh::Real mu = this->_mu(T);
350 libMesh::RealGradient U( context.fixed_interior_value(this->_flow_vars.u(), qp),
351 context.fixed_interior_value(this->_flow_vars.v(), qp) );
353 libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_flow_vars.u(), qp);
354 libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_flow_vars.v(), qp);
355 libMesh::RealGradient grad_w;
357 if( this->_flow_vars.dim() == 3 )
359 U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
360 grad_w = context.fixed_interior_gradient(this->_flow_vars.w(), qp);
363 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
365 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
366 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
368 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
369 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
371 libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
372 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
373 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
375 for (
unsigned int i=0; i != n_u_dofs; i++)
377 Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
378 + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp] )*JxW[qp];
380 Fv(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
381 + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
383 if( this->_flow_vars.dim() == 3 )
385 (*Fw)(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
386 + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
394 template<
class Mu,
class SH,
class TC>
399 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
402 const std::vector<libMesh::Real> &JxW =
403 context.get_element_fe(this->_temp_vars.T())->get_JxW();
406 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
407 context.get_element_fe(this->_temp_vars.T())->get_dphi();
409 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
411 unsigned int n_qpoints = context.get_element_qrule().n_points();
413 for (
unsigned int qp=0; qp != n_qpoints; qp++)
415 libMesh::Number u, v;
416 u = context.fixed_interior_value(this->_flow_vars.u(), qp);
417 v = context.fixed_interior_value(this->_flow_vars.v(), qp);
419 libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_temp_vars.T(), qp);
421 libMesh::NumberVectorValue U(u,v);
422 if (this->_flow_vars.dim() == 3)
423 U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
425 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
426 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
428 libMesh::Real k = this->_k(T);
429 libMesh::Real cp = this->_cp(T);
431 libMesh::Number rho_cp = rho*this->_cp(T);
433 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
435 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
436 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
438 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp,
false );
440 libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
442 for (
unsigned int i=0; i != n_T_dofs; i++)
444 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.
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...
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
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)
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context)
Time dependent part(s) of physics for element interiors.