GRINS-0.7.0
List of all members | Public Member Functions | Private Member Functions
GRINS::IncompressibleNavierStokesAdjointStabilization< Viscosity > Class Template Reference

Adds VMS-based stabilization to LowMachNavierStokes physics class. More...

#include <inc_navier_stokes_adjoint_stab.h>

Inheritance diagram for GRINS::IncompressibleNavierStokesAdjointStabilization< Viscosity >:
Inheritance graph
[legend]
Collaboration diagram for GRINS::IncompressibleNavierStokesAdjointStabilization< Viscosity >:
Collaboration graph
[legend]

Public Member Functions

 IncompressibleNavierStokesAdjointStabilization (const GRINS::PhysicsName &physics_name, const GetPot &input)
 
virtual ~IncompressibleNavierStokesAdjointStabilization ()
 
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...
 
- Public Member Functions inherited from GRINS::IncompressibleNavierStokesStabilizationBase< Viscosity >
 IncompressibleNavierStokesStabilizationBase (const GRINS::PhysicsName &physics_name, const GetPot &input)
 
virtual ~IncompressibleNavierStokesStabilizationBase ()
 
virtual void init_context (AssemblyContext &context)
 Initialize context for added physics variables. More...
 
- Public Member Functions inherited from GRINS::IncompressibleNavierStokesBase< Viscosity >
 IncompressibleNavierStokesBase (const std::string &my_physics_name, const std::string &core_physics_name, const GetPot &input)
 
 ~IncompressibleNavierStokesBase ()
 
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 register_parameter (const std::string &param_name, libMesh::ParameterMultiAccessor< libMesh::Number > &param_pointer) const
 Each subclass will register its copy of an independent. More...
 
libMesh::Real get_viscosity_value (AssemblyContext &context, unsigned int qp) const
 
- 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 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 name of postprocessed quantity with PostProcessedQuantities. 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_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_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)
 
