30 #include "grins_config.h"
37 #include "libmesh/getpot.h"
38 #include "libmesh/quadrature.h"
39 #include "libmesh/boundary_info.h"
58 if( input.have_variable(section) )
60 unsigned int n_vars = input.vector_variable_size(section);
62 for(
unsigned int v = 0; v < n_vars; v++ )
64 std::string name = input(section,
"DIE!",v);
66 if( name == std::string(
"k") )
73 <<
" Found " << name << std::endl
74 <<
" Acceptable values are: k" << std::endl;
88 #ifdef GRINS_USE_GRVY_TIMERS
89 this->_timer->BeginTimer(
"HeatTransfer::element_time_derivative");
93 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
94 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
102 const std::vector<libMesh::Real> &JxW =
103 context.get_element_fe(this->_temp_vars.T())->get_JxW();
106 const std::vector<std::vector<libMesh::Real> >& T_phi =
107 context.get_element_fe(this->_temp_vars.T())->get_phi();
110 const std::vector<std::vector<libMesh::Real> >& vel_phi =
111 context.get_element_fe(this->_flow_vars.u())->get_phi();
115 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
116 context.get_element_fe(this->_temp_vars.T())->get_dphi();
118 const std::vector<libMesh::Point>& u_qpoint =
119 context.get_element_fe(this->_flow_vars.u())->get_xyz();
121 libMesh::DenseSubMatrix<libMesh::Number> &KTT = context.get_elem_jacobian(this->_temp_vars.T(), this->_temp_vars.T());
123 libMesh::DenseSubMatrix<libMesh::Number> &KTu = context.get_elem_jacobian(this->_temp_vars.T(), this->_flow_vars.u());
124 libMesh::DenseSubMatrix<libMesh::Number> &KTv = context.get_elem_jacobian(this->_temp_vars.T(), this->_flow_vars.v());
125 libMesh::DenseSubMatrix<libMesh::Number>* KTw = NULL;
127 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
129 if( this->_dim == 3 )
131 KTw = &context.get_elem_jacobian(this->_temp_vars.T(), this->_flow_vars.w());
140 unsigned int n_qpoints = context.get_element_qrule().n_points();
142 for (
unsigned int qp=0; qp != n_qpoints; qp++)
145 libMesh::Number u, v;
146 u = context.interior_value(this->_flow_vars.u(), qp);
147 v = context.interior_value(this->_flow_vars.v(), qp);
149 libMesh::Gradient grad_T;
150 grad_T = context.interior_gradient(this->_temp_vars.T(), qp);
152 libMesh::NumberVectorValue U (u,v);
154 U(2) = context.interior_value(this->_flow_vars.w(), qp);
156 const libMesh::Number r = u_qpoint[qp](0);
158 libMesh::Real jac = JxW[qp];
161 libMesh::Real _k_qp = this->_k(context, qp);
169 for (
unsigned int i=0; i != n_T_dofs; i++)
172 (-this->_rho*this->_Cp*T_phi[i][qp]*(U*grad_T)
173 -_k_qp*(T_gradphi[i][qp]*grad_T) );
175 if (compute_jacobian)
177 for (
unsigned int j=0; j != n_T_dofs; j++)
182 KTT(i,j) += jac * context.get_elem_solution_derivative() *
183 (-this->_rho*this->_Cp*T_phi[i][qp]*(U*T_gradphi[j][qp])
184 -_k_qp*(T_gradphi[i][qp]*T_gradphi[j][qp]));
188 for (
unsigned int j=0; j != n_u_dofs; j++)
190 KTu(i,j) += jac * context.get_elem_solution_derivative()*(-this->_rho*this->_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(0)));
191 KTv(i,j) += jac * context.get_elem_solution_derivative()*(-this->_rho*this->_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(1)));
193 (*KTw)(i,j) += jac * context.get_elem_solution_derivative()*(-this->_rho*this->_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(2)));
201 #ifdef GRINS_USE_GRVY_TIMERS
202 this->_timer->EndTimer(
"HeatTransfer::element_time_derivative");
213 #ifdef GRINS_USE_GRVY_TIMERS
214 this->_timer->BeginTimer(
"HeatTransfer::mass_residual");
221 const std::vector<libMesh::Real> &JxW =
222 context.get_element_fe(this->_temp_vars.T())->get_JxW();
225 const std::vector<std::vector<libMesh::Real> >& phi =
226 context.get_element_fe(this->_temp_vars.T())->get_phi();
228 const std::vector<libMesh::Point>& u_qpoint =
229 context.get_element_fe(this->_flow_vars.u())->get_xyz();
232 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
235 libMesh::DenseSubVector<libMesh::Real> &F = context.get_elem_residual(this->_temp_vars.T());
237 libMesh::DenseSubMatrix<libMesh::Real> &M = context.get_elem_jacobian(this->_temp_vars.T(), this->_temp_vars.T());
239 unsigned int n_qpoints = context.get_element_qrule().n_points();
241 for (
unsigned int qp = 0; qp != n_qpoints; ++qp)
249 context.interior_rate(this->_temp_vars.T(), qp, T_dot);
251 const libMesh::Number r = u_qpoint[qp](0);
253 libMesh::Real jac = JxW[qp];
260 for (
unsigned int i = 0; i != n_T_dofs; ++i)
262 F(i) -= this->_rho*this->_Cp*T_dot*phi[i][qp]*jac;
264 if( compute_jacobian )
266 for (
unsigned int j=0; j != n_T_dofs; j++)
269 M(i,j) -= this->_rho*this->_Cp*phi[j][qp]*phi[i][qp]*jac * context.get_elem_solution_rate_derivative();
277 #ifdef GRINS_USE_GRVY_TIMERS
278 this->_timer->EndTimer(
"HeatTransfer::mass_residual");
287 const libMesh::Point& point,
288 libMesh::Real& value )
290 if( quantity_index == this->_k_index )
292 value = this->_k(point, context.get_time());
static bool is_axisymmetric()
static PhysicsName heat_transfer()
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
GRINS::ICHandlingBase * _ic_handler
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...
Physics class for Heat Transfer.
Base class for reading and handling initial conditions for physics classes.
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
Compute value of postprocessed quantities at libMesh::Point.
INSTANTIATE_HEAT_TRANSFER_SUBCLASS(HeatTransfer)
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for HeatTransfer.