26 #include "grins_config.h"
34 #include "libmesh/getpot.h"
35 #include "libmesh/fem_system.h"
36 #include "libmesh/quadrature.h"
43 _flow_stab_helper(physics_name+
"FlowStabHelper", input),
44 _temp_stab_helper(physics_name+
"TempStabHelper", input),
72 #ifdef GRINS_USE_GRVY_TIMERS
73 this->_timer->BeginTimer(
"BoussinesqBuoyancySPGSMStabilization::element_time_derivative");
77 const unsigned int n_u_dofs = context.get_dof_indices(_flow_vars.u_var()).size();
80 const std::vector<libMesh::Real> &JxW =
81 context.get_element_fe(_flow_vars.u_var())->get_JxW();
88 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
89 context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
92 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(_flow_vars.u_var());
93 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(_flow_vars.v_var());
94 libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
97 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
106 unsigned int n_qpoints = context.get_element_qrule().n_points();
108 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
110 for (
unsigned int qp=0; qp != n_qpoints; qp++)
112 libMesh::RealGradient g = this->_flow_stab_helper.compute_g( fe, context, qp );
113 libMesh::RealTensor G = this->_flow_stab_helper.compute_G( fe, context, qp );
115 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
116 context.interior_value( this->_flow_vars.v_var(), qp ) );
117 if( this->_dim == 3 )
119 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
123 libMesh::Real mu_qp = this->_mu(context, qp);
125 libMesh::Real tau_M = this->_flow_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, mu_qp, this->_is_steady );
133 T = context.interior_value(_temp_vars.T_var(), qp);
135 libMesh::RealGradient residual = _rho_ref*_beta_T*(T-_T_ref)*_g;
137 for (
unsigned int i=0; i != n_u_dofs; i++)
139 Fu(i) += ( -tau_M*residual(0)*_rho*U*u_gradphi[i][qp] )*JxW[qp];
142 Fv(i) += ( -tau_M*residual(1)*_rho*U*u_gradphi[i][qp] )*JxW[qp];
147 (*Fw)(i) += ( -tau_M*residual(2)*_rho*U*u_gradphi[i][qp] )*JxW[qp];
151 if (compute_jacobian)
153 libmesh_not_implemented();
159 #ifdef GRINS_USE_GRVY_TIMERS
160 this->_timer->EndTimer(
"BoussinesqBuoyancySPGSMStabilization::element_time_derivative");
171 #ifdef GRINS_USE_GRVY_TIMERS
172 this->_timer->BeginTimer(
"BoussinesqBuoyancySPGSMStabilization::element_constraint");
176 const unsigned int n_p_dofs = context.get_dof_indices(_flow_vars.p_var()).size();
179 const std::vector<libMesh::Real> &JxW =
180 context.get_element_fe(_flow_vars.u_var())->get_JxW();
182 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
183 context.get_element_fe(this->_flow_vars.p_var())->get_dphi();
185 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_flow_vars.p_var());
193 unsigned int n_qpoints = context.get_element_qrule().n_points();
195 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u_var());
197 for (
unsigned int qp=0; qp != n_qpoints; qp++)
199 libMesh::RealGradient g = this->_flow_stab_helper.compute_g( fe, context, qp );
200 libMesh::RealTensor G = this->_flow_stab_helper.compute_G( fe, context, qp );
202 libMesh::RealGradient U( context.interior_value( this->_flow_vars.u_var(), qp ),
203 context.interior_value( this->_flow_vars.v_var(), qp ) );
204 if( this->_dim == 3 )
206 U(2) = context.interior_value( this->_flow_vars.w_var(), qp );
210 libMesh::Real mu_qp = this->_mu(context, qp);
212 libMesh::Real tau_M = this->_flow_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, mu_qp, this->_is_steady );
216 T = context.interior_value(_temp_vars.T_var(), qp);
218 libMesh::RealGradient residual = _rho_ref*_beta_T*(T-_T_ref)*_g;
223 for (
unsigned int i=0; i != n_p_dofs; i++)
225 Fp(i) += -tau_M*residual*p_dphi[i][qp]*JxW[qp];
227 if (compute_jacobian)
229 libmesh_not_implemented();
236 #ifdef GRINS_USE_GRVY_TIMERS
237 this->_timer->EndTimer(
"BoussinesqBuoyancySPGSMStabilization::element_constraint");
328 (
const std::string & param_name,
333 _mu.register_parameter(param_name, param_pointer);
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
~BoussinesqBuoyancySPGSMStabilization()
INSTANTIATE_INC_NS_SUBCLASS(BoussinesqBuoyancySPGSMStabilization)
virtual void register_parameter(const std::string ¶m_name, libMesh::ParameterMultiPointer< libMesh::Number > ¶m_pointer) const
Each subclass will register its copy of an independent.
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 void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
virtual void register_parameter(const std::string ¶m_name, libMesh::ParameterMultiPointer< libMesh::Number > ¶m_pointer) const
Each subclass will register its copy of an independent.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
const PhysicsName heat_transfer
BoussinesqBuoyancySPGSMStabilization()