GRINS-0.7.0
List of all members | Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions
GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity > Class Template Reference

Physics class for Incompressible Navier-Stokes. More...

#include <low_mach_navier_stokes.h>

Inheritance diagram for GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >:
Inheritance graph
[legend]
Collaboration diagram for GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >:
Collaboration graph
[legend]

Public Member Functions

 LowMachNavierStokes (const PhysicsName &physics_name, const GetPot &input)
 
 ~LowMachNavierStokes ()
 
virtual void auxiliary_init (MultiphysicsSystem &system)
 Any auxillary initialization a Physics class may need. More...
 
virtual void register_postprocessing_vars (const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
 Register postprocessing variables for LowMachNavierStokes. More...
 
virtual void init_context (AssemblyContext &context)
 Initialize context for added physics variables. More...
 
virtual void element_time_derivative (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Time dependent part(s) of physics for element interiors. More...
 
virtual void element_constraint (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Constraint part(s) of physics for element interiors. More...
 
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. More...
 
virtual void compute_element_time_derivative_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_postprocessed_quantity (unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
 
- Public Member Functions inherited from GRINS::LowMachNavierStokesBase< Viscosity, SpecificHeat, ThermalConductivity >
 LowMachNavierStokesBase (const PhysicsName &physics_name, const std::string &core_physics_name, const GetPot &input)
 
 ~LowMachNavierStokesBase ()
 
virtual void init_variables (libMesh::FEMSystem *system)
 Initialize variables for this physics. More...
 
virtual void set_time_evolving_vars (libMesh::FEMSystem *system)
 Sets velocity variables to be time-evolving. More...
 
libMesh::Real T (const libMesh::Point &p, const AssemblyContext &c) const
 
libMesh::Real rho (libMesh::Real T, libMesh::Real p0) const
 
libMesh::Real d_rho_dT (libMesh::Real T, libMesh::Real p0) const
 
libMesh::Real get_p0_steady (const AssemblyContext &c, unsigned int qp) const
 
libMesh::Real get_p0_steady_side (const AssemblyContext &c, unsigned int qp) const
 
libMesh::Real get_p0_steady (const AssemblyContext &c, const libMesh::Point &p) const
 
libMesh::Real get_p0_transient (AssemblyContext &c, unsigned int qp) const
 
virtual void register_parameter (const std::string &param_name, libMesh::ParameterMultiAccessor< libMesh::Number > &param_pointer) const
 Each subclass will register its copy of an independent. More...
 
- Public Member Functions inherited from GRINS::Physics
 Physics (const GRINS::PhysicsName &physics_name, const GetPot &input)
 
virtual ~Physics ()
 
virtual bool enabled_on_elem (const libMesh::Elem *elem)
 Find if current physics is active on supplied element. More...
 
void set_is_steady (bool is_steady)
 Sets whether this physics is to be solved with a steady solver or not. More...
 
bool is_steady () const
 Returns whether or not this physics is being solved with a steady solver. More...
 
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. More...
 
virtual void nonlocal_time_derivative (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Time dependent part(s) of physics for scalar variables. More...
 
virtual void side_constraint (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Constraint part(s) of physics for boundaries of elements on the domain boundary. More...
 
virtual void nonlocal_constraint (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Constraint part(s) of physics for scalar variables. More...
 
virtual void damping_residual (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Damping matrix part(s) for element interiors. All boundary terms lie within the time_derivative part. More...
 
virtual void nonlocal_mass_residual (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Mass matrix part(s) for scalar variables. More...
 
void init_ics (libMesh::FEMSystem *system, libMesh::CompositeFunction< libMesh::Number > &all_ics)
 
virtual void compute_side_time_derivative_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_nonlocal_time_derivative_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_element_constraint_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_side_constraint_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_nonlocal_constraint_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_damping_residual_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_mass_residual_cache (const AssemblyContext &context, CachedValues &cache)
 
virtual void compute_nonlocal_mass_residual_cache (const AssemblyContext &context, CachedValues &cache)
 
ICHandlingBaseget_ic_handler ()
 
- Public Member Functions inherited from GRINS::ParameterUser
 ParameterUser (const std::string &user_name)
 
virtual ~ParameterUser ()
 
virtual void set_parameter (libMesh::Number &param_variable, const GetPot &input, const std::string &param_name, libMesh::Number param_default)
 Each subclass can simultaneously read a parameter value from. More...
 
virtual void set_parameter (libMesh::ParsedFunction< libMesh::Number, libMesh::Gradient > &func, const GetPot &input, const std::string &func_param_name, const std::string &param_default)
 Each subclass can simultaneously read a parsed function from. More...
 
virtual void set_parameter (libMesh::ParsedFEMFunction< libMesh::Number > &func, const GetPot &input, const std::string &func_param_name, const std::string &param_default)
 Each subclass can simultaneously read a parsed function from. More...
 
virtual void move_parameter (const libMesh::Number &old_parameter, libMesh::Number &new_parameter)
 When cloning an object, we need to update parameter pointers. More...
 
virtual void move_parameter (const libMesh::ParsedFunction< libMesh::Number, libMesh::Gradient > &old_func, libMesh::ParsedFunction< libMesh::Number, libMesh::Gradient > &new_func)
 When cloning an object, we need to update parameter pointers. More...
 
virtual void move_parameter (const libMesh::ParsedFEMFunction< libMesh::Number > &old_func, libMesh::ParsedFEMFunction< libMesh::Number > &new_func)
 When cloning an object, we need to update parameter pointers. More...
 

Protected Member Functions

void assemble_mass_time_deriv (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Helper function. More...
 
void assemble_momentum_time_deriv (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Helper function. More...
 
void assemble_energy_time_deriv (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
 Helper function. More...
 
void assemble_continuity_mass_residual (bool compute_jacobian, AssemblyContext &c)
 Helper function. More...
 
void assemble_momentum_mass_residual (bool compute_jacobian, AssemblyContext &c)
 Helper function. More...
 
void assemble_energy_mass_residual (bool compute_jacobian, AssemblyContext &c)
 Helper function. More...
 
void assemble_thermo_press_elem_time_deriv (bool compute_jacobian, AssemblyContext &c)
 
void assemble_thermo_press_side_time_deriv (bool compute_jacobian, AssemblyContext &c)
 
void assemble_thermo_press_mass_residual (bool compute_jacobian, AssemblyContext &c)
 
- Protected Member Functions inherited from GRINS::Physics
libMesh::UniquePtr< libMesh::FEGenericBase< libMesh::Real > > build_new_fe (const libMesh::Elem *elem, const libMesh::FEGenericBase< libMesh::Real > *fe, const libMesh::Point p)
 
void parse_enabled_subdomains (const GetPot &input, const std::string &physics_name)
 

Protected Attributes

bool _pin_pressure
 Enable pressure pinning. More...
 
PressurePinning _p_pinning
 
unsigned int _rho_index
 Cache index for density post-processing. More...
 
- Protected Attributes inherited from GRINS::LowMachNavierStokesBase< Viscosity, SpecificHeat, ThermalConductivity >
libMesh::Number _p0_over_R
 Thermodynamic pressure divided by gas constant. More...
 
libMesh::Number _p0
 
libMesh::Number _R
 
libMesh::Number _T0
 
unsigned int _dim
 Physical dimension of problem. More...
 
VelocityFEVariables _flow_vars
 
PressureFEVariable _press_var
 
PrimitiveTempFEVariables _temp_vars
 
libMesh::UniquePtr< ThermoPressureFEVariable_p0_var
 
Viscosity _mu
 Viscosity object. More...
 
SpecificHeat _cp
 Specific heat object. More...
 
ThermalConductivity _k
 Thermal conductivity object. More...
 
libMesh::Point _g
 Gravity vector. More...
 
bool _enable_thermo_press_calc
 Flag to enable thermodynamic pressure calculation. More...
 
- Protected Attributes inherited from GRINS::Physics
const PhysicsName _physics_name
 Name of the physics object. Used for reading physics specific inputs. More...
 
GRINS::ICHandlingBase_ic_handler
 
std::set< libMesh::subdomain_id_type > _enabled_subdomains
 Subdomains on which the current Physics class is enabled. More...
 

Private Member Functions

 LowMachNavierStokes ()
 

Additional Inherited Members

- Static Public Member Functions inherited from GRINS::Physics
static void set_is_axisymmetric (bool is_axisymmetric)
 Set whether we should treat the problem as axisymmetric. More...
 
static bool is_axisymmetric ()
 
- Static Public Attributes inherited from GRINS::ParameterUser
static std::string zero_vector_function = std::string("{0}")
 A parseable function string with LIBMESH_DIM components, all 0. More...
 
- Static Protected Attributes inherited from GRINS::Physics
static bool _is_steady = false
 Caches whether or not the solver that's being used is steady or not. More...
 
static bool _is_axisymmetric = false
 Caches whether we are solving an axisymmetric problem or not. More...
 

Detailed Description

template<class Viscosity, class SpecificHeat, class ThermalConductivity>
class GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >

Physics class for Incompressible Navier-Stokes.

This physics class implements the classical Incompressible Navier-Stokes equations.

Definition at line 41 of file low_mach_navier_stokes.h.

Constructor & Destructor Documentation

template<class Mu , class SH , class TC >
GRINS::LowMachNavierStokes< Mu, SH, TC >::LowMachNavierStokes ( const PhysicsName physics_name,
const GetPot &  input 
)

Definition at line 45 of file low_mach_navier_stokes.C.

References GRINS::Physics::_ic_handler, GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >::_pin_pressure, and GRINS::PhysicsNaming::low_mach_navier_stokes().

46  : LowMachNavierStokesBase<Mu,SH,TC>(physics_name,PhysicsNaming::low_mach_navier_stokes(),input),
47  _p_pinning(input,physics_name),
48  _rho_index(0) // Initialize to zero
49  {
50  this->_ic_handler = new GenericICHandler( physics_name, input );
51 
52  this->_pin_pressure = input("Physics/"+PhysicsNaming::low_mach_navier_stokes()+"/pin_pressure", false );
53  }
GRINS::ICHandlingBase * _ic_handler
Definition: physics.h:269
static PhysicsName low_mach_navier_stokes()
unsigned int _rho_index
Cache index for density post-processing.
bool _pin_pressure
Enable pressure pinning.
template<class Viscosity , class SpecificHeat , class ThermalConductivity >
GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >::~LowMachNavierStokes ( )
inline

Definition at line 47 of file low_mach_navier_stokes.h.

47 {};
template<class Viscosity , class SpecificHeat , class ThermalConductivity >
GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >::LowMachNavierStokes ( )
private

Member Function Documentation

template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_continuity_mass_residual ( bool  compute_jacobian,
AssemblyContext c 
)
protected

Helper function.

Definition at line 708 of file low_mach_navier_stokes.C.

710  {
711  // Element Jacobian * quadrature weights for interior integration
712  const std::vector<libMesh::Real> &JxW =
713  context.get_element_fe(this->_flow_vars.u())->get_JxW();
714 
715  // The shape functions at interior quadrature points.
716  const std::vector<std::vector<libMesh::Real> >& p_phi =
717  context.get_element_fe(this->_press_var.p())->get_phi();
718 
719  // The number of local degrees of freedom in each variable
720  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
721 
722  // The subvectors and submatrices we need to fill:
723  libMesh::DenseSubVector<libMesh::Real> &F_p = context.get_elem_residual(this->_press_var.p());
724 
725  unsigned int n_qpoints = context.get_element_qrule().n_points();
726 
727  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
728  {
729  // For the mass residual, we need to be a little careful.
730  // The time integrator is handling the time-discretization
731  // for us so we need to supply M(u_fixed)*u' for the residual.
732  // u_fixed will be given by the fixed_interior_value function
733  // while u' will be given by the interior_rate function.
734  libMesh::Real T_dot;
735  context.interior_rate(this->_temp_vars.T(), qp, T_dot);
736 
737  libMesh::Real T = context.fixed_interior_value(this->_temp_vars.T(), qp);
738 
739  for (unsigned int i = 0; i != n_p_dofs; ++i)
740  {
741  F_p(i) -= T_dot/T*p_phi[i][qp]*JxW[qp];
742  } // End DoF loop i
743 
744  } // End quadrature loop qp
745 
746  return;
747  }
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
VariableIndex p() const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_energy_mass_residual ( bool  compute_jacobian,
AssemblyContext c 
)
protected

Helper function.

Definition at line 830 of file low_mach_navier_stokes.C.

832  {
833  // Element Jacobian * quadrature weights for interior integration
834  const std::vector<libMesh::Real> &JxW =
835  context.get_element_fe(this->_temp_vars.T())->get_JxW();
836 
837  // The shape functions at interior quadrature points.
838  const std::vector<std::vector<libMesh::Real> >& T_phi =
839  context.get_element_fe(this->_temp_vars.T())->get_phi();
840 
841  // The number of local degrees of freedom in each variable
842  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
843 
844  // The subvectors and submatrices we need to fill:
845  libMesh::DenseSubVector<libMesh::Real> &F_T = context.get_elem_residual(this->_temp_vars.T());
846 
847  unsigned int n_qpoints = context.get_element_qrule().n_points();
848 
849  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
850  {
851  // For the mass residual, we need to be a little careful.
852  // The time integrator is handling the time-discretization
853  // for us so we need to supply M(u_fixed)*u' for the residual.
854  // u_fixed will be given by the fixed_interior_value function
855  // while u will be given by the interior_rate function.
856  libMesh::Real T_dot;
857  context.interior_rate(this->_temp_vars.T(), qp, T_dot);
858 
859  libMesh::Real T = context.fixed_interior_value(this->_temp_vars.T(), qp);
860 
861  libMesh::Real cp = this->_cp(T);
862 
863  libMesh::Number rho = this->rho(T, this->get_p0_transient(context, qp));
864 
865  for (unsigned int i = 0; i != n_T_dofs; ++i)
866  {
867  F_T(i) -= rho*cp*T_dot*T_phi[i][qp]*JxW[qp];
868  } // End DoF loop i
869 
870  } // End quadrature loop qp
871 
872  return;
873  }
SpecificHeat _cp
Specific heat object.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
libMesh::Real get_p0_transient(AssemblyContext &c, unsigned int qp) const
libMesh::Real rho(libMesh::Real T, libMesh::Real p0) const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_energy_time_deriv ( bool  compute_jacobian,
AssemblyContext context,
CachedValues cache 
)
protected

Helper function.

Definition at line 596 of file low_mach_navier_stokes.C.

References GRINS::CachedValues::get_cached_gradient_values(), GRINS::CachedValues::get_cached_values(), GRINS::Cache::TEMPERATURE, GRINS::Cache::TEMPERATURE_GRAD, GRINS::Cache::THERMO_PRESSURE, GRINS::Cache::X_VELOCITY, GRINS::Cache::Y_VELOCITY, and GRINS::Cache::Z_VELOCITY.

599  {
600  // The number of local degrees of freedom in each variable.
601  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
602  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
603 
604  // Element Jacobian * quadrature weights for interior integration.
605  const std::vector<libMesh::Real> &JxW =
606  context.get_element_fe(this->_temp_vars.T())->get_JxW();
607 
608  // The temperature shape functions at interior quadrature points.
609  const std::vector<std::vector<libMesh::Real> >& T_phi =
610  context.get_element_fe(this->_temp_vars.T())->get_phi();
611  const std::vector<std::vector<libMesh::Real> >& u_phi =
612  context.get_element_fe(this->_flow_vars.u())->get_phi();
613 
614  // The temperature shape functions gradients at interior quadrature points.
615  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
616  context.get_element_fe(this->_temp_vars.T())->get_dphi();
617 
618  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T()); // R_{T}
619 
620  unsigned int n_qpoints = context.get_element_qrule().n_points();
621  for (unsigned int qp=0; qp != n_qpoints; qp++)
622  {
623  libMesh::Number u, v, T, p0;
624  u = cache.get_cached_values(Cache::X_VELOCITY)[qp];
625  v = cache.get_cached_values(Cache::Y_VELOCITY)[qp];
626  T = cache.get_cached_values(Cache::TEMPERATURE)[qp];
627  p0 = cache.get_cached_values(Cache::THERMO_PRESSURE)[qp];
628 
629  libMesh::Gradient grad_T = cache.get_cached_gradient_values(Cache::TEMPERATURE_GRAD)[qp];
630 
631  libMesh::NumberVectorValue U(u,v);
632  if (this->_dim == 3)
633  U(2) = cache.get_cached_values(Cache::Z_VELOCITY)[qp]; // w
634 
635  libMesh::DenseSubMatrix<libMesh::Number> &KTu = context.get_elem_jacobian(this->_temp_vars.T(), this->_flow_vars.u()); // R_{u},{u}
636  libMesh::DenseSubMatrix<libMesh::Number> &KTv = context.get_elem_jacobian(this->_temp_vars.T(), this->_flow_vars.v()); // R_{u},{u}
637  libMesh::DenseSubMatrix<libMesh::Number>* KTw = NULL;
638 
639  libMesh::DenseSubMatrix<libMesh::Number> &KTT = context.get_elem_jacobian(this->_temp_vars.T(), this->_temp_vars.T()); // R_{u},{u}
640 
641  if( this->_dim == 3 )
642  {
643  KTw = &context.get_elem_jacobian(this->_temp_vars.T(), this->_flow_vars.w()); // R_{u},{w}
644  }
645 
646  libMesh::Number k = this->_k(T);
647  libMesh::Number dk_dT = this->_k.deriv(T);
648  libMesh::Number cp = this->_cp(T);
649  libMesh::Number d_cp = this->_cp.deriv(T);
650  libMesh::Number rho = this->rho( T, p0 );
651  libMesh::Number d_rho = this->d_rho_dT( T, p0 );
652 
653  // Now a loop over the pressure degrees of freedom. This
654  // computes the contributions of the continuity equation.
655  for (unsigned int i=0; i != n_T_dofs; i++)
656  {
657  FT(i) += ( -rho*cp*U*grad_T*T_phi[i][qp] // convection term
658  - k*grad_T*T_gradphi[i][qp] // diffusion term
659  )*JxW[qp];
660 
661 
662  if(compute_jacobian)
663  {
664  for (unsigned int j=0; j!=n_u_dofs; j++)
665  {
666  //pre-compute repeated term
667  libMesh::Number r0 = rho*cp*T_phi[i][qp]*u_phi[j][qp];
668 
669  KTu(i,j) += JxW[qp]*
670  -r0*grad_T(0);
671  //-rho*cp*u_phi[j][qp]*grad_T(0)*T_phi[i][qp];
672 
673  KTv(i,j) += JxW[qp]*
674  -r0*grad_T(1);
675  //-rho*cp*u_phi[j][qp]*grad_T(1)*T_phi[i][qp];
676 
677  if (this->_dim == 3)
678  {
679  (*KTw)(i,j) += JxW[qp]*
680  -r0*grad_T(2);
681  //-rho*cp*u_phi[j][qp]*grad_T(2)*T_phi[i][qp];
682  }
683 
684  } // end u_dofs loop (j)
685 
686 
687  for (unsigned int j=0; j!=n_T_dofs; j++)
688  {
689  KTT(i,j) += JxW[qp]* (
690  -rho*(
691  cp*U*T_phi[i][qp]*T_gradphi[j][qp]
692  + U*grad_T*T_phi[i][qp]*d_cp*T_phi[j][qp]
693  )
694  -cp*U*grad_T*T_phi[i][qp]*d_rho*T_phi[j][qp]
695  -k*T_gradphi[i][qp]*T_gradphi[j][qp]
696  -grad_T*T_gradphi[i][qp]*dk_dT*T_phi[j][qp]
697  );
698  } // end T_dofs loop (j)
699 
700  } // end if compute_jacobian
701  } // end outer T_dofs loop (i)
702  } //end qp loop
703 
704  return;
705  }
unsigned int _dim
Physical dimension of problem.
SpecificHeat _cp
Specific heat object.
ThermalConductivity _k
Thermal conductivity object.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
libMesh::Real d_rho_dT(libMesh::Real T, libMesh::Real p0) const
libMesh::Real rho(libMesh::Real T, libMesh::Real p0) const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_mass_time_deriv ( bool  compute_jacobian,
AssemblyContext context,
CachedValues cache 
)
protected

Helper function.

Definition at line 169 of file low_mach_navier_stokes.C.

References GRINS::CachedValues::get_cached_gradient_values(), GRINS::CachedValues::get_cached_values(), GRINS::Cache::TEMPERATURE, GRINS::Cache::TEMPERATURE_GRAD, GRINS::Cache::X_VELOCITY, GRINS::Cache::X_VELOCITY_GRAD, GRINS::Cache::Y_VELOCITY, GRINS::Cache::Y_VELOCITY_GRAD, GRINS::Cache::Z_VELOCITY, and GRINS::Cache::Z_VELOCITY_GRAD.

172  {
173  // The number of local degrees of freedom in each variable.
174  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
175 
176  // Element Jacobian * quadrature weights for interior integration.
177  const std::vector<libMesh::Real> &JxW =
178  context.get_element_fe(this->_flow_vars.u())->get_JxW();
179 
180  // The pressure shape functions at interior quadrature points.
181  const std::vector<std::vector<libMesh::Real> >& p_phi =
182  context.get_element_fe(this->_press_var.p())->get_phi();
183 
184 
185 
186  // The velocity shape functions at interior quadrature points.
187  const std::vector<std::vector<libMesh::Real> >& u_phi =
188  context.get_element_fe(this->_flow_vars.u())->get_phi();
189 
190  // The velocity shape function gradients (in global coords.)
191  // at interior quadrature points.
192  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
193  context.get_element_fe(this->_flow_vars.u())->get_dphi();
194 
195  // The temperature shape functions at interior quadrature points.
196  const std::vector<std::vector<libMesh::Real> >& T_phi =
197  context.get_element_fe(this->_temp_vars.T())->get_phi();
198 
199  // The temperature shape function gradients (in global coords.)
200  // at interior quadrature points.
201  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
202  context.get_element_fe(this->_temp_vars.T())->get_dphi();
203 
204 
205 
206  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
207 
208  unsigned int n_qpoints = context.get_element_qrule().n_points();
209 
210  // The number of local degrees of freedom in each variable.
211  const unsigned int n_t_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
212 
213  // The number of local degrees of freedom in each variable.
214  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
215 
216  // Check number of dofs is same for _flow_vars.u(), v_var and w_var.
217  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
218  if (this->_dim == 3)
219  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
220 
221  for (unsigned int qp=0; qp != n_qpoints; qp++)
222  {
223  libMesh::Number u, v, T;
224  u = cache.get_cached_values(Cache::X_VELOCITY)[qp];
225  v = cache.get_cached_values(Cache::Y_VELOCITY)[qp];
226 
227  T = cache.get_cached_values(Cache::TEMPERATURE)[qp];
228 
229  libMesh::Gradient grad_u = cache.get_cached_gradient_values(Cache::X_VELOCITY_GRAD)[qp];
230  libMesh::Gradient grad_v = cache.get_cached_gradient_values(Cache::Y_VELOCITY_GRAD)[qp];
231 
232  libMesh::Gradient grad_T = cache.get_cached_gradient_values(Cache::TEMPERATURE_GRAD)[qp];
233 
234  libMesh::NumberVectorValue U(u,v);
235  if (this->_dim == 3)
236  U(2) = cache.get_cached_values(Cache::Z_VELOCITY)[qp]; // w
237 
238  libMesh::Number divU = grad_u(0) + grad_v(1);
239  if (this->_dim == 3)
240  {
241  libMesh::Gradient grad_w = cache.get_cached_gradient_values(Cache::Z_VELOCITY_GRAD)[qp];
242  divU += grad_w(2);
243  }
244 
245 
246  libMesh::DenseSubMatrix<libMesh::Number> &KPu = context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u());
247  libMesh::DenseSubMatrix<libMesh::Number> &KPv = context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v());
248  libMesh::DenseSubMatrix<libMesh::Number> &KPT = context.get_elem_jacobian(this->_press_var.p(), this->_temp_vars.T());
249 
250  libMesh::DenseSubMatrix<libMesh::Number>* KPw = NULL;
251 
252  if( this->_dim == 3 )
253  {
254  KPw = &context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.w());
255  }
256 
257  // Now a loop over the pressure degrees of freedom. This
258  // computes the contributions of the continuity equation.
259  for (unsigned int i=0; i != n_p_dofs; i++)
260  {
261  Fp(i) += (-U*grad_T/T + divU)*p_phi[i][qp]*JxW[qp];
262 
263 
264  if (compute_jacobian)
265  {
266 
267  for (unsigned int j=0; j!=n_u_dofs; j++)
268  {
269  KPu(i,j) += JxW[qp]*(
270  +u_gradphi[j][qp](0)*p_phi[i][qp]
271  -u_phi[j][qp]*p_phi[i][qp]*grad_T(0)/T
272  );
273 
274  KPv(i,j) += JxW[qp]*(
275  +u_gradphi[j][qp](1)*p_phi[i][qp]
276  -u_phi[j][qp]*p_phi[i][qp]*grad_T(1)/T
277  );
278 
279  if (this->_dim == 3)
280  {
281  (*KPw)(i,j) += JxW[qp]*(
282  +u_gradphi[j][qp](2)*p_phi[i][qp]
283  -u_phi[j][qp]*p_phi[i][qp]*grad_T(2)/T
284  );
285  }
286 
287  }
288 
289  for (unsigned int j=0; j!=n_t_dofs; j++)
290  {
291  KPT(i,j) += JxW[qp]*(
292  -T_gradphi[j][qp]*U*p_phi[i][qp]/T
293  +U*p_phi[i][qp]*grad_T*T_phi[j][qp])/(T*T);
294  }
295  } // end if compute_jacobian
296  } // end p_dofs loop
297  } // end qp loop
298 
299  return;
300  }
unsigned int _dim
Physical dimension of problem.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
VariableIndex p() const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_momentum_mass_residual ( bool  compute_jacobian,
AssemblyContext c 
)
protected

Helper function.

Definition at line 750 of file low_mach_navier_stokes.C.

752  {
753  // Element Jacobian * quadrature weights for interior integration
754  const std::vector<libMesh::Real> &JxW =
755  context.get_element_fe(this->_flow_vars.u())->get_JxW();
756 
757  // The shape functions at interior quadrature points.
758  const std::vector<std::vector<libMesh::Real> >& u_phi =
759  context.get_element_fe(this->_flow_vars.u())->get_phi();
760 
761  // The number of local degrees of freedom in each variable
762  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
763 
764  // The subvectors and submatrices we need to fill:
765  libMesh::DenseSubVector<libMesh::Real> &F_u = context.get_elem_residual(this->_flow_vars.u());
766  libMesh::DenseSubVector<libMesh::Real> &F_v = context.get_elem_residual(this->_flow_vars.v());
767  libMesh::DenseSubVector<libMesh::Real>* F_w = NULL;
768 
769  if( this->_dim == 3 )
770  {
771  F_w = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
772  }
773 
774  unsigned int n_qpoints = context.get_element_qrule().n_points();
775 
776  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
777  {
778  // For the mass residual, we need to be a little careful.
779  // The time integrator is handling the time-discretization
780  // for us so we need to supply M(u_fixed)*u' for the residual.
781  // u_fixed will be given by the fixed_interior_value function
782  // while u' will be given by the interior_rate function.
783  libMesh::Real u_dot, v_dot, w_dot = 0.0;
784  context.interior_rate(this->_flow_vars.u(), qp, u_dot);
785  context.interior_rate(this->_flow_vars.v(), qp, v_dot);
786 
787  if( this->_dim == 3 )
788  context.interior_rate(this->_flow_vars.w(), qp, w_dot);
789 
790  libMesh::Real T = context.fixed_interior_value(this->_temp_vars.T(), qp);
791 
792  libMesh::Number rho = this->rho(T, this->get_p0_transient(context, qp));
793 
794  for (unsigned int i = 0; i != n_u_dofs; ++i)
795  {
796  F_u(i) -= rho*u_dot*u_phi[i][qp]*JxW[qp];
797  F_v(i) -= rho*v_dot*u_phi[i][qp]*JxW[qp];
798 
799  if( this->_dim == 3 )
800  (*F_w)(i) -= rho*w_dot*u_phi[i][qp]*JxW[qp];
801 
802  /*
803  if( compute_jacobian )
804  {
805  for (unsigned int j=0; j != n_u_dofs; j++)
806  {
807  // Assuming rho is constant w.r.t. u, v, w
808  // and T (if Boussinesq added).
809  libMesh::Real value = JxW[qp]*_rho*u_phi[i][qp]*u_phi[j][qp];
810 
811  M_uu(i,j) += value;
812  M_vv(i,j) += value;
813 
814  if( _dim == 3)
815  {
816  M_ww(i,j) += value;
817  }
818 
819  } // End DoF loop j
820  } // End Jacobian check
821  */
822 
823  } // End DoF loop i
824  } // End quadrature loop qp
825 
826  return;
827  }
unsigned int _dim
Physical dimension of problem.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
libMesh::Real get_p0_transient(AssemblyContext &c, unsigned int qp) const
libMesh::Real rho(libMesh::Real T, libMesh::Real p0) const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_momentum_time_deriv ( bool  compute_jacobian,
AssemblyContext context,
CachedValues cache 
)
protected

Helper function.

Definition at line 303 of file low_mach_navier_stokes.C.

References GRINS::CachedValues::get_cached_gradient_values(), GRINS::CachedValues::get_cached_values(), GRINS::Cache::PRESSURE, GRINS::Cache::TEMPERATURE, GRINS::Cache::THERMO_PRESSURE, GRINS::Cache::X_VELOCITY, GRINS::Cache::X_VELOCITY_GRAD, GRINS::Cache::Y_VELOCITY, GRINS::Cache::Y_VELOCITY_GRAD, GRINS::Cache::Z_VELOCITY, and GRINS::Cache::Z_VELOCITY_GRAD.

306  {
307  // The number of local degrees of freedom in each variable.
308  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
309  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
310  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
311 
312  // Check number of dofs is same for _flow_vars.u(), v_var and w_var.
313  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
314  if (this->_dim == 3)
315  libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
316 
317  // Element Jacobian * quadrature weights for interior integration.
318  const std::vector<libMesh::Real> &JxW =
319  context.get_element_fe(this->_flow_vars.u())->get_JxW();
320 
321  // The pressure shape functions at interior quadrature points.
322  const std::vector<std::vector<libMesh::Real> >& u_phi =
323  context.get_element_fe(this->_flow_vars.u())->get_phi();
324  const std::vector<std::vector<libMesh::Real> >& p_phi =
325  context.get_element_fe(this->_press_var.p())->get_phi();
326  const std::vector<std::vector<libMesh::Real> >& T_phi =
327  context.get_element_fe(this->_temp_vars.T())->get_phi();
328 
329  // The velocity shape function gradients at interior quadrature points.
330  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
331  context.get_element_fe(this->_flow_vars.u())->get_dphi();
332 
333  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{u}
334  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{v}
335  libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
336 
337  if( this->_dim == 3 )
338  {
339  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
340  }
341 
342  unsigned int n_qpoints = context.get_element_qrule().n_points();
343  for (unsigned int qp=0; qp != n_qpoints; qp++)
344  {
345  libMesh::Number u, v, p, p0, T;
346  u = cache.get_cached_values(Cache::X_VELOCITY)[qp];
347  v = cache.get_cached_values(Cache::Y_VELOCITY)[qp];
348 
349  T = cache.get_cached_values(Cache::TEMPERATURE)[qp];
350  p = cache.get_cached_values(Cache::PRESSURE)[qp];
351  p0 = cache.get_cached_values(Cache::THERMO_PRESSURE)[qp];
352 
353  libMesh::Gradient grad_u = cache.get_cached_gradient_values(Cache::X_VELOCITY_GRAD)[qp];
354  libMesh::Gradient grad_v = cache.get_cached_gradient_values(Cache::Y_VELOCITY_GRAD)[qp];
355 
356  libMesh::Gradient grad_w;
357  if (this->_dim == 3)
358  grad_w = cache.get_cached_gradient_values(Cache::Z_VELOCITY_GRAD)[qp];
359 
360  libMesh::NumberVectorValue grad_uT( grad_u(0), grad_v(0) );
361  libMesh::NumberVectorValue grad_vT( grad_u(1), grad_v(1) );
362  libMesh::NumberVectorValue grad_wT;
363  if( this->_dim == 3 )
364  {
365  grad_uT(2) = grad_w(0);
366  grad_vT(2) = grad_w(1);
367  grad_wT = libMesh::NumberVectorValue( grad_u(2), grad_v(2), grad_w(2) );
368  }
369 
370  libMesh::NumberVectorValue U(u,v);
371  if (this->_dim == 3)
372  U(2) = cache.get_cached_values(Cache::Z_VELOCITY)[qp]; // w
373 
374  libMesh::Number divU = grad_u(0) + grad_v(1);
375  if (this->_dim == 3)
376  divU += grad_w(2);
377 
378  libMesh::Number rho = this->rho( T, p0 );
379  libMesh::Number d_rho = this->d_rho_dT( T, p0 );
380  libMesh::Number d_mu = this->_mu.deriv(T);
381 
382  libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u()); // R_{u},{u}
383  libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v()); // R_{u},{v}
384  libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
385 
386  libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u()); // R_{v},{u}
387  libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v()); // R_{v},{v}
388  libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
389 
390  libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
391  libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
392  libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
393 
394  libMesh::DenseSubMatrix<libMesh::Number> &Kup = context.get_elem_jacobian(this->_flow_vars.u(), this->_press_var.p()); // R_{u},{p}
395  libMesh::DenseSubMatrix<libMesh::Number> &Kvp = context.get_elem_jacobian(this->_flow_vars.v(), this->_press_var.p()); // R_{v},{p}
396  libMesh::DenseSubMatrix<libMesh::Number>* Kwp = NULL;
397 
398  libMesh::DenseSubMatrix<libMesh::Number> &KuT = context.get_elem_jacobian(this->_flow_vars.u(), this->_temp_vars.T()); // R_{u},{p}
399  libMesh::DenseSubMatrix<libMesh::Number> &KvT = context.get_elem_jacobian(this->_flow_vars.v(), this->_temp_vars.T()); // R_{v},{p}
400  libMesh::DenseSubMatrix<libMesh::Number>* KwT = NULL;
401 
402  if( this->_dim == 3 )
403  {
404  Kuw = &context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.w()); // R_{u},{w}
405  Kvw = &context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.w()); // R_{v},{w}
406  Kwu = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.u()); // R_{w},{u};
407  Kwv = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.v()); // R_{w},{v};
408  Kww = &context.get_elem_jacobian(this->_flow_vars.w(), this->_flow_vars.w()); // R_{w},{w}
409  Kwp = &context.get_elem_jacobian(this->_flow_vars.w(), this->_press_var.p()); // R_{w},{p}
410  KwT = &context.get_elem_jacobian(this->_flow_vars.w(), this->_temp_vars.T()); // R_{w},{T}
411  }
412 
413  // Now a loop over the pressure degrees of freedom. This
414  // computes the contributions of the continuity equation.
415  for (unsigned int i=0; i != n_u_dofs; i++)
416  {
417  Fu(i) += ( -rho*U*grad_u*u_phi[i][qp] // convection term
418  + p*u_gradphi[i][qp](0) // pressure term
419  - this->_mu(T)*(u_gradphi[i][qp]*grad_u + u_gradphi[i][qp]*grad_uT
420  - 2.0/3.0*divU*u_gradphi[i][qp](0) ) // diffusion term
421  + rho*this->_g(0)*u_phi[i][qp] // hydrostatic term
422  )*JxW[qp];
423 
424  Fv(i) += ( -rho*U*grad_v*u_phi[i][qp] // convection term
425  + p*u_gradphi[i][qp](1) // pressure term
426  - this->_mu(T)*(u_gradphi[i][qp]*grad_v + u_gradphi[i][qp]*grad_vT
427  - 2.0/3.0*divU*u_gradphi[i][qp](1) ) // diffusion term
428  + rho*this->_g(1)*u_phi[i][qp] // hydrostatic term
429  )*JxW[qp];
430  if (this->_dim == 3)
431  {
432  (*Fw)(i) += ( -rho*U*grad_w*u_phi[i][qp] // convection term
433  + p*u_gradphi[i][qp](2) // pressure term
434  - this->_mu(T)*(u_gradphi[i][qp]*grad_w + u_gradphi[i][qp]*grad_wT
435  - 2.0/3.0*divU*u_gradphi[i][qp](2) ) // diffusion term
436  + rho*this->_g(2)*u_phi[i][qp] // hydrostatic term
437  )*JxW[qp];
438  }
439 
440  if (compute_jacobian && context.get_elem_solution_derivative())
441  {
442  libmesh_assert (context.get_elem_solution_derivative() == 1.0);
443 
444  for (unsigned int j=0; j != n_u_dofs; j++)
445  {
446  //precompute repeated terms
447  libMesh::Number r0 = rho*U*u_phi[i][qp]*u_gradphi[j][qp];
448  libMesh::Number r1 = u_gradphi[i][qp]*u_gradphi[j][qp];
449  libMesh::Number r2 = rho*u_phi[i][qp]*u_phi[j][qp];
450 
451 
452  Kuu(i,j) += JxW[qp]*(
453  -r0
454  //-rho*U*u_gradphi[j][qp]*u_phi[i][qp]
455  -r2*grad_u(0)
456  //-rho*u_phi[i][qp]*grad_u(0)*u_phi[j][qp]
457  -this->_mu(T)*(
458  r1
459  //u_gradphi[i][qp]*u_gradphi[j][qp]
460  + u_gradphi[i][qp](0)*u_gradphi[j][qp](0) // transpose
461  - 2.0/3.0*u_gradphi[i][qp](0)*u_gradphi[j][qp](0)
462  ));
463  Kvv(i,j) += JxW[qp]*(
464  -r0
465  //-rho*U*u_gradphi[j][qp]*u_phi[i][qp]
466  -r2*grad_v(1)
467  //-rho*u_phi[i][qp]*grad_v(1)*u_phi[j][qp]
468  -this->_mu(T)*(
469  r1
470  //u_gradphi[i][qp]*u_gradphi[j][qp]
471  + u_gradphi[i][qp](1)*u_gradphi[j][qp](1) // transpose
472  - 2.0/3.0*u_gradphi[i][qp](1)*u_gradphi[j][qp](1)
473  ));
474 
475  Kuv(i,j) += JxW[qp]*(
476  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](1)
477  -this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](0)
478  -r2*grad_u(1)
479  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_u(1));
480  );
481 
482  Kvu(i,j) += JxW[qp]*(
483  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](0)
484  -this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](1)
485  -r2*grad_v(0)
486  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_v(0));
487  );
488 
489 
490 
491  if (this->_dim == 3)
492  {
493  (*Kuw)(i,j) += JxW[qp]*(
494  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](2)
495  -this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](0)
496  -r2*grad_u(2)
497  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_u(2)
498  );
499 
500  (*Kvw)(i,j) += JxW[qp]*(
501  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](2)
502  -this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](1)
503  -r2*grad_v(2)
504  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_v(2)
505  );
506 
507  (*Kwu)(i,j) += JxW[qp]*(
508  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](0)
509  -this->_mu(T)*u_gradphi[i][qp](0)*u_gradphi[j][qp](2)
510  -r2*grad_w(0)
511  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_w(0)
512  );
513 
514  (*Kwv)(i,j) += JxW[qp]*(
515  +2.0/3.0*this->_mu(T)*u_gradphi[i][qp](2)*u_gradphi[j][qp](1)
516  -this->_mu(T)*u_gradphi[i][qp](1)*u_gradphi[j][qp](2)
517  -r2*grad_w(1)
518  //-rho*u_phi[i][qp]*u_phi[j][qp]*grad_w(1)
519  );
520 
521  (*Kww)(i,j) += JxW[qp]*(
522  -r0
523  //-rho*U*u_gradphi[j][qp]*u_phi[i][qp]
524  -r2*grad_w(2)
525  //-rho*u_phi[i][qp]*grad_w(2)*u_phi[j][qp]
526  -this->_mu(T)*(
527  r1
528  //u_gradphi[i][qp]*u_gradphi[j][qp]
529  + u_gradphi[i][qp](2)*u_gradphi[j][qp](2) // transpose
530  - 2.0/3.0*u_gradphi[i][qp](2)*u_gradphi[j][qp](2)
531  ));
532 
533  } // end if _dim==3
534  } // end of the inner dof (j) loop
535 
536  for (unsigned int j=0; j!=n_T_dofs; j++)
537  {
538 
539  //precompute repeated term
540  libMesh:: Number r3 = d_rho*u_phi[i][qp]*T_phi[j][qp];
541 
542  // Analytical Jacobains
543  KuT(i,j) += JxW[qp]*(
544  -r3*U*grad_u
545  //-d_rho*T_phi[j][qp]*U*grad_u*u_phi[i][qp]
546  -d_mu*T_phi[j][qp]*grad_u*u_gradphi[i][qp]
547  -d_mu*T_phi[j][qp]*grad_u(0)*u_gradphi[i][qp](0) // transpose
548  +2.0/3.0*d_mu*T_phi[j][qp]*divU*u_gradphi[i][qp](0)
549  +r3*this->_g(0)
550  //+d_rho*T_phi[j][qp]*this->_g(0)*u_phi[i][qp]
551  );
552 
553  KvT(i,j) += JxW[qp]*(
554  -r3*U*grad_v
555  //-d_rho*T_phi[j][qp]*U*grad_v*u_phi[i][qp]
556  -d_mu*T_phi[j][qp]*grad_v*u_gradphi[i][qp]
557  -d_mu*T_phi[j][qp]*grad_v(1)*u_gradphi[i][qp](1) // transpose
558  +2.0/3.0*d_mu*T_phi[j][qp]*divU*u_gradphi[i][qp](1)
559  +r3*this->_g(1)
560  //+d_rho*T_phi[j][qp]*this->_g(1)*u_phi[i][qp]
561  );
562 
563  if (this->_dim == 3)
564  {
565  (*KwT)(i,j) += JxW[qp]*(
566  -r3*U*grad_w
567  //-d_rho*T_phi[j][qp]*U*grad_v*u_phi[i][qp]
568  -d_mu*T_phi[j][qp]*grad_w*u_gradphi[i][qp]
569  -d_mu*T_phi[j][qp]*grad_w(2)*u_gradphi[i][qp](2) // transpose
570  +2.0/3.0*d_mu*T_phi[j][qp]*divU*u_gradphi[i][qp](2)
571  +r3*this->_g(2)
572  //+d_rho*T_phi[j][qp]*this->_g(2)*u_phi[i][qp]
573  );
574 
575  } // end if _dim==3
576  } // end T_dofs loop
577 
578  // Matrix contributions for the up, vp and wp couplings
579  for (unsigned int j=0; j != n_p_dofs; j++)
580  {
581  Kup(i,j) += JxW[qp]*p_phi[j][qp]*u_gradphi[i][qp](0);
582  Kvp(i,j) += JxW[qp]*p_phi[j][qp]*u_gradphi[i][qp](1);
583  if (this->_dim == 3)
584  (*Kwp)(i,j) += JxW[qp]*p_phi[j][qp]*u_gradphi[i][qp](2);
585  } // end of the inner dof (j) loop
586 
587  } // end - if (compute_jacobian && context.get_elem_solution_derivative())
588 
589  } // end of the outer dof (i) loop
590  } // end of the quadrature point (qp) loop
591 
592  return;
593  }
unsigned int _dim
Physical dimension of problem.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
VariableIndex p() const
libMesh::Point _g
Gravity vector.
libMesh::Real d_rho_dT(libMesh::Real T, libMesh::Real p0) const
libMesh::Real rho(libMesh::Real T, libMesh::Real p0) const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_thermo_press_elem_time_deriv ( bool  compute_jacobian,
AssemblyContext c 
)
protected

