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

#include <inc_navier_stokes.h>

Inheritance diagram for GRINS::IncompressibleNavierStokes< Viscosity >:

Collaboration diagram for GRINS::IncompressibleNavierStokes< Viscosity >:

Public Member Functions
	IncompressibleNavierStokes (const std::string &physics_name, const GetPot &input)

	~IncompressibleNavierStokes ()

virtual void	read_input_options (const GetPot &input)
	Read options from GetPot input file. More...

virtual void	register_postprocessing_vars (const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
	Register postprocessing variables for IncompressibleNavierStokes. 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	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	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_postprocessed_quantity (unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
	Compute value of postprocessed quantities at libMesh::Point. More...

virtual void	init_variables (libMesh::FEMSystem *system)
	Initialization of Navier-Stokes variables. More...

virtual void	set_time_evolving_vars (libMesh::FEMSystem *system)
	Sets velocity variables to be time-evolving. More...

virtual void	init_context (AssemblyContext &context)
	Initialize context for added physics variables. More...

virtual void	register_parameter (const std::string &param_name, libMesh::ParameterMultiPointer< libMesh::Number > &param_pointer) const
	Each subclass will register its copy of an independent. More...

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	auxiliary_init (MultiphysicsSystem &system)
	Any auxillary initialization a Physics class may need. 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	nonlocal_constraint (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
	Constraint part(s) of physics for scalar variables. More...

virtual void	nonlocal_mass_residual (bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
	Mass matrix part(s) for scalar variables. More...

void	init_bcs (libMesh::FEMSystem *system)

void	init_ics (libMesh::FEMSystem *system, libMesh::CompositeFunction< libMesh::Number > &all_ics)

void	attach_neumann_bound_func (GRINS::NBCContainer &neumann_bcs)

void	attach_dirichlet_bound_func (const GRINS::DBCContainer &dirichlet_bc)

virtual void	compute_element_time_derivative_cache (const AssemblyContext &context, CachedValues &cache)

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_mass_residual_cache (const AssemblyContext &context, CachedValues &cache)

virtual void	compute_nonlocal_mass_residual_cache (const AssemblyContext &context, CachedValues &cache)

BCHandlingBase *	get_bc_handler ()

ICHandlingBase *	get_ic_handler ()

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...

Protected Attributes
PressurePinning	_p_pinning

bool	_pin_pressure
	Enable pressure pinning. More...

unsigned int	_mu_index
	Index from registering this quantity. More...

unsigned int	_dim
	Physical dimension of problem. More...

PrimitiveFlowFEVariables	_flow_vars

libMesh::Number	_rho
	Material parameters, read from input. More...

Viscosity	_mu
	Viscosity object. More...

const PhysicsName	_physics_name
	Name of the physics object. Used for reading physics specific inputs. More...

GRINS::BCHandlingBase *	_bc_handler

GRINS::ICHandlingBase *	_ic_handler

std::set< libMesh::subdomain_id_type >	_enabled_subdomains
	Subdomains on which the current Physics class is enabled. More...

bool	_is_axisymmetric

Static Protected Attributes
static bool	_is_steady = false
	Caches whether or not the solver that's being used is steady or not. More...

Private Member Functions
	IncompressibleNavierStokes ()

Detailed Description

template<class Viscosity>
class GRINS::IncompressibleNavierStokes< Viscosity >

Physics class for Incompressible Navier-Stokes.

This physics class implements the classical Incompressible Navier-Stokes equations. This is a templated class, the class Viscosity can be instantiated as a specific type (right now:ConstantViscosity or SpatiallyVaryingViscosity) to allow the user to specify a constant or spatially varying viscosity in the input file

Definition at line 44 of file inc_navier_stokes.h.

Constructor & Destructor Documentation

template<class Mu >

GRINS::IncompressibleNavierStokes< Mu >::IncompressibleNavierStokes	(	const std::string &	physics_name,
		const GetPot &	input
	)

Definition at line 43 of file inc_navier_stokes.C.

References GRINS::Physics::_bc_handler, GRINS::Physics::_ic_handler, GRINS::Physics::_is_axisymmetric, GRINS::BCHandlingBase::is_axisymmetric(), and GRINS::IncompressibleNavierStokes< Viscosity >::read_input_options().

     : IncompressibleNavierStokesBase<Mu>(physics_name,
                                          incompressible_navier_stokes, /* "core" Physics name */
                                          input),
     _p_pinning(input,physics_name),
     _mu_index(0)
   {
     this->read_input_options(input);
 
     // This is deleted in the base class
     this->_bc_handler = new IncompressibleNavierStokesBCHandling( physics_name, input );
 
     if( this->_bc_handler->is_axisymmetric() )
       {
         this->_is_axisymmetric = true;
       }
 
     this->_ic_handler = new GenericICHandler( physics_name, input );
 
     return;
   }

template<class Mu >

GRINS::IncompressibleNavierStokes< Mu >::~IncompressibleNavierStokes ( )

Definition at line 66 of file inc_navier_stokes.C.

   {
     return;
   }

template<class Viscosity >

GRINS::IncompressibleNavierStokes< Viscosity >::IncompressibleNavierStokes ( )

private

Member Function Documentation

void GRINS::Physics::attach_dirichlet_bound_func ( const GRINS::DBCContainer & dirichlet_bc )

inherited

Definition at line 150 of file physics.C.

References GRINS::Physics::_bc_handler, and GRINS::BCHandlingBase::attach_dirichlet_bound_func().

   {
     _bc_handler->attach_dirichlet_bound_func( dirichlet_bc );
     return;
   }

void GRINS::Physics::attach_neumann_bound_func ( GRINS::NBCContainer & neumann_bcs )

inherited

Definition at line 144 of file physics.C.

References GRINS::Physics::_bc_handler, and GRINS::BCHandlingBase::attach_neumann_bound_func().

   {
     _bc_handler->attach_neumann_bound_func( neumann_bcs );
     return;
   }

void GRINS::Physics::auxiliary_init ( MultiphysicsSystem & system )

virtualinherited

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.

Definition at line 113 of file physics.C.

   {
     return;
   }

void GRINS::Physics::compute_element_constraint_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Definition at line 185 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::element_constraint().

   {
     return;
   }

void GRINS::Physics::compute_element_time_derivative_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Reimplemented in GRINS::LowMachNavierStokes< Viscosity, SpecificHeat, ThermalConductivity >, and GRINS::ReactingLowMachNavierStokes< Mixture, Evaluator >.

Definition at line 167 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::element_time_derivative().

   {
     return;
   }

void GRINS::Physics::compute_mass_residual_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Definition at line 203 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::mass_residual().

   {
     return;
   }

void GRINS::Physics::compute_nonlocal_constraint_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Definition at line 197 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::nonlocal_constraint().

   {
     return;
   }

void GRINS::Physics::compute_nonlocal_mass_residual_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Definition at line 209 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::nonlocal_mass_residual().

   {
     return;
   }

void GRINS::Physics::compute_nonlocal_time_derivative_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Definition at line 179 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::nonlocal_time_derivative().

   {
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokes< Mu >::compute_postprocessed_quantity	(	unsigned int	quantity_index,
		const AssemblyContext &	context,
		const libMesh::Point &	point,
		libMesh::Real &	value
	)

virtual

Compute value of postprocessed quantities at libMesh::Point.

Reimplemented from GRINS::Physics.

Definition at line 575 of file inc_navier_stokes.C.

   {
     if( quantity_index == this->_mu_index )
       {
         value = this->_mu(point, context.get_time());
       }
 
     return;
   }

void GRINS::Physics::compute_side_constraint_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Definition at line 191 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::side_constraint().

   {
     return;
   }

void GRINS::Physics::compute_side_time_derivative_cache	(	const AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Reimplemented in GRINS::ReactingLowMachNavierStokes< Mixture, Evaluator >.

Definition at line 173 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::side_time_derivative().

   {
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokes< Mu >::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 351 of file inc_navier_stokes.C.

   {
 #ifdef GRINS_USE_GRVY_TIMERS
     this->_timer->BeginTimer("IncompressibleNavierStokes::element_constraint");
 #endif
 
     // The number of local degrees of freedom in each variable.
     const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
     const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
 
     // We get some references to cell-specific data that
     // will be used to assemble the linear system.
 
     // Element Jacobian * quadrature weights for interior integration.
     const std::vector<libMesh::Real> &JxW =
       context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
 
     // The velocity shape function gradients (in global coords.)
     // at interior quadrature points.
     const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
       context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
 
     // The velocity shape function gradients (in global coords.)
     // at interior quadrature points.
     const std::vector<std::vector<libMesh::Real> >& u_phi =
       context.get_element_fe(this->_flow_vars.u_var())->get_phi();
 
     // The pressure shape functions at interior quadrature points.
     const std::vector<std::vector<libMesh::Real> >& p_phi =
       context.get_element_fe(this->_flow_vars.p_var())->get_phi();
 
     const std::vector<libMesh::Point>& u_qpoint = 
       context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
     
     // The subvectors and submatrices we need to fill:
     //
     // Kpu, Kpv, Kpw, Fp
 
     libMesh::DenseSubMatrix<libMesh::Number> &Kpu = context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.u_var()); // R_{p},{u}
     libMesh::DenseSubMatrix<libMesh::Number> &Kpv = context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.v_var()); // R_{p},{v}
     libMesh::DenseSubMatrix<libMesh::Number>* Kpw = NULL;
 
     libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_flow_vars.p_var()); // R_{p}
 
     if( this->_dim == 3 )
       {
         Kpw = &context.get_elem_jacobian(this->_flow_vars.p_var(), this->_flow_vars.w_var()); // R_{p},{w}
       }
 
     // Add the constraint given by the continuity equation.
     unsigned int n_qpoints = context.get_element_qrule().n_points();
     for (unsigned int qp=0; qp != n_qpoints; qp++)
       {
         // Compute the velocity gradient at the old Newton iterate.
         libMesh::Gradient grad_u, grad_v, grad_w;
         grad_u = context.interior_gradient(this->_flow_vars.u_var(), qp);
         grad_v = context.interior_gradient(this->_flow_vars.v_var(), qp);
         if (this->_dim == 3)
           grad_w = context.interior_gradient(this->_flow_vars.w_var(), qp);
 
         libMesh::Number divU = grad_u(0) + grad_v(1);
         if (this->_dim == 3)
           divU += grad_w(2);
 
         const libMesh::Number r = u_qpoint[qp](0);
 
         libMesh::Real jac = JxW[qp];
 
         if( this->_is_axisymmetric )
           {
             libMesh::Number u = context.interior_value( this->_flow_vars.u_var(), qp );
             divU += u/r;
             jac *= r;
           }
 
         // Now a loop over the pressure degrees of freedom.  This
         // computes the contributions of the continuity equation.
         for (unsigned int i=0; i != n_p_dofs; i++)
           {
             Fp(i) += p_phi[i][qp]*divU*jac;
 
             if (compute_jacobian)
               {
                 libmesh_assert_equal_to (context.get_elem_solution_derivative(), 1.0);
 
                 for (unsigned int j=0; j != n_u_dofs; j++)
                   {
                     Kpu(i,j) += p_phi[i][qp]*u_gradphi[j][qp](0)*jac * context.get_elem_solution_derivative();
                     Kpv(i,j) += p_phi[i][qp]*u_gradphi[j][qp](1)*jac * context.get_elem_solution_derivative();
                     if (this->_dim == 3)
                       (*Kpw)(i,j) += p_phi[i][qp]*u_gradphi[j][qp](2)*jac * context.get_elem_solution_derivative();
 
                     if( this->_is_axisymmetric )
                       {
                         Kpu(i,j) += p_phi[i][qp]*u_phi[j][qp]/r*jac * context.get_elem_solution_derivative();
                       }
                   } // end of the inner dof (j) loop
 
               } // end - if (compute_jacobian)
 
           } // end of the outer dof (i) loop
       } // end of the quadrature point (qp) loop
 
 
     
   
 
 #ifdef GRINS_USE_GRVY_TIMERS
     this->_timer->EndTimer("IncompressibleNavierStokes::element_constraint");
 #endif
 
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokes< Mu >::element_time_derivative	(	bool	compute_jacobian,
		AssemblyContext &	context,
		CachedValues &	cache
	)

virtual

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

Todo:: Would it be better to put this in its own DoF loop and do the if check once?

Reimplemented from GRINS::Physics.

Definition at line 114 of file inc_navier_stokes.C.

   {
 #ifdef GRINS_USE_GRVY_TIMERS
     this->_timer->BeginTimer("IncompressibleNavierStokes::element_time_derivative");
 #endif
 
     // The number of local degrees of freedom in each variable.
     const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
     const unsigned int n_p_dofs = context.get_dof_indices(this->_flow_vars.p_var()).size();
 
     // Check number of dofs is same for this->_flow_vars.u_var(), v_var and w_var.
     libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v_var()).size());
     if (this->_dim == 3)
       libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w_var()).size());
 
     // We get some references to cell-specific data that
     // will be used to assemble the linear system.
 
     // Element Jacobian * quadrature weights for interior integration.
     const std::vector<libMesh::Real> &JxW =
       context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
 
     // The velocity shape functions at interior quadrature points.
     const std::vector<std::vector<libMesh::Real> >& u_phi =
       context.get_element_fe(this->_flow_vars.u_var())->get_phi();
 
     // The velocity shape function gradients (in global coords.)
     // at interior quadrature points.
     const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
       context.get_element_fe(this->_flow_vars.u_var())->get_dphi();
 
     // The pressure shape functions at interior quadrature points.
     const std::vector<std::vector<libMesh::Real> >& p_phi =
       context.get_element_fe(this->_flow_vars.p_var())->get_phi();
 
     const std::vector<libMesh::Point>& u_qpoint = 
       context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
 
     // The subvectors and submatrices we need to fill:
     //
     // K_{\alpha \beta} = R_{\alpha},{\beta} = \partial{ R_{\alpha} } / \partial{ {\beta} } (where R denotes residual)
     // e.g., for \alpha = v and \beta = u we get: K{vu} = R_{v},{u}
     // Note that Kpu, Kpv, Kpw and Fp comes as constraint.
 
     libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var()); // R_{u},{u}
     libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.v_var()); // R_{u},{v}
     libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
 
     libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.u_var()); // R_{v},{u}
     libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var()); // R_{v},{v}
     libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
 
     libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
     libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
     libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
 
     libMesh::DenseSubMatrix<libMesh::Number> &Kup = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.p_var()); // R_{u},{p}
     libMesh::DenseSubMatrix<libMesh::Number> &Kvp = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.p_var()); // R_{v},{p}
     libMesh::DenseSubMatrix<libMesh::Number>* Kwp = NULL;
 
     libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var()); // R_{u}
     libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var()); // R_{v}
     libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
 
     if( this->_dim == 3 )
       {
         Kuw = &context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.w_var()); // R_{u},{w}
         Kvw = &context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.w_var()); // R_{v},{w}
         Kwu = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.u_var()); // R_{w},{u};
         Kwv = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.v_var()); // R_{w},{v};
         Kww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var()); // R_{w},{w}
         Kwp = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.p_var()); // R_{w},{p}
         Fw  = &context.get_elem_residual(this->_flow_vars.w_var()); // R_{w}
       }
 
     // Now we will build the element Jacobian and residual.
     // Constructing the residual requires the solution and its
     // gradient from the previous timestep.  This must be
     // calculated at each quadrature point by summing the
     // solution degree-of-freedom values by the appropriate
     // weight functions.
     unsigned int n_qpoints = context.get_element_qrule().n_points();
 
     for (unsigned int qp=0; qp != n_qpoints; qp++)
       {
         // Compute the solution & its gradient at the old Newton iterate.
         libMesh::Number p, u, v;
         p = context.interior_value(this->_flow_vars.p_var(), qp);
         u = context.interior_value(this->_flow_vars.u_var(), qp);
         v = context.interior_value(this->_flow_vars.v_var(), qp);
 
         libMesh::Gradient grad_u, grad_v, grad_w;
         grad_u = context.interior_gradient(this->_flow_vars.u_var(), qp);
         grad_v = context.interior_gradient(this->_flow_vars.v_var(), qp);
         if (this->_dim == 3)
           grad_w = context.interior_gradient(this->_flow_vars.w_var(), qp);
 
         libMesh::NumberVectorValue U(u,v);
         if (this->_dim == 3)
           U(2) = context.interior_value(this->_flow_vars.w_var(), qp); // w
 
         const libMesh::Number  grad_u_x = grad_u(0);
         const libMesh::Number  grad_u_y = grad_u(1);
         const libMesh::Number  grad_u_z = (this->_dim == 3)?grad_u(2):0;
         const libMesh::Number  grad_v_x = grad_v(0);
         const libMesh::Number  grad_v_y = grad_v(1);
         const libMesh::Number  grad_v_z = (this->_dim == 3)?grad_v(2):0;
         const libMesh::Number  grad_w_x = (this->_dim == 3)?grad_w(0):0;
         const libMesh::Number  grad_w_y = (this->_dim == 3)?grad_w(1):0;
         const libMesh::Number  grad_w_z = (this->_dim == 3)?grad_w(2):0;
 
         const libMesh::Number r = u_qpoint[qp](0);
 
         libMesh::Real jac = JxW[qp];
 
         // Compute the viscosity at this qp
         libMesh::Real _mu_qp = this->_mu(context, qp);
 
         if( this->_is_axisymmetric )
           {
             jac *= r;
           }
 
         // First, an i-loop over the velocity degrees of freedom.
         // We know that n_u_dofs == n_v_dofs so we can compute contributions
         // for both at the same time.
         for (unsigned int i=0; i != n_u_dofs; i++)
           {
             Fu(i) += jac *
               (-this->_rho*u_phi[i][qp]*(U*grad_u)        // convection term
                +p*u_gradphi[i][qp](0)              // pressure term
                -_mu_qp*(u_gradphi[i][qp]*grad_u) ); // diffusion term
 
             if( this->_is_axisymmetric )
               {
                 Fu(i) += u_phi[i][qp]*( p/r - _mu_qp*U(0)/(r*r) )*jac;
               }
 
             Fv(i) += jac *
               (-this->_rho*u_phi[i][qp]*(U*grad_v)        // convection term
                +p*u_gradphi[i][qp](1)              // pressure term
                -_mu_qp*(u_gradphi[i][qp]*grad_v) ); // diffusion term
 
             if (this->_dim == 3)
               {
                 (*Fw)(i) += jac *
                   (-this->_rho*u_phi[i][qp]*(U*grad_w)        // convection term
                    +p*u_gradphi[i][qp](2)              // pressure term
                    -_mu_qp*(u_gradphi[i][qp]*grad_w) ); // diffusion term
               }
 
             if (compute_jacobian)
               {
                 for (unsigned int j=0; j != n_u_dofs; j++)
                   {
                     // TODO: precompute some terms like:
                     //   (Uvec*u_gradphi[j][qp]),
                     //   u_phi[i][qp]*u_phi[j][qp],
                     //   (u_gradphi[i][qp]*u_gradphi[j][qp])
 
                     Kuu(i,j) += jac * context.get_elem_solution_derivative() *
                       (-this->_rho*u_phi[i][qp]*(U*u_gradphi[j][qp])       // convection term
                        -this->_rho*u_phi[i][qp]*grad_u_x*u_phi[j][qp]             // convection term
                        -_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp])); // diffusion term
 
                     
                     if( this->_is_axisymmetric )
                       {
                         Kuu(i,j) -= u_phi[i][qp]*_mu_qp*u_phi[j][qp]/(r*r)*jac * context.get_elem_solution_derivative();
                       }
 
                     Kuv(i,j) += jac * context.get_elem_solution_derivative() *
                       (-this->_rho*u_phi[i][qp]*grad_u_y*u_phi[j][qp]);           // convection term
 
                     Kvv(i,j) += jac * context.get_elem_solution_derivative() *
                       (-this->_rho*u_phi[i][qp]*(U*u_gradphi[j][qp])       // convection term
                        -this->_rho*u_phi[i][qp]*grad_v_y*u_phi[j][qp]             // convection term
                        -_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp])); // diffusion term
 
                     Kvu(i,j) += jac * context.get_elem_solution_derivative() *
                       (-this->_rho*u_phi[i][qp]*grad_v_x*u_phi[j][qp]);           // convection term
 
                     if (this->_dim == 3)
                       {
                         (*Kuw)(i,j) += jac * context.get_elem_solution_derivative() *
                           (-this->_rho*u_phi[i][qp]*grad_u_z*u_phi[j][qp]);           // convection term
 
                         (*Kvw)(i,j) += jac * context.get_elem_solution_derivative() *
                           (-this->_rho*u_phi[i][qp]*grad_v_z*u_phi[j][qp]);           // convection term
 
                         (*Kww)(i,j) += jac * context.get_elem_solution_derivative() *
                           (-this->_rho*u_phi[i][qp]*(U*u_gradphi[j][qp])       // convection term
                            -this->_rho*u_phi[i][qp]*grad_w_z*u_phi[j][qp]             // convection term
                            -_mu_qp*(u_gradphi[i][qp]*u_gradphi[j][qp])); // diffusion term
                         (*Kwu)(i,j) += jac * context.get_elem_solution_derivative() *
                           (-this->_rho*u_phi[i][qp]*grad_w_x*u_phi[j][qp]);           // convection term
                         (*Kwv)(i,j) += jac * context.get_elem_solution_derivative() *
                           (-this->_rho*u_phi[i][qp]*grad_w_y*u_phi[j][qp]);           // convection term
                       }
                   } // end of the inner dof (j) loop
 
                 // Matrix contributions for the up, vp and wp couplings
                 for (unsigned int j=0; j != n_p_dofs; j++)
                   {
                     Kup(i,j) += u_gradphi[i][qp](0)*p_phi[j][qp]*jac * context.get_elem_solution_derivative();
                     Kvp(i,j) += u_gradphi[i][qp](1)*p_phi[j][qp]*jac * context.get_elem_solution_derivative();
 
                     if (this->_dim == 3)
                       {
                         (*Kwp)(i,j) += u_gradphi[i][qp](2)*p_phi[j][qp]*jac * context.get_elem_solution_derivative();
                       }
 
                     if( this->_is_axisymmetric )
                       {
                         Kup(i,j) += u_phi[i][qp]*p_phi[j][qp]/r*jac * context.get_elem_solution_derivative();
                       }
 
                   } // end of the inner dof (j) loop
 
                 
 
               } // end - if (compute_jacobian)
 
           } // end of the outer dof (i) loop
       } // end of the quadrature point (qp) loop
 
 #ifdef GRINS_USE_GRVY_TIMERS
     this->_timer->EndTimer("IncompressibleNavierStokes::element_time_derivative");
 #endif
 
     return;
   }

