33 #include "libmesh/quadrature.h"
49 #ifdef GRINS_USE_GRVY_TIMERS
50 this->_timer->BeginTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_time_derivative");
54 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
57 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
60 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
64 const std::vector<libMesh::Real> &JxW =
65 context.get_element_fe(this->_flow_vars.u())->get_JxW();
68 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
69 context.get_element_fe(this->_flow_vars.u())->get_dphi();
71 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
72 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
73 libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
76 Fw = &context.get_elem_residual(this->_flow_vars.w());
79 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
81 unsigned int n_qpoints = context.get_element_qrule().n_points();
83 for (
unsigned int qp=0; qp != n_qpoints; qp++)
85 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
86 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
88 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
89 context.interior_value( this->_flow_vars.v(), qp ) );
92 U(2) = context.interior_value( this->_flow_vars.w(), qp );
96 libMesh::Real _mu_qp = this->_mu(context, qp);
98 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
99 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
101 libMesh::RealGradient RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
102 libMesh::Real RC = this->_stab_helper.compute_res_continuity( context, qp );
104 for (
unsigned int i=0; i != n_u_dofs; i++)
106 Fu(i) += ( - tau_C*RC*u_gradphi[i][qp](0)
107 - tau_M*RM_s(0)*this->_rho*U*u_gradphi[i][qp] )*JxW[qp];
109 Fv(i) += ( - tau_C*RC*u_gradphi[i][qp](1)
110 - tau_M*RM_s(1)*this->_rho*U*u_gradphi[i][qp] )*JxW[qp];
112 if( this->_dim == 3 )
114 (*Fw)(i) += ( - tau_C*RC*u_gradphi[i][qp](2)
115 - tau_M*RM_s(2)*this->_rho*U*u_gradphi[i][qp] )*JxW[qp];
119 if( compute_jacobian )
121 libmesh_not_implemented();
126 #ifdef GRINS_USE_GRVY_TIMERS
127 this->_timer->EndTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_time_derivative");
138 #ifdef GRINS_USE_GRVY_TIMERS
139 this->_timer->BeginTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_constraint");
143 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
146 const std::vector<libMesh::Real> &JxW =
147 context.get_element_fe(this->_flow_vars.u())->get_JxW();
150 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
151 context.get_element_fe(this->_press_var.p())->get_dphi();
153 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
155 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
157 unsigned int n_qpoints = context.get_element_qrule().n_points();
159 for (
unsigned int qp=0; qp != n_qpoints; qp++)
161 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
162 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
164 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
165 context.interior_value( this->_flow_vars.v(), qp ) );
166 if( this->_dim == 3 )
168 U(2) = context.interior_value( this->_flow_vars.w(), qp );
172 libMesh::Real _mu_qp = this->_mu(context, qp);
174 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
176 libMesh::RealGradient RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
178 for (
unsigned int i=0; i != n_p_dofs; i++)
180 Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
183 if( compute_jacobian )
185 libmesh_not_implemented();
190 #ifdef GRINS_USE_GRVY_TIMERS
191 this->_timer->EndTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_constraint");
202 #ifdef GRINS_USE_GRVY_TIMERS
203 this->_timer->BeginTimer(
"IncompressibleNavierStokesSPGSMStabilization::mass_residual");
207 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
208 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
211 const std::vector<libMesh::Real> &JxW =
212 context.get_element_fe(this->_flow_vars.u())->get_JxW();
215 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
216 context.get_element_fe(this->_press_var.p())->get_dphi();
218 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
219 context.get_element_fe(this->_flow_vars.u())->get_dphi();
221 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
222 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
223 libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
225 Fw = &context.get_elem_residual(this->_flow_vars.w());
227 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
229 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
231 unsigned int n_qpoints = context.get_element_qrule().n_points();
233 for (
unsigned int qp=0; qp != n_qpoints; qp++)
235 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
236 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
238 libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
239 context.fixed_interior_value( this->_flow_vars.v(), qp ) );
241 libMesh::Real _mu_qp = this->_mu(context, qp);
243 if( this->_dim == 3 )
245 U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
248 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp,
false );
250 libMesh::RealGradient RM_t = this->_stab_helper.compute_res_momentum_transient( context, qp, this->_rho );
252 for (
unsigned int i=0; i != n_p_dofs; i++)
254 Fp(i) -= tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
257 for (
unsigned int i=0; i != n_u_dofs; i++)
259 Fu(i) -= tau_M*RM_t(0)*this->_rho*U*u_gradphi[i][qp]*JxW[qp];
261 Fv(i) -= tau_M*RM_t(1)*this->_rho*U*u_gradphi[i][qp]*JxW[qp];
263 if( this->_dim == 3 )
265 (*Fw)(i) -= tau_M*RM_t(2)*this->_rho*U*u_gradphi[i][qp]*JxW[qp];
269 if( compute_jacobian )
271 libmesh_not_implemented();
276 #ifdef GRINS_USE_GRVY_TIMERS
277 this->_timer->EndTimer(
"IncompressibleNavierStokesSPGSMStabilization::mass_residual");
IncompressibleNavierStokesSPGSMStabilization()
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
INSTANTIATE_INC_NS_SUBCLASS(IncompressibleNavierStokesSPGSMStabilization)
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
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...