Definition at line 876 of file low_mach_navier_stokes.C.

878  {
879  // Element Jacobian * quadrature weights for interior integration
880  const std::vector<libMesh::Real> &JxW =
881  context.get_element_fe(this->_temp_vars.T())->get_JxW();
882 
883  // The number of local degrees of freedom in each variable
884  const unsigned int n_p0_dofs = context.get_dof_indices(this->_p0_var->p0()).size();
885 
886  // The subvectors and submatrices we need to fill:
887  libMesh::DenseSubVector<libMesh::Real> &F_p0 = context.get_elem_residual(this->_p0_var->p0());
888 
889  unsigned int n_qpoints = context.get_element_qrule().n_points();
890 
891  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
892  {
893  libMesh::Number T;
894  T = context.interior_value(this->_temp_vars.T(), qp);
895 
896  libMesh::Gradient grad_u, grad_v, grad_w;
897  grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
898  grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
899  if (this->_dim == 3)
900  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
901 
902  libMesh::Number divU = grad_u(0) + grad_v(1);
903  if(this->_dim==3)
904  divU += grad_w(2);
905 
906  //libMesh::Number cp = this->_cp(T);
907  //libMesh::Number cv = cp + this->_R;
908  //libMesh::Number gamma = cp/cv;
909  //libMesh::Number gamma_ratio = gamma/(gamma-1.0);
910 
911  libMesh::Number p0 = context.interior_value( this->_p0_var->p0(), qp );
912 
913  for (unsigned int i = 0; i != n_p0_dofs; ++i)
914  {
915  F_p0(i) += (p0/T - this->_p0/this->_T0)*JxW[qp];
916  //F_p0(i) -= p0*gamma_ratio*divU*JxW[qp];
917  } // End DoF loop i
918  }
919 
920  return;
921  }
libMesh::UniquePtr< ThermoPressureFEVariable > _p0_var
unsigned int _dim
Physical dimension of problem.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_thermo_press_mass_residual ( bool  compute_jacobian,
AssemblyContext c 
)
protected

