GRINS-0.8.0
reacting_low_mach_navier_stokes_spgsm_stab.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // GRINS - General Reacting Incompressible Navier-Stokes
5 //
6 // Copyright (C) 2014-2017 Paul T. Bauman, Roy H. Stogner
7 // Copyright (C) 2010-2013 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 // This class
26 #include "grins/reacting_low_mach_navier_stokes_spgsm_stab.h"
27 
28 // GRINS
29 #include "grins/assembly_context.h"
30 
31 // libMesh
32 #include "libmesh/quadrature.h"
33 
34 namespace GRINS
35 {
36  template<typename Mixture, typename Evaluator>
37  ReactingLowMachNavierStokesSPGSMStabilization<Mixture,Evaluator>::ReactingLowMachNavierStokesSPGSMStabilization
38  ( const GRINS::PhysicsName& physics_name, const GetPot& input,libMesh::UniquePtr<Mixture> & gas_mix )
39  : ReactingLowMachNavierStokesStabilizationBase<Mixture,Evaluator>(physics_name,input,gas_mix)
40  {}
41 
42  template<typename Mixture, typename Evaluator>
43  void ReactingLowMachNavierStokesSPGSMStabilization<Mixture,Evaluator>::element_time_derivative
44  ( bool compute_jacobian,
45  AssemblyContext & context )
46  {
47  // The number of local degrees of freedom in each variable.
48  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
49  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
50  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
51  const VariableIndex s0_var = this->_species_vars.species(0);
52  const unsigned int n_s_dofs = context.get_dof_indices(s0_var).size();
53 
54  // Element Jacobian * quadrature weights for interior integration.
55  const std::vector<libMesh::Real> &JxW =
56  context.get_element_fe(this->_flow_vars.u())->get_JxW();
57 
58  // The pressure shape functions at interior quadrature points.
59  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
60  context.get_element_fe(this->_press_var.p())->get_dphi();
61 
62  const std::vector<std::vector<libMesh::Real> >& u_phi =
63  context.get_element_fe(this->_flow_vars.u())->get_phi();
64 
65  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
66  context.get_element_fe(this->_flow_vars.u())->get_dphi();
67 
68  // The temperature shape functions gradients at interior quadrature points.
69  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
70  context.get_element_fe(this->_temp_vars.T())->get_dphi();
71 
72  const std::vector<std::vector<libMesh::Gradient> >& s_gradphi = context.get_element_fe(s0_var)->get_dphi();
73 
74  // We're assuming the quadrature rule is the same for all variables
75  const std::vector<libMesh::Point>& u_qpoint =
76  context.get_element_fe(this->_flow_vars.u())->get_xyz();
77 
78  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
79 
80  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{u}
81 
82  libMesh::DenseSubVector<libMesh::Number>* Fv = NULL;
83  if( this->_flow_vars.dim() > 1 )
84  Fv = &context.get_elem_residual(this->_flow_vars.v()); // R_{v}
85 
86  libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
87  if( this->_flow_vars.dim() == 3 )
88  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
89 
90  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T()); // R_{T}
91 
92  unsigned int n_qpoints = context.get_element_qrule().n_points();
93 
94 
95  libMesh::FEBase* u_fe = context.get_element_fe(this->_flow_vars.u());
96  for (unsigned int qp=0; qp != n_qpoints; qp++)
97  {
98  libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
99 
100  libMesh::RealGradient U( context.interior_value( this->_flow_vars.u(), qp ) );
101  if( this->_flow_vars.dim() > 1 )
102  U(1) = context.interior_value( this->_flow_vars.v(), qp );
103  if( this->_flow_vars.dim() == 3 )
104  U(2) = context.interior_value( this->_flow_vars.w(), qp );
105 
106  std::vector<libMesh::Real> ws(this->n_species());
107  for(unsigned int s=0; s < this->_n_species; s++ )
108  {
109  ws[s] = context.fixed_interior_value(this->_species_vars.species(s), qp);
110  }
111 
112  Evaluator gas_evaluator( *(this->_gas_mixture) );
113  const libMesh::Real R_mix = gas_evaluator.R_mix(ws);
114  const libMesh::Real p0 = this->get_p0_steady(context,qp);
115  libMesh::Real rho = this->rho(T, p0, R_mix);
116 
117  const libMesh::Real cp = gas_evaluator.cp(T,p0,ws);
118 
119  std::vector<libMesh::Real> D( this->n_species() );
120  libMesh::Real mu, k;
121 
122  gas_evaluator.mu_and_k_and_D( T, rho, cp, ws, mu, k, D );
123 
124  libMesh::RealGradient g = this->_stab_helper.compute_g( u_fe, context, qp );
125  libMesh::RealTensor G = this->_stab_helper.compute_G( u_fe, context, qp );
126 
127  // Taus
128  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
129 
130  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
131 
132  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
133 
134 
135 
136  // Strong form residuals
137  libMesh::RealGradient RM_s = 0.0;
138  libMesh::Real RC_s = 0.0;
139  libMesh::Real RE_s = 0.0;
140  std::vector<libMesh::Real> Rs_s;
141 
142  this->compute_res_steady( context, qp, RC_s, RM_s, RE_s, Rs_s );
143 
144  const libMesh::Number r = u_qpoint[qp](0);
145 
146  libMesh::Real jac = JxW[qp];
147 
148  if( this->_is_axisymmetric )
149  jac *= r;
150 
151  // Pressure PSPG term
152  for (unsigned int i=0; i != n_p_dofs; i++)
153  Fp(i) += tau_M*RM_s*p_dphi[i][qp]*jac;
154 
155  // Momentum SUPG + div-div terms
156  for (unsigned int i=0; i != n_u_dofs; i++)
157  {
158  Fu(i) += ( - tau_C*RC_s*u_gradphi[i][qp](0)
159  - tau_M*RM_s(0)*rho*U*u_gradphi[i][qp] )*jac;
160 
161  if( this->_is_axisymmetric )
162  Fu(i) += (-tau_C*RC_s/r)*u_phi[i][qp]*jac;
163 
164  if( this->_flow_vars.dim() > 1 )
165  (*Fv)(i) += ( - tau_C*RC_s*u_gradphi[i][qp](1)
166  - tau_M*RM_s(1)*rho*U*u_gradphi[i][qp] )*jac;
167 
168  if( this->_flow_vars.dim() == 3 )
169  (*Fw)(i) += ( - tau_C*RC_s*u_gradphi[i][qp](2)
170  - tau_M*RM_s(2)*rho*U*u_gradphi[i][qp] )*jac;
171  }
172 
173  // Energy SUPG terms
174  for (unsigned int i=0; i != n_T_dofs; i++)
175  FT(i) -= rho*cp*tau_E*RE_s*U*T_gradphi[i][qp]*jac;
176 
177  // Species SUPG terms
178  for(unsigned int s=0; s < this->n_species(); s++)
179  {
180  libMesh::DenseSubVector<libMesh::Number> &Fs =
181  context.get_elem_residual(this->_species_vars.species(s));
182 
183  libMesh::Real tau_s = this->_stab_helper.compute_tau_species( context, qp, g, G, rho, U, D[s], this->_is_steady );
184 
185  for (unsigned int i=0; i != n_s_dofs; i++)
186  Fs(i) -= rho*tau_s*Rs_s[s]*U*s_gradphi[i][qp]*jac;
187  }
188 
189  if(compute_jacobian)
190  libmesh_not_implemented();
191  }
192  }
193 
194  template<typename Mixture, typename Evaluator>
195  void ReactingLowMachNavierStokesSPGSMStabilization<Mixture,Evaluator>::mass_residual
196  ( bool compute_jacobian, AssemblyContext & context )
197  {
198  const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
199  const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
200  const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
201  const VariableIndex s0_var = this->_species_vars.species(0);
202  const unsigned int n_s_dofs = context.get_dof_indices(s0_var).size();
203 
204  const std::vector<libMesh::Real> &JxW =
205  context.get_element_fe(this->_flow_vars.u())->get_JxW();
206 
207  const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
208  context.get_element_fe(this->_press_var.p())->get_dphi();
209 
210  const std::vector<std::vector<libMesh::Real> >& u_phi =
211  context.get_element_fe(this->_flow_vars.u())->get_phi();
212 
213  const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
214  context.get_element_fe(this->_flow_vars.u())->get_dphi();
215 
216  const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
217  context.get_element_fe(this->_temp_vars.T())->get_dphi();
218 
219  const std::vector<std::vector<libMesh::Gradient> >& s_gradphi = context.get_element_fe(s0_var)->get_dphi();
220 
221  libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p()); // R_{p}
222  libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u()); // R_{u}
223 
224  libMesh::DenseSubVector<libMesh::Number>* Fv = NULL;
225  if( this->_flow_vars.dim() > 1 )
226  Fv = &context.get_elem_residual(this->_flow_vars.v()); // R_{v}
227 
228  libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
229  if( this->_flow_vars.dim() == 3 )
230  Fw = &context.get_elem_residual(this->_flow_vars.w()); // R_{w}
231 
232  libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
233 
234  libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
235 
236  unsigned int n_qpoints = context.get_element_qrule().n_points();
237 
238  const std::vector<libMesh::Point>& u_qpoint =
239  context.get_element_fe(this->_flow_vars.u())->get_xyz();
240 
241  for (unsigned int qp=0; qp != n_qpoints; qp++)
242  {
243  libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
244  libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
245 
246  libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
247 
248  libMesh::Number u = context.fixed_interior_value(this->_flow_vars.u(), qp);
249 
250  libMesh::NumberVectorValue U(u);
251  if (this->_flow_vars.dim() > 1)
252  U(1) = context.fixed_interior_value(this->_flow_vars.v(), qp);
253  if (this->_flow_vars.dim() == 3)
254  U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
255 
256  std::vector<libMesh::Real> ws(this->n_species());
257  for(unsigned int s=0; s < this->_n_species; s++ )
258  ws[s] = context.fixed_interior_value(this->_species_vars.species(s), qp);
259 
260  Evaluator gas_evaluator( *(this->_gas_mixture) );
261  const libMesh::Real R_mix = gas_evaluator.R_mix(ws);
262  const libMesh::Real p0 = this->get_p0_steady(context,qp);
263  libMesh::Real rho = this->rho(T, p0, R_mix);
264 
265  const libMesh::Real cp = gas_evaluator.cp(T,p0,ws);
266 
267  std::vector<libMesh::Real> D( this->n_species() );
268  libMesh::Real mu, k;
269 
270  gas_evaluator.mu_and_k_and_D( T, rho, cp, ws, mu, k, D );
271 
272  libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, false );
273  libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
274  libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, false );
275 
276  // Strong residuals
277  libMesh::Real RC_t;
278  libMesh::RealGradient RM_t;
279  libMesh::Real RE_t;
280  std::vector<libMesh::Real> Rs_t(this->n_species());
281 
282  this->compute_res_transient( context, qp, RC_t, RM_t, RE_t, Rs_t );
283 
284  libMesh::Real jac = JxW[qp];
285  const libMesh::Number r = u_qpoint[qp](0);
286 
287  if( this->_is_axisymmetric )
288  jac *= r;
289 
290  for (unsigned int i=0; i != n_p_dofs; i++)
291  Fp(i) -= tau_M*RM_t*p_dphi[i][qp]*jac;
292 
293  for(unsigned int s=0; s < this->n_species(); s++)
294  {
295  libMesh::DenseSubVector<libMesh::Number> &Fs =
296  context.get_elem_residual(this->_species_vars.species(s));
297 
298  libMesh::Real tau_s = this->_stab_helper.compute_tau_species( context, qp, g, G, rho, U, D[s], false );
299 
300  for (unsigned int i=0; i != n_s_dofs; i++)
301  Fs(i) -= rho*tau_s*Rs_t[s]*U*s_gradphi[i][qp]*jac;
302  }
303 
304  for (unsigned int i=0; i != n_u_dofs; i++)
305  {
306  Fu(i) -= ( tau_C*RC_t*u_gradphi[i][qp](0)
307  + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp] )*jac;
308 
309  if( this->_is_axisymmetric )
310  Fu(i) += (-tau_C*RC_t/r)*u_phi[i][qp]*jac;
311 
312  if( this->_flow_vars.dim() > 1 )
313  (*Fv)(i) -= ( tau_C*RC_t*u_gradphi[i][qp](1)
314  + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp] )*jac;
315 
316  if( this->_flow_vars.dim() == 3 )
317  (*Fw)(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
318  + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp] )*jac;
319  }
320 
321  for (unsigned int i=0; i != n_T_dofs; i++)
322  FT(i) -= rho*cp*tau_E*RE_t*U*T_gradphi[i][qp]*jac;
323 
324  if(compute_jacobian)
325  libmesh_not_implemented();
326  }
327 
328  return;
329  }
330 } // end namespace GRINS
GRINS::ReactingLowMachNavierStokesSPGSMStabilization::mass_residual
virtual void mass_residual(bool compute_jacobian, AssemblyContext &context)
Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part...
Definition: reacting_low_mach_navier_stokes_spgsm_stab.C:196
GRINS::VariableIndex
unsigned int VariableIndex
More descriptive name of the type used for variable indices.
Definition: var_typedefs.h:42
assembly_context.h
GRINS
GRINS namespace.
Definition: arrhenius_catalycity.h:31
GRINS::ReactingLowMachNavierStokesSPGSMStabilization::element_time_derivative
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context)
Time dependent part(s) of physics for element interiors.
Definition: reacting_low_mach_navier_stokes_spgsm_stab.C:44
GRINS::ReactingLowMachNavierStokesStabilizationBase
Adds VMS-based stabilization to LowMachNavierStokes physics class.
Definition: reacting_low_mach_navier_stokes_stab_base.h:37
GRINS::PhysicsName
std::string PhysicsName
Definition: physics_naming.h:33
reacting_low_mach_navier_stokes_spgsm_stab.h
Generated on Tue Dec 19 2017 12:47:28 for GRINS-0.8.0 by  doxygen 1.8.9.1