From 82f229437fdfae0684558847f13f142567db3e98 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 18 Jul 2024 21:08:35 -1000 Subject: [PATCH 1/5] start --- src/equations/polytropic_euler_2d.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index e900fd64073..73eb2e48fe0 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -113,9 +113,9 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat phi = atan(y_norm, x_norm) # 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) + rho = r > 0.5f0 ? 1.0 : 1.1691 + v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos(phi) + v2 = r > 0.5f0 ? 0.0 : 0.1882 * sin(phi) return prim2cons(SVector(rho, v1, v2), equations) end @@ -186,12 +186,12 @@ For details see Section 3.2 of the following reference else # equations.gamma > 1 # polytropic gas rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) end - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] @@ -212,16 +212,16 @@ end else # equations.gamma > 1 # polytropic gas rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) end - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) if orientation == 1 # x-direction - f1 = rho_mean * 0.5 * (v1_ll + v1_rr) + f1 = rho_mean * 0.5f0 * (v1_ll + v1_rr) f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg else # y-direction - f1 = rho_mean * 0.5 * (v2_ll + v2_rr) + f1 = rho_mean * 0.5f0 * (v2_ll + v2_rr) f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg end @@ -367,7 +367,7 @@ end (equations.gamma - 1.0) end - w1 = internal_energy + p / rho - 0.5 * v_square + w1 = internal_energy + p / rho - 0.5f0 * v_square w2 = v1 w3 = v2 From d1511220b3d3b630231441ed78df940bcaf0e2b7 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sun, 28 Jul 2024 14:45:16 -1000 Subject: [PATCH 2/5] fix flux tests --- test/test_type.jl | 314 +++++++++++++--------------------------------- 1 file changed, 86 insertions(+), 228 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index 6c51460e6d9..39f8895fdd3 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -45,39 +45,17 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_wall(u_inner, - orientation, - direction, x, - t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_wall(u_inner, orientation, - direction, x, t, - surface_flux_function, - equations)) == RealT - end + @test eltype(@inferred boundary_condition_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations)) == RealT end end - - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, - normal_direction, - x, t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - normal_direction, x, t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, x, t, + surface_flux_function, + equations)) == + RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, @@ -850,25 +828,13 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT end @test eltype(@inferred flux(u, orientation, equations)) == RealT @@ -911,25 +877,13 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT end end @@ -985,25 +939,13 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT end end @@ -1631,24 +1573,13 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, - orientation, - direction, x, - t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, x, + t, + surface_flux_function, + equations)) == + RealT end @test eltype(@inferred flux(u, orientation, equations)) == RealT @@ -1701,58 +1632,30 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - @test_broken eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - @test eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT end end @@ -1806,26 +1709,14 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT end end @@ -1904,24 +1795,12 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == RealT - end - + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == RealT @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, direction, equations)) == RealT @@ -2006,22 +1885,12 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == RealT + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == RealT - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, - normal_direction, - x, t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - normal_direction, - x, t, - surface_flux_function, - equations)) == RealT - end @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, normal_direction_ll, @@ -2063,23 +1932,12 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - if RealT == Float32 - # check `surface_flux_function` (test broken) - @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT - else - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - orientation, - direction, x, t, - surface_flux_function, - equations)) == - RealT - end + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, x, t, + surface_flux_function, + equations)) == + RealT end @test eltype(@inferred flux(u, orientation, equations)) == RealT From c6fecc4e6c221e7cf40a6eff75cbbb75f9d150d0 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sun, 28 Jul 2024 15:08:46 -1000 Subject: [PATCH 3/5] keep conflicts --- test/test_type.jl | 314 +++++++++++++++++++++++++++++++++------------- 1 file changed, 228 insertions(+), 86 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index 39f8895fdd3..6c51460e6d9 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -45,17 +45,39 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - @test eltype(@inferred boundary_condition_wall(u_inner, orientation, - direction, x, t, - surface_flux_function, - equations)) == RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_wall(u_inner, + orientation, + direction, x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations)) == RealT + end end end - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - normal_direction, x, t, - surface_flux_function, - equations)) == - RealT + + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, x, t, + surface_flux_function, + equations)) == + RealT + end @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, @@ -828,13 +850,25 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for direction in directions - @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end end @test eltype(@inferred flux(u, orientation, equations)) == RealT @@ -877,13 +911,25 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end end end @@ -939,13 +985,25 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == - RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end end end @@ -1573,13 +1631,24 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for direction in directions - @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, - orientation, - direction, x, - t, - surface_flux_function, - equations)) == - RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, x, + t, + surface_flux_function, + equations)) == + RealT + end end @test eltype(@inferred flux(u, orientation, equations)) == RealT @@ -1632,30 +1701,58 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - @test eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT - @test eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test_broken eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + end end end @@ -1709,14 +1806,26 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - @test eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, - orientation, - direction, - x, - t, - surface_flux_function, - equations)) == - RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + end end end @@ -1795,12 +1904,24 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for direction in directions - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - orientation, - direction, - x, t, - surface_flux_function, - equations)) == RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == RealT + end + @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, direction, equations)) == RealT @@ -1885,12 +2006,22 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == RealT - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - normal_direction, - x, t, - surface_flux_function, - equations)) == RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == RealT + end @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, normal_direction_ll, @@ -1932,12 +2063,23 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations for direction in directions - @test eltype(@inferred boundary_condition_slip_wall(u_inner, - orientation, - direction, x, t, - surface_flux_function, - equations)) == - RealT + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, x, t, + surface_flux_function, + equations)) == + RealT + end end @test eltype(@inferred flux(u, orientation, equations)) == RealT From 21b20124b1fa1ed731c759b0fb1778d8d8fd2fe0 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Sun, 28 Jul 2024 15:38:18 -1000 Subject: [PATCH 4/5] complete equations --- src/equations/polytropic_euler_2d.jl | 32 +++++++++++++++------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 73eb2e48fe0..571d4723f6f 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -62,7 +62,7 @@ in combination with [`source_terms_convergence_test`](@ref). function initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D) # manufactured initial condition from Winters (2019) [0.1007/s10543-019-00789-w] # domain must be set to [0, 1] x [0, 1] - h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) + h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t) return SVector(h, h / 2, 3 * h / 2) end @@ -78,19 +78,20 @@ Source terms used for convergence tests in combination with rho, v1, v2 = cons2prim(u, equations) # Residual from Winters (2019) [0.1007/s10543-019-00789-w] eq. (5.2). - h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) - h_t = -2 * pi * cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * sin(2 * pi * t) - h_x = -2 * pi * sin(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) - h_y = 2 * pi * cos(2 * pi * x[1]) * cos(2 * pi * x[2]) * cos(2 * pi * t) + RealT = eltype(u) + h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t) + h_t = -2 * convert(RealT, pi) * cospi(2 * x[1]) * sinpi(2 * x[2]) * sinpi(2 * t) + h_x = -2 * convert(RealT, pi) * sinpi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t) + h_y = 2 * convert(RealT, pi) * cospi(2 * x[1]) * cospi(2 * x[2]) * cospi(2 * t) rho_x = h_x rho_y = h_y b = equations.kappa * equations.gamma * h^(equations.gamma - 1) - r_1 = h_t + h_x / 2 + 3 / 2 * h_y - r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 / 4 * h_y - r_3 = 3 / 2 * h_t + 3 / 4 * h_x + 9 / 4 * h_y + b * rho_y + r_1 = h_t + h_x / 2 + 3 * h_y / 2 + r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 * h_y / 4 + r_3 = 3 * h_t / 2 + 3 * h_x / 4 + 9 * h_y / 4 + b * rho_y return SVector(r_1, r_2, r_3) end @@ -113,9 +114,10 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat phi = atan(y_norm, x_norm) # Calculate primitive variables - rho = r > 0.5f0 ? 1.0 : 1.1691 - v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos(phi) - v2 = r > 0.5f0 ? 0.0 : 0.1882 * sin(phi) + RealT = eltype(x) + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) return prim2cons(SVector(rho, v1, v2), equations) end @@ -181,7 +183,7 @@ For details see Section 3.2 of the following reference v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # Compute the necessary mean values - if equations.gamma == 1.0 # isothermal gas + if equations.gamma == 1 # isothermal gas rho_mean = ln_mean(rho_ll, rho_rr) else # equations.gamma > 1 # polytropic gas rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) @@ -207,7 +209,7 @@ end p_rr = equations.kappa * rho_rr^equations.gamma # Compute the necessary mean values - if equations.gamma == 1.0 # isothermal gas + if equations.gamma == 1 # isothermal gas rho_mean = ln_mean(rho_ll, rho_rr) else # equations.gamma > 1 # polytropic gas rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) @@ -360,11 +362,11 @@ end v_square = v1^2 + v2^2 p = pressure(u, equations) # Form of the internal energy depends on gas type - if equations.gamma == 1.0 # isothermal gas + if equations.gamma == 1 # isothermal gas internal_energy = equations.kappa * log(rho) else # equations.gamma > 1 # polytropic gas internal_energy = equations.kappa * rho^(equations.gamma - 1) / - (equations.gamma - 1.0) + (equations.gamma - 1) end w1 = internal_energy + p / rho - 0.5f0 * v_square From eff115adb8cd6253829950b39a62b82405d3c548 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Thu, 15 Aug 2024 15:49:41 -1000 Subject: [PATCH 5/5] complete unit test --- test/test_type.jl | 74 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index 9b8c3366cfb..d22afa65c0a 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1926,6 +1926,80 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Polytropic Euler 2D" begin + for RealT in (Float32, Float64) + equations1 = @inferred PolytropicEulerEquations2D(RealT(1), + RealT(1)) # equations.gamma == 1 + equations2 = @inferred PolytropicEulerEquations2D(RealT(1.4), RealT(0.5)) # equations.gamma > 1 + + for equations in (equations1, equations2) + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = SVector(one(RealT), zero(RealT)) + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `stolarsky_mean` (test broken) + @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + else + @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction, + equations)) == + RealT + end + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `stolarsky_mean` (test broken) + @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation, + equations)) == + RealT + end + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == + RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + end + end + end + @timed_testset "Shallow Water 1D" begin for RealT in (Float32, Float64) equations = @inferred ShallowWaterEquations1D(gravity_constant = RealT(9.81))