Definition at line 974 of file low_mach_navier_stokes.C.

976  {
977  // The number of local degrees of freedom in each variable.
978  const unsigned int n_p0_dofs = context.get_dof_indices(this->_p0_var->p0()).size();
979  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
980  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
981 
982  // Element Jacobian * quadrature weights for interior integration
983  const std::vector<libMesh::Real> &JxW =
984  context.get_element_fe(this->_temp_vars.T())->get_JxW();
985 
986  // The temperature shape functions at interior quadrature points.
987  const std::vector<std::vector<libMesh::Real> >& T_phi =
988  context.get_element_fe(this->_temp_vars.T())->get_phi();
989 
990  // The temperature shape functions at interior quadrature points.
991  const std::vector<std::vector<libMesh::Real> >& p_phi =
992  context.get_element_fe(this->_press_var.p())->get_phi();
993 
994  // The subvectors and submatrices we need to fill:
995  libMesh::DenseSubVector<libMesh::Real> &F_p0 = context.get_elem_residual(this->_p0_var->p0());
996  libMesh::DenseSubVector<libMesh::Real> &F_T = context.get_elem_residual(this->_temp_vars.T());
997  libMesh::DenseSubVector<libMesh::Real> &F_p = context.get_elem_residual(this->_press_var.p());
998 
999  unsigned int n_qpoints = context.get_element_qrule().n_points();
1000 
1001  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
1002  {
1003  libMesh::Number T;
1004  T = context.fixed_interior_value(this->_temp_vars.T(), qp);
1005 
1006  libMesh::Number cp = this->_cp(T);
1007  libMesh::Number cv = cp + this->_R;
1008  libMesh::Number gamma = cp/cv;
1009  libMesh::Number one_over_gamma = 1.0/(gamma-1.0);
1010 
1011  libMesh::Number p0_dot;
1012  context.interior_rate(this->_p0_var->p0(), qp, p0_dot);
1013 
1014  libMesh::Number p0 = context.fixed_interior_value(this->_p0_var->p0(), qp );
1015 
1016  for (unsigned int i=0; i != n_p0_dofs; i++)
1017  {
1018  F_p0(i) -= p0_dot*one_over_gamma*JxW[qp];
1019  }
1020 
1021  for (unsigned int i=0; i != n_T_dofs; i++)
1022  {
1023  F_T(i) += p0_dot*T_phi[i][qp]*JxW[qp];
1024  }
1025 
1026  for (unsigned int i=0; i != n_p_dofs; i++)
1027  {
1028  F_p(i) += p0_dot/p0*p_phi[i][qp]*JxW[qp];
1029  }
1030 
1031  }
1032  return;
1033  }
libMesh::UniquePtr< ThermoPressureFEVariable > _p0_var
SpecificHeat _cp
Specific heat object.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
VariableIndex p() const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::assemble_thermo_press_side_time_deriv ( bool  compute_jacobian,
AssemblyContext c 
)
protected

