36 #include "libmesh/quadrature.h"
37 #include "libmesh/fem_system.h"
41 template<
typename Mixture,
typename Evaluator>
44 _p_pinning(input,physics_name),
55 template<
typename Mixture,
typename Evaluator>
59 _p_pinning.check_pin_location(system.get_mesh());
62 template<
typename Mixture,
typename Evaluator>
68 if( input.have_variable(section) )
70 unsigned int n_vars = input.vector_variable_size(section);
72 for(
unsigned int v = 0; v < n_vars; v++ )
74 std::string name = input(section,
"DIE!",v);
76 if( name == std::string(
"rho") )
80 else if( name == std::string(
"mu") )
84 else if( name == std::string(
"k") )
88 else if( name == std::string(
"cp") )
92 else if( name == std::string(
"mole_fractions") )
94 this->_mole_fractions_index.resize(this->n_species());
96 for(
unsigned int s = 0; s < this->n_species(); s++ )
98 this->_mole_fractions_index[s] = postprocessing.
register_quantity(
"X_"+this->_gas_mixture.species_name(s) );
101 else if( name == std::string(
"h_s") )
103 this->_h_s_index.resize(this->n_species());
105 for(
unsigned int s = 0; s < this->n_species(); s++ )
107 this->_h_s_index[s] = postprocessing.
register_quantity(
"h_"+this->_gas_mixture.species_name(s) );
110 else if( name == std::string(
"omega_dot") )
112 this->_omega_dot_index.resize(this->n_species());
114 for(
unsigned int s = 0; s < this->n_species(); s++ )
116 this->_omega_dot_index[s] = postprocessing.
register_quantity(
"omega_dot_"+this->_gas_mixture.species_name(s) );
119 std::cout <<
"omega_dot size = " << _omega_dot_index.size() << std::endl;
124 <<
" Found " << name << std::endl
125 <<
" Acceptable values are: rho" << std::endl
126 <<
" mu" << std::endl
128 <<
" cp" << std::endl
129 <<
" mole_fractions" << std::endl
130 <<
" omega_dot" << std::endl;
139 template<
typename Mixture,
typename Evaluator>
146 context.get_side_fe(this->_flow_vars.u())->get_JxW();
147 context.get_side_fe(this->_flow_vars.u())->get_phi();
148 context.get_side_fe(this->_flow_vars.u())->get_dphi();
149 context.get_side_fe(this->_flow_vars.u())->get_xyz();
151 context.get_side_fe(this->_temp_vars.T())->get_JxW();
152 context.get_side_fe(this->_temp_vars.T())->get_phi();
153 context.get_side_fe(this->_temp_vars.T())->get_dphi();
154 context.get_side_fe(this->_temp_vars.T())->get_xyz();
157 template<
typename Mixture,
typename Evaluator>
162 if( compute_jacobian )
164 libmesh_not_implemented();
170 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
171 const unsigned int n_s_dofs = context.get_dof_indices(s0_var).size();
172 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
173 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
176 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
178 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
181 const std::vector<libMesh::Real>& JxW =
182 context.get_element_fe(this->_flow_vars.u())->get_JxW();
185 const std::vector<std::vector<libMesh::Real> >& p_phi =
186 context.get_element_fe(this->_press_var.p())->get_phi();
189 const std::vector<std::vector<libMesh::Real> >& s_phi = context.get_element_fe(s0_var)->get_phi();
192 const std::vector<std::vector<libMesh::Gradient> >& s_grad_phi = context.get_element_fe(s0_var)->get_dphi();
195 const std::vector<std::vector<libMesh::Real> >& u_phi =
196 context.get_element_fe(this->_flow_vars.u())->get_phi();
199 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
200 context.get_element_fe(this->_flow_vars.u())->get_dphi();
203 const std::vector<std::vector<libMesh::Real> >& T_phi =
204 context.get_element_fe(this->_temp_vars.T())->get_phi();
207 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
208 context.get_element_fe(this->_temp_vars.T())->get_dphi();
210 const std::vector<libMesh::Point>& u_qpoint =
211 context.get_element_fe(this->_flow_vars.u())->get_xyz();
213 libMesh::DenseSubVector<libMesh::Number>& Fp = context.get_elem_residual(this->_press_var.p());
215 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
216 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
217 libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
219 if( this->_dim == 3 )
221 Fw = &context.get_elem_residual(this->_flow_vars.w());
224 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
226 unsigned int n_qpoints = context.get_element_qrule().n_points();
227 for (
unsigned int qp=0; qp != n_qpoints; qp++)
231 libMesh::Number u, v, T;
237 const libMesh::Gradient& grad_T =
240 libMesh::NumberVectorValue U(u,v);
247 libMesh::Gradient grad_w;
251 libMesh::Number divU = grad_u(0) + grad_v(1);
255 libMesh::NumberVectorValue grad_uT( grad_u(0), grad_v(0) );
256 libMesh::NumberVectorValue grad_vT( grad_u(1), grad_v(1) );
257 libMesh::NumberVectorValue grad_wT;
258 if( this->_dim == 3 )
260 grad_uT(2) = grad_w(0);
261 grad_vT(2) = grad_w(1);
262 grad_wT = libMesh::NumberVectorValue( grad_u(2), grad_v(2), grad_w(2) );
268 const std::vector<libMesh::Real>& D =
277 const std::vector<libMesh::Real>& omega_dot =
280 const std::vector<libMesh::Real>& h =
283 const libMesh::Number r = u_qpoint[qp](0);
285 libMesh::Real jac = JxW[qp];
296 libmesh_assert_equal_to( grad_ws.size(), this->_n_species );
299 libMesh::Gradient mass_term(0.0,0.0,0.0);
300 for(
unsigned int s=0; s < this->_n_species; s++ )
302 mass_term += grad_ws[s]/this->_gas_mixture.M(s);
306 for (
unsigned int i=0; i != n_p_dofs; i++)
308 Fp(i) += (-U*(mass_term + grad_T/T) + divU)*jac*p_phi[i][qp];
312 for(
unsigned int s=0; s < this->_n_species; s++ )
314 libMesh::DenseSubVector<libMesh::Number> &Fs =
315 context.get_elem_residual(this->_species_vars.species(s));
317 const libMesh::Real term1 = -rho*(U*grad_ws[s]) + omega_dot[s];
318 const libMesh::Gradient term2 = -rho*D[s]*grad_ws[s];
320 for (
unsigned int i=0; i != n_s_dofs; i++)
323 Fs(i) += ( term1*s_phi[i][qp] + term2*s_grad_phi[i][qp] )*jac;
328 for (
unsigned int i=0; i != n_u_dofs; i++)
330 Fu(i) += ( -rho*U*grad_u*u_phi[i][qp]
331 + p*u_gradphi[i][qp](0)
332 - mu*(u_gradphi[i][qp]*grad_u + u_gradphi[i][qp]*grad_uT
333 - 2.0/3.0*divU*u_gradphi[i][qp](0) )
334 + rho*this->_g(0)*u_phi[i][qp]
341 Fu(i) += u_phi[i][qp]*( p/r - 2*mu*U(0)/(r*r) - 2.0/3.0*mu*divU/r )*jac;
344 Fv(i) += ( -rho*U*grad_v*u_phi[i][qp]
345 + p*u_gradphi[i][qp](1)
346 - mu*(u_gradphi[i][qp]*grad_v + u_gradphi[i][qp]*grad_vT
347 - 2.0/3.0*divU*u_gradphi[i][qp](1) )
348 + rho*this->_g(1)*u_phi[i][qp]
353 (*Fw)(i) += ( -rho*U*grad_w*u_phi[i][qp]
354 + p*u_gradphi[i][qp](2)
355 - mu*(u_gradphi[i][qp]*grad_w + u_gradphi[i][qp]*grad_wT
356 - 2.0/3.0*divU*u_gradphi[i][qp](2) )
357 + rho*this->_g(2)*u_phi[i][qp]
363 libMesh::Real chem_term = 0.0;
364 for(
unsigned int s=0; s < this->_n_species; s++ )
366 chem_term += h[s]*omega_dot[s];
369 for (
unsigned int i=0; i != n_T_dofs; i++)
371 FT(i) += ( ( -rho*cp*U*grad_T - chem_term )*T_phi[i][qp]
372 - k*grad_T*T_gradphi[i][qp] )*jac;
380 template<
typename Mixture,
typename Evaluator>
386 if( this->_pin_pressure )
388 _p_pinning.pin_value( context, compute_jacobian, this->_press_var.p() );
394 template<
typename Mixture,
typename Evaluator>
399 libmesh_not_implemented();
404 template<
typename Mixture,
typename Evaluator>
408 Evaluator gas_evaluator( this->_gas_mixture );
410 const unsigned int n_qpoints = context.get_element_qrule().n_points();
412 std::vector<libMesh::Real> u, v, w, T, p, p0;
420 p0.resize(n_qpoints);
422 std::vector<libMesh::Gradient> grad_u, grad_v, grad_w, grad_T;
423 grad_u.resize(n_qpoints);
424 grad_v.resize(n_qpoints);
426 grad_w.resize(n_qpoints);
428 grad_T.resize(n_qpoints);
430 std::vector<std::vector<libMesh::Real> > mass_fractions;
431 std::vector<std::vector<libMesh::Gradient> > grad_mass_fractions;
432 mass_fractions.resize(n_qpoints);
433 grad_mass_fractions.resize(n_qpoints);
435 std::vector<libMesh::Real> M;
438 std::vector<libMesh::Real> R;
441 std::vector<libMesh::Real> rho;
442 rho.resize(n_qpoints);
444 std::vector<libMesh::Real> cp;
445 cp.resize(n_qpoints);
447 std::vector<libMesh::Real> mu;
448 mu.resize(n_qpoints);
450 std::vector<libMesh::Real> k;
453 std::vector<std::vector<libMesh::Real> > h_s;
454 h_s.resize(n_qpoints);
456 std::vector<std::vector<libMesh::Real> > D_s;
457 D_s.resize(n_qpoints);
459 std::vector<std::vector<libMesh::Real> > omega_dot_s;
460 omega_dot_s.resize(n_qpoints);
462 for (
unsigned int qp = 0; qp != n_qpoints; ++qp)
464 u[qp] = context.interior_value(this->_flow_vars.u(), qp);
465 v[qp] = context.interior_value(this->_flow_vars.v(), qp);
467 grad_u[qp] = context.interior_gradient(this->_flow_vars.u(), qp);
468 grad_v[qp] = context.interior_gradient(this->_flow_vars.v(), qp);
471 w[qp] = context.interior_value(this->_flow_vars.w(), qp);
472 grad_w[qp] = context.interior_gradient(this->_flow_vars.w(), qp);
475 T[qp] = context.interior_value(this->_temp_vars.T(), qp);
476 grad_T[qp] = context.interior_gradient(this->_temp_vars.T(), qp);
478 p[qp] = context.interior_value(this->_press_var.p(), qp);
479 p0[qp] = this->get_p0_steady(context, qp);
481 mass_fractions[qp].resize(this->_n_species);
482 grad_mass_fractions[qp].resize(this->_n_species);
483 h_s[qp].resize(this->_n_species);
485 for(
unsigned int s = 0; s < this->_n_species; s++ )
489 mass_fractions[qp][s] = std::max( context.interior_value(this->_species_vars.species(s),qp), 0.0 );
490 grad_mass_fractions[qp][s] = context.interior_gradient(this->_species_vars.species(s),qp);
491 h_s[qp][s] = gas_evaluator.h_s( T[qp], s );
494 M[qp] = gas_evaluator.M_mix( mass_fractions[qp] );
496 R[qp] = gas_evaluator.R_mix( mass_fractions[qp] );
498 rho[qp] = this->rho( T[qp], p0[qp], R[qp] );
500 cp[qp] = gas_evaluator.cp(T[qp], p0[qp], mass_fractions[qp]);
502 D_s[qp].resize(this->_n_species);
504 gas_evaluator.mu_and_k_and_D( T[qp], rho[qp], cp[qp], mass_fractions[qp],
505 mu[qp], k[qp], D_s[qp] );
507 omega_dot_s[qp].resize(this->_n_species);
508 gas_evaluator.omega_dot( T[qp], rho[qp], mass_fractions[qp], omega_dot_s[qp] );
539 template<
typename Mixture,
typename Evaluator>
542 const libMesh::Point& point,
543 libMesh::Real& value )
545 Evaluator gas_evaluator( this->_gas_mixture );
547 if( quantity_index == this->_rho_index )
549 std::vector<libMesh::Real> Y( this->_n_species );
550 libMesh::Real T = this->T(point,context);
551 libMesh::Real p0 = this->get_p0_steady(context,point);
552 this->mass_fractions( point, context, Y );
554 value = this->rho( T, p0, gas_evaluator.R_mix(Y) );
556 else if( quantity_index == this->_mu_index )
558 std::vector<libMesh::Real> Y( this->_n_species );
559 libMesh::Real T = this->T(point,context);
560 this->mass_fractions( point, context, Y );
561 libMesh::Real p0 = this->get_p0_steady(context,point);
563 value = gas_evaluator.mu( T, p0, Y );
565 else if( quantity_index == this->_k_index )
567 std::vector<libMesh::Real> Y( this->_n_species );
568 libMesh::Real T = this->T(point,context);
569 this->mass_fractions( point, context, Y );
570 libMesh::Real p0 = this->get_p0_steady(context,point);
572 value = gas_evaluator.k( T, p0, Y );
574 else if( quantity_index == this->_cp_index )
576 std::vector<libMesh::Real> Y( this->_n_species );
577 libMesh::Real T = this->T(point,context);
578 this->mass_fractions( point, context, Y );
579 libMesh::Real p0 = this->get_p0_steady(context,point);
581 value = gas_evaluator.cp( T, p0, Y );
586 if( !this->_h_s_index.empty() )
588 libmesh_assert_equal_to( _h_s_index.size(), this->n_species() );
590 for(
unsigned int s = 0; s < this->n_species(); s++ )
592 if( quantity_index == this->_h_s_index[s] )
594 libMesh::Real T = this->T(point,context);
596 value = gas_evaluator.h_s( T, s );
602 if( !this->_mole_fractions_index.empty() )
604 libmesh_assert_equal_to( _mole_fractions_index.size(), this->n_species() );
606 for(
unsigned int s = 0; s < this->n_species(); s++ )
608 if( quantity_index == this->_mole_fractions_index[s] )
610 std::vector<libMesh::Real> Y( this->_n_species );
611 this->mass_fractions( point, context, Y );
613 libMesh::Real M = gas_evaluator.M_mix(Y);
615 value = gas_evaluator.X( s, M, Y[s] );
621 if( !this->_omega_dot_index.empty() )
623 libmesh_assert_equal_to( _omega_dot_index.size(), this->n_species() );
625 for(
unsigned int s = 0; s < this->n_species(); s++ )
627 if( quantity_index == this->_omega_dot_index[s] )
629 std::vector<libMesh::Real> Y( this->n_species() );
630 this->mass_fractions( point, context, Y );
632 libMesh::Real T = this->T(point,context);
634 libMesh::Real p0 = this->get_p0_steady(context,point);
636 libMesh::Real rho = this->rho( T, p0, gas_evaluator.R_mix(Y) );
638 std::vector<libMesh::Real> omega_dot( this->n_species() );
639 gas_evaluator.omega_dot( T, rho, Y, omega_dot );
641 value = omega_dot[s];
static bool is_axisymmetric()
unsigned int VariableIndex
More descriptive name of the type used for variable indices.
unsigned int register_quantity(std::string name)
Register quantity to be postprocessed.
GRINS::ICHandlingBase * _ic_handler
virtual void register_postprocessing_vars(const GetPot &input, PostProcessedQuantities< libMesh::Real > &postprocessing)
Register postprocessing variables for ReactingLowMachNavierStokes.
void set_values(unsigned int quantity, std::vector< libMesh::Number > &values)
static PhysicsName reacting_low_mach_navier_stokes()
void set_vector_gradient_values(unsigned int quantity, std::vector< std::vector< libMesh::Gradient > > &values)
Base class for reading and handling initial conditions for physics classes.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
const std::vector< std::vector< libMesh::Number > > & get_cached_vector_values(unsigned int quantity) const
virtual void compute_element_time_derivative_cache(const AssemblyContext &context, CachedValues &cache)
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
bool _pin_pressure
Enable pressure pinning.
const std::vector< std::vector< libMesh::Gradient > > & get_cached_vector_gradient_values(unsigned int quantity) const
Interface with libMesh for solving Multiphysics problems.
const std::vector< libMesh::Gradient > & get_cached_gradient_values(unsigned int quantity) const
ReactingLowMachNavierStokes()
void set_vector_values(unsigned int quantity, std::vector< std::vector< libMesh::Number > > &values)
virtual void element_constraint(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Constraint part(s) of physics for element interiors.
void set_gradient_values(unsigned int quantity, std::vector< libMesh::Gradient > &values)
virtual void mass_residual(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part...
virtual void auxiliary_init(MultiphysicsSystem &system)
Any auxillary initialization a Physics class may need.
virtual void compute_postprocessed_quantity(unsigned int quantity_index, const AssemblyContext &context, const libMesh::Point &point, libMesh::Real &value)
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
const std::vector< libMesh::Number > & get_cached_values(unsigned int quantity) const