36 #include "libmesh/quadrature.h" 
   40   template<
class Mu, 
class SH, 
class TC>
 
   47   template<
class Mu, 
class SH, 
class TC>
 
   53   template<
class Mu, 
class SH, 
class TC>
 
   58 #ifdef GRINS_USE_GRVY_TIMERS 
   59     this->_timer->BeginTimer(
"LowMachNavierStokesSPGSMStabilization::element_time_derivative");
 
   62     this->assemble_continuity_time_deriv( compute_jacobian, context );
 
   63     this->assemble_momentum_time_deriv( compute_jacobian, context );
 
   64     this->assemble_energy_time_deriv( compute_jacobian, context );
 
   66 #ifdef GRINS_USE_GRVY_TIMERS 
   67     this->_timer->EndTimer(
"LowMachNavierStokesSPGSMStabilization::element_time_derivative");
 
   72   template<
class Mu, 
class SH, 
class TC>
 
   77 #ifdef GRINS_USE_GRVY_TIMERS 
   78     this->_timer->BeginTimer(
"LowMachNavierStokesSPGSMStabilization::mass_residual");
 
   81     this->assemble_continuity_mass_residual( compute_jacobian, context );
 
   82     this->assemble_momentum_mass_residual( compute_jacobian, context );
 
   83     this->assemble_energy_mass_residual( compute_jacobian, context );
 
   85 #ifdef GRINS_USE_GRVY_TIMERS 
   86     this->_timer->EndTimer(
"LowMachNavierStokesSPGSMStabilization::mass_residual");
 
   91   template<
class Mu, 
class SH, 
class TC>
 
   96     const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
 
   99     const std::vector<libMesh::Real> &JxW =
 
  100       context.get_element_fe(this->_flow_vars.u())->get_JxW();
 
  103     const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
 
  104       context.get_element_fe(this->_press_var.p())->get_dphi();
 
  106     libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); 
 
  108     unsigned int n_qpoints = context.get_element_qrule().n_points();
 
  110     for (
unsigned int qp=0; qp != n_qpoints; qp++)
 
  112         libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
 
  114         libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
 
  115         libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
 
  117         libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
 
  119         libMesh::Real mu = this->_mu(T);
 
  121         libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
 
  123         libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
 
  124                                  context.interior_value( this->_flow_vars.v(), qp ) );
 
  125         if( this->_dim == 3 )
 
  126           U(2) = context.interior_value( this->_flow_vars.w(), qp );
 
  128         libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
 
  130         libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
 
  134         for (
unsigned int i=0; i != n_p_dofs; i++)
 
  136             Fp(i) += tau_M*RM_s*p_dphi[i][qp]*JxW[qp];
 
  144   template<
class Mu, 
class SH, 
class TC>
 
  149     const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
 
  152     libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
 
  154       libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
 
  157     const std::vector<libMesh::Real> &JxW =
 
  158       context.get_element_fe(this->_flow_vars.u())->get_JxW();
 
  161     const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
 
  162       context.get_element_fe(this->_flow_vars.u())->get_dphi();
 
  164     libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); 
 
  165     libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); 
 
  166     libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
 
  168     if( this->_dim == 3 )
 
  170         Fw  = &context.get_elem_residual(this->_flow_vars.w()); 
 
  173     unsigned int n_qpoints = context.get_element_qrule().n_points();
 
  175     for (
unsigned int qp=0; qp != n_qpoints; qp++)
 
  177         libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
 
  178         libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
 
  180         libMesh::RealGradient U( context.interior_value(this->_flow_vars.u(), qp),
 
  181                                  context.interior_value(this->_flow_vars.v(), qp) );
 
  183         libMesh::RealGradient grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
 
  184         libMesh::RealGradient grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
 
  185         libMesh::RealGradient grad_w;
 
  187         if( this->_dim == 3 )
 
  189             U(2) = context.interior_value(this->_flow_vars.w(), qp);
 
  190             grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
 
  193         libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
 
  195         libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
 
  196         libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
 
  197         libMesh::Real mu = this->_mu(T);
 
  199         libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
 
  200         libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
 
  202         libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
 
  203         libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
 
  205         for (
unsigned int i=0; i != n_u_dofs; i++)
 
  207             Fu(i) += ( - tau_C*RC_s*u_gradphi[i][qp](0)
 
  208                        - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp]  )*JxW[qp];
 
  210             Fv(i) += ( - tau_C*RC_s*u_gradphi[i][qp](1)
 
  211                        - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
 
  213             if( this->_dim == 3 )
 
  215                 (*Fw)(i) += ( - tau_C*RC_s*u_gradphi[i][qp](2)
 
  216                               - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
 
  224   template<
class Mu, 
class SH, 
class TC>
 
  229     const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
 
  232     const std::vector<libMesh::Real> &JxW =
 
  233       context.get_element_fe(this->_temp_vars.T())->get_JxW();
 
  236     const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
 
  237       context.get_element_fe(this->_temp_vars.T())->get_dphi();
 
  239     libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T()); 
 
  241     unsigned int n_qpoints = context.get_element_qrule().n_points();
 
  243     for (
unsigned int qp=0; qp != n_qpoints; qp++)
 
  245         libMesh::Number u, v;
 
  246         u = context.interior_value(this->_flow_vars.u(), qp);
 
  247         v = context.interior_value(this->_flow_vars.v(), qp);
 
  249         libMesh::Gradient grad_T = context.interior_gradient(this->_temp_vars.T(), qp);
 
  251         libMesh::NumberVectorValue U(u,v);
 
  253           U(2) = context.interior_value(this->_flow_vars.w(), qp); 
 
  255         libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
 
  256         libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
 
  258         libMesh::Real k = this->_k(T);
 
  259         libMesh::Real cp = this->_cp(T);
 
  261         libMesh::Number rho_cp = rho*this->_cp(T);
 
  263         libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
 
  265         libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
 
  266         libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
 
  268         libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
 
  270         libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
 
  272         for (
unsigned int i=0; i != n_T_dofs; i++)
 
  274             FT(i) -= rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]*JxW[qp];
 
  282   template<