Definition at line 924 of file low_mach_navier_stokes.C.

926  {
927  // The number of local degrees of freedom in each variable.
928  const unsigned int n_p0_dofs = context.get_dof_indices(this->_p0_var->p0()).size();
929 
930  // Element Jacobian * quadrature weight for side integration.
931  //const std::vector<libMesh::Real> &JxW_side = context.get_side_fe(this->_temp_vars.T())->get_JxW();
932 
933  //const std::vector<Point> &normals = context.get_side_fe(this->_temp_vars.T())->get_normals();
934 
935  //libMesh::DenseSubVector<libMesh::Number> &F_p0 = context.get_elem_residual(this->_p0_var->p0()); // residual
936 
937  // Physical location of the quadrature points
938  //const std::vector<libMesh::Point>& qpoint = context.get_side_fe(this->_temp_vars.T())->get_xyz();
939 
940  unsigned int n_qpoints = context.get_side_qrule().n_points();
941  for (unsigned int qp=0; qp != n_qpoints; qp++)
942  {
943  /*
944  libMesh::Number T = context.side_value( this->_temp_vars.T(), qp );
945  libMesh::Gradient U = ( context.side_value( this->_flow_vars.u(), qp ),
946  context.side_value( this->_flow_vars.v(), qp ) );
947  libMesh::Gradient grad_T = context.side_gradient( this->_temp_vars.T(), qp );
948 
949 
950  libMesh::Number p0 = context.side_value( this->_p0_var->p0(), qp );
951 
952  libMesh::Number k = this->_k(T);
953  libMesh::Number cp = this->_cp(T);
954 
955  libMesh::Number cv = cp + this->_R;
956  libMesh::Number gamma = cp/cv;
957  libMesh::Number gamma_ratio = gamma/(gamma-1.0);
958  */
959 
960  //std::cout << "U = " << U << std::endl;
961 
962  //std::cout << "x = " << qpoint[qp] << ", grad_T = " << grad_T << std::endl;
963 
964  for (unsigned int i=0; i != n_p0_dofs; i++)
965  {
966  //F_p0(i) += (k*grad_T*normals[qp] - p0*gamma_ratio*U*normals[qp] )*JxW_side[qp];
967  }
968  }
969 
970  return;
971  }
libMesh::UniquePtr< ThermoPressureFEVariable > _p0_var
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::auxiliary_init ( MultiphysicsSystem system)
virtual

