12 template <
int dim,
int nstate,
typename real>
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,
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)
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)
40 static_assert(nstate==dim+2,
"Physics::Euler() should be created with nstate=dim+2");
43 temperature_inf = gam*pressure_inf/density_inf * mach_inf_sqr;
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;
52 velocities_inf[0] = 1.0;
54 velocities_inf[0] = cos(angle_of_attack);
55 velocities_inf[1] = sin(angle_of_attack);
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);
62 assert(std::abs(velocities_inf.norm() - 1.0) < 1e-14);
64 double velocity_inf_sqr = 1.0;
65 dynamic_pressure_inf = 0.5 * density_inf * velocity_inf_sqr;
68 template <
int dim,
int nstate,
typename real>
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 )
const 76 return source_term(pos,conservative_soln,current_time);
79 template <
int dim,
int nstate,
typename real>
82 const dealii::Point<dim,real> &pos,
83 const std::array<real,nstate> &,
86 std::array<real,nstate> source_term = convective_source_term(pos);
90 template <
int dim,
int nstate,
typename real>
93 const dealii::Point<dim,real> &pos)
const 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);
99 assert(manufactured_solution[s] > 0);
102 return manufactured_solution;
105 template <
int dim,
int nstate,
typename real>
108 const dealii::Point<dim,real> &pos)
const 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];
118 return manufactured_solution_gradient;
121 template <
int dim,
int nstate,
typename real>
124 const dealii::Point<dim,real> &pos)
const 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);
129 dealii::Tensor<1,nstate,real> convective_flux_divergence;
130 for (
int d=0;d<dim;d++) {
131 dealii::Tensor<1,dim,real> normal;
133 const dealii::Tensor<2,nstate,real> jacobian = convective_flux_directional_jacobian(manufactured_solution, normal);
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];
141 convective_flux_divergence[sr] += jac_grad_row;
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];
149 return convective_source_term;
152 template <
int dim,
int nstate,
typename real>
153 template<
typename real2>
156 bool qty_is_positive;
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) {
162 qty = this->
template handle_non_physical_result<real2>(qty_name +
" is negative.");
163 qty_is_positive =
false;
166 qty_is_positive =
true;
169 qty_is_positive =
true;
172 return qty_is_positive;
175 template <
int dim,
int nstate,
typename real>
176 template<
typename real2>
180 std::array<real2, nstate> primitive_soln;
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);
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];
192 primitive_soln[nstate-1] = pressure;
194 return primitive_soln;
197 template <
int dim,
int nstate,
typename real>
202 const real density = primitive_soln[0];
203 const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive<real>(primitive_soln);
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];
210 conservative_soln[nstate-1] = compute_total_energy(primitive_soln);
212 return conservative_soln;
222 template <
int dim,
int nstate,
typename real>
223 template<
typename real2>
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; }
233 template <
int dim,
int nstate,
typename real>
234 template <
typename real2>
239 for (
int d=0; d<dim; d++) {
240 vel2 = vel2 + velocities[d]*velocities[d];
246 template <
int dim,
int nstate,
typename real>
247 template<
typename real2>
251 dealii::Tensor<1,dim,real2> velocities;
252 for (
int d=0; d<dim; d++) { velocities[d] = primitive_soln[1+d]; }
256 template <
int dim,
int nstate,
typename real>
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;
266 template <
int dim,
int nstate,
typename real>
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;
277 template <
int dim,
int nstate,
typename real>
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;
286 template <
int dim,
int nstate,
typename real>
290 real density = conservative_soln[0];
291 const real pressure = compute_pressure<real>(conservative_soln);
292 return compute_entropy_measure(density, pressure);
295 template <
int dim,
int nstate,
typename real>
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;
307 template <
int dim,
int nstate,
typename real>
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;
317 template <
int dim,
int nstate,
typename real>
321 const real pressure = compute_pressure<real>(conservative_soln);
322 const real density = conservative_soln[0];
324 const real entropy = compute_entropy<real>(density, pressure);
326 const real numerical_entropy_function = - density * entropy;
328 return numerical_entropy_function;
331 template <
int dim,
int nstate,
typename real>
332 template<
typename real2>
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);
342 template <
int dim,
int nstate,
typename real>
346 const real density = gam*mach_inf_sqr*(pressure/temperature);
350 template <
int dim,
int nstate,
typename real>
354 const real temperature = gam*mach_inf_sqr*(pressure/density);
358 template <
int dim,
int nstate,
typename real>
362 const real pressure = density*temperature/(gam*mach_inf_sqr);
366 template <
int dim,
int nstate,
typename real>
367 template<
typename real2>
371 const real2 density = conservative_soln[0];
373 const real2 tot_energy = conservative_soln[nstate-1];
375 const dealii::Tensor<1,dim,real2> vel = compute_velocities<real2>(conservative_soln);
377 const real2 vel2 = compute_velocity_squared<real2>(vel);
378 real2 pressure = gamm1*(tot_energy - 0.5*density*vel2);
380 check_positive_quantity<real2>(pressure,
"pressure");
385 template <
int dim,
int nstate,
typename real>
386 template<
typename real2>
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);
400 return (real2)this->BIG_NUMBER;
405 template <
int dim,
int nstate,
typename real>
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);
416 template <
int dim,
int nstate,
typename real>
421 const real sound = sqrt(pressure*gam/density);
425 template <
int dim,
int nstate,
typename real>
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;
438 template <
int dim,
int nstate,
typename real>
441 const std::array<real,nstate> &conservative_soln2)
const 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);
454 return conv_num_split_flux;
457 template <
int dim,
int nstate,
typename real>
460 const std::array<real,nstate> &conservative_soln2)
const 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);
468 for (
int flux_dim = 0; flux_dim < dim; ++flux_dim)
471 conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];
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];
476 conv_num_split_flux[1+flux_dim][flux_dim] += mean_pressure;
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];
481 return conv_num_split_flux;
484 template <
int dim,
int nstate,
typename real>
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];
494 ismail_roe_parameter_vector[nstate-1] = sqrt(primitive_soln[0]*primitive_soln[nstate-1]);
496 return ismail_roe_parameter_vector;
499 template <
int dim,
int nstate,
typename real>
505 const real zeta = val1/val2;
506 const real f = (zeta-1.0)/(zeta+1.0);
510 if(u<1.0e-2){ F = 1.0 + u/3.0 + u*u/5.0 + u*u*u/7.0; }
512 if constexpr(std::is_same<real,double>::value) F = std::log(zeta)/2.0/f;
515 const real log_mean_val = (val1+val2)/(2.0*F);
520 template <
int dim,
int nstate,
typename real>
523 const std::array<real,nstate> &conservative_soln2)
const 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));
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]);
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]);
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];
549 const real mean_pressure = avg_parameter_vector[nstate-1]/avg_parameter_vector[0];
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);
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];
560 mean_enthalpy += 0.5*mean_velocities_sqr_sum;
563 std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
564 for (
int flux_dim = 0; flux_dim < dim; ++flux_dim)
567 conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];
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];
572 conv_num_split_flux[1+flux_dim][flux_dim] += mean_pressure;
574 conv_num_split_flux[nstate-1][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_enthalpy;
577 return conv_num_split_flux;
580 template <
int dim,
int nstate,
typename real>
583 const std::array<real,nstate> &conservative_soln2)
const 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);
591 const real beta1 = conservative_soln1[0]/(2.0*pressure1);
592 const real beta2 = conservative_soln2[0]/(2.0*pressure2);
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);
598 const real pressure_hat = 0.5*(conservative_soln1[0] + conservative_soln2[0])/(2.0*0.5*(beta1+beta2));
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]));
607 real enthalpy_hat = 1.0/(2.0*beta_log*gamm1) + vel_square_avg + pressure_hat/rho_log;
609 for(
int idim=0; idim<dim; idim++){
610 enthalpy_hat -= 0.5*(0.5*(vel1[idim]*vel1[idim] + vel2[idim]*vel2[idim]));
613 for(
int flux_dim=0; flux_dim<dim; flux_dim++){
615 conv_num_split_flux[0][flux_dim] = rho_log * vel_avg[flux_dim];
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];
620 conv_num_split_flux[1+flux_dim][flux_dim] += pressure_hat;
623 conv_num_split_flux[nstate-1][flux_dim] = rho_log * vel_avg[flux_dim] * enthalpy_hat;
626 return conv_num_split_flux;
629 template <
int dim,
int nstate,
typename real>
632 const std::array<real,nstate> &conservative_soln2)
const 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);
640 const real beta1 = conservative_soln1[0]/(pressure1);
641 const real beta2 = conservative_soln2[0]/(pressure2);
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);
647 const real pressure_hat = 0.5*(pressure1+pressure2);
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]));
656 real enthalpy_hat = 1.0/(beta_log*gamm1) + vel_square_avg + 2.0*pressure_hat/rho_log;
658 for(
int idim=0; idim<dim; idim++){
659 enthalpy_hat -= 0.5*(0.5*(vel1[idim]*vel1[idim] + vel2[idim]*vel2[idim]));
662 for(
int flux_dim=0; flux_dim<dim; flux_dim++){
664 conv_num_split_flux[0][flux_dim] = rho_log * vel_avg[flux_dim];
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];
669 conv_num_split_flux[1+flux_dim][flux_dim] += pressure_hat;
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]));
676 return conv_num_split_flux;
680 template <
int dim,
int nstate,
typename real>
683 const std::array<real,nstate> &conservative_soln)
const 685 std::array<real,nstate> entropy_var;
686 const real density = conservative_soln[0];
687 const real pressure = compute_pressure<real>(conservative_soln);
689 const real entropy = compute_entropy<real>(density, pressure);
691 const real rho_theta = pressure / gamm1;
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;
697 entropy_var[nstate-1] = - density / rho_theta;
702 template <
int dim,
int nstate,
typename real>
705 const std::array<real,nstate> &entropy_var)
const 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];
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);
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];
722 conservative_var[nstate-1] = rho_theta * (1.0 - 0.5 * entropy_var_vel_squared / entropy_var[nstate-1]);
723 return conservative_var;
726 template <
int dim,
int nstate,
typename real>
729 const std::array<real,nstate> &conservative_soln)
const 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);
735 kin_energy_var[0] = - 0.5 * vel2;
736 for(
int idim=0; idim<dim; idim++){
737 kin_energy_var[idim+1] = vel[idim];
739 kin_energy_var[nstate-1] = 0;
741 return kin_energy_var;
744 template <
int dim,
int nstate,
typename real>
747 const std::array<real,nstate> &conservative_soln2)
const 749 return (conservative_soln1[0] + conservative_soln2[0])/2.;
752 template <
int dim,
int nstate,
typename real>
755 const std::array<real,nstate> &conservative_soln2)
const 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.;
762 template <
int dim,
int nstate,
typename real>
765 const std::array<real,nstate> &conservative_soln2)
const 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]);
776 template <
int dim,
int nstate,
typename real>
779 const std::array<real,nstate> &conservative_soln2)
const 781 return ((conservative_soln1[nstate-1]/conservative_soln1[0]) + (conservative_soln2[nstate-1]/conservative_soln2[0]))/2.;
785 template <
int dim,
int nstate,
typename real>
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;
796 for (
int flux_dim=0; flux_dim<dim; ++flux_dim) {
798 conv_flux[0][flux_dim] = conservative_soln[1+flux_dim];
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];
803 conv_flux[1+flux_dim][flux_dim] += pressure;
805 conv_flux[nstate-1][flux_dim] = density*vel[flux_dim]*specific_total_enthalpy;
810 template <
int dim,
int nstate,
typename real>
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];
822 const real total_energy = conservative_soln[nstate-1];
823 const real specific_total_enthalpy = (total_energy + pressure) / density;
825 const real rhoV = density*normal_vel;
827 conv_normal_flux[0] = rhoV;
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;
833 conv_normal_flux[nstate-1] = rhoV*specific_total_enthalpy;
834 return conv_normal_flux;
837 template <
int dim,
int nstate,
typename real>
840 const std::array<real,nstate> &conservative_soln,
841 const dealii::Tensor<1,dim,real> &normal)
const 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]; }
850 const real vel2 = compute_velocity_squared<real>(vel);
851 const real phi = 0.5*gamm1 * vel2;
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;
860 dealii::Tensor<2,nstate,real> jacobian;
861 for (
int d=0; d<dim; ++d) {
862 jacobian[0][1+d] = normal[d];
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];
870 jacobian[1+row_dim][1+col_dim] = normal[col_dim]*vel[row_dim] - a2*normal[row_dim]*vel[col_dim];
873 jacobian[1+row_dim][nstate-1] = normal[row_dim]*a2;
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;
879 jacobian[nstate-1][nstate-1] = gam*vel_normal;
884 template <
int dim,
int nstate,
typename real>
887 const std::array<real,nstate> &conservative_soln,
888 const dealii::Tensor<1,dim,real> &normal)
const 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++) {
904 template <
int dim,
int nstate,
typename real>
908 const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
910 const real sound = compute_sound (conservative_soln);
912 real vel2 = compute_velocity_squared<real>(vel);
914 const real max_eig = sqrt(vel2) + sound;
919 template <
int dim,
int nstate,
typename real>
922 const std::array<real,nstate> &conservative_soln,
923 const dealii::Tensor<1,dim,real> &normal)
const 925 const dealii::Tensor<1,dim,real> vel = compute_velocities<real>(conservative_soln);
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;
932 return max_normal_eig;
935 template <
int dim,
int nstate,
typename real>
942 template <
int dim,
int nstate,
typename real>
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 )
const 949 return dissipative_flux(conservative_soln, solution_gradient);
952 template <
int dim,
int nstate,
typename real>
955 const std::array<real,nstate> &,
956 const std::array<dealii::Tensor<1,dim,real>,nstate> &)
const 958 std::array<dealii::Tensor<1,dim,real>,nstate> diss_flux;
960 for (
int i=0; i<nstate; i++) {
966 template <
int dim,
int nstate,
typename real>
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 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;
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);
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] );
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];
993 const real out_riemann_invariant = vel_int_dot_normal + 2.0/gamm1*sound_int,
994 inc_riemann_invariant = vel_ext_dot_normal - 2.0/gamm1*sound_ext;
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);
999 std::array<real,nstate> primitive_bc;
1000 if (abs(normal_velocity_bc) >= abs(sound_bc)) {
1001 if (normal_velocity_bc < 0.0) {
1002 primitive_bc = primitive_ext;
1004 primitive_bc = primitive_int;
1009 dealii::Tensor<1,dim,real> velocities_bc;
1012 dealii::Tensor<1,dim,real> velocities_tangential;
1013 if (normal_velocity_bc < 0.0) {
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];
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];
1026 for (
int d=0; d<dim; ++d) {
1027 velocities_bc[d] = velocities_tangential[d] + normal_velocity_bc*normal_int[d];
1030 pressure_bc = 1.0/gam * sound_bc * sound_bc * density_bc;
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;
1037 soln_bc = convert_primitive_to_conservative(primitive_bc);
1040 template <
int dim,
int nstate,
typename real>
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 1054 const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive<real>(soln_int);
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];
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);
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];
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];
1074 for (
int d=0; d<dim; ++d) {
1075 primitive_boundary_values[1+d] = velocities_bc[d];
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];
1083 for (
int istate=0; istate<nstate; ++istate) {
1084 soln_grad_bc[istate] = -soln_grad_int[istate];
1088 template <
int dim,
int nstate,
typename real>
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 1098 boundary_slip_wall(normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1101 template <
int dim,
int nstate,
typename real>
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 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);
1118 std::array<real,nstate> primitive_boundary_values = convert_conservative_to_primitive<real>(conservative_boundary_values);
1119 for (
int istate=0; istate<nstate; ++istate) {
1121 std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(conservative_boundary_values, normal_int);
1122 const bool inflow = (characteristic_dot_n[istate] <= 0.);
1126 soln_bc[istate] = conservative_boundary_values[istate];
1127 soln_grad_bc[istate] = soln_grad_int[istate];
1134 const std::array<real,nstate> modified_conservative_boundary_values = convert_primitive_to_conservative(primitive_boundary_values);
1135 (void) modified_conservative_boundary_values;
1137 soln_bc[istate] = conservative_boundary_values[istate];
1141 soln_bc[istate] = soln_int[istate];
1148 soln_grad_bc[istate] = soln_grad_int[istate];
1154 soln_bc[istate] = conservative_boundary_values[istate];
1158 template <
int dim,
int nstate,
typename real>
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 1169 const real mach_int = compute_mach_number(soln_int);
1170 if (mach_int > 1.0) {
1172 for (
int istate=0; istate<nstate; ++istate) {
1173 soln_bc[istate] = soln_int[istate];
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];
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);
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;
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];
1198 template <
int dim,
int nstate,
typename real>
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 1210 const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive<real>(soln_int);
1212 const dealii::Tensor<1,dim,real> normal = -normal_int;
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];
1219 real normal_vel_i = 0.0;
1220 for (
int d=0; d<dim; ++d) {
1221 normal_vel_i += velocities_i[d]*normal[d];
1223 const real sound_i = compute_sound(soln_int);
1231 if(mach_inf < 1.0) {
1238 const real riemann_pos = normal_vel_i + 2.0*sound_i/gamm1;
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;
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);
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;
1252 const real sound_bc = std::max(sound_bc1, sound_bc2);
1255 const real velocity_magnitude_bc = riemann_pos - 2.0*sound_bc/gamm1;
1256 const real mach_bc = velocity_magnitude_bc/sound_bc;
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);
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];
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);
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;
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]; }
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];
1305 template <
int dim,
int nstate,
typename real>
1308 std::array<real,nstate> &soln_bc)
const 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]; }
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];
1322 template <
int dim,
int nstate,
typename real>
1325 const int boundary_type,
1326 const dealii::Point<dim, real> &pos,
1327 const dealii::Tensor<1,dim,real> &normal_int,
1328 const std::array<real,nstate> &soln_int,
1329 const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
1330 std::array<real,nstate> &soln_bc,
1331 std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc)
const 1334 const real total_inlet_pressure = pressure_inf*pow(1.0+0.5*gamm1*mach_inf_sqr, gam/gamm1);
1335 const real total_inlet_temperature = temperature_inf*pow(total_inlet_pressure/pressure_inf, gamm1/gam);
1337 if (boundary_type == 1000) {
1339 boundary_manufactured_solution (pos, normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1341 else if (boundary_type == 1001) {
1343 boundary_wall (normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1345 else if (boundary_type == 1002) {
1347 const real back_pressure = 0.99;
1348 boundary_pressure_outflow (total_inlet_pressure, back_pressure, soln_int, soln_bc);
1350 else if (boundary_type == 1003) {
1352 boundary_inflow (total_inlet_pressure, total_inlet_temperature, normal_int, soln_int, soln_bc);
1354 else if (boundary_type == 1004) {
1356 boundary_riemann (normal_int, soln_int, soln_bc);
1358 else if (boundary_type == 1005) {
1360 boundary_farfield(soln_bc);
1362 else if (boundary_type == 1006) {
1364 boundary_slip_wall (normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
1367 this->pcout <<
"Invalid boundary_type: " << boundary_type << std::endl;
1372 template <
int dim,
int nstate,
typename real>
1374 const dealii::Vector<double> &uh,
1375 const std::vector<dealii::Tensor<1,dim> > &duh,
1376 const std::vector<dealii::Tensor<2,dim> > &dduh,
1377 const dealii::Tensor<1,dim> &normals,
1378 const dealii::Point<dim> &evaluation_points)
const 1380 std::vector<std::string> names = post_get_names ();
1382 unsigned int current_data_index = computed_quantities.size() - 1;
1383 computed_quantities.grow_or_shrink(names.size());
1384 if constexpr (std::is_same<real,double>::value) {
1386 std::array<double, nstate> conservative_soln;
1387 for (
unsigned int s=0; s<nstate; ++s) {
1388 conservative_soln[s] = uh(s);
1390 const std::array<double, nstate> primitive_soln = convert_conservative_to_primitive<real>(conservative_soln);
1394 computed_quantities(++current_data_index) = primitive_soln[0];
1396 for (
unsigned int d=0; d<dim; ++d) {
1397 computed_quantities(++current_data_index) = primitive_soln[1+d];
1400 for (
unsigned int d=0; d<dim; ++d) {
1401 computed_quantities(++current_data_index) = conservative_soln[1+d];
1404 computed_quantities(++current_data_index) = conservative_soln[nstate-1];
1406 computed_quantities(++current_data_index) = primitive_soln[nstate-1];
1408 computed_quantities(++current_data_index) = (primitive_soln[nstate-1] - pressure_inf) / dynamic_pressure_inf;
1410 computed_quantities(++current_data_index) = compute_temperature<real>(primitive_soln);
1412 computed_quantities(++current_data_index) = compute_entropy_measure(conservative_soln) - entropy_inf;
1414 computed_quantities(++current_data_index) = compute_mach_number(conservative_soln);
1417 if (computed_quantities.size()-1 != current_data_index) {
1418 this->pcout <<
" Did not assign a value to all the data. Missing " << computed_quantities.size() - current_data_index <<
" variables." 1419 <<
" If you added a new output variable, make sure the names and DataComponentInterpretation match the above. " 1423 return computed_quantities;
1426 template <
int dim,
int nstate,
typename real>
1430 namespace DCI = dealii::DataComponentInterpretation;
1432 interpretation.push_back (DCI::component_is_scalar);
1433 for (
unsigned int d=0; d<dim; ++d) {
1434 interpretation.push_back (DCI::component_is_part_of_vector);
1436 for (
unsigned int d=0; d<dim; ++d) {
1437 interpretation.push_back (DCI::component_is_part_of_vector);
1439 interpretation.push_back (DCI::component_is_scalar);
1440 interpretation.push_back (DCI::component_is_scalar);
1441 interpretation.push_back (DCI::component_is_scalar);
1442 interpretation.push_back (DCI::component_is_scalar);
1443 interpretation.push_back (DCI::component_is_scalar);
1444 interpretation.push_back (DCI::component_is_scalar);
1446 std::vector<std::string> names = post_get_names();
1447 if (names.size() != interpretation.size()) {
1448 this->pcout <<
"Number of DataComponentInterpretation is not the same as number of names for output file" << std::endl;
1450 return interpretation;
1454 template <
int dim,
int nstate,
typename real>
1459 names.push_back (
"density");
1460 for (
unsigned int d=0; d<dim; ++d) {
1461 names.push_back (
"velocity");
1463 for (
unsigned int d=0; d<dim; ++d) {
1464 names.push_back (
"momentum");
1466 names.push_back (
"energy");
1467 names.push_back (
"pressure");
1468 names.push_back (
"pressure_coeffcient");
1469 names.push_back (
"temperature");
1471 names.push_back (
"entropy_generation");
1472 names.push_back (
"mach_number");
1476 template <
int dim,
int nstate,
typename real>
1481 return dealii::update_values
1482 | dealii::update_quadrature_points
LimiterType
Limiter type to be applied on the solution.
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
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.
Sacado::Fad::DFad< double > FadType
Sacado AD type for first derivatives.
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
Manufactured solution used for grid studies to check convergence orders.
codi_JacobianComputationType RadType
CoDiPaco reverse-AD type for first derivatives.
Files for the baseline physics.
Main parameter class that contains the various other sub-parameter classes.
Euler equations. Derived from PhysicsBase.
codi_HessianComputationType RadFadType
Nested reverse-forward mode type for Jacobian and Hessian computation using TapeHelper.