From 3c1a1a99be4dee4cdcaba5ec14c922cd15dbaf74 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 16 Jul 2024 19:24:25 -1000 Subject: [PATCH 1/2] fix equations --- test/test_unit.jl | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 09cf506576c..aa66ce2556f 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -439,14 +439,17 @@ end # Neither Mach number nor velocity set @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = nothing, Re = 1000) # Both Mach number and velocity set - @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, u0 = 1) + @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, + u0 = 1.0) # Neither Reynolds number nor viscosity set @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = nothing) # Both Reynolds number and viscosity set - @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, nu = 1) + @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, + nu = 1.0) # No non-dimensional values set - @test LatticeBoltzmannEquations2D(Ma = nothing, Re = nothing, u0 = 1, nu = 1) isa + @test LatticeBoltzmannEquations2D(Ma = nothing, Re = nothing, u0 = 1.0, + nu = 1.0) isa LatticeBoltzmannEquations2D end @@ -454,14 +457,17 @@ end # Neither Mach number nor velocity set @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = nothing, Re = 1000) # Both Mach number and velocity set - @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, u0 = 1) + @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, + u0 = 1.0) # Neither Reynolds number nor viscosity set @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = nothing) # Both Reynolds number and viscosity set - @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, nu = 1) + @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, + nu = 1.0) # No non-dimensional values set - @test LatticeBoltzmannEquations3D(Ma = nothing, Re = nothing, u0 = 1, nu = 1) isa + @test LatticeBoltzmannEquations3D(Ma = nothing, Re = nothing, u0 = 1.0, + nu = 1.0) isa LatticeBoltzmannEquations3D end From 4d8be9031d32be44c11575bda94e35699221f552 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 16 Jul 2024 22:17:10 -1000 Subject: [PATCH 2/2] partial complete --- src/equations/shallow_water_1d.jl | 56 +++++++-------- src/equations/shallow_water_2d.jl | 90 +++++++++++++------------ src/equations/shallow_water_quasi_1d.jl | 26 +++---- 3 files changed, 89 insertions(+), 83 deletions(-) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 7007bea887a..a48d14d181f 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -74,13 +74,14 @@ A smooth initial condition used for convergence tests in combination with function initial_condition_convergence_test(x, t, equations::ShallowWaterEquations1D) # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] - c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi + RealT = eltype(x) + c = 7 + omega_x = 2 * convert(RealT, pi) * sqrt(2.0) + omega_t = 2 * convert(RealT, pi) H = c + cos(omega_x * x[1]) * cos(omega_t * t) - v = 0.5 - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x[1]) + v = 0.5f0 + b = 2.0 + 0.5f0 * sinpi(sqrt(2.0) * x[1]) return prim2cons(SVector(H, v, b), equations) end @@ -92,7 +93,7 @@ Source terms used for convergence tests in combination with (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). This manufactured solution source term is specifically designed for the bottom topography function -`b(x) = 2.0 + 0.5 * sin(sqrt(2.0)*pi*x[1])` +`b(x) = 2.0 + 0.5 * sinpi(sqrt(2.0)*x[1])` as defined in [`initial_condition_convergence_test`](@ref). """ @@ -100,11 +101,12 @@ as defined in [`initial_condition_convergence_test`](@ref). equations::ShallowWaterEquations1D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocity is taken to be constant + RealT = eltype(u) c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi - omega_b = sqrt(2.0) * pi - v = 0.5 + omega_x = 2.0 * convert(RealT, pi) * sqrt(2.0) + omega_t = 2.0 * convert(RealT, pi) + omega_b = sqrt(2.0) * convert(RealT, pi) + v = 0.5f0 sinX, cosX = sincos(omega_x * x[1]) sinT, cosT = sincos(omega_t * t) @@ -116,8 +118,8 @@ as defined in [`initial_condition_convergence_test`](@ref). H_t = -omega_t * cosX * sinT # bottom topography and its spatial derivative - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x[1]) - b_x = 0.5 * omega_b * cos(omega_b * x[1]) + b = 2.0 + 0.5f0 * sinpi(sqrt(2.0) * x[1]) + b_x = 0.5f0 * omega_b * cos(omega_b * x[1]) du1 = H_t + v * (H_x - b_x) du2 = v * du1 + equations.gravity * (H - b) * H_x @@ -136,8 +138,8 @@ function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquation r = abs(x_norm) # Calculate primitive variables - H = r > 0.5 ? 3.25 : 4.0 - v = r > 0.5 ? 0.0 : 0.1882 + H = r > 0.5f0 ? 3.25 : 4.0 + v = r > 0.5f0 ? 0.0 : 0.1882 b = sin(x[1]) # arbitrary continuous function return prim2cons(SVector(H, v, b), equations) @@ -185,7 +187,7 @@ end h, h_v, _ = u v = velocity(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 f1 = h_v f2 = h_v * v + p @@ -249,7 +251,7 @@ and for curvilinear 2D case in the paper: h_ll, _, b_ll = u_ll h_rr, _, b_rr = u_rr - h_average = 0.5 * (h_ll + h_rr) + h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll # Includes two parts: @@ -367,8 +369,8 @@ Details are available in Eq. (4.1) in the paper: v_rr = velocity(u_rr, equations) # Average each factor of products in flux - h_avg = 0.5 * (h_ll + h_rr) - v_avg = 0.5 * (v_ll + v_rr) + h_avg = 0.5f0 * (h_ll + h_rr) + v_avg = 0.5f0 * (v_ll + v_rr) p_avg = 0.25 * equations.gravity * (h_ll^2 + h_rr^2) # Calculate fluxes depending on orientation @@ -403,11 +405,11 @@ Further details are available in Theorem 1 of the paper: v_rr = velocity(u_rr, equations) # Average each factor of products in flux - v_avg = 0.5 * (v_ll + v_rr) - p_avg = 0.5 * equations.gravity * h_ll * h_rr + v_avg = 0.5f0 * (v_ll + v_rr) + p_avg = 0.5f0 * equations.gravity * h_ll * h_rr # Calculate fluxes depending on orientation - f1 = 0.5 * (h_v_ll + h_v_rr) + f1 = 0.5f0 * (h_v_ll + h_v_rr) f2 = f1 * v_avg + p_avg return SVector(f1, f2, zero(eltype(u_ll))) @@ -471,7 +473,7 @@ end equations::ShallowWaterEquations1D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) + diss = -0.5f0 * λ * (u_rr - u_ll) return SVector(diss[1], diss[2], zero(eltype(u_ll))) end @@ -581,7 +583,7 @@ end v = velocity(u, equations) - w1 = equations.gravity * (h + b) - 0.5 * v^2 + w1 = equations.gravity * (h + b) - 0.5f0 * v^2 w2 = v return SVector(w1, w2, b) @@ -591,7 +593,7 @@ end @inline function entropy2cons(w, equations::ShallowWaterEquations1D) w1, w2, b = w - h = (w1 + 0.5 * w2^2) / equations.gravity - b + h = (w1 + 0.5f0 * w2^2) / equations.gravity - b h_v = h * w2 return SVector(h, h_v, b) end @@ -612,7 +614,7 @@ end @inline function pressure(u, equations::ShallowWaterEquations1D) h = waterheight(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 return p end @@ -638,7 +640,7 @@ Or equation (9.17) in [this lecture notes](https://metaphor.ethz.ch/x/2019/hs/40 h_rr = waterheight(u_rr, equations) v_rr = velocity(u_rr, equations) - h_roe = 0.5 * (h_ll + h_rr) + h_roe = 0.5f0 * (h_ll + h_rr) c_roe = sqrt(equations.gravity * h_roe) h_ll_sqrt = sqrt(h_ll) @@ -658,7 +660,7 @@ end @inline function energy_total(cons, equations::ShallowWaterEquations1D) h, h_v, b = cons - e = (h_v^2) / (2 * h) + 0.5 * equations.gravity * h^2 + equations.gravity * h * b + e = (h_v^2) / (2 * h) + 0.5f0 * equations.gravity * h^2 + equations.gravity * h * b return e end diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 73fae89a0fa..327c5ebd2bd 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -77,16 +77,17 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::ShallowWaterEquations2D) # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2 + RealT = eltype(x) c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi + omega_x = 2.0 * convert(RealT, pi) * sqrt(2.0) + omega_t = 2.0 * convert(RealT, pi) x1, x2 = x H = c + cos(omega_x * x1) * sin(omega_x * x2) * cos(omega_t * t) - v1 = 0.5 + v1 = 0.5f0 v2 = 1.5 - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2) + b = 2.0 + 0.5f0 * sinpi(sqrt(2.0) * x1) + 0.5f0 * sinpi(sqrt(2.0) * x2) return prim2cons(SVector(H, v1, v2, b), equations) end @@ -98,18 +99,19 @@ Source terms used for convergence tests in combination with (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). This manufactured solution source term is specifically designed for the bottom topography function -`b(x,y) = 2 + 0.5 * sin(sqrt(2)*pi*x) + 0.5 * sin(sqrt(2)*pi*y)` +`b(x,y) = 2 + 0.5 * sinpi(sqrt(2)*x) + 0.5 * sinpi(sqrt(2)*y)` as defined in [`initial_condition_convergence_test`](@ref). """ @inline function source_terms_convergence_test(u, x, t, equations::ShallowWaterEquations2D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocities are taken to be constants + RealT = eltype(u) c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi - omega_b = sqrt(2.0) * pi - v1 = 0.5 + omega_x = 2.0 * convert(RealT, pi) * sqrt(2.0) + omega_t = 2.0 * convert(RealT, pi) + omega_b = sqrt(2.0) * convert(RealT, pi) + v1 = 0.5f0 v2 = 1.5 x1, x2 = x @@ -126,8 +128,8 @@ as defined in [`initial_condition_convergence_test`](@ref). H_t = -omega_t * cosX * sinY * sinT # bottom topography and its gradient - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2) - tmp1 = 0.5 * omega_b + b = 2.0 + 0.5f0 * sinpi(sqrt(2.0) * x1) + 0.5f0 * sinpi(sqrt(2.0) * x2) + tmp1 = 0.5f0 * omega_b b_x = tmp1 * cos(omega_b * x1) b_y = tmp1 * cos(omega_b * x2) @@ -153,9 +155,9 @@ function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquation sin_phi, cos_phi = sincos(phi) # Calculate primitive variables - H = r > 0.5 ? 3.25 : 4.0 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi + H = r > 0.5f0 ? 3.25 : 4.0 + v1 = r > 0.5f0 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5f0 ? 0.0 : 0.1882 * sin_phi b = 0.0 # by default assume there is no bottom topography return prim2cons(SVector(H, v1, v2, b), equations) @@ -228,7 +230,7 @@ end h, h_v1, h_v2, _ = u v1, v2 = velocity(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 if orientation == 1 f1 = h_v1 f2 = h_v1 * v1 + p @@ -250,7 +252,7 @@ end v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] h_v_normal = h * v_normal - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 f1 = h_v_normal f2 = h_v_normal * v1 + p * normal_direction[1] @@ -349,7 +351,7 @@ and for curvilinear 2D case in the paper: h_ll, _, _, b_ll = u_ll h_rr, _, _, b_rr = u_rr - h_average = 0.5 * (h_ll + h_rr) + h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll # Includes two parts: @@ -389,7 +391,7 @@ end # (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump` # to handle discontinuous bathymetry - h_average = 0.5 * (h_ll + h_rr) + h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll f2 += normal_direction_ll[1] * equations.gravity * h_average * b_jump @@ -611,9 +613,9 @@ Details are available in Eq. (4.1) in the paper: v1_rr, v2_rr = velocity(u_rr, equations) # Average each factor of products in flux - h_avg = 0.5 * (h_ll + h_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) + h_avg = 0.5f0 * (h_ll + h_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) p_avg = 0.25 * equations.gravity * (h_ll^2 + h_rr^2) # Calculate fluxes depending on orientation @@ -642,12 +644,12 @@ end v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # Average each factor of products in flux - h_avg = 0.5 * (h_ll + h_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - h2_avg = 0.5 * (h_ll^2 + h_rr^2) - p_avg = 0.5 * equations.gravity * h2_avg - v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr) + h_avg = 0.5f0 * (h_ll + h_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + h2_avg = 0.5f0 * (h_ll^2 + h_rr^2) + p_avg = 0.5f0 * equations.gravity * h2_avg + v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) # Calculate fluxes depending on normal_direction f1 = h_avg * v_dot_n_avg @@ -682,17 +684,17 @@ Further details are available in Theorem 1 of the paper: v1_rr, v2_rr = velocity(u_rr, equations) # Average each factor of products in flux - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * equations.gravity * h_ll * h_rr + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * equations.gravity * h_ll * h_rr # Calculate fluxes depending on orientation if orientation == 1 - f1 = 0.5 * (h_v1_ll + h_v1_rr) + f1 = 0.5f0 * (h_v1_ll + h_v1_rr) f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg else - f1 = 0.5 * (h_v2_ll + h_v2_rr) + f1 = 0.5f0 * (h_v2_ll + h_v2_rr) f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg end @@ -711,11 +713,11 @@ end v1_rr, v2_rr = velocity(u_rr, equations) # Average each factor of products in flux - h_v1_avg = 0.5 * (h_v1_ll + h_v1_rr) - h_v2_avg = 0.5 * (h_v2_ll + h_v2_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * equations.gravity * h_ll * h_rr + h_v1_avg = 0.5f0 * (h_v1_ll + h_v1_rr) + h_v2_avg = 0.5f0 * (h_v2_ll + h_v2_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * equations.gravity * h_ll * h_rr # Calculate fluxes depending on normal_direction f1 = h_v1_avg * normal_direction[1] + h_v2_avg * normal_direction[2] @@ -771,7 +773,7 @@ end equations::ShallowWaterEquations2D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) + diss = -0.5f0 * λ * (u_rr - u_ll) return SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll))) end @@ -958,7 +960,7 @@ end v1, v2 = velocity(u, equations) v_square = v1^2 + v2^2 - w1 = equations.gravity * (h + b) - 0.5 * v_square + w1 = equations.gravity * (h + b) - 0.5f0 * v_square w2 = v1 w3 = v2 return SVector(w1, w2, w3, b) @@ -968,7 +970,7 @@ end @inline function entropy2cons(w, equations::ShallowWaterEquations2D) w1, w2, w3, b = w - h = (w1 + 0.5 * (w2^2 + w3^2)) / equations.gravity - b + h = (w1 + 0.5f0 * (w2^2 + w3^2)) / equations.gravity - b h_v1 = h * w2 h_v2 = h * w3 return SVector(h, h_v1, h_v2, b) @@ -990,7 +992,7 @@ end @inline function pressure(u, equations::ShallowWaterEquations2D) h = waterheight(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 return p end @@ -1017,7 +1019,7 @@ slides 8 and 9. h_rr = waterheight(u_rr, equations) v1_rr, v2_rr = velocity(u_rr, equations) - h_roe = 0.5 * (h_ll + h_rr) + h_roe = 0.5f0 * (h_ll + h_rr) c_roe = sqrt(equations.gravity * h_roe) h_ll_sqrt = sqrt(h_ll) @@ -1041,7 +1043,7 @@ end norm_ = norm(normal_direction) - h_roe = 0.5 * (h_ll + h_rr) + h_roe = 0.5f0 * (h_ll + h_rr) c_roe = sqrt(equations.gravity * h_roe) * norm_ h_ll_sqrt = sqrt(h_ll) @@ -1064,7 +1066,7 @@ end @inline function energy_total(cons, equations::ShallowWaterEquations2D) h, h_v1, h_v2, b = cons - e = (h_v1^2 + h_v2^2) / (2 * h) + 0.5 * equations.gravity * h^2 + + e = (h_v1^2 + h_v2^2) / (2 * h) + 0.5f0 * equations.gravity * h^2 + equations.gravity * h * b return e end diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index 08620021254..9e2069c4833 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/src/equations/shallow_water_quasi_1d.jl @@ -72,8 +72,9 @@ function initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D) # generates a manufactured solution. # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] - Omega = sqrt(2) * pi - H = 2.0 + 0.5 * sin(Omega * x[1] - t) + RealT = eltype(x) + Omega = sqrt(2) * convert(RealT, pi) + H = 2.0 + 0.5f0 * sin(Omega * x[1] - t) v = 0.25 b = 0.2 - 0.05 * sin(Omega * x[1]) a = 1 + 0.1 * cos(Omega * x[1]) @@ -88,17 +89,18 @@ Source terms used for convergence tests in combination with (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). This manufactured solution source term is specifically designed for the bottom topography function -`b(x) = 0.2 - 0.05 * sin(sqrt(2) * pi *x[1])` and channel width 'a(x)= 1 + 0.1 * cos(sqrt(2) * pi * x[1])' +`b(x) = 0.2 - 0.05 * sinpi(sqrt(2) *x[1])` and channel width 'a(x)= 1 + 0.1 * cospi(sqrt(2) * x[1])' as defined in [`initial_condition_convergence_test`](@ref). """ @inline function source_terms_convergence_test(u, x, t, equations::ShallowWaterEquationsQuasi1D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocity is taken to be constant - Omega = sqrt(2) * pi - H = 2.0 + 0.5 * sin(Omega * x[1] - t) - H_x = 0.5 * cos(Omega * x[1] - t) * Omega - H_t = -0.5 * cos(Omega * x[1] - t) + RealT = eltype(u) + Omega = sqrt(2) * convert(RealT, pi) + H = 2.0 + 0.5f0 * sin(Omega * x[1] - t) + H_x = 0.5f0 * cos(Omega * x[1] - t) * Omega + H_t = -0.5f0 * cos(Omega * x[1] - t) v = 0.25 @@ -200,8 +202,8 @@ Further details are available in the paper: v_ll = velocity(u_ll, equations) v_rr = velocity(u_rr, equations) - f1 = 0.5 * (a_h_v_ll + a_h_v_rr) - f2 = f1 * 0.5 * (v_ll + v_rr) + f1 = 0.5f0 * (a_h_v_ll + a_h_v_rr) + f2 = f1 * 0.5f0 * (v_ll + v_rr) return SVector(f1, f2, zero(eltype(u_ll)), zero(eltype(u_ll))) end @@ -240,7 +242,7 @@ end equations::ShallowWaterEquationsQuasi1D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) + diss = -0.5f0 * λ * (u_rr - u_ll) return SVector(diss[1], diss[2], zero(eltype(u_ll)), zero(eltype(u_ll))) end @@ -278,7 +280,7 @@ end h = waterheight(u, equations) v = velocity(u, equations) #entropy variables are the same as ones in standard shallow water equations - w1 = equations.gravity * (h + b) - 0.5 * v^2 + w1 = equations.gravity * (h + b) - 0.5f0 * v^2 w2 = v return SVector(w1, w2, b, a) @@ -306,7 +308,7 @@ end # Calculate total energy for a conservative state `cons` @inline function energy_total(cons, equations::ShallowWaterEquationsQuasi1D) a_h, a_h_v, b, a = cons - e = (a_h_v^2) / (2 * a * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + + e = (a_h_v^2) / (2 * a * a_h) + 0.5f0 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b return e end