Any auxillary initialization a Physics class may need.

This is called after all variables are added, so this method can safely query the MultiphysicsSystem about variable information.

Reimplemented from GRINS::Physics.

Definition at line 56 of file low_mach_navier_stokes.C.

57  {
58  if( _pin_pressure )
59  _p_pinning.check_pin_location(system.get_mesh());
60  }
void check_pin_location(const libMesh::MeshBase &mesh)
Check the mesh to ensure pin location is found.
bool _pin_pressure
Enable pressure pinning.
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::compute_element_time_derivative_cache ( const AssemblyContext context,
CachedValues cache 
)
virtual

Reimplemented from GRINS::Physics.

Definition at line 1036 of file low_mach_navier_stokes.C.

References GRINS::Cache::PRESSURE, GRINS::CachedValues::set_gradient_values(), GRINS::CachedValues::set_values(), GRINS::Cache::TEMPERATURE, GRINS::Cache::TEMPERATURE_GRAD, GRINS::Cache::THERMO_PRESSURE, GRINS::Cache::X_VELOCITY, GRINS::Cache::X_VELOCITY_GRAD, GRINS::Cache::Y_VELOCITY, GRINS::Cache::Y_VELOCITY_GRAD, GRINS::Cache::Z_VELOCITY, and GRINS::Cache::Z_VELOCITY_GRAD.

