27 #include "grins_config.h"
35 #include "libmesh/getpot.h"
36 #include "libmesh/quadrature.h"
45 _stab_helper( physics_name+
"StabHelper", input )
58 context.get_element_fe(this->_flow_vars.p_var())->get_dphi();
60 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
61 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
62 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
63 context.get_element_fe(this->_flow_vars.u_var())->get_d2phi();
73 #ifdef GRINS_USE_GRVY_TIMERS
74 this->_timer->BeginTimer(
"ParsedVelocitySourceAdjointStabilization::element_time_derivative");
78 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
81 const std::vector<libMesh::Real> &JxW =
82 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
84 const std::vector<libMesh::Point>& u_qpoint =
85 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
87 const std::vector<std::vector<libMesh::Real> >& u_phi =
88 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
90 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
91 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
93 const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
94 context.get_element_fe(this->_flow_vars.u_var())->get_d2phi();
97 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
98 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
99 libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
101 libMesh::DenseSubMatrix<libMesh::Number> &Kuu =
102 context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
103 libMesh::DenseSubMatrix<libMesh::Number> &Kuv =
104 context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.v_var());
105 libMesh::DenseSubMatrix<libMesh::Number> &Kvu =
106 context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.u_var());
107 libMesh::DenseSubMatrix<libMesh::Number> &Kvv =
108 context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
110 libMesh::DenseSubMatrix<libMesh::Number> *Kuw = NULL;
111 libMesh::DenseSubMatrix<libMesh::Number> *Kvw = NULL;
112 libMesh::DenseSubMatrix<libMesh::Number> *Kwu = NULL;
113 libMesh::DenseSubMatrix<libMesh::Number> *Kwv = NULL;
114 libMesh::DenseSubMatrix<libMesh::Number> *Kww = NULL;
118 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
119 Kuw = &context.get_elem_jacobian
120 (this->_flow_vars.u_var(), this->_flow_vars.w_var());
121 Kvw = &context.get_elem_jacobian
122 (this->_flow_vars.v_var(), this->_flow_vars.w_var());
123 Kwu = &context.get_elem_jacobian
124 (this->_flow_vars.w_var(), this->_flow_vars.u_var());
125 Kwv = &context.get_elem_jacobian
126 (this->_flow_vars.w_var(), this->_flow_vars.v_var());
127 Kww = &context.get_elem_jacobian
128 (this->_flow_vars.w_var(), this->_flow_vars.w_var());
137 unsigned int n_qpoints = context.get_element_qrule().n_points();
139 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
141 for (
unsigned int qp=0; qp != n_qpoints; qp++)
143 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
144 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
146 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
147 context.interior_value( this->_flow_vars.v_var(), qp ) );
148 if( this->_dim == 3 )
150 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
154 libMesh::Real mu_qp = this->_mu(context, qp);
157 libMesh::Real d_tau_M_d_rho;
158 libMesh::Gradient d_tau_M_dU;
160 if (compute_jacobian)
161 this->_stab_helper.compute_tau_momentum_and_derivs
162 ( context, qp, g, G, this->_rho, U, mu_qp,
163 tau_M, d_tau_M_d_rho, d_tau_M_dU,
166 tau_M = this->_stab_helper.compute_tau_momentum
167 ( context, qp, g, G, this->_rho, U, mu_qp,
170 libMesh::NumberVectorValue F;
171 libMesh::NumberTensorValue dFdU;
172 libMesh::NumberTensorValue* dFdU_ptr =
173 compute_jacobian ? &dFdU : NULL;
174 if (!this->compute_force(u_qpoint[qp], context.time, context, F, dFdU_ptr))
177 for (
unsigned int i=0; i != n_u_dofs; i++)
179 libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
180 mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
181 Fu(i) += tau_M*F(0)*test_func*JxW[qp];
183 Fv(i) += tau_M*F(1)*test_func*JxW[qp];
187 (*Fw)(i) += tau_M*F(2)*test_func*JxW[qp];
190 if (compute_jacobian)
192 libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
194 for (
unsigned int j=0; j != n_u_dofs; ++j)
196 Kuu(i,j) += tau_M*F(0)*d_test_func_dU(0)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
197 Kuu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(0)*test_func*JxW[qp]*context.get_elem_solution_derivative();
198 Kuu(i,j) += tau_M*dFdU(0,0)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
199 Kuv(i,j) += tau_M*F(0)*d_test_func_dU(1)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
200 Kuv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(0)*test_func*JxW[qp]*context.get_elem_solution_derivative();
201 Kuv(i,j) += tau_M*dFdU(0,1)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
202 Kvu(i,j) += tau_M*F(1)*d_test_func_dU(0)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
203 Kvu(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(1)*test_func*JxW[qp]*context.get_elem_solution_derivative();
204 Kvu(i,j) += tau_M*dFdU(1,0)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
205 Kvv(i,j) += tau_M*F(1)*d_test_func_dU(1)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
206 Kvv(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(1)*test_func*JxW[qp]*context.get_elem_solution_derivative();
207 Kvv(i,j) += tau_M*dFdU(1,1)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
212 for (
unsigned int j=0; j != n_u_dofs; ++j)
214 (*Kuw)(i,j) += tau_M*F(0)*d_test_func_dU(2)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
215 (*Kuw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(0)*test_func*JxW[qp]*context.get_elem_solution_derivative();
216 (*Kuw)(i,j) += tau_M*dFdU(0,2)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
217 (*Kvw)(i,j) += tau_M*F(1)*d_test_func_dU(2)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
218 (*Kvw)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(1)*test_func*JxW[qp]*context.get_elem_solution_derivative();
219 (*Kvw)(i,j) += tau_M*dFdU(1,2)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
220 (*Kwu)(i,j) += tau_M*F(2)*d_test_func_dU(0)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
221 (*Kwu)(i,j) += d_tau_M_dU(0)*u_phi[j][qp]*F(2)*test_func*JxW[qp]*context.get_elem_solution_derivative();
222 (*Kwu)(i,j) += tau_M*dFdU(2,0)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
223 (*Kwv)(i,j) += tau_M*F(2)*d_test_func_dU(1)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
224 (*Kwv)(i,j) += d_tau_M_dU(1)*u_phi[j][qp]*F(2)*test_func*JxW[qp]*context.get_elem_solution_derivative();
225 (*Kwv)(i,j) += tau_M*dFdU(2,1)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
226 (*Kww)(i,j) += tau_M*F(2)*d_test_func_dU(2)*u_phi[j][qp]*JxW[qp]*context.get_elem_solution_derivative();
227 (*Kww)(i,j) += d_tau_M_dU(2)*u_phi[j][qp]*F(2)*test_func*JxW[qp]*context.get_elem_solution_derivative();
228 (*Kww)(i,j) += tau_M*dFdU(2,2)*u_phi[j][qp]*test_func*JxW[qp]*context.get_elem_solution_derivative();
237 #ifdef GRINS_USE_GRVY_TIMERS
238 this->_timer->EndTimer(
"BoussinesqBuoyancyAdjointStabilization::element_time_derivative");
249 #ifdef GRINS_USE_GRVY_TIMERS
250 this->_timer->BeginTimer(
"ParsedVelocitySourceAdjointStabilization::element_constraint");
254 const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
255 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
258 const std::vector<libMesh::Real> &JxW =
259 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
261 const std::vector<libMesh::Point>& u_qpoint =
262 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
264 const std::vector<std::vector<libMesh::Real> >& u_phi =
265 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
267 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
268 context.get_element_fe(this->_flow_vars.p_var())->get_dphi();
270 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_flow_vars.p_var());
272 libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
273 context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.u_var());
274 libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
275 context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.v_var());
276 libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
280 Kpw = &context.get_elem_jacobian
281 (this->_flow_vars.p_var(), this->_flow_vars.w_var());
290 unsigned int n_qpoints = context.get_element_qrule().n_points();
292 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
294 for (
unsigned int qp=0; qp != n_qpoints; qp++)
296 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
297 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
299 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
300 context.interior_value( this->_flow_vars.v_var(), qp ) );
301 if( this->_dim == 3 )
303 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
307 libMesh::Real mu_qp = this->_mu(context, qp);
310 libMesh::Real d_tau_M_d_rho;
311 libMesh::Gradient d_tau_M_dU;
313 if (compute_jacobian)
314 this->_stab_helper.compute_tau_momentum_and_derivs
315 ( context, qp, g, G, this->_rho, U, mu_qp,
316 tau_M, d_tau_M_d_rho, d_tau_M_dU,
319 tau_M = this->_stab_helper.compute_tau_momentum
320 ( context, qp, g, G, this->_rho, U, mu_qp,
323 libMesh::NumberVectorValue F;
324 libMesh::NumberTensorValue dFdU;
325 libMesh::NumberTensorValue* dFdU_ptr =
326 compute_jacobian ? &dFdU : NULL;
327 if (!this->compute_force(u_qpoint[qp], context.time, context, F, dFdU_ptr))
333 for (
unsigned int i=0; i != n_p_dofs; i++)
335 Fp(i) += -tau_M*F*p_dphi[i][qp]*JxW[qp];
337 if (compute_jacobian)
339 for (
unsigned int j=0; j != n_u_dofs; ++j)
341 Kpu(i,j) += -d_tau_M_dU(0)*u_phi[j][qp]*F*p_dphi[i][qp]*JxW[qp]*context.get_elem_solution_derivative();
342 Kpv(i,j) += -d_tau_M_dU(1)*u_phi[j][qp]*F*p_dphi[i][qp]*JxW[qp]*context.get_elem_solution_derivative();
343 for (
unsigned int d=0; d != 3; ++d)
345 Kpu(i,j) += -tau_M*dFdU(d,0)*u_phi[j][qp]*p_dphi[i][qp](d)*JxW[qp]*context.get_elem_solution_derivative();
346 Kpv(i,j) += -tau_M*dFdU(d,1)*u_phi[j][qp]*p_dphi[i][qp](d)*JxW[qp]*context.get_elem_solution_derivative();
349 if( this->_dim == 3 )
350 for (
unsigned int j=0; j != n_u_dofs; ++j)
352 (*Kpw)(i,j) += -d_tau_M_dU(2)*u_phi[j][qp]*F*p_dphi[i][qp]*JxW[qp]*context.get_elem_solution_derivative();
353 for (
unsigned int d=0; d != 3; ++d)
355 (*Kpw)(i,j) += -tau_M*dFdU(d,2)*u_phi[j][qp]*p_dphi[i][qp](d)*JxW[qp]*context.get_elem_solution_derivative();
362 #ifdef GRINS_USE_GRVY_TIMERS
363 this->_timer->EndTimer(
"ParsedVelocitySourceAdjointStabilization::element_constraint");
INSTANTIATE_INC_NS_SUBCLASS(ParsedVelocitySourceAdjointStabilization)
~ParsedVelocitySourceAdjointStabilization()
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
ParsedVelocitySourceAdjointStabilization()
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.