From ab8ecad4249fa718fab0ba6acba60bc5d92b8a75 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Wed, 11 Dec 2024 18:52:40 -1000 Subject: [PATCH] Type stability of functions from examples tree_2d_dgsem (#2145) * Start * Add more * Add more * Add more * Add more * Add more * Add more * Add more * Mark * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * Fix * Apply pi * Fix on second review --------- Co-authored-by: Hendrik Ranocha --- .../elixir_acoustics_gauss_wall.jl | 7 +- .../elixir_acoustics_gaussian_source.jl | 15 ++- .../elixir_acoustics_monopole.jl | 22 ++-- .../elixir_advection_amr_coarsen_twice.jl | 2 +- .../elixir_advection_amr_refine_twice.jl | 2 +- ...ixir_advection_amr_solution_independent.jl | 11 +- .../elixir_advection_callbacks.jl | 2 +- .../elixir_advection_diffusion.jl | 9 +- .../elixir_advection_diffusion_nonperiodic.jl | 7 +- .../elixir_euler_astro_jet_amr.jl | 11 +- .../tree_2d_dgsem/elixir_euler_blast_wave.jl | 11 +- .../elixir_euler_blast_wave_amr.jl | 11 +- .../elixir_euler_blast_wave_pure_fv.jl | 11 +- ...euler_blast_wave_sc_subcell_nonperiodic.jl | 11 +- .../tree_2d_dgsem/elixir_euler_blob_amr.jl | 20 +-- .../tree_2d_dgsem/elixir_euler_blob_mortar.jl | 20 +-- .../elixir_euler_colliding_flow.jl | 13 +- .../elixir_euler_colliding_flow_amr.jl | 13 +- ...uler_colliding_flow_amr_entropy_bounded.jl | 13 +- ...ixir_euler_kelvin_helmholtz_instability.jl | 12 +- ..._euler_kelvin_helmholtz_instability_amr.jl | 12 +- ...in_helmholtz_instability_fjordholm_etal.jl | 56 ++++---- ...kelvin_helmholtz_instability_sc_subcell.jl | 12 +- .../tree_2d_dgsem/elixir_euler_positivity.jl | 17 +-- .../elixir_euler_sedov_blast_wave.jl | 17 +-- ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 17 +-- .../elixir_euler_shockcapturing_subcell.jl | 11 +- ...r_euler_source_terms_amr_refine_coarsen.jl | 11 +- examples/tree_2d_dgsem/elixir_euler_vortex.jl | 17 +-- .../tree_2d_dgsem/elixir_euler_vortex_amr.jl | 21 +-- .../elixir_euler_vortex_mortar.jl | 17 +-- ...ixir_euler_vortex_mortar_shockcapturing.jl | 17 +-- .../elixir_euler_vortex_mortar_split.jl | 17 +-- .../elixir_euler_vortex_shockcapturing.jl | 17 +-- .../tree_2d_dgsem/elixir_euler_warm_bubble.jl | 21 +-- ..._euleracoustics_co-rotating_vortex_pair.jl | 13 +- .../elixir_eulermulti_shock_bubble.jl | 27 ++-- ...ck_bubble_shockcapturing_subcell_minmax.jl | 27 ++-- ...ubble_shockcapturing_subcell_positivity.jl | 27 ++-- .../tree_2d_dgsem/elixir_hypdiff_godunov.jl | 16 ++- .../elixir_hypdiff_harmonic_nonperiodic.jl | 27 ++-- .../elixir_hypdiff_lax_friedrichs.jl | 16 ++- examples/tree_2d_dgsem/elixir_kpp.jl | 16 ++- examples/tree_2d_dgsem/elixir_lbm_couette.jl | 3 +- .../tree_2d_dgsem/elixir_mhd_blast_wave.jl | 29 +++-- .../tree_2d_dgsem/elixir_mhd_orszag_tang.jl | 18 +-- examples/tree_2d_dgsem/elixir_mhd_rotor.jl | 41 +++--- .../elixir_mhd_shockcapturing_subcell.jl | 31 ++--- .../tree_2d_dgsem/elixir_mhdmulti_rotor.jl | 47 +++---- .../elixir_navierstokes_convergence.jl | 120 +++++++++--------- .../elixir_navierstokes_lid_driven_cavity.jl | 17 ++- .../elixir_navierstokes_shearlayer_amr.jl | 13 +- ...elixir_navierstokes_taylor_green_vortex.jl | 9 +- ...erstokes_taylor_green_vortex_sutherland.jl | 20 +-- .../tree_2d_dgsem/elixir_shallowwater_ec.jl | 19 +-- .../tree_2d_dgsem/elixir_shallowwater_wall.jl | 8 +- .../elixir_shallowwater_well_balanced.jl | 17 +-- .../elixir_shallowwater_well_balanced_wall.jl | 17 +-- 58 files changed, 575 insertions(+), 506 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall.jl b/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall.jl index 918c0831fcd..4027bce76d7 100644 --- a/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall.jl +++ b/examples/tree_2d_dgsem/elixir_acoustics_gauss_wall.jl @@ -26,9 +26,10 @@ A Gaussian pulse, used in the `gauss_wall` example elixir in combination with [`boundary_condition_wall`](@ref). Uses the global mean values from `equations`. """ function initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D) - v1_prime = 0.0 - v2_prime = 0.0 - p_prime = exp(-log(2) * (x[1]^2 + (x[2] - 25)^2) / 25) + RealT = eltype(x) + v1_prime = 0 + v2_prime = 0 + p_prime = exp(-log(convert(RealT, 2)) * (x[1]^2 + (x[2] - 25)^2) / 25) prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...) diff --git a/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl b/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl index 71d4f1a9f68..522b5844546 100644 --- a/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl +++ b/examples/tree_2d_dgsem/elixir_acoustics_gaussian_source.jl @@ -3,18 +3,19 @@ using Trixi # Oscillating Gaussian-shaped source terms function source_terms_gauss(u, x, t, equations::AcousticPerturbationEquations2D) - r = 0.1 - A = 1.0 - f = 2.0 + RealT = eltype(u) + r = convert(RealT, 0.1) + A = 1 + f = 2 # Velocity sources - s1 = 0.0 - s2 = 0.0 + s1 = 0 + s2 = 0 # Pressure source - s3 = exp(-(x[1]^2 + x[2]^2) / (2 * r^2)) * A * sin(2 * pi * f * t) + s3 = exp(-(x[1]^2 + x[2]^2) / (2 * r^2)) * A * sinpi(2 * f * t) # Mean sources - s4 = s5 = s6 = s7 = 0.0 + s4 = s5 = s6 = s7 = 0 return SVector(s1, s2, s3, s4, s5, s6, s7) end diff --git a/examples/tree_2d_dgsem/elixir_acoustics_monopole.jl b/examples/tree_2d_dgsem/elixir_acoustics_monopole.jl index d7265775114..416fe051d90 100644 --- a/examples/tree_2d_dgsem/elixir_acoustics_monopole.jl +++ b/examples/tree_2d_dgsem/elixir_acoustics_monopole.jl @@ -20,16 +20,17 @@ Initial condition for the monopole in a boundary layer setup, used in combinatio [`boundary_condition_monopole`](@ref). """ function initial_condition_monopole(x, t, equations::AcousticPerturbationEquations2D) - m = 0.3 # Mach number + RealT = eltype(x) + m = convert(RealT, 0.3) # Mach number - v1_prime = 0.0 - v2_prime = 0.0 - p_prime = 0.0 + v1_prime = 0 + v2_prime = 0 + p_prime = 0 v1_mean = x[2] > 1 ? m : m * (2 * x[2] - 2 * x[2]^2 + x[2]^4) - v2_mean = 0.0 - c_mean = 1.0 - rho_mean = 1.0 + v2_mean = 0 + c_mean = 1 + rho_mean = 1 prim = SVector(v1_prime, v2_prime, p_prime, v1_mean, v2_mean, c_mean, rho_mean) @@ -48,6 +49,7 @@ with [`initial_condition_monopole`](@ref). function boundary_condition_monopole(u_inner, orientation, direction, x, t, surface_flux_function, equations::AcousticPerturbationEquations2D) + RealT = eltype(u_inner) if direction != 3 error("expected direction = 3, got $direction instead") end @@ -56,9 +58,9 @@ function boundary_condition_monopole(u_inner, orientation, direction, x, t, # we use a sinusoidal boundary state for the perturbed variables. For the rest of the -y boundary # we set the boundary state to the inner state and multiply the perturbed velocity in the # y-direction by -1. - if -0.05 <= x[1] <= 0.05 # Monopole - v1_prime = 0.0 - v2_prime = p_prime = sin(2 * pi * t) + if RealT(-0.05) <= x[1] <= RealT(0.05) # Monopole + v1_prime = 0 + v2_prime = p_prime = sinpi(2 * t) prim_boundary = SVector(v1_prime, v2_prime, p_prime, u_inner[4], u_inner[5], u_inner[6], u_inner[7]) diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_coarsen_twice.jl b/examples/tree_2d_dgsem/elixir_advection_amr_coarsen_twice.jl index aa042e7500e..70c392ee21a 100644 --- a/examples/tree_2d_dgsem/elixir_advection_amr_coarsen_twice.jl +++ b/examples/tree_2d_dgsem/elixir_advection_amr_coarsen_twice.jl @@ -28,7 +28,7 @@ function (indicator::IndicatorAlwaysCoarsen)(u::AbstractArray{<:Any, 4}, alpha = indicator.cache.alpha resize!(alpha, nelements(dg, cache)) - alpha .= -1.0 + fill!(alpha, -1) return alpha end diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_refine_twice.jl b/examples/tree_2d_dgsem/elixir_advection_amr_refine_twice.jl index 7b441775204..900ce5cc059 100644 --- a/examples/tree_2d_dgsem/elixir_advection_amr_refine_twice.jl +++ b/examples/tree_2d_dgsem/elixir_advection_amr_refine_twice.jl @@ -28,7 +28,7 @@ function (indicator::IndicatorAlwaysRefine)(u::AbstractArray{<:Any, 4}, alpha = indicator.cache.alpha resize!(alpha, nelements(dg, cache)) - alpha .= 1.0 + fill!(alpha, 1) return alpha end diff --git a/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl index c7412660b0c..b9fbbe1b37c 100644 --- a/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/tree_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -19,16 +19,17 @@ end function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, mesh, equations, dg, cache; t, kwargs...) + RealT = eltype(u) mesh = indicator.cache.mesh alpha = indicator.cache.alpha resize!(alpha, nelements(dg, cache)) #Predict the theoretical center. - advection_velocity = (0.2, -0.7) + advection_velocity = (convert(RealT, 0.2), convert(RealT, -0.7)) center = t .* advection_velocity inner_distance = 1 - outer_distance = 1.85 + outer_distance = convert(RealT, 1.85) #Iterate over all elements for element in eachindex(alpha) @@ -38,16 +39,16 @@ function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, #The geometric shape of the amr should be preserved when the base_level is increased. #This is done by looking at the original coordinates of each cell. - cell_coordinates = original_coordinates(coordinates, 5 / 8) + cell_coordinates = original_coordinates(coordinates, 0.625f0) cell_distance = periodic_distance_2d(cell_coordinates, center, 10) if cell_distance < (inner_distance + outer_distance) / 2 - cell_coordinates = original_coordinates(coordinates, 5 / 16) + cell_coordinates = original_coordinates(coordinates, 0.3125f0) cell_distance = periodic_distance_2d(cell_coordinates, center, 10) end #Set alpha according to cells position inside the circles. target_level = (cell_distance < inner_distance) + (cell_distance < outer_distance) - alpha[element] = target_level / 2 + alpha[element] = convert(RealT, target_level) / 2 end return alpha end diff --git a/examples/tree_2d_dgsem/elixir_advection_callbacks.jl b/examples/tree_2d_dgsem/elixir_advection_callbacks.jl index 5b58f92fe29..9f9c087b1eb 100644 --- a/examples/tree_2d_dgsem/elixir_advection_callbacks.jl +++ b/examples/tree_2d_dgsem/elixir_advection_callbacks.jl @@ -130,7 +130,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2cons) -example_callback = TrixiExtensionExample.ExampleStepCallback(message = "안녕하세요?") +example_callback = TrixiExtensionExample.ExampleStepCallback(message = "Initializing callback") stepsize_callback = StepsizeCallback(cfl = 1.6) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl index 1f765ff3564..a5c2ca37b80 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl @@ -25,14 +25,15 @@ mesh = TreeMesh(coordinates_min, coordinates_max, function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t nu = diffusivity() - c = 1.0 - A = 0.5 + c = 1 + A = 0.5f0 L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) return SVector(scalar) end diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl index a3fc820af70..06e2cd22b3e 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -28,14 +28,15 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # to numerical partial differential equations. # [DOI](https://doi.org/10.1007/978-3-319-41640-3_6). function initial_condition_eriksson_johnson(x, t, equations) + RealT = eltype(x) l = 4 epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) - r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) - s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + r1 = (1 + sqrt(1 + 4 * convert(RealT, pi)^2 * epsilon^2)) / (2 * epsilon) + s1 = (1 - sqrt(1 + 4 * convert(RealT, pi)^2 * epsilon^2)) / (2 * epsilon) u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) + - cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) + cospi(x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) return SVector{1}(u) end initial_condition = initial_condition_eriksson_johnson diff --git a/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl b/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl index 1393f4a6cc8..7fbaa8c533b 100644 --- a/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl +++ b/examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl @@ -13,18 +13,19 @@ equations = CompressibleEulerEquations2D(gamma) # https://tinyurl.com/c76fjtx4 # Mach = 2000 jet function initial_condition_astro_jet(x, t, equations::CompressibleEulerEquations2D) + RealT = eltype(x) @unpack gamma = equations - rho = 0.5 + rho = 0.5f0 v1 = 0 v2 = 0 - p = 0.4127 + p = convert(RealT, 0.4127) # add inflow for t>0 at x=-0.5 # domain size is [-0.5,+0.5]^2 - if (t > 0) && (x[1] ≈ -0.5) && (abs(x[2]) < 0.05) - rho = 5 + if (t > 0) && (x[1] ≈ -0.5f0) && (abs(x[2]) < RealT(0.05)) + rho = 5.0f0 v1 = 800 # about Mach number Ma = 2000 v2 = 0 - p = 0.4127 + p = convert(RealT, 0.4127) end return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave.jl index ccd7b54086b..005ba442854 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave.jl @@ -18,7 +18,8 @@ A medium blast wave taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -26,10 +27,10 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation sin_phi, cos_phi = sincos(phi) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-3 : 1.245 + rho = r > 0.5f0 ? one(RealT) : RealT(1.1691) + v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi + p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245) return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr.jl index d32c2e51b06..6ee73277e76 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_amr.jl @@ -18,7 +18,8 @@ A medium blast wave taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -26,10 +27,10 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation sin_phi, cos_phi = sincos(phi) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-3 : 1.245 + rho = r > 0.5f0 ? one(RealT) : RealT(1.1691) + v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi + p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245) return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_pure_fv.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_pure_fv.jl index a2392d05e5a..c12f775008f 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_pure_fv.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_pure_fv.jl @@ -18,7 +18,8 @@ A medium blast wave taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -26,10 +27,10 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation sin_phi, cos_phi = sincos(phi) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-3 : 1.245 + rho = r > 0.5f0 ? one(RealT) : RealT(1.1691) + v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi + p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245) return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl index 21a89fdcbe6..47c31aceac9 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl @@ -18,7 +18,8 @@ A medium blast wave taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -26,10 +27,10 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation sin_phi, cos_phi = sincos(phi) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-3 : 1.245 + rho = r > 0.5f0 ? one(RealT) : RealT(1.1691) + v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi + p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245) return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_euler_blob_amr.jl b/examples/tree_2d_dgsem/elixir_euler_blob_amr.jl index 2b3659017a3..7ad035e0a9e 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blob_amr.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blob_amr.jl @@ -23,17 +23,18 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D) # resolution 128^2, 256^2 # domain size is [-20.0,20.0]^2 # gamma = 5/3 for this test case - R = 1.0 # radius of the blob + RealT = eltype(x) + R = 1 # radius of the blob # background density - dens0 = 1.0 - Chi = 10.0 # density contrast + dens0 = 1 + Chi = convert(RealT, 10) # density contrast # reference time of characteristic growth of KH instability equal to 1.0 - tau_kh = 1.0 - tau_cr = tau_kh / 1.6 # crushing time + tau_kh = 1 + tau_cr = tau_kh / convert(RealT, 1.6) # crushing time # determine background velocity velx0 = 2 * R * sqrt(Chi) / tau_cr - vely0 = 0.0 - Ma0 = 2.7 # background flow Mach number Ma=v/c + vely0 = 0 + Ma0 = convert(RealT, 2.7) # background flow Mach number Ma=v/c c = velx0 / Ma0 # sound speed # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density p0 = c * c * dens0 / equations.gamma @@ -45,9 +46,10 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D) slope = 2 # density blob dens = dens0 + - (Chi - 1) * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) + (Chi - 1) * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) # velocity blob is zero - velx = velx0 - velx0 * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) + velx = velx0 - + velx0 * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) return prim2cons(SVector(dens, velx, vely0, p0), equations) end initial_condition = initial_condition_blob diff --git a/examples/tree_2d_dgsem/elixir_euler_blob_mortar.jl b/examples/tree_2d_dgsem/elixir_euler_blob_mortar.jl index 8bd5db00c9a..1d363a241d2 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blob_mortar.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blob_mortar.jl @@ -23,17 +23,18 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D) # resolution 128^2, 256^2 # domain size is [-20.0,20.0]^2 # gamma = 5/3 for this test case - R = 1.0 # radius of the blob + RealT = eltype(x) + R = 1 # radius of the blob # background density - dens0 = 1.0 - Chi = 10.0 # density contrast + dens0 = 1 + Chi = convert(RealT, 10) # density contrast # reference time of characteristic growth of KH instability equal to 1.0 - tau_kh = 1.0 - tau_cr = tau_kh / 1.6 # crushing time + tau_kh = 1 + tau_cr = tau_kh / convert(RealT, 1.6) # crushing time # determine background velocity velx0 = 2 * R * sqrt(Chi) / tau_cr - vely0 = 0.0 - Ma0 = 2.7 # background flow Mach number Ma=v/c + vely0 = 0 + Ma0 = convert(RealT, 2.7) # background flow Mach number Ma=v/c c = velx0 / Ma0 # sound speed # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density p0 = c * c * dens0 / equations.gamma @@ -45,9 +46,10 @@ function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D) slope = 2 # density blob dens = dens0 + - (Chi - 1) * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) + (Chi - 1) * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) # velocity blob is zero - velx = velx0 - velx0 * 0.5 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) + velx = velx0 - + velx0 * 0.5f0 * (1 + (tanh(slope * (r + R)) - (tanh(slope * (r - R)) + 1))) return prim2cons(SVector(dens, velx, vely0, p0), equations) end initial_condition = initial_condition_blob diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl index 984ac3ff1f6..7396f11a684 100644 --- a/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl +++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow.jl @@ -15,20 +15,21 @@ function initial_condition_colliding_flow_astro(x, t, # change discontinuity to tanh # resolution 128^2 elements (refined close to the interface) and polydeg=3 (total of 512^2 DOF) # domain size is [-64,+64]^2 + RealT = eltype(x) @unpack gamma = equations # the quantities are chosen such, that they are as close as possible to the astro examples # keep in mind, that in the astro example, the physical units are weird (parsec, mega years, ...) - rho = 0.0247 - c = 0.2 + rho = convert(RealT, 0.0247) + c = convert(RealT, 0.2) p = c^2 / gamma * rho - vel = 13.907432274789372 - slope = 1.0 + vel = convert(RealT, 13.907432274789372) + slope = 1 v1 = -vel * tanh(slope * x[1]) # add small initial disturbance to the field, but only close to the interface if abs(x[1]) < 10 - v1 = v1 * (1 + 0.01 * sin(pi * x[2])) + v1 = v1 * (1 + RealT(0.01) * sinpi(x[2])) end - v2 = 0.0 + v2 = 0 return prim2cons(SVector(rho, v1, v2, p), equations) end initial_condition = initial_condition_colliding_flow_astro diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl index a9eb671929f..a68fee19178 100644 --- a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl +++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl @@ -15,20 +15,21 @@ function initial_condition_colliding_flow_astro(x, t, # change discontinuity to tanh # resolution 128^2 elements (refined close to the interface) and polydeg=3 (total of 512^2 DOF) # domain size is [-64,+64]^2 + RealT = eltype(x) @unpack gamma = equations # the quantities are chosen such, that they are as close as possible to the astro examples # keep in mind, that in the astro example, the physical units are weird (parsec, mega years, ...) - rho = 0.0247 - c = 0.2 + rho = convert(RealT, 0.0247) + c = convert(RealT, 0.2) p = c^2 / gamma * rho - vel = 13.907432274789372 - slope = 1.0 + vel = convert(RealT, 13.907432274789372) + slope = 1 v1 = -vel * tanh(slope * x[1]) # add small initial disturbance to the field, but only close to the interface if abs(x[1]) < 10 - v1 = v1 * (1 + 0.01 * sin(pi * x[2])) + v1 = v1 * (1 + RealT(0.01) * sinpi(x[2])) end - v2 = 0.0 + v2 = 0 return prim2cons(SVector(rho, v1, v2, p), equations) end initial_condition = initial_condition_colliding_flow_astro diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl index 6e7cac9a175..61474efd82c 100644 --- a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl +++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl @@ -15,20 +15,21 @@ function initial_condition_colliding_flow_astro(x, t, # change discontinuity to tanh # resolution 128^2 elements (refined close to the interface) and polydeg=3 (total of 512^2 DOF) # domain size is [-64,+64]^2 + RealT = eltype(x) @unpack gamma = equations # the quantities are chosen such, that they are as close as possible to the astro examples # keep in mind, that in the astro example, the physical units are weird (parsec, mega years, ...) - rho = 0.0247 - c = 0.2 + rho = convert(RealT, 0.0247) + c = convert(RealT, 0.2) p = c^2 / gamma * rho - vel = 13.907432274789372 - slope = 1.0 + vel = convert(RealT, 13.907432274789372) + slope = 1 v1 = -vel * tanh(slope * x[1]) # add small initial disturbance to the field, but only close to the interface if abs(x[1]) < 10 - v1 = v1 * (1 + 0.01 * sin(pi * x[2])) + v1 = v1 * (1 + RealT(0.01) * sinpi(x[2])) end - v2 = 0.0 + v2 = 0 return prim2cons(SVector(rho, v1, v2, p), equations) end initial_condition = initial_condition_colliding_flow_astro diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl index 5e6b1e0cc0d..dfeb0c7f13b 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl @@ -21,13 +21,13 @@ function initial_condition_kelvin_helmholtz_instability(x, t, # change discontinuity to tanh # typical resolution 128^2, 256^2 # domain size is [-1,+1]^2 + RealT = eltype(x) slope = 15 - amplitude = 0.02 - B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5) - rho = 0.5 + 0.75 * B - v1 = 0.5 * (B - 1) - v2 = 0.1 * sin(2 * pi * x[1]) - p = 1.0 + B = tanh(slope * x[2] + 7.5f0) - tanh(slope * x[2] - 7.5f0) + rho = 0.5f0 + 0.75f0 * B + v1 = 0.5f0 * (B - 1) + v2 = convert(RealT, 0.1) * sinpi(2 * x[1]) + p = 1 return prim2cons(SVector(rho, v1, v2, p), equations) end initial_condition = initial_condition_kelvin_helmholtz_instability diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr.jl index 5c237835cc5..3a4ff051a06 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr.jl @@ -21,13 +21,13 @@ function initial_condition_kelvin_helmholtz_instability(x, t, # change discontinuity to tanh # typical resolution 128^2, 256^2 # domain size is [-1,+1]^2 + RealT = eltype(x) slope = 15 - amplitude = 0.02 - B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5) - rho = 0.5 + 0.75 * B - v1 = 0.5 * (B - 1) - v2 = 0.1 * sin(2 * pi * x[1]) - p = 1.0 + B = tanh(slope * x[2] + 7.5f0) - tanh(slope * x[2] - 7.5f0) + rho = 0.5f0 + 0.75f0 * B + v1 = 0.5f0 * (B - 1) + v2 = convert(RealT, 0.1) * sinpi(2 * x[1]) + p = 1 return prim2cons(SVector(rho, v1, v2, p), equations) end initial_condition = initial_condition_kelvin_helmholtz_instability diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl index 5b2b80d84f3..13388da6269 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl @@ -33,43 +33,49 @@ function initial_condition_kelvin_helmholtz_instability_fjordholm_etal(x, t, # b1 = (rand(rng, m) .- 0.5) .* pi # b2 = (rand(rng, m) .- 0.5) .* pi + RealT = eltype(x) m = 10 - a1 = [0.04457096674422902, 0.03891512410182607, 0.0030191053979293433, - 0.0993913172320319, - 0.1622302137588842, 0.1831383653456182, 0.11758003014101702, 0.07964318348142958, - 0.0863245324711805, 0.18518716132585408] - a2 = [0.061688440856337096, 0.23000237877135882, 0.04453793881833177, - 0.19251530387370916, - 0.11107917357941084, 0.05898041974649702, 0.09949312336096268, 0.07022276346006465, - 0.10670366489014596, 0.02477679264318211] - b1 = [0.06582340543754152, 0.9857886297001535, 0.8450452205037154, -1.279648120993805, - 0.45454198915209526, -0.13359370986823993, 0.07062615913363897, -1.0097986278512623, - 1.0810669017430343, -0.14207309803877177] - b2 = [-1.1376882185131414, -1.4798197129947765, 0.6139290513283818, -0.3319087388365522, - 0.14633328999192285, -0.06373231463100072, -0.6270101051216724, 0.13941252226261905, - -1.0337526453303645, 1.0441408867083155] - Y1 = 0.0 - Y2 = 0.0 + a1 = RealT[0.04457096674422902, 0.03891512410182607, 0.0030191053979293433, + 0.0993913172320319, 0.1622302137588842, 0.1831383653456182, + 0.11758003014101702, 0.07964318348142958, 0.0863245324711805, + 0.18518716132585408] + + a2 = RealT[0.061688440856337096, 0.23000237877135882, 0.04453793881833177, + 0.19251530387370916, 0.11107917357941084, 0.05898041974649702, + 0.09949312336096268, 0.07022276346006465, 0.10670366489014596, + 0.02477679264318211] + + b1 = RealT[0.06582340543754152, 0.9857886297001535, 0.8450452205037154, + -1.279648120993805, 0.45454198915209526, -0.13359370986823993, + 0.07062615913363897, -1.0097986278512623, 1.0810669017430343, + -0.14207309803877177] + + b2 = RealT[-1.1376882185131414, -1.4798197129947765, 0.6139290513283818, + -0.3319087388365522, 0.14633328999192285, -0.06373231463100072, + -0.6270101051216724, 0.13941252226261905, -1.0337526453303645, + 1.0441408867083155] + Y1 = zero(RealT) + Y2 = zero(RealT) for n in 1:m - Y1 += a1[n] * cos(b1[n] + 2 * n * pi * x[1]) - Y2 += a2[n] * cos(b2[n] + 2 * n * pi * x[1]) + Y1 += a1[n] * cos(b1[n] + 2 * n * convert(RealT, pi) * x[1]) + Y2 += a2[n] * cos(b2[n] + 2 * n * convert(RealT, pi) * x[1]) end - J1 = 0.25 - J2 = 0.75 - epsilon = 0.01 + J1 = 0.25f0 + J2 = 0.75f0 + epsilon = convert(RealT, 0.01) I1 = J1 + epsilon * Y1 I2 = J2 + epsilon * Y2 if (x[2] > I1) && (x[2] < I2) rho = 2 - v1 = -0.5 + v1 = -0.5f0 else rho = 1 - v1 = 0.5 + v1 = 0.5f0 end - v2 = 0 - p = 2.5 + v2 = zero(RealT) + p = 2.5f0 return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl index 9e9fb45e7d1..578e51e53dc 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl @@ -21,13 +21,13 @@ function initial_condition_kelvin_helmholtz_instability(x, t, # change discontinuity to tanh # typical resolution 128^2, 256^2 # domain size is [-1,+1]^2 + RealT = eltype(x) slope = 15 - amplitude = 0.02 - B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5) - rho = 0.5 + 0.75 * B - v1 = 0.5 * (B - 1) - v2 = 0.1 * sin(2 * pi * x[1]) - p = 1.0 + B = tanh(slope * x[2] + 7.5f0) - tanh(slope * x[2] - 7.5f0) + rho = 0.5f0 + 0.75f0 * B + v1 = 0.5f0 * (B - 1) + v2 = convert(RealT, 0.1) * sinpi(2 * x[1]) + p = 1 return prim2cons(SVector(rho, v1, v2, p), equations) end initial_condition = initial_condition_kelvin_helmholtz_instability diff --git a/examples/tree_2d_dgsem/elixir_euler_positivity.jl b/examples/tree_2d_dgsem/elixir_euler_positivity.jl index 6fec4c1bf9b..73393e33c7f 100644 --- a/examples/tree_2d_dgsem/elixir_euler_positivity.jl +++ b/examples/tree_2d_dgsem/elixir_euler_positivity.jl @@ -15,23 +15,24 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 - r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup - E = 1.0 - p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0^2) + p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables - rho = 1.0 - v1 = 0.0 - v2 = 0.0 + rho = 1 + v1 = 0 + v2 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, v2, p), equations) diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl index 5b8959b97d1..018f81c21f3 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -15,23 +15,24 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 - r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup - E = 1.0 - p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0^2) + p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables - rho = 1.0 - v1 = 0.0 - v2 = 0.0 + rho = 1 + v1 = 0 + v2 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, v2, p), equations) diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 1dc0586ec7e..5148b5c5d2f 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -15,23 +15,24 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 - r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup - E = 1.0 - p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0^2) + p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables - rho = 1.0 - v1 = 0.0 - v2 = 0.0 + rho = 1 + v1 = 0 + v2 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, v2, p), equations) diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 44e63a0872e..a50689a2102 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -18,7 +18,8 @@ A medium blast wave (modified to lower density and higher pressure) taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> modified to lower density, higher pressure # Set up polar coordinates - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -26,10 +27,10 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation sin_phi, cos_phi = sincos(phi) # Calculate primitive variables "normal" medium blast wave - rho = r > 0.5 ? 0.1 : 0.2691 # rho = r > 0.5 ? 1 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-1 : 1.245 # p = r > 0.5 ? 1.0E-3 : 1.245 + rho = r > 0.5f0 ? RealT(0.1) : RealT(0.2691) # rho = r > 0.5 ? 1 : 1.1691 + v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * sin_phi + p = r > 0.5f0 ? RealT(1.0E-1) : RealT(1.245) # p = r > 0.5 ? 1.0E-3 : 1.245 return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_euler_source_terms_amr_refine_coarsen.jl b/examples/tree_2d_dgsem/elixir_euler_source_terms_amr_refine_coarsen.jl index 28ab4cec1d3..f2984f49cac 100644 --- a/examples/tree_2d_dgsem/elixir_euler_source_terms_amr_refine_coarsen.jl +++ b/examples/tree_2d_dgsem/elixir_euler_source_terms_amr_refine_coarsen.jl @@ -25,17 +25,18 @@ end function (indicator::IndicatorRefineCoarsen)(u::AbstractArray{<:Any, 4}, mesh, equations, dg, cache; t, kwargs...) + RealT = eltype(u) alpha = indicator.cache.alpha resize!(alpha, nelements(dg, cache)) - if t >= 0.7 && t < 1.0 + if t >= RealT(0.7) && t < 1 # Refine to max level - alpha .= 1.0 - elseif t >= 1.0 + fill!(alpha, 1) + elseif t >= 1 # Coarsen to base level - alpha .= -1.0 + fill!(alpha, -1) else - alpha .= 0.0 + fill!(alpha, 0) end return alpha diff --git a/examples/tree_2d_dgsem/elixir_euler_vortex.jl b/examples/tree_2d_dgsem/elixir_euler_vortex.jl index c87d6f49ba8..67fa302429c 100644 --- a/examples/tree_2d_dgsem/elixir_euler_vortex.jl +++ b/examples/tree_2d_dgsem/elixir_euler_vortex.jl @@ -21,17 +21,18 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # for error convergence: make sure that the end time is such that the vortex is back at the initial state!! # for the current velocity and domain size: t_end should be a multiple of 20s # initial center of the vortex - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) # size and strength of the vortex - iniamplitude = 5.0 + iniamplitude = 5 # base flow - rho = 1.0 - v1 = 1.0 - v2 = 1.0 + rho = 1 + v1 = 1 + v2 = 1 vel = SVector(v1, v2) - p = 25.0 + p = convert(RealT, 25) rt = p / rho # ideal gas equation - t_loc = 0.0 + t_loc = 0 cent = inicenter + vel * t_loc # advection of center # ATTENTION: handle periodic BC, but only for v1 = v2 = 1.0 (!!!!) cent = x - cent # distance to center point @@ -39,7 +40,7 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # cross product with iniaxis = [0, 0, 1] cent = SVector(-cent[2], cent[1]) r2 = cent[1]^2 + cent[2]^2 - du = iniamplitude / (2 * π) * exp(0.5 * (1 - r2)) # vel. perturbation + du = iniamplitude / (2 * convert(RealT, pi)) * exp(0.5f0 * (1 - r2)) # vel. perturbation dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1)) vel = vel + du * cent diff --git a/examples/tree_2d_dgsem/elixir_euler_vortex_amr.jl b/examples/tree_2d_dgsem/elixir_euler_vortex_amr.jl index fd25defd417..82ada17a93e 100644 --- a/examples/tree_2d_dgsem/elixir_euler_vortex_amr.jl +++ b/examples/tree_2d_dgsem/elixir_euler_vortex_amr.jl @@ -29,10 +29,10 @@ function (indicator_vortex::IndicatorVortex)(u::AbstractArray{<:Any, 4}, alpha = indicator_vortex.cache.alpha resize!(alpha, nelements(dg, cache)) - # get analytical vortex center (based on assumption that center=[0.0,0.0] + # get analytical vortex center (based on assumption that center=[0.0, 0.0] # at t=0.0 and that we stop after one period) domain_length = mesh.tree.length_level_0 - if t < 0.5 * domain_length + if t < 0.5f0 * domain_length center = (t, t) else center = (t - domain_length, t - domain_length) @@ -84,17 +84,18 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # for error convergence: make sure that the end time is such that the vortex is back at the initial state!! # for the current velocity and domain size: t_end should be a multiple of 20s # initial center of the vortex - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) # size and strength of the vortex - iniamplitude = 5.0 + iniamplitude = 5 # base flow - rho = 1.0 - v1 = 1.0 - v2 = 1.0 + rho = 1 + v1 = 1 + v2 = 1 vel = SVector(v1, v2) - p = 25.0 + p = convert(RealT, 25) rt = p / rho # ideal gas equation - t_loc = 0.0 + t_loc = 0 cent = inicenter + vel * t_loc # advection of center # ATTENTION: handle periodic BC, but only for v1 = v2 = 1.0 (!!!!) @@ -104,7 +105,7 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # cross product with iniaxis = [0, 0, 1] cent = SVector(-cent[2], cent[1]) r2 = cent[1]^2 + cent[2]^2 - du = iniamplitude / (2 * π) * exp(0.5 * (1 - r2)) # vel. perturbation + du = iniamplitude / (2 * convert(RealT, pi)) * exp(0.5f0 * (1 - r2)) # vel. perturbation dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1)) vel = vel + du * cent diff --git a/examples/tree_2d_dgsem/elixir_euler_vortex_mortar.jl b/examples/tree_2d_dgsem/elixir_euler_vortex_mortar.jl index 858799d2d3d..ba4d76aabc1 100644 --- a/examples/tree_2d_dgsem/elixir_euler_vortex_mortar.jl +++ b/examples/tree_2d_dgsem/elixir_euler_vortex_mortar.jl @@ -21,17 +21,18 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # for error convergence: make sure that the end time is such that the vortex is back at the initial state!! # for the current velocity and domain size: t_end should be a multiple of 20s # initial center of the vortex - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) # size and strength of the vortex - iniamplitude = 5.0 + iniamplitude = 5 # base flow - rho = 1.0 - v1 = 1.0 - v2 = 1.0 + rho = 1 + v1 = 1 + v2 = 1 vel = SVector(v1, v2) - p = 25.0 + p = convert(RealT, 25) rt = p / rho # ideal gas equation - t_loc = 0.0 + t_loc = 0 cent = inicenter + vel * t_loc # advection of center # ATTENTION: handle periodic BC, but only for v1 = v2 = 1.0 (!!!!) cent = x - cent # distance to center point @@ -39,7 +40,7 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # cross product with iniaxis = [0, 0, 1] cent = SVector(-cent[2], cent[1]) r2 = cent[1]^2 + cent[2]^2 - du = iniamplitude / (2 * π) * exp(0.5 * (1 - r2)) # vel. perturbation + du = iniamplitude / (2 * convert(RealT, pi)) * exp(0.5f0 * (1 - r2)) # vel. perturbation dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1)) vel = vel + du * cent diff --git a/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_shockcapturing.jl b/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_shockcapturing.jl index 026f6d1462c..3d210f103ee 100644 --- a/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_shockcapturing.jl +++ b/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_shockcapturing.jl @@ -21,17 +21,18 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # for error convergence: make sure that the end time is such that the vortex is back at the initial state!! # for the current velocity and domain size: t_end should be a multiple of 20s # initial center of the vortex - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) # size and strength of the vortex - iniamplitude = 5.0 + iniamplitude = 5 # base flow - rho = 1.0 - v1 = 1.0 - v2 = 1.0 + rho = 1 + v1 = 1 + v2 = 1 vel = SVector(v1, v2) - p = 25.0 + p = convert(RealT, 25) rt = p / rho # ideal gas equation - t_loc = 0.0 + t_loc = 0 cent = inicenter + vel * t_loc # advection of center # ATTENTION: handle periodic BC, but only for v1 = v2 = 1.0 (!!!!) cent = x - cent # distance to center point @@ -39,7 +40,7 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # cross product with iniaxis = [0, 0, 1] cent = SVector(-cent[2], cent[1]) r2 = cent[1]^2 + cent[2]^2 - du = iniamplitude / (2 * π) * exp(0.5 * (1 - r2)) # vel. perturbation + du = iniamplitude / (2 * convert(RealT, pi)) * exp(0.5f0 * (1 - r2)) # vel. perturbation dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1)) vel = vel + du * cent diff --git a/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_split.jl b/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_split.jl index d719e01fd7c..69635b4d2cb 100644 --- a/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_split.jl +++ b/examples/tree_2d_dgsem/elixir_euler_vortex_mortar_split.jl @@ -21,17 +21,18 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # for error convergence: make sure that the end time is such that the vortex is back at the initial state!! # for the current velocity and domain size: t_end should be a multiple of 20s # initial center of the vortex - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) # size and strength of the vortex - iniamplitude = 5.0 + iniamplitude = 5 # base flow - rho = 1.0 - v1 = 1.0 - v2 = 1.0 + rho = 1 + v1 = 1 + v2 = 1 vel = SVector(v1, v2) - p = 25.0 + p = convert(RealT, 25) rt = p / rho # ideal gas equation - t_loc = 0.0 + t_loc = 0 cent = inicenter + vel * t_loc # advection of center # ATTENTION: handle periodic BC, but only for v1 = v2 = 1.0 (!!!!) cent = x - cent # distance to center point @@ -39,7 +40,7 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # cross product with iniaxis = [0, 0, 1] cent = SVector(-cent[2], cent[1]) r2 = cent[1]^2 + cent[2]^2 - du = iniamplitude / (2 * π) * exp(0.5 * (1 - r2)) # vel. perturbation + du = iniamplitude / (2 * convert(RealT, pi)) * exp(0.5f0 * (1 - r2)) # vel. perturbation dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1)) vel = vel + du * cent diff --git a/examples/tree_2d_dgsem/elixir_euler_vortex_shockcapturing.jl b/examples/tree_2d_dgsem/elixir_euler_vortex_shockcapturing.jl index 99e4b090633..8280bfbf01d 100644 --- a/examples/tree_2d_dgsem/elixir_euler_vortex_shockcapturing.jl +++ b/examples/tree_2d_dgsem/elixir_euler_vortex_shockcapturing.jl @@ -21,17 +21,18 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # for error convergence: make sure that the end time is such that the vortex is back at the initial state!! # for the current velocity and domain size: t_end should be a multiple of 20s # initial center of the vortex - inicenter = SVector(0.0, 0.0) + RealT = eltype(x) + inicenter = SVector(0, 0) # size and strength of the vortex - iniamplitude = 5.0 + iniamplitude = 5 # base flow - rho = 1.0 - v1 = 1.0 - v2 = 1.0 + rho = 1 + v1 = 1 + v2 = 1 vel = SVector(v1, v2) - p = 25.0 + p = convert(RealT, 25) rt = p / rho # ideal gas equation - t_loc = 0.0 + t_loc = 0 cent = inicenter + vel * t_loc # advection of center # ATTENTION: handle periodic BC, but only for v1 = v2 = 1.0 (!!!!) cent = x - cent # distance to center point @@ -39,7 +40,7 @@ function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerE # cross product with iniaxis = [0, 0, 1] cent = SVector(-cent[2], cent[1]) r2 = cent[1]^2 + cent[2]^2 - du = iniamplitude / (2 * π) * exp(0.5 * (1 - r2)) # vel. perturbation + du = iniamplitude / (2 * convert(RealT, pi)) * exp(0.5f0 * (1 - r2)) # vel. perturbation dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1)) vel = vel + du * cent diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl index 2632e80db71..64f222df2f9 100644 --- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -28,21 +28,22 @@ end # Initial condition function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + RealT = eltype(x) @unpack g, c_p, c_v = setup # center of perturbation - center_x = 10000.0 - center_z = 2000.0 + center_x = 10000 + center_z = 2000 # radius of perturbation - radius = 2000.0 + radius = 2000 # distance of current x to center of perturbation r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) # perturbation in potential temperature - potential_temperature_ref = 300.0 - potential_temperature_perturbation = 0.0 + potential_temperature_ref = 300 + potential_temperature_perturbation = zero(RealT) if r <= radius - potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 + potential_temperature_perturbation = 2 * cospi(0.5f0 * r / radius)^2 end potential_temperature = potential_temperature_ref + potential_temperature_perturbation @@ -50,7 +51,7 @@ function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) exner = 1 - g / (c_p * potential_temperature) * x[2] # pressure - p_0 = 100_000.0 # reference pressure + p_0 = 100_000 # reference pressure R = c_p - c_v # gas constant (dry air) p = p_0 * exner^(c_p / R) @@ -60,9 +61,9 @@ function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) # density rho = p / (R * T) - v1 = 20.0 - v2 = 0.0 - E = c_v * T + 0.5 * (v1^2 + v2^2) + v1 = 20 + v2 = 0 + E = c_v * T + 0.5f0 * (v1^2 + v2^2) return SVector(rho, rho * v1, rho * v2, rho * E) end diff --git a/examples/tree_2d_dgsem/elixir_euleracoustics_co-rotating_vortex_pair.jl b/examples/tree_2d_dgsem/elixir_euleracoustics_co-rotating_vortex_pair.jl index ea81bd049e4..ec1be83357f 100644 --- a/examples/tree_2d_dgsem/elixir_euleracoustics_co-rotating_vortex_pair.jl +++ b/examples/tree_2d_dgsem/elixir_euleracoustics_co-rotating_vortex_pair.jl @@ -30,9 +30,10 @@ end # Analytical flow solution, used for the initial condition of the flow simulation function velocity(x, t, vortex_pair::VortexPair) + RealT = eltype(x) @unpack r0, rc, circulation = vortex_pair - omega = circulation / (4 * pi * r0^2) + omega = circulation / (4 * convert(RealT, pi) * r0^2) si, co = sincos(omega * t) b = SVector(r0 * co, r0 * si) # vortex centers are b and -b z_plus = x - b @@ -47,9 +48,9 @@ function velocity(x, t, vortex_pair::VortexPair) si_plus, co_plus = sincos(theta_plus) si_minus, co_minus = sincos(theta_minus) - v1 = -circulation / (2 * pi) * (r_plus / (rc^2 + r_plus^2) * si_plus + + v1 = -circulation / (2 * convert(RealT, pi)) * (r_plus / (rc^2 + r_plus^2) * si_plus + r_minus / (rc^2 + r_minus^2) * si_minus) - v2 = circulation / (2 * pi) * (r_plus / (rc^2 + r_plus^2) * co_plus + + v2 = circulation / (2 * convert(RealT, pi)) * (r_plus / (rc^2 + r_plus^2) * co_plus + r_minus / (rc^2 + r_minus^2) * co_minus) return SVector(v1, v2) @@ -69,7 +70,7 @@ function (initial_condition::InitialCondition)(x, t, v = velocity(x, t, vortex_pair) p0 = rho0 * c0^2 / gamma - p = p0 - 0.5 * (gamma - 1) / gamma * sum(v .^ 2) # Bernoulli's principle + p = p0 - 0.5f0 * (gamma - 1) / gamma * sum(v .^ 2) # Bernoulli's principle prim = SVector(rho0, v[1], v[2], p) return prim2cons(prim, equations) @@ -341,9 +342,9 @@ summary_callback() ############################################################################### # set up coupled semidiscretization -source_region(x) = sum(abs2, x) < 6.0^2 # calculate sources within radius 6 around origin +source_region(x) = sum(abs2, x) < 6^2 # calculate sources within radius 6 around origin # gradually reduce acoustic source term amplitudes to zero, starting at radius 5 -weights(x) = sum(abs2, x) < 5.0^2 ? 1.0 : cospi(0.5 * (norm(x) - 5.0)) +weights(x) = sum(abs2, x) < 5^2 ? one(eltype(x)) : cospi(0.5f0 * (norm(x) - 5)) semi = SemidiscretizationEulerAcoustics(semi_acoustics, semi_euler, source_region = source_region, weights = weights) diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble.jl index c6ed07dcda1..08346a39af8 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble.jl @@ -22,46 +22,47 @@ function initial_condition_shock_bubble(x, t, # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972 # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf # typical domain is rectangular, we change it to a square, as Trixi can only do squares + RealT = eltype(x) @unpack gas_constants = equations # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving - delta = 0.03 + delta = convert(RealT, 0.03) # Region I rho1_1 = delta - rho2_1 = 1.225 * gas_constants[1] / gas_constants[2] - delta - v1_1 = zero(delta) - v2_1 = zero(delta) + rho2_1 = RealT(1.225) * gas_constants[1] / gas_constants[2] - delta + v1_1 = zero(RealT) + v2_1 = zero(RealT) p_1 = 101325 # Region II - rho1_2 = 1.225 - delta + rho1_2 = RealT(1.225) - delta rho2_2 = delta - v1_2 = zero(delta) - v2_2 = zero(delta) + v1_2 = zero(RealT) + v2_2 = zero(RealT) p_2 = 101325 # Region III - rho1_3 = 1.6861 - delta + rho1_3 = RealT(1.6861) - delta rho2_3 = delta - v1_3 = -113.5243 - v2_3 = zero(delta) + v1_3 = -RealT(113.5243) + v2_3 = zero(RealT) p_3 = 159060 # Set up Region I & II: - inicenter = SVector(zero(delta), zero(delta)) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - if (x[1] > 0.50) + if (x[1] > 0.5f0) # Set up Region III rho1 = rho1_3 rho2 = rho2_3 v1 = v1_3 v2 = v2_3 p = p_3 - elseif (r < 0.25) + elseif (r < 0.25f0) # Set up Region I rho1 = rho1_1 rho2 = rho2_1 diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl index 9bf310d7536..20fec4e83bd 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl @@ -22,46 +22,47 @@ function initial_condition_shock_bubble(x, t, # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972 # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf # typical domain is rectangular, we change it to a square, as Trixi can only do squares + RealT = eltype(x) @unpack gas_constants = equations # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving - delta = 0.03 + delta = convert(RealT, 0.03) # Region I rho1_1 = delta - rho2_1 = 1.225 * gas_constants[1] / gas_constants[2] - delta - v1_1 = zero(delta) - v2_1 = zero(delta) + rho2_1 = RealT(1.225) * gas_constants[1] / gas_constants[2] - delta + v1_1 = zero(RealT) + v2_1 = zero(RealT) p_1 = 101325 # Region II - rho1_2 = 1.225 - delta + rho1_2 = RealT(1.225) - delta rho2_2 = delta - v1_2 = zero(delta) - v2_2 = zero(delta) + v1_2 = zero(RealT) + v2_2 = zero(RealT) p_2 = 101325 # Region III - rho1_3 = 1.6861 - delta + rho1_3 = RealT(1.6861) - delta rho2_3 = delta - v1_3 = -113.5243 - v2_3 = zero(delta) + v1_3 = -RealT(113.5243) + v2_3 = zero(RealT) p_3 = 159060 # Set up Region I & II: - inicenter = SVector(zero(delta), zero(delta)) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - if (x[1] > 0.50) + if (x[1] > 0.5f0) # Set up Region III rho1 = rho1_3 rho2 = rho2_3 v1 = v1_3 v2 = v2_3 p = p_3 - elseif (r < 0.25) + elseif (r < 0.25f0) # Set up Region I rho1 = rho1_1 rho2 = rho2_1 diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl index be14c448e4d..ffdfdaec8b5 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl @@ -22,46 +22,47 @@ function initial_condition_shock_bubble(x, t, # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972 # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf # typical domain is rectangular, we change it to a square, as Trixi can only do squares + RealT = eltype(x) @unpack gas_constants = equations # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving - delta = 0.03 + delta = convert(RealT, 0.03) # Region I rho1_1 = delta - rho2_1 = 1.225 * gas_constants[1] / gas_constants[2] - delta - v1_1 = zero(delta) - v2_1 = zero(delta) + rho2_1 = RealT(1.225) * gas_constants[1] / gas_constants[2] - delta + v1_1 = zero(RealT) + v2_1 = zero(RealT) p_1 = 101325 # Region II - rho1_2 = 1.225 - delta + rho1_2 = RealT(1.225) - delta rho2_2 = delta - v1_2 = zero(delta) - v2_2 = zero(delta) + v1_2 = zero(RealT) + v2_2 = zero(RealT) p_2 = 101325 # Region III - rho1_3 = 1.6861 - delta + rho1_3 = RealT(1.6861) - delta rho2_3 = delta - v1_3 = -113.5243 - v2_3 = zero(delta) + v1_3 = -RealT(113.5243) + v2_3 = zero(RealT) p_3 = 159060 # Set up Region I & II: - inicenter = SVector(zero(delta), zero(delta)) + inicenter = SVector(0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - if (x[1] > 0.50) + if (x[1] > 0.5f0) # Set up Region III rho1 = rho1_3 rho2 = rho2_3 v1 = v1_3 v2 = v2_3 p = p_3 - elseif (r < 0.25) + elseif (r < 0.25f0) # Set up Region I rho1 = rho1_1 rho2 = rho2_1 diff --git a/examples/tree_2d_dgsem/elixir_hypdiff_godunov.jl b/examples/tree_2d_dgsem/elixir_hypdiff_godunov.jl index a9b8334c324..6260c3a5c6a 100644 --- a/examples/tree_2d_dgsem/elixir_hypdiff_godunov.jl +++ b/examples/tree_2d_dgsem/elixir_hypdiff_godunov.jl @@ -10,14 +10,15 @@ equations = HyperbolicDiffusionEquations2D() function initial_condition_poisson_periodic(x, t, equations::HyperbolicDiffusionEquations2D) # elliptic equation: -νΔϕ = f # depending on initial constant state, c, for phi this converges to the solution ϕ + c + RealT = eltype(x) if iszero(t) - phi = 0.0 - q1 = 0.0 - q2 = 0.0 + phi = zero(RealT) + q1 = zero(RealT) + q2 = zero(RealT) else - phi = sin(2.0 * pi * x[1]) * sin(2.0 * pi * x[2]) - q1 = 2 * pi * cos(2.0 * pi * x[1]) * sin(2.0 * pi * x[2]) - q2 = 2 * pi * sin(2.0 * pi * x[1]) * cos(2.0 * pi * x[2]) + phi = sinpi(2 * x[1]) * sinpi(2 * x[2]) + q1 = 2 * convert(RealT, pi) * cospi(2 * x[1]) * sinpi(2 * x[2]) + q2 = 2 * convert(RealT, pi) * sinpi(2 * x[1]) * cospi(2 * x[2]) end return SVector(phi, q1, q2) end @@ -27,8 +28,9 @@ initial_condition = initial_condition_poisson_periodic equations::HyperbolicDiffusionEquations2D) # elliptic equation: -νΔϕ = f # analytical solution: phi = sin(2πx)*sin(2πy) and f = -8νπ^2 sin(2πx)*sin(2πy) + RealT = eltype(u) @unpack inv_Tr = equations - C = -8 * equations.nu * pi^2 + C = -8 * equations.nu * convert(RealT, pi)^2 x1, x2 = x tmp1 = sinpi(2 * x1) diff --git a/examples/tree_2d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl index 65bc02187fc..71041044a4f 100644 --- a/examples/tree_2d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl @@ -10,21 +10,22 @@ equations = HyperbolicDiffusionEquations2D() @inline function initial_condition_harmonic_nonperiodic(x, t, equations::HyperbolicDiffusionEquations2D) # elliptic equation: -ν Δϕ = 0 in Ω, u = g on ∂Ω - if t == 0.0 - phi = 1.0 - q1 = 1.0 - q2 = 1.0 + RealT = eltype(x) + if t == 0 + phi = one(RealT) + q1 = one(RealT) + q2 = one(RealT) else - C = inv(sinh(pi)) - sinpi_x1, cospi_x1 = sincos(pi * x[1]) - sinpi_x2, cospi_x2 = sincos(pi * x[2]) - sinh_pix1 = sinh(pi * x[1]) - cosh_pix1 = cosh(pi * x[1]) - sinh_pix2 = sinh(pi * x[2]) - cosh_pix2 = cosh(pi * x[2]) + C = inv(sinh(convert(RealT, pi))) + sinpi_x1, cospi_x1 = sincospi(x[1]) + sinpi_x2, cospi_x2 = sincospi(x[2]) + sinh_pix1 = sinh(convert(RealT, pi) * x[1]) + cosh_pix1 = cosh(convert(RealT, pi) * x[1]) + sinh_pix2 = sinh(convert(RealT, pi) * x[2]) + cosh_pix2 = cosh(convert(RealT, pi) * x[2]) phi = C * (sinh_pix1 * sinpi_x2 + sinh_pix2 * sinpi_x1) - q1 = C * pi * (cosh_pix1 * sinpi_x2 + sinh_pix2 * cospi_x1) - q2 = C * pi * (sinh_pix1 * cospi_x2 + cosh_pix2 * sinpi_x1) + q1 = C * convert(RealT, pi) * (cosh_pix1 * sinpi_x2 + sinh_pix2 * cospi_x1) + q2 = C * convert(RealT, pi) * (sinh_pix1 * cospi_x2 + cosh_pix2 * sinpi_x1) end return SVector(phi, q1, q2) end diff --git a/examples/tree_2d_dgsem/elixir_hypdiff_lax_friedrichs.jl b/examples/tree_2d_dgsem/elixir_hypdiff_lax_friedrichs.jl index c8680dec1ed..901fe173bc8 100644 --- a/examples/tree_2d_dgsem/elixir_hypdiff_lax_friedrichs.jl +++ b/examples/tree_2d_dgsem/elixir_hypdiff_lax_friedrichs.jl @@ -10,14 +10,15 @@ equations = HyperbolicDiffusionEquations2D() function initial_condition_poisson_periodic(x, t, equations::HyperbolicDiffusionEquations2D) # elliptic equation: -νΔϕ = f # depending on initial constant state, c, for phi this converges to the solution ϕ + c + RealT = eltype(x) if iszero(t) - phi = 0.0 - q1 = 0.0 - q2 = 0.0 + phi = zero(RealT) + q1 = zero(RealT) + q2 = zero(RealT) else - phi = sin(2.0 * pi * x[1]) * sin(2.0 * pi * x[2]) - q1 = 2 * pi * cos(2.0 * pi * x[1]) * sin(2.0 * pi * x[2]) - q2 = 2 * pi * sin(2.0 * pi * x[1]) * cos(2.0 * pi * x[2]) + phi = sinpi(2 * x[1]) * sinpi(2 * x[2]) + q1 = 2 * convert(RealT, pi) * cospi(2 * x[1]) * sinpi(2 * x[2]) + q2 = 2 * convert(RealT, pi) * sinpi(2 * x[1]) * cospi(2 * x[2]) end return SVector(phi, q1, q2) end @@ -27,8 +28,9 @@ initial_condition = initial_condition_poisson_periodic equations::HyperbolicDiffusionEquations2D) # elliptic equation: -νΔϕ = f # analytical solution: phi = sin(2πx)*sin(2πy) and f = -8νπ^2 sin(2πx)*sin(2πy) + RealT = eltype(u) @unpack inv_Tr = equations - C = -8 * equations.nu * pi^2 + C = -8 * equations.nu * convert(RealT, pi)^2 x1, x2 = x tmp1 = sinpi(2 * x1) diff --git a/examples/tree_2d_dgsem/elixir_kpp.jl b/examples/tree_2d_dgsem/elixir_kpp.jl index b48bde0d155..8f155c82f90 100644 --- a/examples/tree_2d_dgsem/elixir_kpp.jl +++ b/examples/tree_2d_dgsem/elixir_kpp.jl @@ -24,11 +24,12 @@ end # Since the KPP problem is a scalar equation, the entropy-conservative flux is uniquely determined @inline function Trixi.flux_ec(u_ll, u_rr, orientation::Integer, ::KPPEquation2D) # The tolerance of 1e-12 is based on experience and somewhat arbitrarily chosen - if abs(u_ll[1] - u_rr[1]) < 1e-12 - return 0.5 * (flux(u_ll, orientation, KPPEquation2D()) + + RealT = eltype(u_ll) + if abs(u_ll[1] - u_rr[1]) < RealT(1e-12) + return 0.5f0 * (flux(u_ll, orientation, KPPEquation2D()) + flux(u_rr, orientation, KPPEquation2D())) else - factor = 1.0 / (u_rr[1] - u_ll[1]) + factor = 1 / (u_rr[1] - u_ll[1]) if orientation == 1 return SVector(factor * (-cos(u_rr[1]) + cos(u_ll[1]))) else @@ -38,6 +39,8 @@ end end # Wavespeeds +# Please note that the return type should be modified if other +# floating point types should be used. @inline wavespeed(::KPPEquation2D) = 1.0 @inline Trixi.max_abs_speeds(u, equation::KPPEquation2D) = (wavespeed(equation), wavespeed(equation)) @@ -46,7 +49,7 @@ end norm(normal_direction) # Compute entropy: we use the square entropy -@inline Trixi.entropy(u::Real, ::KPPEquation2D) = 0.5 * u^2 +@inline Trixi.entropy(u::Real, ::KPPEquation2D) = 0.5f0 * u^2 @inline Trixi.entropy(u, ::KPPEquation2D) = entropy(u[1], equation) # Convert between conservative, primitive, and entropy variables. The conserved quantity "u" is also @@ -59,10 +62,11 @@ Trixi.varnames(::Any, ::KPPEquation2D) = ("u",) # Standard KPP test problem with discontinuous initial condition function initial_condition_kpp(x, t, ::KPPEquation2D) + RealT = eltype(x) if x[1]^2 + x[2]^2 < 1 - return SVector(0.25 * 14.0 * pi) + return SVector(0.25f0 * 14 * convert(RealT, pi)) else - return SVector(0.25 * pi) + return SVector(0.25f0 * convert(RealT, pi)) end end diff --git a/examples/tree_2d_dgsem/elixir_lbm_couette.jl b/examples/tree_2d_dgsem/elixir_lbm_couette.jl index 1ba040405d1..5c651bad15f 100644 --- a/examples/tree_2d_dgsem/elixir_lbm_couette.jl +++ b/examples/tree_2d_dgsem/elixir_lbm_couette.jl @@ -16,12 +16,13 @@ incompressible Navier-Stokes equations. To be used in combination with this setup will converge to the state set in [`initial_condition_couette_steady`](@ref). """ function initial_condition_couette_unsteady(x, t, equations::LatticeBoltzmannEquations2D) + RealT = eltype(x) @unpack L, u0, rho0, nu = equations x1, x2 = x v1 = u0 * x2 / L for m in 1:100 - lambda_m = m * pi / L + lambda_m = m * convert(RealT, pi) / L v1 += 2 * u0 * (-1)^m / (lambda_m * L) * exp(-nu * lambda_m^2 * t) * sin(lambda_m * x2) end diff --git a/examples/tree_2d_dgsem/elixir_mhd_blast_wave.jl b/examples/tree_2d_dgsem/elixir_mhd_blast_wave.jl index a0909ca7580..e016115b0c3 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_blast_wave.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_blast_wave.jl @@ -18,23 +18,24 @@ An MHD blast wave taken from function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) # setup taken from Derigs et al. DMV article (2018) # domain must be [-0.5, 0.5] x [-0.5, 0.5], γ = 1.4 + RealT = eltype(x) r = sqrt(x[1]^2 + x[2]^2) - f = (0.1 - r) / 0.01 - if r <= 0.09 - p = 1000.0 - elseif r >= 0.1 - p = 0.1 + f = (convert(RealT, 0.1) - r) / convert(RealT, 0.01) + if r <= RealT(0.09) + p = convert(RealT, 1000) + elseif r >= RealT(0.1) + p = convert(RealT, 0.1) else - p = 0.1 + 999.9 * f + p = convert(RealT, 0.1) + convert(RealT, 999.9) * f end - rho = 1.0 - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - B1 = 100.0 / sqrt(4.0 * pi) - B2 = 0.0 - B3 = 0.0 - psi = 0.0 + rho = 1 + v1 = 0 + v2 = 0 + v3 = 0 + B1 = 100 / sqrt(4 * convert(RealT, pi)) + B2 = 0 + B3 = 0 + psi = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end initial_condition = initial_condition_blast_wave diff --git a/examples/tree_2d_dgsem/elixir_mhd_orszag_tang.jl b/examples/tree_2d_dgsem/elixir_mhd_orszag_tang.jl index 7f26f270d6e..a8e5ce9828e 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_orszag_tang.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_orszag_tang.jl @@ -18,15 +18,15 @@ The classical Orszag-Tang vortex test case. Here, the setup is taken from function initial_condition_orszag_tang(x, t, equations::IdealGlmMhdEquations2D) # setup taken from Derigs et al. DMV article (2018) # domain must be [0, 1] x [0, 1], γ = 5/3 - rho = 1.0 - v1 = -sin(2.0 * pi * x[2]) - v2 = sin(2.0 * pi * x[1]) - v3 = 0.0 - p = 1.0 / equations.gamma - B1 = -sin(2.0 * pi * x[2]) / equations.gamma - B2 = sin(4.0 * pi * x[1]) / equations.gamma - B3 = 0.0 - psi = 0.0 + rho = 1 + v1 = -sinpi(2 * x[2]) + v2 = sinpi(2 * x[1]) + v3 = 0 + p = 1 / equations.gamma + B1 = -sinpi(2 * x[2]) / equations.gamma + B2 = sinpi(4 * x[1]) / equations.gamma + B3 = 0 + psi = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end initial_condition = initial_condition_orszag_tang diff --git a/examples/tree_2d_dgsem/elixir_mhd_rotor.jl b/examples/tree_2d_dgsem/elixir_mhd_rotor.jl index 3109b1ce303..7a54a71ff26 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_rotor.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_rotor.jl @@ -17,29 +17,30 @@ The classical MHD rotor test case. Here, the setup is taken from function initial_condition_rotor(x, t, equations::IdealGlmMhdEquations2D) # setup taken from Derigs et al. DMV article (2018) # domain must be [0, 1] x [0, 1], γ = 1.4 - dx = x[1] - 0.5 - dy = x[2] - 0.5 + RealT = eltype(x) + dx = x[1] - 0.5f0 + dy = x[2] - 0.5f0 r = sqrt(dx^2 + dy^2) - f = (0.115 - r) / 0.015 - if r <= 0.1 - rho = 10.0 - v1 = -20.0 * dy - v2 = 20.0 * dx - elseif r >= 0.115 - rho = 1.0 - v1 = 0.0 - v2 = 0.0 + f = (convert(RealT, 0.115) - r) / convert(RealT, 0.015) + if r <= RealT(0.1) + rho = convert(RealT, 10) + v1 = -20 * dy + v2 = 20 * dx + elseif r >= RealT(0.115) + rho = one(RealT) + v1 = zero(RealT) + v2 = zero(RealT) else - rho = 1.0 + 9.0 * f - v1 = -20.0 * f * dy - v2 = 20.0 * f * dx + rho = 1 + 9 * f + v1 = -20 * f * dy + v2 = 20 * f * dx end - v3 = 0.0 - p = 1.0 - B1 = 5.0 / sqrt(4.0 * pi) - B2 = 0.0 - B3 = 0.0 - psi = 0.0 + v3 = 0 + p = 1 + B1 = 5 / sqrt(4 * convert(RealT, pi)) + B2 = 0 + B3 = 0 + psi = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end initial_condition = initial_condition_rotor diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl index 74d0370647a..6d878f6e7ff 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -19,29 +19,30 @@ This setup needs a positivity limiter for the density. function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) # setup taken from Derigs et al. DMV article (2018) # domain must be [-0.5, 0.5] x [-0.5, 0.5], γ = 1.4 + RealT = eltype(x) r = sqrt(x[1]^2 + x[2]^2) - pmax = 10.0 - pmin = 0.01 - rhomax = 1.0 - rhomin = 0.01 - if r <= 0.09 + pmax = convert(RealT, 10) + pmin = convert(RealT, 0.01) + rhomax = one(RealT) + rhomin = convert(RealT, 0.01) + if r <= RealT(0.09) p = pmax rho = rhomax - elseif r >= 0.1 + elseif r >= RealT(0.1) p = pmin rho = rhomin else - p = pmin + (0.1 - r) * (pmax - pmin) / 0.01 - rho = rhomin + (0.1 - r) * (rhomax - rhomin) / 0.01 + p = pmin + (convert(RealT, 0.1) - r) * (pmax - pmin) / convert(RealT, 0.01) + rho = rhomin + (convert(RealT, 0.1) - r) * (rhomax - rhomin) / convert(RealT, 0.01) end - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - B1 = 1.0 / sqrt(4.0 * pi) - B2 = 0.0 - B3 = 0.0 - psi = 0.0 + v1 = 0 + v2 = 0 + v3 = 0 + B1 = 1 / sqrt(4 * convert(RealT, pi)) + B2 = 0 + B3 = 0 + psi = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end initial_condition = initial_condition_blast_wave diff --git a/examples/tree_2d_dgsem/elixir_mhdmulti_rotor.jl b/examples/tree_2d_dgsem/elixir_mhdmulti_rotor.jl index 5ca21cc5e9c..4f0c07ddd3f 100644 --- a/examples/tree_2d_dgsem/elixir_mhdmulti_rotor.jl +++ b/examples/tree_2d_dgsem/elixir_mhdmulti_rotor.jl @@ -15,32 +15,33 @@ The classical MHD rotor test case adapted to two density components. function initial_condition_rotor(x, t, equations::IdealGlmMhdMulticomponentEquations2D) # setup taken from Derigs et al. DMV article (2018) # domain must be [0, 1] x [0, 1], γ = 1.4 - dx = x[1] - 0.5 - dy = x[2] - 0.5 + RealT = eltype(x) + dx = x[1] - 0.5f0 + dy = x[2] - 0.5f0 r = sqrt(dx^2 + dy^2) - f = (0.115 - r) / 0.015 - if r <= 0.1 - rho1 = 10.0 - rho2 = 5.0 - v1 = -20.0 * dy - v2 = 20.0 * dx - elseif r >= 0.115 - rho1 = 1.0 - rho2 = 0.5 - v1 = 0.0 - v2 = 0.0 + f = (convert(RealT, 0.115) - r) / convert(RealT, 0.015) + if r <= RealT(0.1) + rho1 = convert(RealT, 10) + rho2 = convert(RealT, 5) + v1 = -20 * dy + v2 = 20 * dx + elseif r >= RealT(0.115) + rho1 = one(RealT) + rho2 = convert(RealT, 0.5) + v1 = zero(RealT) + v2 = zero(RealT) else - rho1 = 1.0 + 9.0 * f - rho2 = 0.5 + 4.5 * f - v1 = -20.0 * f * dy - v2 = 20.0 * f * dx + rho1 = 1 + 9 * f + rho2 = 0.5f0 + 4.5f0 * f + v1 = -20 * f * dy + v2 = 20 * f * dx end - v3 = 0.0 - p = 1.0 - B1 = 5.0 / sqrt(4.0 * pi) - B2 = 0.0 - B3 = 0.0 - psi = 0.0 + v3 = 0 + p = 1 + B1 = 5 / sqrt(4 * convert(RealT, pi)) + B2 = 0 + B3 = 0 + psi = 0 return prim2cons(SVector(v1, v2, v3, p, B1, B2, B3, psi, rho1, rho2), equations) end initial_condition = initial_condition_rotor diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index c76e7290e01..a257514d455 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -31,16 +31,17 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) function initial_condition_navier_stokes_convergence_test(x, t, equations) # Amplitude and shift - A = 0.5 - c = 2.0 + RealT = eltype(x) + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_y = convert(RealT, pi) * x[2] + pi_t = convert(RealT, pi) * t rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0))) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2) * (1 - exp(-A * (x[2] - 1))) * cos(pi_t) v2 = v1 p = rho^2 @@ -48,6 +49,7 @@ function initial_condition_navier_stokes_convergence_test(x, t, equations) end @inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + RealT = eltype(x) y = x[2] # TODO: parabolic @@ -59,36 +61,38 @@ end # Same settings as in `initial_condition` # Amplitude and shift - A = 0.5 - c = 2.0 + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_y = convert(RealT, pi) * x[2] + pi_t = convert(RealT, pi) * t # compute the manufactured solution and all necessary derivatives rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) - rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) - rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) - rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - - v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) - v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + rho_t = -convert(RealT, pi) * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = convert(RealT, pi) * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -convert(RealT, pi) * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -convert(RealT, pi)^2 * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -convert(RealT, pi)^2 * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * cos(pi_t) + v1_t = -convert(RealT, pi) * sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * + sin(pi_t) + v1_x = convert(RealT, pi) * cos(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * cos(pi_t) v1_y = sin(pi_x) * - (A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) - v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_xy = pi * cos(pi_x) * - (A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + (A * log(y + 2) * exp(-A * (y - 1)) + + (1 - exp(-A * (y - 1))) / (y + 2)) * cos(pi_t) + v1_xx = -convert(RealT, pi)^2 * sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * + cos(pi_t) + v1_xy = convert(RealT, pi) * cos(pi_x) * + (A * log(y + 2) * exp(-A * (y - 1)) + + (1 - exp(-A * (y - 1))) / (y + 2)) * cos(pi_t) v1_yy = (sin(pi_x) * - (2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) - - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) - - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + (2 * A * exp(-A * (y - 1)) / (y + 2) - + A * A * log(y + 2) * exp(-A * (y - 1)) - + (1 - exp(-A * (y - 1))) / ((y + 2) * (y + 2))) * cos(pi_t)) v2 = v1 v2_t = v1_t v2_x = v1_x @@ -98,21 +102,21 @@ end v2_yy = v1_yy p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_y = 2.0 * rho * rho_y - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x - p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_y = 2 * rho * rho_y + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x + p_yy = 2 * rho * rho_yy + 2 * rho_y * rho_y # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) - E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + E = p * inv_gamma_minus_one + 0.5f0 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2 * rho * v1 * v1_y # Some convenience constants T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_rho_cubed = 1.0 / (rho^3) + inv_rho_cubed = 1 / (rho^3) # compute the source terms # density equation @@ -120,13 +124,13 @@ end # x-momentum equation du2 = (rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 - + 2.0 * rho * v1 * v1_x + + 2 * rho * v1 * v1_x + rho_y * v1 * v2 + rho * v1_y * v2 + rho * v1 * v2_y - # stress tensor from x-direction - 4.0 / 3.0 * v1_xx * mu_ + - 2.0 / 3.0 * v2_xy * mu_ - + RealT(4) / 3 * v1_xx * mu_ + + RealT(2) / 3 * v2_xy * mu_ - v1_yy * mu_ - v2_xy * mu_) # y-momentum equation @@ -134,42 +138,42 @@ end + rho * v1_x * v2 + rho * v1 * v2_x + rho_y * v2^2 - + 2.0 * rho * v2 * v2_y - + + 2 * rho * v2 * v2_y - # stress tensor from y-direction v1_xy * mu_ - v2_xx * mu_ - - 4.0 / 3.0 * v2_yy * mu_ + - 2.0 / 3.0 * v1_xy * mu_) + RealT(4) / 3 * v2_yy * mu_ + + RealT(2) / 3 * v1_xy * mu_) # total energy equation du4 = (E_t + v1_x * (E + p) + v1 * (E_x + p_x) + v2_y * (E + p) + v2 * (E_y + p_y) - # stress tensor and temperature gradient terms from x-direction - 4.0 / 3.0 * v1_xx * v1 * mu_ + - 2.0 / 3.0 * v2_xy * v1 * mu_ - - 4.0 / 3.0 * v1_x * v1_x * mu_ + - 2.0 / 3.0 * v2_y * v1_x * mu_ - + RealT(4) / 3 * v1_xx * v1 * mu_ + + RealT(2) / 3 * v2_xy * v1 * mu_ - + RealT(4) / 3 * v1_x * v1_x * mu_ + + RealT(2) / 3 * v2_y * v1_x * mu_ - v1_xy * v2 * mu_ - v2_xx * v2 * mu_ - v1_y * v2_x * mu_ - v2_x * v2_x * mu_ - T_const * inv_rho_cubed * (p_xx * rho * rho - - 2.0 * p_x * rho * rho_x + - 2.0 * p * rho_x * rho_x - + 2 * p_x * rho * rho_x + + 2 * p * rho_x * rho_x - p * rho * rho_xx) * mu_ - # stress tensor and temperature gradient terms from y-direction v1_yy * v1 * mu_ - v2_xy * v1 * mu_ - v1_y * v1_y * mu_ - v2_x * v1_y * mu_ - - 4.0 / 3.0 * v2_yy * v2 * mu_ + - 2.0 / 3.0 * v1_xy * v2 * mu_ - - 4.0 / 3.0 * v2_y * v2_y * mu_ + - 2.0 / 3.0 * v1_x * v2_y * mu_ - + RealT(4) / 3 * v2_yy * v2 * mu_ + + RealT(2) / 3 * v1_xy * v2 * mu_ - + RealT(4) / 3 * v2_y * v2_y * mu_ + + RealT(2) / 3 * v1_x * v2_y * mu_ - T_const * inv_rho_cubed * (p_yy * rho * rho - - 2.0 * p_y * rho * rho_y + - 2.0 * p * rho_y * rho_y - + 2 * p_y * rho * rho_y + + 2 * p * rho_y * rho_y - p * rho * rho_yy) * mu_) return SVector(du1, du2, du3, du4) @@ -183,7 +187,7 @@ velocity_bc_top_bottom = NoSlip() do x, t, equations_parabolic # This may be an `SVector` or simply a `Tuple` return (u_cons[2] / u_cons[1], u_cons[3] / u_cons[1]) end -heat_bc_top_bottom = Adiabatic((x, t, equations_parabolic) -> 0.0) +heat_bc_top_bottom = Adiabatic((x, t, equations_parabolic) -> zero(eltype(x))) boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index a923844a6e5..10de811a0ff 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -24,18 +24,21 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 30_000) # set maximum capacity of tree data structure function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D) - Ma = 0.1 - rho = 1.0 - u, v = 0.0, 0.0 - p = 1.0 / (Ma^2 * equations.gamma) + RealT = eltype(x) + Ma = convert(RealT, 0.1) + rho = 1 + u, v = 0, 0 + p = 1 / (Ma^2 * equations.gamma) return prim2cons(SVector(rho, u, v, p), equations) end initial_condition = initial_condition_cavity # BC types -velocity_bc_lid = NoSlip((x, t, equations_parabolic) -> SVector(1.0, 0.0)) -velocity_bc_cavity = NoSlip((x, t, equations_parabolic) -> SVector(0.0, 0.0)) -heat_bc = Adiabatic((x, t, equations_parabolic) -> 0.0) +velocity_bc_lid = NoSlip((x, t, equations_parabolic) -> SVector(one(eltype(x)), + zero(eltype(x)))) +velocity_bc_cavity = NoSlip((x, t, equations_parabolic) -> SVector(zero(eltype(x)), + zero(eltype(x)))) +heat_bc = Adiabatic((x, t, equations_parabolic) -> zero(eltype(x))) boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc) boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index a7492bafb47..9f9608de316 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -22,15 +22,16 @@ Brown and Minion (1995). """ function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D) # Shear layer parameters + RealT = eltype(x) k = 80 - delta = 0.05 - u0 = 1.0 + delta = convert(RealT, 0.05) + u0 = 1 - Ms = 0.1 # maximum Mach number + Ms = convert(RealT, 0.1) # maximum Mach number - rho = 1.0 - v1 = x[2] <= 0.5 ? u0 * tanh(k * (x[2] - 0.25)) : u0 * tanh(k * (0.75 - x[2])) - v2 = u0 * delta * sin(2 * pi * (x[1] + 0.25)) + rho = 1 + v1 = x[2] <= 0.5f0 ? u0 * tanh(k * (x[2] - 0.25f0)) : u0 * tanh(k * (0.75f0 - x[2])) + v2 = u0 * delta * sinpi(2 * (x[1] + 0.25f0)) p = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms return prim2cons(SVector(rho, v1, v2, p), equations) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index c6e5f0bc40a..95c3f3ee252 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -23,14 +23,15 @@ This forms the basis behind the 3D case found for instance in """ function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D) - A = 1.0 # magnitude of speed - Ms = 0.1 # maximum Mach number + RealT = eltype(x) + A = 1 # magnitude of speed + Ms = convert(RealT, 0.1) # maximum Mach number - rho = 1.0 + rho = 1 v1 = A * sin(x[1]) * cos(x[2]) v2 = -A * cos(x[1]) * sin(x[2]) p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms - p = p + 1.0 / 4.0 * A^2 * rho * (cos(2 * x[1]) + cos(2 * x[2])) + p = p + 0.25f0 * A^2 * rho * (cos(2 * x[1]) + cos(2 * x[2])) return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl index 9598ae184cd..f48cff6a710 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl @@ -13,15 +13,16 @@ prandtl_number() = 0.72 # 1991, McGraw-Hill, ISBN, 0-07-069712-4 # Pages 28 and 29. @inline function mu(u, equations) - T_ref = 291.15 + RealT = eltype(u) + T_ref = convert(RealT, 291.15) - R_specific_air = 287.052874 + R_specific_air = convert(RealT, 287.052874) T = R_specific_air * Trixi.temperature(u, equations) - C_air = 120.0 - mu_ref_air = 1.827e-5 + C_air = 120 + mu_ref_air = convert(RealT, 1.827e-5) - return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5 + return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5f0 end equations = CompressibleEulerEquations2D(1.4) @@ -39,14 +40,15 @@ This forms the basis behind the 3D case found for instance in """ function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D) - A = 1.0 # magnitude of speed - Ms = 0.1 # maximum Mach number + RealT = eltype(x) + A = 1 # magnitude of speed + Ms = convert(RealT, 0.1) # maximum Mach number - rho = 1.0 + rho = 1 v1 = A * sin(x[1]) * cos(x[2]) v2 = -A * cos(x[1]) * sin(x[2]) p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms - p = p + 1.0 / 4.0 * A^2 * rho * (cos(2 * x[1]) + cos(2 * x[2])) + p = p + 0.25f0 * A^2 * rho * (cos(2 * x[1]) + cos(2 * x[2])) return prim2cons(SVector(rho, v1, v2, p), equations) end diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl index 8221dfebe39..97e9065dc50 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl @@ -50,7 +50,8 @@ ode = semidiscretize(semi, tspan) function initial_condition_ec_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations2D) # Set up polar coordinates - inicenter = SVector(0.7, 0.7) + RealT = eltype(x) + inicenter = SVector(RealT(0.7), RealT(0.7)) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -58,21 +59,21 @@ function initial_condition_ec_discontinuous_bottom(x, t, element_id, sin_phi, cos_phi = sincos(phi) # Set the background values - H = 4.25 - v1 = 0.0 - v2 = 0.0 - b = 0.0 + H = 4.25f0 + v1 = zero(RealT) + v2 = zero(RealT) + b = zero(RealT) # setup the discontinuous water height and velocities if element_id == 10 - H = 5.0 - v1 = 0.1882 * cos_phi - v2 = 0.1882 * sin_phi + H = 5.0f0 + v1 = RealT(0.1882) * cos_phi + v2 = RealT(0.1882) * sin_phi end # Setup a discontinuous bottom topography using the element id number if element_id == 7 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[2]) + b = 2 + 0.5f0 * sinpi(2 * x[1]) + 0.5f0 * cospi(2 * x[2]) end return prim2cons(SVector(H, v1, v2, b), equations) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl index f8f601d4120..96beade8669 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -11,13 +11,13 @@ equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) function initial_condition_perturbation(x, t, equations::ShallowWaterEquations2D) # Set the background values H = equations.H0 - v1 = 0.0 - v2 = 0.0 + v1 = 0 + v2 = 0 # Bottom topography - b = 1.5 * exp(-0.5 * ((x[1])^2 + (x[2])^2)) + b = 1.5f0 * exp(-0.5f0 * ((x[1])^2 + (x[2])^2)) # Waterheight perturbation - H = H + 0.5 * exp(-10.0 * ((x[1])^2 + (x[2])^2)) + H = H + 0.5f0 * exp(-10 * ((x[1])^2 + (x[2])^2)) return prim2cons(SVector(H, v1, v2, b), equations) end diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl index 22043392b2a..50948aa44b2 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -14,12 +14,12 @@ equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D) # Set the background values H = equations.H0 - v1 = 0.0 - v2 = 0.0 + v1 = 0 + v2 = 0 # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) x1, x2 = x - b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + - 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + b = (1.5f0 / exp(0.5f0 * ((x1 - 1)^2 + (x2 - 1)^2)) + + 0.75f0 / exp(0.5f0 * ((x1 + 1)^2 + (x2 + 1)^2))) return prim2cons(SVector(H, v1, v2, b), equations) end @@ -64,14 +64,15 @@ ode = semidiscretize(semi, tspan) function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations2D) # Set the background values + RealT = eltype(x) H = equations.H0 - v1 = 0.0 - v2 = 0.0 - b = 0.0 + v1 = 0 + v2 = 0 + b = zero(RealT) # Setup a discontinuous bottom topography using the element id number if element_id == 7 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[2]) + b = 2 + 0.5f0 * sinpi(2 * x[1]) + 0.5f0 * cospi(2 * x[2]) end return prim2cons(SVector(H, v1, v2, b), equations) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl index 19073b0504a..bc1200db6fe 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl @@ -13,12 +13,12 @@ equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D) # Set the background values H = equations.H0 - v1 = 0.0 - v2 = 0.0 + v1 = 0 + v2 = 0 # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) x1, x2 = x - b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + - 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + b = (1.5f0 / exp(0.5f0 * ((x1 - 1)^2 + (x2 - 1)^2)) + + 0.75f0 / exp(0.5f0 * ((x1 + 1)^2 + (x2 + 1)^2))) return prim2cons(SVector(H, v1, v2, b), equations) end @@ -67,14 +67,15 @@ ode = semidiscretize(semi, tspan) function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations2D) # Set the background values + RealT = eltype(x) H = equations.H0 - v1 = 0.0 - v2 = 0.0 - b = 0.0 + v1 = 0 + v2 = 0 + b = zero(RealT) # Setup a discontinuous bottom topography using the element id number if element_id == 7 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[2]) + b = 2 + 0.5f0 * sinpi(2 * x[1]) + 0.5f0 * cospi(2 * x[2]) end return prim2cons(SVector(H, v1, v2, b), equations)