32 #include "libmesh/quadrature.h"
56 #ifdef GRINS_USE_GRVY_TIMERS
57 this->_timer->BeginTimer(
"HeatTransferAdjointStabilization::element_time_derivative");
61 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T_var()).size();
62 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
65 const std::vector<libMesh::Real> &JxW =
66 context.get_element_fe(this->_temp_vars.T_var())->get_JxW();
68 const std::vector<std::vector<libMesh::Real> >& u_phi =
69 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
71 const std::vector<std::vector<libMesh::Real> >& T_phi =
72 context.get_element_fe(this->_temp_vars.T_var())->get_phi();
74 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
75 context.get_element_fe(this->_temp_vars.T_var())->get_dphi();
77 const std::vector<std::vector<libMesh::RealTensor> >& T_hessphi =
78 context.get_element_fe(this->_temp_vars.T_var())->get_d2phi();
95 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T_var());
96 libMesh::DenseSubMatrix<libMesh::Number> &KTT =
97 context.get_elem_jacobian(this->_temp_vars.T_var(), this->_temp_vars.T_var());
98 libMesh::DenseSubMatrix<libMesh::Number> &KTu =
99 context.get_elem_jacobian(this->_temp_vars.T_var(), this->_flow_vars.u_var());
100 libMesh::DenseSubMatrix<libMesh::Number> &KTv =
101 context.get_elem_jacobian(this->_temp_vars.T_var(), this->_flow_vars.v_var());
102 libMesh::DenseSubMatrix<libMesh::Number> *KTw = NULL;
106 KTw = &context.get_elem_jacobian
107 (this->_temp_vars.T_var(), this->_flow_vars.w_var());
110 unsigned int n_qpoints = context.get_element_qrule().n_points();
112 for (
unsigned int qp=0; qp != n_qpoints; qp++)
114 libMesh::FEBase* fe = context.get_element_fe(this->_temp_vars.T_var());
116 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
117 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
119 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
120 context.interior_value( this->_flow_vars.v_var(), qp ) );
121 if( this->_dim == 3 )
123 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
128 libMesh::Real tau_E, RE_s;
129 libMesh::Real d_tau_E_dT_rho, d_RE_s_dT;
130 libMesh::Gradient d_tau_E_dU, d_RE_s_dgradT, d_RE_s_dU;
131 libMesh::Tensor d_RE_s_dhessT;
134 libMesh::Real _k_qp = this->_k(context, qp);
136 if (compute_jacobian)
138 this->_stab_helper.compute_tau_energy_and_derivs
139 ( context, G, this->_rho, this->_Cp, _k_qp, U,
140 tau_E, d_tau_E_dT_rho, d_tau_E_dU, this->_is_steady );
141 this->_stab_helper.compute_res_energy_steady_and_derivs
142 ( context, qp, this->_rho, this->_Cp, _k_qp,
143 RE_s, d_RE_s_dT, d_RE_s_dgradT, d_RE_s_dhessT,
148 tau_E = this->_stab_helper.compute_tau_energy
149 ( context, G, this->_rho, this->_Cp, _k_qp, U, this->_is_steady );
151 RE_s = this->_stab_helper.compute_res_energy_steady
152 ( context, qp, this->_rho, this->_Cp, _k_qp );
167 for (
unsigned int i=0; i != n_T_dofs; i++)
169 FT(i) += -tau_E*RE_s*( this->_rho*this->_Cp*U*T_gradphi[i][qp]
170 + _k_qp*(T_hessphi[i][qp](0,0) + T_hessphi[i][qp](1,1) + T_hessphi[i][qp](2,2))
172 if (compute_jacobian)
174 for (
unsigned int j=0; j != n_T_dofs; ++j)
177 (d_RE_s_dT*T_phi[j][qp] +
178 d_RE_s_dgradT*T_gradphi[j][qp] +
179 d_RE_s_dhessT.contract(T_hessphi[j][qp])
181 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
182 + _k_qp*(T_hessphi[i][qp](0,0) +
183 T_hessphi[i][qp](1,1) +
184 T_hessphi[i][qp](2,2))
186 * context.get_fixed_solution_derivative();
188 for (
unsigned int j=0; j != n_u_dofs; ++j)
190 KTu(i,j) += -tau_E*RE_s*
191 ( this->_rho*this->_Cp*u_phi[j][qp]*T_gradphi[i][qp](0) )*JxW[qp]
192 * context.get_fixed_solution_derivative();
194 -(tau_E*d_RE_s_dU(0)+d_tau_E_dU(0)*RE_s)*u_phi[j][qp]*
195 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
196 + _k_qp*(T_hessphi[i][qp](0,0) +
197 T_hessphi[i][qp](1,1) +
198 T_hessphi[i][qp](2,2))
200 * context.get_fixed_solution_derivative();
201 KTv(i,j) += -tau_E*RE_s*
202 ( this->_rho*this->_Cp*u_phi[j][qp]*T_gradphi[i][qp](1) )*JxW[qp]
203 * context.get_fixed_solution_derivative();
205 -(tau_E*d_RE_s_dU(1)+d_tau_E_dU(1)*RE_s)*u_phi[j][qp]*
206 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
207 + _k_qp*(T_hessphi[i][qp](0,0) +
208 T_hessphi[i][qp](1,1) +
209 T_hessphi[i][qp](2,2))
211 * context.get_fixed_solution_derivative();
215 for (
unsigned int j=0; j != n_u_dofs; ++j)
217 (*KTw)(i,j) += -tau_E*RE_s*
218 ( this->_rho*this->_Cp*u_phi[j][qp]*T_gradphi[i][qp](2) )*JxW[qp]
219 * context.get_fixed_solution_derivative();
221 -(tau_E*d_RE_s_dU(2)+d_tau_E_dU(2)*RE_s)*u_phi[j][qp]*
222 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
223 + _k_qp*(T_hessphi[i][qp](0,0) +
224 T_hessphi[i][qp](1,1) +
225 T_hessphi[i][qp](2,2))
227 * context.get_fixed_solution_derivative();
234 #ifdef GRINS_USE_GRVY_TIMERS
235 this->_timer->EndTimer(
"HeatTransferAdjointStabilization::element_time_derivative");
245 #ifdef GRINS_USE_GRVY_TIMERS
246 this->_timer->BeginTimer(
"HeatTransferAdjointStabilization::mass_residual");
250 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T_var()).size();
251 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
254 const std::vector<libMesh::Real> &JxW =
255 context.get_element_fe(this->_temp_vars.T_var())->get_JxW();
257 const std::vector<std::vector<libMesh::Real> >& u_phi =
258 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
260 const std::vector<std::vector<libMesh::Real> >& T_phi =
261 context.get_element_fe(this->_temp_vars.T_var())->get_phi();
263 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
264 context.get_element_fe(this->_temp_vars.T_var())->get_dphi();
266 const std::vector<std::vector<libMesh::RealTensor> >& T_hessphi =
267 context.get_element_fe(this->_temp_vars.T_var())->get_d2phi();
284 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T_var());
285 libMesh::DenseSubMatrix<libMesh::Number> &KTT =
286 context.get_elem_jacobian(this->_temp_vars.T_var(), this->_temp_vars.T_var());
287 libMesh::DenseSubMatrix<libMesh::Number> &KTu =
288 context.get_elem_jacobian(this->_temp_vars.T_var(), this->_flow_vars.u_var());
289 libMesh::DenseSubMatrix<libMesh::Number> &KTv =
290 context.get_elem_jacobian(this->_temp_vars.T_var(), this->_flow_vars.v_var());
291 libMesh::DenseSubMatrix<libMesh::Number> *KTw = NULL;
295 KTw = &context.get_elem_jacobian
296 (this->_temp_vars.T_var(), this->_flow_vars.w_var());
299 unsigned int n_qpoints = context.get_element_qrule().n_points();
301 for (
unsigned int qp=0; qp != n_qpoints; qp++)
303 libMesh::FEBase* fe = context.get_element_fe(this->_temp_vars.T_var());
305 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
306 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
308 libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u_var(), qp ),
309 context.fixed_interior_value( this->_flow_vars.v_var(), qp ) );
310 if( this->_dim == 3 )
311 U(2) = context.fixed_interior_value( this->_flow_vars.w_var(), qp );
315 libMesh::Real tau_E, RE_t;
316 libMesh::Real d_tau_E_d_rho, d_RE_t_dT;
317 libMesh::Gradient d_tau_E_dU;
320 libMesh::Real _k_qp = this->_k(context, qp);
322 if (compute_jacobian)
324 this->_stab_helper.compute_tau_energy_and_derivs
325 ( context, G, this->_rho, this->_Cp, _k_qp, U,
326 tau_E, d_tau_E_d_rho, d_tau_E_dU,
false );
327 this->_stab_helper.compute_res_energy_transient_and_derivs
328 ( context, qp, this->_rho, this->_Cp,
333 tau_E = this->_stab_helper.compute_tau_energy
334 ( context, G, this->_rho, this->_Cp, _k_qp, U,
false );
336 RE_t = this->_stab_helper.compute_res_energy_transient
337 ( context, qp, this->_rho, this->_Cp );
353 for (
unsigned int i=0; i != n_T_dofs; i++)
355 FT(i) -= tau_E*RE_t*( this->_rho*this->_Cp*U*T_gradphi[i][qp]
356 + _k_qp*(T_hessphi[i][qp](0,0) + T_hessphi[i][qp](1,1) + T_hessphi[i][qp](2,2))
358 if (compute_jacobian)
360 const libMesh::Real fixed_deriv =
361 context.get_fixed_solution_derivative();
363 for (
unsigned int j=0; j != n_T_dofs; ++j)
366 (tau_E*d_RE_t_dT)*T_phi[j][qp]*
367 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
368 + _k_qp*(T_hessphi[i][qp](0,0) +
369 T_hessphi[i][qp](1,1) +
370 T_hessphi[i][qp](2,2))
373 for (
unsigned int j=0; j != n_u_dofs; ++j)
376 d_tau_E_dU(0)*u_phi[j][qp]*RE_t*
377 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
378 + _k_qp*(T_hessphi[i][qp](0,0) +
379 T_hessphi[i][qp](1,1) +
380 T_hessphi[i][qp](2,2))
381 )*fixed_deriv*JxW[qp];
384 ( this->_rho*this->_Cp*u_phi[j][qp]*T_gradphi[i][qp](0)*fixed_deriv
387 d_tau_E_dU(1)*u_phi[j][qp]*RE_t*
388 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
389 + _k_qp*(T_hessphi[i][qp](0,0) +
390 T_hessphi[i][qp](1,1) +
391 T_hessphi[i][qp](2,2))
392 )*fixed_deriv*JxW[qp];
395 ( this->_rho*this->_Cp*u_phi[j][qp]*T_gradphi[i][qp](1)*fixed_deriv
400 for (
unsigned int j=0; j != n_u_dofs; ++j)
403 d_tau_E_dU(2)*u_phi[j][qp]*RE_t*
404 ( this->_rho*this->_Cp*U*T_gradphi[i][qp]
405 + _k_qp*(T_hessphi[i][qp](0,0) +
406 T_hessphi[i][qp](1,1) +
407 T_hessphi[i][qp](2,2))
408 )*fixed_deriv*JxW[qp];
411 ( this->_rho*this->_Cp*u_phi[j][qp]*T_gradphi[i][qp](2)*fixed_deriv
420 #ifdef GRINS_USE_GRVY_TIMERS
421 this->_timer->EndTimer(
"HeatTransferAdjointStabilization::mass_residual");
INSTANTIATE_HEAT_TRANSFER_SUBCLASS(HeatTransferAdjointStabilization)
HeatTransferAdjointStabilization()
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...
virtual ~HeatTransferAdjointStabilization()