diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index dca846cac1e..74d672ce7ae 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -1,3 +1,10 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + @doc raw""" CompressibleNavierStokesDiffusion1D(equations; mu, Pr, gradient_variables=GradientVariablesPrimitive()) @@ -77,7 +84,8 @@ w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p} This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, - E <: AbstractCompressibleEulerEquations{1}} <: + E <: AbstractCompressibleEulerEquations{1} + } <: AbstractCompressibleNavierStokesDiffusion{1, 3} # TODO: parabolic # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations @@ -109,7 +117,8 @@ function CompressibleNavierStokesDiffusion1D(equations::CompressibleEulerEquatio CompressibleNavierStokesDiffusion1D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, μ, Pr, kappa, - equations, gradient_variables) + equations, + gradient_variables) end # TODO: parabolic @@ -263,7 +272,8 @@ end u_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesPrimitive @@ -278,7 +288,8 @@ end u_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesPrimitive @@ -299,7 +310,8 @@ end u_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesPrimitive @@ -316,7 +328,8 @@ end u_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesPrimitive @@ -337,7 +350,8 @@ end w_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesEntropy @@ -354,7 +368,8 @@ end w_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesEntropy @@ -374,7 +389,8 @@ end w_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesEntropy @@ -394,10 +410,12 @@ end w_inner, orientation::Integer, direction, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion1D{ GradientVariablesEntropy }) return SVector(flux_inner[1], flux_inner[2], flux_inner[3]) end +end # @muladd diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index f762fe5d5ee..b10ffa3b9d3 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -1,3 +1,10 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + @doc raw""" CompressibleNavierStokesDiffusion2D(equations; mu, Pr, gradient_variables=GradientVariablesPrimitive()) @@ -77,7 +84,8 @@ w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, - E <: AbstractCompressibleEulerEquations{2}} <: + E <: AbstractCompressibleEulerEquations{2} + } <: AbstractCompressibleNavierStokesDiffusion{2, 4} # TODO: parabolic # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations @@ -109,7 +117,8 @@ function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquatio CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, μ, Pr, kappa, - equations, gradient_variables) + equations, + gradient_variables) end # TODO: parabolic @@ -301,12 +310,14 @@ end <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesPrimitive }) - v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, + t, equations) return SVector(u_inner[1], v1, v2, u_inner[4]) end @@ -315,7 +326,8 @@ end <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesPrimitive @@ -324,7 +336,8 @@ end normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) - v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, + t, equations) _, tau_1n, tau_2n, _ = flux_inner # extract fluxes for 2nd and 3rd equations normal_energy_flux = v1 * tau_1n + v2 * tau_2n + normal_heat_flux @@ -335,12 +348,14 @@ end <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesPrimitive }) - v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, + t, equations) T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations) @@ -351,7 +366,8 @@ end <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesPrimitive @@ -371,12 +387,14 @@ end <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesEntropy }) - v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, + t, equations) negative_rho_inv_p = w_inner[4] # w_4 = -rho / p return SVector(w_inner[1], -v1 * negative_rho_inv_p, -v2 * negative_rho_inv_p, @@ -388,7 +406,8 @@ end <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesEntropy @@ -396,7 +415,8 @@ end normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations) - v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, + t, equations) _, tau_1n, tau_2n, _ = flux_inner # extract fluxes for 2nd and 3rd equations normal_energy_flux = v1 * tau_1n + v2 * tau_2n + normal_heat_flux @@ -407,12 +427,14 @@ end <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesEntropy }) - v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, + v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, + t, equations) T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations) @@ -426,10 +448,12 @@ end <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion2D{ GradientVariablesEntropy }) return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4]) end +end # @muladd diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 166b53bf615..de2cad99ea8 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -1,3 +1,10 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + @doc raw""" CompressibleNavierStokesDiffusion3D(equations; mu, Pr, gradient_variables=GradientVariablesPrimitive()) @@ -77,7 +84,8 @@ w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p} This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, - E <: AbstractCompressibleEulerEquations{3}} <: + E <: AbstractCompressibleEulerEquations{3} + } <: AbstractCompressibleNavierStokesDiffusion{3, 5} # TODO: parabolic # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations @@ -109,7 +117,8 @@ function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquatio CompressibleNavierStokesDiffusion3D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one, μ, Pr, kappa, - equations, gradient_variables) + equations, + gradient_variables) end # TODO: parabolic @@ -319,9 +328,12 @@ end @inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion3D) # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. - _, dv1dx, dv2dx, dv3dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) - _, dv1dy, dv2dy, dv3dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) - _, dv1dz, dv2dz, dv3dz, _ = convert_derivative_to_primitive(u, gradients[3], equations) + _, dv1dx, dv2dx, dv3dx, _ = convert_derivative_to_primitive(u, gradients[1], + equations) + _, dv1dy, dv2dy, dv3dy, _ = convert_derivative_to_primitive(u, gradients[2], + equations) + _, dv1dz, dv2dz, dv3dz, _ = convert_derivative_to_primitive(u, gradients[3], + equations) return SVector(dv3dy - dv2dz, dv1dz - dv3dx, dv2dx - dv1dy) end @@ -330,7 +342,8 @@ end <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesPrimitive @@ -345,7 +358,8 @@ end <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesPrimitive @@ -367,7 +381,8 @@ end <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesPrimitive @@ -384,7 +399,8 @@ end <:Isothermal})(flux_inner, u_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesPrimitive @@ -404,7 +420,8 @@ end <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesEntropy @@ -422,7 +439,8 @@ end <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesEntropy @@ -443,7 +461,8 @@ end <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Gradient, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesEntropy @@ -463,7 +482,8 @@ end <:Isothermal})(flux_inner, w_inner, normal::AbstractVector, - x, t, + x, + t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion3D{ GradientVariablesEntropy @@ -471,3 +491,4 @@ end return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4], flux_inner[5]) end +end # @muladd diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 3c2b8855ce8..f00138c698c 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -80,8 +80,8 @@ isdir(outdir) && rm(outdir, recursive=true) @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls_amr.jl"), equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), Prandtl=prandtl_number()), - l2 = [2.527877257772131e-5, 2.5539911566937718e-5, 0.0001211860451244785], - linf = [0.00014663867588948776, 0.00019422448348348196, 0.0009556439394007299] + l2 = [2.5278824700860636e-5, 2.5540078777006958e-5, 0.00012118655083858043], + linf = [0.0001466387075579334, 0.00019422427462629705, 0.0009556446847707178] ) end @@ -90,8 +90,8 @@ isdir(outdir) && rm(outdir, recursive=true) equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), Prandtl=prandtl_number(), gradient_variables = GradientVariablesEntropy()), - l2 = [2.4593699163175966e-5, 2.392863645712634e-5, 0.00011252526651714956], - linf = [0.00011850555445525046, 0.0001898777490968537, 0.0009597561467877824] + l2 = [2.459359632523962e-5, 2.3928390718460263e-5, 0.00011252414117082376], + linf = [0.0001185052018830568, 0.00018987717854305393, 0.0009597503607920999] ) end end