virtual void compute_postprocessed_quantity (unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
 
ICHandlingBase * get_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...
 

Private Member Functions

 IncompressibleNavierStokesAdjointStabilization ()
 

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...
 
- 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 inherited from GRINS::IncompressibleNavierStokesStabilizationBase< Viscosity >
IncompressibleNavierStokesStabilizationHelper _stab_helper
 
- Protected Attributes inherited from GRINS::IncompressibleNavierStokesBase< Viscosity >
unsigned int _dim
 Physical dimension of problem. More...
 
VelocityFEVariables _flow_vars
 
PressureFEVariable _press_var
 
libMesh::Number _rho
 Material parameters, read from input. More...
 
Viscosity _mu
 Viscosity object. 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...
 
- 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 GRINS::IncompressibleNavierStokesAdjointStabilization< Viscosity >

Adds VMS-based stabilization to LowMachNavierStokes physics class.

Definition at line 35 of file inc_navier_stokes_adjoint_stab.h.

Constructor & Destructor Documentation

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

Definition at line 40 of file inc_navier_stokes_adjoint_stab.C.

42  : IncompressibleNavierStokesStabilizationBase<Mu>(physics_name,input)
43  {}
template<class Viscosity >
virtual GRINS::IncompressibleNavierStokesAdjointStabilization< Viscosity >::~IncompressibleNavierStokesAdjointStabilization ( )
inlinevirtual

Definition at line 41 of file inc_navier_stokes_adjoint_stab.h.

41 {};
template<class Viscosity >
GRINS::IncompressibleNavierStokesAdjointStabilization< Viscosity >::IncompressibleNavierStokesAdjointStabilization ( )
private

Member Function Documentation

template<class Mu >
void GRINS::IncompressibleNavierStokesAdjointStabilization< 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 376 of file inc_navier_stokes_adjoint_stab.C.

379  {
380 #ifdef GRINS_USE_GRVY_TIMERS
381  this->_timer->BeginTimer("IncompressibleNavierStokesAdjointStabilization::element_constraint");
382 #endif
383 
384  // The number of local degrees of freedom in each variable.
385  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
386  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
387 
388  // Element Jacobian * quadrature weights for interior integration.
389  const std::vector<libMesh::Real> &JxW =
390  context.get_element_fe(this->_flow_vars.u())->get_JxW();
391 
392  // The pressure shape functions at interior quadrature points.
393  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
394  context.get_element_fe(this->_press_var.p())->get_dphi();
395 
396  // The velocity shape functions at interior quadrature points.
397  const std::vector<std::vector<libMesh::Real> >& u_phi =
398  context.get_element_fe(this->_flow_vars.u())->get_phi();
399 
400  const std::vector<std::vector<libMesh::RealGradient> >& u_dphi =
401  context.get_element_fe(this->_flow_vars.u())->get_dphi();
402 
403  const std::vector<std::vector<libMesh::RealTensor> >& u_d2phi =
404  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
405 
406  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
407 
408  libMesh::DenseSubMatrix<libMesh::Number> &Kpp =
409  context.get_elem_jacobian(this->_press_var.p(), this->_press_var.p()); // J_{pp}
410  libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
411  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u()); // J_{pu}
412  libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
413  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v()); // J_{pv}
414  libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
415 
416  if(this->_dim == 3)
417  {
418  Kpw = &context.get_elem_jacobian
419  (this->_press_var.p(), this->_flow_vars.w()); // J_{pw}
420  }
421 
422  unsigned int n_qpoints = context.get_element_qrule().n_points();
423 
424  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
425 
426  for (unsigned int qp=0; qp != n_qpoints; qp++)
427  {
428  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
429  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
430 
431  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
432  context.interior_value( this->_flow_vars.v(), qp ) );
433  if( this->_dim == 3 )
434  {
435  U(2) = context.interior_value( this->_flow_vars.w(), qp );
436  }
437 
438  libMesh::Real tau_M;
439  libMesh::Real d_tau_M_d_rho;
440  libMesh::Gradient d_tau_M_dU;
441  libMesh::Gradient RM_s, d_RM_s_uvw_dgraduvw;
442  libMesh::Tensor d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dhessuvw;
443 
444  // Compute the viscosity at this qp
445  libMesh::Real _mu_qp = this->_mu(context, qp);
446 
447  if (compute_jacobian)
448  {
449  this->_stab_helper.compute_tau_momentum_and_derivs
450  ( context, qp, g, G, this->_rho, U, _mu_qp,
451  tau_M, d_tau_M_d_rho, d_tau_M_dU,
452  this->_is_steady );
453  this->_stab_helper.compute_res_momentum_steady_and_derivs
454  ( context, qp, this->_rho, _mu_qp,
455  RM_s, d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dgraduvw,
456  d_RM_s_uvw_dhessuvw);
457  }
458  else
459  {
460  tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
461  RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
462  }
463 
464  // Now a loop over the pressure degrees of freedom. This
465  // computes the contributions of the continuity equation.
466  for (unsigned int i=0; i != n_p_dofs; i++)
467  {
468  Fp(i) += tau_M*(RM_s*p_dphi[i][qp])*JxW[qp];
469 
470  if (compute_jacobian)
471  {
472  const libMesh::Real fixed_deriv =
473  context.get_fixed_solution_derivative();
474 
475  const libMesh::TypeVector<libMesh::Number>
476  p_dphiiT_d_RM_s_dgradp =
477  d_RM_s_dgradp.transpose() * p_dphi[i][qp];
478 
479  for (unsigned int j=0; j != n_p_dofs; j++)
480  {
481  Kpp(i,j) += tau_M*(p_dphiiT_d_RM_s_dgradp*p_dphi[j][qp])*JxW[qp];
482  }
483 
484  for (unsigned int j=0; j != n_u_dofs; j++)
485  {
486  Kpu(i,j) += (
487  d_tau_M_dU(0)*u_phi[j][qp]*(RM_s*p_dphi[i][qp])
488  + tau_M*
489  (
490  d_RM_s_dU(0,0)*u_phi[j][qp]*p_dphi[i][qp](0) +
491  d_RM_s_dU(1,0)*u_phi[j][qp]*p_dphi[i][qp](1) +
492  d_RM_s_dU(2,0)*u_phi[j][qp]*p_dphi[i][qp](2) +
493  d_RM_s_uvw_dgraduvw*u_dphi[j][qp]*p_dphi[i][qp](0) +
494  d_RM_s_uvw_dhessuvw.contract(u_d2phi[j][qp])*p_dphi[i][qp](0)
495  )*fixed_deriv
496  )*JxW[qp];
497  Kpv(i,j) += (
498  d_tau_M_dU(1)*u_phi[j][qp]*(RM_s*p_dphi[i][qp])
499  + tau_M*
500  (
501  d_RM_s_dU(0,1)*u_phi[j][qp]*p_dphi[i][qp](0) +
502  d_RM_s_dU(1,1)*u_phi[j][qp]*p_dphi[i][qp](1) +
503  d_RM_s_dU(2,1)*u_phi[j][qp]*p_dphi[i][qp](2) +
504  d_RM_s_uvw_dgraduvw*u_dphi[j][qp]*p_dphi[i][qp](1) +
505  d_RM_s_uvw_dhessuvw.contract(u_d2phi[j][qp])*p_dphi[i][qp](0)
506  )*fixed_deriv
507  )*JxW[qp];
508  }
509 
510  if(this->_dim == 3)
511  {
512  for (unsigned int j=0; j != n_u_dofs; j++)
513  {
514  (*Kpw)(i,j) += (
515  d_tau_M_dU(2)*u_phi[j][qp]*(RM_s*p_dphi[i][qp])
516  + tau_M*
517  (
518  d_RM_s_dU(0,2)*u_phi[j][qp]*p_dphi[i][qp](0) +
519  d_RM_s_dU(1,2)*u_phi[j][qp]*p_dphi[i][qp](1) +
520  d_RM_s_dU(2,2)*u_phi[j][qp]*p_dphi[i][qp](2) +
521  d_RM_s_uvw_dgraduvw*u_dphi[j][qp]*p_dphi[i][qp](2) +
522  d_RM_s_uvw_dhessuvw.contract(u_d2phi[j][qp])*p_dphi[i][qp](0)
523  )*fixed_deriv
524  )*JxW[qp];
525  }
526  }
527  }
528  }
529  }
530 
531 #ifdef GRINS_USE_GRVY_TIMERS
532  this->_timer->EndTimer("IncompressibleNavierStokesAdjointStabilization::element_constraint");
533 #endif
534 
535  return;
536  }
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_momentum_and_derivs
void compute_tau_momentum_and_derivs(AssemblyContext &c, unsigned int qp, libMesh::RealGradient &g, libMesh::RealTensor &G, libMesh::Real rho, libMesh::Gradient U, libMesh::Real T, libMesh::Real &tau_M, libMesh::Real &d_tau_M_d_rho, libMesh::Gradient &d_tau_M_d_U, bool is_steady) const
Definition: inc_navier_stokes_stab_helper.h:231
GRINS::StabilizationHelper::compute_G
libMesh::RealTensor compute_G(libMesh::FEBase *fe, AssemblyContext &c, unsigned int qp) const
Definition: stab_helper.C:64
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_momentum_steady
libMesh::RealGradient compute_res_momentum_steady(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, const libMesh::Real mu) const
Definition: inc_navier_stokes_stab_helper.C:212
GRINS::IncompressibleNavierStokesBase::_rho
libMesh::Number _rho
Material parameters, read from input.
Definition: inc_navier_stokes_base.h:87
GRINS::IncompressibleNavierStokesBase::_dim
unsigned int _dim
Physical dimension of problem.
Definition: inc_navier_stokes_base.h:81
GRINS::IncompressibleNavierStokesStabilizationBase::_stab_helper
IncompressibleNavierStokesStabilizationHelper _stab_helper
Definition: inc_navier_stokes_stab_base.h:50
GRINS::StabilizationHelper::compute_g
libMesh::RealGradient compute_g(libMesh::FEBase *fe, AssemblyContext &c, unsigned int qp) const
Definition: stab_helper.C:47
GRINS::PressureFEVariable::p
VariableIndex p() const
Definition: pressure_fe_variable.h:71
GRINS::Physics::_is_steady
static bool _is_steady
Caches whether or not the solver that's being used is steady or not.
Definition: physics.h:277
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_momentum_steady_and_derivs
void compute_res_momentum_steady_and_derivs(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, const libMesh::Real mu, libMesh::Gradient &res_M, libMesh::Tensor &d_res_M_dgradp, libMesh::Tensor &d_res_M_dU, libMesh::Gradient &d_res_Muvw_dgraduvw, libMesh::Tensor &d_res_Muvw_dhessuvw) const
Definition: inc_navier_stokes_stab_helper.C:250
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_momentum
libMesh::Real compute_tau_momentum(AssemblyContext &c, unsigned int qp, libMesh::RealGradient &g, libMesh::RealTensor &G, libMesh::Real rho, libMesh::Gradient U, libMesh::Real mu, bool is_steady) const
Definition: inc_navier_stokes_stab_helper.h:217
template<class Mu >
void GRINS::IncompressibleNavierStokesAdjointStabilization< Mu >::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 46 of file inc_navier_stokes_adjoint_stab.C.

