30 #include "grins_config.h"
38 #include "libmesh/getpot.h"
39 #include "libmesh/quadrature.h"
40 #include "libmesh/boundary_info.h"
81 std::string section =
"Physics/"+
heat_transfer+
"/output_vars";
83 if( input.have_variable(section) )
85 unsigned int n_vars = input.vector_variable_size(section);
87 for(
unsigned int v = 0; v < n_vars; v++ )
89 std::string name = input(section,
"DIE!",v);
91 if( name == std::string(
"k") )
97 std::cerr <<
"Error: Invalid output_vars value for "+
heat_transfer << std::endl
98 <<
" Found " << name << std::endl
99 <<
" Acceptable values are: k" << std::endl;
113 #ifdef GRINS_USE_GRVY_TIMERS
114 this->_timer->BeginTimer(
"HeatTransfer::element_time_derivative");
118 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T_var()).size();
119 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
127 const std::vector<libMesh::Real> &JxW =
128 context.get_element_fe(this->_temp_vars.T_var())->get_JxW();
131 const std::vector<std::vector<libMesh::Real> >& T_phi =
132 context.get_element_fe(this->_temp_vars.T_var())->get_phi();
135 const std::vector<std::vector<libMesh::Real> >& vel_phi =
136 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
140 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
141 context.get_element_fe(this->_temp_vars.T_var())->get_dphi();
143 const std::vector<libMesh::Point>& u_qpoint =
144 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
146 libMesh::DenseSubMatrix<libMesh::Number> &KTT = context.get_elem_jacobian(this->_temp_vars.T_var(), this->_temp_vars.T_var());
148 libMesh::DenseSubMatrix<libMesh::Number> &KTu = context.get_elem_jacobian(this->_temp_vars.T_var(), this->_flow_vars.u_var());
149 libMesh::DenseSubMatrix<libMesh::Number> &KTv = context.get_elem_jacobian(this->_temp_vars.T_var(), this->_flow_vars.v_var());
150 libMesh::DenseSubMatrix<libMesh::Number>* KTw = NULL;
152 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T_var());
154 if( this->_dim == 3 )
156 KTw = &context.get_elem_jacobian(this->_temp_vars.T_var(), this->_flow_vars.w_var());
165 unsigned int n_qpoints = context.get_element_qrule().n_points();
167 for (
unsigned int qp=0; qp != n_qpoints; qp++)
170 libMesh::Number u, v;
171 u = context.interior_value(this->_flow_vars.u_var(), qp);
172 v = context.interior_value(this->_flow_vars.v_var(), qp);
174 libMesh::Gradient grad_T;
175 grad_T = context.interior_gradient(this->_temp_vars.T_var(), qp);
177 libMesh::NumberVectorValue U (u,v);
179 U(2) = context.interior_value(this->_flow_vars.w_var(), qp);
181 const libMesh::Number r = u_qpoint[qp](0);
183 libMesh::Real jac = JxW[qp];
186 libMesh::Real _k_qp = this->_k(context, qp);
188 if( this->_is_axisymmetric )
194 for (
unsigned int i=0; i != n_T_dofs; i++)
197 (-this->_rho*this->_Cp*T_phi[i][qp]*(U*grad_T)
198 -_k_qp*(T_gradphi[i][qp]*grad_T) );
200 if (compute_jacobian)
202 for (
unsigned int j=0; j != n_T_dofs; j++)
207 KTT(i,j) += jac * context.get_elem_solution_derivative() *
208 (-this->_rho*this->_Cp*T_phi[i][qp]*(U*T_gradphi[j][qp])
209 -_k_qp*(T_gradphi[i][qp]*T_gradphi[j][qp]));
213 for (
unsigned int j=0; j != n_u_dofs; j++)
215 KTu(i,j) += jac * context.get_elem_solution_derivative()*(-this->_rho*this->_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(0)));
216 KTv(i,j) += jac * context.get_elem_solution_derivative()*(-this->_rho*this->_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(1)));
218 (*KTw)(i,j) += jac * context.get_elem_solution_derivative()*(-this->_rho*this->_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(2)));
226 #ifdef GRINS_USE_GRVY_TIMERS
227 this->_timer->EndTimer(
"HeatTransfer::element_time_derivative");
238 #ifdef GRINS_USE_GRVY_TIMERS
239 this->_timer->BeginTimer(
"HeatTransfer::side_time_derivative");
242 std::vector<BoundaryID> ids = context.side_boundary_ids();
244 for( std::vector<BoundaryID>::const_iterator it = ids.begin();
245 it != ids.end(); it++ )
247 libmesh_assert (*it != libMesh::BoundaryInfo::invalid_id);
249 this->_bc_handler->apply_neumann_bcs( context, cache, compute_jacobian, *it );
252 #ifdef GRINS_USE_GRVY_TIMERS
253 this->_timer->EndTimer(
"HeatTransfer::side_time_derivative");
264 #ifdef GRINS_USE_GRVY_TIMERS
265 this->_timer->BeginTimer(
"HeatTransfer::mass_residual");
272 const std::vector<libMesh::Real> &JxW =
273 context.get_element_fe(this->_temp_vars.T_var())->get_JxW();
276 const std::vector<std::vector<libMesh::Real> >& phi =
277 context.get_element_fe(this->_temp_vars.T_var())->get_phi();
279 const std::vector<libMesh::Point>& u_qpoint =
280 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
283 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T_var()).size();
286 libMesh::DenseSubVector<libMesh::Real> &F = context.get_elem_residual(this->_temp_vars.T_var());
288 libMesh::DenseSubMatrix<libMesh::Real> &M = context.get_elem_jacobian(this->_temp_vars.T_var(), this->_temp_vars.T_var());
290 unsigned int n_qpoints = context.get_element_qrule().n_points();
292 for (
unsigned int qp = 0; qp != n_qpoints; ++qp)
300 context.interior_rate(this->_temp_vars.T_var(), qp, T_dot);
302 const libMesh::Number r = u_qpoint[qp](0);
304 libMesh::Real jac = JxW[qp];
306 if( this->_is_axisymmetric )
311 for (
unsigned int i = 0; i != n_T_dofs; ++i)
313 F(i) -= this->_rho*this->_Cp*T_dot*phi[i][qp]*jac;
315 if( compute_jacobian )
317 for (
unsigned int j=0; j != n_T_dofs; j++)
320 M(i,j) -= this->_rho*this->_Cp*phi[j][qp]*phi[i][qp]*jac * context.get_elem_solution_rate_derivative();
328 #ifdef GRINS_USE_GRVY_TIMERS
329 this->_timer->EndTimer(
"HeatTransfer::mass_residual");
338 const libMesh::Point& point,
339 libMesh::Real& value )
341 if( quantity_index == this->_k_index )
343 value = this->_k(point, context.get_time());
Base class for reading and handling boundary conditions for physics classes.
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.
GRINS::BCHandlingBase * _bc_handler
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 read_input_options(const GetPot &input)
Read options from GetPot input file.
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.
virtual void side_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for boundaries of elements on the domain boundary.
bool is_axisymmetric() const
const PhysicsName heat_transfer
INSTANTIATE_HEAT_TRANSFER_SUBCLASS(HeatTransfer)
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for HeatTransfer.