33 #include "libmesh/quadrature.h"
59 #ifdef GRINS_USE_GRVY_TIMERS
60 this->_timer->BeginTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_time_derivative");
64 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
67 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v_var()).size());
70 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w_var()).size());
74 const std::vector<libMesh::Real> &JxW =
75 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
78 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
79 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
81 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
82 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
83 libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
86 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
89 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
91 unsigned int n_qpoints = context.get_element_qrule().n_points();
93 for (
unsigned int qp=0; qp != n_qpoints; qp++)
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::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
99 context.interior_value( this->_flow_vars.v_var(), qp ) );
100 if( this->_dim == 3 )
102 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
106 libMesh::Real _mu_qp = this->_mu(context, qp);
108 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
109 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
111 libMesh::RealGradient RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
112 libMesh::Real RC = this->_stab_helper.compute_res_continuity( context, qp );
114 for (
unsigned int i=0; i != n_u_dofs; i++)
116 Fu(i) += ( - tau_C*RC*u_gradphi[i][qp](0)
117 - tau_M*RM_s(0)*this->_rho*U*u_gradphi[i][qp] )*JxW[qp];
119 Fv(i) += ( - tau_C*RC*u_gradphi[i][qp](1)
120 - tau_M*RM_s(1)*this->_rho*U*u_gradphi[i][qp] )*JxW[qp];
122 if( this->_dim == 3 )
124 (*Fw)(i) += ( - tau_C*RC*u_gradphi[i][qp](2)
125 - tau_M*RM_s(2)*this->_rho*U*u_gradphi[i][qp] )*JxW[qp];
129 if( compute_jacobian )
131 libmesh_not_implemented();
136 #ifdef GRINS_USE_GRVY_TIMERS
137 this->_timer->EndTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_time_derivative");
148 #ifdef GRINS_USE_GRVY_TIMERS
149 this->_timer->BeginTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_constraint");
153 const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
156 const std::vector<libMesh::Real> &JxW =
157 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
160 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
161 context.get_element_fe(this->_flow_vars.p_var())->get_dphi();
163 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_flow_vars.p_var());
165 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
167 unsigned int n_qpoints = context.get_element_qrule().n_points();
169 for (
unsigned int qp=0; qp != n_qpoints; qp++)
171 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
172 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
174 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
175 context.interior_value( this->_flow_vars.v_var(), qp ) );
176 if( this->_dim == 3 )
178 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
182 libMesh::Real _mu_qp = this->_mu(context, qp);
184 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
186 libMesh::RealGradient RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
188 for (
unsigned int i=0; i != n_p_dofs; i++)
190 Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
193 if( compute_jacobian )
195 libmesh_not_implemented();
200 #ifdef GRINS_USE_GRVY_TIMERS
201 this->_timer->EndTimer(
"IncompressibleNavierStokesSPGSMStabilization::element_constraint");
212 #ifdef GRINS_USE_GRVY_TIMERS
213 this->_timer->BeginTimer(
"IncompressibleNavierStokesSPGSMStabilization::mass_residual");
217 const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
218 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
221 const std::vector<libMesh::Real> &JxW =
222 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
225 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
226 context.get_element_fe(this->_flow_vars.p_var())->get_dphi();
228 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
229 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
231 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
232 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
233 libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
235 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
237 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_flow_vars.p_var());
239 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
241 unsigned int n_qpoints = context.get_element_qrule().n_points();
243 for (
unsigned int qp=0; qp != n_qpoints; qp++)
245 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
246 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
248 libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u_var(), qp ),
249 context.fixed_interior_value( this->_flow_vars.v_var(), qp ) );
251 libMesh::Real _mu_qp = this->_mu(context, qp);
253 if( this->_dim == 3 )
255 U(2) = context.fixed_interior_value( this->_flow_vars.w_var(), qp );
258 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp,
false );
260 libMesh::RealGradient RM_t = this->_stab_helper.compute_res_momentum_transient( context, qp, this->_rho );
262 for (
unsigned int i=0; i != n_p_dofs; i++)
264 Fp(i) -= tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
267 for (
unsigned int i=0; i != n_u_dofs; i++)
269 Fu(i) -= tau_M*RM_t(0)*this->_rho*U*u_gradphi[i][qp]*JxW[qp];
271 Fv(i) -= tau_M*RM_t(1)*this->_rho*U*u_gradphi[i][qp]*JxW[qp];
273 if( this->_dim == 3 )
275 (*Fw)(i) -= tau_M*RM_t(2)*this->_rho*U*u_gradphi[i][qp]*JxW[qp];
279 if( compute_jacobian )
281 libmesh_not_implemented();
286 #ifdef GRINS_USE_GRVY_TIMERS
287 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.
virtual ~IncompressibleNavierStokesSPGSMStabilization()
virtual void read_input_options(const GetPot &input)
Read options from GetPot input file. By default, nothing is read.
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...