1038  {
1039  const unsigned int n_qpoints = context.get_element_qrule().n_points();
1040 
1041  std::vector<libMesh::Real> u, v, w, T, p, p0;
1042  u.resize(n_qpoints);
1043  v.resize(n_qpoints);
1044  if( this->_dim > 2 )
1045  w.resize(n_qpoints);
1046 
1047  T.resize(n_qpoints);
1048  p.resize(n_qpoints);
1049  p0.resize(n_qpoints);
1050 
1051  std::vector<libMesh::Gradient> grad_u, grad_v, grad_w, grad_T;
1052  grad_u.resize(n_qpoints);
1053  grad_v.resize(n_qpoints);
1054  if( this->_dim > 2 )
1055  grad_w.resize(n_qpoints);
1056 
1057  grad_T.resize(n_qpoints);
1058 
1059  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
1060  {
1061  u[qp] = context.interior_value(this->_flow_vars.u(), qp);
1062  v[qp] = context.interior_value(this->_flow_vars.v(), qp);
1063 
1064  grad_u[qp] = context.interior_gradient(this->_flow_vars.u(), qp);
1065  grad_v[qp] = context.interior_gradient(this->_flow_vars.v(), qp);
1066  if( this->_dim > 2 )
1067  {
1068  w[qp] = context.interior_value(this->_flow_vars.w(), qp);
1069  grad_w[qp] = context.interior_gradient(this->_flow_vars.w(), qp);
1070  }
1071  T[qp] = context.interior_value(this->_temp_vars.T(), qp);
1072  grad_T[qp] = context.interior_gradient(this->_temp_vars.T(), qp);
1073 
1074  p[qp] = context.interior_value(this->_press_var.p(), qp);
1075  p0[qp] = this->get_p0_steady(context, qp);
1076  }
1077 
1078  cache.set_values(Cache::X_VELOCITY, u);
1079  cache.set_values(Cache::Y_VELOCITY, v);
1080 
1081  cache.set_gradient_values(Cache::X_VELOCITY_GRAD, grad_u);
1082  cache.set_gradient_values(Cache::Y_VELOCITY_GRAD, grad_v);
1083 
1084  if(this->_dim > 2)
1085  {
1086  cache.set_values(Cache::Z_VELOCITY, w);
1087  cache.set_gradient_values(Cache::Z_VELOCITY_GRAD, grad_w);
1088  }
1089 
1090  cache.set_values(Cache::TEMPERATURE, T);
1091  cache.set_gradient_values(Cache::TEMPERATURE_GRAD, grad_T);
1092 
1093  cache.set_values(Cache::PRESSURE, p);
1094  cache.set_values(Cache::THERMO_PRESSURE, p0);
1095 
1096  return;
1097  }
libMesh::Real get_p0_steady(const AssemblyContext &c, unsigned int qp) const
unsigned int _dim
Physical dimension of problem.
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
VariableIndex p() const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::compute_postprocessed_quantity ( unsigned int  quantity_index,
const AssemblyContext context,
const libMesh::Point &  point,
libMesh::Real &  value 
)
virtual