bool GRINS::Physics::enabled_on_elem ( const libMesh::Elem * elem )

virtualinherited

Find if current physics is active on supplied element.

Definition at line 83 of file physics.C.

References GRINS::Physics::_enabled_subdomains.

   {
     // Check if enabled_subdomains flag has been set and if we're
     // looking at a real element (rather than a nonlocal evaluation)
     if( !elem || _enabled_subdomains.empty() ) 
       return true;
   
     // Check if current physics is enabled on elem
     if( _enabled_subdomains.find( elem->subdomain_id() ) == _enabled_subdomains.end() )
       return false;
 
     return true;
   }

BCHandlingBase * GRINS::Physics::get_bc_handler ( )

inlineinherited

Definition at line 282 of file physics.h.

References GRINS::Physics::_bc_handler.

   {
     return _bc_handler;
   }

ICHandlingBase * GRINS::Physics::get_ic_handler ( )

inlineinherited

Definition at line 288 of file physics.h.

References GRINS::Physics::_ic_handler.

   {
     return _ic_handler;
   }

void GRINS::Physics::init_bcs ( libMesh::FEMSystem * system )

inherited

Definition at line 118 of file physics.C.

References GRINS::Physics::_bc_handler, GRINS::BCHandlingBase::init_bc_data(), GRINS::BCHandlingBase::init_dirichlet_bc_func_objs(), GRINS::BCHandlingBase::init_dirichlet_bcs(), and GRINS::BCHandlingBase::init_periodic_bcs().

   {
     // Only need to init BC's if the physics actually created a handler
     if( _bc_handler )
       {
         _bc_handler->init_bc_data( *system );
         _bc_handler->init_dirichlet_bcs( system );
         _bc_handler->init_dirichlet_bc_func_objs( system );
         _bc_handler->init_periodic_bcs( system );
       }
 
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokesBase< Mu >::init_context ( AssemblyContext & context )

virtualinherited

Initialize context for added physics variables.

Reimplemented from GRINS::Physics.

Reimplemented in GRINS::AveragedTurbine< Viscosity >, GRINS::VelocityDrag< Viscosity >, GRINS::AveragedFan< Viscosity >, GRINS::ParsedVelocitySource< Viscosity >, GRINS::VelocityPenalty< Viscosity >, GRINS::AveragedTurbineAdjointStabilization< Viscosity >, GRINS::VelocityDragAdjointStabilization< Viscosity >, GRINS::AveragedFanAdjointStabilization< Viscosity >, GRINS::ParsedVelocitySourceAdjointStabilization< Viscosity >, GRINS::VelocityPenaltyAdjointStabilization< Viscosity >, and GRINS::IncompressibleNavierStokesStabilizationBase< Viscosity >.

Definition at line 91 of file inc_navier_stokes_base.C.

Referenced by GRINS::IncompressibleNavierStokesStabilizationBase< Viscosity >::init_context().

   {
     // We should prerequest all the data
     // we will need to build the linear system
     // or evaluate a quantity of interest.
     context.get_element_fe(_flow_vars.u_var())->get_JxW();
     context.get_element_fe(_flow_vars.u_var())->get_phi();
     context.get_element_fe(_flow_vars.u_var())->get_dphi();
     context.get_element_fe(_flow_vars.u_var())->get_xyz();
 
     context.get_element_fe(_flow_vars.p_var())->get_phi();
     context.get_element_fe(_flow_vars.p_var())->get_xyz();
 
     context.get_side_fe(_flow_vars.u_var())->get_JxW();
     context.get_side_fe(_flow_vars.u_var())->get_phi();
     context.get_side_fe(_flow_vars.u_var())->get_dphi();
     context.get_side_fe(_flow_vars.u_var())->get_xyz();
 
     return;
   }

void GRINS::Physics::init_ics	(	libMesh::FEMSystem *	system,
		libMesh::CompositeFunction< libMesh::Number > &	all_ics
	)

inherited

Definition at line 133 of file physics.C.

References GRINS::Physics::_ic_handler, and GRINS::ICHandlingBase::init_ic_data().

   {
     if( _ic_handler )
       {
         _ic_handler->init_ic_data( *system, all_ics );
       }
 
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokesBase< Mu >::init_variables ( libMesh::FEMSystem * system )

virtualinherited

Initialization of Navier-Stokes variables.

Add velocity and pressure variables to system.

Implements GRINS::Physics.

Reimplemented in GRINS::AveragedTurbineBase< Viscosity >, and GRINS::IncompressibleNavierStokesStabilizationBase< Viscosity >.

Definition at line 63 of file inc_navier_stokes_base.C.

Referenced by GRINS::IncompressibleNavierStokesStabilizationBase< Viscosity >::init_variables(), and GRINS::AveragedTurbineBase< Viscosity >::init_variables().

   {
     this->_dim = system->get_mesh().mesh_dimension();
 
     this->_flow_vars.init(system);
     
     this->_mu.init(system); 
    
     return;
   }

bool GRINS::Physics::is_steady ( ) const

inherited

Returns whether or not this physics is being solved with a steady solver.

Definition at line 103 of file physics.C.

References GRINS::Physics::_is_steady.

Referenced by GRINS::Physics::set_is_steady().

   {
     return _is_steady;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokes< Mu >::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 482 of file inc_navier_stokes.C.

   {
     // Element Jacobian * quadrature weights for interior integration
     // We assume the same for each flow variable
     const std::vector<libMesh::Real> &JxW = 
       context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
 
     // The shape functions at interior quadrature points.
     // We assume the same for each flow variable
     const std::vector<std::vector<libMesh::Real> >& u_phi = 
       context.get_element_fe(this->_flow_vars.u_var())->get_phi();
 
     const std::vector<libMesh::Point>& u_qpoint = 
       context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
 
     // The number of local degrees of freedom in each variable
     const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
 
     // The subvectors and submatrices we need to fill:
     libMesh::DenseSubVector<libMesh::Real> &F_u = context.get_elem_residual(this->_flow_vars.u_var());
     libMesh::DenseSubVector<libMesh::Real> &F_v = context.get_elem_residual(this->_flow_vars.v_var());
     libMesh::DenseSubVector<libMesh::Real>* F_w = NULL;
 
     libMesh::DenseSubMatrix<libMesh::Real> &M_uu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
     libMesh::DenseSubMatrix<libMesh::Real> &M_vv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
     libMesh::DenseSubMatrix<libMesh::Real>* M_ww = NULL;
 
     if( this->_dim == 3 )
       {
         F_w  = &context.get_elem_residual(this->_flow_vars.w_var()); // R_{w}
         M_ww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
       }
 
     unsigned int n_qpoints = context.get_element_qrule().n_points();
 
     for (unsigned int qp = 0; qp != n_qpoints; ++qp)
       {
         // For the mass residual, we need to be a little careful.
         // The time integrator is handling the time-discretization
         // for us so we need to supply M(u_fixed)*u' for the residual.
         // u_fixed will be given by the fixed_interior_value function
         // while u' will be given by the interior_rate function.
         libMesh::Real u_dot, v_dot, w_dot = 0.0; 
         context.interior_rate(this->_flow_vars.u_var(), qp, u_dot);
         context.interior_rate(this->_flow_vars.v_var(), qp, v_dot);
 
         if(this->_dim == 3 )
           context.interior_rate(this->_flow_vars.w_var(), qp, w_dot);
       
         const libMesh::Number r = u_qpoint[qp](0);
 
         libMesh::Real jac = JxW[qp];
 
         if( this->_is_axisymmetric )
           {
             jac *= r;
           }
 
         for (unsigned int i = 0; i != n_u_dofs; ++i)
           {
             F_u(i) -= this->_rho*u_dot*u_phi[i][qp]*jac;
             F_v(i) -= this->_rho*v_dot*u_phi[i][qp]*jac;
 
             if( this->_dim == 3 )
               (*F_w)(i) -= this->_rho*w_dot*u_phi[i][qp]*jac;
           
             if( compute_jacobian )
               {
                 for (unsigned int j=0; j != n_u_dofs; j++)
                   {
                     // Assuming rho is constant w.r.t. u, v, w
                     // and T (if Boussinesq added).
                     libMesh::Real value = this->_rho*u_phi[i][qp]*u_phi[j][qp]*jac * context.get_elem_solution_rate_derivative();
 
                     M_uu(i,j) -= value;
                     M_vv(i,j) -= value;
 
                     if( this->_dim == 3)
                       {
                         (*M_ww)(i,j) -= value;
                       }
 
                   } // End dof loop
               } // End Jacobian check
           } // End dof loop
       } // End quadrature loop
 
     return;
   }

void GRINS::Physics::nonlocal_constraint	(	bool	compute_jacobian,
		AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Constraint part(s) of physics for scalar variables.

Reimplemented in GRINS::ScalarODE.

Definition at line 250 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::nonlocal_constraint().

   {
     return;
   }   

void GRINS::Physics::nonlocal_mass_residual	(	bool	compute_jacobian,
		AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Mass matrix part(s) for scalar variables.

Reimplemented in GRINS::ScalarODE, and GRINS::AveragedTurbine< Viscosity >.

Definition at line 264 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::nonlocal_mass_residual().

   {
     return;
   }

void GRINS::Physics::nonlocal_time_derivative	(	bool	compute_jacobian,
		AssemblyContext &	context,
		CachedValues &	cache
	)

virtualinherited

Time dependent part(s) of physics for scalar variables.

Reimplemented in GRINS::AveragedTurbine< Viscosity >, and GRINS::ScalarODE.

Definition at line 229 of file physics.C.

Referenced by GRINS::MultiphysicsSystem::nonlocal_time_derivative().

   {
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokes< Mu >::read_input_options ( const GetPot & input )

virtual

Read options from GetPot input file.

Reimplemented from GRINS::Physics.

Definition at line 72 of file inc_navier_stokes.C.

References GRINS::incompressible_navier_stokes.

Referenced by GRINS::IncompressibleNavierStokes< Viscosity >::IncompressibleNavierStokes().

   {
     // Other quantities read in base class
 
     // Read pressure pinning information
     this->_pin_pressure = input("Physics/"+incompressible_navier_stokes+"/pin_pressure", false );
 
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokesBase< Mu >::register_parameter	(	const std::string &	param_name,
		libMesh::ParameterMultiPointer< libMesh::Number > &	param_pointer
	)		const

virtualinherited

Each subclass will register its copy of an independent.

Reimplemented from GRINS::ParameterUser.

Definition at line 114 of file inc_navier_stokes_base.C.

References GRINS::ParameterUser::register_parameter().

   {
     ParameterUser::register_parameter(param_name, param_pointer);
     _mu.register_parameter(param_name, param_pointer);
   }

template<class Mu >

void GRINS::IncompressibleNavierStokes< Mu >::register_postprocessing_vars	(	const GetPot &	input,
		PostProcessedQuantities< libMesh::Real > &	postprocessing
	)

virtual

Register postprocessing variables for IncompressibleNavierStokes.

Reimplemented from GRINS::Physics.

Definition at line 83 of file inc_navier_stokes.C.

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

   {
     std::string section = "Physics/"+incompressible_navier_stokes+"/output_vars";
 
     if( input.have_variable(section) )
       {
         unsigned int n_vars = input.vector_variable_size(section);
 
         for( unsigned int v = 0; v < n_vars; v++ )
           {
             std::string name = input(section,"DIE!",v);
 
             if( name == std::string("mu") )
               {
                 this->_mu_index = postprocessing.register_quantity( name );
               }
             else
               {
                 std::cerr << "Error: Invalid output_vars value for "+incompressible_navier_stokes << std::endl
                           << "       Found " << name << std::endl
                           << "       Acceptable values are: mu" << std::endl;
                 libmesh_error();
               }
           }
       }
 
     return;
   }

void GRINS::Physics::set_is_steady ( bool is_steady )

inherited

Sets whether this physics is to be solved with a steady solver or not.

Since the member variable is static, only needs to be called on a single physics.

Definition at line 97 of file physics.C.

References GRINS::Physics::_is_steady, and GRINS::Physics::is_steady().

   {
     _is_steady = is_steady;
     return;
   }

void GRINS::ParameterUser::set_parameter	(	libMesh::Number &	param_variable,
		const GetPot &	input,
		const std::string &	param_name,
		libMesh::Number	param_default
	)

virtualinherited

Each subclass can simultaneously read a parameter value from.

Definition at line 35 of file parameter_user.C.

References GRINS::ParameterUser::_my_name, and GRINS::ParameterUser::_my_parameters.

Referenced by GRINS::AveragedFanAdjointStabilization< Viscosity >::AveragedFanAdjointStabilization(), GRINS::AveragedTurbineAdjointStabilization< Viscosity >::AveragedTurbineAdjointStabilization(), GRINS::BoussinesqBuoyancyAdjointStabilization< Viscosity >::BoussinesqBuoyancyAdjointStabilization(), GRINS::BoussinesqBuoyancyBase::BoussinesqBuoyancyBase(), GRINS::BoussinesqBuoyancySPGSMStabilization< Viscosity >::BoussinesqBuoyancySPGSMStabilization(), GRINS::ConstantConductivity::ConstantConductivity(), GRINS::ConstantPrandtlConductivity::ConstantPrandtlConductivity(), GRINS::ConstantSourceFunction::ConstantSourceFunction(), GRINS::ConstantSourceTerm::ConstantSourceTerm(), GRINS::ConstantSpecificHeat::ConstantSpecificHeat(), GRINS::ConstantViscosity::ConstantViscosity(), GRINS::ElasticCable< StressStrainLaw >::ElasticCable(), GRINS::ElasticCableConstantGravity::ElasticCableConstantGravity(), GRINS::ElasticMembrane< StressStrainLaw >::ElasticMembrane(), GRINS::ElasticMembraneConstantPressure::ElasticMembraneConstantPressure(), GRINS::HeatConduction< Conductivity >::HeatConduction(), GRINS::HeatTransferBase< Conductivity >::HeatTransferBase(), GRINS::IncompressibleNavierStokesBase< Viscosity >::IncompressibleNavierStokesBase(), GRINS::AverageNusseltNumber::init(), GRINS::MooneyRivlin::MooneyRivlin(), GRINS::ReactingLowMachNavierStokesBase< Mixture, Evaluator >::ReactingLowMachNavierStokesBase(), GRINS::HookesLaw1D::read_input_options(), GRINS::HookesLaw::read_input_options(), GRINS::AxisymmetricBoussinesqBuoyancy::read_input_options(), and GRINS::VelocityDragAdjointStabilization< Viscosity >::VelocityDragAdjointStabilization().

   {
     param_variable = input(param_name, param_default);
 
     libmesh_assert_msg(!_my_parameters.count(param_name),
       "ERROR: " << _my_name << " double-registered parameter " <<
       param_name);
 
     _my_parameters[param_name] = &param_variable;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokesBase< Mu >::set_time_evolving_vars ( libMesh::FEMSystem * system )

virtualinherited

Sets velocity variables to be time-evolving.

Reimplemented from GRINS::Physics.

Reimplemented in GRINS::AveragedTurbineBase< Viscosity >, and GRINS::ParsedVelocitySourceBase< Viscosity >.

Definition at line 75 of file inc_navier_stokes_base.C.

Referenced by GRINS::AveragedTurbineBase< Viscosity >::set_time_evolving_vars().

   {
     const unsigned int dim = system->get_mesh().mesh_dimension();
 
     // Tell the system to march velocity forward in time, but
     // leave p as a constraint only
     system->time_evolving(_flow_vars.u_var());
     system->time_evolving(_flow_vars.v_var());
 
     if (dim == 3)
       system->time_evolving(_flow_vars.w_var());
 
     return;
   }

template<class Mu >

void GRINS::IncompressibleNavierStokes< Mu >::side_constraint	(	bool	compute_jacobian,
		AssemblyContext &	context,
		CachedValues &	cache
	)

virtual

Constraint part(s) of physics for boundaries of elements on the domain boundary.

Reimplemented from GRINS::Physics.

Definition at line 468 of file inc_navier_stokes.C.

   {
     // Pin p = p_value at p_point
     if( _pin_pressure )
       {
         _p_pinning.pin_value( context, compute_jacobian, this->_flow_vars.p_var() );
       }
 
     return;
   }