27 #include "grins_config.h"
35 #include "libmesh/quadrature.h"
45 _stab_helper( physics_name+
"StabHelper", input )
60 context.get_element_fe(this->_flow_vars.p_var())->get_dphi();
62 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
63 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
64 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
65 context.get_element_fe(this->_flow_vars.u_var())->get_d2phi();
72 (
bool compute_jacobian,
76 #ifdef GRINS_USE_GRVY_TIMERS
77 this->_timer->BeginTimer(
"VelocityDragAdjointStabilization::element_time_derivative");
80 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
83 const std::vector<libMesh::Real> &JxW = fe->get_JxW();
86 const std::vector<std::vector<libMesh::Real> >& u_phi = fe->get_phi();
88 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
89 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
91 const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
92 context.get_element_fe(this->_flow_vars.u_var())->get_d2phi();
94 const std::vector<libMesh::Point>& u_qpoint = fe->get_xyz();
97 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
100 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
101 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.v_var());
102 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.u_var());
103 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
105 libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
106 libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
107 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
108 libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
109 libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
111 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
112 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
113 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
115 if( this->_dim == 3 )
117 Kuw = &context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.w_var());
118 Kvw = &context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.w_var());
120 Kwu = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.u_var());
121 Kwv = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.v_var());
122 Kww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
123 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
126 unsigned int n_qpoints = context.get_element_qrule().n_points();
128 for (
unsigned int qp=0; qp != n_qpoints; qp++)
130 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
131 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
133 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
134 context.interior_value( this->_flow_vars.v_var(), qp ) );
135 if( this->_dim == 3 )
137 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
141 libMesh::Real mu_qp = this->_mu(context, qp);
144 libMesh::Real d_tau_M_d_rho;
145 libMesh::Gradient d_tau_M_dU;
147 if (compute_jacobian)
148 this->_stab_helper.compute_tau_momentum_and_derivs
149 ( context, qp, g, G, this->_rho, U, mu_qp,
150 tau_M, d_tau_M_d_rho, d_tau_M_dU,
153 tau_M = this->_stab_helper.compute_tau_momentum
154 ( context, qp, g, G, this->_rho, U, mu_qp,
157 libMesh::NumberVectorValue F;
158 libMesh::NumberTensorValue dFdU;
159 libMesh::NumberTensorValue* dFdU_ptr =
160 compute_jacobian ? &dFdU : NULL;
161 if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
164 for (
unsigned int i=0; i != n_u_dofs; i++)
166 libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
167 mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
168 Fu(i) += tau_M*F(0)*test_func*JxW[qp];
170 Fv(i) += tau_M*F(1)*test_func*JxW[qp];
174 (*Fw)(i) += tau_M*F(2)*test_func*JxW[qp];
177 if (compute_jacobian)
179 const libMesh::Real sol_deriv =
180 context.get_elem_solution_derivative();
182 const libMesh::Real JxWxS = JxW[qp] * sol_deriv;
184 const libMesh::Real fixed_deriv =
185 context.get_fixed_solution_derivative();
187 const libMesh::Real JxWxF = JxW[qp] * fixed_deriv;
189 libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
191 for (
unsigned int j=0; j != n_u_dofs; ++j)
193 Kuu(i,j) += tau_M*F(0)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
194 Kuu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(0)*test_func*JxWxF;
195 Kuu(i,j) += tau_M*dFdU(0,0)*u_phi[j][qp]*test_func*JxWxS;
196 Kuv(i,j) += tau_M*F(0)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
197 Kuv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(0)*test_func*JxWxF;
198 Kuv(i,j) += tau_M*dFdU(0,1)*u_phi[j][qp]*test_func*JxWxS;
199 Kvu(i,j) += tau_M*F(1)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
200 Kvu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(1)*test_func*JxWxF;
201 Kvu(i,j) += tau_M*dFdU(1,0)*u_phi[j][qp]*test_func*JxWxS;
202 Kvv(i,j) += tau_M*F(1)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
203 Kvv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(1)*test_func*JxWxF;
204 Kvv(i,j) += tau_M*dFdU(1,1)*u_phi[j][qp]*test_func*JxWxS;
209 for (
unsigned int j=0; j != n_u_dofs; ++j)
211 (*Kuw)(i,j) += tau_M*F(0)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
212 (*Kuw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(0)*test_func*JxWxF;
213 (*Kuw)(i,j) += tau_M*dFdU(0,2)*u_phi[j][qp]*test_func*JxWxS;
214 (*Kvw)(i,j) += tau_M*F(1)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
215 (*Kvw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(1)*test_func*JxWxF;
216 (*Kvw)(i,j) += tau_M*dFdU(1,2)*u_phi[j][qp]*test_func*JxWxS;
217 (*Kwu)(i,j) += tau_M*F(2)*d_test_func_dU(0)*u_phi[j][qp]*JxWxS;
218 (*Kwu)(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(2)*test_func*JxWxF;
219 (*Kwu)(i,j) += tau_M*dFdU(2,0)*u_phi[j][qp]*test_func*JxWxS;
220 (*Kwv)(i,j) += tau_M*F(2)*d_test_func_dU(1)*u_phi[j][qp]*JxWxS;
221 (*Kwv)(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(2)*test_func*JxWxF;
222 (*Kwv)(i,j) += tau_M*dFdU(2,1)*u_phi[j][qp]*test_func*JxWxS;
223 (*Kww)(i,j) += tau_M*F(2)*d_test_func_dU(2)*u_phi[j][qp]*JxWxS;
224 (*Kww)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(2)*test_func*JxWxF;
225 (*Kww)(i,j) += tau_M*dFdU(2,2)*u_phi[j][qp]*test_func*JxWxS;
234 #ifdef GRINS_USE_GRVY_TIMERS
235 this->_timer->EndTimer(
"VelocityDragAdjointStabilization::element_time_derivative");
246 #ifdef GRINS_USE_GRVY_TIMERS
247 this->_timer->BeginTimer(
"VelocityDragAdjointStabilization::element_constraint");
251 const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
252 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
255 const std::vector<libMesh::Real> &JxW =
256 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
258 const std::vector<libMesh::Point>& u_qpoint =
259 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
261 const std::vector<std::vector<libMesh::Real> >& u_phi =
262 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
264 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
265 context.get_element_fe(this->_flow_vars.p_var())->get_dphi();
267 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_flow_vars.p_var());
269 libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
270 context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.u_var());
271 libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
272 context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.v_var());
273 libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
277 Kpw = &context.get_elem_jacobian
278 (this->_flow_vars.p_var(), this->_flow_vars.w_var());
287 unsigned int n_qpoints = context.get_element_qrule().n_points();
289 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
291 for (
unsigned int qp=0; qp != n_qpoints; qp++)
293 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
294 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
296 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
297 context.interior_value( this->_flow_vars.v_var(), qp ) );
298 if( this->_dim == 3 )
300 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
304 libMesh::Real mu_qp = this->_mu(context, qp);
307 libMesh::Real d_tau_M_d_rho;
308 libMesh::Gradient d_tau_M_dU;
310 if (compute_jacobian)
311 this->_stab_helper.compute_tau_momentum_and_derivs
312 ( context, qp, g, G, this->_rho, U, mu_qp,
313 tau_M, d_tau_M_d_rho, d_tau_M_dU,
316 tau_M = this->_stab_helper.compute_tau_momentum
317 ( context, qp, g, G, this->_rho, U, mu_qp,
320 libMesh::NumberVectorValue F;
321 libMesh::NumberTensorValue dFdU;
322 libMesh::NumberTensorValue* dFdU_ptr =
323 compute_jacobian ? &dFdU : NULL;
324 if (!this->compute_force(u_qpoint[qp], context.time, U, F, dFdU_ptr))
330 for (
unsigned int i=0; i != n_p_dofs; i++)
332 Fp(i) += -tau_M*F*p_dphi[i][qp]*JxW[qp];
334 if (compute_jacobian)
336 const libMesh::Real sol_deriv =
337 context.get_elem_solution_derivative();
339 const libMesh::Real JxWxS = JxW[qp] * sol_deriv;
341 const libMesh::Real fixed_deriv =
342 context.get_fixed_solution_derivative();
344 const libMesh::Real JxWxF = JxW[qp] * fixed_deriv;
346 for (
unsigned int j=0; j != n_u_dofs; ++j)
348 Kpu(i,j) += -d_tau_M_dU(0)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
349 Kpv(i,j) += -d_tau_M_dU(1)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
350 for (
unsigned int d=0; d != 3; ++d)
352 Kpu(i,j) += -tau_M*dFdU(d,0)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
353 Kpv(i,j) += -tau_M*dFdU(d,1)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
356 if( this->_dim == 3 )
357 for (
unsigned int j=0; j != n_u_dofs; ++j)
359 (*Kpw)(i,j) += -d_tau_M_dU(2)*u_phi[j][qp]*F*p_dphi[i][qp]*JxWxF;
360 for (
unsigned int d=0; d != 3; ++d)
362 (*Kpw)(i,j) += -tau_M*dFdU(d,2)*u_phi[j][qp]*p_dphi[i][qp](d)*JxWxS;
369 #ifdef GRINS_USE_GRVY_TIMERS
370 this->_timer->EndTimer(
"VelocityDragAdjointStabilization::element_constraint");
virtual void set_parameter(libMesh::Number ¶m_variable, const GetPot &input, const std::string ¶m_name, libMesh::Number param_default)
Each subclass can simultaneously read a parameter value from.
const PhysicsName incompressible_navier_stokes
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.