From c919da0d86af29053473beb829569ba2ce515203 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 18 Oct 2024 20:52:35 -1000 Subject: [PATCH 1/8] Start --- .../elixir_euler_blast_wave_entropy_bounded.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl index fa81defb6ca..1bd589a368f 100644 --- a/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl +++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl @@ -18,7 +18,8 @@ A medium blast wave taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # The following code is equivalent to @@ -28,9 +29,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_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 + p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245) return prim2cons(SVector(rho, v1, p), equations) end From c2aa97ce2e7f8448cf813c740f25a826aead8b78 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Fri, 25 Oct 2024 23:06:48 -1000 Subject: [PATCH 2/8] Start again --- examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl index 843d8a7aa4c..d73c08d03b1 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl @@ -35,7 +35,7 @@ mesh = TreeMesh(coordinate_min, coordinate_max, # Discontinuous initial condition (Riemann Problem) leading to a rarefaction fan. function initial_condition_rarefaction(x, t, equation::InviscidBurgersEquation1D) - scalar = x[1] < 0.5 ? 0.5 : 1.5 + scalar = x[1] < 0.5f0 ? 0.5f0 : 1.5f0 return SVector(scalar) end From be916fb588fd79034b9a124a9251c57afd8fc31a Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 5 Nov 2024 13:56:41 -1000 Subject: [PATCH 3/8] Add more --- .../tree_1d_dgsem/elixir_advection_diffusion.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl index b75e3622ac1..b3144d0940d 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl @@ -22,10 +22,12 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 30_000, # set maximum capacity of tree data structure periodicity = true) -function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0)) +function x_trans_periodic(x, domain_length = SVector(oftype(x, 2 * pi)), + center = SVector(oftype(x, 0))) x_normalized = x .- center x_shifted = x_normalized .% domain_length - x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* + x_offset = ((x_shifted .< -0.5f0 * domain_length) - + (x_shifted .> 0.5f0 * domain_length)) .* domain_length return center + x_shifted + x_offset end @@ -37,11 +39,11 @@ function initial_condition_diffusive_convergence_test(x, t, x_trans = x_trans_periodic(x - equation.advection_velocity * t) nu = diffusivity() - c = 0.0 - A = 1.0 + c = 0 + A = 1 L = 2 - f = 1 / L - omega = 1.0 + f = 1.0f0 / L + omega = 1 scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t) return SVector(scalar) end From 9c94c645b1ff4e96c9653689478ae5f21e5ec19c Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 5 Nov 2024 17:45:30 -1000 Subject: [PATCH 4/8] Complete --- .../elixir_burgers_linear_stability.jl | 3 +- .../tree_1d_dgsem/elixir_burgers_shock.jl | 2 +- .../tree_1d_dgsem/elixir_euler_blast_wave.jl | 9 +-- .../tree_1d_dgsem/elixir_euler_positivity.jl | 15 ++--- .../elixir_euler_quasi_1d_discontinuous.jl | 9 +-- .../tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl | 9 +-- .../elixir_euler_sedov_blast_wave.jl | 15 ++--- .../elixir_euler_sedov_blast_wave_pure_fv.jl | 15 ++--- ..._eulermulti_two_interacting_blast_waves.jl | 13 ++-- .../elixir_hypdiff_harmonic_nonperiodic.jl | 9 +-- .../elixir_linearizedeuler_gauss_wall.jl | 2 +- .../elixir_maxwell_E_excitation.jl | 5 +- .../elixir_mhd_briowu_shock_tube.jl | 17 ++--- .../elixir_mhd_ryujones_shock_tube.jl | 15 ++--- .../elixir_mhd_shu_osher_shock_tube.jl | 38 ++++++------ .../elixir_mhd_torrilhon_shock_tube.jl | 17 ++--- .../elixir_mhdmulti_briowu_shock_tube.jl | 35 +++++------ ...lixir_navierstokes_convergence_periodic.jl | 50 ++++++++------- .../elixir_navierstokes_convergence_walls.jl | 62 ++++++++++--------- ...ixir_navierstokes_convergence_walls_amr.jl | 62 ++++++++++--------- .../tree_1d_dgsem/elixir_shallowwater_ec.jl | 15 ++--- ...xir_shallowwater_quasi_1d_discontinuous.jl | 13 ++-- ...xir_shallowwater_quasi_1d_well_balanced.jl | 4 +- .../elixir_shallowwater_shock_capturing.jl | 18 +++--- .../elixir_shallowwater_well_balanced.jl | 8 +-- ..._shallowwater_well_balanced_nonperiodic.jl | 4 +- .../elixir_traffic_flow_lwr_trafficjam.jl | 2 +- 27 files changed, 244 insertions(+), 222 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl b/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl index ae2039edde6..438356144cc 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl @@ -8,8 +8,9 @@ using Trixi equations = InviscidBurgersEquation1D() function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D) + RealT = eltype(x) k = 1 - 2 + sinpi(k * (x[1] - 0.7)) |> SVector + SVector(2 + sinpi(k * (x[1] - convert(RealT, 0.7)))) end volume_flux = flux_ec diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl index 09c6cee9f98..89939e4328a 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl @@ -35,7 +35,7 @@ mesh = TreeMesh(coordinate_min, coordinate_max, # Discontinuous initial condition (Riemann Problem) leading to a shock to test e.g. correct shock speed. function initial_condition_shock(x, t, equation::InviscidBurgersEquation1D) - scalar = x[1] < 0.5 ? 1.5 : 0.5 + scalar = x[1] < 0.5f0 ? 1.5f0 : 0.5f0 return SVector(scalar) end diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl index 9cba4936d22..5b5eb6fdd33 100644 --- a/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl +++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl @@ -18,7 +18,8 @@ A medium blast wave taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # The following code is equivalent to @@ -28,9 +29,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_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 + p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245) return prim2cons(SVector(rho, v1, p), equations) end diff --git a/examples/tree_1d_dgsem/elixir_euler_positivity.jl b/examples/tree_1d_dgsem/elixir_euler_positivity.jl index 3fb96fb807b..34a6ad16a43 100644 --- a/examples/tree_1d_dgsem/elixir_euler_positivity.jl +++ b/examples/tree_1d_dgsem/elixir_euler_positivity.jl @@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # 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 = 6 * (equations.gamma - 1) * E / (3 * pi * r0) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0) + 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 + rho = 1 + v1 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, p), equations) diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl index cc4535be028..0af5b210ee9 100644 --- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl +++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl @@ -18,10 +18,11 @@ A discontinuous initial condition taken from """ function initial_condition_discontinuity(x, t, equations::CompressibleEulerEquationsQuasi1D) - rho = (x[1] < 0) ? 3.4718 : 2.0 - v1 = (x[1] < 0) ? -2.5923 : -3.0 - p = (x[1] < 0) ? 5.7118 : 2.639 - a = (x[1] < 0) ? 1.0 : 1.5 + RealT = eltype(x) + rho = (x[1] < 0) ? RealT(3.4718) : RealT(2.0) + v1 = (x[1] < 0) ? RealT(-2.5923) : RealT(-3.0) + p = (x[1] < 0) ? RealT(5.7118) : RealT(2.639) + a = (x[1] < 0) ? 1.0f0 : 1.5f0 return prim2cons(SVector(rho, v1, p, a), equations) end diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl index ae1b2b24b62..ebd3680cef4 100644 --- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl +++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl @@ -15,10 +15,11 @@ equations = CompressibleEulerEquationsQuasi1D(1.4) # refinement level is changed the initial condition below may need changed as well to # ensure that the discontinuities lie on an element interface. function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D) - v1 = 0.1 - rho = 2.0 + 0.1 * x[1] - p = 3.0 - a = 2.0 + x[1] + RealT = eltype(x) + v1 = convert(RealT, 0.1) + rho = 2 + convert(RealT, 0.1) * x[1] + p = 3 + a = 2 + x[1] return prim2cons(SVector(rho, v1, p, a), equations) end diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl index b67b2bc602e..a89727da147 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # 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 = 6 * (equations.gamma - 1) * E / (3 * pi * r0) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0) + 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 + rho = 1 + v1 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, p), equations) diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl index 8a0241680b9..cf33f7f5b64 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl @@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # 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 = 6 * (equations.gamma - 1) * E / (3 * pi * r0) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0) + 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 + rho = 1 + v1 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, p), equations) diff --git a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl index c093420cc7c..908554e0e39 100644 --- a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl +++ b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl @@ -18,18 +18,19 @@ A multicomponent two interacting blast wave test taken from """ function initial_condition_two_interacting_blast_waves(x, t, equations::CompressibleEulerMulticomponentEquations1D) - rho1 = 0.5 * x[1]^2 - rho2 = 0.5 * (sin(20 * x[1]))^2 + RealT = eltype(x) + rho1 = 0.5f0 * x[1]^2 + rho2 = 0.5f0 * (sin(20 * x[1]))^2 rho3 = 1 - rho1 - rho2 prim_rho = SVector{3, real(equations)}(rho1, rho2, rho3) - v1 = 0.0 + v1 = 0 - if x[1] <= 0.1 + if x[1] <= RealT(0.1) p = 1000 - elseif x[1] < 0.9 - p = 0.01 + elseif x[1] < RealT(0.9) + p = convert(RealT, 0.01) else p = 100 end diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl index b9173ec9f49..b0f68c73689 100644 --- a/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl @@ -19,12 +19,13 @@ A non-priodic harmonic function used in combination with function initial_condition_harmonic_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D) # elliptic equation: -νΔϕ = f - if t == 0.0 - phi = 5.0 - q1 = 0.0 + RealT = eltype(x) + if t == 0 + phi = 5 + q1 = 0 else A = 3 - B = exp(1) + B = exp(one(RealT)) phi = A + B * x[1] q1 = B end diff --git a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl index 0884249559a..7a7b83dfe84 100644 --- a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -22,7 +22,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # that is advected to left with v - c and to the right with v + c. # Correspondigly, the bump splits in half. function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations1D) - v1_prime = 0.0 + v1_prime = 0 rho_prime = p_prime = 2 * exp(-(x[1] - 45)^2 / 25) return SVector(rho_prime, v1_prime, p_prime) end diff --git a/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl b/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl index b5478331a7d..da44be30da2 100644 --- a/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl +++ b/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl @@ -20,9 +20,10 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # Excite the electric field which causes a standing wave. # The solution is an undamped exchange between electric and magnetic energy. function initial_condition_E_excitation(x, t, equations::MaxwellEquations1D) + RealT = eltype(x) c = equations.speed_of_light - E = -c * sin(2 * pi * x[1]) - B = 0.0 + E = -c * sinpi(2 * x[1]) + B = 0 return SVector(E, B) end diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl index 4d537508b47..8f33ad77cfc 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl @@ -18,14 +18,15 @@ MHD extension of the Sod shock tube. Taken from Section V of the article """ function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [0, 1], γ = 2, final time = 0.12 - rho = x[1] < 0.5 ? 1.0 : 0.125 - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - p = x[1] < 0.5 ? 1.0 : 0.1 - B1 = 0.75 - B2 = x[1] < 0.5 ? 1.0 : -1.0 - B3 = 0.0 + RealT = eltype(x) + rho = x[1] < 0.5f0 ? 1.0f0 : 0.125f0 + v1 = 0 + v2 = 0 + v3 = 0 + p = x[1] < 0.5f0 ? one(RealT) : RealT(0.1) + B1 = 0.75f0 + B2 = x[1] < 0.5f0 ? 1 : -1 + B3 = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end initial_condition = initial_condition_briowu_shock_tube diff --git a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl index a7d7689a806..fc25ef87df7 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl @@ -23,14 +23,15 @@ present in the one dimensional MHD equations. It is the second test from Section """ function initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [0, 1], γ = 5/3, final time = 0.2 - rho = x[1] <= 0.5 ? 1.08 : 1.0 - v1 = x[1] <= 0.5 ? 1.2 : 0.0 - v2 = x[1] <= 0.5 ? 0.01 : 0.0 - v3 = x[1] <= 0.5 ? 0.5 : 0.0 - p = x[1] <= 0.5 ? 0.95 : 1.0 - inv_sqrt4pi = 1.0 / sqrt(4 * pi) + RealT = eltype(x) + rho = x[1] <= 0.5f0 ? RealT(1.08) : one(RealT) + v1 = x[1] <= 0.5f0 ? RealT(1.2) : zero(RealT) + v2 = x[1] <= 0.5f0 ? RealT(0.01) : zero(RealT) + v3 = x[1] <= 0.5f0 ? 0.5f0 : 0.0f0 + p = x[1] <= 0.5f0 ? RealT(0.95) : one(RealT) + inv_sqrt4pi = 1 / sqrt(4 * convert(RealT, pi)) B1 = 2 * inv_sqrt4pi - B2 = x[1] <= 0.5 ? 3.6 * inv_sqrt4pi : 4.0 * inv_sqrt4pi + B2 = x[1] <= 0.5f0 ? RealT(3.6) * inv_sqrt4pi : 4 * inv_sqrt4pi B3 = B1 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) diff --git a/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl index 689537ebdd4..1ab4c54985b 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl @@ -19,15 +19,16 @@ Taken from Section 4.1 of function initial_condition_shu_osher_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [-5, 5], γ = 5/3, final time = 0.7 # initial shock location is taken to be at x = -4 - x_0 = -4.0 - rho = x[1] <= x_0 ? 3.5 : 1.0 + 0.2 * sin(5.0 * x[1]) - v1 = x[1] <= x_0 ? 5.8846 : 0.0 - v2 = x[1] <= x_0 ? 1.1198 : 0.0 - v3 = 0.0 - p = x[1] <= x_0 ? 42.0267 : 1.0 - B1 = 1.0 - B2 = x[1] <= x_0 ? 3.6359 : 1.0 - B3 = 0.0 + RealT = eltype(x) + x_0 = -4 + rho = x[1] <= x_0 ? RealT(3.5) : 1 + RealT(0.2) * sin(5 * x[1]) + v1 = x[1] <= x_0 ? RealT(5.8846) : zero(RealT) + v2 = x[1] <= x_0 ? RealT(1.1198) : zero(RealT) + v3 = 0 + p = x[1] <= x_0 ? RealT(42.0267) : one(RealT) + B1 = 1 + B2 = x[1] <= x_0 ? RealT(3.6359) : one(RealT) + B3 = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end @@ -45,15 +46,16 @@ function initial_condition_shu_osher_shock_tube_flipped(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [-5, 5], γ = 5/3, final time = 0.7 # initial shock location is taken to be at x = 4 - x_0 = 4.0 - rho = x[1] <= x_0 ? 1.0 + 0.2 * sin(5.0 * x[1]) : 3.5 - v1 = x[1] <= x_0 ? 0.0 : -5.8846 - v2 = x[1] <= x_0 ? 0.0 : -1.1198 - v3 = 0.0 - p = x[1] <= x_0 ? 1.0 : 42.0267 - B1 = 1.0 - B2 = x[1] <= x_0 ? 1.0 : 3.6359 - B3 = 0.0 + RealT = eltype(x) + x_0 = 4 + rho = x[1] <= x_0 ? 1 + RealT(0.2) * sin(5 * x[1]) : RealT(3.5) + v1 = x[1] <= x_0 ? zero(RealT) : RealT(-5.8846) + v2 = x[1] <= x_0 ? zero(RealT) : RealT(-1.1198) + v3 = 0 + p = x[1] <= x_0 ? one(RealT) : RealT(42.0267) + B1 = 1 + B2 = x[1] <= x_0 ? one(RealT) : RealT(3.6359) + B3 = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end diff --git a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl index 3b366c35e0f..0d2a4418b35 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl @@ -17,14 +17,15 @@ Torrilhon's shock tube test case for one dimensional ideal MHD equations. """ function initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [-1, 1.5], γ = 5/3, final time = 0.4 - rho = x[1] <= 0 ? 3.0 : 1.0 - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - p = x[1] <= 0 ? 3.0 : 1.0 - B1 = 1.5 - B2 = x[1] <= 0 ? 1.0 : cos(1.5) - B3 = x[1] <= 0 ? 0.0 : sin(1.5) + RealT = eltype(x) + rho = x[1] <= 0 ? 3 : 1 + v1 = 0 + v2 = 0 + v3 = 0 + p = x[1] <= 0 ? 3 : 1 + B1 = 1.5f0 + B2 = x[1] <= 0 ? zero(RealT) : cos(RealT(1.5)) + B3 = x[1] <= 0 ? zero(RealT) : sin(RealT(1.5)) return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end initial_condition = initial_condition_torrilhon_shock_tube diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl index 831fa7afedb..999e6c42b29 100644 --- a/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl @@ -19,28 +19,27 @@ MHD extension of the Sod shock tube. Taken from Section V of the article function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdMulticomponentEquations1D) # domain must be set to [0, 1], γ = 2, final time = 0.12 - if x[1] < 0.5 - rho = 1.0 - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + RealT = eltype(x) + if x[1] < 0.5f0 + rho = one(RealT) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) else - rho = 0.125 - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + rho = RealT(0.125) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) end - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - p = x[1] < 0.5 ? 1.0 : 0.1 - B1 = 0.75 - B2 = x[1] < 0.5 ? 1.0 : -1.0 - B3 = 0.0 + v1 = 0 + v2 = 0 + v3 = 0 + p = x[1] < 0.5f0 ? one(RealT) : RealT(0.1) + B1 = 0.75f0 + B2 = x[1] < 0.5f0 ? 1 : -1 + B3 = 0 prim_other = SVector(v1, v2, v3, p, B1, B2, B3) return prim2cons(vcat(prim_other, prim_rho), equations) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl index eab0840f385..1089cf4bb85 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl @@ -16,12 +16,13 @@ equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu = mu(), # (Simplified version of the 2D) 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_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t rho = c + A * sin(pi_x) * cos(pi_t) v1 = sin(pi_x) * cos(pi_t) @@ -34,42 +35,43 @@ initial_condition = initial_condition_navier_stokes_convergence_test @inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) # we currently need to hardcode these parameters until we fix the "combined equation" issue # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + RealT = eltype(x) inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() mu_ = mu() # 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_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t # compute the manufactured solution and all necessary derivatives rho = c + A * sin(pi_x) * cos(pi_t) - rho_t = -pi * A * sin(pi_x) * sin(pi_t) - rho_x = pi * A * cos(pi_x) * cos(pi_t) - rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_t) + rho_t = -convert(RealT, pi) * A * sin(pi_x) * sin(pi_t) + rho_x = convert(RealT, pi) * A * cos(pi_x) * cos(pi_t) + rho_xx = -convert(RealT, pi) * convert(RealT, pi) * A * sin(pi_x) * cos(pi_t) v1 = sin(pi_x) * cos(pi_t) - v1_t = -pi * sin(pi_x) * sin(pi_t) - v1_x = pi * cos(pi_x) * cos(pi_t) - v1_xx = -pi * pi * sin(pi_x) * cos(pi_t) + v1_t = -convert(RealT, pi) * sin(pi_x) * sin(pi_t) + v1_x = convert(RealT, pi) * cos(pi_x) * cos(pi_t) + v1_xx = -convert(RealT, pi) * convert(RealT, pi) * sin(pi_x) * cos(pi_t) p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x - E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 - E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + E = p * inv_gamma_minus_one + 0.5f0 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5f0 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5f0 * rho_x * v1^2 + rho * v1 * v1_x # 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 @@ -77,7 +79,7 @@ initial_condition = initial_condition_navier_stokes_convergence_test # x-momentum equation du2 = (rho_t * v1 + rho * v1_t - + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - + + p_x + rho_x * v1^2 + 2 * rho * v1 * v1_x - # stress tensor from x-direction v1_xx * mu_) @@ -88,8 +90,8 @@ initial_condition = initial_condition_navier_stokes_convergence_test v1_x * v1_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_) return SVector(du1, du2, du3) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl index 40030d53345..d534b7b189b 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl @@ -31,21 +31,23 @@ 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_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t rho = c + A * cos(pi_x) * cos(pi_t) - v1 = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0))) * cos(pi_t) + v1 = log(x[1] + 2) * (1 - exp(-A * (x[1] - 1))) * cos(pi_t) p = rho^2 return prim2cons(SVector(rho, v1, p), equations) end @inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + RealT = eltype(x) x = x[1] # TODO: parabolic @@ -57,40 +59,40 @@ 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 - pi_t = pi * t + pi_x = convert(RealT, pi) * x + pi_t = convert(RealT, pi) * t # compute the manufactured solution and all necessary derivatives rho = c + A * cos(pi_x) * cos(pi_t) - rho_t = -pi * A * cos(pi_x) * sin(pi_t) - rho_x = -pi * A * sin(pi_x) * cos(pi_t) - rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t) - - v1 = log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t) - v1_t = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t) - v1_x = (A * log(x + 2.0) * exp(-A * (x - 1.0)) + - (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t) - v1_xx = ((2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0) - - A * A * log(x + 2.0) * exp(-A * (x - 1.0)) - - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t)) + rho_t = -convert(RealT, pi) * A * cos(pi_x) * sin(pi_t) + rho_x = -convert(RealT, pi) * A * sin(pi_x) * cos(pi_t) + rho_xx = -convert(RealT, pi) * convert(RealT, pi) * A * cos(pi_x) * cos(pi_t) + + v1 = log(x + 2) * (1 - exp(-A * (x - 1))) * cos(pi_t) + v1_t = -convert(RealT, pi) * log(x + 2) * (1 - exp(-A * (x - 1))) * sin(pi_t) + v1_x = (A * log(x + 2) * exp(-A * (x - 1)) + + (1 - exp(-A * (x - 1))) / (x + 2)) * cos(pi_t) + v1_xx = ((2 * A * exp(-A * (x - 1)) / (x + 2) - + A * A * log(x + 2) * exp(-A * (x - 1)) - + (1 - exp(-A * (x - 1))) / ((x + 2) * (x + 2))) * cos(pi_t)) p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 - E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + E = p * inv_gamma_minus_one + 0.5f0 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5f0 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5f0 * rho_x * v1^2 + rho * v1 * v1_x # 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 @@ -98,7 +100,7 @@ end # y-momentum equation du2 = (rho_t * v1 + rho * v1_t - + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - + + p_x + rho_x * v1^2 + 2 * rho * v1 * v1_x - # stress tensor from y-direction v1_xx * mu_) @@ -109,8 +111,8 @@ end v1_x * v1_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_) return SVector(du1, du2, du3) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl index e833155a68e..4f9cd6365fa 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl @@ -31,21 +31,23 @@ 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_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t rho = c + A * cos(pi_x) * cos(pi_t) - v1 = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0))) * cos(pi_t) + v1 = log(x[1] + 2) * (1 - exp(-A * (x[1] - 1))) * cos(pi_t) p = rho^2 return prim2cons(SVector(rho, v1, p), equations) end @inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + RealT = eltype(x) x = x[1] # TODO: parabolic @@ -57,40 +59,40 @@ 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 - pi_t = pi * t + pi_x = convert(RealT, pi) * x + pi_t = convert(RealT, pi) * t # compute the manufactured solution and all necessary derivatives rho = c + A * cos(pi_x) * cos(pi_t) - rho_t = -pi * A * cos(pi_x) * sin(pi_t) - rho_x = -pi * A * sin(pi_x) * cos(pi_t) - rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t) - - v1 = log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t) - v1_t = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t) - v1_x = (A * log(x + 2.0) * exp(-A * (x - 1.0)) + - (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t) - v1_xx = ((2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0) - - A * A * log(x + 2.0) * exp(-A * (x - 1.0)) - - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t)) + rho_t = -convert(RealT, pi) * A * cos(pi_x) * sin(pi_t) + rho_x = -convert(RealT, pi) * A * sin(pi_x) * cos(pi_t) + rho_xx = -convert(RealT, pi) * convert(RealT, pi) * A * cos(pi_x) * cos(pi_t) + + v1 = log(x + 2) * (1 - exp(-A * (x - 1))) * cos(pi_t) + v1_t = -convert(RealT, pi) * log(x + 2) * (1 - exp(-A * (x - 1))) * sin(pi_t) + v1_x = (A * log(x + 2) * exp(-A * (x - 1)) + + (1 - exp(-A * (x - 1))) / (x + 2)) * cos(pi_t) + v1_xx = ((2 * A * exp(-A * (x - 1)) / (x + 2) - + A * A * log(x + 2) * exp(-A * (x - 1)) - + (1 - exp(-A * (x - 1))) / ((x + 2) * (x + 2))) * cos(pi_t)) p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 - E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + E = p * inv_gamma_minus_one + 0.5f0 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5f0 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5f0 * rho_x * v1^2 + rho * v1 * v1_x # 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 @@ -98,7 +100,7 @@ end # y-momentum equation du2 = (rho_t * v1 + rho * v1_t - + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - + + p_x + rho_x * v1^2 + 2 * rho * v1 * v1_x - # stress tensor from y-direction v1_xx * mu_) @@ -109,8 +111,8 @@ end v1_x * v1_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_) return SVector(du1, du2, du3) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index af0da5d1768..c3a94a27d83 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -17,19 +17,20 @@ equations = ShallowWaterEquations1D(gravity_constant = 9.81) # ensure that the discontinuities lie on an element interface. function initial_condition_ec_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D) # Set the background values - H = 4.25 - v = 0.0 + RealT = eltype(x) + H = 4.25f0 + v = 0 b = sin(x[1]) # arbitrary continuous function # Setup the discontinuous water height and velocity - if x[1] >= 0.125 && x[1] <= 0.25 - H = 5.0 - v = 0.1882 + if x[1] >= 0.125f0 && x[1] <= 0.25f0 + H = 5 + v = convert(RealT, 0.1882) end # Setup a discontinuous bottom topography - if x[1] >= -0.25 && x[1] <= -0.125 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + if x[1] >= -0.25f0 && x[1] <= -0.125f0 + b = 2 + 0.5f0 * sinpi(2 * x[1]) end return prim2cons(SVector(H, v, b), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl index 76c04759389..15c2f46ae7d 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl @@ -8,15 +8,16 @@ using Trixi equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81) function initial_condition_discontinuity(x, t, equations::ShallowWaterEquationsQuasi1D) - H = 2 + 0.1 * exp(-25 * x[1]^2) - v = 0.0 + RealT = eltype(x) + H = 2 + convert(RealT, 0.1) * exp(-25 * x[1]^2) + v = 0 if x[1] > 0 - b = 0.1 - a = 1.0 + b = convert(RealT, 0.1) + a = 1 else - b = 0.0 - a = 1.1 + b = 0 + a = convert(RealT, 1.1) end return prim2cons(SVector(H, v, b, a), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl index a4f4b0189ba..4a777dc732c 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl @@ -17,11 +17,11 @@ equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81, H0 = 2.0) function initial_condition_discontinuous_well_balancedness(x, t, equations::ShallowWaterEquationsQuasi1D) H = equations.H0 - v = 0.0 + v = 0 # for a periodic domain, this choice of `b` and `a` mimic # discontinuity across the periodic boundary. - b = 0.5 * (x[1] + 1) + b = 0.5f0 * (x[1] + 1) a = 2 + x[1] return prim2cons(SVector(H, v, b, a), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl index 511f33d1101..fba76f6fa99 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl @@ -20,19 +20,19 @@ function initial_condition_stone_throw_discontinuous_bottom(x, t, H = equations.H0 # Discontinuous velocity - v = 0.0 - if x[1] >= -0.75 && x[1] <= 0.0 - v = -1.0 - elseif x[1] >= 0.0 && x[1] <= 0.75 - v = 1.0 + v = 0 + if x[1] >= -0.75f0 && x[1] <= 0 + v = -1 + elseif x[1] >= 0 && x[1] <= 0.75f0 + v = 1 end - b = (1.5 / exp(0.5 * ((x[1] - 1.0)^2)) + - 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + b = (1.5f0 / exp(0.5f0 * ((x[1] - 1)^2)) + + 0.75f0 / exp(0.5f0 * ((x[1] + 1)^2))) # Force a discontinuous bottom topography - if x[1] >= -1.5 && x[1] <= 0.0 - b = 0.5 + if x[1] >= -1.5f0 && x[1] <= 0 + b = 0.5f0 end return prim2cons(SVector(H, v, b), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index 5851530e230..219f0f9b6dc 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -19,12 +19,12 @@ function initial_condition_discontinuous_well_balancedness(x, t, equations::ShallowWaterEquations1D) # Set the background values H = equations.H0 - v = 0.0 - b = 0.0 + v = 0 + b = 0 # Setup a discontinuous bottom topography - if x[1] >= 0.5 && x[1] <= 0.75 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + if x[1] >= 0.5f0 && x[1] <= 0.75f0 + b = 2 + 0.5f0 * sinpi(2 * x[1]) end return prim2cons(SVector(H, v, b), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl index 9ed02c0e378..a95b7892c4e 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl @@ -11,9 +11,9 @@ equations = ShallowWaterEquations1D(gravity_constant = 1.0, H0 = 3.0) function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations1D) # Set the background values H = equations.H0 - v = 0.0 + v = 0 - b = (1.5 / exp(0.5 * ((x[1] - 1.0)^2)) + 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + b = (1.5f0 / exp(0.5f0 * ((x[1] - 1)^2)) + 0.75f0 / exp(0.5f0 * ((x[1] + 1)^2))) return prim2cons(SVector(H, v, b), equations) end diff --git a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl index d3a17b513fc..4c164a08ec4 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -21,7 +21,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # Discontinuous initial condition (Riemann Problem) leading to a shock that moves to the left. # The shock corresponds to the traffic congestion. function initial_condition_traffic_jam(x, t, equation::TrafficFlowLWREquations1D) - scalar = x[1] < 0.0 ? 0.5 : 1.0 + scalar = x[1] < 0 ? 0.5f0 : 1.0f0 return SVector(scalar) end From a89e7450a5d6e75209ae56efd2910eeb96d26edd Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 5 Nov 2024 17:50:45 -1000 Subject: [PATCH 5/8] Fix --- examples/tree_1d_dgsem/elixir_advection_diffusion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl index b3144d0940d..28e47794c38 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl @@ -22,8 +22,8 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 30_000, # set maximum capacity of tree data structure periodicity = true) -function x_trans_periodic(x, domain_length = SVector(oftype(x, 2 * pi)), - center = SVector(oftype(x, 0))) +function x_trans_periodic(x, domain_length = SVector(oftype(x[1], 2 * pi)), + center = SVector(oftype(x[1], 0))) x_normalized = x .- center x_shifted = x_normalized .% domain_length x_offset = ((x_shifted .< -0.5f0 * domain_length) - From 8c25c476d88e274c053139a09ab75dea73d4eaec Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Thu, 7 Nov 2024 13:55:47 -1000 Subject: [PATCH 6/8] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index 219f0f9b6dc..0bb46129ba4 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -20,7 +20,7 @@ function initial_condition_discontinuous_well_balancedness(x, t, # Set the background values H = equations.H0 v = 0 - b = 0 + b = zero(eltype(x)) # Setup a discontinuous bottom topography if x[1] >= 0.5f0 && x[1] <= 0.75f0 From 1324e982b9977048be5902b5a1f98eb971c9bc02 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Thu, 7 Nov 2024 13:58:23 -1000 Subject: [PATCH 7/8] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/tree_1d_dgsem/elixir_advection_diffusion.jl | 2 -- examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl | 3 ++- examples/tree_1d_dgsem/elixir_burgers_shock.jl | 3 ++- .../elixir_eulermulti_two_interacting_blast_waves.jl | 4 ++-- examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl | 4 ++-- examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl | 4 ++-- examples/tree_1d_dgsem/elixir_shallowwater_ec.jl | 4 ++-- .../elixir_shallowwater_quasi_1d_discontinuous.jl | 4 ++-- examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl | 3 ++- 9 files changed, 16 insertions(+), 15 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl index 28e47794c38..d72e66f8330 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl @@ -41,8 +41,6 @@ function initial_condition_diffusive_convergence_test(x, t, nu = diffusivity() c = 0 A = 1 - L = 2 - f = 1.0f0 / L omega = 1 scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t) return SVector(scalar) diff --git a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl index d73c08d03b1..5e40baf49cc 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl @@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max, # Discontinuous initial condition (Riemann Problem) leading to a rarefaction fan. function initial_condition_rarefaction(x, t, equation::InviscidBurgersEquation1D) - scalar = x[1] < 0.5f0 ? 0.5f0 : 1.5f0 + RealT = eltype(x) + scalar = x[1] < 0.5f0 ? convert(RealT, 0.5f0) : convert(RealT, 1.5f0) return SVector(scalar) end diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl index 89939e4328a..b795e8623fa 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl @@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max, # Discontinuous initial condition (Riemann Problem) leading to a shock to test e.g. correct shock speed. function initial_condition_shock(x, t, equation::InviscidBurgersEquation1D) - scalar = x[1] < 0.5f0 ? 1.5f0 : 0.5f0 + RealT = eltype(x) + scalar = x[1] < 0.5f0 ? convert(RealT, 1.5f0) : convert(RealT, 0.5f0) return SVector(scalar) end diff --git a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl index 908554e0e39..3089e0c6809 100644 --- a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl +++ b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl @@ -28,11 +28,11 @@ function initial_condition_two_interacting_blast_waves(x, t, v1 = 0 if x[1] <= RealT(0.1) - p = 1000 + p = convert(RealT, 1000) elseif x[1] < RealT(0.9) p = convert(RealT, 0.01) else - p = 100 + p = convert(RealT, 100) end prim_other = SVector{2, real(equations)}(v1, p) diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl index b0f68c73689..dda1009cfa8 100644 --- a/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl @@ -21,8 +21,8 @@ function initial_condition_harmonic_nonperiodic(x, t, # elliptic equation: -νΔϕ = f RealT = eltype(x) if t == 0 - phi = 5 - q1 = 0 + phi = convert(RealT, 5) + q1 = zero(RealT) else A = 3 B = exp(one(RealT)) diff --git a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl index 0d2a4418b35..64641679f4f 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl @@ -24,8 +24,8 @@ function initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEqua v3 = 0 p = x[1] <= 0 ? 3 : 1 B1 = 1.5f0 - B2 = x[1] <= 0 ? zero(RealT) : cos(RealT(1.5)) - B3 = x[1] <= 0 ? zero(RealT) : sin(RealT(1.5)) + B2 = x[1] <= 0 ? zero(RealT) : cos(RealT(1.5f0)) + B3 = x[1] <= 0 ? zero(RealT) : sin(RealT(1.5f0)) return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end initial_condition = initial_condition_torrilhon_shock_tube diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index c3a94a27d83..018ce0999eb 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -19,12 +19,12 @@ function initial_condition_ec_discontinuous_bottom(x, t, equations::ShallowWater # Set the background values RealT = eltype(x) H = 4.25f0 - v = 0 + v = zero(RealT) b = sin(x[1]) # arbitrary continuous function # Setup the discontinuous water height and velocity if x[1] >= 0.125f0 && x[1] <= 0.25f0 - H = 5 + H = 5.0f0 v = convert(RealT, 0.1882) end diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl index 15c2f46ae7d..0f9998b60a7 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl @@ -14,9 +14,9 @@ function initial_condition_discontinuity(x, t, equations::ShallowWaterEquationsQ if x[1] > 0 b = convert(RealT, 0.1) - a = 1 + a = one(RealT) else - b = 0 + b = zero(RealT) a = convert(RealT, 1.1) end diff --git a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl index 4c164a08ec4..67d3d56c7af 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -21,7 +21,8 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # Discontinuous initial condition (Riemann Problem) leading to a shock that moves to the left. # The shock corresponds to the traffic congestion. function initial_condition_traffic_jam(x, t, equation::TrafficFlowLWREquations1D) - scalar = x[1] < 0 ? 0.5f0 : 1.0f0 + RealT = eltype(x) + scalar = x[1] < 0 ? RealT(0.5f0) : one(RealT) return SVector(scalar) end From e82cf3e8e9ca12d1b08ca186f78c9496c2162ee3 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 7 Nov 2024 18:55:10 -1000 Subject: [PATCH 8/8] Fix --- examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl index 64641679f4f..c35212b4074 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl @@ -24,7 +24,7 @@ function initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEqua v3 = 0 p = x[1] <= 0 ? 3 : 1 B1 = 1.5f0 - B2 = x[1] <= 0 ? zero(RealT) : cos(RealT(1.5f0)) + B2 = x[1] <= 0 ? one(RealT) : cos(RealT(1.5f0)) B3 = x[1] <= 0 ? zero(RealT) : sin(RealT(1.5f0)) return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end