From ed0e8dcdd2aefc2b9715ac936f843fc47ca9009b Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 18 Dec 2023 13:37:13 +0100 Subject: [PATCH 01/41] add LMARS flux in 2D --- src/equations/compressible_euler_2d.jl | 87 ++++++++++++++++++++++++++ src/equations/compressible_euler_3d.jl | 5 -- 2 files changed, 87 insertions(+), 5 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index a992f99eaf4..0e56f934bbf 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -809,6 +809,93 @@ end return SVector(f1m, f2m, f3m, f4m) end +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquations2D) + +Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using +an estimate `c` of the speed of sound. + +References: +- Xi Chen et al. (2013) + A Control-Volume Model of the Compressible Euler Equations with a Vertical + Lagrangian Coordinate + [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) +""" +struct FluxLMARS{SpeedOfSound} + # Estimate for the speed of sound + speed_of_sound::SpeedOfSound +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquations2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + + if v >= 0 + f1, f2, f3, f4 = v * u_ll + else + f1, f2, f3, f4 = v * u_rr + end + + if orientation == 1 + f2 += p + else # orientation == 2 + f3 += p + end + f4 += p * v + + return SVector(f1, f2, f3, f4) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Note that this is the same as computing v_ll and v_rr with a normalized normal vector + # and then multiplying v by `norm_` again, but this version is slightly faster. + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + + if v >= 0 + f1, f2, f3, f4 = u_ll * v + # f4 += p_ll * v ??? + else + f1, f2, f3, f4 = u_rr * v + # f4 += p_rr * v ??? + end + + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4 + p * v) +end + """ splitting_vanleer_haenel(u, orientation::Integer, equations::CompressibleEulerEquations2D) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index fc56f58025b..66d35edef43 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -944,11 +944,6 @@ References: Lagrangian Coordinate [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) """ -struct FluxLMARS{SpeedOfSound} - # Estimate for the speed of sound - speed_of_sound::SpeedOfSound -end - @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations3D) c = flux_lmars.speed_of_sound From 77f0d461bbac221b864f0bffd913ccce37bdcfb2 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 18 Dec 2023 14:00:47 +0100 Subject: [PATCH 02/41] add dry air warm bubble test case --- .../elixir_euler_warm_bubble.jl | 124 ++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl new file mode 100644 index 00000000000..abe96ecd70a --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -0,0 +1,124 @@ +using OrdinaryDiffEq +using Trixi + +# Physical constants +g::Float64 = 9.81 # gravity of earth +p_0::Float64 = 100_000.0 # reference pressure +c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) +c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) +R = c_p - c_v # gas constant (dry air) +gamma = c_p / c_v # heat capacity ratio (dry air) + +# Warm bubble test from +# Wicker, L. J., and Skamarock, W. C., 1998: A time-splitting scheme +# for the elastic equations incorporating second-order Runge–Kutta +# time differencing. Mon. Wea. Rev., 126, 1992–1999. +function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D) + # center of perturbation + xc = 10000.0 + zc = 2000.0 + # radius of perturbation + rc = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) + + # perturbation in potential temperature + θ_ref = 300.0 + Δθ = 0.0 + if r <= rc + Δθ = 2 * cospi(0.5*r/rc)^2 + end + θ = θ_ref + Δθ # potential temperature + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * θ) * x[2] + + # pressure + p = p_0 * exner^(c_p/R) + + # temperature + T = θ * exner + + # density + rho = p / (R*T) + + v1 = 20.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + +@inline function source_terms_gravity(u, x, t, + equations::CompressibleEulerEquations2D) + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(gamma) + +boundary_conditions = (x_neg=boundary_condition_periodic, + x_pos=boundary_condition_periodic, + y_neg=boundary_condition_slip_wall, + y_pos=boundary_condition_slip_wall) + +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = FluxLMARS(340.0) +volume_flux = flux_ranocha +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = ( 0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +cells_per_dimension = (64, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver, + source_terms = source_terms_gravity, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + extra_analysis_errors=(:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=analysis_interval, + save_initial_solution=true, + save_final_solution=true, + output_directory="out.struct_lmars_ra", + solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl=0.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + maxiters=1.0e7, + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +summary_callback() From c8b53606e7c94acb420eba2bbedd697545033ae8 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 18 Dec 2023 15:43:24 +0100 Subject: [PATCH 03/41] get formatting right --- .../elixir_euler_warm_bubble.jl | 61 +++++++++---------- 1 file changed, 30 insertions(+), 31 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index abe96ecd70a..63991d45b88 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -2,12 +2,12 @@ using OrdinaryDiffEq using Trixi # Physical constants -g::Float64 = 9.81 # gravity of earth -p_0::Float64 = 100_000.0 # reference pressure -c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) -c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) -R = c_p - c_v # gas constant (dry air) -gamma = c_p / c_v # heat capacity ratio (dry air) +g::Float64 = 9.81 # gravity of earth +p_0::Float64 = 100_000.0 # reference pressure +c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) +c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) +R = c_p - c_v # gas constant (dry air) +gamma = c_p / c_v # heat capacity ratio (dry air) # Warm bubble test from # Wicker, L. J., and Skamarock, W. C., 1998: A time-splitting scheme @@ -26,7 +26,7 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio θ_ref = 300.0 Δθ = 0.0 if r <= rc - Δθ = 2 * cospi(0.5*r/rc)^2 + Δθ = 2 * cospi(0.5 * r / rc)^2 end θ = θ_ref + Δθ # potential temperature @@ -34,17 +34,17 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio exner = 1 - g / (c_p * θ) * x[2] # pressure - p = p_0 * exner^(c_p/R) + p = p_0 * exner^(c_p / R) # temperature T = θ * exner - + # density - rho = p / (R*T) - + rho = p / (R * T) + v1 = 20.0 v2 = 0.0 - E = c_v * T + 0.5 * (v1^2 + v2^2) + E = c_v * T + 0.5 * (v1^2 + v2^2) return SVector(rho, rho * v1, rho * v2, rho * E) end @@ -54,16 +54,15 @@ end return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) end - ############################################################################### # semidiscretization of the compressible Euler equations equations = CompressibleEulerEquations2D(gamma) -boundary_conditions = (x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, - y_neg=boundary_condition_slip_wall, - y_pos=boundary_condition_slip_wall) +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) polydeg = 4 basis = LobattoLegendreBasis(polydeg) @@ -74,7 +73,7 @@ volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) -coordinates_min = ( 0.0, 0.0) +coordinates_min = (0.0, 0.0) coordinates_max = (20_000.0, 10_000.0) cells_per_dimension = (64, 32) @@ -95,18 +94,18 @@ summary_callback = SummaryCallback() analysis_interval = 1000 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, - extra_analysis_errors=(:entropy_conservation_error,)) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) -alive_callback = AliveCallback(analysis_interval=analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval=analysis_interval, - save_initial_solution=true, - save_final_solution=true, - output_directory="out.struct_lmars_ra", - solution_variables=cons2prim) +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out.struct_lmars_ra", + solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl=0.2) +stepsize_callback = StepsizeCallback(cfl = 0.2) callbacks = CallbackSet(summary_callback, analysis_callback, @@ -116,9 +115,9 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - maxiters=1.0e7, - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() From 786ccd990780d28a7171775a89864029b0f2410b Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 19 Dec 2023 13:02:11 +0100 Subject: [PATCH 04/41] remove += --- src/equations/compressible_euler_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 0e56f934bbf..8abf365840b 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -854,11 +854,11 @@ end end if orientation == 1 - f2 += p + f2 = f2 + p else # orientation == 2 - f3 += p + f3 = f3 + p end - f4 += p * v + f4 = f4 + p * v return SVector(f1, f2, f3, f4) end From 62512dc395171ceca6cce9dbdc3db3b99a4ba141 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 19 Dec 2023 13:07:08 +0100 Subject: [PATCH 05/41] add DOI --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index 63991d45b88..5e3bba2e77b 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -9,10 +9,11 @@ c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) R = c_p - c_v # gas constant (dry air) gamma = c_p / c_v # heat capacity ratio (dry air) -# Warm bubble test from -# Wicker, L. J., and Skamarock, W. C., 1998: A time-splitting scheme -# for the elastic equations incorporating second-order Runge–Kutta -# time differencing. Mon. Wea. Rev., 126, 1992–1999. +# Warm bubble test case from +# Wicker, L. J., and Skamarock, W. C. +# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta +# time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D) # center of perturbation xc = 10000.0 From 8043aafcae960eed263ebfb95c21fa80a52d58fa Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 19 Dec 2023 13:18:07 +0100 Subject: [PATCH 06/41] cleaned up variables (naming, scope) --- .../elixir_euler_warm_bubble.jl | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index 5e3bba2e77b..a962aade6e1 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -2,12 +2,10 @@ using OrdinaryDiffEq using Trixi # Physical constants -g::Float64 = 9.81 # gravity of earth -p_0::Float64 = 100_000.0 # reference pressure -c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) -c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) -R = c_p - c_v # gas constant (dry air) -gamma = c_p / c_v # heat capacity ratio (dry air) +g::Float64 = 9.81 # gravity of earth +c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) +c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) +gamma::Float64 = c_p / c_v # heat capacity ratio (dry air) # Warm bubble test case from # Wicker, L. J., and Skamarock, W. C. @@ -16,29 +14,31 @@ gamma = c_p / c_v # heat capacity ratio (dry air) # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D) # center of perturbation - xc = 10000.0 - zc = 2000.0 + center_x = 10000.0 + center_z = 2000.0 # radius of perturbation - rc = 2000.0 + radius = 2000.0 # distance of current x to center of perturbation - r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) # perturbation in potential temperature - θ_ref = 300.0 - Δθ = 0.0 - if r <= rc - Δθ = 2 * cospi(0.5 * r / rc)^2 + potential_temperature_ref = 300.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 end - θ = θ_ref + Δθ # potential temperature + potential_temperature = potential_temperature_ref + potential_temperature_perturbation # Exner pressure, solves hydrostatic equation for x[2] - exner = 1 - g / (c_p * θ) * x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) p = p_0 * exner^(c_p / R) # temperature - T = θ * exner + T = potential_temperature * exner # density rho = p / (R * T) From b6e527bd63767e0a124a274d57e986df9a88b766 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 19 Dec 2023 13:18:34 +0100 Subject: [PATCH 07/41] reduce run time of test case --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index a962aade6e1..d250136d5d3 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -77,7 +77,8 @@ solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (20_000.0, 10_000.0) -cells_per_dimension = (64, 32) +# increase resolution to (64, 32) +cells_per_dimension = (32, 16) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver, @@ -87,7 +88,8 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubb ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 1000.0) +# change final time to 1000.0 +tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan) @@ -103,7 +105,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) save_solution = SaveSolutionCallback(interval = analysis_interval, save_initial_solution = true, save_final_solution = true, - output_directory = "out.struct_lmars_ra", + output_directory = "out", solution_variables = cons2prim) stepsize_callback = StepsizeCallback(cfl = 0.2) From 7cfe952baa099842d8a359c8ff797c1bd06558a0 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 09:14:33 +0100 Subject: [PATCH 08/41] Revert "reduce run time of test case" This reverts commit b6e527bd63767e0a124a274d57e986df9a88b766. --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index d250136d5d3..a962aade6e1 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -77,8 +77,7 @@ solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (20_000.0, 10_000.0) -# increase resolution to (64, 32) -cells_per_dimension = (32, 16) +cells_per_dimension = (64, 32) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver, @@ -88,8 +87,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubb ############################################################################### # ODE solvers, callbacks etc. -# change final time to 1000.0 -tspan = (0.0, 10.0) +tspan = (0.0, 1000.0) ode = semidiscretize(semi, tspan) @@ -105,7 +103,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) save_solution = SaveSolutionCallback(interval = analysis_interval, save_initial_solution = true, save_final_solution = true, - output_directory = "out", + output_directory = "out.struct_lmars_ra", solution_variables = cons2prim) stepsize_callback = StepsizeCallback(cfl = 0.2) From 03db8f134a7d6beb7668b3ef87840592e4424eec Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 09:15:24 +0100 Subject: [PATCH 09/41] change output folder --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index a962aade6e1..24ba7d1aad0 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -103,7 +103,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) save_solution = SaveSolutionCallback(interval = analysis_interval, save_initial_solution = true, save_final_solution = true, - output_directory = "out.struct_lmars_ra", + output_directory = "out", solution_variables = cons2prim) stepsize_callback = StepsizeCallback(cfl = 0.2) From 78f22c31874b3037462d20d8f7896cfaf44e8162 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 09:16:00 +0100 Subject: [PATCH 10/41] change energy term in LMARS solver to use p_l/r --- src/equations/compressible_euler_2d.jl | 13 +++++++++---- src/equations/compressible_euler_3d.jl | 20 +++++++++++++------- 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 8abf365840b..f7c3b5c7fa2 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -847,10 +847,14 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + # We treat the energy term in analogy to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4 = v * u_ll + f4 = f4 + p_ll * v else f1, f2, f3, f4 = v * u_rr + f4 = f4 + p_rr * v end if orientation == 1 @@ -858,7 +862,6 @@ end else # orientation == 2 f3 = f3 + p end - f4 = f4 + p * v return SVector(f1, f2, f3, f4) end @@ -882,18 +885,20 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + # We treat the energy term in analogy to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4 = u_ll * v - # f4 += p_ll * v ??? + f4 = f4 + p_ll * v else f1, f2, f3, f4 = u_rr * v - # f4 += p_rr * v ??? + f4 = f4 + p_rr * v end return SVector(f1, f2 + p * normal_direction[1], f3 + p * normal_direction[2], - f4 + p * v) + f4) end """ diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 66d35edef43..2a2b0119d03 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -967,10 +967,14 @@ References: p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + # We treat the energy term in analogy to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll + f5 = f5 + p_ll * v else f1, f2, f3, f4, f5 = v * u_rr + f5 = f5 + p_rr * v end if orientation == 1 @@ -980,7 +984,6 @@ References: else # orientation == 3 f4 += p end - f5 += p * v return SVector(f1, f2, f3, f4, f5) end @@ -1006,18 +1009,21 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + # We treat the energy term in analogy to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll + f5 = f5 + p_ll * v else f1, f2, f3, f4, f5 = v * u_rr + f5 = f5 + p_rr * v end - f2 += p * normal_direction[1] - f3 += p * normal_direction[2] - f4 += p * normal_direction[3] - f5 += p * v - - return SVector(f1, f2, f3, f4, f5) + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4 + p * normal_direction[3], + f5) end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the From 700dd30b9236b2b900fee4b83fd16f617f480730 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 09:29:31 +0100 Subject: [PATCH 11/41] add lmars consistency checks --- test/test_unit.jl | 48 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index b3ed29d38e3..abb97c0b9ac 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1247,6 +1247,49 @@ end end end +@timed_testset "Consistency check for LMARS flux" begin + equations = CompressibleEulerEquations2D(1.4) + flux_lmars = FluxLMARS(340) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + orientations = [1, 2] + u_values = [SVector(1.0, 0.5, -0.7, 1.0), + SVector(1.5, -0.2, 0.1, 5.0)] + + for u in u_values, orientation in orientations + @test flux_lmars(u, u, orientation, equations) ≈ + flux_lmars(u, orientation, equations) + end + + for u in u_values, normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux_lmars(u, normal_direction, equations) + end + + equations = CompressibleEulerEquations3D(1.4) + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + orientations = [1, 2, 3] + u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), + SVector(1.5, -0.2, 0.1, 0.2, 5.0)] + + for u in u_values, orientation in orientations + @test flux_lmars(u, u, orientation, equations) ≈ + flux_lmars(u, orientation, equations) + end + + for u in u_values, normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux_lmars(u, normal_direction, equations) + end +end + @testset "FluxRotated vs. direct implementation" begin @timed_testset "CompressibleEulerMulticomponentEquations2D" begin equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), @@ -1280,7 +1323,7 @@ end u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] for f_std in fluxes f_rot = FluxRotated(f_std) @@ -1303,8 +1346,7 @@ end u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - FluxLMARS(340), - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] for f_std in fluxes f_rot = FluxRotated(f_std) From e25f2dc9a67823a84d160fa4f15b2e8b982b5bf8 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 12:25:37 +0100 Subject: [PATCH 12/41] switched to kennedy gruber flux --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index 24ba7d1aad0..f73e5681392 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -69,7 +69,7 @@ polydeg = 4 basis = LobattoLegendreBasis(polydeg) surface_flux = FluxLMARS(340.0) -volume_flux = flux_ranocha +volume_flux = flux_kennedy_gruber volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) From 115178960c0d78e645c151611cc84bed915e08cf Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 12:26:04 +0100 Subject: [PATCH 13/41] add euler warm bubble elixir to tests --- test/test_structured_2d.jl | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 1addc29e3e6..3a5b1c39f47 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -611,6 +611,33 @@ end end end +@trixi_testset "elixir_euler_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_warm_bubble.jl"), + l2=[ + 0.00019389374978502348, + 0.030867307039472824, + 0.04541800001263111, + 43.89530131084914, + ], + linf=[ + 0.0015906758532284737, + 0.17545022864421256, + 0.3729842228441566, + 307.76174927831744, + ], + cells_per_dimension=(32, 16), + tspan=(0.0, 10.)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), l2=[ From d91e0f61a28ea7f54d4bfe1132f4c84ac96f6185 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 13:16:58 +0100 Subject: [PATCH 14/41] adapt errors due to change flux --- test/test_p4est_3d.jl | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index f2467f30204..2ec8defbcdc 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -352,18 +352,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_circular_wind_nonconforming.jl"), l2=[ - 1.573832094977477e-7, - 3.863090659429634e-5, - 3.867293305754584e-5, - 3.686550296950078e-5, - 0.05508968493733932, + 1.5737711609657832e-7, + 3.8630261900166194e-5, + 3.8672287531936816e-5, + 3.6865116098660796e-5, + 0.05508620970403884, ], linf=[ - 2.2695202613887133e-6, - 0.0005314968179916946, - 0.0005314969614147458, - 0.0005130280733059617, - 0.7944959432352334, + 2.268845333053271e-6, + 0.000531462302113539, + 0.0005314624461298934, + 0.0005129931254772464, + 0.7942778058932163, ], tspan=(0.0, 2e2), coverage_override=(trees_per_cube_face = (1, 1), polydeg = 3)) # Prevent long compile time in CI @@ -381,18 +381,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_baroclinic_instability.jl"), l2=[ - 6.725065410642336e-7, - 0.00021710117340245454, - 0.000438679759422352, - 0.00020836356588024185, - 0.07602006689579247, + 6.725093801700048e-7, + 0.00021710076010951073, + 0.0004386796338203878, + 0.00020836270267103122, + 0.07601887903440395, ], linf=[ - 1.9101671995258585e-5, - 0.029803626911022396, - 0.04847630924006063, - 0.022001371349740104, - 4.847761006938526, + 1.9107530539574924e-5, + 0.02980358831035801, + 0.048476331898047564, + 0.02200137344113612, + 4.848310144356219, ], tspan=(0.0, 1e2), # Decrease tolerance of adaptive time stepping to get similar results across different systems From b38d57e74518cd5e2e80b0954658ded020ce286d Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 13:17:40 +0100 Subject: [PATCH 15/41] add warm bubble test with TreeMesh --- .../tree_2d_dgsem/elixir_euler_warm_bubble.jl | 126 ++++++++++++++++++ test/test_tree_2d_euler.jl | 26 ++++ 2 files changed, 152 insertions(+) create mode 100644 examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl new file mode 100644 index 00000000000..8395656872a --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -0,0 +1,126 @@ +using OrdinaryDiffEq +using Trixi + +# Physical constants +g::Float64 = 9.81 # gravity of earth +c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) +c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) +gamma::Float64 = c_p / c_v # heat capacity ratio (dry air) + +# Warm bubble test case from +# Wicker, L. J., and Skamarock, W. C. +# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta +# time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D) + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + +@inline function source_terms_gravity(u, x, t, + equations::CompressibleEulerEquations2D) + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(gamma) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = FluxLMARS(340.0) +volume_flux = flux_kennedy_gruber +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000, + periodicity = (true, false)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver, + source_terms = source_terms_gravity, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 65899cd5263..39aa82320e1 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -834,6 +834,32 @@ end end end +@trixi_testset "elixir_euler_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_warm_bubble.jl"), + l2=[ + 0.00013820060399128458, + 0.020833682516209616, + 0.03327372067821428, + 31.39884011181286, + ], + linf=[ + 0.0016677536924397662, + 0.15617627577194781, + 0.3368120143100738, + 331.40871663423604, + ], + tspan=(0.0, 10.0), + initial_refinement_level=4) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + # Coverage test for all initial conditions @testset "Compressible Euler: Tests for initial conditions" begin @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin From e373c08591902f9e4d0bc60c9bce54c2cba9299c Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 13:18:09 +0100 Subject: [PATCH 16/41] fix unit test --- test/test_unit.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index abb97c0b9ac..45e625cb760 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1261,12 +1261,12 @@ end for u in u_values, orientation in orientations @test flux_lmars(u, u, orientation, equations) ≈ - flux_lmars(u, orientation, equations) + flux(u, orientation, equations) end for u in u_values, normal_direction in normal_directions @test flux_lmars(u, u, normal_direction, equations) ≈ - flux_lmars(u, normal_direction, equations) + flux(u, normal_direction, equations) end equations = CompressibleEulerEquations3D(1.4) @@ -1281,12 +1281,12 @@ end for u in u_values, orientation in orientations @test flux_lmars(u, u, orientation, equations) ≈ - flux_lmars(u, orientation, equations) + flux(u, orientation, equations) end for u in u_values, normal_direction in normal_directions @test flux_lmars(u, u, normal_direction, equations) ≈ - flux_lmars(u, normal_direction, equations) + flux(u, normal_direction, equations) end end From 285f09b3008f5f212b00ab528cc11e57d4d1efef Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 21 Dec 2023 13:19:24 +0100 Subject: [PATCH 17/41] fix format --- test/test_structured_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 3a5b1c39f47..f6857a3df8a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -627,7 +627,7 @@ end 307.76174927831744, ], cells_per_dimension=(32, 16), - tspan=(0.0, 10.)) + tspan=(0.0, 10.0)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 63d66ce937b78b21b73cc865c7886e2d4f274f74 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 9 Jan 2024 13:41:50 +0100 Subject: [PATCH 18/41] adapt polynomial degree and CFL number --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 4 ++-- examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index f73e5681392..be52f674a4e 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -65,7 +65,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) -polydeg = 4 +polydeg = 3 basis = LobattoLegendreBasis(polydeg) surface_flux = FluxLMARS(340.0) @@ -106,7 +106,7 @@ save_solution = SaveSolutionCallback(interval = analysis_interval, output_directory = "out", solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.2) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl index 8395656872a..15c6e6eb52b 100644 --- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -65,7 +65,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) -polydeg = 4 +polydeg = 3 basis = LobattoLegendreBasis(polydeg) surface_flux = FluxLMARS(340.0) @@ -108,7 +108,7 @@ save_solution = SaveSolutionCallback(interval = analysis_interval, output_directory = "out", solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.2) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, From 82d1d7617929ae817a873b4a38f33b67f8cdabc4 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 9 Jan 2024 13:46:18 +0100 Subject: [PATCH 19/41] fix format --- test/test_unit.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index fff6d915033..e8a8effbe29 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1363,7 +1363,8 @@ end u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc, + ] for f_std in fluxes f_rot = FluxRotated(f_std) @@ -1386,7 +1387,8 @@ end u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc, + ] for f_std in fluxes f_rot = FluxRotated(f_std) From 17e11a21be992651731516e21e10c99774934a74 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 9 Jan 2024 14:23:54 +0100 Subject: [PATCH 20/41] adapt tests due to changed parameters --- test/test_structured_2d.jl | 16 ++++++++-------- test/test_tree_2d_euler.jl | 16 ++++++++-------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index f6857a3df8a..e5d45ebcc07 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -615,16 +615,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_warm_bubble.jl"), l2=[ - 0.00019389374978502348, - 0.030867307039472824, - 0.04541800001263111, - 43.89530131084914, + 0.00019387402388722496, + 0.03086514388623955, + 0.04541427917165, + 43.892826583444716, ], linf=[ - 0.0015906758532284737, - 0.17545022864421256, - 0.3729842228441566, - 307.76174927831744, + 0.0015942305974430138, + 0.17449778969139373, + 0.3729704262394843, + 307.6706958565337, ], cells_per_dimension=(32, 16), tspan=(0.0, 10.0)) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 39aa82320e1..3e3f69d36c3 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -837,16 +837,16 @@ end @trixi_testset "elixir_euler_warm_bubble.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_warm_bubble.jl"), l2=[ - 0.00013820060399128458, - 0.020833682516209616, - 0.03327372067821428, - 31.39884011181286, + 0.0001379946769624388, + 0.02078779689715382, + 0.033237241571263176, + 31.36068872331705, ], linf=[ - 0.0016677536924397662, - 0.15617627577194781, - 0.3368120143100738, - 331.40871663423604, + 0.0016286690573188434, + 0.15623770697198225, + 0.3341371832270615, + 331.5373488726036, ], tspan=(0.0, 10.0), initial_refinement_level=4) From 0c97c9f64118420d36c6df7b348039af14d68cb6 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Wed, 10 Jan 2024 09:05:53 +0100 Subject: [PATCH 21/41] Update src/equations/compressible_euler_2d.jl Co-authored-by: Andrew Winters --- src/equations/compressible_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index fa2c9f5e94f..2d100b23bf0 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -847,7 +847,7 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) - # We treat the energy term in analogy to the potential temperature term in the paper by + # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4 = v * u_ll From fc4da198db06bc01da53db43e801c40f912fd57a Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Wed, 10 Jan 2024 09:06:38 +0100 Subject: [PATCH 22/41] Update src/equations/compressible_euler_2d.jl Co-authored-by: Andrew Winters --- src/equations/compressible_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 2d100b23bf0..3c6f759db2b 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -885,7 +885,7 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ - # We treat the energy term in analogy to the potential temperature term in the paper by + # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4 = u_ll * v From 8b8f1e21827482c5e1a0301bb36525e8a79bd2f2 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Wed, 10 Jan 2024 09:08:04 +0100 Subject: [PATCH 23/41] Update src/equations/compressible_euler_3d.jl Co-authored-by: Andrew Winters --- src/equations/compressible_euler_3d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index db74f1f5485..54b36e0b832 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -967,7 +967,7 @@ References: p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) - # We treat the energy term in analogy to the potential temperature term in the paper by + # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll From 854eefda0bc2955938c6f7b79d77c610706b0fac Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Wed, 10 Jan 2024 09:08:19 +0100 Subject: [PATCH 24/41] Update src/equations/compressible_euler_3d.jl Co-authored-by: Andrew Winters --- src/equations/compressible_euler_3d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 54b36e0b832..292b912f009 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1009,7 +1009,7 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ - # We treat the energy term in analogy to the potential temperature term in the paper by + # We treat the energy term analogous to the potential temperature term in the paper by # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll From 2d55dc9c60033aad03784fa01fd3313418cbae04 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 10 Jan 2024 12:49:48 +0100 Subject: [PATCH 25/41] correct test result --- test/test_tree_2d_euler.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 3e3f69d36c3..61b5c54b5e9 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -846,7 +846,7 @@ end 0.0016286690573188434, 0.15623770697198225, 0.3341371832270615, - 331.5373488726036, + 334.5373488726036, ], tspan=(0.0, 10.0), initial_refinement_level=4) From 54d7de00213ec096fd0a910709d076ddcfd54737 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 10 Jan 2024 14:53:23 +0100 Subject: [PATCH 26/41] use callable struct to hold parameters Thanks sloede! --- .../elixir_euler_warm_bubble.jl | 35 ++++++++++++------ .../tree_2d_dgsem/elixir_euler_warm_bubble.jl | 37 +++++++++++++------ 2 files changed, 48 insertions(+), 24 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index be52f674a4e..ba1f23c8171 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -1,18 +1,27 @@ using OrdinaryDiffEq using Trixi -# Physical constants -g::Float64 = 9.81 # gravity of earth -c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) -c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) -gamma::Float64 = c_p / c_v # heat capacity ratio (dry air) - # Warm bubble test case from # Wicker, L. J., and Skamarock, W. C. # A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta # time differencing # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) -function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + @unpack g, c_p, c_v = setup + # center of perturbation center_x = 10000.0 center_z = 2000.0 @@ -49,16 +58,18 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio return SVector(rho, rho * v1, rho * v2, rho * E) end -@inline function source_terms_gravity(u, x, t, - equations::CompressibleEulerEquations2D) +# Source terms +@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack g = setup rho, _, rho_v2, _ = u return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) end ############################################################################### # semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() -equations = CompressibleEulerEquations2D(gamma) +equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) boundary_conditions = (x_neg = boundary_condition_periodic, x_pos = boundary_condition_periodic, @@ -80,8 +91,8 @@ coordinates_max = (20_000.0, 10_000.0) cells_per_dimension = (64, 32) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver, - source_terms = source_terms_gravity, +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + source_terms = warm_bubble_setup, boundary_conditions = boundary_conditions) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl index 15c6e6eb52b..1064da50232 100644 --- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -1,18 +1,27 @@ using OrdinaryDiffEq using Trixi -# Physical constants -g::Float64 = 9.81 # gravity of earth -c_p::Float64 = 1004.0 # heat capacity for constant pressure (dry air) -c_v::Float64 = 717.0 # heat capacity for constant volume (dry air) -gamma::Float64 = c_p / c_v # heat capacity ratio (dry air) - # Warm bubble test case from # Wicker, L. J., and Skamarock, W. C. # A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta # time differencing # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) -function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + @unpack g, c_p, c_v = setup + # center of perturbation center_x = 10000.0 center_z = 2000.0 @@ -49,16 +58,18 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio return SVector(rho, rho * v1, rho * v2, rho * E) end -@inline function source_terms_gravity(u, x, t, - equations::CompressibleEulerEquations2D) +# Source terms +@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack g = setup rho, _, rho_v2, _ = u return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) end ############################################################################### # semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() -equations = CompressibleEulerEquations2D(gamma) +equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) boundary_conditions = (x_neg = boundary_condition_periodic, x_pos = boundary_condition_periodic, @@ -77,13 +88,15 @@ solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (20_000.0, 10_000.0) +# Same coordinates as in examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +# However TreeMesh will generate a 20_000 x 20_000 square domain instead mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 6, n_cells_max = 10_000, periodicity = (true, false)) -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver, - source_terms = source_terms_gravity, +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + source_terms = warm_bubble_setup, boundary_conditions = boundary_conditions) ############################################################################### From 843173b344921d7388fe701112ce6a1d1ce19472 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Wed, 10 Jan 2024 16:50:25 +0100 Subject: [PATCH 27/41] Update examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl Co-authored-by: Hendrik Ranocha --- .../elixir_euler_warm_bubble.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index ba1f23c8171..06d187568ae 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -2,10 +2,18 @@ using OrdinaryDiffEq using Trixi # Warm bubble test case from -# Wicker, L. J., and Skamarock, W. C. -# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta -# time differencing -# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# - Wicker, L. J., and Skamarock, W. C. +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) struct WarmBubbleSetup # Physical constants g::Float64 # gravity of earth From 623dae344196b62f8ad43bc46769429cd85b26db Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Wed, 10 Jan 2024 16:53:32 +0100 Subject: [PATCH 28/41] Update examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl Co-authored-by: Hendrik Ranocha --- .../tree_2d_dgsem/elixir_euler_warm_bubble.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl index 1064da50232..867c27ab944 100644 --- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -2,10 +2,18 @@ using OrdinaryDiffEq using Trixi # Warm bubble test case from -# Wicker, L. J., and Skamarock, W. C. -# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta -# time differencing -# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# - Wicker, L. J., and Skamarock, W. C. +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) struct WarmBubbleSetup # Physical constants g::Float64 # gravity of earth From bed6975aa40cf9a03f947569ea1ef2d63283c695 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 11 Jan 2024 11:33:06 +0100 Subject: [PATCH 29/41] year of Wicker paper, comment on tspan [no ci] --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 4 ++-- examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index 06d187568ae..b26b3249e9e 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -2,7 +2,7 @@ using OrdinaryDiffEq using Trixi # Warm bubble test case from -# - Wicker, L. J., and Skamarock, W. C. +# - Wicker, L. J., and Skamarock, W. C. (1998) # A time-splitting scheme for the elastic equations incorporating # second-order Runge–Kutta time differencing # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) @@ -106,7 +106,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 1000.0) +tspan = (0.0, 1000.0) # 1000 seconds final time ode = semidiscretize(semi, tspan) diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl index 867c27ab944..e3054090291 100644 --- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -2,7 +2,7 @@ using OrdinaryDiffEq using Trixi # Warm bubble test case from -# - Wicker, L. J., and Skamarock, W. C. +# - Wicker, L. J., and Skamarock, W. C. (1998) # A time-splitting scheme for the elastic equations incorporating # second-order Runge–Kutta time differencing # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) @@ -110,7 +110,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 1000.0) +tspan = (0.0, 1000.0) # 1000 seconds final time ode = semidiscretize(semi, tspan) From cc43a5074cfb74fc9105be036f75315e8bea7c37 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 11 Jan 2024 18:12:42 +0100 Subject: [PATCH 30/41] add comment on speed of sound --- examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 3 +++ examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl | 3 +++ 2 files changed, 6 insertions(+) diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index b26b3249e9e..05c09d57530 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -87,7 +87,10 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 3 basis = LobattoLegendreBasis(polydeg) +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. surface_flux = FluxLMARS(340.0) + volume_flux = flux_kennedy_gruber volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl index e3054090291..f2e14273ae7 100644 --- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -87,7 +87,10 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 3 basis = LobattoLegendreBasis(polydeg) +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. surface_flux = FluxLMARS(340.0) + volume_flux = flux_kennedy_gruber volume_integral = VolumeIntegralFluxDifferencing(volume_flux) From 863f8b95397521f71a1ef7b5bc24d51c9d5436b7 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jan 2024 11:36:14 +0100 Subject: [PATCH 31/41] restructure CompressibleEulerMulticomponentEquations - introduce AbstractCompressibleEulerMulticomponentEquations{2} holding all common methods - introducde CompressibleMoistEulerEquations2D with methods that need to take into account latent heat - CompressibleEulerMulticomponentEquations{2} holds those methods in their original formulation --- .../elixir_eulermultimoist_warm_bubble.jl | 384 ++++++++++ src/Trixi.jl | 3 +- .../compressible_euler_multicomponent_2d.jl | 724 ++++-------------- ...ssible_euler_multicomponent_abstract_2d.jl | 641 ++++++++++++++++ ...pressible_euler_multicomponent_moist_2d.jl | 390 ++++++++++ src/equations/equations.jl | 32 +- test/test_structured_2d.jl | 26 + 7 files changed, 1612 insertions(+), 588 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl create mode 100644 src/equations/compressible_euler_multicomponent_abstract_2d.jl create mode 100644 src/equations/compressible_euler_multicomponent_moist_2d.jl diff --git a/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl new file mode 100644 index 00000000000..c0d8966af5e --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl @@ -0,0 +1,384 @@ + +using OrdinaryDiffEq +using Trixi +using NLsolve: nlsolve +using Infiltrator + +# Warm moist bubble test case from +# Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_pd::Float64 # heat capacity for constant pressure (dry air) + c_vd::Float64 # heat capacity for constant volume (dry air) + R_d::Float64 # gas constant (dry air) + c_pv::Float64 # heat capacity for constant pressure (water vapor) + c_vv::Float64 # heat capacity for constant volume (water vapor) + R_v::Float64 # gas constant (water vapor) + c_pl::Float64 # heat capacity for constant pressure (liquid water) + p_0::Float64 # reference pressure 1000 hPa + L_00::Float64 # latent heat of evaporation at 0 K + # TODO value ??? + # L00 = parameter("L00",2.5000e6 + (c_pl - c_pv) * 273.15) + + function WarmBubbleSetup(; g = 9.81, + c_pd = 1004.0, c_vd = 717.0, R_d = c_pd - c_vd, + c_pv = 1885, c_vv = 1424.0, R_v = c_pv - c_vv, + c_pl = 4186.0, p_0 = 100_000.0, L_00 = 3147620.0) + new(g, c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_0, L_00) + end +end + +function get_c_p(setup::WarmBubbleSetup) + return (setup.c_pd, setup.c_pv, setup.c_pl) +end + +function get_c_v(setup::WarmBubbleSetup) + return (setup.c_vd, setup.c_vv, setup.c_pl) # c_pl = c_vl for liquid phase +end + +############################################################################### +# source terms +@inline function source_terms_geopotential(u, setup::WarmBubbleSetup, + equations::CompressibleMoistEulerEquations2D) + @unpack g = setup + _, rho_v2 = u + rho = density(u, equations) + + return SVector(zero(eltype(u)), -g * rho, -g * rho_v2, + zero(eltype(u)), zero(eltype(u)), zero(eltype(u))) +end + +@inline function source_terms_phase_change(u, setup::WarmBubbleSetup, + equations::CompressibleMoistEulerEquations2D) + + # This source term models the phase chance between could water and vapor. + @unpack R_v = setup + rho_v1, rho_v2, rho_e, rho_qd, rho_qv, rho_ql = u + rho = density(u, equations) + T = temperature(u, equations) + + T_C = T - 273.15 + # saturation vapor pressure + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + + # saturation density of vapor + rho_star_qv = p_vs / (R_v * T) + + # Fisher-Burgmeister-Function + a = rho_star_qv - rho_qv + b = rho - rho_qv - rho_qd + + # saturation control factor + # < 1: stronger saturation effect + # > 1: weaker saturation effect + C = 1 + + Q_ph = (a + b - sqrt(a^2 + b^2)) * C + + return SVector(zero(eltype(u)), zero(eltype(u)), zero(eltype(u)), + zero(eltype(u)), Q_ph, -Q_ph) +end + +@inline function (setup::WarmBubbleSetup)(u, x, t, + equations::CompressibleMoistEulerEquations2D) + return source_terms_geopotential(u, setup, equations) + + source_terms_phase_change(u, setup, equations) +end + +############################################################################### +# Initial condition +struct AtmossphereLayers{RealT <: Real} + setup::WarmBubbleSetup + # structure: 1--> i-layer (z = total_hight/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql + LayerData::Matrix{RealT} + total_hight::RealT + preciseness::Int + layers::Int + ground_state::NTuple{2, RealT} + equivalentpotential_temperature::RealT + mixing_ratios::NTuple{2, RealT} +end + +function AtmossphereLayers(setup; total_hight = 10010.0, preciseness = 10, + ground_state = (1.4, 100000.0), + equivalentpotential_temperature = 320, + mixing_ratios = (0.02, 0.02), RealT = Float64) + @unpack c_pd, R_d, c_pv, R_v, c_pl = setup + rho0, p0 = ground_state + r_t0, r_v0 = mixing_ratios + theta_e0 = equivalentpotential_temperature + + rho_qv0 = rho0 * r_v0 + T0 = theta_e0 + y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] + + n = convert(Int, total_hight / preciseness) + dz = 0.01 + LayerData = zeros(RealT, n + 1, 4) + + F = generate_function_of_y(dz, y0, r_t0, theta_e0, setup) + sol = nlsolve(F, y0) + p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero + + rho_d = rho / (1 + r_t) + rho_ql = rho - rho_d - rho_qv + kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + LayerData[1, :] = [rho, rho_theta, rho_qv, rho_ql] + for i in (1:n) + #println("Layer", i) + y0 = deepcopy(sol.zero) + dz = preciseness + F = generate_function_of_y(dz, y0, r_t0, theta_e0, setup) + sol = nlsolve(F, y0) + p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero + rho_d = rho / (1 + r_t) + rho_ql = rho - rho_d - rho_qv + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + LayerData[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] + end + + return AtmossphereLayers{RealT}(setup, LayerData, total_hight, dz, n, ground_state, + theta_e0, mixing_ratios) +end + +function moist_state(y, dz, y0, r_t0, theta_e0, setup::WarmBubbleSetup) + @unpack g, c_pd, R_d, c_pv, R_v, c_pl, p_0, L_00 = setup + (p, rho, T, r_t, r_v, rho_qv, theta_e) = y + p0 = y0[1] + + F = zeros(7, 1) + rho_d = rho / (1 + r_t) + p_d = R_d * rho_d * T + T_C = T - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + L = L_00 - (c_pl - c_pv) * T + + F[1] = (p - p0) / dz + g * rho + F[2] = p - (R_d * rho_d + R_v * rho_qv) * T + # H = 1 is assumed + F[3] = (theta_e - + T * (p_d / p_0)^(-R_d / (c_pd + c_pl * r_t)) * + exp(L * r_v / ((c_pd + c_pl * r_t) * T))) + F[4] = r_t - r_t0 + F[5] = rho_qv - rho_d * r_v + F[6] = theta_e - theta_e0 + a = p_vs / (R_v * T) - rho_qv + b = rho - rho_qv - rho_d + # H=1 => phi=0 + F[7] = a + b - sqrt(a * a + b * b) + + return F +end + +function generate_function_of_y(dz, y0, r_t0, theta_e0, setup::WarmBubbleSetup) + function function_of_y(y) + return moist_state(y, dz, y0, r_t0, theta_e0, setup) + end +end + +# Moist bubble test case from paper: +# G.H. Bryan, J.M. Fritsch, A Benchmark Simulation for Moist Nonhydrostatic Numerical +# Models, MonthlyWeather Review Vol.130, 2917–2928, 2002, +# https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml. +function initial_condition_moist_bubble(x, t, setup::WarmBubbleSetup, + AtmosphereLayers::AtmossphereLayers) + @unpack LayerData, preciseness, total_hight = AtmosphereLayers + dz = preciseness + z = x[2] + if (z > total_hight && !(isapprox(z, total_hight))) + error("The atmossphere does not match the simulation domain") + end + n = convert(Int, floor((z + eps()) / dz)) + 1 + z_l = (n - 1) * dz + (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :] + z_r = n * dz + if (z_l == total_hight) + z_r = z_l + dz + n = n - 1 + end + (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = LayerData[n + 1, :] + rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz + rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / + dz + rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz + rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz + + rho, rho_e, rho_qv, rho_ql = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, + setup) + + v1 = 0.0 + v2 = 0.0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2) + + rho_qd = rho - rho_qv - rho_ql + + return SVector(rho_v1, rho_v2, rho_E, rho_qd, rho_qv, rho_ql) +end + +function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbleSetup) + @unpack c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_0, L_00 = setup + xc = 10000.0 + zc = 2000.0 + rc = 2000.0 + Δθ = 2.0 + + r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + p_loc = p_0 * (R_d * rho_theta / p_0)^(1 / (1 - kappa_M)) + T_loc = p_loc / (R_d * rho_d + R_v * rho_qv) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + + p_v = rho_qv * R_v * T_loc + p_d = p_loc - p_v + T_C = T_loc - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + H = p_v / p_vs + r_v = rho_qv / rho_d + r_l = rho_ql / rho_d + r_t = r_v + r_l + + # Aequivalentpotential temperature + a = T_loc * (p_0 / p_d)^(R_d / (c_pd + r_t * c_pl)) + b = H^(-r_v * R_v / c_pd) + L_v = L_00 + (c_pv - c_pl) * T_loc + c = exp(L_v * r_v / ((c_pd + r_t * c_pl) * T_loc)) + aeq_pot = (a * b * c) + + # Assume pressure stays constant + if (r < rc && Δθ > 0) + kappa = 1 - c_vd / c_pd + # Calculate background density potential temperature + θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) + # Calculate perturbed density potential temperature + θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5 * r / rc)^2 / 300) + rt = (rho_qv + rho_ql) / rho_d + rv = rho_qv / rho_d + # Calculate moist potential temperature + θ_loc = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rv) + # Adjust varuables until the temperature is met + if rt > 0 + while true + T_loc = θ_loc * (p_loc / p_0)^kappa + T_C = T_loc - 273.15 + # SaturVapor + pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + rho_d_new = (p_loc - pvs) / (R_d * T_loc) + rvs = pvs / (R_v * rho_d_new * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + if abs(θ_new - θ_loc) <= θ_loc * 1.0e-12 + break + else + θ_loc = θ_new + end + end + else + rvs = 0 + T_loc = θ_loc * (p_loc / p_0)^kappa + rho_d_new = p_loc / (R_d * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + end + rho_qv = rvs * rho_d_new + rho_ql = (rt - rvs) * rho_d_new + rho = rho_d_new * (1 + rt) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + end + return SVector(rho, rho_e, rho_qv, rho_ql) +end + +############################################################################### +# semidiscretization of the compressible moist Euler equations +warm_bubble_setup = WarmBubbleSetup() + +AtmossphereData = AtmossphereLayers(warm_bubble_setup) + +# Create the initial condition with the initial data set +function initial_condition_moist(x, t, equations) + return initial_condition_moist_bubble(x, t, warm_bubble_setup, AtmossphereData) +end + +equations = CompressibleMoistEulerEquations2D(; c_p = get_c_p(warm_bubble_setup), + c_v = get_c_v(warm_bubble_setup), + L_00 = warm_bubble_setup.L_00) + +boundary_condition = (x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = FluxLMARS(360.0) +volume_flux = flux_chandrashekar +#volume_flux = flux_ranocha +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20000.0, 10000.0) + +cells_per_dimension = (64, 32) + +# Create curved mesh with 64 x 32 elements +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_moist, solver, + boundary_conditions = boundary_condition, + source_terms = warm_bubble_setup) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +#solution_variables = cons2aeqpot +solution_variables = cons2prim + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) +#analysis_callback = AnalysisCallback(semi, interval=analysis_interval, extra_analysis_errors=(:entropy_conservation_error, ), extra_analysis_integrals=(entropy, energy_total, Trixi.saturation_pressure)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = solution_variables) + +stepsize_callback = StepsizeCallback(cfl = 0.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index e18b2f6415c..1b46224f777 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -141,6 +141,7 @@ export AcousticPerturbationEquations2D, CompressibleEulerEquations3D, CompressibleEulerMulticomponentEquations1D, CompressibleEulerMulticomponentEquations2D, + CompressibleMoistEulerEquations2D, CompressibleEulerEquationsQuasi1D, IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D, IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D, @@ -213,7 +214,7 @@ export initial_condition_eoc_test_coupled_euler_gravity, export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, cons2entropy, entropy2cons export density, pressure, density_pressure, velocity, global_mean_vars, - equilibrium_distribution, waterheight_pressure + equilibrium_distribution, waterheight_pressure, temperature export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, enstrophy diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 60fce222f21..1e7a6459600 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -102,191 +102,164 @@ end RealT end -function varnames(::typeof(cons2cons), - equations::CompressibleEulerMulticomponentEquations2D) - cons = ("rho_v1", "rho_v2", "rho_e") - rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) - return (cons..., rhos...) -end - -function varnames(::typeof(cons2prim), - equations::CompressibleEulerMulticomponentEquations2D) - prim = ("v1", "v2", "p") - rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) - return (prim..., rhos...) -end - -# Set initial conditions at physical location `x` for time `t` - """ - initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D) + temperature(u, equations::CompressibleEulerMulticomponentEquations2D) -A smooth initial condition used for convergence tests in combination with -[`source_terms_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +Calculate temperature. """ -function initial_condition_convergence_test(x, t, - equations::CompressibleEulerMulticomponentEquations2D) - c = 2 - A = 0.1 - L = 2 - f = 1 / L - omega = 2 * pi * f - ini = c + A * sin(omega * (x[1] + x[2] - t)) - - v1 = 1.0 - v2 = 1.0 - - rho = ini - - # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1) - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho - for i in eachcomponent(equations)) - - prim1 = rho * v1 - prim2 = rho * v2 - prim3 = rho^2 +@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations - prim_other = SVector{3, real(equations)}(prim1, prim2, prim3) + rho_v1, rho_v2, rho_e, rho_d, rho_v = u - return vcat(prim_other, prim_rho) -end + rho = density(u, equations) + help1 = zero(rho) -""" - source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D) + # compute weighted average of cv + # normalization by rho not required, cancels below + for i in eachcomponent(equations) + help1 += u[i + 3] * cv[i] + end -Source terms used for convergence tests in combination with -[`initial_condition_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). -""" -@inline function source_terms_convergence_test(u, x, t, - equations::CompressibleEulerMulticomponentEquations2D) - # Same settings as in `initial_condition` - c = 2 - A = 0.1 - L = 2 - f = 1 / L - omega = 2 * pi * f - - gamma = totalgamma(u, equations) - - x1, x2 = x - si, co = sincos((x1 + x2 - t) * omega) - tmp1 = co * A * omega - tmp2 = si * A - tmp3 = gamma - 1 - tmp4 = (2 * c - 1) * tmp3 - tmp5 = (2 * tmp2 * gamma - 2 * tmp2 + tmp4 + 1) * tmp1 - tmp6 = tmp2 + c - - # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1 - du_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - tmp1 - for i in eachcomponent(equations)) - - du1 = tmp5 - du2 = tmp5 - du3 = 2 * ((tmp6 - 1.0) * tmp3 + tmp6 * gamma) * tmp1 - - du_other = SVector{3, real(equations)}(du1, du2, du3) - - return vcat(du_other, du_rho) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + T = (rho_e - 0.5 * rho * v_square) / help1 + return T end """ - initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D) + totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) -A for multicomponent adapted weak blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +Function that calculates the total gamma out of all partial gammas using the +partial density fractions as well as the partial specific heats at constant volume. """ -function initial_condition_weak_blast_wave(x, t, - equations::CompressibleEulerMulticomponentEquations2D) - # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) - # Set up polar coordinates - inicenter = SVector(0.0, 0.0) - x_norm = x[1] - inicenter[1] - y_norm = x[2] - inicenter[2] - r = sqrt(x_norm^2 + y_norm^2) - phi = atan(y_norm, x_norm) - sin_phi, cos_phi = sincos(phi) - - prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? - 2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - 1.0 : - 2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - 1.1691 - for i in eachcomponent(equations)) +@inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, cp = equations - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0 : 1.245 + help1 = zero(u[1]) + help2 = zero(u[1]) - prim_other = SVector{3, real(equations)}(v1, v2, p) + # compute weighted averages of cp and cv + # normalization by total rho not required, would cancel below + for i in eachcomponent(equations) + help1 += u[i + 3] * cp[i] + help2 += u[i + 3] * cv[i] + end - return prim2cons(vcat(prim_other, prim_rho), equations) + return help1 / help2 end -# Calculate 1D flux for a single point -@inline function flux(u, orientation::Integer, - equations::CompressibleEulerMulticomponentEquations2D) +""" +cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + +Convert conservative variables to entropy. +""" +@inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations rho_v1, rho_v2, rho_e = u rho = density(u, equations) + # Multicomponent stuff + help1 = zero(rho) + gas_constant = zero(rho) + for i in eachcomponent(equations) + help1 += u[i + 3] * cv[i] + gas_constant += gas_constants[i] * (u[i + 3] / rho) + end + v1 = rho_v1 / rho v2 = rho_v2 / rho - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + v_square = v1^2 + v2^2 + p = pressure(u, equations) + rho_p = rho / p + T = temperature(u, equations) - if orientation == 1 - f_rho = densities(u, v1, equations) - f1 = rho_v1 * v1 + p - f2 = rho_v2 * v1 - f3 = (rho_e + p) * v1 - else - f_rho = densities(u, v2, equations) - f1 = rho_v1 * v2 - f2 = rho_v2 * v2 + p - f3 = (rho_e + p) * v2 - end + entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 3])) - + v_square / (2 * T)) + for i in eachcomponent(equations)) - f_other = SVector{3, real(equations)}(f1, f2, f3) + w1 = gas_constant * v1 * rho_p + w2 = gas_constant * v2 * rho_p + w3 = gas_constant * (-rho_p) - return vcat(f_other, f_rho) + entrop_other = SVector{3, real(equations)}(w1, w2, w3) + + return vcat(entrop_other, entrop_rho) end -# Calculate 1D flux for a single point -@inline function flux(u, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u +""" +entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) - rho = density(u, equations) +Convert entropy variables to conservative variables +""" +@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) + @unpack gammas, gas_constants, cp, cv = equations + T = -1 / w[3] + v1 = w[1] * T + v2 = w[2] * T + v_squared = v1^2 + v2^2 + cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - + cv[i] * + (1 - log(T)) + + v_squared / + (2 * T)) / + gas_constants[i] - + 1) + for i in eachcomponent(equations)) - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + u2 = rho * v2 + gamma = help1 / help2 + u3 = p / (gamma - 1) + 0.5 * rho * v_squared + cons_other = SVector{3, real(equations)}(u1, u2, u3) + return vcat(cons_other, cons_rho) +end - f_rho = densities(u, v_normal, equations) - f1 = rho_v1 * v_normal + p * normal_direction[1] - f2 = rho_v2 * v_normal + p * normal_direction[2] - f3 = (rho_e + p) * v_normal +""" + total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) - f_other = SVector{3, real(equations)}(f1, f2, f3) +Calculate total entropy. +""" +@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations + T = temperature(u, equations) - return vcat(f_other, f_rho) + total_entropy = zero(u[1]) + for i in eachcomponent(equations) + total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) + end + + return total_entropy +end + +""" + density_gas_constant(u, equations::CompressibleEulerMulticomponentEquations2D{2}) + +Function that calculates overall density times overall gas constant. +""" +@inline function density_gas_constant(u, + equations::CompressibleEulerMulticomponentEquations2D) + @unpack gas_constants = equations + help = zero(u[1]) + for i in eachcomponent(equations) + help += u[i + 3] * gas_constant[i] + end + return help end """ @@ -326,17 +299,12 @@ Adaption of the entropy conserving two-point flux by v_sum = v1_avg + v2_avg enth = zero(v_sum) - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - for i in eachcomponent(equations) enth += rhok_avg[i] * gas_constants[i] - help1_ll += u_ll[i + 3] * cv[i] - help1_rr += u_rr[i + 3] * cv[i] end - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) @@ -371,27 +339,12 @@ Adaption of the entropy conserving two-point flux by return vcat(f_other, f_rho) end -""" - flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, - equations::CompressibleEulerMulticomponentEquations2D) - -Adaption of the entropy conserving and kinetic energy preserving two-point flux by -- Hendrik Ranocha (2018) - Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods - for Hyperbolic Balance Laws - [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743) -See also -- Hendrik Ranocha (2020) - Entropy Conserving and Kinetic Energy Preserving Numerical Methods for - the Euler Equations Using Summation-by-Parts Operators - [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42) -""" -@inline function flux_ranocha(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMulticomponentEquations2D) +@inline function flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) # Unpack left and right state @unpack gammas, gas_constants, cv = equations - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], u_rr[i + 3]) for i in eachcomponent(equations)) @@ -403,414 +356,45 @@ See also rho_ll = density(u_ll, equations) rho_rr = density(u_rr, equations) - # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) - inv_gamma_minus_one = 1 / (gamma - 1) - # extract velocities v1_ll = rho_v1_ll / rho_ll - v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) v2_ll = rho_v2_ll / rho_ll - v2_rr = rho_v2_rr / rho_rr - v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) - - # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) - for i in eachcomponent(equations) - enth_ll += u_ll[i + 3] * gas_constants[i] - enth_rr += u_rr[i + 3] * gas_constants[i] - help1_ll += u_ll[i + 3] * cv[i] - help1_rr += u_rr[i + 3] * cv[i] - end - - # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr - p_ll = T_ll * enth_ll - p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) - inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - - f_rho_sum = zero(T_rr) - if orientation == 1 - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg - for i in eachcomponent(equations)) - for i in eachcomponent(equations) - f_rho_sum += f_rho[i] - end - f1 = f_rho_sum * v1_avg + p_avg - f2 = f_rho_sum * v2_avg - f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v1_rr + p_rr * v1_ll) - else - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg - for i in eachcomponent(equations)) - for i in eachcomponent(equations) - f_rho_sum += f_rho[i] - end - f1 = f_rho_sum * v1_avg - f2 = f_rho_sum * v2_avg + p_avg - f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v2_rr + p_rr * v2_ll) - end - - # momentum and energy flux - f_other = SVector{3, real(equations)}(f1, f2, f3) - - return vcat(f_other, f_rho) -end - -@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - # Unpack left and right state - @unpack gammas, gas_constants, cv = equations - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], - u_rr[i + 3]) - for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + - u_rr[i + 3]) - for i in eachcomponent(equations)) - - # Iterating over all partial densities - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - - # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) - inv_gamma_minus_one = 1 / (gamma - 1) - - # extract velocities - v1_ll = rho_v1_ll / rho_ll v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_ll = rho_v2_ll / rho_ll v2_rr = rho_v2_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) - v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] - v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] - - # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) - for i in eachcomponent(equations) - enth_ll += u_ll[i + 3] * gas_constants[i] - enth_rr += u_rr[i + 3] * gas_constants[i] - help1_ll += u_ll[i + 3] * cv[i] - help1_rr += u_rr[i + 3] * cv[i] - end - - # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr - p_ll = T_ll * enth_ll - p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) - inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - - f_rho_sum = zero(T_rr) - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * - (v_dot_n_ll + v_dot_n_rr) - for i in eachcomponent(equations)) - for i in eachcomponent(equations) - f_rho_sum += f_rho[i] - end - f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] - f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] - f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) - - # momentum and energy flux - f_other = SVector(f1, f2, f3) - - return vcat(f_other, f_rho) -end - -# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation -@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - - # Get the density and gas gamma - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - gamma_ll = totalgamma(u_ll, equations) - gamma_rr = totalgamma(u_rr, equations) - - # Get the velocities based on direction - if orientation == 1 - v_ll = rho_v1_ll / rho_ll - v_rr = rho_v1_rr / rho_rr - else # orientation == 2 - v_ll = rho_v2_ll / rho_ll - v_rr = rho_v2_rr / rho_rr - end - - # Compute the sound speeds on the left and right - p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) - c_ll = sqrt(gamma_ll * p_ll / rho_ll) - p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) - c_rr = sqrt(gamma_rr * p_rr / rho_rr) - - λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) -end - -@inline function max_abs_speeds(u, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - v1 = rho_v1 / rho - v2 = rho_v2 / rho - - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 1 / 2 * rho * (v1^2 + v2^2)) - c = sqrt(gamma * p / rho) - - return (abs(v1) + c, abs(v2) + c) -end - -@inline function rotate_to_x(u, normal_vector, - equations::CompressibleEulerMulticomponentEquations2D) - # cos and sin of the angle between the x-axis and the normalized normal_vector are - # the normalized vector's x and y coordinates respectively (see unit circle). - c = normal_vector[1] - s = normal_vector[2] - - # Apply the 2D rotation matrix with normal and tangent directions of the form - # [ n_1 n_2 0 0; - # t_1 t_2 0 0; - # 0 0 1 0 - # 0 0 0 1] - # where t_1 = -n_2 and t_2 = n_1 - - densities = @view u[4:end] - return SVector(c * u[1] + s * u[2], - -s * u[1] + c * u[2], - u[3], - densities...) -end - -# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction -# has been normalized prior to this back-rotation of the state vector -@inline function rotate_from_x(u, normal_vector, - equations::CompressibleEulerMulticomponentEquations2D) - # cos and sin of the angle between the x-axis and the normalized normal_vector are - # the normalized vector's x and y coordinates respectively (see unit circle). - c = normal_vector[1] - s = normal_vector[2] - - # Apply the 2D back-rotation matrix with normal and tangent directions of the form - # [ n_1 t_1 0 0; - # n_2 t_2 0 0; - # 0 0 1 0; - # 0 0 0 1 ] - # where t_1 = -n_2 and t_2 = n_1 - - densities = @view u[4:end] - return SVector(c * u[1] - s * u[2], - s * u[1] + c * u[2], - u[3], - densities...) -end - -# Convert conservative variables to primitive -@inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u - - prim_rho = SVector{ncomponents(equations), real(equations)}(u[i + 3] - for i in eachcomponent(equations)) - - rho = density(u, equations) - v1 = rho_v1 / rho - v2 = rho_v2 / rho - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) - prim_other = SVector{3, real(equations)}(v1, v2, p) - - return vcat(prim_other, prim_rho) -end - -# Convert conservative variables to entropy -@inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas, gas_constants = equations - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - - # Multicomponent stuff - help1 = zero(rho) - gas_constant = zero(rho) - for i in eachcomponent(equations) - help1 += u[i + 3] * cv[i] - gas_constant += gas_constants[i] * (u[i + 3] / rho) - end - - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_square = v1^2 + v2^2 - gamma = totalgamma(u, equations) - - p = (gamma - 1) * (rho_e - 0.5 * rho * v_square) - s = log(p) - gamma * log(rho) - log(gas_constant) - rho_p = rho / p - T = (rho_e - 0.5 * rho * v_square) / (help1) - - entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * - (1 - log(T)) + - gas_constants[i] * - (1 + log(u[i + 3])) - - v_square / (2 * T)) - for i in eachcomponent(equations)) - - w1 = gas_constant * v1 * rho_p - w2 = gas_constant * v2 * rho_p - w3 = gas_constant * (-rho_p) - - entrop_other = SVector{3, real(equations)}(w1, w2, w3) - - return vcat(entrop_other, entrop_rho) -end - -# Convert entropy variables to conservative variables -@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) - @unpack gammas, gas_constants, cp, cv = equations - T = -1 / w[3] - v1 = w[1] * T - v2 = w[2] * T - v_squared = v1^2 + v2^2 - cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - - cv[i] * - (1 - log(T)) + - v_squared / - (2 * T)) / - gas_constants[i] - - 1) - for i in eachcomponent(equations)) - - rho = zero(cons_rho[1]) - help1 = zero(cons_rho[1]) - help2 = zero(cons_rho[1]) - p = zero(cons_rho[1]) - for i in eachcomponent(equations) - rho += cons_rho[i] - help1 += cons_rho[i] * cv[i] * gammas[i] - help2 += cons_rho[i] * cv[i] - p += cons_rho[i] * gas_constants[i] * T - end - u1 = rho * v1 - u2 = rho * v2 - gamma = help1 / help2 - u3 = p / (gamma - 1) + 0.5 * rho * v_squared - cons_other = SVector{3, real(equations)}(u1, u2, u3) - return vcat(cons_other, cons_rho) -end - -# Convert primitive to conservative variables -@inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas = equations - v1, v2, p = prim - - cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i + 3] - for i in eachcomponent(equations)) - rho = density(prim, equations) - gamma = totalgamma(prim, equations) - - rho_v1 = rho * v1 - rho_v2 = rho * v2 - rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2) - - cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e) - - return vcat(cons_other, cons_rho) -end - -@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas, gas_constants = equations - rho = density(u, equations) - T = temperature(u, equations) - - total_entropy = zero(u[1]) - for i in eachcomponent(equations) - total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) - end - - return total_entropy -end - -@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas, gas_constants = equations - - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - help1 = zero(rho) + v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg + v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v_sum = v1_avg + v2_avg + enth = zero(v_sum) for i in eachcomponent(equations) - help1 += u[i + 3] * cv[i] + enth += rhok_avg[i] * gas_constants[i] end - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_square = v1^2 + v2^2 - T = (rho_e - 0.5 * rho * v_square) / help1 - - return T -end - -""" - totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) - -Function that calculates the total gamma out of all partial gammas using the -partial density fractions as well as the partial specific heats at constant volume. -""" -@inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas = equations + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) - help1 = zero(u[1]) - help2 = zero(u[1]) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v_dot_n_avg + for i in eachcomponent(equations)) + help1 = zero(T_ll) + help2 = zero(T_rr) for i in eachcomponent(equations) - help1 += u[i + 3] * cv[i] * gammas[i] - help2 += u[i + 3] * cv[i] + help1 += f_rho[i] * cv[i] + help2 += f_rho[i] end - return help1 / help2 -end - -@inline function density_pressure(u, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - gamma = totalgamma(u, equations) - rho_times_p = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) - - return rho_times_p -end - -@inline function density(u, equations::CompressibleEulerMulticomponentEquations2D) - rho = zero(u[1]) + f1 = (help2) * v1_avg + normal_direction[1] * enth / T + f2 = (help2) * v2_avg + normal_direction[2] * enth / T + f3 = ((help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + + v1_avg * f1 + v2_avg * f2) - for i in eachcomponent(equations) - rho += u[i + 3] - end - - return rho -end + f_other = SVector{3, real(equations)}(f1, f2, f3) -@inline function densities(u, v, equations::CompressibleEulerMulticomponentEquations2D) - return SVector{ncomponents(equations), real(equations)}(u[i + 3] * v - for i in eachcomponent(equations)) + return vcat(f_other, f_rho) end end # @muladd diff --git a/src/equations/compressible_euler_multicomponent_abstract_2d.jl b/src/equations/compressible_euler_multicomponent_abstract_2d.jl new file mode 100644 index 00000000000..889b380b582 --- /dev/null +++ b/src/equations/compressible_euler_multicomponent_abstract_2d.jl @@ -0,0 +1,641 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function varnames(::typeof(cons2cons), + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + cons = ("rho_v1", "rho_v2", "rho_e") + rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) + return (cons..., rhos...) +end + +function varnames(::typeof(cons2prim), + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + prim = ("v1", "v2", "p") + rhos = ntuple(n -> "rho" * string(n), Val(ncomponents(equations))) + return (prim..., rhos...) +end + +# Calculate 1D flux for a single point +@inline function flux(u, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = pressure(u, equations) + + if orientation == 1 + f_rho = u[4:end] .* v1 + f1 = rho_v1 * v1 + p + f2 = rho_v2 * v1 + f3 = (rho_e + p) * v1 + else + f_rho = u[4:end] .* v2 + f1 = rho_v1 * v2 + f2 = rho_v2 * v2 + p + f3 = (rho_e + p) * v2 + end + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +# Calculate 1D flux for a single point +@inline function flux(u, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + p = pressure(u, equations) + + f_rho = u[4:end] .* v_normal + f1 = rho_v1 * v_normal + p * normal_direction[1] + f2 = rho_v2 * v_normal + p * normal_direction[2] + f3 = (rho_e + p) * v_normal + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function rotate_to_x(u, normal_vector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D rotation matrix with normal and tangent directions of the form + # [ n_1 n_2 0 0; + # t_1 t_2 0 0; + # 0 0 1 0 + # 0 0 0 1] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] + s * u[2], + -s * u[1] + c * u[2], + u[3], + densities...) +end + +# Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction +# has been normalized prior to this back-rotation of the state vector +@inline function rotate_from_x(u, normal_vector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D back-rotation matrix with normal and tangent directions of the form + # [ n_1 t_1 0 0; + # n_2 t_2 0 0; + # 0 0 1 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + densities = @view u[4:end] + return SVector(c * u[1] - s * u[2], + s * u[1] + c * u[2], + u[3], + densities...) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Determine the boundary numerical surface flux for a slip wall condition. +Imposes a zero normal velocity at the wall. +Density is taken from the internal solution state and pressure is computed as an +exact solution of a 1D Riemann problem. Further details about this boundary state +are available in the paper: +- J. J. W. van der Vegt and H. van der Ven (2002) + Slip flow boundary conditions in discontinuous Galerkin discretizations of + the Euler equations of gas dynamics + [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1) + +Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of the book +- Eleuterio F. Toro (2009) + Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + 3rd edition + [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + +Should be used together with [`UnstructuredMesh2D`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, surface_flux_function, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal = normal_direction / norm_ + + # rotate the internal solution state + u_local = rotate_to_x(u_inner, normal, equations) + + # compute the primitive variables + v_normal, v_tangent, p_local = cons2prim(u_local, equations) + + rho_local = density(u_local, equations) + gamma = totalgamma(u_inner, equations) + + # Get the solution of the pressure Riemann problem + if v_normal <= 0.0 + sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2 * gamma * + inv(gamma - 1)) + else # v_normal > 0.0 + A = 2 / ((gamma + 1) * rho_local) + B = p_local * (gamma - 1) / (gamma + 1) + p_star = p_local + + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + return SVector(p_star * normal[1], + p_star * normal[2], + zeros(eltype(u_inner), ncomponents(equations) + 1)...) * norm_ +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, + surface_flux_function, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Should be used together with [`StructuredMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x and -y sides due to the orientation convention used by StructuredMesh + if isodd(direction) + boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, + equations) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, + equations) + end + + return boundary_flux +end + +# Convert conservative variables to primitive +@inline function cons2prim(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + prim_rho = SVector{ncomponents(equations), real(equations)}(u[i + 3] + for i in eachcomponent(equations)) + + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = pressure(u, equations) + prim_other = SVector{3, real(equations)}(v1, v2, p) + + return vcat(prim_other, prim_rho) +end + +# Convert primitive to conservative variables +@inline function prim2cons(prim, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + v1, v2, p = prim + + cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i + 3] + for i in eachcomponent(equations)) + rho = density(prim, equations) + gamma = totalgamma(prim, equations) + + rho_v1 = rho * v1 + rho_v2 = rho * v2 + + rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + + cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e) + + return vcat(cons_other, cons_rho) +end + +""" + density(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Function that calculates the overall density. +""" +@inline function density(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho = zero(u[1]) + for i in eachcomponent(equations) + rho += u[i + 3] + end + return rho +end + +""" + totalgamma(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Function that calculates the total gamma out of all partial gammas. +This function has to be implemented for subtypes! +""" +function totalgamma end + +""" + temperature(u, equations::CompressibleEulerMulticomponentEquations2D) + +Calculate temperature. +This function has to be implemented for subtypes! +""" +function temperature end + +""" + pressure(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Function that calculates the overall pressure. +""" +@inline function pressure(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + gamma = totalgamma(u, equations) + + p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + return p +end + +""" + density_gas_constant(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Function that calculates overall density times overall gas constant. +This has to be implemented for subtypes! +""" +function density_gas_constant end + +""" + initial_condition_convergence_test(x, t, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function initial_condition_convergence_test(x, t, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + c = 2 + A = 0.1 + L = 2 + f = 1 / L + omega = 2 * pi * f + ini = c + A * sin(omega * (x[1] + x[2] - t)) + + v1 = 1.0 + v2 = 1.0 + + rho = ini + + # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + rho + for i in eachcomponent(equations)) + + prim1 = rho * v1 + prim2 = rho * v2 + prim3 = rho^2 + + prim_other = SVector{3, real(equations)}(prim1, prim2, prim3) + + return vcat(prim_other, prim_rho) +end + +""" + source_terms_convergence_test(u, x, t, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +@inline function source_terms_convergence_test(u, x, t, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Same settings as in `initial_condition` + c = 2 + A = 0.1 + L = 2 + f = 1 / L + omega = 2 * pi * f + + gamma = totalgamma(u, equations) + + x1, x2 = x + si, co = sincos((x1 + x2 - t) * omega) + tmp1 = co * A * omega + tmp2 = si * A + tmp3 = gamma - 1 + tmp4 = (2 * c - 1) * tmp3 + tmp5 = (2 * tmp2 * gamma - 2 * tmp2 + tmp4 + 1) * tmp1 + tmp6 = tmp2 + c + + # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1 + du_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + tmp1 + for i in eachcomponent(equations)) + + du1 = tmp5 + du2 = tmp5 + du3 = 2 * ((tmp6 - 1.0) * tmp3 + tmp6 * gamma) * tmp1 + + du_other = SVector{3, real(equations)}(du1, du2, du3) + + return vcat(du_other, du_rho) +end + +""" + initial_condition_weak_blast_wave(x, t, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +A for multicomponent adapted weak blast wave taken from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_weak_blast_wave(x, t, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? + 2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + 1.0 : + 2^(i - 1) * (1 - 2) / + (1 - + 2^ncomponents(equations)) * + 1.1691 + for i in eachcomponent(equations)) + + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi + p = r > 0.5 ? 1.0 : 1.245 + + prim_other = SVector{3, real(equations)}(v1, v2, p) + + return prim2cons(vcat(prim_other, prim_rho), equations) +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + + # Get the density and gas gamma + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + gamma_ll = totalgamma(u_ll, equations) + gamma_rr = totalgamma(u_rr, equations) + + # Get the velocities based on direction + if orientation == 1 + v_ll = rho_v1_ll / rho_ll + v_rr = rho_v1_rr / rho_rr + else # orientation == 2 + v_ll = rho_v2_ll / rho_ll + v_rr = rho_v2_rr / rho_rr + end + + # Compute the sound speeds on the left and right + p_ll = pressure(u_ll, equations) + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + p_rr = pressure(u_rr, equations) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function max_abs_speeds(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2 = u + + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + gamma = totalgamma(u, equations) + p = pressure(u, equations) + c = sqrt(gamma * p / rho) + + return (abs(v1) + c, abs(v2) + c) +end + +""" +flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Adaption of the entropy conserving and kinetic energy preserving two-point flux by +- Hendrik Ranocha (2018) +Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods +for Hyperbolic Balance Laws +[PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743) +See also +- Hendrik Ranocha (2020) +Entropy Conserving and Kinetic Energy Preserving Numerical Methods for +the Euler Equations Using Summation-by-Parts Operators +[Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42) +""" +@inline function flux_ranocha(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + v2_avg = 0.5 * (v2_ll + v2_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # helpful variables + enth_ll = density_gas_constant(u_ll) + enth_rr = density_gas_constant(u_rr) + + # temperature and pressure + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + p_ll = T_ll * enth_ll + p_rr = T_rr * enth_rr + p_avg = 0.5 * (p_ll + p_rr) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + + f_rho_sum = zero(T_rr) + if orientation == 1 + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + p_avg + f2 = f_rho_sum * v2_avg + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + else + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + f2 = f_rho_sum * v2_avg + p_avg + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + end + + # momentum and energy flux + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + v2_avg = 0.5 * (v2_ll + v2_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # helpful variables + enth_ll = density_gas_constant(u_ll) + enth_rr = density_gas_constant(u_rr) + + # temperature and pressure + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + p_ll = T_ll * enth_ll + p_rr = T_rr * enth_rr + p_avg = 0.5 * (p_ll + p_rr) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + + f_rho_sum = zero(T_rr) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * + (v_dot_n_ll + v_dot_n_rr) + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] + f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Note that this is the same as computing v_ll and v_rr with a normalized normal vector + # and then multiplying v by `norm_` again, but this version is slightly faster. + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + + # We treat the energy term analogous to the potential temperature term in the paper + # by Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3 = u_ll * v + f3 = f3 + p_ll * v + # density fluxes + f_rho = SVector{ncomponents(equations), real(equations)}(u_ll[ncomponents(equations) + i] * + v + for i in eachcomponent(equations)) + else + f1, f2, f3 = u_rr * v + f3 = f3 + p_rr * v + # density fluxes + f_rho = SVector{ncomponents(equations), real(equations)}(u_rr[ncomponents(equations) + i] * + v + for i in eachcomponent(equations)) + end + f1 = f1 + p * normal_direction[1] + f2 = f2 + p * normal_direction[2] + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end +end # @muladd diff --git a/src/equations/compressible_euler_multicomponent_moist_2d.jl b/src/equations/compressible_euler_multicomponent_moist_2d.jl new file mode 100644 index 00000000000..3406ea1a0e6 --- /dev/null +++ b/src/equations/compressible_euler_multicomponent_moist_2d.jl @@ -0,0 +1,390 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" +CompressibleMoistEulerEquations2D(;gammas, gas_constants, c_p, c_v, L_00) + +Similar to `CompressibleEulerMulticomponentEquations2D` but including latent heat of +vaporization +```math +... +``` + +Either the specific heat ratios `gammas` and the gas constants `gas_constants` in should be +passed as tuples, or the specific heats at constant volume `cv` and the specific heats at +constant pressure `cp`. The other quantities are then calculated considering a calorically +perfect gas. +""" +struct CompressibleMoistEulerEquations2D{NVARS, NCOMP, RealT <: Real} <: + AbstractCompressibleEulerMulticomponentEquations{2, NVARS, NCOMP} + gammas::SVector{NCOMP, RealT} + gas_constants::SVector{NCOMP, RealT} + c_p::SVector{NCOMP, RealT} + c_v::SVector{NCOMP, RealT} + L_00::RealT # latent heat of evaporation at 0 K + + function CompressibleMoistEulerEquations2D{NVARS, NCOMP, RealT}(gammas::SVector{NCOMP, + RealT}, + gas_constants::SVector{NCOMP, + RealT}, + c_p::SVector{NCOMP, + RealT}, + c_v::SVector{NCOMP, + RealT}, + L_00::RealT) where { + NVARS, + NCOMP, + RealT <: + Real + } + new(gammas, gas_constants, c_p, c_v, L_00) + end +end + +function CompressibleMoistEulerEquations2D(; + gammas = nothing, gas_constants = nothing, + c_p = nothing, c_v = nothing, L_00) + if gammas !== nothing && gas_constants !== nothing + _gammas = promote(gammas...) + _gas_constants = promote(gas_constants...) + RealT = promote_type(eltype(_gammas), eltype(_gas_constants), + typeof(gas_constants[1] / (gammas[1] - 1))) + _c_v = _gas_constants ./ (_gammas .- 1) + _c_p = _gas_constants + _gas_constants ./ (_gammas .- 1) + elseif c_p !== nothing && c_v !== nothing + _c_p = promote(c_p...) + _c_v = promote(c_v...) + RealT = promote_type(eltype(_c_p), eltype(_c_v), typeof(c_p[1] / c_v[1])) + _gas_constants = _c_p .- _c_v + _gammas = _c_p ./ _c_v + else + throw(DimensionMismatch("Either `gammas` and `gas_constants` or `c_p` and `c_v` \ + have to be filled with at least one value")) + end + + NVARS = length(_gammas) + 3 + NCOMP = length(_gammas) + + __gammas = SVector(map(RealT, _gammas)) + __gas_constants = SVector(map(RealT, _gas_constants)) + __c_p = SVector(map(RealT, _c_p)) + __c_v = SVector(map(RealT, _c_v)) + + return CompressibleMoistEulerEquations2D{NVARS, NCOMP, RealT}(__gammas, + __gas_constants, + __c_p, __c_v, L_00) +end + +@inline function Base.real(::CompressibleMoistEulerEquations2D{NVARS, NCOMP, + RealT}) where {NVARS, + NCOMP, + RealT} + RealT +end + +""" + totalgamma(u, equations::CompressibleMoistEulerEquations2D) + +Function that calculates the total gamma out of all partial gammas using the +partial density fractions as well as the partial specific heats at constant volume. +""" +@inline function totalgamma(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_v, c_p = equations + + help1 = zero(u[1]) + help2 = zero(u[1]) + + # compute weighted averages of cp and cv + # normalization by total rho not required, would cancel below + for i in eachcomponent(equations) + help1 += u[i + 3] * c_p[i] + help2 += u[i + 3] * c_v[i] + end + + return help1 / help2 +end + +""" +pressure(u, equations::CompressibleMoistEulerEquations2D) + +Calculate pressure. This differs from the calculation in `AbstractCompressibleEulerMulticomponentEquations` in that the latent heat is accounted for. +""" +@inline function pressure(u, equations::CompressibleMoistEulerEquations2D) + @unpack L_00 = equations + rho_v1, rho_v2, rho_e, rho_d, rho_v = u + rho = density(u, equations) + v1 = rho_v1 / rho + v2 = rho_v2 / rho + gamma = totalgamma(u, equations) + p = (gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - L_00 * rho_v) + return p +end + +""" +temperature(u, equations::CompressibleMoistEulerEquations2D) + +Calculate temperature. Account for latent heat. +""" +@inline function temperature(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_v, gammas, gas_constants, L_00 = equations + rho_v1, rho_v2, rho_e, rho_d, rho_v = u + + rho = density(u, equations) + help1 = zero(rho) + + # compute weighted average of cv + # normalization by rho not required, cancels below + for i in eachcomponent(equations) + help1 += u[i + 3] * c_v[i] + end + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + T = (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2) - L_00 * rho_v) / help1 + return T +end + +""" + density_gas_constant(u, equations::CompressibleMoistEulerEquations2D) + +Function that calculates overall density times overall gas constant. +""" +@inline function density_gas_constant(u, + equations::CompressibleMoistEulerEquations2D) + @unpack gas_constants = equations + help = zero(u[1]) + for i in eachcomponent(equations) + help += u[i + 3] * gas_constant[i] + end + return help +end + + +""" +cons2entropy(u, equations::CompressibleMoistEulerEquations2D) + +Convert conservative variables to entropy. +""" +@inline function cons2entropy(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_v, gammas, gas_constants = equations + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + # Multicomponent stuff + help1 = zero(rho) + gas_constant = zero(rho) + for i in eachcomponent(equations) + help1 += u[i + 3] * c_v[i] + gas_constant += gas_constants[i] * (u[i + 3] / rho) + end + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + p = pressure(u, equations) + rho_p = rho / p + T = temperature(u, equations) + + entrop_rho = SVector{ncomponents(equations), real(equations)}((c_v[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 3])) - + v_square / (2 * T)) + for i in eachcomponent(equations)) + + w1 = gas_constant * v1 * rho_p + w2 = gas_constant * v2 * rho_p + w3 = gas_constant * (-rho_p) + + entrop_other = SVector{3, real(equations)}(w1, w2, w3) + + return vcat(entrop_other, entrop_rho) +end + +""" +entropy2cons(w, equations::CompressibleMoistEulerEquations2D) + +Convert entropy variables to conservative variables +""" +@inline function entropy2cons(w, equations::CompressibleMoistEulerEquations2D) + @unpack gammas, gas_constants, cp, cv = equations + T = -1 / w[3] + v1 = w[1] * T + v2 = w[2] * T + v_squared = v1^2 + v2^2 + cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - + cv[i] * + (1 - log(T)) + + v_squared / + (2 * T)) / + gas_constants[i] - + 1) + for i in eachcomponent(equations)) + + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + u2 = rho * v2 + gamma = help1 / help2 + u3 = p / (gamma - 1) + 0.5 * rho * v_squared + L_00 * cons_rho[2] + cons_other = SVector{3, real(equations)}(u1, u2, u3) + return vcat(cons_other, cons_rho) +end + +""" + flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleMoistEulerEquations2D) + +Adaption of the entropy conserving two-point flux by +- Ayoub Gouasmi, Karthik Duraisamy (2020) + "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" + [arXiv:1904.00972v3](https://arxiv.org/abs/1904.00972) [math.NA] 4 Feb 2020 +""" +@inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, + equations::CompressibleMoistEulerEquations2D) + # Unpack left and right state + @unpack gammas, gas_constants, cv, L_00 = equations + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v_sum = v1_avg + v2_avg + + enth = zero(v_sum) + for i in eachcomponent(equations) + enth += rhok_avg[i] * gas_constants[i] + end + + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + + # Calculate fluxes depending on orientation + help1 = zero(T_ll) + help2 = zero(T_rr) + if orientation == 1 + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + help1 += f_rho[i] * cv[i] + help2 += f_rho[i] + end + f1 = (help2) * v1_avg + enth / T + f2 = (help2) * v2_avg + # Account for latent heat. This is the only difference compared to + # AbstractCompressibleEulerMulticomponentEquations + f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + v2_avg * f2 + L_00 * f_rho[2] + else + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + help1 += f_rho[i] * cv[i] + help2 += f_rho[i] + end + f1 = (help2) * v1_avg + f2 = (help2) * v2_avg + enth / T + # Account for latent heat. This is the only difference compared to + # AbstractCompressibleEulerMulticomponentEquations + f3 = (help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + + v2_avg * f2 + L_00 * f_rho[2] + end + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + # Unpack left and right state + @unpack gammas, gas_constants, c_v, L_00 = equations + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg + v1_square = 0.5 * (v1_ll^2 + v1_rr^2) + v2_square = 0.5 * (v2_ll^2 + v2_rr^2) + v_sum = v1_avg + v2_avg + + enth = zero(v_sum) + for i in eachcomponent(equations) + enth += rhok_avg[i] * gas_constants[i] + end + + T_ll = temperature(u_ll, equations) + T_rr = temperature(u_rr, equations) + T = 0.5 * (1.0 / T_ll + 1.0 / T_rr) + T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr) + + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v_dot_n_avg + for i in eachcomponent(equations)) + + help1 = zero(T_ll) + help2 = zero(T_rr) + for i in eachcomponent(equations) + help1 += f_rho[i] * c_v[i] + help2 += f_rho[i] + end + + f1 = (help2) * v1_avg + normal_direction[1] * enth / T + f2 = (help2) * v2_avg + normal_direction[2] * enth / T + + # Account for latent heat. This is the only difference compared to + # AbstractCompressibleEulerMulticomponentEquations + f3 = ((help1) / T_log - 0.5 * (v1_square + v2_square) * (help2) + + v1_avg * f1 + v2_avg * f2 + L_00 * f_rho[2]) + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 7a3c326984d..ec73e8ef2c1 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -414,24 +414,15 @@ include("compressible_euler_quasi_1d.jl") # CompressibleEulerMulticomponentEquations abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <: AbstractEquations{NDIMS, NVARS} end -include("compressible_euler_multicomponent_1d.jl") -include("compressible_euler_multicomponent_2d.jl") - -# PolytropicEulerEquations -abstract type AbstractPolytropicEulerEquations{NDIMS, NVARS} <: - AbstractEquations{NDIMS, NVARS} end -include("polytropic_euler_2d.jl") # Retrieve number of components from equation instance for the multicomponent case @inline function ncomponents(::AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, - NCOMP}) where { - NDIMS, - NVARS, - NCOMP - } + NCOMP}) where + {NDIMS, NVARS, NCOMP} NCOMP end + """ eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations) @@ -443,6 +434,16 @@ In particular, not the components themselves are returned. Base.OneTo(ncomponents(equations)) end +include("compressible_euler_multicomponent_1d.jl") +include("compressible_euler_multicomponent_abstract_2d.jl") +include("compressible_euler_multicomponent_2d.jl") +include("compressible_euler_multicomponent_moist_2d.jl") + +# PolytropicEulerEquations +abstract type AbstractPolytropicEulerEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end +include("polytropic_euler_2d.jl") + # Ideal MHD abstract type AbstractIdealGlmMhdEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end @@ -458,11 +459,8 @@ include("ideal_glm_mhd_multicomponent_2d.jl") # Retrieve number of components from equation instance for the multicomponent case @inline function ncomponents(::AbstractIdealGlmMhdMulticomponentEquations{NDIMS, NVARS, - NCOMP}) where { - NDIMS, - NVARS, - NCOMP - } + NCOMP}) where + {NDIMS, NVARS, NCOMP} NCOMP end """ diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index e5d45ebcc07..d969cbfc055 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -273,6 +273,32 @@ end end end +@trixi_testset "elixir_eulermultimoist_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermultimoist_warm_bubble.jl"), + l2=[ + 1.5123651627525257e-5, + 1.51236516273878e-5, + 2.4544918394022538e-5, + 5.904791661362391e-6, + 1.1809583322724782e-5, + ], + linf=[ + 8.393471747591974e-5, + 8.393471748258108e-5, + 0.00015028562494778797, + 3.504466610437795e-5, + 7.00893322087559e-5, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), # Expected errors are exactly the same as with TreeMesh! From 405bf74513e4684c7c784b9d2297c8aa12b5534d Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jan 2024 11:58:17 +0100 Subject: [PATCH 32/41] fix spelling --- .../elixir_eulermultimoist_warm_bubble.jl | 16 ++++++++-------- ...compressible_euler_multicomponent_moist_2d.jl | 1 - test/test_structured_2d.jl | 4 +++- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl index c0d8966af5e..3fef2e044ea 100644 --- a/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl @@ -92,9 +92,9 @@ end # Initial condition struct AtmossphereLayers{RealT <: Real} setup::WarmBubbleSetup - # structure: 1--> i-layer (z = total_hight/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql + # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql LayerData::Matrix{RealT} - total_hight::RealT + total_height::RealT preciseness::Int layers::Int ground_state::NTuple{2, RealT} @@ -102,7 +102,7 @@ struct AtmossphereLayers{RealT <: Real} mixing_ratios::NTuple{2, RealT} end -function AtmossphereLayers(setup; total_hight = 10010.0, preciseness = 10, +function AtmossphereLayers(setup; total_height = 10010.0, preciseness = 10, ground_state = (1.4, 100000.0), equivalentpotential_temperature = 320, mixing_ratios = (0.02, 0.02), RealT = Float64) @@ -115,7 +115,7 @@ function AtmossphereLayers(setup; total_hight = 10010.0, preciseness = 10, T0 = theta_e0 y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] - n = convert(Int, total_hight / preciseness) + n = convert(Int, total_height / preciseness) dz = 0.01 LayerData = zeros(RealT, n + 1, 4) @@ -145,7 +145,7 @@ function AtmossphereLayers(setup; total_hight = 10010.0, preciseness = 10, LayerData[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] end - return AtmossphereLayers{RealT}(setup, LayerData, total_hight, dz, n, ground_state, + return AtmossphereLayers{RealT}(setup, LayerData, total_height, dz, n, ground_state, theta_e0, mixing_ratios) end @@ -190,17 +190,17 @@ end # https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml. function initial_condition_moist_bubble(x, t, setup::WarmBubbleSetup, AtmosphereLayers::AtmossphereLayers) - @unpack LayerData, preciseness, total_hight = AtmosphereLayers + @unpack LayerData, preciseness, total_height = AtmosphereLayers dz = preciseness z = x[2] - if (z > total_hight && !(isapprox(z, total_hight))) + if (z > total_height && !(isapprox(z, total_height))) error("The atmossphere does not match the simulation domain") end n = convert(Int, floor((z + eps()) / dz)) + 1 z_l = (n - 1) * dz (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :] z_r = n * dz - if (z_l == total_hight) + if (z_l == total_height) z_r = z_l + dz n = n - 1 end diff --git a/src/equations/compressible_euler_multicomponent_moist_2d.jl b/src/equations/compressible_euler_multicomponent_moist_2d.jl index 3406ea1a0e6..c7bdf98bb1c 100644 --- a/src/equations/compressible_euler_multicomponent_moist_2d.jl +++ b/src/equations/compressible_euler_multicomponent_moist_2d.jl @@ -164,7 +164,6 @@ Function that calculates overall density times overall gas constant. return help end - """ cons2entropy(u, equations::CompressibleMoistEulerEquations2D) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index d969cbfc055..73f236df946 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -288,7 +288,9 @@ end 0.00015028562494778797, 3.504466610437795e-5, 7.00893322087559e-5, - ]) + ], + cells_per_dimension=(32, 16), + tspan=(0.0, 10.0)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From c52da4e973d9638d20df2038d3d8a803d6b847c5 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jan 2024 13:48:09 +0100 Subject: [PATCH 33/41] fix density_gas_constant calls --- src/equations/compressible_euler_multicomponent_2d.jl | 2 +- .../compressible_euler_multicomponent_abstract_2d.jl | 8 ++++---- .../compressible_euler_multicomponent_moist_2d.jl | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 1e7a6459600..11290349269 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -263,7 +263,7 @@ Function that calculates overall density times overall gas constant. end """ - flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D) + flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerMulticomponentEquations2D) Adaption of the entropy conserving two-point flux by - Ayoub Gouasmi, Karthik Duraisamy (2020) diff --git a/src/equations/compressible_euler_multicomponent_abstract_2d.jl b/src/equations/compressible_euler_multicomponent_abstract_2d.jl index 889b380b582..6dc923fea78 100644 --- a/src/equations/compressible_euler_multicomponent_abstract_2d.jl +++ b/src/equations/compressible_euler_multicomponent_abstract_2d.jl @@ -493,8 +493,8 @@ the Euler Equations Using Summation-by-Parts Operators velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) # helpful variables - enth_ll = density_gas_constant(u_ll) - enth_rr = density_gas_constant(u_rr) + enth_ll = density_gas_constant(u_ll, equations) + enth_rr = density_gas_constant(u_rr, equations) # temperature and pressure T_ll = temperature(u_ll, equations) @@ -562,8 +562,8 @@ end v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # helpful variables - enth_ll = density_gas_constant(u_ll) - enth_rr = density_gas_constant(u_rr) + enth_ll = density_gas_constant(u_ll, equations) + enth_rr = density_gas_constant(u_rr, equations) # temperature and pressure T_ll = temperature(u_ll, equations) diff --git a/src/equations/compressible_euler_multicomponent_moist_2d.jl b/src/equations/compressible_euler_multicomponent_moist_2d.jl index c7bdf98bb1c..9d1f00898d4 100644 --- a/src/equations/compressible_euler_multicomponent_moist_2d.jl +++ b/src/equations/compressible_euler_multicomponent_moist_2d.jl @@ -245,7 +245,7 @@ Convert entropy variables to conservative variables end """ - flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleMoistEulerEquations2D) + flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleMoistEulerEquations2D) Adaption of the entropy conserving two-point flux by - Ayoub Gouasmi, Karthik Duraisamy (2020) From d9efc2aa8fa472a94ca7943a6b824a276d9e1bb5 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jan 2024 14:49:20 +0100 Subject: [PATCH 34/41] some cleanup --- src/equations/compressible_euler_multicomponent_2d.jl | 6 ++---- .../compressible_euler_multicomponent_moist_2d.jl | 10 +++++----- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 11290349269..28dec1fa711 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -108,9 +108,7 @@ end Calculate temperature. """ @inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) - @unpack cv, gammas, gas_constants = equations - - rho_v1, rho_v2, rho_e, rho_d, rho_v = u + rho_v1, rho_v2, rho_e = u rho = density(u, equations) help1 = zero(rho) @@ -257,7 +255,7 @@ Function that calculates overall density times overall gas constant. @unpack gas_constants = equations help = zero(u[1]) for i in eachcomponent(equations) - help += u[i + 3] * gas_constant[i] + help += u[i + 3] * gas_constants[i] end return help end diff --git a/src/equations/compressible_euler_multicomponent_moist_2d.jl b/src/equations/compressible_euler_multicomponent_moist_2d.jl index 9d1f00898d4..6e64f50444b 100644 --- a/src/equations/compressible_euler_multicomponent_moist_2d.jl +++ b/src/equations/compressible_euler_multicomponent_moist_2d.jl @@ -130,8 +130,8 @@ temperature(u, equations::CompressibleMoistEulerEquations2D) Calculate temperature. Account for latent heat. """ @inline function temperature(u, equations::CompressibleMoistEulerEquations2D) - @unpack c_v, gammas, gas_constants, L_00 = equations - rho_v1, rho_v2, rho_e, rho_d, rho_v = u + @unpack c_v, L_00 = equations + rho_v1, rho_v2, rho_e = u rho = density(u, equations) help1 = zero(rho) @@ -159,7 +159,7 @@ Function that calculates overall density times overall gas constant. @unpack gas_constants = equations help = zero(u[1]) for i in eachcomponent(equations) - help += u[i + 3] * gas_constant[i] + help += u[i + 3] * gas_constants[i] end return help end @@ -170,7 +170,7 @@ cons2entropy(u, equations::CompressibleMoistEulerEquations2D) Convert conservative variables to entropy. """ @inline function cons2entropy(u, equations::CompressibleMoistEulerEquations2D) - @unpack c_v, gammas, gas_constants = equations + @unpack c_v, gas_constants = equations rho_v1, rho_v2, rho_e = u rho = density(u, equations) @@ -212,7 +212,7 @@ entropy2cons(w, equations::CompressibleMoistEulerEquations2D) Convert entropy variables to conservative variables """ @inline function entropy2cons(w, equations::CompressibleMoistEulerEquations2D) - @unpack gammas, gas_constants, cp, cv = equations + @unpack gammas, gas_constants, cv = equations T = -1 / w[3] v1 = w[1] * T v2 = w[2] * T From f95024de8c38eeaf6e0d312c81a7ff6ec81bcd12 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jan 2024 17:50:19 +0100 Subject: [PATCH 35/41] missing cv --- src/equations/compressible_euler_multicomponent_2d.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 28dec1fa711..d02ceccf56b 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -109,6 +109,7 @@ Calculate temperature. """ @inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) rho_v1, rho_v2, rho_e = u + @unpack cv = equations rho = density(u, equations) help1 = zero(rho) From f32fdff9582a046b5ab77b494edccaabe4ac9e7e Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jan 2024 17:50:48 +0100 Subject: [PATCH 36/41] add shima, kenndy_gruber, and LMARS fluxes with tests --- ...ssible_euler_multicomponent_abstract_2d.jl | 261 +++++++++++++++++- test/test_structured_2d.jl | 2 + test/test_unit.jl | 34 ++- 3 files changed, 269 insertions(+), 28 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_abstract_2d.jl b/src/equations/compressible_euler_multicomponent_abstract_2d.jl index 6dc923fea78..02915eb56c6 100644 --- a/src/equations/compressible_euler_multicomponent_abstract_2d.jl +++ b/src/equations/compressible_euler_multicomponent_abstract_2d.jl @@ -506,8 +506,7 @@ the Euler Equations Using Summation-by-Parts Operators f_rho_sum = zero(T_rr) if orientation == 1 - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg - for i in eachcomponent(equations)) + f_rho = rhok_mean .* v1_avg for i in eachcomponent(equations) f_rho_sum += f_rho[i] end @@ -516,8 +515,7 @@ the Euler Equations Using Summation-by-Parts Operators f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + 0.5 * (p_ll * v1_rr + p_rr * v1_ll) else - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v2_avg - for i in eachcomponent(equations)) + f_rho = rhok_mean .* v2_avg for i in eachcomponent(equations) f_rho_sum += f_rho[i] end @@ -574,9 +572,7 @@ end inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) f_rho_sum = zero(T_rr) - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * - (v_dot_n_ll + v_dot_n_rr) - for i in eachcomponent(equations)) + f_rho = rhok_mean .* (0.5 * (v_dot_n_ll + v_dot_n_rr)) for i in eachcomponent(equations) f_rho_sum += f_rho[i] end @@ -591,6 +587,249 @@ end return vcat(f_other, f_rho) end +""" + flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +This flux is is a modification of the original kinetic energy preserving two-point flux by +- Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) + Kinetic energy and entropy preserving schemes for compressible flows + by split convective forms + [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058) + +The modification is in the energy flux to guarantee pressure equilibrium and was developed by +- Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) + Preventing spurious pressure oscillations in split convective form discretizations for + compressible flows + [DOI: 10.1016/j.jcp.2020.110060](https://doi.org/10.1016/j.jcp.2020.110060) +""" +@inline function flux_shima_etal(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # Average each factor of products in flux + f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + kin_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on orientation + if orientation == 1 + f_rho = f_rho .* v1_avg + rho_avg_v = rho_avg * v1_avg + pv1_avg = 0.5 * (p_ll * v1_rr + p_rr * v1_ll) + f1 = rho_avg_v * v1_avg + p_avg + f2 = rho_avg_v * v2_avg + f3 = p_avg * v1_avg * inv_gamma_minus_one + rho_avg_v * kin_avg + pv1_avg + else + f_rho = f_rho .* v2_avg + rho_avg_v = rho_avg * v2_avg + pv2_avg = 0.5 * (p_ll * v2_rr + p_rr * v2_ll) + f1 = rho_avg_v * v1_avg + f2 = rho_avg_v * v2_avg + p_avg + f3 = p_avg * v2_avg * inv_gamma_minus_one + rho_avg_v * kin_avg + pv2_avg + end + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # Average each factor of products in flux + f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr) + p_avg = 0.5 * (p_ll + p_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + + # Calculate fluxes depending on normal_direction + f_rho = f_rho .* v_dot_n_avg + rho_avg_v = rho_avg * v_dot_n_avg + f1 = rho_avg_v * v1_avg + p_avg * normal_direction[1] + f2 = rho_avg_v * v2_avg + p_avg * normal_direction[2] + f3 = (rho_avg_v * velocity_square_avg + + p_avg * v_dot_n_avg * inv_gamma_minus_one + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll)) + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +""" + flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Kinetic energy preserving two-point flux by +- Kennedy and Gruber (2008) + Reduced aliasing formulations of the convective terms within the + Navier-Stokes equations for a compressible fluid + [DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020) +""" +@inline function flux_kennedy_gruber(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_e_ll = u_ll[3] + rho_e_rr = u_rr[3] + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Average each factor of products in flux + f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on orientation + if orientation == 1 + f_rho = f_rho .* v1_avg + f1 = rho_avg * v1_avg * v1_avg + p_avg + f2 = rho_avg * v1_avg * v2_avg + f3 = (rho_avg * e_avg + p_avg) * v1_avg + else + f_rho = f_rho .* v2_avg + f1 = rho_avg * v2_avg * v1_avg + f2 = rho_avg * v2_avg * v2_avg + p_avg + f3 = (rho_avg * e_avg + p_avg) * v2_avg + end + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +@inline function flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + # Unpack left and right state + rho_e_ll = u_ll[3] + rho_e_rr = u_rr[3] + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Average each factor of products in flux + f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + rho_avg = 0.5 * (rho_ll + rho_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2] + p_avg = 0.5 * (p_ll + p_rr) + e_avg = 0.5 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr) + + # Calculate fluxes depending on normal_direction + f_rho = f_rho .* v_dot_n_avg + rho_avg_v = rho_avg * v_dot_n_avg + f1 = rho_avg_v * v1_avg + p_avg * normal_direction[1] + f2 = rho_avg_v * v2_avg + p_avg * normal_direction[2] + f3 = rho_avg_v * e_avg + p_avg * v_dot_n_avg + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using +an estimate `c` of the speed of sound. + +References: +- Xi Chen et al. (2013) + A Control-Volume Model of the Compressible Euler Equations with a Vertical + Lagrangian Coordinate + [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) +""" + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3 = v * u_ll + f3 = f3 + p_ll * v + # density fluxes + f_rho = u_ll[4:end] .* v + else + f1, f2, f3 = v * u_rr + f3 = f3 + p_rr * v + # density fluxes + f_rho = u_rr[4:end] .* v + end + + if orientation == 1 + f1 = f1 + p + else # orientation == 2 + f2 = f2 + p + end + + # momentum and energy flux + f_other = SVector(f1, f2, f3) + + return vcat(f_other, f_rho) +end + @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, equations::AbstractCompressibleEulerMulticomponentEquations{2}) c = flux_lmars.speed_of_sound @@ -619,16 +858,12 @@ end f1, f2, f3 = u_ll * v f3 = f3 + p_ll * v # density fluxes - f_rho = SVector{ncomponents(equations), real(equations)}(u_ll[ncomponents(equations) + i] * - v - for i in eachcomponent(equations)) + f_rho = u_ll[4:end] .* v else f1, f2, f3 = u_rr * v f3 = f3 + p_rr * v # density fluxes - f_rho = SVector{ncomponents(equations), real(equations)}(u_rr[ncomponents(equations) + i] * - v - for i in eachcomponent(equations)) + f_rho = u_rr[4:end] .* v end f1 = f1 + p * normal_direction[1] f2 = f2 + p * normal_direction[2] diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 73f236df946..e4d9717a8fa 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -273,6 +273,7 @@ end end end +#= @trixi_testset "elixir_eulermultimoist_warm_bubble.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermultimoist_warm_bubble.jl"), l2=[ @@ -300,6 +301,7 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end +=# @trixi_testset "elixir_euler_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), diff --git a/test/test_unit.jl b/test/test_unit.jl index e8a8effbe29..31cfa73254f 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -614,14 +614,17 @@ end end @timed_testset "Consistency check for single point flux: CEMCE" begin - equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), - gas_constants = (0.4, 0.4)) + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.2), + gas_constants = (0.4, 0.6)) u = SVector(0.1, -0.5, 1.0, 1.0, 2.0) orientations = [1, 2] - for orientation in orientations - @test flux(u, orientation, equations) ≈ - flux_ranocha(u, u, orientation, equations) + fluxes = [flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340)] + for f_test in fluxes + for orientation in orientations + @test flux(u, orientation, equations) ≈ + f_test(u, u, orientation, equations) + end end end @@ -1332,25 +1335,26 @@ end @testset "FluxRotated vs. direct implementation" begin @timed_testset "CompressibleEulerMulticomponentEquations2D" begin - equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), + equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.2), gas_constants = (0.4, - 0.4)) + 0.6)) normal_directions = [SVector(1.0, 0.0), SVector(0.0, 1.0), SVector(0.5, -0.5), SVector(-1.2, 0.3)] u_values = [SVector(0.1, -0.5, 1.0, 1.0, 2.0), SVector(-0.1, -0.3, 1.2, 1.3, 1.4)] + fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, + FluxLMARS(340)] - f_std = flux - f_rot = FluxRotated(f_std) - println(typeof(f_std)) - println(typeof(f_rot)) - for u in u_values, - normal_direction in normal_directions + for f_std in fluxes + f_rot = FluxRotated(f_std) + for u_ll in u_values, u_rr in u_values, + normal_direction in normal_directions - @test f_rot(u, normal_direction, equations) ≈ - f_std(u, normal_direction, equations) + @test f_rot(u_ll, u_rr, normal_direction, equations) ≈ + f_std(u_ll, u_rr, normal_direction, equations) + end end end From 778df22d9848b163225e75f4489af11ea84dafba Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jan 2024 17:52:12 +0100 Subject: [PATCH 37/41] format --- test/test_unit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 31cfa73254f..dd7235a3957 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1345,7 +1345,7 @@ end u_values = [SVector(0.1, -0.5, 1.0, 1.0, 2.0), SVector(-0.1, -0.3, 1.2, 1.3, 1.4)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - FluxLMARS(340)] + FluxLMARS(340)] for f_std in fluxes f_rot = FluxRotated(f_std) From 01abb39e93a87eec99d563a1625b3e3d6d42ce01 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 26 Jan 2024 16:59:51 +0100 Subject: [PATCH 38/41] fight allocations --- ...ssible_euler_multicomponent_abstract_2d.jl | 50 +++++++++++++++---- ...pressible_euler_multicomponent_moist_2d.jl | 14 +++++- 2 files changed, 54 insertions(+), 10 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_abstract_2d.jl b/src/equations/compressible_euler_multicomponent_abstract_2d.jl index 02915eb56c6..5bf2acb415c 100644 --- a/src/equations/compressible_euler_multicomponent_abstract_2d.jl +++ b/src/equations/compressible_euler_multicomponent_abstract_2d.jl @@ -265,7 +265,7 @@ function temperature end """ pressure(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) -Function that calculates the overall pressure. +Calculates overall pressure. """ @inline function pressure(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) @@ -279,6 +279,22 @@ Function that calculates the overall pressure. return p end +""" + density_pressure(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) + +Calculates overall density times overall pressure. +""" +@inline function density_pressure(u, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + gamma = totalgamma(u, equations) + + rho_times_p = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) + return rho_times_p +end + """ density_gas_constant(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) @@ -617,7 +633,9 @@ The modification is in the energy flux to guarantee pressure equilibrium and was inv_gamma_minus_one = 1 / (gamma - 1) # Average each factor of products in flux - f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) rho_avg = 0.5 * (rho_ll + rho_rr) v1_avg = 0.5 * (v1_ll + v1_rr) v2_avg = 0.5 * (v2_ll + v2_rr) @@ -663,7 +681,9 @@ end inv_gamma_minus_one = 1 / (gamma - 1) # Average each factor of products in flux - f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) rho_avg = 0.5 * (rho_ll + rho_rr) v1_avg = 0.5 * (v1_ll + v1_rr) v2_avg = 0.5 * (v2_ll + v2_rr) @@ -708,7 +728,9 @@ Kinetic energy preserving two-point flux by rho_rr = density(u_rr, equations) # Average each factor of products in flux - f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) rho_avg = 0.5 * (rho_ll + rho_rr) v1_avg = 0.5 * (v1_ll + v1_rr) v2_avg = 0.5 * (v2_ll + v2_rr) @@ -746,7 +768,9 @@ end rho_rr = density(u_rr, equations) # Average each factor of products in flux - f_rho = 0.5 * (u_ll[4:end] + u_rr[4:end]) + f_rho = SVector{ncomponents(equations), + real(equations)}(0.5 * (u_ll[i + 3] + u_rr[i + 3]) + for i in eachcomponent(equations)) rho_avg = 0.5 * (rho_ll + rho_rr) v1_avg = 0.5 * (v1_ll + v1_rr) v2_avg = 0.5 * (v2_ll + v2_rr) @@ -810,12 +834,16 @@ References: f1, f2, f3 = v * u_ll f3 = f3 + p_ll * v # density fluxes - f_rho = u_ll[4:end] .* v + f_rho = SVector{ncomponents(equations), + real(equations)}(u_ll[i + 3] * v + for i in eachcomponent(equations)) else f1, f2, f3 = v * u_rr f3 = f3 + p_rr * v # density fluxes - f_rho = u_rr[4:end] .* v + f_rho = SVector{ncomponents(equations), + real(equations)}(u_rr[i + 3] * v + for i in eachcomponent(equations)) end if orientation == 1 @@ -858,12 +886,16 @@ end f1, f2, f3 = u_ll * v f3 = f3 + p_ll * v # density fluxes - f_rho = u_ll[4:end] .* v + f_rho = SVector{ncomponents(equations), + real(equations)}(u_ll[i + 3] * v + for i in eachcomponent(equations)) else f1, f2, f3 = u_rr * v f3 = f3 + p_rr * v # density fluxes - f_rho = u_rr[4:end] .* v + f_rho = SVector{ncomponents(equations), + real(equations)}(u_rr[i + 3] * v + for i in eachcomponent(equations)) end f1 = f1 + p * normal_direction[1] f2 = f2 + p * normal_direction[2] diff --git a/src/equations/compressible_euler_multicomponent_moist_2d.jl b/src/equations/compressible_euler_multicomponent_moist_2d.jl index 6e64f50444b..1e3e2b28d9f 100644 --- a/src/equations/compressible_euler_multicomponent_moist_2d.jl +++ b/src/equations/compressible_euler_multicomponent_moist_2d.jl @@ -124,6 +124,18 @@ Calculate pressure. This differs from the calculation in `AbstractCompressibleEu return p end +""" + density_pressure(u, equations::CompressibleMoistEulerEquations2D) + +Calculates overall density times overall pressure. +""" +@inline function density_pressure(u, + equations::CompressibleMoistEulerEquations2D) + rho = density(u, equations) + p = pressure(u, equations) + return rho * p +end + """ temperature(u, equations::CompressibleMoistEulerEquations2D) @@ -131,7 +143,7 @@ Calculate temperature. Account for latent heat. """ @inline function temperature(u, equations::CompressibleMoistEulerEquations2D) @unpack c_v, L_00 = equations - rho_v1, rho_v2, rho_e = u + rho_v1, rho_v2, rho_e, rho_d, rho_v = u rho = density(u, equations) help1 = zero(rho) From 393e69e52923b3ffe8efc8be0cec99d8bd300c30 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 26 Feb 2024 12:39:46 +0100 Subject: [PATCH 39/41] avoid allocation? --- .../compressible_euler_multicomponent_abstract_2d.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_abstract_2d.jl b/src/equations/compressible_euler_multicomponent_abstract_2d.jl index 5bf2acb415c..46971dd27aa 100644 --- a/src/equations/compressible_euler_multicomponent_abstract_2d.jl +++ b/src/equations/compressible_euler_multicomponent_abstract_2d.jl @@ -165,9 +165,11 @@ Should be used together with [`UnstructuredMesh2D`](@ref). (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) end - return SVector(p_star * normal[1], - p_star * normal[2], - zeros(eltype(u_inner), ncomponents(equations) + 1)...) * norm_ + f_v = SVector{2, real(equations)}(p_star * normal_direction[1], + p_star * normal_direction[2]) + f_other = zeros(SVector{ncomponents(equations) + 1, real(equations)}) + + return vcat(f_v, f_other) end """ From ce6b975bf73eb0ba6500c6687111581b4abf89e5 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 26 Feb 2024 12:40:09 +0100 Subject: [PATCH 40/41] implement max_abs_speed_naive --- ...ssible_euler_multicomponent_abstract_2d.jl | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/src/equations/compressible_euler_multicomponent_abstract_2d.jl b/src/equations/compressible_euler_multicomponent_abstract_2d.jl index 46971dd27aa..1b80367ea10 100644 --- a/src/equations/compressible_euler_multicomponent_abstract_2d.jl +++ b/src/equations/compressible_euler_multicomponent_abstract_2d.jl @@ -454,6 +454,36 @@ end λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) end +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::AbstractCompressibleEulerMulticomponentEquations{2}) + + rho_v1_ll, rho_v2_ll = u_ll + rho_v1_rr, rho_v2_rr = u_rr + + # Get the density and gas gamma + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + gamma_ll = totalgamma(u_ll, equations) + gamma_rr = totalgamma(u_rr, equations) + + # Calculate normal velocities and sound speed + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + v_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2]) + v_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2]) + + # Compute the sound speeds on the left and right + p_ll = pressure(u_ll, equations) + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + p_rr = pressure(u_rr, equations) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + @inline function max_abs_speeds(u, equations::AbstractCompressibleEulerMulticomponentEquations{2}) rho_v1, rho_v2 = u From cadf32b2f617acc239a53542ba0e5597e7c5cc0b Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 24 Jul 2024 09:51:28 +0200 Subject: [PATCH 41/41] alternative initial condition --- .../elixir_eulermultimoist_warm_bubble.jl | 313 +++++++++++++++--- 1 file changed, 271 insertions(+), 42 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl index 3fef2e044ea..80a42697d51 100644 --- a/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_eulermultimoist_warm_bubble.jl @@ -8,26 +8,144 @@ using Infiltrator # Bryan and Fritsch (2002) # A Benchmark Simulation for Moist Nonhydrostatic Numerical Models # [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) + +@inline function getqvs(eps, p, t) + # Equation for saturation vapor pressure + # Bolton (1980, MWR, p. 1047) + es = 611.2 * exp(17.67 * (t - 273.15) / (t - 29.65)) + return eps * es / (p - es) +end + +# Compute a base state in hydrostatic equilibrium +# The hydrostatic ODE is integrated upwards from surface (z=0) to `height` +# A one dimensional distribution (in z) with `nz` points is stored +# Adapted from https://www2.mmm.ucar.edu/people/bryan/Code/mbm.F +function calculate_hydrostatic_base_state(g, c_pd, R_d, c_pv, R_v, c_pl, L_00, + p_ref, theta_ref, r_t, nz, height) + eps = R_d / R_v + tolerance = 0.0001 # tolerance for iterative scheme + relaxation = 0.3 # relaxation factor for iterative scheme + + theta_0 = Vector{Float64}(undef, nz) # potential temperature + theta_rho_0 = Vector{Float64}(undef, nz) # density potential temperature + p_0 = Vector{Float64}(undef, nz) # pressure + q_v0 = Vector{Float64}(undef, nz) # vapor mixing ratio + q_l0 = Vector{Float64}(undef, nz) # liquid water mixing ratio + + # Start at the surface + p_0[1] = p_ref # ==> exner = 1 + q_v0[1] = getqvs(eps, p_ref, theta_ref) # T = theta_ref * exner = theta_ref + q_l0[1] = r_t - q_v0[1] + theta_0[1] = theta_ref + theta_rho_0[1] = theta_ref * (1.0 + inv(eps) * q_v0[1]) / (1.0 + q_v0[1] + q_l0[1]) + + T_last = theta_0[1] # T = theta_ref * exner = theta_ref + exner_last = 1 + + # Integrate upwards + dz = height / (nz - 1) + for k in 2:nz + theta_new = theta_0[k - 1] + theta_tmp = theta_0[k - 1] + theta_rho_tmp = theta_rho_0[k - 1] + + exner_tmp = 1.0 + p_tmp = 1.0 + T_tmp = 1.0 + q_v_tmp = 1.0 + q_l_tmp = 1.0 + steps = 0 + converged = false + while !converged + steps = steps + 1 + + # Trapezoidal rule? Only applied to theta_rho? + exner_tmp = exner_last - + dz * g / (c_pd * 0.5 * (theta_rho_0[k - 1] + theta_rho_tmp)) + p_tmp = p_ref * (exner_tmp^(c_pd / R_d)) + T_tmp = theta_new * exner_tmp + q_v_tmp = getqvs(eps, p_tmp, T_tmp) + q_l_tmp = r_t - q_v_tmp + theta_rho_tmp = theta_tmp * (1.0 + inv(eps) * q_v_tmp) / + (1.0 + q_v_tmp + q_l_tmp) + + # Trapezoidal rule? + T_bar = 0.5 * (T_last + T_tmp) + q_vbar = 0.5 * (q_v0[k - 1] + q_v_tmp) + q_lbar = 0.5 * (q_l0[k - 1] + q_l_tmp) + + lhv = L_00 - (c_pl - c_pv) * T_bar # latent heat + R_m = R_d + R_v * q_vbar # gas constant for moist air + c_pm = c_pd + c_pv * q_vbar + c_pl * q_lbar # c_p for moist air + + # An exact thermodynamic equation (e.g., for CM1) + theta_tmp = theta_0[k - 1] * + exp(-lhv * (q_v_tmp - q_v0[k - 1]) / (c_pm * T_bar) + + (R_m / c_pm - R_d / c_pd) * log(p_tmp / p_0[k - 1])) + + if abs(theta_tmp - theta_new) > tolerance + theta_new = theta_new + relaxation * (theta_tmp - theta_new) + if steps > 100 + error("calculate_hydrostatic_base_state: not converging!") + end + else + converged = true + if q_l_tmp < 0.0 + error("calculate_hydrostatic_base_state: liquid ratio negative!") + end + end + end + + theta_0[k] = theta_tmp + theta_rho_0[k] = theta_tmp * (1.0 + inv(eps) * q_v_tmp) / (1.0 + q_v_tmp + q_l_tmp) + p_0[k] = p_tmp + q_v0[k] = q_v_tmp + q_l0[k] = q_l_tmp + + T_last = T_tmp + exner_last = exner_tmp + end + + return theta_0, theta_rho_0, p_0, q_v0, q_l0 +end + struct WarmBubbleSetup - # Physical constants - g::Float64 # gravity of earth - c_pd::Float64 # heat capacity for constant pressure (dry air) - c_vd::Float64 # heat capacity for constant volume (dry air) - R_d::Float64 # gas constant (dry air) - c_pv::Float64 # heat capacity for constant pressure (water vapor) - c_vv::Float64 # heat capacity for constant volume (water vapor) - R_v::Float64 # gas constant (water vapor) - c_pl::Float64 # heat capacity for constant pressure (liquid water) - p_0::Float64 # reference pressure 1000 hPa - L_00::Float64 # latent heat of evaporation at 0 K - # TODO value ??? - # L00 = parameter("L00",2.5000e6 + (c_pl - c_pv) * 273.15) + g::Float64 # gravity of earth + c_pd::Float64 # heat capacity for constant pressure (dry air) + c_vd::Float64 # heat capacity for constant volume (dry air) + R_d::Float64 # gas constant (dry air) + c_pv::Float64 # heat capacity for constant pressure (water vapor) + c_vv::Float64 # heat capacity for constant volume (water vapor) + R_v::Float64 # gas constant (water vapor) + c_pl::Float64 # heat capacity for constant pressure (liquid water) + p_ref::Float64 # reference pressure + theta_ref::Float64 # potential temperature at surface + L_00::Float64 # latent heat of evaporation at 0 K + r_t::Float64 # total water mixing ratio + theta_0::Vector{Float64} # potential temperature + theta_rho_0::Vector{Float64} # density potential temperature + p_0::Vector{Float64} # pressure + q_v0::Vector{Float64} # vapor mixing ratio + q_l0::Vector{Float64} # liquid water mixing ratio + nz::Int64 # resolution for vertical hydrostatic distribution + height::Float64 # maximal height in vertical hydrostatic distribution function WarmBubbleSetup(; g = 9.81, c_pd = 1004.0, c_vd = 717.0, R_d = c_pd - c_vd, c_pv = 1885, c_vv = 1424.0, R_v = c_pv - c_vv, - c_pl = 4186.0, p_0 = 100_000.0, L_00 = 3147620.0) - new(g, c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_0, L_00) + c_pl = 4186.0, + p_ref = 100_000.0, theta_ref = 289.8486, + L_00 = 2.5e6 + (c_pl - c_pv) * 273.15, + r_t = 0.020, + nz, height) + theta_0, + theta_rho_0, + p_0, + q_v0, + q_l0 = calculate_hydrostatic_base_state(g, c_pd, R_d, c_pv, R_v, c_pl, L_00, p_ref, + theta_ref, r_t, nz, height) + new(g, c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_ref, theta_ref, L_00, r_t, theta_0, + theta_rho_0, p_0, q_v0, q_l0, nz, height) end end @@ -51,8 +169,27 @@ end zero(eltype(u)), zero(eltype(u)), zero(eltype(u))) end -@inline function source_terms_phase_change(u, setup::WarmBubbleSetup, - equations::CompressibleMoistEulerEquations2D) +# condensation, equation by Rutledge and Hobbs (1983) +@inline function source_terms_phase_change_rh(u, setup::WarmBubbleSetup, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, R_d, c_pv, R_v, c_pl, L_00 = setup + eps = R_d / R_v + rho_v1, rho_v2, rho_e, rho_qd, rho_qv, rho_ql = u + rho = density(u, equations) + T = temperature(u, equations) + p = pressure(u, equations) + + q_vs = getqvs(eps, p, T) + q_v = rho_qv / rho + lhv = L_00 - (c_pl - c_pv) * T + rcond = (q_v - q_vs) / (1.0 + q_vs * lhv^2 / (c_pd * R_v * T^2)) + #qcond(mgs) = Min( Max( 0.0, tmp ), (qvap(mgs)-qvs(mgs)) ) + return SVector(zero(eltype(u)), zero(eltype(u)), zero(eltype(u)), + zero(eltype(u)), -rcond, rcond) +end + +@inline function source_terms_phase_change_fb(u, setup::WarmBubbleSetup, + equations::CompressibleMoistEulerEquations2D) # This source term models the phase chance between could water and vapor. @unpack R_v = setup @@ -74,7 +211,7 @@ end # saturation control factor # < 1: stronger saturation effect # > 1: weaker saturation effect - C = 1 + C = 94.0 Q_ph = (a + b - sqrt(a^2 + b^2)) * C @@ -85,11 +222,99 @@ end @inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleMoistEulerEquations2D) return source_terms_geopotential(u, setup, equations) + - source_terms_phase_change(u, setup, equations) + source_terms_phase_change_fb(u, setup, equations) end ############################################################################### # Initial condition + +# Use the precomputed vertical distribution in `WarmBubbleSetup` to interpolate and perturb +# the initial solution +@inline function (setup::WarmBubbleSetup)(x, t, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_vd, R_d, c_vv, R_v, c_pl, L_00, r_t, p_ref, theta_0, + theta_rho_0, p_0, q_v0, q_l0, nz, height = setup + eps = R_d / R_v + tolerance = 0.0001 # tolerance for iterative scheme + + # get value at current z via interpolation + dz = height / (nz - 1) + k_l = convert(Int, floor(x[2] / dz)) + 1 + z_l = (k_l - 1) * dz + z_u = k_l * dz + if k_l == nz + k_u = nz + else + k_u = k_l + 1 + end + theta, theta_rho, p, q_v, q_l = map((a, b) -> (z_u - x[2]) * a / dz + + (x[2] - z_l) * b / dz, + [ + theta_0[k_l], + theta_rho_0[k_l], + p_0[k_l], + q_v0[k_l], + q_l0[k_l], + ], + [ + theta_0[k_u], + theta_rho_0[k_u], + p_0[k_u], + q_v0[k_u], + q_l0[k_u], + ]) + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + exner = (p / p_ref)^(R_d / c_pd) + T = theta * exner + + if r <= radius + theta_perturbation = 2 * cospi(0.5 * r / radius)^2 + theta_tmp = theta + converged = false + steps = 0 + while !converged + steps = steps + 1 + theta = ((theta_perturbation / 300.0) + (1.0 + r_t) / (1.0 + q_v)) * + theta_rho * (1.0 + q_v) / (1.0 + inv(eps) * q_v) + T = theta * exner + q_v = getqvs(eps, p, T) + q_l = r_t - q_v + if abs(theta - theta_tmp) < tolerance + converged = true + if q_l < 0.0 + error("calculate_hydrostatic_base_state: liquid ratio negative!") + end + else + theta_tmp = theta + if steps > 100 + error("initial condition perturbation: not converging!") + end + end + end + end + + rho_d = p / (R_d * T * (1.0 + q_v * inv(eps))) # p = rho_d * R_d * T * (1 + q_v / eps) + rho_v = rho_d * q_v + rho_l = rho_d * q_l + rho = rho_d + rho_l + rho_v + v1 = 0.0 + v2 = 0.0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_e = (c_vd * rho_d + c_vv * rho_v + c_pl * rho_l) * T + + L_00 * rho_v + + 0.5 * rho * (v1^2 + v2^2) + return SVector(rho_v1, rho_v2, rho_e, rho_d, rho_v, rho_l) +end + +# Initial condition Knoth struct AtmossphereLayers{RealT <: Real} setup::WarmBubbleSetup # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql @@ -150,7 +375,7 @@ function AtmossphereLayers(setup; total_height = 10010.0, preciseness = 10, end function moist_state(y, dz, y0, r_t0, theta_e0, setup::WarmBubbleSetup) - @unpack g, c_pd, R_d, c_pv, R_v, c_pl, p_0, L_00 = setup + @unpack g, c_pd, R_d, c_pv, R_v, c_pl, p_ref, theta_ref, L_00 = setup (p, rho, T, r_t, r_v, rho_qv, theta_e) = y p0 = y0[1] @@ -165,7 +390,7 @@ function moist_state(y, dz, y0, r_t0, theta_e0, setup::WarmBubbleSetup) F[2] = p - (R_d * rho_d + R_v * rho_qv) * T # H = 1 is assumed F[3] = (theta_e - - T * (p_d / p_0)^(-R_d / (c_pd + c_pl * r_t)) * + T * (p_d / p_ref)^(-R_d / (c_pd + c_pl * r_t)) * exp(L * r_v / ((c_pd + c_pl * r_t) * T))) F[4] = r_t - r_t0 F[5] = rho_qv - rho_d * r_v @@ -184,10 +409,6 @@ function generate_function_of_y(dz, y0, r_t0, theta_e0, setup::WarmBubbleSetup) end end -# Moist bubble test case from paper: -# G.H. Bryan, J.M. Fritsch, A Benchmark Simulation for Moist Nonhydrostatic Numerical -# Models, MonthlyWeather Review Vol.130, 2917–2928, 2002, -# https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml. function initial_condition_moist_bubble(x, t, setup::WarmBubbleSetup, AtmosphereLayers::AtmossphereLayers) @unpack LayerData, preciseness, total_height = AtmosphereLayers @@ -226,7 +447,7 @@ function initial_condition_moist_bubble(x, t, setup::WarmBubbleSetup, end function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbleSetup) - @unpack c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_0, L_00 = setup + @unpack c_pd, c_vd, R_d, c_pv, c_vv, R_v, c_pl, p_ref, L_00 = setup xc = 10000.0 zc = 2000.0 rc = 2000.0 @@ -235,7 +456,7 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbl r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) rho_d = rho - rho_qv - rho_ql kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) - p_loc = p_0 * (R_d * rho_theta / p_0)^(1 / (1 - kappa_M)) + p_loc = p_ref * (R_d * rho_theta / p_ref)^(1 / (1 - kappa_M)) T_loc = p_loc / (R_d * rho_d + R_v * rho_qv) rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv @@ -249,7 +470,7 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbl r_t = r_v + r_l # Aequivalentpotential temperature - a = T_loc * (p_0 / p_d)^(R_d / (c_pd + r_t * c_pl)) + a = T_loc * (p_ref / p_d)^(R_d / (c_pd + r_t * c_pl)) b = H^(-r_v * R_v / c_pd) L_v = L_00 + (c_pv - c_pl) * T_loc c = exp(L_v * r_v / ((c_pd + r_t * c_pl) * T_loc)) @@ -259,7 +480,7 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbl if (r < rc && Δθ > 0) kappa = 1 - c_vd / c_pd # Calculate background density potential temperature - θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) + θ_dens = rho_theta / rho * (p_loc / p_ref)^(kappa_M - kappa) # Calculate perturbed density potential temperature θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5 * r / rc)^2 / 300) rt = (rho_qv + rho_ql) / rho_d @@ -269,7 +490,7 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbl # Adjust varuables until the temperature is met if rt > 0 while true - T_loc = θ_loc * (p_loc / p_0)^kappa + T_loc = θ_loc * (p_loc / p_ref)^kappa T_C = T_loc - 273.15 # SaturVapor pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) @@ -284,7 +505,7 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbl end else rvs = 0 - T_loc = θ_loc * (p_loc / p_0)^kappa + T_loc = θ_loc * (p_loc / p_ref)^kappa rho_d_new = p_loc / (R_d * T_loc) θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) end @@ -294,7 +515,7 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, setup::WarmBubbl rho_d = rho - rho_qv - rho_ql kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) - rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M) + rho_theta = rho * θ_dens_new * (p_loc / p_ref)^(kappa - kappa_M) rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv end return SVector(rho, rho_e, rho_qv, rho_ql) @@ -302,7 +523,11 @@ end ############################################################################### # semidiscretization of the compressible moist Euler equations -warm_bubble_setup = WarmBubbleSetup() +nelements_z = 100 +polydeg = 1 +height = 10_000.0 + +warm_bubble_setup = WarmBubbleSetup(; nz = nelements_z * (polydeg + 1), height = height) AtmossphereData = AtmossphereLayers(warm_bubble_setup) @@ -320,20 +545,21 @@ boundary_condition = (x_neg = boundary_condition_slip_wall, y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) -polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = FluxLMARS(360.0) -volume_flux = flux_chandrashekar +#surface_flux = FluxLMARS(340.0) +surface_flux = flux_lax_friedrichs +#volume_flux = flux_chandrashekar #volume_flux = flux_ranocha +volume_flux = flux_kennedy_gruber volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) -coordinates_max = (20000.0, 10000.0) +coordinates_max = (20_000.0, height) -cells_per_dimension = (64, 32) +cells_per_dimension = (2 * nelements_z, nelements_z) # Create curved mesh with 64 x 32 elements mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) @@ -341,14 +567,14 @@ mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) ############################################################################### # create the semi discretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_moist, solver, +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, boundary_conditions = boundary_condition, source_terms = warm_bubble_setup) ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 1000.0) +tspan = (0.0, 1200.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -365,12 +591,13 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true, + output_directory = "out_bf_struct200_p1_cfl02_lf_kennedy_gruber", solution_variables = solution_variables) stepsize_callback = StepsizeCallback(cfl = 0.2) callbacks = CallbackSet(summary_callback, - analysis_callback, + #analysis_callback, alive_callback, save_solution, stepsize_callback) @@ -380,5 +607,7 @@ callbacks = CallbackSet(summary_callback, sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); + save_everystep = false, + maxiters = 1e6, + callback = callbacks); summary_callback() # print the timer summary