[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
euler.cpp
1 #include <cmath>
2 #include <vector>
3 
4 #include "ADTypes.hpp"
5 
6 #include "physics.h"
7 #include "euler.h"
8 
9 namespace PHiLiP {
10 namespace Physics {
11 
12 template <int dim, int nstate, typename real>
13 Euler<dim,nstate,real>::Euler (
14  const Parameters::AllParameters *const parameters_input,
15  const double ref_length,
16  const double gamma_gas,
17  const double mach_inf,
18  const double angle_of_attack,
19  const double side_slip_angle,
20  std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function,
21  const two_point_num_flux_enum two_point_num_flux_type_input,
22  const bool has_nonzero_diffusion,
23  const bool has_nonzero_physical_source)
24  : PhysicsBase<dim,nstate,real>(parameters_input, has_nonzero_diffusion,has_nonzero_physical_source,manufactured_solution_function)
25  , ref_length(ref_length)
26  , gam(gamma_gas)
27  , gamm1(gam-1.0)
28  , density_inf(1.0) // Nondimensional - Free stream values
29  , mach_inf(mach_inf)
30  , mach_inf_sqr(mach_inf*mach_inf)
31  , angle_of_attack(angle_of_attack)
32  , side_slip_angle(side_slip_angle)
33  , sound_inf(1.0/(mach_inf))
34  , pressure_inf(1.0/(gam*mach_inf_sqr))
35  , entropy_inf(pressure_inf*pow(density_inf,-gam))
36  , two_point_num_flux_type(two_point_num_flux_type_input)
37  //, internal_energy_inf(1.0/(gam*(gam-1.0)*mach_inf_sqr))
38  // Note: Eq.(3.11.18) has a typo in internal_energy_inf expression, mach_inf_sqr should be in denominator.
39 {
40  static_assert(nstate==dim+2, "Physics::Euler() should be created with nstate=dim+2");
41 
42  // Nondimensional temperature at infinity
43  temperature_inf = gam*pressure_inf/density_inf * mach_inf_sqr; // Note by JB: this can simply be set = 1
44 
45  // For now, don't allow side-slip angle
46  if (std::abs(side_slip_angle) >= 1e-14) {
47  this->pcout << "Side slip angle = " << side_slip_angle << ". Side_slip_angle must be zero. " << std::endl;
48  this->pcout << "I have not figured out the side slip angles just yet." << std::endl;
49  std::abort();
50  }
51  if(dim==1) {
52  velocities_inf[0] = 1.0;
53  } else if(dim==2) {
54  velocities_inf[0] = cos(angle_of_attack);
55  velocities_inf[1] = sin(angle_of_attack); // Maybe minus?? -- Clarify with Doug
56  } else if (dim==3) {
57  velocities_inf[0] = cos(angle_of_attack)*cos(side_slip_angle);
58  velocities_inf[1] = sin(angle_of_attack)*cos(side_slip_angle);
59  velocities_inf[2] = sin(side_slip_angle);
60  }
61 
62  assert(std::abs(velocities_inf.norm() - 1.0) < 1e-14);
63 
64  double velocity_inf_sqr = 1.0;
65  dynamic_pressure_inf = 0.5 * density_inf * velocity_inf_sqr;
66 }
67 
68 template <int dim, int nstate, typename real>
69 std::array<real,nstate> Euler<dim,nstate,real>
70 ::source_term (
71  const dealii::Point<dim,real> &pos,
72  const std::array<real,nstate> &conservative_soln,
73  const real current_time,
74  const dealii::types::global_dof_index /*cell_index*/) const
75 {
76  return source_term(pos,conservative_soln,current_time);
77 }
78 
79 template <int dim, int nstate, typename real>
80 std::array<real,nstate> Euler<dim,nstate,real>
81 ::source_term (
82  const dealii::Point<dim,real> &pos,
83  const std::array<real,nstate> &/*conservative_soln*/,
84  const real /*current_time*/) const
85 {
86  std::array<real,nstate> source_term = convective_source_term(pos);
87  return source_term;
88 }
89 
90 template <int dim, int nstate, typename real>
91 std::array<real,nstate> Euler<dim,nstate,real>
92 ::get_manufactured_solution_value (
93  const dealii::Point<dim,real> &pos) const
94 {
95  std::array<real,nstate> manufactured_solution;
96  for (int s=0; s<nstate; s++) {
97  manufactured_solution[s] = this->manufactured_solution_function->value (pos, s);
98  if (s==0) {
99  assert(manufactured_solution[s] > 0);
100  }
101  }
102  return manufactured_solution;
103 }
104 
105 template <int dim, int nstate, typename real>
106 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim,nstate,real>
107 ::get_manufactured_solution_gradient (
108  const dealii::Point<dim,real> &pos) const
109 {
110  std::vector<dealii::Tensor<1,dim,real>> manufactured_solution_gradient_dealii(nstate);
111  this->manufactured_solution_function->vector_gradient(pos,manufactured_solution_gradient_dealii);
112  std::array<dealii::Tensor<1,dim,real>,nstate> manufactured_solution_gradient;
113  for (int d=0;d<dim;d++) {
114  for (int s=0; s<nstate; s++) {
115  manufactured_solution_gradient[s][d] = manufactured_solution_gradient_dealii[s][d];
116  }
117  }
118  return manufactured_solution_gradient;
119 }
120 
121 template <int dim, int nstate, typename real>
122 std::array<real,nstate> Euler<dim,nstate,real>
123 ::convective_source_term (
124  const dealii::Point<dim,real> &pos) const
125 {
126  const std::array<real,nstate> manufactured_solution = get_manufactured_solution_value(pos);
127  const std::array<dealii::Tensor<1,dim,real>,nstate> manufactured_solution_gradient = get_manufactured_solution_gradient(pos);
128 
129  dealii::Tensor<1,nstate,real> convective_flux_divergence;
130  for (int d=0;d<dim;d++) {
131  dealii::Tensor<1,dim,real> normal;
132  normal[d] = 1.0;
133  const dealii::Tensor<2,nstate,real> jacobian = convective_flux_directional_jacobian(manufactured_solution, normal);
134 
135  //convective_flux_divergence += jacobian*manufactured_solution_gradient[d];
136  for (int sr = 0; sr < nstate; ++sr) {
137  real jac_grad_row = 0.0;
138  for (int sc = 0; sc < nstate; ++sc) {
139  jac_grad_row += jacobian[sr][sc]*manufactured_solution_gradient[sc][d];
140  }
141  convective_flux_divergence[sr] += jac_grad_row;
142  }
143  }
144  std::array<real,nstate> convective_source_term;
145  for (int s=0; s<nstate; s++) {
146  convective_source_term[s] = convective_flux_divergence[s];
147  }
148 
149  return convective_source_term;
150 }
151 
152 template <int dim, int nstate, typename real>
153 template<typename real2>
154 bool Euler<dim,nstate,real>::check_positive_quantity(real2 &qty, const std::string qty_name) const {
155  using limiter_enum = Parameters::LimiterParam::LimiterType;
156  bool qty_is_positive;
157 
158  if (this->all_parameters->limiter_param.bound_preserving_limiter != limiter_enum::positivity_preservingZhang2010
159  && this->all_parameters->limiter_param.bound_preserving_limiter != limiter_enum::positivity_preservingWang2012) {
160  if (qty < 0.0) {
161  // Refer to base class for non-physical results handling
162  qty = this->template handle_non_physical_result<real2>(qty_name + " is negative.");
163  qty_is_positive = false;
164  }
165  else {
166  qty_is_positive = true;
167  }
168  } else {
169  qty_is_positive = true;
170  }
171 
172  return qty_is_positive;
173 }
174 
175 template <int dim, int nstate, typename real>
176 template<typename real2>
177 inline std::array<real2,nstate> Euler<dim,nstate,real>
178 ::convert_conservative_to_primitive ( const std::array<real2,nstate> &conservative_soln ) const
179 {
180  std::array<real2, nstate> primitive_soln;
181 
182  real2 density = conservative_soln[0];
183  dealii::Tensor<1,dim,real2> vel = compute_velocities<real2>(conservative_soln);
184  real2 pressure = compute_pressure<real2>(conservative_soln);
185 
186  check_positive_quantity<real2>(density, "density");
187  check_positive_quantity<real2>(pressure, "pressure");
188  primitive_soln[0] = density;
189  for (int d=0; d<dim; ++d) {
190  primitive_soln[1+d] = vel[d];
191  }
192  primitive_soln[nstate-1] = pressure;
193 
194  return primitive_soln;
195 }
196 
197 template <int dim, int nstate, typename real>
198 inline std::array<real,nstate> Euler<dim,nstate,real>
199 ::convert_primitive_to_conservative ( const std::array<real,nstate> &primitive_soln ) const
200 {
201 
202  const real density = primitive_soln[0];
203  const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive<real>(primitive_soln);
204 
205  std::array<real, nstate> conservative_soln;
206  conservative_soln[0] = density;
207  for (int d=0; d<dim; ++d) {
208  conservative_soln[1+d] = density*velocities[d];
209  }
210  conservative_soln[nstate-1] = compute_total_energy(primitive_soln);
211 
212  return conservative_soln;
213 }
214 
215 //template <int dim, int nstate, typename real>
216 //inline dealii::Tensor<1,dim,double> Euler<dim,nstate,real>::compute_velocities_inf() const
217 //{
218 // dealii::Tensor<1,dim,double> velocities;
219 // return velocities;
220 //}
221 
222 template <int dim, int nstate, typename real>
223 template<typename real2>
224 inline dealii::Tensor<1,dim,real2> Euler<dim,nstate,real>
225 ::compute_velocities ( const std::array<real2,nstate> &conservative_soln ) const
226 {
227  const real2 density = conservative_soln[0];
228  dealii::Tensor<1,dim,real2> vel;
229  for (int d=0; d<dim; ++d) { vel[d] = conservative_soln[1+d]/density; }
230  return vel;
231 }
232 
233 template <int dim, int nstate, typename real>
234 template <typename real2>
235 inline real2 Euler<dim,nstate,real>
236 ::compute_velocity_squared ( const dealii::Tensor<1,dim,real2> &velocities ) const
237 {
238  real2 vel2 = 0.0;
239  for (int d=0; d<dim; d++) {
240  vel2 = vel2 + velocities[d]*velocities[d];
241  }
242 
243  return vel2;
244 }
245 
246 template <int dim, int nstate, typename real>
247 template<typename real2>
248 inline dealii::Tensor<1,dim,real2> Euler<dim,nstate,real>
249 ::extract_velocities_from_primitive ( const std::array<real2,nstate> &primitive_soln ) const
250 {
251  dealii::Tensor<1,dim,real2> velocities;
252  for (int d=0; d<dim; d++) { velocities[d] = primitive_soln[1+d]; }
253  return velocities;
254 }
255 
256 template <int dim, int nstate, typename real>
257 inline real Euler<dim,nstate,real>
258 ::compute_total_energy ( const std::array<real,nstate> &primitive_soln ) const
259 {
260  const real pressure = primitive_soln[nstate-1];
261  const real kinetic_energy = compute_kinetic_energy_from_primitive_solution(primitive_soln);
262  const real tot_energy = pressure / this->gamm1 + kinetic_energy;
263  return tot_energy;
264 }
265 
266 template <int dim, int nstate, typename real>
267 inline real Euler<dim,nstate,real>
268 ::compute_kinetic_energy_from_primitive_solution ( const std::array<real,nstate> &primitive_soln ) const
269 {
270  const real density = primitive_soln[0];
271  const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive<real>(primitive_soln);
272  const real vel2 = compute_velocity_squared<real>(velocities);
273  const real kinetic_energy = 0.5*density*vel2;
274  return kinetic_energy;
275 }
276 
277 template <int dim, int nstate, typename real>
278 inline real Euler<dim,nstate,real>
279 ::compute_kinetic_energy_from_conservative_solution ( const std::array<real,nstate> &conservative_soln ) const
280 {
281  const std::array<real,nstate> primitive_soln = convert_conservative_to_primitive<real>(conservative_soln);
282  const real kinetic_energy = compute_kinetic_energy_from_primitive_solution(primitive_soln);
283  return kinetic_energy;
284 }
285 
286 template <int dim, int nstate, typename real>
287 inline real Euler<dim,nstate,real>
288 ::compute_entropy_measure ( const std::array<real,nstate> &conservative_soln ) const
289 {
290  real density = conservative_soln[0];
291  const real pressure = compute_pressure<real>(conservative_soln);
292  return compute_entropy_measure(density, pressure);
293 }
294 
295 template <int dim, int nstate, typename real>
296 inline real Euler<dim,nstate,real>
297 ::compute_entropy_measure ( const real density, const real pressure ) const
298 {
299  //Copy such that we don't modify the original density that is passed
300  real density_check = density;
301  const bool density_is_positive = check_positive_quantity<real>(density_check, "density");
302  if (density_is_positive) return pressure*pow(density,-gam);
303  else return (real)this->BIG_NUMBER;
304 }
305 
306 
307 template <int dim, int nstate, typename real>
308 inline real Euler<dim,nstate,real>
309 ::compute_specific_enthalpy ( const std::array<real,nstate> &conservative_soln, const real pressure ) const
310 {
311  const real density = conservative_soln[0];
312  const real total_energy = conservative_soln[nstate-1];
313  const real specific_enthalpy = (total_energy+pressure)/density;
314  return specific_enthalpy;
315 }
316 
317 template <int dim, int nstate, typename real>
318 inline real Euler<dim,nstate,real>
319 ::compute_numerical_entropy_function ( const std::array<real,nstate> &conservative_soln ) const
320 {
321  const real pressure = compute_pressure<real>(conservative_soln);
322  const real density = conservative_soln[0];
323 
324  const real entropy = compute_entropy<real>(density, pressure);
325 
326  const real numerical_entropy_function = - density * entropy;
327 
328  return numerical_entropy_function;
329 }
330 
331 template <int dim, int nstate, typename real>
332 template<typename real2>
333 inline real2 Euler<dim,nstate,real>
334 ::compute_temperature ( const std::array<real2,nstate> &primitive_soln ) const
335 {
336  const real2 density = primitive_soln[0];
337  const real2 pressure = primitive_soln[nstate-1];
338  const real2 temperature = gam*mach_inf_sqr*(pressure/density);
339  return temperature;
340 }
341 
342 template <int dim, int nstate, typename real>
343 inline real Euler<dim,nstate,real>
344 ::compute_density_from_pressure_temperature ( const real pressure, const real temperature ) const
345 {
346  const real density = gam*mach_inf_sqr*(pressure/temperature);
347  return density;
348 }
349 
350 template <int dim, int nstate, typename real>
351 inline real Euler<dim,nstate,real>
352 ::compute_temperature_from_density_pressure ( const real density, const real pressure ) const
353 {
354  const real temperature = gam*mach_inf_sqr*(pressure/density);
355  return temperature;
356 }
357 
358 template <int dim, int nstate, typename real>
359 inline real Euler<dim,nstate,real>
360 ::compute_pressure_from_density_temperature ( const real density, const real temperature ) const
361 {
362  const real pressure = density*temperature/(gam*mach_inf_sqr);
363  return pressure;
364 }
365 
366 template <int dim, int nstate, typename real>
367 template<typename real2>
368 inline real2 Euler<dim,nstate,real>
369 ::compute_pressure ( const std::array<real2,nstate> &conservative_soln ) const
370 {
371  const real2 density = conservative_soln[0];
372 
373  const real2 tot_energy = conservative_soln[nstate-1];
374 
375  const dealii::Tensor<1,dim,real2> vel = compute_velocities<real2>(conservative_soln);
376 
377  const real2 vel2 = compute_velocity_squared<real2>(vel);
378  real2 pressure = gamm1*(tot_energy - 0.5*density*vel2);
379 
380  check_positive_quantity<real2>(pressure, "pressure");
381  return pressure;
382 }
383 
384 
385 template <int dim, int nstate, typename real>
386 template<typename real2>
387 inline real2 Euler<dim,nstate,real>
388 :: compute_entropy (const real2 density, const real2 pressure) const
389 {
390  // copy density and pressure such that the check will not modify originals
391  real2 density_check = density;
392  real2 pressure_check = pressure;
393  const bool density_is_positive = check_positive_quantity(density_check, "density");
394  const bool pressure_is_positive = check_positive_quantity(pressure_check, "pressure");
395  if (density_is_positive && pressure_is_positive) {
396  real2 entropy = pressure * pow(density, -gam);
397  entropy = log(entropy);
398  return entropy;
399  } else {
400  return (real2)this->BIG_NUMBER;
401  }
402 
403 }
404 
405 template <int dim, int nstate, typename real>
406 inline real Euler<dim,nstate,real>
407 ::compute_sound ( const std::array<real,nstate> &conservative_soln ) const
408 {
409  real density = conservative_soln[0];
410  check_positive_quantity<real>(density, "density");
411  const real pressure = compute_pressure<real>(conservative_soln);
412  const real sound = sqrt(pressure*gam/density);
413  return sound;
414 }
415 
416 template <int dim, int nstate, typename real>
417 inline real Euler<dim,nstate,real>
418 ::compute_sound ( const real density, const real pressure ) const
419 {
420  //assert(density > 0);
421  const real sound = sqrt(pressure*gam/density);
422  return sound;
423 }
424 
425 template <int dim, int nstate, typename real>
426 inline real Euler<dim,nstate,real>
427 ::compute_mach_number ( const std::array<real,nstate> &conservative_soln ) const
428 {
429  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
430  const real velocity = sqrt(compute_velocity_squared<real>(vel));
431  const real sound = compute_sound (conservative_soln);
432  const real mach_number = velocity/sound;
433  return mach_number;
434 }
435 
436 // Split form functions:
437 
438 template <int dim, int nstate, typename real>
439 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim, nstate, real>
440 ::convective_numerical_split_flux(const std::array<real,nstate> &conservative_soln1,
441  const std::array<real,nstate> &conservative_soln2) const
442 {
443  std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
444  if(two_point_num_flux_type == two_point_num_flux_enum::KG) {
445  conv_num_split_flux = convective_numerical_split_flux_kennedy_gruber(conservative_soln1, conservative_soln2);
446  } else if(two_point_num_flux_type == two_point_num_flux_enum::IR) {
447  conv_num_split_flux = convective_numerical_split_flux_ismail_roe(conservative_soln1, conservative_soln2);
448  } else if(two_point_num_flux_type == two_point_num_flux_enum::CH) {
449  conv_num_split_flux = convective_numerical_split_flux_chandrashekar(conservative_soln1, conservative_soln2);
450  } else if(two_point_num_flux_type == two_point_num_flux_enum::Ra) {
451  conv_num_split_flux = convective_numerical_split_flux_ranocha(conservative_soln1, conservative_soln2);
452  }
453 
454  return conv_num_split_flux;
455 }
456 
457 template <int dim, int nstate, typename real>
458 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim, nstate, real>
459 ::convective_numerical_split_flux_kennedy_gruber(const std::array<real,nstate> &conservative_soln1,
460  const std::array<real,nstate> &conservative_soln2) const
461 {
462  std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
463  const real mean_density = compute_mean_density(conservative_soln1, conservative_soln2);
464  const real mean_pressure = compute_mean_pressure(conservative_soln1, conservative_soln2);
465  const dealii::Tensor<1,dim,real> mean_velocities = compute_mean_velocities(conservative_soln1,conservative_soln2);
466  const real mean_specific_total_energy = compute_mean_specific_total_energy(conservative_soln1, conservative_soln2);
467 
468  for (int flux_dim = 0; flux_dim < dim; ++flux_dim)
469  {
470  // Density equation
471  conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];
472  // Momentum equation
473  for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
474  conv_num_split_flux[1+velocity_dim][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_velocities[velocity_dim];
475  }
476  conv_num_split_flux[1+flux_dim][flux_dim] += mean_pressure; // Add diagonal of pressure
477  // Energy equation
478  conv_num_split_flux[nstate-1][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_specific_total_energy + mean_pressure * mean_velocities[flux_dim];
479  }
480 
481  return conv_num_split_flux;
482 }
483 
484 template <int dim, int nstate, typename real>
485 std::array<real,nstate> Euler<dim, nstate, real>
486 ::compute_ismail_roe_parameter_vector_from_primitive(const std::array<real,nstate> &primitive_soln) const
487 {
488  // Ismail-Roe parameter vector; Eq (3.14) [Gassner, Winters, and Kopriva, 2016, SBP]
489  std::array<real,nstate> ismail_roe_parameter_vector;
490  ismail_roe_parameter_vector[0] = sqrt(primitive_soln[0]/primitive_soln[nstate-1]);
491  for(int d=0; d<dim; ++d){
492  ismail_roe_parameter_vector[1+d] = ismail_roe_parameter_vector[0]*primitive_soln[1+d];
493  }
494  ismail_roe_parameter_vector[nstate-1] = sqrt(primitive_soln[0]*primitive_soln[nstate-1]);
495 
496  return ismail_roe_parameter_vector;
497 }
498 
499 template <int dim, int nstate, typename real>
500 real Euler<dim, nstate, real>
501 ::compute_ismail_roe_logarithmic_mean(const real val1, const real val2) const
502 {
503  // See Appendix B [Ismail and Roe, 2009, Entropy-Consistent Euler Flux Functions II]
504  // -- Numerically stable algorithm for computing the logarithmic mean
505  const real zeta = val1/val2;
506  const real f = (zeta-1.0)/(zeta+1.0);
507  const real u = f*f;
508 
509  real F;
510  if(u<1.0e-2){ F = 1.0 + u/3.0 + u*u/5.0 + u*u*u/7.0; }
511  else {
512  if constexpr(std::is_same<real,double>::value) F = std::log(zeta)/2.0/f;
513  }
514 
515  const real log_mean_val = (val1+val2)/(2.0*F);
516 
517  return log_mean_val;
518 }
519 
520 template <int dim, int nstate, typename real>
521 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim, nstate, real>
522 ::convective_numerical_split_flux_ismail_roe(const std::array<real,nstate> &conservative_soln1,
523  const std::array<real,nstate> &conservative_soln2) const
524 {
525  // Get Ismail Roe parameter vectors
526  const std::array<real,nstate> parameter_vector1 = compute_ismail_roe_parameter_vector_from_primitive(
527  convert_conservative_to_primitive<real>(conservative_soln1));
528  const std::array<real,nstate> parameter_vector2 = compute_ismail_roe_parameter_vector_from_primitive(
529  convert_conservative_to_primitive<real>(conservative_soln2));
530 
531  // Compute mean (average) parameter vector
532  std::array<real,nstate> avg_parameter_vector;
533  for(int s=0; s<nstate; ++s){
534  avg_parameter_vector[s] = 0.5*(parameter_vector1[s] + parameter_vector2[s]);
535  }
536 
537  // Compute logarithmic mean parameter vector
538  std::array<real,nstate> log_mean_parameter_vector;
539  for(int s=0; s<nstate; ++s){
540  log_mean_parameter_vector[s] = compute_ismail_roe_logarithmic_mean(parameter_vector1[s], parameter_vector2[s]);
541  }
542 
543  // Compute Ismail Roe mean primitive variables; Eq (3.15) [Gassner, Winters, and Kopriva, 2016, SBP]
544  std::array<real,dim> mean_velocities;
545  const real mean_density = avg_parameter_vector[0]*log_mean_parameter_vector[nstate-1];
546  for(int d=0; d<dim; ++d){
547  mean_velocities[d] = avg_parameter_vector[1+d]/avg_parameter_vector[0];
548  }
549  const real mean_pressure = avg_parameter_vector[nstate-1]/avg_parameter_vector[0];
550  // -- enthalpy
551  real mean_enthalpy = (gam+1.0)*(log_mean_parameter_vector[nstate-1]/log_mean_parameter_vector[0]) + gamm1*mean_pressure;
552  mean_enthalpy /= 2.0*gam;
553  mean_enthalpy *= gam/(mean_density*gamm1);
554  // -- get sum of mean velocities squared
555  real mean_velocities_sqr_sum = 0.0;
556  for(int d=0; d<dim; ++d){
557  mean_velocities_sqr_sum += mean_velocities[d]*mean_velocities[d];
558  }
559  // -- add to enthalpy
560  mean_enthalpy += 0.5*mean_velocities_sqr_sum;
561 
562  // Compute Ismail Roe convective numerical split flux
563  std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
564  for (int flux_dim = 0; flux_dim < dim; ++flux_dim)
565  {
566  // Density equation
567  conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];
568  // Momentum equation
569  for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
570  conv_num_split_flux[1+velocity_dim][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_velocities[velocity_dim];
571  }
572  conv_num_split_flux[1+flux_dim][flux_dim] += mean_pressure; // Add diagonal of pressure
573  // Energy equation
574  conv_num_split_flux[nstate-1][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_enthalpy;
575  }
576 
577  return conv_num_split_flux;
578 }
579 
580 template <int dim, int nstate, typename real>
581 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim, nstate, real>
582 ::convective_numerical_split_flux_chandrashekar(const std::array<real,nstate> &conservative_soln1,
583  const std::array<real,nstate> &conservative_soln2) const
584 {
585 
586  std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
587  const real rho_log = compute_ismail_roe_logarithmic_mean(conservative_soln1[0], conservative_soln2[0]);
588  const real pressure1 = compute_pressure<real>(conservative_soln1);
589  const real pressure2 = compute_pressure<real>(conservative_soln2);
590 
591  const real beta1 = conservative_soln1[0]/(2.0*pressure1);
592  const real beta2 = conservative_soln2[0]/(2.0*pressure2);
593 
594  const real beta_log = compute_ismail_roe_logarithmic_mean(beta1, beta2);
595  const dealii::Tensor<1,dim,real> vel1 = compute_velocities<real>(conservative_soln1);
596  const dealii::Tensor<1,dim,real> vel2 = compute_velocities<real>(conservative_soln2);
597 
598  const real pressure_hat = 0.5*(conservative_soln1[0] + conservative_soln2[0])/(2.0*0.5*(beta1+beta2));
599 
600  dealii::Tensor<1,dim,real> vel_avg;
601  real vel_square_avg = 0.0;;
602  for(int idim=0; idim<dim; idim++){
603  vel_avg[idim] = 0.5*(vel1[idim]+vel2[idim]);
604  vel_square_avg += (0.5 *(vel1[idim]+vel2[idim])) * (0.5 *(vel1[idim]+vel2[idim]));
605  }
606 
607  real enthalpy_hat = 1.0/(2.0*beta_log*gamm1) + vel_square_avg + pressure_hat/rho_log;
608 
609  for(int idim=0; idim<dim; idim++){
610  enthalpy_hat -= 0.5*(0.5*(vel1[idim]*vel1[idim] + vel2[idim]*vel2[idim]));
611  }
612 
613  for(int flux_dim=0; flux_dim<dim; flux_dim++){
614  // Density equation
615  conv_num_split_flux[0][flux_dim] = rho_log * vel_avg[flux_dim];
616  // Momentum equation
617  for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
618  conv_num_split_flux[1+velocity_dim][flux_dim] = rho_log*vel_avg[flux_dim]*vel_avg[velocity_dim];
619  }
620  conv_num_split_flux[1+flux_dim][flux_dim] += pressure_hat; // Add diagonal of pressure
621 
622  // Energy equation
623  conv_num_split_flux[nstate-1][flux_dim] = rho_log * vel_avg[flux_dim] * enthalpy_hat;
624  }
625 
626  return conv_num_split_flux;
627 }
628 
629 template <int dim, int nstate, typename real>
630 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim, nstate, real>
631 ::convective_numerical_split_flux_ranocha(const std::array<real,nstate> &conservative_soln1,
632  const std::array<real,nstate> &conservative_soln2) const
633 {
634 
635  std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
636  const real rho_log = compute_ismail_roe_logarithmic_mean(conservative_soln1[0], conservative_soln2[0]);
637  const real pressure1 = compute_pressure<real>(conservative_soln1);
638  const real pressure2 = compute_pressure<real>(conservative_soln2);
639 
640  const real beta1 = conservative_soln1[0]/(pressure1);
641  const real beta2 = conservative_soln2[0]/(pressure2);
642 
643  const real beta_log = compute_ismail_roe_logarithmic_mean(beta1, beta2);
644  const dealii::Tensor<1,dim,real> vel1 = compute_velocities<real>(conservative_soln1);
645  const dealii::Tensor<1,dim,real> vel2 = compute_velocities<real>(conservative_soln2);
646 
647  const real pressure_hat = 0.5*(pressure1+pressure2);
648 
649  dealii::Tensor<1,dim,real> vel_avg;
650  real vel_square_avg = 0.0;;
651  for(int idim=0; idim<dim; idim++){
652  vel_avg[idim] = 0.5*(vel1[idim]+vel2[idim]);
653  vel_square_avg += (0.5 *(vel1[idim]+vel2[idim])) * (0.5 *(vel1[idim]+vel2[idim]));
654  }
655 
656  real enthalpy_hat = 1.0/(beta_log*gamm1) + vel_square_avg + 2.0*pressure_hat/rho_log;
657 
658  for(int idim=0; idim<dim; idim++){
659  enthalpy_hat -= 0.5*(0.5*(vel1[idim]*vel1[idim] + vel2[idim]*vel2[idim]));
660  }
661 
662  for(int flux_dim=0; flux_dim<dim; flux_dim++){
663  // Density equation
664  conv_num_split_flux[0][flux_dim] = rho_log * vel_avg[flux_dim];
665  // Momentum equation
666  for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
667  conv_num_split_flux[1+velocity_dim][flux_dim] = rho_log*vel_avg[flux_dim]*vel_avg[velocity_dim];
668  }
669  conv_num_split_flux[1+flux_dim][flux_dim] += pressure_hat; // Add diagonal of pressure
670 
671  // Energy equation
672  conv_num_split_flux[nstate-1][flux_dim] = rho_log * vel_avg[flux_dim] * enthalpy_hat;
673  conv_num_split_flux[nstate-1][flux_dim] -= ( 0.5 *(pressure1*vel1[flux_dim] + pressure2*vel2[flux_dim]));
674  }
675 
676  return conv_num_split_flux;
677 
678 }
679 
680 template <int dim, int nstate, typename real>
681 std::array<real,nstate> Euler<dim, nstate, real>
682 ::compute_entropy_variables (
683  const std::array<real,nstate> &conservative_soln) const
684 {
685  std::array<real,nstate> entropy_var;
686  const real density = conservative_soln[0];
687  const real pressure = compute_pressure<real>(conservative_soln);
688 
689  const real entropy = compute_entropy<real>(density, pressure);
690 
691  const real rho_theta = pressure / gamm1;
692 
693  entropy_var[0] = (rho_theta *(gam + 1.0 - entropy) - conservative_soln[nstate-1])/rho_theta;
694  for(int idim=0; idim<dim; idim++){
695  entropy_var[idim+1] = conservative_soln[idim+1] / rho_theta;
696  }
697  entropy_var[nstate-1] = - density / rho_theta;
698 
699  return entropy_var;
700 }
701 
702 template <int dim, int nstate, typename real>
703 std::array<real,nstate> Euler<dim, nstate, real>
704 ::compute_conservative_variables_from_entropy_variables (
705  const std::array<real,nstate> &entropy_var) const
706 {
707  //Eq. 119 and 120 from Chan, Jesse. "On discretely entropy conservative and entropy stable discontinuous Galerkin methods." Journal of Computational Physics 362 (2018): 346-374.
708  //Extrapolated for 3D
709  std::array<real,nstate> conservative_var;
710  real entropy_var_vel_squared = 0.0;
711  for(int idim=0; idim<dim; idim++){
712  entropy_var_vel_squared += entropy_var[idim + 1] * entropy_var[idim + 1];
713  }
714  const real entropy = gam - entropy_var[0] + 0.5 * entropy_var_vel_squared / entropy_var[nstate-1];
715  const real rho_theta = pow( (gamm1/ pow(- entropy_var[nstate-1], gam)), 1.0 /gamm1)
716  * exp( - entropy / gamm1);
717 
718  conservative_var[0] = - rho_theta * entropy_var[nstate-1];
719  for(int idim=0; idim<dim; idim++){
720  conservative_var[idim+1] = rho_theta * entropy_var[idim+1];
721  }
722  conservative_var[nstate-1] = rho_theta * (1.0 - 0.5 * entropy_var_vel_squared / entropy_var[nstate-1]);
723  return conservative_var;
724 }
725 
726 template <int dim, int nstate, typename real>
727 std::array<real,nstate> Euler<dim, nstate, real>
728 ::compute_kinetic_energy_variables (
729  const std::array<real,nstate> &conservative_soln) const
730 {
731  std::array<real,nstate> kin_energy_var;
732  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
733  const real vel2 = compute_velocity_squared<real>(vel);
734 
735  kin_energy_var[0] = - 0.5 * vel2;
736  for(int idim=0; idim<dim; idim++){
737  kin_energy_var[idim+1] = vel[idim];
738  }
739  kin_energy_var[nstate-1] = 0;
740 
741  return kin_energy_var;
742 }
743 
744 template <int dim, int nstate, typename real>
745 inline real Euler<dim,nstate,real>::
746 compute_mean_density(const std::array<real,nstate> &conservative_soln1,
747  const std::array<real,nstate> &conservative_soln2) const
748 {
749  return (conservative_soln1[0] + conservative_soln2[0])/2.;
750 }
751 
752 template <int dim, int nstate, typename real>
753 inline real Euler<dim,nstate,real>::
754 compute_mean_pressure(const std::array<real,nstate> &conservative_soln1,
755  const std::array<real,nstate> &conservative_soln2) const
756 {
757  real pressure_1 = compute_pressure<real>(conservative_soln1);
758  real pressure_2 = compute_pressure<real>(conservative_soln2);
759  return (pressure_1 + pressure_2)/2.;
760 }
761 
762 template <int dim, int nstate, typename real>
763 inline dealii::Tensor<1,dim,real> Euler<dim,nstate,real>::
764 compute_mean_velocities(const std::array<real,nstate> &conservative_soln1,
765  const std::array<real,nstate> &conservative_soln2) const
766 {
767  dealii::Tensor<1,dim,real> vel_1 = compute_velocities<real>(conservative_soln1);
768  dealii::Tensor<1,dim,real> vel_2 = compute_velocities<real>(conservative_soln2);
769  dealii::Tensor<1,dim,real> mean_vel;
770  for (int d=0; d<dim; ++d) {
771  mean_vel[d] = 0.5*(vel_1[d]+vel_2[d]);
772  }
773  return mean_vel;
774 }
775 
776 template <int dim, int nstate, typename real>
777 inline real Euler<dim,nstate,real>::
778 compute_mean_specific_total_energy(const std::array<real,nstate> &conservative_soln1,
779  const std::array<real,nstate> &conservative_soln2) const
780 {
781  return ((conservative_soln1[nstate-1]/conservative_soln1[0]) + (conservative_soln2[nstate-1]/conservative_soln2[0]))/2.;
782 }
783 
784 
785 template <int dim, int nstate, typename real>
786 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim,nstate,real>
787 ::convective_flux (const std::array<real,nstate> &conservative_soln) const
788 {
789  std::array<dealii::Tensor<1,dim,real>,nstate> conv_flux;
790  const real density = conservative_soln[0];
791  const real pressure = compute_pressure<real>(conservative_soln);
792  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
793  const real specific_total_energy = conservative_soln[nstate-1]/conservative_soln[0];
794  const real specific_total_enthalpy = specific_total_energy + pressure/density;
795 
796  for (int flux_dim=0; flux_dim<dim; ++flux_dim) {
797  // Density equation
798  conv_flux[0][flux_dim] = conservative_soln[1+flux_dim];
799  // Momentum equation
800  for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
801  conv_flux[1+velocity_dim][flux_dim] = density*vel[flux_dim]*vel[velocity_dim];
802  }
803  conv_flux[1+flux_dim][flux_dim] += pressure; // Add diagonal of pressure
804  // Energy equation
805  conv_flux[nstate-1][flux_dim] = density*vel[flux_dim]*specific_total_enthalpy;
806  }
807  return conv_flux;
808 }
809 
810 template <int dim, int nstate, typename real>
811 std::array<real,nstate> Euler<dim,nstate,real>
812 ::convective_normal_flux (const std::array<real,nstate> &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const
813 {
814  std::array<real, nstate> conv_normal_flux;
815  const real density = conservative_soln[0];
816  const real pressure = compute_pressure<real>(conservative_soln);
817  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
818  real normal_vel = 0.0;
819  for (int d=0; d<dim; ++d) {
820  normal_vel += vel[d]*normal[d];
821  }
822  const real total_energy = conservative_soln[nstate-1];
823  const real specific_total_enthalpy = (total_energy + pressure) / density;
824 
825  const real rhoV = density*normal_vel;
826  // Density equation
827  conv_normal_flux[0] = rhoV;
828  // Momentum equation
829  for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
830  conv_normal_flux[1+velocity_dim] = rhoV*vel[velocity_dim] + normal[velocity_dim] * pressure;
831  }
832  // Energy equation
833  conv_normal_flux[nstate-1] = rhoV*specific_total_enthalpy;
834  return conv_normal_flux;
835 }
836 
837 template <int dim, int nstate, typename real>
838 dealii::Tensor<2,nstate,real> Euler<dim,nstate,real>
839 ::convective_flux_directional_jacobian (
840  const std::array<real,nstate> &conservative_soln,
841  const dealii::Tensor<1,dim,real> &normal) const
842 {
843  // See Blazek (year?) Appendix A.9 p. 429-430
844  // For Blazek (2001), see Appendix A.7 p. 419-420
845  // Alternatively, see Masatsuka 2018 "I do like CFD", p.77, eq.(3.6.8)
846  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
847  real vel_normal = 0.0;
848  for (int d=0;d<dim;d++) { vel_normal += vel[d] * normal[d]; }
849 
850  const real vel2 = compute_velocity_squared<real>(vel);
851  const real phi = 0.5*gamm1 * vel2;
852 
853  const real density = conservative_soln[0];
854  const real tot_energy = conservative_soln[nstate-1];
855  const real E = tot_energy / density;
856  const real a1 = gam*E-phi;
857  const real a2 = gam-1.0;
858  const real a3 = gam-2.0;
859 
860  dealii::Tensor<2,nstate,real> jacobian;
861  for (int d=0; d<dim; ++d) {
862  jacobian[0][1+d] = normal[d];
863  }
864  for (int row_dim=0; row_dim<dim; ++row_dim) {
865  jacobian[1+row_dim][0] = normal[row_dim]*phi - vel[row_dim] * vel_normal;
866  for (int col_dim=0; col_dim<dim; ++col_dim){
867  if (row_dim == col_dim) {
868  jacobian[1+row_dim][1+col_dim] = vel_normal - a3*normal[row_dim]*vel[row_dim];
869  } else {
870  jacobian[1+row_dim][1+col_dim] = normal[col_dim]*vel[row_dim] - a2*normal[row_dim]*vel[col_dim];
871  }
872  }
873  jacobian[1+row_dim][nstate-1] = normal[row_dim]*a2;
874  }
875  jacobian[nstate-1][0] = vel_normal*(phi-a1);
876  for (int d=0; d<dim; ++d){
877  jacobian[nstate-1][1+d] = normal[d]*a1 - a2*vel[d]*vel_normal;
878  }
879  jacobian[nstate-1][nstate-1] = gam*vel_normal;
880 
881  return jacobian;
882 }
883 
884 template <int dim, int nstate, typename real>
885 std::array<real,nstate> Euler<dim,nstate,real>
886 ::convective_eigenvalues (
887  const std::array<real,nstate> &conservative_soln,
888  const dealii::Tensor<1,dim,real> &normal) const
889 {
890  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
891  std::array<real,nstate> eig;
892  real vel_dot_n = 0.0;
893  for (int d=0;d<dim;++d) { vel_dot_n += vel[d]*normal[d]; };
894  for (int i=0; i<nstate; i++) {
895  eig[i] = vel_dot_n;
896  //eig[i] = advection_speed*normal;
897 
898  //eig[i] = 1.0;
899  //eig[i] = -1.0;
900  }
901  return eig;
902 }
903 
904 template <int dim, int nstate, typename real>
905 real Euler<dim,nstate,real>
906 ::max_convective_eigenvalue (const std::array<real,nstate> &conservative_soln) const
907 {
908  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
909 
910  const real sound = compute_sound (conservative_soln);
911 
912  real vel2 = compute_velocity_squared<real>(vel);
913 
914  const real max_eig = sqrt(vel2) + sound;
915 
916  return max_eig;
917 }
918 
919 template <int dim, int nstate, typename real>
920 real Euler<dim,nstate,real>
921 ::max_convective_normal_eigenvalue (
922  const std::array<real,nstate> &conservative_soln,
923  const dealii::Tensor<1,dim,real> &normal) const
924 {
925  const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
926 
927  const real sound = compute_sound (conservative_soln);
928  real vel_dot_n = 0.0;
929  for (int d=0;d<dim;++d) { vel_dot_n += vel[d]*normal[d]; };
930  const real max_normal_eig = abs(vel_dot_n) + sound;
931 
932  return max_normal_eig;
933 }
934 
935 template <int dim, int nstate, typename real>
936 real Euler<dim,nstate,real>
937 ::max_viscous_eigenvalue (const std::array<real,nstate> &/*conservative_soln*/) const
938 {
939  return 0.0;
940 }
941 
942 template <int dim, int nstate, typename real>
943 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim,nstate,real>
944 ::dissipative_flux (
945  const std::array<real,nstate> &conservative_soln,
946  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
947  const dealii::types::global_dof_index /*cell_index*/) const
948 {
949  return dissipative_flux(conservative_soln, solution_gradient);
950 }
951 
952 template <int dim, int nstate, typename real>
953 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim,nstate,real>
954 ::dissipative_flux (
955  const std::array<real,nstate> &/*conservative_soln*/,
956  const std::array<dealii::Tensor<1,dim,real>,nstate> &/*solution_gradient*/) const
957 {
958  std::array<dealii::Tensor<1,dim,real>,nstate> diss_flux;
959  // No dissipation for Euler
960  for (int i=0; i<nstate; i++) {
961  diss_flux[i] = 0;
962  }
963  return diss_flux;
964 }
965 
966 template <int dim, int nstate, typename real>
967 void Euler<dim,nstate,real>
968 ::boundary_riemann (
969  const dealii::Tensor<1,dim,real> &normal_int,
970  const std::array<real,nstate> &soln_int,
971  std::array<real,nstate> &soln_bc) const
972 {
973  std::array<real,nstate> primitive_int = convert_conservative_to_primitive<real>(soln_int);
974  std::array<real,nstate> primitive_ext;
975  primitive_ext[0] = density_inf;
976  for (int d=0;d<dim;d++) { primitive_ext[1+d] = velocities_inf[d]; }
977  primitive_ext[nstate-1] = pressure_inf;
978 
979  const dealii::Tensor<1,dim,real> velocities_int = extract_velocities_from_primitive<real>(primitive_int);
980  const dealii::Tensor<1,dim,real> velocities_ext = extract_velocities_from_primitive<real>(primitive_ext);
981 
982  const real sound_int = compute_sound ( primitive_int[0], primitive_int[nstate-1] );
983  const real sound_ext = compute_sound ( primitive_ext[0], primitive_ext[nstate-1] );
984 
985  real vel_int_dot_normal = 0.0;
986  real vel_ext_dot_normal = 0.0;
987  for (int d=0; d<dim; d++) {
988  vel_int_dot_normal = vel_int_dot_normal + velocities_int[d]*normal_int[d];
989  vel_ext_dot_normal = vel_ext_dot_normal + velocities_ext[d]*normal_int[d];
990  }
991 
992  // Riemann invariants
993  const real out_riemann_invariant = vel_int_dot_normal + 2.0/gamm1*sound_int, // Outgoing
994  inc_riemann_invariant = vel_ext_dot_normal - 2.0/gamm1*sound_ext; // Incoming
995 
996  const real normal_velocity_bc = 0.5*(out_riemann_invariant+inc_riemann_invariant),
997  sound_bc = 0.25*gamm1*(out_riemann_invariant-inc_riemann_invariant);
998 
999  std::array<real,nstate> primitive_bc;
1000  if (abs(normal_velocity_bc) >= abs(sound_bc)) { // Supersonic
1001  if (normal_velocity_bc < 0.0) { // Inlet
1002  primitive_bc = primitive_ext;
1003  } else { // Outlet
1004  primitive_bc = primitive_int;
1005  }
1006  } else { // Subsonic
1007 
1008  real density_bc;
1009  dealii::Tensor<1,dim,real> velocities_bc;
1010  real pressure_bc;
1011 
1012  dealii::Tensor<1,dim,real> velocities_tangential;
1013  if (normal_velocity_bc < 0.0) { // Inlet
1014  const real entropy_ext = compute_entropy_measure(primitive_ext[0], primitive_ext[nstate-1]);
1015  density_bc = pow( 1.0/gam * sound_bc * sound_bc / entropy_ext, 1.0/gamm1 );
1016  for (int d=0; d<dim; ++d) {
1017  velocities_tangential[d] = velocities_ext[d] - vel_ext_dot_normal * normal_int[d];
1018  }
1019  } else { // Outlet
1020  const real entropy_int = compute_entropy_measure(primitive_int[0], primitive_int[nstate-1]);
1021  density_bc = pow( 1.0/gam * sound_bc * sound_bc / entropy_int, 1.0/gamm1 );
1022  for (int d=0; d<dim; ++d) {
1023  velocities_tangential[d] = velocities_int[d] - vel_int_dot_normal * normal_int[d];
1024  }
1025  }
1026  for (int d=0; d<dim; ++d) {
1027  velocities_bc[d] = velocities_tangential[d] + normal_velocity_bc*normal_int[d];
1028  }
1029 
1030  pressure_bc = 1.0/gam * sound_bc * sound_bc * density_bc;
1031 
1032  primitive_bc[0] = density_bc;
1033  for (int d=0;d<dim;d++) { primitive_bc[1+d] = velocities_bc[d]; }
1034  primitive_bc[nstate-1] = pressure_bc;
1035  }
1036 
1037  soln_bc = convert_primitive_to_conservative(primitive_bc);
1038 }
1039 
1040 template <int dim, int nstate, typename real>
1041 void Euler<dim,nstate,real>
1042 ::boundary_slip_wall (
1043  const dealii::Tensor<1,dim,real> &normal_int,
1044  const std::array<real,nstate> &soln_int,
1045  const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
1046  std::array<real,nstate> &soln_bc,
1047  std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
1048 {
1049  // Slip wall boundary conditions (No penetration)
1050  // Given by Algorithm II of the following paper
1051  // Krivodonova, L., and Berger, M.,
1052  // “High-order accurate implementation of solid wall boundary conditions in curved geometries,”
1053  // Journal of Computational Physics, vol. 211, 2006, pp. 492–512.
1054  const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive<real>(soln_int);
1055 
1056  // Copy density and pressure
1057  std::array<real,nstate> primitive_boundary_values;
1058  primitive_boundary_values[0] = primitive_interior_values[0];
1059  primitive_boundary_values[nstate-1] = primitive_interior_values[nstate-1];
1060 
1061  const dealii::Tensor<1,dim,real> surface_normal = -normal_int;
1062  const dealii::Tensor<1,dim,real> velocities_int = extract_velocities_from_primitive<real>(primitive_interior_values);
1063  //const dealii::Tensor<1,dim,real> velocities_bc = velocities_int - 2.0*(velocities_int*surface_normal)*surface_normal;
1064  real vel_int_dot_normal = 0.0;
1065  for (int d=0; d<dim; d++) {
1066  vel_int_dot_normal = vel_int_dot_normal + velocities_int[d]*surface_normal[d];
1067  }
1068  dealii::Tensor<1,dim,real> velocities_bc;
1069  for (int d=0; d<dim; d++) {
1070  velocities_bc[d] = velocities_int[d] - 2.0*(vel_int_dot_normal)*surface_normal[d];
1071  //velocities_bc[d] = velocities_int[d] - (vel_int_dot_normal)*surface_normal[d];
1072  //velocities_bc[d] += velocities_int[d] * surface_normal.norm_square();
1073  }
1074  for (int d=0; d<dim; ++d) {
1075  primitive_boundary_values[1+d] = velocities_bc[d];
1076  }
1077 
1078  const std::array<real,nstate> modified_conservative_boundary_values = convert_primitive_to_conservative(primitive_boundary_values);
1079  for (int istate=0; istate<nstate; ++istate) {
1080  soln_bc[istate] = modified_conservative_boundary_values[istate];
1081  }
1082 
1083  for (int istate=0; istate<nstate; ++istate) {
1084  soln_grad_bc[istate] = -soln_grad_int[istate];
1085  }
1086 }
1087 
1088 template <int dim, int nstate, typename real>
1089 void Euler<dim,nstate,real>
1090 ::boundary_wall (
1091  const dealii::Tensor<1,dim,real> &normal_int,
1092  const std::array<real,nstate> &soln_int,
1093  const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
1094  std::array<real,nstate> &soln_bc,
1095  std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
1096 {
1097  // Slip wall boundary for Euler
1098  boundary_slip_wall(normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1099 }
1100 
1101 template <int dim, int nstate, typename real>
1102 void Euler<dim,nstate,real>
1103 ::boundary_manufactured_solution (
1104  const dealii::Point<dim, real> &pos,
1105  const dealii::Tensor<1,dim,real> &normal_int,
1106  const std::array<real,nstate> &soln_int,
1107  const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
1108  std::array<real,nstate> &soln_bc,
1109  std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
1110 {
1111  // Manufactured solution
1112  std::array<real,nstate> conservative_boundary_values;
1113  std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
1114  for (int s=0; s<nstate; s++) {
1115  conservative_boundary_values[s] = this->manufactured_solution_function->value (pos, s);
1116  boundary_gradients[s] = this->manufactured_solution_function->gradient (pos, s);
1117  }
1118  std::array<real,nstate> primitive_boundary_values = convert_conservative_to_primitive<real>(conservative_boundary_values);
1119  for (int istate=0; istate<nstate; ++istate) {
1120 
1121  std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(conservative_boundary_values, normal_int);
1122  const bool inflow = (characteristic_dot_n[istate] <= 0.);
1123 
1124  if (inflow) { // Dirichlet boundary condition
1125 
1126  soln_bc[istate] = conservative_boundary_values[istate];
1127  soln_grad_bc[istate] = soln_grad_int[istate];
1128 
1129  // Only set the pressure and velocity
1130  // primitive_boundary_values[0] = soln_int[0];;
1131  // for(int d=0;d<dim;d++){
1132  // primitive_boundary_values[1+d] = soln_int[1+d]/soln_int[0];;
1133  //}
1134  const std::array<real,nstate> modified_conservative_boundary_values = convert_primitive_to_conservative(primitive_boundary_values);
1135  (void) modified_conservative_boundary_values;
1136  //conservative_boundary_values[nstate-1] = soln_int[nstate-1];
1137  soln_bc[istate] = conservative_boundary_values[istate];
1138 
1139  } else { // Neumann boundary condition
1140  // soln_bc[istate] = -soln_int[istate]+2*conservative_boundary_values[istate];
1141  soln_bc[istate] = soln_int[istate];
1142 
1143  // **************************************************************************************************************
1144  // Note I don't know how to properly impose the soln_grad_bc to obtain an adjoint consistent scheme
1145  // Currently, Neumann boundary conditions are only imposed for the linear advection
1146  // Therefore, soln_grad_bc does not affect the solution
1147  // **************************************************************************************************************
1148  soln_grad_bc[istate] = soln_grad_int[istate];
1149  //soln_grad_bc[istate] = boundary_gradients[istate];
1150  //soln_grad_bc[istate] = -soln_grad_int[istate]+2*boundary_gradients[istate];
1151  }
1152 
1153  // HARDCODE DIRICHLET BC
1154  soln_bc[istate] = conservative_boundary_values[istate];
1155  }
1156 }
1157 
1158 template <int dim, int nstate, typename real>
1159 void Euler<dim,nstate,real>
1160 ::boundary_pressure_outflow (
1161  const real total_inlet_pressure,
1162  const real back_pressure,
1163  const std::array<real,nstate> &soln_int,
1164  std::array<real,nstate> &soln_bc) const
1165 {
1166  // Pressure Outflow Boundary Condition (back pressure)
1167  // Reference: Carlson 2011, sec. 2.4
1168 
1169  const real mach_int = compute_mach_number(soln_int);
1170  if (mach_int > 1.0) {
1171  // Supersonic, simply extrapolate
1172  for (int istate=0; istate<nstate; ++istate) {
1173  soln_bc[istate] = soln_int[istate];
1174  }
1175  }
1176  else {
1177  const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive<real>(soln_int);
1178  const real pressure_int = primitive_interior_values[nstate-1];
1179 
1180  const real radicant = 1.0+0.5*gamm1*mach_inf_sqr;
1181  const real pressure_inlet = total_inlet_pressure * pow(radicant, -gam/gamm1);
1182  const real pressure_bc = (mach_int >= 1) * pressure_int + (1-(mach_int >= 1)) * back_pressure*pressure_inlet;
1183  const real temperature_int = compute_temperature<real>(primitive_interior_values);
1184 
1185  // Assign primitive boundary values
1186  std::array<real,nstate> primitive_boundary_values;
1187  primitive_boundary_values[0] = compute_density_from_pressure_temperature(pressure_bc, temperature_int);
1188  for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = primitive_interior_values[1+d]; }
1189  primitive_boundary_values[nstate-1] = pressure_bc;
1190 
1191  const std::array<real,nstate> conservative_bc = convert_primitive_to_conservative(primitive_boundary_values);
1192  for (int istate=0; istate<nstate; ++istate) {
1193  soln_bc[istate] = conservative_bc[istate];
1194  }
1195  }
1196 }
1197 
1198 template <int dim, int nstate, typename real>
1199 void Euler<dim,nstate,real>
1200 ::boundary_inflow (
1201  const real total_inlet_pressure,
1202  const real total_inlet_temperature,
1203  const dealii::Tensor<1,dim,real> &normal_int,
1204  const std::array<real,nstate> &soln_int,
1205  std::array<real,nstate> &soln_bc) const
1206 {
1207  // Inflow boundary conditions (both subsonic and supersonic)
1208  // Carlson 2011, sec. 2.2 & sec 2.9
1209 
1210  const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive<real>(soln_int);
1211 
1212  const dealii::Tensor<1,dim,real> normal = -normal_int;
1213 
1214  const real density_i = primitive_interior_values[0];
1215  const dealii::Tensor<1,dim,real> velocities_i = extract_velocities_from_primitive<real>(primitive_interior_values);
1216  const real pressure_i = primitive_interior_values[nstate-1];
1217 
1218  //const real normal_vel_i = velocities_i*normal;
1219  real normal_vel_i = 0.0;
1220  for (int d=0; d<dim; ++d) {
1221  normal_vel_i += velocities_i[d]*normal[d];
1222  }
1223  const real sound_i = compute_sound(soln_int);
1224  //const real mach_i = std::abs(normal_vel_i)/sound_i;
1225 
1226  //const dealii::Tensor<1,dim,real> velocities_o = velocities_inf;
1227  //const real normal_vel_o = velocities_o*normal;
1228  //const real sound_o = sound_inf;
1229  //const real mach_o = mach_inf;
1230 
1231  if(mach_inf < 1.0) {
1232  // Subsonic inflow, sec 2.7
1233 
1234  //this->pcout << "Subsonic inflow, mach=" << mach_i << std::endl;
1235 
1236  // Want to solve for c_b (sound_bc), to then solve for U (velocity_magnitude_bc) and M_b (mach_bc)
1237  // Eq. 37
1238  const real riemann_pos = normal_vel_i + 2.0*sound_i/gamm1;
1239  // Could evaluate enthalpy from primitive like eq.36, but easier to use the following
1240  const real specific_total_energy = soln_int[nstate-1]/density_i;
1241  const real specific_total_enthalpy = specific_total_energy + pressure_i/density_i;
1242  // Eq. 43
1243  const real a = 1.0+2.0/gamm1;
1244  const real b = -2.0*riemann_pos;
1245  const real c = 0.5*gamm1 * (riemann_pos*riemann_pos - 2.0*specific_total_enthalpy);
1246  // Eq. 42
1247  const real term1 = -0.5*b/a;
1248  const real term2= 0.5*sqrt(b*b-4.0*a*c)/a;
1249  const real sound_bc1 = term1 + term2;
1250  const real sound_bc2 = term1 - term2;
1251  // Eq. 44
1252  const real sound_bc = std::max(sound_bc1, sound_bc2);
1253  // Eq. 45
1254  //const real velocity_magnitude_bc = 2.0*sound_bc/gamm1 - riemann_pos;
1255  const real velocity_magnitude_bc = riemann_pos - 2.0*sound_bc/gamm1;
1256  const real mach_bc = velocity_magnitude_bc/sound_bc;
1257  // Eq. 46
1258  const real radicant = 1.0+0.5*gamm1*mach_bc*mach_bc;
1259  const real pressure_bc = total_inlet_pressure * pow(radicant, -gam/gamm1);
1260  const real temperature_bc = total_inlet_temperature * pow(radicant, -1.0);
1261  //this->pcout << " pressure_bc " << pressure_bc << "pressure_inf" << pressure_inf << std::endl;
1262  //this->pcout << " temperature_bc " << temperature_bc << "temperature_inf" << temperature_inf << std::endl;
1263 
1264  const real density_bc = compute_density_from_pressure_temperature(pressure_bc, temperature_bc);
1265  std::array<real,nstate> primitive_boundary_values;
1266  primitive_boundary_values[0] = density_bc;
1267  for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = velocity_magnitude_bc*normal[d]; }
1268  primitive_boundary_values[nstate-1] = pressure_bc;
1269  const std::array<real,nstate> conservative_bc = convert_primitive_to_conservative(primitive_boundary_values);
1270  for (int istate=0; istate<nstate; ++istate) {
1271  soln_bc[istate] = conservative_bc[istate];
1272  }
1273 
1274  //this->pcout << " entropy_bc " << compute_entropy_measure(soln_bc) << "entropy_inf" << entropy_inf << std::endl;
1275 
1276  }
1277  else {
1278  // Supersonic inflow, sec 2.9
1279 
1280  // Specify all quantities through
1281  // total_inlet_pressure, total_inlet_temperature, mach_inf & angle_of_attack
1282  //this->pcout << "Supersonic inflow, mach=" << mach_i << std::endl;
1283  const real radicant = 1.0+0.5*gamm1*mach_inf_sqr;
1284  const real static_inlet_pressure = total_inlet_pressure * pow(radicant, -gam/gamm1);
1285  const real static_inlet_temperature = total_inlet_temperature * pow(radicant, -1.0);
1286 
1287  const real pressure_bc = static_inlet_pressure;
1288  const real temperature_bc = static_inlet_temperature;
1289  const real density_bc = compute_density_from_pressure_temperature(pressure_bc, temperature_bc);
1290  const real sound_bc = sqrt(gam * pressure_bc / density_bc);
1291  const real velocity_magnitude_bc = mach_inf * sound_bc;
1292 
1293  // Assign primitive boundary values
1294  std::array<real,nstate> primitive_boundary_values;
1295  primitive_boundary_values[0] = density_bc;
1296  for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = -velocity_magnitude_bc*normal_int[d]; } // minus since it's inflow
1297  primitive_boundary_values[nstate-1] = pressure_bc;
1298  const std::array<real,nstate> conservative_bc = convert_primitive_to_conservative(primitive_boundary_values);
1299  for (int istate=0; istate<nstate; ++istate) {
1300  soln_bc[istate] = conservative_bc[istate];
1301  }
1302  }
1303 }
1304 
1305 template <int dim, int nstate, typename real>
1306 void Euler<dim,nstate,real>
1307 ::boundary_farfield (
1308  std::array<real,nstate> &soln_bc) const
1309 {
1310  const real density_bc = density_inf;
1311  const real pressure_bc = 1.0/(gam*mach_inf_sqr);
1312  std::array<real,nstate> primitive_boundary_values;
1313  primitive_boundary_values[0] = density_bc;
1314  for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = velocities_inf[d]; } // minus since it's inflow
1315  primitive_boundary_values[nstate-1] = pressure_bc;
1316  const std::array<real,nstate> conservative_bc = convert_primitive_to_conservative(primitive_boundary_values);
1317  for (int istate=0; istate<nstate; ++istate) {
1318  soln_bc[istate] = conservative_bc[istate];
1319  }
1320 }
1321 
1322 template <int dim, int nstate, typename real>
1323 void Euler<dim, nstate, real>
1324 ::boundary_p0_extrapolation(
1325  const std::array<real, nstate>& soln_int,
1326  std::array<real, nstate>& soln_bc,
1327  std::array<dealii::Tensor<1, dim, real>, nstate>& soln_grad_bc) const
1328 {
1329  for (int istate = 0; istate < nstate; ++istate) {
1330  soln_bc[istate] = soln_int[istate];
1331  soln_grad_bc[istate] = 0;
1332  }
1333 
1334 }
1335 
1336 template <int dim, int nstate, typename real>
1337 void Euler<dim, nstate, real>
1338 ::boundary_custom(
1339  std::array<real, nstate>& soln_bc) const
1340 {
1341  std::array<real, nstate> primitive_boundary_values;
1342  for (int istate = 0; istate < nstate; ++istate) {
1343  primitive_boundary_values[istate] = this->all_parameters->euler_param.custom_boundary_for_each_state[istate];
1344  }
1345 
1346  const std::array<real, nstate> conservative_bc = convert_primitive_to_conservative(primitive_boundary_values);
1347  for (int istate = 0; istate < nstate; ++istate) {
1348  soln_bc[istate] = conservative_bc[istate];
1349  }
1350 }
1351 
1352 template <int dim, int nstate, typename real>
1353 void Euler<dim, nstate, real>
1354 ::boundary_astrophysical_inflow(
1355  std::array<real, nstate>& soln_bc) const
1356 {
1357  std::array<real, nstate> primitive_boundary_values;
1358  for (int istate = 0; istate < nstate; ++istate) {
1359  if(istate == 0)
1360  primitive_boundary_values[istate] = 0.5;
1361  if(istate == 1)
1362  primitive_boundary_values[istate] = 0.0;
1363  if(istate == 2)
1364  primitive_boundary_values[istate] = 0.0;
1365  if(istate == 3)
1366  primitive_boundary_values[istate] = 0.4127;
1367  if(istate == 4)
1368  primitive_boundary_values[istate] = 0.0;
1369  }
1370 
1371  const std::array<real, nstate> conservative_bc = convert_primitive_to_conservative(primitive_boundary_values);
1372  for (int istate = 0; istate < nstate; ++istate) {
1373  soln_bc[istate] = conservative_bc[istate];
1374  }
1375 }
1376 
1377 template <int dim, int nstate, typename real>
1378 void Euler<dim,nstate,real>
1379 ::boundary_face_values (
1380  const int boundary_type,
1381  const dealii::Point<dim, real> &pos,
1382  const dealii::Tensor<1,dim,real> &normal_int,
1383  const std::array<real,nstate> &soln_int,
1384  const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
1385  std::array<real,nstate> &soln_bc,
1386  std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
1387 {
1388  // NEED TO PROVIDE AS INPUT ************************************** (ask Doug where this should be moved to, protected member?)
1389  const real total_inlet_pressure = pressure_inf*pow(1.0+0.5*gamm1*mach_inf_sqr, gam/gamm1);
1390  const real total_inlet_temperature = temperature_inf*pow(total_inlet_pressure/pressure_inf, gamm1/gam);
1391 
1392  if (boundary_type == 1000) {
1393  // Manufactured solution boundary condition
1394  boundary_manufactured_solution (pos, normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1395  }
1396  else if (boundary_type == 1001) {
1397  // Wall boundary condition (slip for Euler, no-slip for Navier-Stokes; done through polymorphism)
1398  boundary_wall (normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1399  }
1400  else if (boundary_type == 1002) {
1401  // Pressure outflow boundary condition (back pressure)
1402  const real back_pressure = 0.99;
1403  boundary_pressure_outflow (total_inlet_pressure, back_pressure, soln_int, soln_bc);
1404  }
1405  else if (boundary_type == 1003) {
1406  // Inflow boundary condition
1407  boundary_inflow (total_inlet_pressure, total_inlet_temperature, normal_int, soln_int, soln_bc);
1408  }
1409  else if (boundary_type == 1004) {
1410  // Riemann-based farfield boundary condition
1411  boundary_riemann (normal_int, soln_int, soln_bc);
1412  }
1413  else if (boundary_type == 1005) {
1414  // Simple farfield boundary condition
1415  boundary_farfield(soln_bc);
1416  }
1417  else if (boundary_type == 1006) {
1418  // Slip wall boundary condition
1419  boundary_slip_wall (normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1420  }
1421  else if (boundary_type == 1007) {
1422  // Do nothing bc, p0 interpolation
1423  boundary_p0_extrapolation (soln_int, soln_bc, soln_grad_bc);
1424  }
1425  else if (boundary_type == 1008) {
1426  // Custom boundary condition, user defined in parameters
1427  boundary_custom (soln_bc);
1428  }
1429  else if (boundary_type == 1009) {
1430  // Boundary specific to the the astrophysical jet case
1431  boundary_astrophysical_inflow (soln_bc);
1432  }
1433  else {
1434  this->pcout << "Invalid boundary_type: " << boundary_type << std::endl;
1435  std::abort();
1436  }
1437 }
1438 
1439 template <int dim, int nstate, typename real>
1440 dealii::Vector<double> Euler<dim,nstate,real>::post_compute_derived_quantities_vector (
1441  const dealii::Vector<double> &uh,
1442  const std::vector<dealii::Tensor<1,dim> > &duh,
1443  const std::vector<dealii::Tensor<2,dim> > &dduh,
1444  const dealii::Tensor<1,dim> &normals,
1445  const dealii::Point<dim> &evaluation_points) const
1446 {
1447  std::vector<std::string> names = post_get_names ();
1448  dealii::Vector<double> computed_quantities = PhysicsBase<dim,nstate,real>::post_compute_derived_quantities_vector ( uh, duh, dduh, normals, evaluation_points);
1449  unsigned int current_data_index = computed_quantities.size() - 1;
1450  computed_quantities.grow_or_shrink(names.size());
1451  if constexpr (std::is_same<real,double>::value) {
1452 
1453  std::array<double, nstate> conservative_soln;
1454  for (unsigned int s=0; s<nstate; ++s) {
1455  conservative_soln[s] = uh(s);
1456  }
1457  const std::array<double, nstate> primitive_soln = convert_conservative_to_primitive<real>(conservative_soln);
1458  // if (primitive_soln[0] < 0) this->pcout << evaluation_points << std::endl;
1459 
1460  // Density
1461  computed_quantities(++current_data_index) = primitive_soln[0];
1462  // Velocities
1463  for (unsigned int d=0; d<dim; ++d) {
1464  computed_quantities(++current_data_index) = primitive_soln[1+d];
1465  }
1466  // Momentum
1467  for (unsigned int d=0; d<dim; ++d) {
1468  computed_quantities(++current_data_index) = conservative_soln[1+d];
1469  }
1470  // Energy
1471  computed_quantities(++current_data_index) = conservative_soln[nstate-1];
1472  // Pressure
1473  computed_quantities(++current_data_index) = primitive_soln[nstate-1];
1474  // Pressure coefficient
1475  computed_quantities(++current_data_index) = (primitive_soln[nstate-1] - pressure_inf) / dynamic_pressure_inf;
1476  // Temperature
1477  computed_quantities(++current_data_index) = compute_temperature<real>(primitive_soln);
1478  // Entropy generation
1479  computed_quantities(++current_data_index) = compute_entropy_measure(conservative_soln) - entropy_inf;
1480  // Mach Number
1481  computed_quantities(++current_data_index) = compute_mach_number(conservative_soln);
1482 
1483  }
1484  if (computed_quantities.size()-1 != current_data_index) {
1485  this->pcout << " Did not assign a value to all the data. Missing " << computed_quantities.size() - current_data_index << " variables."
1486  << " If you added a new output variable, make sure the names and DataComponentInterpretation match the above. "
1487  << std::endl;
1488  }
1489 
1490  return computed_quantities;
1491 }
1492 
1493 template <int dim, int nstate, typename real>
1494 std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> Euler<dim,nstate,real>
1495 ::post_get_data_component_interpretation () const
1496 {
1497  namespace DCI = dealii::DataComponentInterpretation;
1498  std::vector<DCI::DataComponentInterpretation> interpretation = PhysicsBase<dim,nstate,real>::post_get_data_component_interpretation (); // state variables
1499  interpretation.push_back (DCI::component_is_scalar); // Density
1500  for (unsigned int d=0; d<dim; ++d) {
1501  interpretation.push_back (DCI::component_is_part_of_vector); // Velocity
1502  }
1503  for (unsigned int d=0; d<dim; ++d) {
1504  interpretation.push_back (DCI::component_is_part_of_vector); // Momentum
1505  }
1506  interpretation.push_back (DCI::component_is_scalar); // Energy
1507  interpretation.push_back (DCI::component_is_scalar); // Pressure
1508  interpretation.push_back (DCI::component_is_scalar); // Pressure coefficient
1509  interpretation.push_back (DCI::component_is_scalar); // Temperature
1510  interpretation.push_back (DCI::component_is_scalar); // Entropy generation
1511  interpretation.push_back (DCI::component_is_scalar); // Mach number
1512 
1513  std::vector<std::string> names = post_get_names();
1514  if (names.size() != interpretation.size()) {
1515  this->pcout << "Number of DataComponentInterpretation is not the same as number of names for output file" << std::endl;
1516  }
1517  return interpretation;
1518 }
1519 
1520 
1521 template <int dim, int nstate, typename real>
1522 std::vector<std::string> Euler<dim,nstate,real>
1523 ::post_get_names () const
1524 {
1525  std::vector<std::string> names = PhysicsBase<dim,nstate,real>::post_get_names ();
1526  names.push_back ("density");
1527  for (unsigned int d=0; d<dim; ++d) {
1528  names.push_back ("velocity");
1529  }
1530  for (unsigned int d=0; d<dim; ++d) {
1531  names.push_back ("momentum");
1532  }
1533  names.push_back ("energy");
1534  names.push_back ("pressure");
1535  names.push_back ("pressure_coeffcient");
1536  names.push_back ("temperature");
1537 
1538  names.push_back ("entropy_generation");
1539  names.push_back ("mach_number");
1540  return names;
1541 }
1542 
1543 template <int dim, int nstate, typename real>
1544 dealii::UpdateFlags Euler<dim,nstate,real>
1545 ::post_get_needed_update_flags () const
1546 {
1547  //return update_values | update_gradients;
1548  return dealii::update_values
1549  | dealii::update_quadrature_points
1550  ;
1551 }
1552 
1553 // Instantiate explicitly
1554 template class Euler < PHILIP_DIM, PHILIP_DIM+2, double >;
1555 template class Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >;
1556 template class Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >;
1557 template class Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >;
1558 template class Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >;
1559 
1560 //==============================================================================
1561 // -> Templated member functions: // could be automated later on using Boost MPL
1562 //------------------------------------------------------------------------------
1563 // -- check_positive_quantity
1564 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, double >::check_positive_quantity< double >(double &qty, const std::string qty_name) const;
1565 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >::check_positive_quantity< FadType >(FadType &qty, const std::string qty_name) const;
1566 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::check_positive_quantity< RadType >(RadType &qty, const std::string qty_name) const;
1567 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::check_positive_quantity< FadFadType >(FadFadType &qty, const std::string qty_name) const;
1568 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::check_positive_quantity< RadFadType >(RadFadType &qty, const std::string qty_name) const;
1569 // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
1570 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, double >::check_positive_quantity< FadType >(FadType &qty, const std::string qty_name) const;
1571 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::check_positive_quantity< FadType >(FadType &qty, const std::string qty_name) const;
1572 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::check_positive_quantity< FadType >(FadType &qty, const std::string qty_name) const;
1573 template bool Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::check_positive_quantity< FadType >(FadType &qty, const std::string qty_name) const;
1574 // -- compute_pressure()
1575 template double Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_pressure< double >(const std::array<double, PHILIP_DIM+2> &conservative_soln) const;
1576 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >::compute_pressure< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1577 template RadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_pressure< RadType >(const std::array<RadType, PHILIP_DIM+2> &conservative_soln) const;
1578 template FadFadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_pressure< FadFadType >(const std::array<FadFadType,PHILIP_DIM+2> &conservative_soln) const;
1579 template RadFadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_pressure< RadFadType >(const std::array<RadFadType,PHILIP_DIM+2> &conservative_soln) const;
1580 // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
1581 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_pressure< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1582 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_pressure< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1583 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_pressure< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1584 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_pressure< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1585 // -- compute_temperature()
1586 template double Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_temperature< double >(const std::array<double, PHILIP_DIM+2> &primitive_soln) const;
1587 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >::compute_temperature< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1588 template RadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_temperature< RadType >(const std::array<RadType, PHILIP_DIM+2> &primitive_soln) const;
1589 template FadFadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_temperature< FadFadType >(const std::array<FadFadType,PHILIP_DIM+2> &primitive_soln) const;
1590 template RadFadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_temperature< RadFadType >(const std::array<RadFadType,PHILIP_DIM+2> &primitive_soln) const;
1591 // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
1592 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_temperature< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1593 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_temperature< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1594 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_temperature< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1595 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_temperature< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1596 // -- compute_velocity_squared()
1597 template double Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_velocity_squared< double >(const dealii::Tensor<1,PHILIP_DIM,double > &velocities) const;
1598 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >::compute_velocity_squared< FadType >(const dealii::Tensor<1,PHILIP_DIM,FadType > &velocities) const;
1599 template RadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_velocity_squared< RadType >(const dealii::Tensor<1,PHILIP_DIM,RadType > &velocities) const;
1600 template FadFadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_velocity_squared< FadFadType >(const dealii::Tensor<1,PHILIP_DIM,FadFadType> &velocities) const;
1601 template RadFadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_velocity_squared< RadFadType >(const dealii::Tensor<1,PHILIP_DIM,RadFadType> &velocities) const;
1602 // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
1603 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_velocity_squared< FadType >(const dealii::Tensor<1,PHILIP_DIM,FadType > &velocities) const;
1604 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_velocity_squared< FadType >(const dealii::Tensor<1,PHILIP_DIM,FadType > &velocities) const;
1605 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_velocity_squared< FadType >(const dealii::Tensor<1,PHILIP_DIM,FadType > &velocities) const;
1606 template FadType Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_velocity_squared< FadType >(const dealii::Tensor<1,PHILIP_DIM,FadType > &velocities) const;
1607 // -- convert_conservative_to_primitive()
1608 template std::array<double, PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, double >::convert_conservative_to_primitive< double >(const std::array<double, PHILIP_DIM+2> &conservative_soln) const;
1609 template std::array<FadType, PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >::convert_conservative_to_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1610 template std::array<RadType, PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::convert_conservative_to_primitive< RadType >(const std::array<RadType, PHILIP_DIM+2> &conservative_soln) const;
1611 template std::array<FadFadType,PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::convert_conservative_to_primitive< FadFadType >(const std::array<FadFadType,PHILIP_DIM+2> &conservative_soln) const;
1612 template std::array<RadFadType,PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::convert_conservative_to_primitive< RadFadType >(const std::array<RadFadType,PHILIP_DIM+2> &conservative_soln) const;
1613 // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
1614 template std::array<FadType, PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, double >::convert_conservative_to_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1615 template std::array<FadType, PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::convert_conservative_to_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1616 template std::array<FadType, PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::convert_conservative_to_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1617 template std::array<FadType, PHILIP_DIM+2> Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::convert_conservative_to_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1618 // -- extract_velocities_from_primitive()
1619 template dealii::Tensor<1,PHILIP_DIM,double > Euler < PHILIP_DIM, PHILIP_DIM+2, double >::extract_velocities_from_primitive< double >(const std::array<double, PHILIP_DIM+2> &primitive_soln) const;
1620 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >::extract_velocities_from_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1621 template dealii::Tensor<1,PHILIP_DIM,RadType > Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::extract_velocities_from_primitive< RadType >(const std::array<RadType, PHILIP_DIM+2> &primitive_soln) const;
1622 template dealii::Tensor<1,PHILIP_DIM,FadFadType> Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::extract_velocities_from_primitive< FadFadType >(const std::array<FadFadType,PHILIP_DIM+2> &primitive_soln) const;
1623 template dealii::Tensor<1,PHILIP_DIM,RadFadType> Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::extract_velocities_from_primitive< RadFadType >(const std::array<RadFadType,PHILIP_DIM+2> &primitive_soln) const;
1624 // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
1625 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, double >::extract_velocities_from_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1626 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::extract_velocities_from_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1627 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::extract_velocities_from_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1628 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::extract_velocities_from_primitive< FadType >(const std::array<FadType, PHILIP_DIM+2> &primitive_soln) const;
1629 // -- compute_velocities()
1630 template dealii::Tensor<1,PHILIP_DIM,double > Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_velocities< double >(const std::array<double, PHILIP_DIM+2> &conservative_soln) const;
1631 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, FadType >::compute_velocities< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1632 template dealii::Tensor<1,PHILIP_DIM,RadType > Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_velocities< RadType >(const std::array<RadType, PHILIP_DIM+2> &conservative_soln) const;
1633 template dealii::Tensor<1,PHILIP_DIM,FadFadType> Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_velocities< FadFadType >(const std::array<FadFadType,PHILIP_DIM+2> &conservative_soln) const;
1634 template dealii::Tensor<1,PHILIP_DIM,RadFadType> Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_velocities< RadFadType >(const std::array<RadFadType,PHILIP_DIM+2> &conservative_soln) const;
1635 // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
1636 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, double >::compute_velocities< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1637 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, RadType >::compute_velocities< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1638 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, FadFadType >::compute_velocities< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1639 template dealii::Tensor<1,PHILIP_DIM,FadType > Euler < PHILIP_DIM, PHILIP_DIM+2, RadFadType >::compute_velocities< FadType >(const std::array<FadType, PHILIP_DIM+2> &conservative_soln) const;
1640 //==============================================================================
1641 
1642 } // Physics namespace
1643 } // PHiLiP namespace
1644 
PHiLiP::Parameters::LimiterParam::LimiterType
LimiterType
Limiter type to be applied on the solution.
Definition: parameters_limiter.h:16
PHiLiP::FadFadType
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
Definition: ADTypes.hpp:12
PHiLiP::Physics::Euler::Euler
Euler(const Parameters::AllParameters *const parameters_input, const double ref_length, const double gamma_gas, const double mach_inf, const double angle_of_attack, const double side_slip_angle, std::shared_ptr< ManufacturedSolutionFunction< dim, real > > manufactured_solution_function=nullptr, const two_point_num_flux_enum two_point_num_flux_type=two_point_num_flux_enum::KG, const bool has_nonzero_diffusion=false, const bool has_nonzero_physical_source=false)
Constructor.
Definition: euler.cpp:13
PHiLiP::FadType
Sacado::Fad::DFad< double > FadType
Sacado AD type for first derivatives.
Definition: ADTypes.hpp:11
PHiLiP::Physics::PhysicsBase
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
Definition: physics.h:34
PHiLiP::ManufacturedSolutionFunction
Manufactured solution used for grid studies to check convergence orders.
Definition: manufactured_solution.h:22
PHiLiP::RadType
codi_JacobianComputationType RadType
CoDiPaco reverse-AD type for first derivatives.
Definition: ADTypes.hpp:27
PHiLiP
Files for the baseline physics.
Definition: ADTypes.hpp:10
PHiLiP::Parameters::AllParameters
Main parameter class that contains the various other sub-parameter classes.
Definition: all_parameters.h:33
PHiLiP::Physics::Euler
Euler equations. Derived from PhysicsBase.
Definition: euler.h:78
PHiLiP::RadFadType
codi_HessianComputationType RadFadType
Nested reverse-forward mode type for Jacobian and Hessian computation using TapeHelper.
Definition: ADTypes.hpp:28