27 #include "grins_config.h"
35 #include "libmesh/quadrature.h"
43 _stab_helper( physics_name+
"StabHelper", input )
55 context.get_element_fe(this->_press_var.p())->get_dphi();
57 context.get_element_fe(this->_flow_vars.u())->get_xyz();
58 context.get_element_fe(this->_flow_vars.u())->get_phi();
59 context.get_element_fe(this->_flow_vars.u())->get_dphi();
60 context.get_element_fe(this->_flow_vars.u())->get_d2phi();
67 (
bool compute_jacobian,
71 #ifdef GRINS_USE_GRVY_TIMERS
72 this->_timer->BeginTimer(
"VelocityDragAdjointStabilization::element_time_derivative");
75 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
78 const std::vector<libMesh::Real> &JxW = fe->get_JxW();
81 const std::vector<std::vector<libMesh::Real> >& u_phi = fe->get_phi();
83 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
84 context.get_element_fe(this->_flow_vars.u())->get_dphi();
86 const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
87 context.get_element_fe(this->_flow_vars.u())->get_d2phi();
89 const std::vector<libMesh::Point>& u_qpoint = fe->get_xyz();
92 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
95 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u());
96 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v());
97 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u());
98 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v());
100 libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
101 libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
102 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
103 libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
104 libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
106 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
107 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
108 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
110 if( this->_dim == 3 )
112 Kuw = &context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.w());
113 Kvw = &context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.w());
115 Kwu = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.u());
116 Kwv = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.v());
117 Kww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w());
118 Fw = &context.get_elem_residual(this->_flow_vars.w());
121 unsigned int n_qpoints = context.get_element_qrule().n_points();
123 for (
unsigned int qp=0; qp != n_qpoints; qp++)
125 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
126 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
128 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
129 context.interior_value( this->_flow_vars.v(), qp ) );
130 if( this->_dim == 3 )
132 U(2) = context.interior_value( this->_flow_vars.w(), qp );
136 libMesh::Real mu_qp = this->_mu(context, qp);
139 libMesh::Real d_tau_M_d_rho;
140 libMesh::Gradient d_tau_M_dU;
142 if (compute_jacobian)
143 this->_stab_helper.compute_tau_momentum_and_derivs
144 ( context, qp, g, G, this->_rho, U, mu_qp,
145 tau_M, d_tau_M_d_rho, d_tau_M_dU,
148 tau_M = this->_stab_helper.compute_tau_momentum
149 ( context, qp, g, G, this->_rho, U, mu_qp,
152 libMesh::NumberVectorValue F;
153 libMesh::NumberTensorValue dFdU;
154 libMesh::NumberTensorValue* dFdU_ptr =
155 compute_jacobian ? &dFdU : NULL;
156 if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
159 for (
unsigned int i=0; i != n_u_dofs; i++)
161 libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
162 mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
163 Fu(i) += tau_M*F(0)*test_func*JxW[qp];
165 Fv(i) += tau_M*F(1)*test_func*JxW[qp];
169 (*Fw)(i) += tau_M*F(2)*test_func*JxW[qp];
172 if (compute_jacobian)
174 const libMesh::Real sol_deriv =
175 context.get_elem_solution_derivative();
177 const libMesh::Real JxWxS = JxW[qp] * sol_deriv;
179 const libMesh::Real fixed_deriv =
180 context.get_fixed_solution_derivative();
182 const libMesh::Real JxWxF = JxW[qp] * fixed_deriv;
184 libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
186 for (
unsigned int j=0; j != n_u_dofs; ++j)
188 Kuu(i,j) += tau_M*F(0)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
189 Kuu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(0)*test_func*JxWxF;
190 Kuu(i,j) += tau_M*dFdU(0,0)*u_phi[j][qp]*test_func*JxWxS;
191 Kuv(i,j) += tau_M*F(0)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
192 Kuv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(0)*test_func*JxWxF;
193 Kuv(i,j) += tau_M*dFdU(0,1)*u_phi[j][qp]*test_func*JxWxS;
194 Kvu(i,j) += tau_M*F(1)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
195 Kvu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(1)*test_func*JxWxF;
196 Kvu(i,j) += tau_M*dFdU(1,0)*u_phi[j][qp]*test_func*JxWxS;
197 Kvv(i,j) += tau_M*F(1)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
198 Kvv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(1)*test_func*JxWxF;
199 Kvv(i,j) += tau_M*dFdU(1,1)*u_phi[j][qp]*test_func*JxWxS;
204 for (
unsigned int j=0; j != n_u_dofs; ++j)
206 (*Kuw)(i,j) += tau_M*F(0)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
207 (*Kuw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(0)*test_func*JxWxF;
208 (*Kuw)(i,j) += tau_M*dFdU(0,2)*u_phi[j][qp]*test_func*JxWxS;
209 (*Kvw)(i,j) += tau_M*F(1)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
210 (*Kvw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(1)*test_func*JxWxF;
211 (*Kvw)(i,j) += tau_M*dFdU(1,2)*u_phi[j][qp]*test_func*JxWxS;
212 (*Kwu)(i,j) += tau_M*F(2)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
213 (*Kwu)(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(2)*test_func*JxWxF;
214 (*Kwu)(i,j) += tau_M*dFdU(2,0)*u_phi[j][qp]*test_func*JxWxS;
215 (*Kwv)(i,j) += tau_M*F(2)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
216 (*Kwv)(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(2)*test_func*JxWxF;
217 (*Kwv)(i,j) += tau_M*dFdU(2,1)*u_phi[j][qp]*test_func*JxWxS;
218 (*Kww)(i,j) += tau_M*F(2)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
219 (*Kww)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(2)*test_func*JxWxF;
220 (*Kww)(i,j) += tau_M*dFdU(2,2)*u_phi[j][qp]*test_func*JxWxS;
229 #ifdef GRINS_USE_GRVY_TIMERS
230 this->_timer->EndTimer(
"VelocityDragAdjointStabilization::element_time_derivative");
241 #ifdef GRINS_USE_GRVY_TIMERS
242 this->_timer->BeginTimer(
"VelocityDragAdjointStabilization::element_constraint");
246 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
247 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
250 const std::vector<libMesh::Real> &JxW =
251 context.get_element_fe(this->_flow_vars.u())->get_JxW();
253 const std::vector<libMesh::Point>& u_qpoint =
254 context.get_element_fe(this->_flow_vars.u())->get_xyz();
256 const std::vector<std::vector<libMesh::Real> >& u_phi =
257 context.get_element_fe(this->_flow_vars.u())->get_phi();
259 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
260 context.get_element_fe(this->_press_var.p())->get_dphi();
262 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
264 libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
265 context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u());
266 libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
267 context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v());
268 libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
272 Kpw = &context.get_elem_jacobian
273 (this->_press_var.p(), this->_flow_vars.w());
282 unsigned int n_qpoints = context.get_element_qrule().n_points();
284 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
286 for (
unsigned int qp=0; qp != n_qpoints; qp++)
288 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
289 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
291 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
292 context.interior_value( this->_flow_vars.v(), qp ) );
293 if( this->_dim == 3 )
295 U(2) = context.interior_value( this->_flow_vars.w(), qp );
299 libMesh::Real mu_qp = this->_mu(context, qp);
302 libMesh::Real d_tau_M_d_rho;
303 libMesh::Gradient d_tau_M_dU;
305 if (compute_jacobian)
306 this->_stab_helper.compute_tau_momentum_and_derivs
307 ( context, qp, g, G, this->_rho, U, mu_qp,
308 tau_M, d_tau_M_d_rho, d_tau_M_dU,
311 tau_M = this->_stab_helper.compute_tau_momentum
312 ( context, qp, g, G, this->_rho, U, mu_qp,
315 libMesh::NumberVectorValue F;
316 libMesh::NumberTensorValue dFdU;
317 libMesh::NumberTensorValue* dFdU_ptr =
318 compute_jacobian ? &dFdU : NULL;
319 if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
325 for (
unsigned int i=0; i != n_p_dofs; i++)
327 Fp(i) += -tau_M*F*p_dphi[i][qp]*JxW[qp];
329 if (compute_jacobian)
331 const libMesh::Real sol_deriv =
332 context.get_elem_solution_derivative();
334 const libMesh::Real JxWxS = JxW[qp] * sol_deriv;
336 const libMesh::Real fixed_deriv =
337 context.get_fixed_solution_derivative();
339 const libMesh::Real JxWxF = JxW[qp] * fixed_deriv;
341 for (
unsigned int j=0; j != n_u_dofs; ++j)
343 Kpu(i,j) += -d_tau_M_dU(0)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
344 Kpv(i,j) += -d_tau_M_dU(1)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
345 for (
unsigned int d=0; d != 3; ++d)
347 Kpu(i,j) += -tau_M*dFdU(d,0)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
348 Kpv(i,j) += -tau_M*dFdU(d,1)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
351 if( this->_dim == 3 )
352 for (
unsigned int j=0; j != n_u_dofs; ++j)
354 (*Kpw)(i,j) += -d_tau_M_dU(2)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
355 for (
unsigned int d=0; d != 3; ++d)
357 (*Kpw)(i,j) += -tau_M*dFdU(d,2)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
364 #ifdef GRINS_USE_GRVY_TIMERS
365 this->_timer->EndTimer(
"VelocityDragAdjointStabilization::element_constraint");
INSTANTIATE_INC_NS_SUBCLASS(VelocityDragAdjointStabilization)
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
~VelocityDragAdjointStabilization()
VelocityDragAdjointStabilization()
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.