3 #include "convective_numerical_flux.hpp" 6 namespace NumericalFlux {
8 using AllParam = Parameters::AllParameters;
11 template<
int nstate,
typename real_tensor>
12 std::array<real_tensor, nstate> array_average(
13 const std::array<real_tensor, nstate> &array1,
14 const std::array<real_tensor, nstate> &array2)
16 std::array<real_tensor,nstate> array_average;
17 for (
int s=0; s<nstate; s++) {
18 array_average[s] = 0.5*(array1[s] + array2[s]);
23 template <
int dim,
int nstate,
typename real>
27 : baseline(
std::move(baseline_input))
28 , riemann_solver_dissipation(
std::move(riemann_solver_dissipation_input))
31 template<
int dim,
int nstate,
typename real>
34 const std::array<real, nstate> &soln_int,
35 const std::array<real, nstate> &soln_ext,
36 const dealii::Tensor<1,dim,real> &normal_int)
const 39 const std::array<real, nstate> baseline_flux_dot_n
40 = this->
baseline->evaluate_flux(soln_int, soln_ext, normal_int);
43 const std::array<real, nstate> riemann_solver_dissipation_dot_n
47 std::array<real, nstate> numerical_flux_dot_n;
48 for (
int s=0; s<nstate; s++) {
49 numerical_flux_dot_n[s] = baseline_flux_dot_n[s] + riemann_solver_dissipation_dot_n[s];
51 return numerical_flux_dot_n;
54 template <
int dim,
int nstate,
typename real>
62 template <
int dim,
int nstate,
typename real>
70 template <
int dim,
int nstate,
typename real>
78 template <
int dim,
int nstate,
typename real>
86 template <
int dim,
int nstate,
typename real>
94 template <
int dim,
int nstate,
typename real>
102 template <
int dim,
int nstate,
typename real>
110 template <
int dim,
int nstate,
typename real>
118 template <
int dim,
int nstate,
typename real>
120 const std::array<real, nstate> &soln_int,
121 const std::array<real, nstate> &soln_ext,
122 const dealii::Tensor<1,dim,real> &normal_int)
const 124 using RealArrayVector = std::array<dealii::Tensor<1,dim,real>,nstate>;
125 RealArrayVector conv_phys_flux_int;
126 RealArrayVector conv_phys_flux_ext;
128 conv_phys_flux_int = pde_physics->convective_flux (soln_int);
129 conv_phys_flux_ext = pde_physics->convective_flux (soln_ext);
131 RealArrayVector flux_avg;
132 for (
int s=0; s<nstate; s++) {
134 for (
int d=0; d<dim; ++d) {
135 flux_avg[s][d] = 0.5*(conv_phys_flux_int[s][d] + conv_phys_flux_ext[s][d]);
139 std::array<real, nstate> numerical_flux_dot_n;
140 for (
int s=0; s<nstate; s++) {
141 real flux_dot_n = 0.0;
142 for (
int d=0; d<dim; ++d) {
143 flux_dot_n += flux_avg[s][d]*normal_int[d];
145 numerical_flux_dot_n[s] = flux_dot_n;
147 return numerical_flux_dot_n;
150 template <
int dim,
int nstate,
typename real>
152 const std::array<real, nstate> &soln_int,
153 const std::array<real, nstate> &soln_ext,
154 const dealii::Tensor<1,dim,real> &normal_int)
const 156 using RealArrayVector = std::array<dealii::Tensor<1,dim,real>,nstate>;
157 RealArrayVector conv_phys_split_flux;
159 conv_phys_split_flux = pde_physics->convective_numerical_split_flux (soln_int,soln_ext);
162 std::array<real, nstate> numerical_flux_dot_n;
163 for (
int s=0; s<nstate; s++) {
164 real flux_dot_n = 0.0;
165 for (
int d=0; d<dim; ++d) {
166 flux_dot_n += conv_phys_split_flux[s][d] * normal_int[d];
168 numerical_flux_dot_n[s] = flux_dot_n;
170 return numerical_flux_dot_n;
173 template<
int dim,
int nstate,
typename real>
176 const std::array<real, nstate> &,
177 const std::array<real, nstate> &,
178 const dealii::Tensor<1,dim,real> &)
const 181 std::array<real, nstate> numerical_flux_dot_n;
182 numerical_flux_dot_n.fill(0.0);
183 return numerical_flux_dot_n;
186 template<
int dim,
int nstate,
typename real>
189 const std::array<real, nstate> &soln_int,
190 const std::array<real, nstate> &soln_ext,
191 const dealii::Tensor<1,dim,real> &normal_int)
const 193 const real conv_max_eig_int = pde_physics->max_convective_normal_eigenvalue(soln_int,normal_int);
194 const real conv_max_eig_ext = pde_physics->max_convective_normal_eigenvalue(soln_ext,normal_int);
198 if (conv_max_eig_int > conv_max_eig_ext) {
199 conv_max_eig = conv_max_eig_int;
201 conv_max_eig = conv_max_eig_ext;
205 std::array<real, nstate> numerical_flux_dot_n;
206 for (
int s=0; s<nstate; s++) {
207 numerical_flux_dot_n[s] = - 0.5 * conv_max_eig * (soln_ext[s]-soln_int[s]);
210 return numerical_flux_dot_n;
213 template <
int dim,
int nstate,
typename real>
216 const std::array<real, 3> &eig_L,
217 const std::array<real, 3> &eig_R,
218 std::array<real, 3> &eig_RoeAvg,
224 for(
int e=0;e<3;e++) {
225 const real eps = std::max(abs(eig_RoeAvg[e] - eig_L[e]), abs(eig_R[e] - eig_RoeAvg[e]));
226 if(eig_RoeAvg[e] < eps) {
227 eig_RoeAvg[e] = 0.5*(eig_RoeAvg[e] * eig_RoeAvg[e]/eps + eps);
232 template <
int dim,
int nstate,
typename real>
235 const std::array<real, nstate> &,
236 const std::array<real, nstate> &,
237 const std::array<real, 3> &,
238 const std::array<real, 3> &,
240 dealii::Tensor<1,dim,real> &
246 template <
int dim,
int nstate,
typename real>
249 const std::array<real, 3> &eig_L,
250 const std::array<real, 3> &eig_R,
252 int &ssw_RIGHT)
const 257 ssw_LEFT=0; ssw_RIGHT=0;
260 if((eig_L[0]>0.0 && eig_R[0]<0.0) || (eig_L[2]>0.0 && eig_R[2]<0.0)) {
265 if((eig_R[0]>0.0 && eig_L[0]<0.0) || (eig_R[2]>0.0 && eig_L[2]<0.0)) {
270 template <
int dim,
int nstate,
typename real>
273 const std::array<real, 3> &eig_L,
274 const std::array<real, 3> &eig_R,
275 std::array<real, 3> &eig_RoeAvg,
276 const real vel2_ravg,
277 const real sound_ravg)
const 281 for(
int e=0;e<3;e++) {
284 const real deig = std::max(static_cast<real>(eig_R[e] - eig_L[e]), static_cast<real>(0.0));
285 if(eig_RoeAvg[e] < 2.0*deig) {
286 eig_RoeAvg[e] = 0.25*(eig_RoeAvg[e] * eig_RoeAvg[e]/deig) + deig;
294 evaluate_shock_indicator(eig_L,eig_R,ssw_L,ssw_R);
295 if(ssw_L!=0 || ssw_R!=0) {
296 eig_RoeAvg[1] = std::max(sound_ravg, static_cast<real>(sqrt(vel2_ravg)));
300 template <
int dim,
int nstate,
typename real>
303 const std::array<real, nstate> &soln_int,
304 const std::array<real, nstate> &soln_ext,
305 const std::array<real, 3> &eig_L,
306 const std::array<real, 3> &eig_R,
308 dealii::Tensor<1,dim,real> &dV_tangent)
const 310 const real mach_number_L = this->euler_physics->compute_mach_number(soln_int);
311 const real mach_number_R = this->euler_physics->compute_mach_number(soln_ext);
315 const real blending_factor = std::min(static_cast<real>(1.0), std::max(mach_number_L,mach_number_R));
318 evaluate_shock_indicator(eig_L,eig_R,ssw_L,ssw_R);
319 if(ssw_L==0 && ssw_R==0)
321 dV_normal *= blending_factor;
322 for (
int d=0;d<dim;d++)
324 dV_tangent[d] *= blending_factor;
329 template <
int dim,
int nstate,
typename real>
332 const std::array<real, nstate> &soln_int,
333 const std::array<real, nstate> &soln_ext,
334 const dealii::Tensor<1,dim,real> &normal_int)
const 344 const std::array<real,nstate> prim_soln_int = euler_physics->convert_conservative_to_primitive(soln_int);
345 const std::array<real,nstate> prim_soln_ext = euler_physics->convert_conservative_to_primitive(soln_ext);
347 const real density_L = prim_soln_int[0];
348 const dealii::Tensor< 1,dim,real > velocities_L = euler_physics->extract_velocities_from_primitive(prim_soln_int);
349 const real pressure_L = prim_soln_int[nstate-1];
352 real normal_vel_L = 0.0;
353 for (
int d=0; d<dim; ++d) {
354 normal_vel_L+= velocities_L[d]*normal_int[d];
356 const real specific_enthalpy_L = euler_physics->compute_specific_enthalpy(soln_int, pressure_L);
359 const real density_R = prim_soln_ext[0];
360 const dealii::Tensor< 1,dim,real > velocities_R = euler_physics->extract_velocities_from_primitive(prim_soln_ext);
361 const real pressure_R = prim_soln_ext[nstate-1];
364 real normal_vel_R = 0.0;
365 for (
int d=0; d<dim; ++d) {
366 normal_vel_R+= velocities_R[d]*normal_int[d];
368 const real specific_enthalpy_R = euler_physics->compute_specific_enthalpy(soln_ext, pressure_R);
371 const real r = sqrt(density_R/density_L);
372 const real rp1 = r+1.0;
374 const real density_ravg = r*density_L;
376 dealii::Tensor< 1,dim,real > velocities_ravg;
377 for (
int d=0; d<dim; ++d) {
378 velocities_ravg[d] = (r*velocities_R[d] + velocities_L[d]) / rp1;
380 const real specific_total_enthalpy_ravg = (r*specific_enthalpy_R + specific_enthalpy_L) / rp1;
382 const real vel2_ravg = euler_physics->compute_velocity_squared (velocities_ravg);
384 real normal_vel_ravg = 0.0;
385 for (
int d=0; d<dim; ++d) {
386 normal_vel_ravg += velocities_ravg[d]*normal_int[d];
389 const real sound2_ravg = euler_physics->gamm1*(specific_total_enthalpy_ravg-0.5*vel2_ravg);
390 real sound_ravg = 1e10;
391 if (sound2_ravg > 0.0) {
392 sound_ravg = sqrt(sound2_ravg);
396 std::array<real, 3> eig_ravg;
397 eig_ravg[0] = abs(normal_vel_ravg-sound_ravg);
398 eig_ravg[1] = abs(normal_vel_ravg);
399 eig_ravg[2] = abs(normal_vel_ravg+sound_ravg);
401 const real sound_L = euler_physics->compute_sound(density_L, pressure_L);
402 std::array<real, 3> eig_L;
403 eig_L[0] = abs(normal_vel_L-sound_L);
404 eig_L[1] = abs(normal_vel_L);
405 eig_L[2] = abs(normal_vel_L+sound_L);
407 const real sound_R = euler_physics->compute_sound(density_R, pressure_R);
408 std::array<real, 3> eig_R;
409 eig_R[0] = abs(normal_vel_R-sound_R);
410 eig_R[1] = abs(normal_vel_R);
411 eig_R[2] = abs(normal_vel_R+sound_R);
414 const real dp = pressure_R - pressure_L;
415 const real drho = density_R - density_L;
418 real dVn = normal_vel_R-normal_vel_L;
421 dealii::Tensor<1,dim,real> dVt;
422 for (
int d=0;d<dim;d++) {
423 dVt[d] = (velocities_R[d] - velocities_L[d]) - dVn*normal_int[d];
427 evaluate_entropy_fix (eig_L, eig_R, eig_ravg, vel2_ravg, sound_ravg);
430 evaluate_additional_modifications (soln_int, soln_ext, eig_L, eig_R, dVn, dVt);
434 coeff[0] = eig_ravg[0]*(dp-density_ravg*sound_ravg*dVn)/(2.0*sound2_ravg);
435 coeff[1] = eig_ravg[1]*(drho - dp/sound2_ravg);
436 coeff[2] = eig_ravg[1]*density_ravg;
437 coeff[3] = eig_ravg[2]*(dp+density_ravg*sound_ravg*dVn)/(2.0*sound2_ravg);
440 std::array<real,nstate> AdW;
443 AdW[0] = coeff[0] * 1.0;
444 for (
int d=0;d<dim;d++) {
445 AdW[1+d] = coeff[0] * (velocities_ravg[d] - sound_ravg * normal_int[d]);
447 AdW[nstate-1] = coeff[0] * (specific_total_enthalpy_ravg - sound_ravg*normal_vel_ravg);
450 AdW[0] += coeff[1] * 1.0;
451 for (
int d=0;d<dim;d++) {
452 AdW[1+d] += coeff[1] * velocities_ravg[d];
454 AdW[nstate-1] += coeff[1] * vel2_ravg * 0.5;
457 AdW[0] += coeff[2] * 0.0;
458 real dVt_dot_vel_ravg = 0.0;
459 for (
int d=0;d<dim;d++) {
460 AdW[1+d] += coeff[2]*dVt[d];
461 dVt_dot_vel_ravg += velocities_ravg[d]*dVt[d];
463 AdW[nstate-1] += coeff[2]*dVt_dot_vel_ravg;
466 AdW[0] += coeff[3] * 1.0;
467 for (
int d=0;d<dim;d++) {
468 AdW[1+d] += coeff[3] * (velocities_ravg[d] + sound_ravg * normal_int[d]);
470 AdW[nstate-1] += coeff[3] * (specific_total_enthalpy_ravg + sound_ravg*normal_vel_ravg);
472 std::array<real, nstate> numerical_flux_dot_n;
473 for (
int s=0; s<nstate; s++) {
474 numerical_flux_dot_n[s] = - 0.5 * AdW[s];
477 return numerical_flux_dot_n;
std::unique_ptr< BaselineNumericalFluxConvective< dim, nstate, real > > baseline
Baseline convective numerical flux object.
EntropyConservingWithRoeDissipation(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
std::unique_ptr< RiemannSolverDissipation< dim, nstate, real > > riemann_solver_dissipation
Upwind convective numerical flux object.
RoePike(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
Central numerical flux. Derived from BaselineNumericalFluxConvective.
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
void evaluate_entropy_fix(const std::array< real, 3 > &eig_L, const std::array< real, 3 > &eig_R, std::array< real, 3 > &eig_RoeAvg, const real vel2_ravg, const real sound_ravg) const
Entropy conserving numerical flux with Roe dissipation. Derived from NumericalFluxConvective.
std::array< real, nstate > evaluate_riemann_solver_dissipation(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const dealii::Tensor< 1, dim, real > &normal1) const
Returns zeros.
Files for the baseline physics.
EntropyConservingWithL2RoeDissipation(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
EntropyConservingWithLaxFriedrichsDissipation(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
Roe-Pike numerical flux. Derived from NumericalFluxConvective.
std::array< real, nstate > evaluate_riemann_solver_dissipation(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const dealii::Tensor< 1, dim, real > &normal1) const
Entropy conserving numerical flux with L2Roe dissipation. Derived from NumericalFluxConvective.
Base class of Roe (Roe-Pike) flux with entropy fix. Derived from RiemannSolverDissipation.
Base class of Riemann solver dissipation (i.e. upwind-term) for numerical flux associated with convec...
void evaluate_additional_modifications(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const std::array< real, 3 > &eig_L, const std::array< real, 3 > &eig_R, real &dV_normal, dealii::Tensor< 1, dim, real > &dV_tangent) const
Entropy conserving numerical flux. Derived from NumericalFluxConvective.
void evaluate_entropy_fix(const std::array< real, 3 > &eig_L, const std::array< real, 3 > &eig_R, std::array< real, 3 > &eig_RoeAvg, const real vel2_ravg, const real sound_ravg) const
NumericalFluxConvective(std::unique_ptr< BaselineNumericalFluxConvective< dim, nstate, real > > baseline_input, std::unique_ptr< RiemannSolverDissipation< dim, nstate, real > > riemann_solver_dissipation_input)
Constructor.
Base class of baseline numerical flux (without upwind term) associated with convection.
Base class of numerical flux associated with convection.
L2Roe numerical flux. Derived from NumericalFluxConvective.
LaxFriedrichs(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
Lax-Friedrichs Riemann solver dissipation. Derived from RiemannSolverDissipation. ...
RoePike flux with entropy fix. Derived from RoeBase.
std::array< real, nstate > evaluate_riemann_solver_dissipation(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const dealii::Tensor< 1, dim, real > &normal1) const
Entropy conserving numerical flux with Lax Friedrichs dissipation. Derived from NumericalFluxConvecti...
Zero Riemann solver dissipation. Derived from RiemannSolverDissipation.
void evaluate_additional_modifications(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const std::array< real, 3 > &eig_L, const std::array< real, 3 > &eig_R, real &dV_normal, dealii::Tensor< 1, dim, real > &dV_tangent) const
Empty function. No additional modifications for the Roe-Pike scheme.
EntropyConserving(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
Central numerical flux. Derived from NumericalFluxConvective.
std::array< real, nstate > evaluate_flux(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const dealii::Tensor< 1, dim, real > &normal1) const
Returns the convective numerical flux at an interface.
void evaluate_shock_indicator(const std::array< real, 3 > &eig_L, const std::array< real, 3 > &eig_R, int &ssw_LEFT, int &ssw_RIGHT) const
std::array< real, nstate > evaluate_flux(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const dealii::Tensor< 1, dim, real > &normal1) const
Returns the convective numerical flux at an interface.
Lax-Friedrichs numerical flux. Derived from NumericalFluxConvective.
std::array< real, nstate > evaluate_flux(const std::array< real, nstate > &soln_int, const std::array< real, nstate > &soln_ext, const dealii::Tensor< 1, dim, real > &normal1) const
Returns the convective numerical flux at an interface.
Central(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
L2Roe(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real >> physics_input)
Constructor.
Entropy Conserving Numerical Flux. Derived from BaselineNumericalFluxConvective.