Reimplemented from GRINS::Physics.

Definition at line 1100 of file low_mach_navier_stokes.C.

1104  {
1105  if( quantity_index == this->_rho_index )
1106  {
1107  libMesh::Real T = this->T(point,context);
1108  libMesh::Real p0 = this->get_p0_steady(context,point);
1109 
1110  value = this->rho( T, p0 );
1111  }
1112 
1113  return;
1114  }
libMesh::Real get_p0_steady(const AssemblyContext &c, unsigned int qp) const
libMesh::Real T(const libMesh::Point &p, const AssemblyContext &c) const
unsigned int _rho_index
Cache index for density post-processing.
libMesh::Real rho(libMesh::Real T, libMesh::Real p0) const
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::element_constraint ( bool  compute_jacobian,
AssemblyContext context,
CachedValues cache 
)
virtual

Constraint part(s) of physics for element interiors.

Reimplemented from GRINS::Physics.

Definition at line 138 of file low_mach_navier_stokes.C.

141  {
142  // Pin p = p_value at p_point
143  if( this->_pin_pressure )
144  {
145  this->_p_pinning.pin_value( context, compute_jacobian, this->_press_var.p());
146  }
147 
148  return;
149  }
VariableIndex p() const
void pin_value(libMesh::DiffContext &context, const bool request_jacobian, const GRINS::VariableIndex var, const double penalty=1.0)
The idea here is to pin a variable to a particular value if there is a null space - e...
bool _pin_pressure
Enable pressure pinning.
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::element_time_derivative ( bool  compute_jacobian,
AssemblyContext context,
CachedValues cache 
)
virtual