class Mu, 
class SH, 
class TC>
 
  287     const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
 
  290     const std::vector<libMesh::Real> &JxW =
 
  291       context.get_element_fe(this->_flow_vars.u())->get_JxW();
 
  294     const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
 
  295       context.get_element_fe(this->_press_var.p())->get_dphi();
 
  297     libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); 
 
  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->_flow_vars.u());
 
  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::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
 
  309         libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
 
  311         libMesh::Real mu = this->_mu(T);
 
  313         libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
 
  314                                  context.fixed_interior_value( this->_flow_vars.v(), qp ) );
 
  315         if( this->_dim == 3 )
 
  316           U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
 
  318         libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, 
false );
 
  319         libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
 
  323         for (
unsigned int i=0; i != n_p_dofs; i++)
 
  325             Fp(i) += tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
 
  332   template<
class Mu, 
class SH, 
class TC>
 
  337     const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
 
  340     libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
 
  342       libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
 
  345     const std::vector<libMesh::Real> &JxW =
 
  346       context.get_element_fe(this->_flow_vars.u())->get_JxW();
 
  349     const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
 
  350       context.get_element_fe(this->_flow_vars.u())->get_dphi();
 
  352     libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); 
 
  353     libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); 
 
  354     libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
 
  356     if( this->_dim == 3 )
 
  358         Fw  = &context.get_elem_residual(this->_flow_vars.w()); 
 
  361     unsigned int n_qpoints = context.get_element_qrule().n_points();
 
  362     for (
unsigned int qp=0; qp != n_qpoints; qp++)
 
  364         libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
 
  365         libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
 
  367         libMesh::Real mu = this->_mu(T);
 
  369         libMesh::RealGradient U( context.fixed_interior_value(this->_flow_vars.u(), qp),
 
  370                                  context.fixed_interior_value(this->_flow_vars.v(), qp) );
 
  372         libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_flow_vars.u(), qp);
 
  373         libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_flow_vars.v(), qp);
 
  374         libMesh::RealGradient grad_w;
 
  376         if( this->_dim == 3 )
 
  378             U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
 
  379             grad_w = context.fixed_interior_gradient(this->_flow_vars.w(), qp);
 
  382         libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
 
  384         libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
 
  385         libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
 
  387         libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, 
false );
 
  388         libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
 
  390         libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
 
  391         libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
 
  392         libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
 
  394         for (
unsigned int i=0; i != n_u_dofs; i++)
 
  396             Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
 
  397                        + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp] )*JxW[qp];
 
  399             Fv(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
 
  400                        + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp] )*JxW[qp];
 
  402             if( this->_dim == 3 )
 
  404                 (*Fw)(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
 
  405                               + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp] )*JxW[qp];
 
  413   template<
class Mu, 
class SH, 
class TC>
 
  418     const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
 
  421     const std::vector<libMesh::Real> &JxW =
 
  422       context.get_element_fe(this->_temp_vars.T())->get_JxW();
 
  425     const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
 
  426       context.get_element_fe(this->_temp_vars.T())->get_dphi();
 
  428     libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T()); 
 
  430     unsigned int n_qpoints = context.get_element_qrule().n_points();
 
  432     for (
unsigned int qp=0; qp != n_qpoints; qp++)
 
  434         libMesh::Number u, v;
 
  435         u = context.fixed_interior_value(this->_flow_vars.u(), qp);
 
  436         v = context.fixed_interior_value(this->_flow_vars.v(), qp);
 
  438         libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_temp_vars.T(), qp);
 
  440         libMesh::NumberVectorValue U(u,v);
 
  442           U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp); 
 
  444         libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
 
  445         libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
 
  447         libMesh::Real k = this->_k(T);
 
  448         libMesh::Real cp = this->_cp(T);
 
  450         libMesh::Number rho_cp = rho*this->_cp(T);
 
  452         libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
 
  454         libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
 
  455         libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
 
  457         libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, 
false );
 
  459         libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
 
  461         for (
unsigned int i=0; i != n_T_dofs; i++)
 
  463             FT(i) -= rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]*JxW[qp];
 
virtual ~LowMachNavierStokesSPGSMStabilization()
 
Adds SPGSM-based stabilization to LowMachNavierStokes physics class. 
 
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
 
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors. 
 
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
 
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...
 
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context)
 
LowMachNavierStokesSPGSMStabilization()
 
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &context)
 
Adds VMS-based stabilization to LowMachNavierStokes physics class. 
 
void assemble_continuity_time_deriv(bool compute_jacobian, AssemblyContext &context)
 
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context)