49  {
50 #ifdef GRINS_USE_GRVY_TIMERS
51  this->_timer->BeginTimer("IncompressibleNavierStokesAdjointStabilization::element_time_derivative");
52 #endif
53 
54  // The number of local degrees of freedom in each variable.
55  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
56  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
57 
58  // Element Jacobian * quadrature weights for interior integration.
59  const std::vector<libMesh::Real> &JxW =
60  context.get_element_fe(this->_flow_vars.u())->get_JxW();
61 
62  const std::vector<std::vector<libMesh::RealGradient> >& p_gradphi =
63  context.get_element_fe(this->_press_var.p())->get_dphi();
64 
65  const std::vector<std::vector<libMesh::Real> >& u_phi =
66  context.get_element_fe(this->_flow_vars.u())->get_phi();
67 
68  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
69  context.get_element_fe(this->_flow_vars.u())->get_dphi();
70 
71  const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
72  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
73 
74  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{p}
75  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{p}
76  libMesh::DenseSubMatrix<libMesh::Number> &Kup =
77  context.get_elem_jacobian(this->_flow_vars.u(), this->_press_var.p()); // J_{up}
78  libMesh::DenseSubMatrix<libMesh::Number> &Kuu =
79  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u()); // J_{uu}
80  libMesh::DenseSubMatrix<libMesh::Number> &Kuv =
81  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v()); // J_{uv}
82  libMesh::DenseSubMatrix<libMesh::Number> &Kvp =
83  context.get_elem_jacobian(this->_flow_vars.v(), this->_press_var.p()); // J_{vp}
84  libMesh::DenseSubMatrix<libMesh::Number> &Kvu =
85  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u()); // J_{vu}
86  libMesh::DenseSubMatrix<libMesh::Number> &Kvv =
87  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v()); // J_{vv}
88 
89  libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
90  libMesh::DenseSubMatrix<libMesh::Number> *Kuw = NULL;
91  libMesh::DenseSubMatrix<libMesh::Number> *Kvw = NULL;
92  libMesh::DenseSubMatrix<libMesh::Number> *Kwp = NULL;
93  libMesh::DenseSubMatrix<libMesh::Number> *Kwu = NULL;
94  libMesh::DenseSubMatrix<libMesh::Number> *Kwv = NULL;
95  libMesh::DenseSubMatrix<libMesh::Number> *Kww = NULL;
96 
97  if(this->_dim == 3)
98  {
99  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
100  Kuw = &context.get_elem_jacobian
101  (this->_flow_vars.u(), this->_flow_vars.w()); // J_{uw}
102  Kvw = &context.get_elem_jacobian
103  (this->_flow_vars.v(), this->_flow_vars.w()); // J_{vw}
104  Kwp = &context.get_elem_jacobian
105  (this->_flow_vars.w(), this->_press_var.p()); // J_{wp}
106  Kwu = &context.get_elem_jacobian
107  (this->_flow_vars.w(), this->_flow_vars.u()); // J_{wu}
108  Kwv = &context.get_elem_jacobian
109  (this->_flow_vars.w(), this->_flow_vars.v()); // J_{wv}
110  Kww = &context.get_elem_jacobian
111  (this->_flow_vars.w(), this->_flow_vars.w()); // J_{ww}
112  }
113 
114  unsigned int n_qpoints = context.get_element_qrule().n_points();
115 
116  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
117 
118  for (unsigned int qp=0; qp != n_qpoints; qp++)
119  {
120  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
121  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
122 
123  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ),
124  context.interior_value( this->_flow_vars.v(), qp ) );
125  if( this->_dim == 3 )
126  {
127  U(2) = context.interior_value( this->_flow_vars.w(), qp );
128  }
129 
130  /*
131  libMesh::Gradient grad_u, grad_v, grad_w;
132  grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
133  grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
134  if (_dim == 3)
135  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
136  */
137 
138  libMesh::Real tau_M, tau_C;
139  libMesh::Real d_tau_M_d_rho, d_tau_C_d_rho;
140  libMesh::Gradient d_tau_M_dU, d_tau_C_dU;
141  libMesh::Gradient RM_s, d_RM_s_uvw_dgraduvw;
142  libMesh::Real RC;
143  libMesh::Tensor d_RC_dgradU,
144  d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dhessuvw;
145 
146  // Compute the viscosity at this qp
147  libMesh::Real _mu_qp = this->_mu(context, qp);
148 
149  if (compute_jacobian)
150  {
151  this->_stab_helper.compute_tau_momentum_and_derivs
152  ( context, qp, g, G, this->_rho, U, _mu_qp,
153  tau_M, d_tau_M_d_rho, d_tau_M_dU,
154  this->_is_steady );
155  this->_stab_helper.compute_tau_continuity_and_derivs
156  ( tau_M, d_tau_M_d_rho, d_tau_M_dU,
157  g,
158  tau_C, d_tau_C_d_rho, d_tau_C_dU );
159  this->_stab_helper.compute_res_momentum_steady_and_derivs
160  ( context, qp, this->_rho, _mu_qp,
161  RM_s, d_RM_s_dgradp, d_RM_s_dU, d_RM_s_uvw_dgraduvw,
162  d_RM_s_uvw_dhessuvw);
163  this->_stab_helper.compute_res_continuity_and_derivs
164  ( context, qp, RC, d_RC_dgradU );
165  }
166  else
167  {
168  tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, this->_is_steady );
169  tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
170  RM_s = this->_stab_helper.compute_res_momentum_steady( context, qp, this->_rho, _mu_qp );
171  RC = this->_stab_helper.compute_res_continuity( context, qp );
172  }
173 
174  for (unsigned int i=0; i != n_u_dofs; i++)
175  {
176  libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
177  _mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
178  libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
179 
180  //libMesh::RealGradient zeroth_order_term = - this->_rho*u_phi[i][qp]*(grad_u + grad_v + grad_w);
181 
182  Fu(i) += ( -tau_M*RM_s(0)*test_func - tau_C*RC*u_gradphi[i][qp](0) )*JxW[qp];
183 
184  Fv(i) += ( -tau_M*RM_s(1)*test_func - tau_C*RC*u_gradphi[i][qp](1) )*JxW[qp];
185 
186  if(this->_dim == 3)
187  {
188  (*Fw)(i) += ( -tau_M*RM_s(2)*test_func - tau_C*RC*u_gradphi[i][qp](2) )*JxW[qp];
189  }
190 
191  if (compute_jacobian)
192  {
193  const libMesh::Real fixed_deriv =
194  context.get_fixed_solution_derivative();
195  for (unsigned int j=0; j != n_p_dofs; j++)
196  {
197  Kup(i,j) += ( -tau_M*
198  ( d_RM_s_dgradp(0,0) *
199  p_gradphi[j][qp](0) +
200  d_RM_s_dgradp(0,1) *
201  p_gradphi[j][qp](1) +
202  d_RM_s_dgradp(0,2) *
203  p_gradphi[j][qp](2)
204  )*test_func)*fixed_deriv*JxW[qp];
205  Kvp(i,j) += ( -tau_M*
206  ( d_RM_s_dgradp(1,0) *
207  p_gradphi[j][qp](0) +
208  d_RM_s_dgradp(1,1) *
209  p_gradphi[j][qp](1) +
210  d_RM_s_dgradp(1,2) *
211  p_gradphi[j][qp](2)
212  )*test_func)*fixed_deriv*JxW[qp];
213  }
214  for (unsigned int j=0; j != n_u_dofs; j++)
215  {
216  Kuu(i,j) += ( -d_tau_M_dU(0)*RM_s(0)*test_func
217  -tau_M*d_RM_s_dU(0,0)*test_func
218  - d_tau_C_dU(0)*RC*u_gradphi[i][qp](0)
219  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
220  Kuu(i,j) += ( -tau_M*RM_s(0)*d_test_func_dU(0)
221  )*u_phi[j][qp]*JxW[qp];
222  Kuu(i,j) += ( - tau_C*
223  (
224  d_RC_dgradU(0,0)*u_gradphi[j][qp](0) +
225  d_RC_dgradU(0,1)*u_gradphi[j][qp](1) +
226  d_RC_dgradU(0,2)*u_gradphi[j][qp](2)
227  )
228  *fixed_deriv*u_gradphi[i][qp](0)
229  )*JxW[qp];
230  Kuu(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dgraduvw*u_gradphi[j][qp])
231  )*fixed_deriv*JxW[qp];
232  Kuu(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dhessuvw.contract(u_hessphi[j][qp]))
233  )*fixed_deriv*JxW[qp];
234  Kuv(i,j) += ( -d_tau_M_dU(1)*RM_s(0)*test_func
235  -tau_M*d_RM_s_dU(0,1)*test_func
236  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
237  Kuv(i,j) += ( -tau_M*RM_s(0)*d_test_func_dU(1)
238  )*u_phi[j][qp]*JxW[qp];
239  Kuv(i,j) += ( - tau_C*
240  (
241  d_RC_dgradU(1,0)*u_gradphi[j][qp](0) +
242  d_RC_dgradU(1,1)*u_gradphi[j][qp](1) +
243  d_RC_dgradU(1,2)*u_gradphi[j][qp](2)
244  )
245  *fixed_deriv*u_gradphi[i][qp](0)
246  )*JxW[qp];
247  Kvu(i,j) += ( -d_tau_M_dU(0)*RM_s(1)*test_func
248  -tau_M*d_RM_s_dU(1,0)*test_func
249  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
250  Kvu(i,j) += ( -tau_M*RM_s(1)*d_test_func_dU(0)
251  )*u_phi[j][qp]*JxW[qp];
252  Kvu(i,j) += ( - tau_C*
253  (
254  d_RC_dgradU(0,0)*u_gradphi[j][qp](0) +
255  d_RC_dgradU(0,1)*u_gradphi[j][qp](1) +
256  d_RC_dgradU(0,2)*u_gradphi[j][qp](2)
257  )
258  *fixed_deriv*u_gradphi[i][qp](1)
259  )*JxW[qp];
260  Kvv(i,j) += ( -d_tau_M_dU(1)*RM_s(1)*test_func
261  -tau_M*d_RM_s_dU(1,1)*test_func
262  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
263  Kvv(i,j) += ( -tau_M*RM_s(1)*d_test_func_dU(1)
264  )*u_phi[j][qp]*JxW[qp];
265  Kvv(i,j) += ( - tau_C*
266  (
267  d_RC_dgradU(1,0)*u_gradphi[j][qp](0) +
268  d_RC_dgradU(1,1)*u_gradphi[j][qp](1) +
269  d_RC_dgradU(1,2)*u_gradphi[j][qp](2)
270  )
271  *fixed_deriv*u_gradphi[i][qp](1)
272  )*JxW[qp];
273  Kvv(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dgraduvw*u_gradphi[j][qp])
274  )*fixed_deriv*JxW[qp];
275  Kvv(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dhessuvw.contract(u_hessphi[j][qp]))
276  )*fixed_deriv*JxW[qp];
277  }
278  if(this->_dim == 3)
279  {
280  for (unsigned int j=0; j != n_p_dofs; j++)
281  {
282  (*Kwp)(i,j) += ( -tau_M*
283  ( d_RM_s_dgradp(2,0) *
284  p_gradphi[j][qp](0) +
285  d_RM_s_dgradp(2,1) *
286  p_gradphi[j][qp](1) +
287  d_RM_s_dgradp(2,2) *
288  p_gradphi[j][qp](2)
289  )*test_func)*fixed_deriv*JxW[qp];
290  }
291  for (unsigned int j=0; j != n_u_dofs; j++)
292  {
293  (*Kuw)(i,j) += ( -d_tau_M_dU(2)*RM_s(0)*test_func
294  -tau_M*d_RM_s_dU(0,2)*test_func
295  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
296  (*Kuw)(i,j) += ( -tau_M*RM_s(0)*d_test_func_dU(2)
297  )*u_phi[j][qp]*JxW[qp];
298  (*Kuw)(i,j) += ( - tau_C*
299  (
300  d_RC_dgradU(2,0)*u_gradphi[j][qp](0) +
301  d_RC_dgradU(2,1)*u_gradphi[j][qp](1) +
302  d_RC_dgradU(2,2)*u_gradphi[j][qp](2)
303  )
304  *fixed_deriv* u_gradphi[i][qp](0)
305  )*JxW[qp];
306  (*Kvw)(i,j) += ( -d_tau_M_dU(2)*RM_s(1)*test_func
307  -tau_M*d_RM_s_dU(1,2)*test_func
308  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
309  (*Kvw)(i,j) += ( -tau_M*RM_s(1)*d_test_func_dU(2)
310  )*u_phi[j][qp]*JxW[qp];
311  (*Kvw)(i,j) += ( - tau_C*
312  (
313  d_RC_dgradU(2,0)*u_gradphi[j][qp](0) +
314  d_RC_dgradU(2,1)*u_gradphi[j][qp](1) +
315  d_RC_dgradU(2,2)*u_gradphi[j][qp](2)
316  )
317  *fixed_deriv* u_gradphi[i][qp](1)
318  )*JxW[qp];
319  (*Kwu)(i,j) += ( -d_tau_M_dU(0)*RM_s(2)*test_func
320  -tau_M*d_RM_s_dU(2,0)*test_func
321  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
322  (*Kwu)(i,j) += ( -tau_M*RM_s(2)*d_test_func_dU(0)
323  )*u_phi[j][qp]*JxW[qp];
324  (*Kwu)(i,j) += ( - tau_C*
325  (
326  d_RC_dgradU(0,0)*u_gradphi[j][qp](0) +
327  d_RC_dgradU(0,1)*u_gradphi[j][qp](1) +
328  d_RC_dgradU(0,2)*u_gradphi[j][qp](2)
329  )
330  *fixed_deriv* u_gradphi[i][qp](2)
331  )*JxW[qp];
332  (*Kwv)(i,j) += ( -d_tau_M_dU(1)*RM_s(2)*test_func
333  -tau_M*d_RM_s_dU(2,1)*test_func
334  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
335  (*Kwv)(i,j) += ( -tau_M*RM_s(2)*d_test_func_dU(1)
336  )*u_phi[j][qp]*JxW[qp];
337  (*Kwv)(i,j) += ( - tau_C*
338  (
339  d_RC_dgradU(1,0)*u_gradphi[j][qp](0) +
340  d_RC_dgradU(1,1)*u_gradphi[j][qp](1) +
341  d_RC_dgradU(1,2)*u_gradphi[j][qp](2)
342  )
343  *fixed_deriv* u_gradphi[i][qp](2)
344  )*JxW[qp];
345  (*Kww)(i,j) += ( -d_tau_M_dU(2)*RM_s(2)*test_func
346  -tau_M*d_RM_s_dU(2,2)*test_func
347  )*fixed_deriv*u_phi[j][qp]*JxW[qp];
348  (*Kww)(i,j) += ( -tau_M*RM_s(2)*d_test_func_dU(2)
349  )*u_phi[j][qp]*JxW[qp];
350  (*Kww)(i,j) += ( - tau_C*
351  (
352  d_RC_dgradU(2,0)*u_gradphi[j][qp](0) +
353  d_RC_dgradU(2,1)*u_gradphi[j][qp](1) +
354  d_RC_dgradU(2,2)*u_gradphi[j][qp](2)
355  )
356  *fixed_deriv* u_gradphi[i][qp](2)
357  )*JxW[qp];
358  (*Kww)(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dgraduvw*u_gradphi[j][qp])
359  )*fixed_deriv*JxW[qp];
360  (*Kww)(i,j) += ( -tau_M*test_func*(d_RM_s_uvw_dhessuvw.contract(u_hessphi[j][qp]))
361  )*fixed_deriv*JxW[qp];
362 
363  }
364  }
365  }
366  }
367  }
368 
369 #ifdef GRINS_USE_GRVY_TIMERS
370  this->_timer->EndTimer("IncompressibleNavierStokesAdjointStabilization::element_time_derivative");
371 #endif
372  return;
373  }
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_momentum_and_derivs
void compute_tau_momentum_and_derivs(AssemblyContext &c, unsigned int qp, libMesh::RealGradient &g, libMesh::RealTensor &G, libMesh::Real rho, libMesh::Gradient U, libMesh::Real T, libMesh::Real &tau_M, libMesh::Real &d_tau_M_d_rho, libMesh::Gradient &d_tau_M_d_U, bool is_steady) const
Definition: inc_navier_stokes_stab_helper.h:231
GRINS::StabilizationHelper::compute_G
libMesh::RealTensor compute_G(libMesh::FEBase *fe, AssemblyContext &c, unsigned int qp) const
Definition: stab_helper.C:64
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_continuity
libMesh::Real compute_tau_continuity(libMesh::Real tau_C, libMesh::RealGradient &g) const
Definition: inc_navier_stokes_stab_helper.h:193
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_momentum_steady
libMesh::RealGradient compute_res_momentum_steady(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, const libMesh::Real mu) const
Definition: inc_navier_stokes_stab_helper.C:212
GRINS::IncompressibleNavierStokesBase::_rho
libMesh::Number _rho
Material parameters, read from input.
Definition: inc_navier_stokes_base.h:87
GRINS::IncompressibleNavierStokesBase::_dim
unsigned int _dim
Physical dimension of problem.
Definition: inc_navier_stokes_base.h:81
GRINS::IncompressibleNavierStokesStabilizationBase::_stab_helper
IncompressibleNavierStokesStabilizationHelper _stab_helper
Definition: inc_navier_stokes_stab_base.h:50
GRINS::StabilizationHelper::compute_g
libMesh::RealGradient compute_g(libMesh::FEBase *fe, AssemblyContext &c, unsigned int qp) const
Definition: stab_helper.C:47
GRINS::PressureFEVariable::p
VariableIndex p() const
Definition: pressure_fe_variable.h:71
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_continuity_and_derivs
void compute_tau_continuity_and_derivs(libMesh::Real tau_M, libMesh::Real d_tau_M_d_rho, libMesh::Gradient d_tau_M_d_U, libMesh::RealGradient &g, libMesh::Real &tau_C, libMesh::Real &d_tau_C_d_rho, libMesh::Gradient &d_tau_C_d_U) const
Definition: inc_navier_stokes_stab_helper.h:201
GRINS::Physics::_is_steady
static bool _is_steady
Caches whether or not the solver that's being used is steady or not.
Definition: physics.h:277
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_momentum_steady_and_derivs
void compute_res_momentum_steady_and_derivs(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, const libMesh::Real mu, libMesh::Gradient &res_M, libMesh::Tensor &d_res_M_dgradp, libMesh::Tensor &d_res_M_dU, libMesh::Gradient &d_res_Muvw_dgraduvw, libMesh::Tensor &d_res_Muvw_dhessuvw) const
Definition: inc_navier_stokes_stab_helper.C:250
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_momentum
libMesh::Real compute_tau_momentum(AssemblyContext &c, unsigned int qp, libMesh::RealGradient &g, libMesh::RealTensor &G, libMesh::Real rho, libMesh::Gradient U, libMesh::Real mu, bool is_steady) const
Definition: inc_navier_stokes_stab_helper.h:217
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_continuity_and_derivs
void compute_res_continuity_and_derivs(AssemblyContext &context, unsigned int qp, libMesh::Real &res_C, libMesh::Tensor &d_res_C_dgradU) const
Definition: inc_navier_stokes_stab_helper.C:187
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_continuity
libMesh::Real compute_res_continuity(AssemblyContext &context, unsigned int qp) const
Definition: inc_navier_stokes_stab_helper.C:168
template<class Mu >
void GRINS::IncompressibleNavierStokesAdjointStabilization< 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 539 of file inc_navier_stokes_adjoint_stab.C.

542  {
543 #ifdef GRINS_USE_GRVY_TIMERS
544  this->_timer->BeginTimer("IncompressibleNavierStokesAdjointStabilization::mass_residual");
545 #endif
546 
547  // The number of local degrees of freedom in each variable.
548  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
549  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
550 
551  // Element Jacobian * quadrature weights for interior integration.
552  const std::vector<libMesh::Real> &JxW =
553  context.get_element_fe(this->_flow_vars.u())->get_JxW();
554 
555  // The pressure shape functions at interior quadrature points.
556  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
557  context.get_element_fe(this->_press_var.p())->get_dphi();
558 
559  const std::vector<std::vector<libMesh::Real> >& u_phi =
560  context.get_element_fe(this->_flow_vars.u())->get_phi();
561 
562  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
563  context.get_element_fe(this->_flow_vars.u())->get_dphi();
564 
565  const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
566  context.get_element_fe(this->_flow_vars.u())->get_d2phi();
567 
568  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{p}
569  libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v()); // R_{p}
570  libMesh::DenseSubVector<libMesh::Number> *Fw = NULL;
571 
572  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
573 
574  libMesh::DenseSubMatrix<libMesh::Number> &Kuu =
575  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.u()); // J_{uu}
576  libMesh::DenseSubMatrix<libMesh::Number> &Kuv =
577  context.get_elem_jacobian(this->_flow_vars.u(), this->_flow_vars.v()); // J_{uv}
578  libMesh::DenseSubMatrix<libMesh::Number> &Kvu =
579  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.u()); // J_{vu}
580  libMesh::DenseSubMatrix<libMesh::Number> &Kvv =
581  context.get_elem_jacobian(this->_flow_vars.v(), this->_flow_vars.v()); // J_{vv}
582 
583  libMesh::DenseSubMatrix<libMesh::Number> *Kuw = NULL;
584  libMesh::DenseSubMatrix<libMesh::Number> *Kvw = NULL;
585  libMesh::DenseSubMatrix<libMesh::Number> *Kwu = NULL;
586  libMesh::DenseSubMatrix<libMesh::Number> *Kwv = NULL;
587  libMesh::DenseSubMatrix<libMesh::Number> *Kww = NULL;
588 
589  libMesh::DenseSubMatrix<libMesh::Number> &Kpu =
590  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.u()); // J_{pu}
591  libMesh::DenseSubMatrix<libMesh::Number> &Kpv =
592  context.get_elem_jacobian(this->_press_var.p(), this->_flow_vars.v()); // J_{pv}
593  libMesh::DenseSubMatrix<libMesh::Number> *Kpw = NULL;
594 
595  if(this->_dim == 3)
596  {
597  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
598 
599  Kuw = &context.get_elem_jacobian
600  (this->_flow_vars.u(), this->_flow_vars.w()); // J_{uw}
601  Kvw = &context.get_elem_jacobian
602  (this->_flow_vars.v(), this->_flow_vars.w()); // J_{vw}
603  Kwu = &context.get_elem_jacobian
604  (this->_flow_vars.w(), this->_flow_vars.u()); // J_{wu}
605  Kwv = &context.get_elem_jacobian
606  (this->_flow_vars.w(), this->_flow_vars.v()); // J_{wv}
607  Kww = &context.get_elem_jacobian
608  (this->_flow_vars.w(), this->_flow_vars.w()); // J_{ww}
609 
610  Kpw = &context.get_elem_jacobian
611  (this->_press_var.p(), this->_flow_vars.w()); // J_{pw}
612  }
613 
614  unsigned int n_qpoints = context.get_element_qrule().n_points();
615 
616  for (unsigned int qp=0; qp != n_qpoints; qp++)
617  {
618  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
619 
620  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
621  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
622 
623  libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
624  context.fixed_interior_value( this->_flow_vars.v(), qp ) );
625  if( this->_dim == 3 )
626  {
627  U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
628  }
629 
630  /*
631  libMesh::Gradient grad_u, grad_v, grad_w;
632  grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
633  grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
634  if (_dim == 3)
635  grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
636  */
637 
638  libMesh::Real tau_M;
639 
640  libMesh::Real d_tau_M_d_rho;
641 
642  libMesh::Gradient d_tau_M_dU;
643 
644  libMesh::RealGradient RM_t;
645 
646  libMesh::Real d_RM_t_uvw_duvw;
647 
648  // Compute the viscosity at this qp
649  libMesh::Real _mu_qp = this->_mu(context, qp);
650 
651  if (compute_jacobian)
652  {
653  this->_stab_helper.compute_tau_momentum_and_derivs
654  ( context, qp, g, G, this->_rho, U, _mu_qp,
655  tau_M, d_tau_M_d_rho, d_tau_M_dU,
656  false );
657  this->_stab_helper.compute_res_momentum_transient_and_derivs
658  ( context, qp, this->_rho,
659  RM_t, d_RM_t_uvw_duvw );
660  }
661  else
662  {
663  tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, this->_rho, U, _mu_qp, false );
664  RM_t = this->_stab_helper.compute_res_momentum_transient( context, qp, this->_rho );
665  }
666 
667 
668  // Now a loop over the pressure degrees of freedom. This
669  // computes the contributions of the continuity equation.
670  for (unsigned int i=0; i != n_p_dofs; i++)
671  {
672  Fp(i) -= tau_M*RM_t*p_dphi[i][qp]*JxW[qp];
673 
674  if (compute_jacobian)
675  {
676  const libMesh::Real fixed_deriv =
677  context.get_fixed_solution_derivative();
678 
679  for (unsigned int j=0; j != n_u_dofs; j++)
680  {
681  Kpu(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t*p_dphi[i][qp]*fixed_deriv*JxW[qp];
682  Kpu(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*p_dphi[i][qp](0)*fixed_deriv*JxW[qp];
683 
684  Kpv(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t*p_dphi[i][qp]*fixed_deriv*JxW[qp];
685  Kpv(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*p_dphi[i][qp](1)*fixed_deriv*JxW[qp];
686  }
687 
688  if(this->_dim == 3)
689  {
690  for (unsigned int j=0; j != n_u_dofs; j++)
691  {
692  (*Kpw)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t*p_dphi[i][qp]*fixed_deriv*JxW[qp];
693  (*Kpw)(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*p_dphi[i][qp](2)*fixed_deriv*JxW[qp];
694  }
695  }
696  }
697  }
698 
699  for (unsigned int i=0; i != n_u_dofs; i++)
700  {
701  libMesh::Real test_func = this->_rho*U*u_gradphi[i][qp] +
702  _mu_qp*( u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2) );
703  libMesh::Gradient d_test_func_dU = this->_rho*u_gradphi[i][qp];
704 
705  //libMesh::RealGradient zeroth_order_term = - this->_rho*u_phi[i][qp]*(grad_u + grad_v + grad_w);
706 
707  Fu(i) -= tau_M*RM_t(0)*test_func*JxW[qp];
708 
709  Fv(i) -= tau_M*RM_t(1)*test_func*JxW[qp];
710 
711  if(this->_dim == 3)
712  {
713  (*Fw)(i) -= tau_M*RM_t(2)*test_func*JxW[qp];
714  }
715 
716  if (compute_jacobian)
717  {
718  const libMesh::Real fixed_deriv =
719  context.get_fixed_solution_derivative();
720 
721  for (unsigned int j=0; j != n_u_dofs; j++)
722  {
723  Kuu(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t(0)*test_func*fixed_deriv*JxW[qp];
724  Kuu(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*test_func*fixed_deriv*JxW[qp];
725  Kuu(i,j) -= tau_M*RM_t(0)*d_test_func_dU(0)*u_phi[j][qp]*fixed_deriv*JxW[qp];
726 
727  Kuv(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t(0)*test_func*fixed_deriv*JxW[qp];
728  Kuv(i,j) -= tau_M*RM_t(0)*d_test_func_dU(1)*u_phi[j][qp]*fixed_deriv*JxW[qp];
729 
730  Kvu(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t(1)*test_func*fixed_deriv*JxW[qp];
731  Kvu(i,j) -= tau_M*RM_t(1)*d_test_func_dU(0)*u_phi[j][qp]*fixed_deriv*JxW[qp];
732 
733  Kvv(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t(1)*test_func*fixed_deriv*JxW[qp];
734  Kvv(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*test_func*fixed_deriv*JxW[qp];
735  Kvv(i,j) -= tau_M*RM_t(1)*d_test_func_dU(1)*u_phi[j][qp]*fixed_deriv*JxW[qp];
736  }
737  if(this->_dim == 3)
738  {
739  for (unsigned int j=0; j != n_u_dofs; j++)
740  {
741  (*Kuw)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t(0)*test_func*fixed_deriv*JxW[qp];
742  (*Kuw)(i,j) -= tau_M*RM_t(0)*d_test_func_dU(2)*u_phi[j][qp]*fixed_deriv*JxW[qp];
743 
744  (*Kvw)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t(1)*test_func*fixed_deriv*JxW[qp];
745  (*Kvw)(i,j) -= tau_M*RM_t(1)*d_test_func_dU(2)*u_phi[j][qp]*fixed_deriv*JxW[qp];
746 
747  (*Kwu)(i,j) -= d_tau_M_dU(0)*u_phi[j][qp]*RM_t(2)*test_func*fixed_deriv*JxW[qp];
748  (*Kwu)(i,j) -= tau_M*RM_t(2)*d_test_func_dU(0)*u_phi[j][qp]*fixed_deriv*JxW[qp];
749 
750  (*Kwv)(i,j) -= d_tau_M_dU(1)*u_phi[j][qp]*RM_t(2)*test_func*fixed_deriv*JxW[qp];
751  (*Kwv)(i,j) -= tau_M*RM_t(2)*d_test_func_dU(1)*u_phi[j][qp]*fixed_deriv*JxW[qp];
752 
753  (*Kww)(i,j) -= d_tau_M_dU(2)*u_phi[j][qp]*RM_t(2)*test_func*fixed_deriv*JxW[qp];
754  (*Kww)(i,j) -= tau_M*d_RM_t_uvw_duvw*u_phi[j][qp]*test_func*fixed_deriv*JxW[qp];
755  (*Kww)(i,j) -= tau_M*RM_t(2)*d_test_func_dU(2)*u_phi[j][qp]*fixed_deriv*JxW[qp];
756  }
757  }
758  }
759  }
760 
761  }
762 
763 #ifdef GRINS_USE_GRVY_TIMERS
764  this->_timer->EndTimer("IncompressibleNavierStokesAdjointStabilization::mass_residual");
765 #endif
766  return;
767  }
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_momentum_and_derivs
void compute_tau_momentum_and_derivs(AssemblyContext &c, unsigned int qp, libMesh::RealGradient &g, libMesh::RealTensor &G, libMesh::Real rho, libMesh::Gradient U, libMesh::Real T, libMesh::Real &tau_M, libMesh::Real &d_tau_M_d_rho, libMesh::Gradient &d_tau_M_d_U, bool is_steady) const
Definition: inc_navier_stokes_stab_helper.h:231
GRINS::StabilizationHelper::compute_G
libMesh::RealTensor compute_G(libMesh::FEBase *fe, AssemblyContext &c, unsigned int qp) const
Definition: stab_helper.C:64
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_momentum_transient
libMesh::RealGradient compute_res_momentum_transient(AssemblyContext &context, unsigned int qp, const libMesh::Real rho) const
Definition: inc_navier_stokes_stab_helper.C:319
GRINS::IncompressibleNavierStokesBase::_rho
libMesh::Number _rho
Material parameters, read from input.
Definition: inc_navier_stokes_base.h:87
GRINS::IncompressibleNavierStokesBase::_dim
unsigned int _dim
Physical dimension of problem.
Definition: inc_navier_stokes_base.h:81
GRINS::IncompressibleNavierStokesStabilizationBase::_stab_helper
IncompressibleNavierStokesStabilizationHelper _stab_helper
Definition: inc_navier_stokes_stab_base.h:50
GRINS::StabilizationHelper::compute_g
libMesh::RealGradient compute_g(libMesh::FEBase *fe, AssemblyContext &c, unsigned int qp) const
Definition: stab_helper.C:47
GRINS::PressureFEVariable::p
VariableIndex p() const
Definition: pressure_fe_variable.h:71
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_res_momentum_transient_and_derivs
void compute_res_momentum_transient_and_derivs(AssemblyContext &context, unsigned int qp, const libMesh::Real rho, libMesh::RealGradient &res_M, libMesh::Real &d_res_Muvw_duvw) const
Definition: inc_navier_stokes_stab_helper.C:332
GRINS::IncompressibleNavierStokesStabilizationHelper::compute_tau_momentum
libMesh::Real compute_tau_momentum(AssemblyContext &c, unsigned int qp, libMesh::RealGradient &g, libMesh::RealTensor &G, libMesh::Real rho, libMesh::Gradient U, libMesh::Real mu, bool is_steady) const
Definition: inc_navier_stokes_stab_helper.h:217
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