Time dependent part(s) of physics for element interiors.

Reimplemented from GRINS::Physics.

Definition at line 115 of file low_mach_navier_stokes.C.

118  {
119 #ifdef GRINS_USE_GRVY_TIMERS
120  this->_timer->BeginTimer("LowMachNavierStokes::element_time_derivative");
121 #endif
122 
123  this->assemble_mass_time_deriv( compute_jacobian, context, cache );
124  this->assemble_momentum_time_deriv( compute_jacobian, context, cache );
125  this->assemble_energy_time_deriv( compute_jacobian, context, cache );
126 
127  if( this->_enable_thermo_press_calc )
128  this->assemble_thermo_press_elem_time_deriv( compute_jacobian, context );
129 
130 #ifdef GRINS_USE_GRVY_TIMERS
131  this->_timer->EndTimer("LowMachNavierStokes::element_time_derivative");
132 #endif
133 
134  return;
135  }
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Helper function.
bool _enable_thermo_press_calc
Flag to enable thermodynamic pressure calculation.
void assemble_thermo_press_elem_time_deriv(bool compute_jacobian, AssemblyContext &c)
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Helper function.
void assemble_mass_time_deriv(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Helper function.
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::init_context ( AssemblyContext context)
virtual

Initialize context for added physics variables.

Reimplemented from GRINS::LowMachNavierStokesBase< Viscosity, SpecificHeat, ThermalConductivity >.

Definition at line 94 of file low_mach_navier_stokes.C.

References GRINS::LowMachNavierStokesBase< Viscosity, SpecificHeat, ThermalConductivity >::init_context().

95  {
96  // First call base class
98 
99  // We also need the side shape functions, etc.
100  context.get_side_fe(this->_flow_vars.u())->get_JxW();
101  context.get_side_fe(this->_flow_vars.u())->get_phi();
102  context.get_side_fe(this->_flow_vars.u())->get_dphi();
103  context.get_side_fe(this->_flow_vars.u())->get_xyz();
104 
105  context.get_side_fe(this->_temp_vars.T())->get_JxW();
106  context.get_side_fe(this->_temp_vars.T())->get_phi();
107  context.get_side_fe(this->_temp_vars.T())->get_dphi();
108  context.get_side_fe(this->_temp_vars.T())->get_xyz();
109 
110  return;
111  }
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::mass_residual ( bool  compute_jacobian,
AssemblyContext context,
CachedValues cache 
)
virtual

Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part.

Reimplemented from GRINS::Physics.

Definition at line 152 of file low_mach_navier_stokes.C.

155  {
156  this->assemble_continuity_mass_residual( compute_jacobian, context );
157 
158  this->assemble_momentum_mass_residual( compute_jacobian, context );
159 
160  this->assemble_energy_mass_residual( compute_jacobian, context );
161 
162  if( this->_enable_thermo_press_calc )
163  this->assemble_thermo_press_mass_residual( compute_jacobian, context );
164 
165  return;
166  }
bool _enable_thermo_press_calc
Flag to enable thermodynamic pressure calculation.
void assemble_thermo_press_mass_residual(bool compute_jacobian, AssemblyContext &c)
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &c)
Helper function.
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &c)
Helper function.
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &c)
Helper function.
template<class Mu , class SH , class TC >
void GRINS::LowMachNavierStokes< Mu, SH, TC >::register_postprocessing_vars ( const GetPot &  input,
PostProcessedQuantities< libMesh::Real > &  postprocessing 
)
virtual

Register postprocessing variables for LowMachNavierStokes.

Reimplemented from GRINS::Physics.

Definition at line 63 of file low_mach_navier_stokes.C.

References GRINS::PhysicsNaming::low_mach_navier_stokes(), and GRINS::PostProcessedQuantities< NumericType >::register_quantity().

65  {
66  std::string section = "Physics/"+PhysicsNaming::low_mach_navier_stokes()+"/output_vars";
67 
68  if( input.have_variable(section) )
69  {
70  unsigned int n_vars = input.vector_variable_size(section);
71 
72  for( unsigned int v = 0; v < n_vars; v++ )
73  {
74  std::string name = input(section,"DIE!",v);
75 
76  if( name == std::string("rho") )
77  {
78  this->_rho_index = postprocessing.register_quantity( name );
79  }
80  else
81  {
82  std::cerr << "Error: Invalue output_vars value for "+PhysicsNaming::low_mach_navier_stokes() << std::endl
83  << " Found " << name << std::endl
84  << " Acceptable values are: rho" << std::endl;
85  libmesh_error();
86  }
87  }
88  }
89 
90  return;
91  }
static PhysicsName low_mach_navier_stokes()
unsigned int _rho_index
Cache index for density post-processing.

Member Data Documentation

template<class Viscosity , class SpecificHeat , class ThermalConductivity >
PressurePinning GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >::_p_pinning
protected

Definition at line 85 of file low_mach_navier_stokes.h.

template<class Viscosity , class SpecificHeat , class ThermalConductivity >
bool GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >::_pin_pressure
protected
template<class Viscosity , class SpecificHeat , class ThermalConductivity >
unsigned int GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >::_rho_index
protected

Cache index for density post-processing.

Definition at line 88 of file low_mach_navier_stokes.h.


The documentation for this class was generated from the following files:

Generated on Thu Jun 2 2016 21:52:31 for GRINS-0.7.0 by  doxygen 1.8.10