From e12a357a041758966ca4904ff694707591900310 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Thu, 19 Sep 2024 22:28:53 +0200 Subject: [PATCH 01/27] Potential Temperature Euler Eq. --- examples/elixir_potential_euler_dry_bubble.jl | 114 +++++++ examples/julia | 0 src/TrixiAtmo.jl | 3 +- .../compressible_potential_euler_2d.jl | 288 ++++++++++++++++++ src/equations/equations.jl | 3 + 5 files changed, 407 insertions(+), 1 deletion(-) create mode 100644 examples/elixir_potential_euler_dry_bubble.jl create mode 100644 examples/julia create mode 100644 src/equations/compressible_potential_euler_2d.jl diff --git a/examples/elixir_potential_euler_dry_bubble.jl b/examples/elixir_potential_euler_dry_bubble.jl new file mode 100644 index 0000000..adcbf34 --- /dev/null +++ b/examples/elixir_potential_euler_dry_bubble.jl @@ -0,0 +1,114 @@ +using OrdinaryDiffEq +using Trixi +using TrixiAtmo +using TrixiAtmo: flux_LMARS, source_terms_bubble, flux_theta +using Plots + +function initial_condition_warm_bubble(x, t, equations::CompressiblePotentialEulerEquations2D) + g = equations.g + c_p = equations.c_p + c_v = equations.c_v + # 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_ref) * 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) + T = potential_temperature * exner + + # density + rho = p / (R * T) + v1 = 20.0 + v2 = 0.0 + + return SVector(rho, rho * v1, rho * v2, rho * potential_temperature) +end + +############################################################################### +# semidiscretization of the compressible Euler equations + + +equations = CompressiblePotentialEulerEquations2D() + +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 = 3 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = flux_LMARS + +volume_flux = flux_theta +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 = (32, 16) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) +initial_condition = initial_condition_warm_bubble + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_bubble, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time + +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 = 1.0) + +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/examples/julia b/examples/julia new file mode 100644 index 0000000..e69de29 diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index e6c4608..49eeb8c 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -24,8 +24,9 @@ include("solvers/solvers.jl") include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D +export CompressiblePotentialEulerEquations2D -export flux_chandrashekar, flux_LMARS +export flux_chandrashekar, flux_LMARS, flux_theta export examples_dir diff --git a/src/equations/compressible_potential_euler_2d.jl b/src/equations/compressible_potential_euler_2d.jl new file mode 100644 index 0000000..3d69c0d --- /dev/null +++ b/src/equations/compressible_potential_euler_2d.jl @@ -0,0 +1,288 @@ +using Trixi +using Trixi: ln_mean, stolarsky_mean +import Trixi: varnames, cons2cons, cons2prim, cons2entropy, entropy + +@muladd begin +#! format: noindent +struct CompressiblePotentialEulerEquations2D{RealT <: Real} <: + AbstractCompressiblePotentialEulerEquations{2, 4} + p_0::RealT + c_p::RealT + c_v::RealT + g::RealT + R::RealT + gamma::RealT + a::RealT +end + +function CompressiblePotentialEulerEquations2D(; g = 9.81, RealT = Float64) + p_0 = 100_000.0 + c_p = 1004.0 + c_v = 717.0 + R = c_p - c_v + gamma = c_p / c_v + a = 340.0 + return CompressiblePotentialEulerEquations2D{RealT}(p_0, c_p, c_v, g, R, gamma, a) +end + +function varnames(::typeof(cons2cons), ::CompressiblePotentialEulerEquations2D) + ("rho", "rho_v1", "rho_v2", "rho_theta") +end + +varnames(::typeof(cons2prim), ::CompressiblePotentialEulerEquations2D) = ("rho", "v1", "v2", "p1") + +# Calculate 1D flux for a single point in the normal direction. +# Note, this directional vector is not normalized. +@inline function flux(u, normal_direction::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = (rho_v_normal) * v1 + p * normal_direction[1] + f3 = (rho_v_normal) * v2 + p * normal_direction[2] + f4 = (rho_theta) * v_normal + return SVector(f1, f2, f3, f4) +end + +@inline function flux(u, orientation::Integer, equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma + + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + f4 = (rho_theta) * v1 + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + f4 = (rho_theta) * v2 + end + + return SVector(f1, f2, f3, f4) +end + +# Slip-wall boundary condition +# Determine the boundary numerical surface flux for a slip wall condition. +# Imposes a zero normal velocity at the wall. +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressiblePotentialEulerEquations2D) + @unpack gamma = equations + 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 + rho_local, v_normal, v_tangent, p_local = cons2prim(u_local, equations) + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0.0 + sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * + inv(gamma - 1)) + else # v_normal > 0.0 + A = 2.0 / ((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.0 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return SVector(zero(eltype(u_inner)), + p_star * normal[1], + p_star * normal[2], + zero(eltype(u_inner))) * norm_ +end + + +# Fix sign for structured mesh. +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressiblePotentialEulerEquations2D) + # 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 + +@inline function source_terms_bubble(u, x, t, equations::CompressiblePotentialEulerEquations2D) + rho, _, _, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -equations.g * rho, zero(eltype(u))) +end + +# Rotate momentum flux. The same as in compressible Euler. +@inline function rotate_to_x(u, normal_vector::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + # 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 + # [ 1 0 0 0; + # 0 n_1 n_2 0; + # 0 t_1 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] + s * u[3], + -s * u[2] + c * u[3], + u[4]) +end + +# Low Mach number approximate Riemann solver (LMARS) from +# X. Chen, N. Andronova, B. Van Leer, J. E. Penner, J. P. Boyd, C. Jablonowski, S. +# Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian +# Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, +# https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. +@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + @unpack a = equations + # 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] + + + # Compute the necessary interface flux components + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + + p_interface = 0.5 * (p_ll + p_rr) - 0.5 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + + if (v_interface > 0) + f1, f2, f3, f4 = u_ll * v_interface + else + f1, f2, f3, f4 = u_rr * v_interface + end + + return SVector(f1, + f2 + p_interface * normal_direction[1], + f3 + p_interface * normal_direction[2], + f4) +end + +## Entropy (total energy) conservative flux for the Compressible Euler with the Potential Formulation +@inline function flux_theta(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + # 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_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_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` + # in exact arithmetic since + # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) + # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = gammamean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + return SVector(f1, f2, f3, f4) +end + +@inline function prim2cons(prim, equations::CompressiblePotentialEulerEquations2D) + rho, v1, v2, p = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_theta = (p / equations.p_0)^(1 / equations.gamma) * equations.p_0 / equations.R + return SVector(rho, rho_v1, rho_v2, rho_theta) +end + +@inline function cons2prim(u, equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma + + return SVector(rho, v1, v2, p) +end + +@inline function cons2cons(u, equations::CompressiblePotentialEulerEquations2D) + return u +end + +@inline function cons2entropy(u, equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + + k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma + w1 = -0.5 * rho_v1^2 / (rho)^2 - 0.5 * rho_v2^2 / (rho)^2 + w2 = rho_v1 / rho + w3 = rho_v2 / rho + w4 = equations.gamma / (equations.gamma - 1) * k * (rho_theta)^(equations.gamma - 1) + + return SVector(w1, w2, w3, w4) +end + +@inline function entropy_math(cons, equations::CompressiblePotentialEulerEquations2D) + # Mathematical entropy + p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma + + U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2 + cons[3]^2) / (cons[1])) + + return U + +end + +# Default entropy is the mathematical entropy +@inline function entropy(cons, equations::CompressiblePotentialEulerEquations2D) + entropy_math(cons, equations) +end + +@inline function energy_total(cons, equations::CompressiblePotentialEulerEquations2D) + entropy(cons, equations) +end + +@inline function energy_kinetic(cons, equations::CompressiblePotentialEulerEquations2D) + return 0.5 * (cons[2]^2 + cons[3]^2) / (cons[1]) +end + +@inline function max_abs_speeds(u, equations::CompressiblePotentialEulerEquations2D) + rho, v1, v2, p = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c +end + +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a833960..9eb34eb 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -2,4 +2,7 @@ using Trixi: AbstractEquations abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +abstract type AbstractCompressiblePotentialEulerEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end include("compressible_moist_euler_2d_lucas.jl") +include("compressible_potential_euler_2d.jl") \ No newline at end of file From 4602b1ce1248911c3260fdbd35525ab96e9662ce Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Thu, 19 Sep 2024 22:34:20 +0200 Subject: [PATCH 02/27] delete backup files --- examples/julia | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 examples/julia diff --git a/examples/julia b/examples/julia deleted file mode 100644 index e69de29..0000000 From 680186eda412cacc4a3aa333d18a0afd29574b8c Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Thu, 19 Sep 2024 22:45:29 +0200 Subject: [PATCH 03/27] formatting --- examples/elixir_potential_euler_dry_bubble.jl | 228 +++---- .../compressible_potential_euler_2d.jl | 576 +++++++++--------- src/equations/equations.jl | 4 +- 3 files changed, 404 insertions(+), 404 deletions(-) diff --git a/examples/elixir_potential_euler_dry_bubble.jl b/examples/elixir_potential_euler_dry_bubble.jl index adcbf34..da28fdb 100644 --- a/examples/elixir_potential_euler_dry_bubble.jl +++ b/examples/elixir_potential_euler_dry_bubble.jl @@ -1,114 +1,114 @@ -using OrdinaryDiffEq -using Trixi -using TrixiAtmo -using TrixiAtmo: flux_LMARS, source_terms_bubble, flux_theta -using Plots - -function initial_condition_warm_bubble(x, t, equations::CompressiblePotentialEulerEquations2D) - g = equations.g - c_p = equations.c_p - c_v = equations.c_v - # 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_ref) * 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) - T = potential_temperature * exner - - # density - rho = p / (R * T) - v1 = 20.0 - v2 = 0.0 - - return SVector(rho, rho * v1, rho * v2, rho * potential_temperature) -end - -############################################################################### -# semidiscretization of the compressible Euler equations - - -equations = CompressiblePotentialEulerEquations2D() - -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 = 3 -basis = LobattoLegendreBasis(polydeg) - -surface_flux = flux_LMARS - -volume_flux = flux_theta -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 = (32, 16) -mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, - periodicity = (true, false)) -initial_condition = initial_condition_warm_bubble - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms = source_terms_bubble, - boundary_conditions = boundary_conditions) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 1000.0) # 1000 seconds final time - -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 = 1.0) - -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() +using OrdinaryDiffEq +using Trixi +using TrixiAtmo +using TrixiAtmo: flux_LMARS, source_terms_bubble, flux_theta +using Plots + +function initial_condition_warm_bubble(x, t, + equations::CompressiblePotentialEulerEquations2D) + g = equations.g + c_p = equations.c_p + c_v = equations.c_v + # 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_ref) * 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) + T = potential_temperature * exner + + # density + rho = p / (R * T) + v1 = 20.0 + v2 = 0.0 + + return SVector(rho, rho * v1, rho * v2, rho * potential_temperature) +end + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressiblePotentialEulerEquations2D() + +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 = 3 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = flux_LMARS + +volume_flux = flux_theta +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 = (32, 16) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) +initial_condition = initial_condition_warm_bubble + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_bubble, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time + +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 = 1.0) + +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/src/equations/compressible_potential_euler_2d.jl b/src/equations/compressible_potential_euler_2d.jl index 3d69c0d..81b031b 100644 --- a/src/equations/compressible_potential_euler_2d.jl +++ b/src/equations/compressible_potential_euler_2d.jl @@ -1,288 +1,288 @@ -using Trixi -using Trixi: ln_mean, stolarsky_mean -import Trixi: varnames, cons2cons, cons2prim, cons2entropy, entropy - -@muladd begin -#! format: noindent -struct CompressiblePotentialEulerEquations2D{RealT <: Real} <: - AbstractCompressiblePotentialEulerEquations{2, 4} - p_0::RealT - c_p::RealT - c_v::RealT - g::RealT - R::RealT - gamma::RealT - a::RealT -end - -function CompressiblePotentialEulerEquations2D(; g = 9.81, RealT = Float64) - p_0 = 100_000.0 - c_p = 1004.0 - c_v = 717.0 - R = c_p - c_v - gamma = c_p / c_v - a = 340.0 - return CompressiblePotentialEulerEquations2D{RealT}(p_0, c_p, c_v, g, R, gamma, a) -end - -function varnames(::typeof(cons2cons), ::CompressiblePotentialEulerEquations2D) - ("rho", "rho_v1", "rho_v2", "rho_theta") -end - -varnames(::typeof(cons2prim), ::CompressiblePotentialEulerEquations2D) = ("rho", "v1", "v2", "p1") - -# Calculate 1D flux for a single point in the normal direction. -# Note, this directional vector is not normalized. -@inline function flux(u, normal_direction::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) - rho, rho_v1, rho_v2, rho_theta = u - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] - rho_v_normal = rho * v_normal - f1 = rho_v_normal - f2 = (rho_v_normal) * v1 + p * normal_direction[1] - f3 = (rho_v_normal) * v2 + p * normal_direction[2] - f4 = (rho_theta) * v_normal - return SVector(f1, f2, f3, f4) -end - -@inline function flux(u, orientation::Integer, equations::CompressiblePotentialEulerEquations2D) - rho, rho_v1, rho_v2, rho_theta = u - v1 = rho_v1 / rho - v2 = rho_v2 / rho - p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma - - if orientation == 1 - f1 = rho_v1 - f2 = rho_v1 * v1 + p - f3 = rho_v1 * v2 - f4 = (rho_theta) * v1 - else - f1 = rho_v2 - f2 = rho_v2 * v1 - f3 = rho_v2 * v2 + p - f4 = (rho_theta) * v2 - end - - return SVector(f1, f2, f3, f4) -end - -# Slip-wall boundary condition -# Determine the boundary numerical surface flux for a slip wall condition. -# Imposes a zero normal velocity at the wall. -@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, - x, t, - surface_flux_function, - equations::CompressiblePotentialEulerEquations2D) - @unpack gamma = equations - 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 - rho_local, v_normal, v_tangent, p_local = cons2prim(u_local, equations) - # Get the solution of the pressure Riemann problem - # See Section 6.3.3 of - # Eleuterio F. Toro (2009) - # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction - # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) - if v_normal <= 0.0 - sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed - p_star = p_local * - (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * - inv(gamma - 1)) - else # v_normal > 0.0 - A = 2.0 / ((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.0 * A * (p_local + B))) - end - - # For the slip wall we directly set the flux as the normal velocity is zero - return SVector(zero(eltype(u_inner)), - p_star * normal[1], - p_star * normal[2], - zero(eltype(u_inner))) * norm_ -end - - -# Fix sign for structured mesh. -@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, - direction, x, t, - surface_flux_function, - equations::CompressiblePotentialEulerEquations2D) - # 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 - -@inline function source_terms_bubble(u, x, t, equations::CompressiblePotentialEulerEquations2D) - rho, _, _, _ = u - return SVector(zero(eltype(u)), zero(eltype(u)), -equations.g * rho, zero(eltype(u))) -end - -# Rotate momentum flux. The same as in compressible Euler. -@inline function rotate_to_x(u, normal_vector::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) - # 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 - # [ 1 0 0 0; - # 0 n_1 n_2 0; - # 0 t_1 t_2 0; - # 0 0 0 1 ] - # where t_1 = -n_2 and t_2 = n_1 - - return SVector(u[1], - c * u[2] + s * u[3], - -s * u[2] + c * u[3], - u[4]) -end - -# Low Mach number approximate Riemann solver (LMARS) from -# X. Chen, N. Andronova, B. Van Leer, J. E. Penner, J. P. Boyd, C. Jablonowski, S. -# Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian -# Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, -# https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. -@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) - @unpack a = equations - # 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] - - - # Compute the necessary interface flux components - norm_ = norm(normal_direction) - - rho = 0.5 * (rho_ll + rho_rr) - - p_interface = 0.5 * (p_ll + p_rr) - 0.5 * a * rho * (v_rr - v_ll) / norm_ - v_interface = 0.5 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ - - if (v_interface > 0) - f1, f2, f3, f4 = u_ll * v_interface - else - f1, f2, f3, f4 = u_rr * v_interface - end - - return SVector(f1, - f2 + p_interface * normal_direction[1], - f3 + p_interface * normal_direction[2], - f4) -end - -## Entropy (total energy) conservative flux for the Compressible Euler with the Potential Formulation -@inline function flux_theta(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) - # 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_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_theta_ll = u_ll - _, _, _, rho_theta_rr = u_rr - # Compute the necessary mean values - rho_mean = ln_mean(rho_ll, rho_rr) - gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) - # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` - # in exact arithmetic since - # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) - # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) - - # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) - f2 = f1 * v1_avg + p_avg * normal_direction[1] - f3 = f1 * v2_avg + p_avg * normal_direction[2] - f4 = gammamean * 0.5 * (v_dot_n_ll + v_dot_n_rr) - return SVector(f1, f2, f3, f4) -end - -@inline function prim2cons(prim, equations::CompressiblePotentialEulerEquations2D) - rho, v1, v2, p = prim - rho_v1 = rho * v1 - rho_v2 = rho * v2 - rho_theta = (p / equations.p_0)^(1 / equations.gamma) * equations.p_0 / equations.R - return SVector(rho, rho_v1, rho_v2, rho_theta) -end - -@inline function cons2prim(u, equations::CompressiblePotentialEulerEquations2D) - rho, rho_v1, rho_v2, rho_theta = u - v1 = rho_v1 / rho - v2 = rho_v2 / rho - p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma - - return SVector(rho, v1, v2, p) -end - -@inline function cons2cons(u, equations::CompressiblePotentialEulerEquations2D) - return u -end - -@inline function cons2entropy(u, equations::CompressiblePotentialEulerEquations2D) - rho, rho_v1, rho_v2, rho_theta = u - - k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma - w1 = -0.5 * rho_v1^2 / (rho)^2 - 0.5 * rho_v2^2 / (rho)^2 - w2 = rho_v1 / rho - w3 = rho_v2 / rho - w4 = equations.gamma / (equations.gamma - 1) * k * (rho_theta)^(equations.gamma - 1) - - return SVector(w1, w2, w3, w4) -end - -@inline function entropy_math(cons, equations::CompressiblePotentialEulerEquations2D) - # Mathematical entropy - p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma - - U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2 + cons[3]^2) / (cons[1])) - - return U - -end - -# Default entropy is the mathematical entropy -@inline function entropy(cons, equations::CompressiblePotentialEulerEquations2D) - entropy_math(cons, equations) -end - -@inline function energy_total(cons, equations::CompressiblePotentialEulerEquations2D) - entropy(cons, equations) -end - -@inline function energy_kinetic(cons, equations::CompressiblePotentialEulerEquations2D) - return 0.5 * (cons[2]^2 + cons[3]^2) / (cons[1]) -end - -@inline function max_abs_speeds(u, equations::CompressiblePotentialEulerEquations2D) - rho, v1, v2, p = cons2prim(u, equations) - c = sqrt(equations.gamma * p / rho) - - return abs(v1) + c, abs(v2) + c -end - -end # @muladd +using Trixi +using Trixi: ln_mean, stolarsky_mean +import Trixi: varnames, cons2cons, cons2prim, cons2entropy, entropy + +@muladd begin +#! format: noindent +struct CompressiblePotentialEulerEquations2D{RealT <: Real} <: + AbstractCompressiblePotentialEulerEquations{2, 4} + p_0::RealT + c_p::RealT + c_v::RealT + g::RealT + R::RealT + gamma::RealT + a::RealT +end + +function CompressiblePotentialEulerEquations2D(; g = 9.81, RealT = Float64) + p_0 = 100_000.0 + c_p = 1004.0 + c_v = 717.0 + R = c_p - c_v + gamma = c_p / c_v + a = 340.0 + return CompressiblePotentialEulerEquations2D{RealT}(p_0, c_p, c_v, g, R, gamma, a) +end + +function varnames(::typeof(cons2cons), ::CompressiblePotentialEulerEquations2D) + ("rho", "rho_v1", "rho_v2", "rho_theta") +end + +varnames(::typeof(cons2prim), ::CompressiblePotentialEulerEquations2D) = ("rho", "v1", + "v2", "p1") + +# Calculate 1D flux for a single point in the normal direction. +# Note, this directional vector is not normalized. +@inline function flux(u, normal_direction::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = (rho_v_normal) * v1 + p * normal_direction[1] + f3 = (rho_v_normal) * v2 + p * normal_direction[2] + f4 = (rho_theta) * v_normal + return SVector(f1, f2, f3, f4) +end + +@inline function flux(u, orientation::Integer, + equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma + + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + f4 = (rho_theta) * v1 + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + f4 = (rho_theta) * v2 + end + + return SVector(f1, f2, f3, f4) +end + +# Slip-wall boundary condition +# Determine the boundary numerical surface flux for a slip wall condition. +# Imposes a zero normal velocity at the wall. +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressiblePotentialEulerEquations2D) + @unpack gamma = equations + 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 + rho_local, v_normal, v_tangent, p_local = cons2prim(u_local, equations) + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0.0 + sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * + inv(gamma - 1)) + else # v_normal > 0.0 + A = 2.0 / ((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.0 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return SVector(zero(eltype(u_inner)), + p_star * normal[1], + p_star * normal[2], + zero(eltype(u_inner))) * norm_ +end + +# Fix sign for structured mesh. +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressiblePotentialEulerEquations2D) + # 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 + +@inline function source_terms_bubble(u, x, t, + equations::CompressiblePotentialEulerEquations2D) + rho, _, _, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -equations.g * rho, + zero(eltype(u))) +end + +# Rotate momentum flux. The same as in compressible Euler. +@inline function rotate_to_x(u, normal_vector::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + # 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 + # [ 1 0 0 0; + # 0 n_1 n_2 0; + # 0 t_1 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] + s * u[3], + -s * u[2] + c * u[3], + u[4]) +end + +# Low Mach number approximate Riemann solver (LMARS) from +# X. Chen, N. Andronova, B. Van Leer, J. E. Penner, J. P. Boyd, C. Jablonowski, S. +# Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian +# Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, +# https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. +@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + @unpack a = equations + # 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] + + # Compute the necessary interface flux components + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + + p_interface = 0.5 * (p_ll + p_rr) - 0.5 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + + if (v_interface > 0) + f1, f2, f3, f4 = u_ll * v_interface + else + f1, f2, f3, f4 = u_rr * v_interface + end + + return SVector(f1, + f2 + p_interface * normal_direction[1], + f3 + p_interface * normal_direction[2], + f4) +end + +## Entropy (total energy) conservative flux for the Compressible Euler with the Potential Formulation +@inline function flux_theta(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressiblePotentialEulerEquations2D) + # 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_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_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` + # in exact arithmetic since + # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) + # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = gammamean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + return SVector(f1, f2, f3, f4) +end + +@inline function prim2cons(prim, equations::CompressiblePotentialEulerEquations2D) + rho, v1, v2, p = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_theta = (p / equations.p_0)^(1 / equations.gamma) * equations.p_0 / equations.R + return SVector(rho, rho_v1, rho_v2, rho_theta) +end + +@inline function cons2prim(u, equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = equations.p_0 * (equations.R * rho_theta / equations.p_0)^equations.gamma + + return SVector(rho, v1, v2, p) +end + +@inline function cons2cons(u, equations::CompressiblePotentialEulerEquations2D) + return u +end + +@inline function cons2entropy(u, equations::CompressiblePotentialEulerEquations2D) + rho, rho_v1, rho_v2, rho_theta = u + + k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma + w1 = -0.5 * rho_v1^2 / (rho)^2 - 0.5 * rho_v2^2 / (rho)^2 + w2 = rho_v1 / rho + w3 = rho_v2 / rho + w4 = equations.gamma / (equations.gamma - 1) * k * (rho_theta)^(equations.gamma - 1) + + return SVector(w1, w2, w3, w4) +end + +@inline function entropy_math(cons, equations::CompressiblePotentialEulerEquations2D) + # Mathematical entropy + p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma + + U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2 + cons[3]^2) / (cons[1])) + + return U +end + +# Default entropy is the mathematical entropy +@inline function entropy(cons, equations::CompressiblePotentialEulerEquations2D) + entropy_math(cons, equations) +end + +@inline function energy_total(cons, equations::CompressiblePotentialEulerEquations2D) + entropy(cons, equations) +end + +@inline function energy_kinetic(cons, equations::CompressiblePotentialEulerEquations2D) + return 0.5 * (cons[2]^2 + cons[3]^2) / (cons[1]) +end + +@inline function max_abs_speeds(u, equations::CompressiblePotentialEulerEquations2D) + rho, v1, v2, p = cons2prim(u, equations) + c = sqrt(equations.gamma * p / rho) + + return abs(v1) + c, abs(v2) + c +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 9eb34eb..a97393c 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -3,6 +3,6 @@ using Trixi: AbstractEquations abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end abstract type AbstractCompressiblePotentialEulerEquations{NDIMS, NVARS} <: - AbstractEquations{NDIMS, NVARS} end + AbstractEquations{NDIMS, NVARS} end include("compressible_moist_euler_2d_lucas.jl") -include("compressible_potential_euler_2d.jl") \ No newline at end of file +include("compressible_potential_euler_2d.jl") From dd5568ad3ced9b647dced0481d5446268842bfb5 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Thu, 19 Sep 2024 22:46:14 +0200 Subject: [PATCH 04/27] formatting --- examples/julia | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 examples/julia diff --git a/examples/julia b/examples/julia deleted file mode 100644 index e69de29..0000000 From 04c304d01574071711348e3ec6f3932cbf8544c7 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 14:36:15 +0200 Subject: [PATCH 05/27] renaming --- examples/elixir_potential_euler_dry_bubble.jl | 6 +-- src/TrixiAtmo.jl | 2 +- ...essible_euler_potential_temperature_2d.jl} | 46 +++++++++---------- src/equations/equations.jl | 2 +- 4 files changed, 28 insertions(+), 28 deletions(-) rename src/equations/{compressible_potential_euler_2d.jl => compressible_euler_potential_temperature_2d.jl} (82%) diff --git a/examples/elixir_potential_euler_dry_bubble.jl b/examples/elixir_potential_euler_dry_bubble.jl index da28fdb..288a797 100644 --- a/examples/elixir_potential_euler_dry_bubble.jl +++ b/examples/elixir_potential_euler_dry_bubble.jl @@ -5,7 +5,7 @@ using TrixiAtmo: flux_LMARS, source_terms_bubble, flux_theta using Plots function initial_condition_warm_bubble(x, t, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) g = equations.g c_p = equations.c_p c_v = equations.c_v @@ -45,7 +45,7 @@ end ############################################################################### # semidiscretization of the compressible Euler equations -equations = CompressiblePotentialEulerEquations2D() +equations = CompressibleEulerPotentialTemperatureEquations2D() boundary_conditions = (x_neg = boundary_condition_periodic, x_pos = boundary_condition_periodic, @@ -65,7 +65,7 @@ solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (20_000.0, 10_000.0) -cells_per_dimension = (32, 16) +cells_per_dimension = (64, 32) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, periodicity = (true, false)) initial_condition = initial_condition_warm_bubble diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 49eeb8c..cdce0a3 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -24,7 +24,7 @@ include("solvers/solvers.jl") include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D -export CompressiblePotentialEulerEquations2D +export CompressibleEulerPotentialTemperatureEquations2D export flux_chandrashekar, flux_LMARS, flux_theta diff --git a/src/equations/compressible_potential_euler_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl similarity index 82% rename from src/equations/compressible_potential_euler_2d.jl rename to src/equations/compressible_euler_potential_temperature_2d.jl index 81b031b..8ddea9f 100644 --- a/src/equations/compressible_potential_euler_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -4,8 +4,8 @@ import Trixi: varnames, cons2cons, cons2prim, cons2entropy, entropy @muladd begin #! format: noindent -struct CompressiblePotentialEulerEquations2D{RealT <: Real} <: - AbstractCompressiblePotentialEulerEquations{2, 4} +struct CompressibleEulerPotentialTemperatureEquations2D{RealT <: Real} <: + AbstractCompressibleEulerPotentialTemperatureEquations{2, 4} p_0::RealT c_p::RealT c_v::RealT @@ -15,27 +15,27 @@ struct CompressiblePotentialEulerEquations2D{RealT <: Real} <: a::RealT end -function CompressiblePotentialEulerEquations2D(; g = 9.81, RealT = Float64) +function CompressibleEulerPotentialTemperatureEquations2D(; g = 9.81, RealT = Float64) p_0 = 100_000.0 c_p = 1004.0 c_v = 717.0 R = c_p - c_v gamma = c_p / c_v a = 340.0 - return CompressiblePotentialEulerEquations2D{RealT}(p_0, c_p, c_v, g, R, gamma, a) + return CompressibleEulerPotentialTemperatureEquations2D{RealT}(p_0, c_p, c_v, g, R, gamma, a) end -function varnames(::typeof(cons2cons), ::CompressiblePotentialEulerEquations2D) +function varnames(::typeof(cons2cons), ::CompressibleEulerPotentialTemperatureEquations2D) ("rho", "rho_v1", "rho_v2", "rho_theta") end -varnames(::typeof(cons2prim), ::CompressiblePotentialEulerEquations2D) = ("rho", "v1", +varnames(::typeof(cons2prim), ::CompressibleEulerPotentialTemperatureEquations2D) = ("rho", "v1", "v2", "p1") # Calculate 1D flux for a single point in the normal direction. # Note, this directional vector is not normalized. @inline function flux(u, normal_direction::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, rho_v1, rho_v2, rho_theta = u v1 = rho_v1 / rho v2 = rho_v2 / rho @@ -49,7 +49,7 @@ varnames(::typeof(cons2prim), ::CompressiblePotentialEulerEquations2D) = ("rho", end @inline function flux(u, orientation::Integer, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, rho_v1, rho_v2, rho_theta = u v1 = rho_v1 / rho v2 = rho_v2 / rho @@ -76,7 +76,7 @@ end @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) @unpack gamma = equations norm_ = norm(normal_direction) # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later @@ -116,7 +116,7 @@ end @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, direction, x, t, surface_flux_function, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) # 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) @@ -133,7 +133,7 @@ end end @inline function source_terms_bubble(u, x, t, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, _, _, _ = u return SVector(zero(eltype(u)), zero(eltype(u)), -equations.g * rho, zero(eltype(u))) @@ -141,7 +141,7 @@ end # Rotate momentum flux. The same as in compressible Euler. @inline function rotate_to_x(u, normal_vector::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) # 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] @@ -166,7 +166,7 @@ end # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. @inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) @unpack a = equations # Unpack left and right state rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) @@ -197,7 +197,7 @@ end ## Entropy (total energy) conservative flux for the Compressible Euler with the Potential Formulation @inline function flux_theta(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressiblePotentialEulerEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) # 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) @@ -224,7 +224,7 @@ end return SVector(f1, f2, f3, f4) end -@inline function prim2cons(prim, equations::CompressiblePotentialEulerEquations2D) +@inline function prim2cons(prim, equations::CompressibleEulerPotentialTemperatureEquations2D) rho, v1, v2, p = prim rho_v1 = rho * v1 rho_v2 = rho * v2 @@ -232,7 +232,7 @@ end return SVector(rho, rho_v1, rho_v2, rho_theta) end -@inline function cons2prim(u, equations::CompressiblePotentialEulerEquations2D) +@inline function cons2prim(u, equations::CompressibleEulerPotentialTemperatureEquations2D) rho, rho_v1, rho_v2, rho_theta = u v1 = rho_v1 / rho v2 = rho_v2 / rho @@ -241,11 +241,11 @@ end return SVector(rho, v1, v2, p) end -@inline function cons2cons(u, equations::CompressiblePotentialEulerEquations2D) +@inline function cons2cons(u, equations::CompressibleEulerPotentialTemperatureEquations2D) return u end -@inline function cons2entropy(u, equations::CompressiblePotentialEulerEquations2D) +@inline function cons2entropy(u, equations::CompressibleEulerPotentialTemperatureEquations2D) rho, rho_v1, rho_v2, rho_theta = u k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma @@ -257,7 +257,7 @@ end return SVector(w1, w2, w3, w4) end -@inline function entropy_math(cons, equations::CompressiblePotentialEulerEquations2D) +@inline function entropy_math(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) # Mathematical entropy p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma @@ -267,19 +267,19 @@ end end # Default entropy is the mathematical entropy -@inline function entropy(cons, equations::CompressiblePotentialEulerEquations2D) +@inline function entropy(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) entropy_math(cons, equations) end -@inline function energy_total(cons, equations::CompressiblePotentialEulerEquations2D) +@inline function energy_total(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) entropy(cons, equations) end -@inline function energy_kinetic(cons, equations::CompressiblePotentialEulerEquations2D) +@inline function energy_kinetic(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) return 0.5 * (cons[2]^2 + cons[3]^2) / (cons[1]) end -@inline function max_abs_speeds(u, equations::CompressiblePotentialEulerEquations2D) +@inline function max_abs_speeds(u, equations::CompressibleEulerPotentialTemperatureEquations2D) rho, v1, v2, p = cons2prim(u, equations) c = sqrt(equations.gamma * p / rho) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a97393c..2bc2984 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -5,4 +5,4 @@ abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: abstract type AbstractCompressiblePotentialEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("compressible_moist_euler_2d_lucas.jl") -include("compressible_potential_euler_2d.jl") +include("compressible_euler_potential_temperature_2d.jl") From 12931b81a94fea1568a6118b9ac98236a60d5a02 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 14:39:27 +0200 Subject: [PATCH 06/27] renaming source term --- examples/elixir_potential_euler_dry_bubble.jl | 2 +- src/equations/compressible_euler_potential_temperature_2d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/elixir_potential_euler_dry_bubble.jl b/examples/elixir_potential_euler_dry_bubble.jl index 288a797..1e65276 100644 --- a/examples/elixir_potential_euler_dry_bubble.jl +++ b/examples/elixir_potential_euler_dry_bubble.jl @@ -71,7 +71,7 @@ mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, initial_condition = initial_condition_warm_bubble semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms = source_terms_bubble, + source_terms = source_terms_gravity, boundary_conditions = boundary_conditions) ############################################################################### diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index 8ddea9f..f318072 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -132,7 +132,7 @@ end return boundary_flux end -@inline function source_terms_bubble(u, x, t, +@inline function source_terms_gravity(u, x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) rho, _, _, _ = u return SVector(zero(eltype(u)), zero(eltype(u)), -equations.g * rho, From 5124bc8455773205e61dfba6978c959105f0ed93 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 14:46:46 +0200 Subject: [PATCH 07/27] renaming --- ...bubble.jl => elixir_euler_potential_temperature_dry_bubble.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/{elixir_potential_euler_dry_bubble.jl => elixir_euler_potential_temperature_dry_bubble.jl} (100%) diff --git a/examples/elixir_potential_euler_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl similarity index 100% rename from examples/elixir_potential_euler_dry_bubble.jl rename to examples/elixir_euler_potential_temperature_dry_bubble.jl From af1eb1f121d0502606c2d103aa2ecf601831bce2 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 15:05:48 +0200 Subject: [PATCH 08/27] renaming and test added --- ..._euler_potential_temperature_dry_bubble.jl | 3 +- ...ressible_euler_potential_temperature_2d.jl | 41 ++++++++++++------- src/equations/equations.jl | 2 +- test/runtests.jl | 4 ++ test/test_2d_euler_potential_temperature.jl | 37 +++++++++++++++++ 5 files changed, 70 insertions(+), 17 deletions(-) create mode 100644 test/test_2d_euler_potential_temperature.jl diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl index 1e65276..a885f87 100644 --- a/examples/elixir_euler_potential_temperature_dry_bubble.jl +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -1,8 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: flux_LMARS, source_terms_bubble, flux_theta -using Plots +using TrixiAtmo: flux_LMARS, source_terms_gravity, flux_theta function initial_condition_warm_bubble(x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index f318072..a7e2658 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -22,15 +22,19 @@ function CompressibleEulerPotentialTemperatureEquations2D(; g = 9.81, RealT = Fl R = c_p - c_v gamma = c_p / c_v a = 340.0 - return CompressibleEulerPotentialTemperatureEquations2D{RealT}(p_0, c_p, c_v, g, R, gamma, a) + return CompressibleEulerPotentialTemperatureEquations2D{RealT}(p_0, c_p, c_v, g, R, + gamma, a) end -function varnames(::typeof(cons2cons), ::CompressibleEulerPotentialTemperatureEquations2D) +function varnames(::typeof(cons2cons), + ::CompressibleEulerPotentialTemperatureEquations2D) ("rho", "rho_v1", "rho_v2", "rho_theta") end -varnames(::typeof(cons2prim), ::CompressibleEulerPotentialTemperatureEquations2D) = ("rho", "v1", - "v2", "p1") +varnames(::typeof(cons2prim), ::CompressibleEulerPotentialTemperatureEquations2D) = ("rho", + "v1", + "v2", + "p1") # Calculate 1D flux for a single point in the normal direction. # Note, this directional vector is not normalized. @@ -133,7 +137,7 @@ end end @inline function source_terms_gravity(u, x, t, - equations::CompressibleEulerPotentialTemperatureEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, _, _, _ = u return SVector(zero(eltype(u)), zero(eltype(u)), -equations.g * rho, zero(eltype(u))) @@ -224,7 +228,8 @@ end return SVector(f1, f2, f3, f4) end -@inline function prim2cons(prim, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function prim2cons(prim, + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, v1, v2, p = prim rho_v1 = rho * v1 rho_v2 = rho * v2 @@ -232,7 +237,8 @@ end return SVector(rho, rho_v1, rho_v2, rho_theta) end -@inline function cons2prim(u, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function cons2prim(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, rho_v1, rho_v2, rho_theta = u v1 = rho_v1 / rho v2 = rho_v2 / rho @@ -241,11 +247,13 @@ end return SVector(rho, v1, v2, p) end -@inline function cons2cons(u, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function cons2cons(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) return u end -@inline function cons2entropy(u, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function cons2entropy(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, rho_v1, rho_v2, rho_theta = u k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma @@ -257,7 +265,8 @@ end return SVector(w1, w2, w3, w4) end -@inline function entropy_math(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function entropy_math(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) # Mathematical entropy p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma @@ -267,19 +276,23 @@ end end # Default entropy is the mathematical entropy -@inline function entropy(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function entropy(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) entropy_math(cons, equations) end -@inline function energy_total(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function energy_total(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) entropy(cons, equations) end -@inline function energy_kinetic(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function energy_kinetic(cons, + equations::CompressibleEulerPotentialTemperatureEquations2D) return 0.5 * (cons[2]^2 + cons[3]^2) / (cons[1]) end -@inline function max_abs_speeds(u, equations::CompressibleEulerPotentialTemperatureEquations2D) +@inline function max_abs_speeds(u, + equations::CompressibleEulerPotentialTemperatureEquations2D) rho, v1, v2, p = cons2prim(u, equations) c = sqrt(equations.gamma * p / rho) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 2bc2984..f2350cc 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -2,7 +2,7 @@ using Trixi: AbstractEquations abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end -abstract type AbstractCompressiblePotentialEulerEquations{NDIMS, NVARS} <: +abstract type AbstractCompressibleEulerPotentialTemperatureEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("compressible_moist_euler_2d_lucas.jl") include("compressible_euler_potential_temperature_2d.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 541cfd8..9635f2e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,4 +20,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection" include("test_spherical_advection.jl") end + + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "euler_potential_temperature" + include("test_2d_euler_potential_temperature.jl") + end end diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl new file mode 100644 index 0000000..0b96e69 --- /dev/null +++ b/test/test_2d_euler_potential_temperature.jl @@ -0,0 +1,37 @@ +module TestExamples2DEulerPotentialTemperature + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") # TODO - This is a repetition from Trixi.jl + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory for examples? + +@trixiatmo_testset "elixir_euler_potential_temperature_dry_bubble" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_potential_temperature_dry_bubble.jl"), + l2=[ + 1.300427456987984e-6, + 2.6010873806861595e-5, + 0.0006660120005093007, + 9.51074191163579e-6, + ], + linf=[ + 1.031131432183141e-5, + 0.00020623855042956052, + 0.006392070867001616, + 7.841424786647622e-5, + ], + polydeg=3, + tspan=(0.0, 0.1)) + # 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 TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +end # module From e382140be0f9a5dcd4b9fd072c9b420f7793771d Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 15:40:49 +0200 Subject: [PATCH 09/27] added ec test --- .../elixir_euler_potential_temperature_ec.jl | 77 +++++++++++++++++++ src/TrixiAtmo.jl | 2 +- test/test_2d_euler_potential_temperature.jl | 29 ++++++- test/test_2d_moist_euler.jl | 28 +++---- test/test_spherical_advection.jl | 4 +- 5 files changed, 121 insertions(+), 19 deletions(-) create mode 100644 examples/elixir_euler_potential_temperature_ec.jl diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl new file mode 100644 index 0000000..eb27c5e --- /dev/null +++ b/examples/elixir_euler_potential_temperature_ec.jl @@ -0,0 +1,77 @@ +using OrdinaryDiffEq +using Trixi +using TrixiAtmo +using TrixiAtmo: flux_theta, prim2cons +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerPotentialTemperatureEquations2D() + +function initial_condition_weak_blast_wave(x, t, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) + # Set up polar coordinates + inicenter = SVector(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) + + # Calculate primitive variables + RealT = eltype(x) + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +volume_flux = flux_theta +solver = DGSEM(polydeg = 3, surface_flux = flux_theta, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) + +cells_per_dimension = (32, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, true)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.4) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +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/TrixiAtmo.jl b/src/TrixiAtmo.jl index cdce0a3..da2cc5b 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -26,7 +26,7 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, flux_LMARS, flux_theta +export flux_chandrashekar, flux_LMARS, flux_theta, prim2cons export examples_dir diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl index 0b96e69..69629af 100644 --- a/test/test_2d_euler_potential_temperature.jl +++ b/test/test_2d_euler_potential_temperature.jl @@ -14,13 +14,13 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 1.300427456987984e-6, 2.6010873806861595e-5, 0.0006660120005093007, - 9.51074191163579e-6, + 9.51074191163579e-6 ], linf=[ 1.031131432183141e-5, 0.00020623855042956052, 0.006392070867001616, - 7.841424786647622e-5, + 7.841424786647622e-5 ], polydeg=3, tspan=(0.0, 0.1)) @@ -34,4 +34,29 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory end end +@trixiatmo_testset "elixir_euler_potential_temperature_ec" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_potential_temperature_ec.jl"), + l2=[ + 0.06174365254988615, + 0.05018008812040643, + 0.05018774429827882, + 0.0057820937341387935 + ], + linf=[ + 0.2932352942992503, + 0.3108954213591686, + 0.31082193857647905, + 0.027465217019157606 + ]) + # 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 TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + end # module diff --git a/test/test_2d_moist_euler.jl b/test/test_2d_moist_euler.jl index d669e96..646ea0e 100644 --- a/test/test_2d_moist_euler.jl +++ b/test/test_2d_moist_euler.jl @@ -16,7 +16,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.0006660124630171347, 0.008969786054960861, 0.0, - 0.0, + 0.0 ], linf=[ 1.0312042909910168e-5, @@ -24,7 +24,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.006392107590872236, 0.07612038028310053, 0.0, - 0.0, + 0.0 ], polydeg=3, tspan=(0.0, 0.1)) @@ -46,7 +46,7 @@ end 7.938812668709457, 4500.747616411578, 0.00015592413050814787, - 0.00014163475049532796, + 0.00014163475049532796 ], linf=[ 0.1427479052298466, @@ -54,7 +54,7 @@ end 91.56822550162855, 49528.383866247605, 0.0019364397602254623, - 0.0013259689889851285, + 0.0013259689889851285 ], polydeg=3, cells_per_dimension=(16, 16), @@ -77,7 +77,7 @@ end 0.0006974588377288118, 1.715668353329522, 8.831269143134121e-7, - 1.025579538944668e-6, + 1.025579538944668e-6 ], linf=[ 8.055695643149896e-5, @@ -85,7 +85,7 @@ end 0.005897639251024437, 19.24776030163048, 1.0043133039065386e-5, - 1.1439046776775402e-5, + 1.1439046776775402e-5 ], polydeg=3, cells_per_dimension=(16, 8), @@ -109,7 +109,7 @@ end 0.01172830034500581, 9.898560584459009, 0.0, - 0.0, + 0.0 ], linf=[ 0.0017602202683439927, @@ -117,7 +117,7 @@ end 0.5945327351674782, 489.89171406268724, 0.0, - 0.0, + 0.0 ], polydeg=3, cells_per_dimension=(10, 8), @@ -140,7 +140,7 @@ end 2.609465609739115e-5, 6.323484379066432e-5, 0.0, - 0.0, + 0.0 ], linf=[ 7.549984224430872e-5, @@ -148,7 +148,7 @@ end 0.00015964938767742964, 0.0005425860570698049, 0.0, - 0.0, + 0.0 ], polydeg=3, cells_per_dimension=(10, 8), @@ -171,7 +171,7 @@ end 0.015104690877128945, 0.5242098451814421, 5.474006801215573e-10, - 1.1103883907226752e-10, + 1.1103883907226752e-10 ], linf=[ 0.00013219484616722177, @@ -179,7 +179,7 @@ end 0.03789645369775574, 3.90888311638264, 3.938382289041286e-9, - 6.892033377287209e-10, + 6.892033377287209e-10 ], polydeg=3, cells_per_dimension=(10, 8), @@ -203,7 +203,7 @@ end 0.07460345535073129, 5.943049264963717, 4.471792794168725e-9, - 7.10320253652373e-10, + 7.10320253652373e-10 ], linf=[ 0.0007084183215528839, @@ -211,7 +211,7 @@ end 0.3697160082709827, 27.843155286573165, 2.1168438904322837e-8, - 3.691699932047233e-9, + 3.691699932047233e-9 ], polydeg=3, cells_per_dimension=(10, 8), diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 31f7a29..c826091 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -15,14 +15,14 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") 8.23414117e-03, 1.84210648e-03, 0.00000000e+00, - 6.44302430e-02, + 6.44302430e-02 ], linf=[ 9.48950488e-02, 9.64811952e-02, 1.37453400e-02, 0.00000000e+00, - 4.09322999e-01, + 4.09322999e-01 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 374cb635af927e266d9f08f39a42b8d65fb53d1a Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 16:19:27 +0200 Subject: [PATCH 10/27] restyling and clean up --- ..._euler_potential_temperature_dry_bubble.jl | 2 +- examples/elixir_moist_euler_dry_bubble.jl | 2 +- examples/elixir_moist_euler_moist_bubble.jl | 2 +- ...oist_euler_nonhydrostatic_gravity_waves.jl | 2 +- ...ressible_euler_potential_temperature_2d.jl | 42 +++++++++---------- .../compressible_moist_euler_2d_lucas.jl | 7 ++-- src/equations/equations.jl | 7 ++++ 7 files changed, 35 insertions(+), 29 deletions(-) diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl index a885f87..55167c6 100644 --- a/examples/elixir_euler_potential_temperature_dry_bubble.jl +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -54,7 +54,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 3 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS +surface_flux = flux_LMARS(340.0) volume_flux = flux_theta volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index d34c63d..30ba16e 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -61,7 +61,7 @@ source_term = source_terms_geopotential polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS +surface_flux = flux_LMARS(340.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index bdf7b4e..9fd2bb3 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -242,7 +242,7 @@ source_term = source_terms_moist_bubble polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS +surface_flux = flux_LMARS(340.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index 02ab73e..82ea26e 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -55,7 +55,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS +surface_flux = flux_LMARS(340.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index a7e2658..b7f6cee 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -12,7 +12,6 @@ struct CompressibleEulerPotentialTemperatureEquations2D{RealT <: Real} <: g::RealT R::RealT gamma::RealT - a::RealT end function CompressibleEulerPotentialTemperatureEquations2D(; g = 9.81, RealT = Float64) @@ -21,9 +20,8 @@ function CompressibleEulerPotentialTemperatureEquations2D(; g = 9.81, RealT = Fl c_v = 717.0 R = c_p - c_v gamma = c_p / c_v - a = 340.0 return CompressibleEulerPotentialTemperatureEquations2D{RealT}(p_0, c_p, c_v, g, R, - gamma, a) + gamma) end function varnames(::typeof(cons2cons), @@ -99,13 +97,13 @@ end if v_normal <= 0.0 sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * + (1.0 + 0.5f0 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * inv(gamma - 1)) else # v_normal > 0.0 A = 2.0 / ((gamma + 1) * rho_local) B = p_local * (gamma - 1) / (gamma + 1) p_star = p_local + - 0.5 * v_normal / A * + 0.5f0 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4.0 * A * (p_local + B))) end @@ -169,9 +167,11 @@ end # Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. -@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, + + +@inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquations2D) - @unpack a = equations + a = 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) @@ -179,13 +179,12 @@ end 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 necessary interface flux components norm_ = norm(normal_direction) - rho = 0.5 * (rho_ll + rho_rr) + rho = 0.5f0 * (rho_ll + rho_rr) - p_interface = 0.5 * (p_ll + p_rr) - 0.5 * a * rho * (v_rr - v_ll) / norm_ - v_interface = 0.5 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + p_interface = 0.5f0 * (p_ll + p_rr) - 0.5f0 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5f0 * (v_ll + v_rr) - 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ if (v_interface > 0) f1, f2, f3, f4 = u_ll * v_interface @@ -200,6 +199,7 @@ end end ## Entropy (total energy) conservative flux for the Compressible Euler with the Potential Formulation + @inline function flux_theta(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquations2D) # Unpack left and right state @@ -211,20 +211,18 @@ end _, _, _, rho_theta_rr = u_rr # Compute the necessary mean values rho_mean = ln_mean(rho_ll, rho_rr) + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) - # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` - # in exact arithmetic since - # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) - # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] - f4 = gammamean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f4 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) return SVector(f1, f2, f3, f4) end @@ -257,7 +255,7 @@ end rho, rho_v1, rho_v2, rho_theta = u k = equations.p_0 * (equations.R / equations.p_0)^equations.gamma - w1 = -0.5 * rho_v1^2 / (rho)^2 - 0.5 * rho_v2^2 / (rho)^2 + w1 = -0.5f0 * rho_v1^2 / (rho)^2 - 0.5f0 * rho_v2^2 / (rho)^2 w2 = rho_v1 / rho w3 = rho_v2 / rho w4 = equations.gamma / (equations.gamma - 1) * k * (rho_theta)^(equations.gamma - 1) @@ -288,7 +286,7 @@ end @inline function energy_kinetic(cons, equations::CompressibleEulerPotentialTemperatureEquations2D) - return 0.5 * (cons[2]^2 + cons[3]^2) / (cons[1]) + return 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1]) end @inline function max_abs_speeds(u, diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 638cc39..18b85e8 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -23,7 +23,6 @@ struct CompressibleMoistEulerEquations2D{RealT <: Real} <: kappa::RealT # ratio of the gas constant R_d gamma::RealT # = inv(kappa- 1); can be used to write slow divisions as fast multiplications L_00::RealT # latent heat of evaporation at 0 K - a::RealT end function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64) @@ -452,9 +451,11 @@ end # Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. -@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, + +@inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleMoistEulerEquations2D) - @unpack a = equations + + a = flux_lmars.speed_of_sound # Unpack left and right state rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, rho_qv_ll, rho_ql_ll = u_ll rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr, rho_qv_rr, rho_ql_rr = u_rr diff --git a/src/equations/equations.jl b/src/equations/equations.jl index f2350cc..201fe53 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -4,5 +4,12 @@ abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end abstract type AbstractCompressibleEulerPotentialTemperatureEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end + + + struct flux_LMARS{SpeedOfSound} + # Estimate for the speed of sound + speed_of_sound::SpeedOfSound + end + include("compressible_moist_euler_2d_lucas.jl") include("compressible_euler_potential_temperature_2d.jl") From b75cfef9f5563fae10830b07f88d77e97b7b6931 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 16:22:12 +0200 Subject: [PATCH 11/27] formatting --- .../compressible_euler_potential_temperature_2d.jl | 5 ++--- src/equations/compressible_moist_euler_2d_lucas.jl | 3 +-- src/equations/equations.jl | 9 ++++----- 3 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index b7f6cee..21ed1d0 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -98,7 +98,7 @@ end sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed p_star = p_local * (1.0 + 0.5f0 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * - inv(gamma - 1)) + inv(gamma - 1)) else # v_normal > 0.0 A = 2.0 / ((gamma + 1) * rho_local) B = p_local * (gamma - 1) / (gamma + 1) @@ -168,9 +168,8 @@ end # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. - @inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerPotentialTemperatureEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) a = flux_lmars.speed_of_sound # Unpack left and right state rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 18b85e8..1bd83df 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -453,8 +453,7 @@ end # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. @inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleMoistEulerEquations2D) - + equations::CompressibleMoistEulerEquations2D) a = flux_lmars.speed_of_sound # Unpack left and right state rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, rho_qv_ll, rho_ql_ll = u_ll diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 201fe53..016d5a9 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -5,11 +5,10 @@ abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: abstract type AbstractCompressibleEulerPotentialTemperatureEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +struct flux_LMARS{SpeedOfSound} + # Estimate for the speed of sound + speed_of_sound::SpeedOfSound +end - struct flux_LMARS{SpeedOfSound} - # Estimate for the speed of sound - speed_of_sound::SpeedOfSound - end - include("compressible_moist_euler_2d_lucas.jl") include("compressible_euler_potential_temperature_2d.jl") From a229219ee1bf55bde7f7247d48a053d1ed128541 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 16:36:17 +0200 Subject: [PATCH 12/27] minor changes --- examples/elixir_moist_euler_dry_bubble.jl | 2 +- examples/elixir_moist_euler_moist_bubble.jl | 2 +- ...oist_euler_nonhydrostatic_gravity_waves.jl | 2 +- test/test_2d_euler_potential_temperature.jl | 8 +++--- test/test_2d_moist_euler.jl | 28 +++++++++---------- test/test_spherical_advection.jl | 4 +-- 6 files changed, 23 insertions(+), 23 deletions(-) diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index 30ba16e..d7ecaa9 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -61,7 +61,7 @@ source_term = source_terms_geopotential polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS(340.0) +surface_flux = flux_LMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 9fd2bb3..209ece8 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -242,7 +242,7 @@ source_term = source_terms_moist_bubble polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS(340.0) +surface_flux = flux_LMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index 82ea26e..d384599 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -55,7 +55,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS(340.0) +surface_flux = flux_LMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl index 69629af..737e6e5 100644 --- a/test/test_2d_euler_potential_temperature.jl +++ b/test/test_2d_euler_potential_temperature.jl @@ -14,13 +14,13 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 1.300427456987984e-6, 2.6010873806861595e-5, 0.0006660120005093007, - 9.51074191163579e-6 + 9.51074191163579e-6, ], linf=[ 1.031131432183141e-5, 0.00020623855042956052, 0.006392070867001616, - 7.841424786647622e-5 + 7.841424786647622e-5, ], polydeg=3, tspan=(0.0, 0.1)) @@ -41,13 +41,13 @@ end 0.06174365254988615, 0.05018008812040643, 0.05018774429827882, - 0.0057820937341387935 + 0.0057820937341387935, ], linf=[ 0.2932352942992503, 0.3108954213591686, 0.31082193857647905, - 0.027465217019157606 + 0.027465217019157606, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_2d_moist_euler.jl b/test/test_2d_moist_euler.jl index 646ea0e..d669e96 100644 --- a/test/test_2d_moist_euler.jl +++ b/test/test_2d_moist_euler.jl @@ -16,7 +16,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.0006660124630171347, 0.008969786054960861, 0.0, - 0.0 + 0.0, ], linf=[ 1.0312042909910168e-5, @@ -24,7 +24,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.006392107590872236, 0.07612038028310053, 0.0, - 0.0 + 0.0, ], polydeg=3, tspan=(0.0, 0.1)) @@ -46,7 +46,7 @@ end 7.938812668709457, 4500.747616411578, 0.00015592413050814787, - 0.00014163475049532796 + 0.00014163475049532796, ], linf=[ 0.1427479052298466, @@ -54,7 +54,7 @@ end 91.56822550162855, 49528.383866247605, 0.0019364397602254623, - 0.0013259689889851285 + 0.0013259689889851285, ], polydeg=3, cells_per_dimension=(16, 16), @@ -77,7 +77,7 @@ end 0.0006974588377288118, 1.715668353329522, 8.831269143134121e-7, - 1.025579538944668e-6 + 1.025579538944668e-6, ], linf=[ 8.055695643149896e-5, @@ -85,7 +85,7 @@ end 0.005897639251024437, 19.24776030163048, 1.0043133039065386e-5, - 1.1439046776775402e-5 + 1.1439046776775402e-5, ], polydeg=3, cells_per_dimension=(16, 8), @@ -109,7 +109,7 @@ end 0.01172830034500581, 9.898560584459009, 0.0, - 0.0 + 0.0, ], linf=[ 0.0017602202683439927, @@ -117,7 +117,7 @@ end 0.5945327351674782, 489.89171406268724, 0.0, - 0.0 + 0.0, ], polydeg=3, cells_per_dimension=(10, 8), @@ -140,7 +140,7 @@ end 2.609465609739115e-5, 6.323484379066432e-5, 0.0, - 0.0 + 0.0, ], linf=[ 7.549984224430872e-5, @@ -148,7 +148,7 @@ end 0.00015964938767742964, 0.0005425860570698049, 0.0, - 0.0 + 0.0, ], polydeg=3, cells_per_dimension=(10, 8), @@ -171,7 +171,7 @@ end 0.015104690877128945, 0.5242098451814421, 5.474006801215573e-10, - 1.1103883907226752e-10 + 1.1103883907226752e-10, ], linf=[ 0.00013219484616722177, @@ -179,7 +179,7 @@ end 0.03789645369775574, 3.90888311638264, 3.938382289041286e-9, - 6.892033377287209e-10 + 6.892033377287209e-10, ], polydeg=3, cells_per_dimension=(10, 8), @@ -203,7 +203,7 @@ end 0.07460345535073129, 5.943049264963717, 4.471792794168725e-9, - 7.10320253652373e-10 + 7.10320253652373e-10, ], linf=[ 0.0007084183215528839, @@ -211,7 +211,7 @@ end 0.3697160082709827, 27.843155286573165, 2.1168438904322837e-8, - 3.691699932047233e-9 + 3.691699932047233e-9, ], polydeg=3, cells_per_dimension=(10, 8), diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index c826091..31f7a29 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -15,14 +15,14 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") 8.23414117e-03, 1.84210648e-03, 0.00000000e+00, - 6.44302430e-02 + 6.44302430e-02, ], linf=[ 9.48950488e-02, 9.64811952e-02, 1.37453400e-02, 0.00000000e+00, - 4.09322999e-01 + 4.09322999e-01, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From bf6566936488732af43e9e0b1b718297956c88e7 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 16:42:58 +0200 Subject: [PATCH 13/27] minor changes --- src/equations/compressible_moist_euler_2d_lucas.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 1bd83df..7a440d2 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -37,9 +37,8 @@ function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64) gamma = c_pd / c_vd # = 1/(1 - kappa) kappa = 1 - inv(gamma) L_00 = 3147620.0 - a = 360.0 return CompressibleMoistEulerEquations2D{RealT}(p_0, c_pd, c_vd, R_d, c_pv, c_vv, - R_v, c_pl, g, kappa, gamma, L_00, a) + R_v, c_pl, g, kappa, gamma, L_00) end # Calculate 1D flux for a single point. From 6b3da0b7fb4db7cfdfd853117b8b32936904f0ca Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 17:04:14 +0200 Subject: [PATCH 14/27] minor changes --- examples/elixir_euler_potential_temperature_ec.jl | 4 ++-- src/TrixiAtmo.jl | 2 +- test/test_trixi_consistency.jl | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl index eb27c5e..2570783 100644 --- a/examples/elixir_euler_potential_temperature_ec.jl +++ b/examples/elixir_euler_potential_temperature_ec.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: flux_theta, prim2cons +using TrixiAtmo: flux_theta ############################################################################### # semidiscretization of the compressible Euler equations equations = CompressibleEulerPotentialTemperatureEquations2D() @@ -24,7 +24,7 @@ function initial_condition_weak_blast_wave(x, t, v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - return prim2cons(SVector(rho, v1, v2, p), equations) + return TrixiAtmo.prim2cons(SVector(rho, v1, v2, p), equations) end initial_condition = initial_condition_weak_blast_wave diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index da2cc5b..cdce0a3 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -26,7 +26,7 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, flux_LMARS, flux_theta, prim2cons +export flux_chandrashekar, flux_LMARS, flux_theta export examples_dir diff --git a/test/test_trixi_consistency.jl b/test/test_trixi_consistency.jl index b13b5d3..72e047f 100644 --- a/test/test_trixi_consistency.jl +++ b/test/test_trixi_consistency.jl @@ -49,7 +49,7 @@ isdir(outdir) && rm(outdir, recursive = true) trixi_include(trixi_elixir, equations = equations_moist, volume_flux = flux_chandrashekar, - surface_flux = flux_LMARS, + surface_flux = flux_LMARS(360.0), maxiters = maxiters) errors_atmo = Main.analysis_callback(Main.sol) From e952915df9117e48c069aab5de6389f233d0efff Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Tue, 24 Sep 2024 21:54:43 +0200 Subject: [PATCH 15/27] renaming, enabling Float32 --- ..._euler_potential_temperature_dry_bubble.jl | 4 +- examples/elixir_moist_euler_dry_bubble.jl | 4 +- examples/elixir_moist_euler_moist_bubble.jl | 4 +- ...oist_euler_nonhydrostatic_gravity_waves.jl | 4 +- src/TrixiAtmo.jl | 2 +- ...ressible_euler_potential_temperature_2d.jl | 4 +- .../compressible_moist_euler_2d_lucas.jl | 50 +++++++++---------- src/equations/equations.jl | 8 +-- test/test_trixi_consistency.jl | 2 +- 9 files changed, 41 insertions(+), 41 deletions(-) diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl index 55167c6..0ec083a 100644 --- a/examples/elixir_euler_potential_temperature_dry_bubble.jl +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: flux_LMARS, source_terms_gravity, flux_theta +using TrixiAtmo: FluxLMARS, source_terms_gravity, flux_theta function initial_condition_warm_bubble(x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) @@ -54,7 +54,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 3 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS(340.0) +surface_flux = FluxLMARS(340.0) volume_flux = flux_theta volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index d7ecaa9..1256a9c 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: flux_LMARS, source_terms_geopotential, cons2drypot +using TrixiAtmo: FluxLMARS, source_terms_geopotential, cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equations @@ -61,7 +61,7 @@ source_term = source_terms_geopotential polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS(360.0) +surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 209ece8..956c4a2 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble, - flux_LMARS + FluxLMARS using NLsolve: nlsolve ############################################################################### @@ -242,7 +242,7 @@ source_term = source_terms_moist_bubble polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS(360.0) +surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index d384599..dd8de4f 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -3,7 +3,7 @@ using Trixi, TrixiAtmo using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential, source_terms_phase_change, source_terms_nonhydrostatic_rayleigh_sponge, - cons2drypot, flux_LMARS + cons2drypot, FluxLMARS ############################################################################### # semidiscretization of the compressible moist Euler equation @@ -55,7 +55,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic, polydeg = 4 basis = LobattoLegendreBasis(polydeg) -surface_flux = flux_LMARS(360.0) +surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index cdce0a3..e4f64c3 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -26,7 +26,7 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, flux_LMARS, flux_theta +export flux_chandrashekar, FluxLMARS, flux_theta export examples_dir diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index 21ed1d0..6aa3ad8 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -1,6 +1,6 @@ using Trixi using Trixi: ln_mean, stolarsky_mean -import Trixi: varnames, cons2cons, cons2prim, cons2entropy, entropy +import Trixi: varnames, cons2cons, cons2prim, cons2entropy, entropy, FluxLMARS @muladd begin #! format: noindent @@ -168,7 +168,7 @@ end # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. -@inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector, +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquations2D) a = flux_lmars.speed_of_sound # Unpack left and right state diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 7a440d2..4bf1676 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -5,7 +5,7 @@ using Trixi using Trixi: ln_mean, inv_ln_mean import Trixi: varnames, flux_chandrashekar, boundary_condition_slip_wall, cons2prim, cons2entropy, max_abs_speed_naive, max_abs_speeds, - entropy, energy_total, flux + entropy, energy_total, flux, FluxLMARS @muladd begin #! format: noindent @@ -118,13 +118,13 @@ end if v_normal <= 0.0 sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed p_star = p_local * - (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * + (1.0 + 0.5f0 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * inv(gamma - 1)) else # v_normal > 0.0 A = 2.0 / ((gamma + 1) * rho_local) B = p_local * (gamma - 1) / (gamma + 1) p_star = p_local + - 0.5 * v_normal / A * + 0.5f0 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4.0 * A * (p_local + B))) end @@ -368,11 +368,11 @@ end # positive even power with default value 2 gamma = 2.0 # relaxation coefficient > 0 - alpha = 0.5 + alpha = 0.5f0 tau_s = zero(eltype(u)) if z > z_s - tau_s = alpha * sin(0.5 * (z - z_s) * inv(z_top - z_s))^(gamma) + tau_s = alpha * sin(0.5f0 * (z - z_s) * inv(z_top - z_s))^(gamma) end return SVector(zero(eltype(u)), @@ -398,7 +398,7 @@ end v2 = rho_v2 / rho # Inner energy - rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + rho_e = (rho_E - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) # Absolute temperature T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) @@ -451,7 +451,7 @@ end # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. -@inline function (flux_lmars::flux_LMARS)(u_ll, u_rr, normal_direction::AbstractVector, +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleMoistEulerEquations2D) a = flux_lmars.speed_of_sound # Unpack left and right state @@ -473,9 +473,9 @@ end # Compute the necessary interface flux components norm_ = norm(normal_direction) - rho = 0.5 * (rho_ll + rho_rr) - p_interface = 0.5 * (p_ll + p_rr) - beta * 0.5 * a * rho * (v_rr - v_ll) / norm_ - v_interface = 0.5 * (v_ll + v_rr) - beta * 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + rho = 0.5f0 * (rho_ll + rho_rr) + p_interface = 0.5f0 * (p_ll + p_rr) - beta * 0.5f0 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5f0 * (v_ll + v_rr) - beta * 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ if (v_interface > 0) f1, f2, f3, f4, f5, f6 = u_ll * v_interface @@ -548,7 +548,7 @@ end g_v = L_00 + (c_pv - s_v) * T g_l = (c_pl - s_l) * T - w1 = g_d - 0.5 * v_square + w1 = g_d - 0.5f0 * v_square w2 = v1 w3 = v2 w4 = -1 @@ -565,7 +565,7 @@ end rho_v2 = rho * v2 rho_qv = rho * qv rho_ql = rho * ql - rho_E = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + rho_E = p * equations.inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2) return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) end @@ -669,7 +669,7 @@ end v2 = rho_v2 / rho # inner energy - rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + rho_e = (rho_E - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) # Absolute Temperature T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) @@ -715,7 +715,7 @@ end v2 = rho_v2 / rho # inner energy - rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + rho_e = (rho_E - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)) # Absolute Temperature T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) @@ -949,8 +949,8 @@ end v2_rr = rho_v2_rr / rho_rr # inner energy - rho_e_ll = (rho_E_ll - 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll)) - rho_e_rr = (rho_E_rr - 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr)) + rho_e_ll = (rho_E_ll - 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll)) + rho_e_rr = (rho_E_rr - 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr)) # Absolute Temperature T_ll = (rho_e_ll - L_00 * rho_qv_ll) / @@ -976,18 +976,18 @@ end inv_T_mean = inv_ln_mean(inv(T_ll), inv(T_rr)) end - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v1_square_avg = 0.5 * (v1_ll^2 + v1_rr^2) - v2_square_avg = 0.5 * (v2_ll^2 + v2_rr^2) - rho_qd_avg = 0.5 * (rho_qd_ll + rho_qd_rr) - rho_qv_avg = 0.5 * (rho_qv_ll + rho_qv_rr) - rho_ql_avg = 0.5 * (rho_ql_ll + rho_ql_rr) - inv_T_avg = 0.5 * (inv(T_ll) + inv(T_rr)) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v1_square_avg = 0.5f0 * (v1_ll^2 + v1_rr^2) + v2_square_avg = 0.5f0 * (v2_ll^2 + v2_rr^2) + rho_qd_avg = 0.5f0 * (rho_qd_ll + rho_qd_rr) + rho_qv_avg = 0.5f0 * (rho_qv_ll + rho_qv_rr) + rho_ql_avg = 0.5f0 * (rho_ql_ll + rho_ql_rr) + inv_T_avg = 0.5f0 * (inv(T_ll) + inv(T_rr)) v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg p_int = inv(inv_T_avg) * (R_d * rho_qd_avg + R_v * rho_qv_avg + R_q * rho_ql_avg) - K_avg = 0.5 * (v1_square_avg + v2_square_avg) + K_avg = 0.5f0 * (v1_square_avg + v2_square_avg) f_1d = rho_qd_mean * v_dot_n_avg f_1v = rho_qv_mean * v_dot_n_avg diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 016d5a9..3df151f 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -5,10 +5,10 @@ abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: abstract type AbstractCompressibleEulerPotentialTemperatureEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end -struct flux_LMARS{SpeedOfSound} - # Estimate for the speed of sound - speed_of_sound::SpeedOfSound -end +# struct flux_LMARS{SpeedOfSound} +# # Estimate for the speed of sound +# speed_of_sound::SpeedOfSound +# end include("compressible_moist_euler_2d_lucas.jl") include("compressible_euler_potential_temperature_2d.jl") diff --git a/test/test_trixi_consistency.jl b/test/test_trixi_consistency.jl index 72e047f..3596d9c 100644 --- a/test/test_trixi_consistency.jl +++ b/test/test_trixi_consistency.jl @@ -49,7 +49,7 @@ isdir(outdir) && rm(outdir, recursive = true) trixi_include(trixi_elixir, equations = equations_moist, volume_flux = flux_chandrashekar, - surface_flux = flux_LMARS(360.0), + surface_flux = FluxLMARS(360.0), maxiters = maxiters) errors_atmo = Main.analysis_callback(Main.sol) From 78d4bcc0a3933f9acbdd5b12125eeccc1448e496 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Tue, 24 Sep 2024 21:57:53 +0200 Subject: [PATCH 16/27] formatting --- .../compressible_euler_potential_temperature_2d.jl | 2 +- src/equations/compressible_moist_euler_2d_lucas.jl | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index 6aa3ad8..1b9aee8 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -169,7 +169,7 @@ end # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerPotentialTemperatureEquations2D) + equations::CompressibleEulerPotentialTemperatureEquations2D) a = flux_lmars.speed_of_sound # Unpack left and right state rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 4bf1676..53557af 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -119,7 +119,7 @@ end sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed p_star = p_local * (1.0 + 0.5f0 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * - inv(gamma - 1)) + inv(gamma - 1)) else # v_normal > 0.0 A = 2.0 / ((gamma + 1) * rho_local) B = p_local * (gamma - 1) / (gamma + 1) @@ -452,7 +452,7 @@ end # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleMoistEulerEquations2D) + equations::CompressibleMoistEulerEquations2D) a = flux_lmars.speed_of_sound # Unpack left and right state rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, rho_qv_ll, rho_ql_ll = u_ll @@ -475,7 +475,8 @@ end rho = 0.5f0 * (rho_ll + rho_rr) p_interface = 0.5f0 * (p_ll + p_rr) - beta * 0.5f0 * a * rho * (v_rr - v_ll) / norm_ - v_interface = 0.5f0 * (v_ll + v_rr) - beta * 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + v_interface = 0.5f0 * (v_ll + v_rr) - + beta * 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ if (v_interface > 0) f1, f2, f3, f4, f5, f6 = u_ll * v_interface From a820c9781f8584dfecaaa9d146544ee31b92e6cf Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Tue, 24 Sep 2024 21:59:32 +0200 Subject: [PATCH 17/27] deleted useless parts --- src/equations/equations.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 3df151f..6fe1256 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -5,10 +5,5 @@ abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: abstract type AbstractCompressibleEulerPotentialTemperatureEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end -# struct flux_LMARS{SpeedOfSound} -# # Estimate for the speed of sound -# speed_of_sound::SpeedOfSound -# end - include("compressible_moist_euler_2d_lucas.jl") include("compressible_euler_potential_temperature_2d.jl") From b2e1b90c9c64cdb26226c318a51909202ff7f424 Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 11:40:37 +0200 Subject: [PATCH 18/27] Update examples/elixir_moist_euler_dry_bubble.jl Co-authored-by: Hendrik Ranocha --- examples/elixir_moist_euler_dry_bubble.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index 1256a9c..8571896 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: FluxLMARS, source_terms_geopotential, cons2drypot +using TrixiAtmo: source_terms_geopotential, cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equations From 9021e28eea68da6fd7a17b148f9fa36ac107cfdd Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 11:40:42 +0200 Subject: [PATCH 19/27] Update examples/elixir_moist_euler_moist_bubble.jl Co-authored-by: Hendrik Ranocha --- examples/elixir_moist_euler_moist_bubble.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 956c4a2..4d18f4c 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble, - FluxLMARS using NLsolve: nlsolve ############################################################################### From c53a3cdacf6eea61c3bd9d48f5668ceab13ea193 Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 11:40:47 +0200 Subject: [PATCH 20/27] Update examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl Co-authored-by: Hendrik Ranocha --- examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index dd8de4f..4646a8b 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -3,7 +3,7 @@ using Trixi, TrixiAtmo using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential, source_terms_phase_change, source_terms_nonhydrostatic_rayleigh_sponge, - cons2drypot, FluxLMARS + cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equation From 3cf8db19daffdeb54b09a4efe9dedaed7482c6bb Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 11:40:52 +0200 Subject: [PATCH 21/27] Update src/TrixiAtmo.jl Co-authored-by: Hendrik Ranocha --- src/TrixiAtmo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index e4f64c3..f466f99 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -26,7 +26,7 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, FluxLMARS, flux_theta +export flux_chandrashekar, flux_theta export examples_dir From d2a23b116201cdb800d9c364ae7a1f088ca88867 Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 11:41:01 +0200 Subject: [PATCH 22/27] Update examples/elixir_moist_euler_moist_bubble.jl Co-authored-by: Hendrik Ranocha --- examples/elixir_moist_euler_moist_bubble.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 4d18f4c..d406bab 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -1,6 +1,6 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo -using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble, +using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble using NLsolve: nlsolve ############################################################################### From dbbd06a4c9542abbc9f16d65feaff5a82178f8c6 Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 11:42:37 +0200 Subject: [PATCH 23/27] Update elixir_euler_potential_temperature_dry_bubble.jl --- examples/elixir_euler_potential_temperature_dry_bubble.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl index 0ec083a..84c2f2b 100644 --- a/examples/elixir_euler_potential_temperature_dry_bubble.jl +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: FluxLMARS, source_terms_gravity, flux_theta +using TrixiAtmo: source_terms_gravity, flux_theta function initial_condition_warm_bubble(x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) From 41ec9bf0597b786974b47e238290e88fcbe6b655 Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 11:56:47 +0200 Subject: [PATCH 24/27] Update test_trixi_consistency.jl --- test/test_trixi_consistency.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_trixi_consistency.jl b/test/test_trixi_consistency.jl index 3596d9c..30981d2 100644 --- a/test/test_trixi_consistency.jl +++ b/test/test_trixi_consistency.jl @@ -1,5 +1,7 @@ module TestTrixiConsistency +using Trixi + include("test_trixiatmo.jl") EXAMPLES_DIR = TrixiAtmo.examples_dir() From de71e98c9271fcb13eac7916cc82a6db4aafba2f Mon Sep 17 00:00:00 2001 From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> Date: Wed, 25 Sep 2024 12:52:10 +0200 Subject: [PATCH 25/27] Update test_trixi_consistency.jl --- test/test_trixi_consistency.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/test_trixi_consistency.jl b/test/test_trixi_consistency.jl index 30981d2..477a3b8 100644 --- a/test/test_trixi_consistency.jl +++ b/test/test_trixi_consistency.jl @@ -1,7 +1,5 @@ module TestTrixiConsistency -using Trixi - include("test_trixiatmo.jl") EXAMPLES_DIR = TrixiAtmo.examples_dir() @@ -51,7 +49,7 @@ isdir(outdir) && rm(outdir, recursive = true) trixi_include(trixi_elixir, equations = equations_moist, volume_flux = flux_chandrashekar, - surface_flux = FluxLMARS(360.0), + surface_flux = Trixi.FluxLMARS(360.0), maxiters = maxiters) errors_atmo = Main.analysis_callback(Main.sol) From 14c3820ead1325d89aad73a827bb329800fce6ce Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 25 Sep 2024 16:29:18 +0200 Subject: [PATCH 26/27] gravity wave test, testing ec fluxes --- ..._euler_potential_temperature_dry_bubble.jl | 1 - .../elixir_euler_potential_temperature_ec.jl | 2 +- ...uler_potential_temperature_gravity_wave.jl | 102 ++++++++++++++++++ examples/elixir_moist_euler_dry_bubble.jl | 2 +- examples/elixir_moist_euler_moist_bubble.jl | 3 +- ...oist_euler_nonhydrostatic_gravity_waves.jl | 2 +- src/TrixiAtmo.jl | 3 +- ...ressible_euler_potential_temperature_2d.jl | 84 ++++++++++++++- test/test_2d_euler_potential_temperature.jl | 27 +++++ 9 files changed, 218 insertions(+), 8 deletions(-) create mode 100644 examples/elixir_euler_potential_temperature_gravity_wave.jl diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl index 0ec083a..1491ee8 100644 --- a/examples/elixir_euler_potential_temperature_dry_bubble.jl +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: FluxLMARS, source_terms_gravity, flux_theta function initial_condition_warm_bubble(x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl index 2570783..c3faf15 100644 --- a/examples/elixir_euler_potential_temperature_ec.jl +++ b/examples/elixir_euler_potential_temperature_ec.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: flux_theta + ############################################################################### # semidiscretization of the compressible Euler equations equations = CompressibleEulerPotentialTemperatureEquations2D() diff --git a/examples/elixir_euler_potential_temperature_gravity_wave.jl b/examples/elixir_euler_potential_temperature_gravity_wave.jl new file mode 100644 index 0000000..64900ef --- /dev/null +++ b/examples/elixir_euler_potential_temperature_gravity_wave.jl @@ -0,0 +1,102 @@ +using OrdinaryDiffEq +using Trixi +using TrixiAtmo + +function initial_condition_gravity_wave(x, t, + equations::CompressibleEulerPotentialTemperatureEquations2D) + g = equations.g + c_p = equations.c_p + c_v = equations.c_v + # center of perturbation + x_c = 100_000.0 + a = 5_000 + H = 10_000 + # perturbation in potential temperature + potential_temperature_ref = 300.0 * exp(0.01^2 / g * x[2]) + potential_temperature_perturbation = 0.01 * sinpi(x[2] / H) / (1 + (x[1] - x_c)^2 / a^2) + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 + g^2 / (c_p * 300.0 * 0.01^2) * (exp(-0.01^2 / g * x[2]) - 1) + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + T = potential_temperature * exner + + # density + rho = p / (R * T) + v1 = 20.0 + v2 = 0.0 + + return SVector(rho, rho * v1, rho * v2, rho * potential_temperature) +end + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerPotentialTemperatureEquations2D() + +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 = 3 +basis = LobattoLegendreBasis(polydeg) + +surface_flux = FluxLMARS(340.0) + +solver = DGSEM(basis, surface_flux) + +coordinates_min = (0.0, 0.0) +coordinates_max = (300_000.0, 10_000.0) + +cells_per_dimension = (600, 20) # Delta x = Delta z = 1 km +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) +initial_condition = initial_condition_gravity_wave + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_gravity, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3000.0) # 1000 seconds final time + +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 = 1.0) + +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 = 1e-1, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index 1256a9c..8571896 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo -using TrixiAtmo: FluxLMARS, source_terms_geopotential, cons2drypot +using TrixiAtmo: source_terms_geopotential, cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equations diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 956c4a2..d406bab 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo -using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble, - FluxLMARS +using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble using NLsolve: nlsolve ############################################################################### diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index dd8de4f..4646a8b 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -3,7 +3,7 @@ using Trixi, TrixiAtmo using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential, source_terms_phase_change, source_terms_nonhydrostatic_rayleigh_sponge, - cons2drypot, FluxLMARS + cons2drypot ############################################################################### # semidiscretization of the compressible moist Euler equation diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index e4f64c3..fcfb4c0 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -26,7 +26,8 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, FluxLMARS, flux_theta +export flux_chandrashekar, flux_theta, flux_theta_rhoAM, flux_theta_physentropy, + flux_theta_physentropy2 export examples_dir diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl index 1b9aee8..926d7bc 100644 --- a/src/equations/compressible_euler_potential_temperature_2d.jl +++ b/src/equations/compressible_euler_potential_temperature_2d.jl @@ -225,6 +225,88 @@ end return SVector(f1, f2, f3, f4) end +@inline function flux_theta_rhoAM(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # 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_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_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + rho_mean = 0.5f0 * (rho_ll + rho_rr) + + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + return SVector(f1, f2, f3, f4) +end + +@inline function flux_theta_physentropy(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # 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_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_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + #rho_mean = ln_mean(rho_ll, rho_rr) + #rho_mean = 0.5f0 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll / rho_theta_ll, rho_rr / rho_theta_rr) + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + #f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f4 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f1 = f4 * rho_mean + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + + return SVector(f1, f2, f3, f4) +end + +@inline function flux_theta_physentropy2(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # 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_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_theta_ll = u_ll + _, _, _, rho_theta_rr = u_rr + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + #rho_mean = 0.5f0 * (rho_ll + rho_rr) + #rho_mean = ln_mean(rho_ll/rho_theta_ll, rho_rr/rho_theta_rr) + gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = inv_ln_mean(rho_ll / rho_theta_ll, rho_rr / rho_theta_rr) * f1 + return SVector(f1, f2, f3, f4) +end + @inline function prim2cons(prim, equations::CompressibleEulerPotentialTemperatureEquations2D) rho, v1, v2, p = prim @@ -267,7 +349,7 @@ end # Mathematical entropy p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma - U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2 + cons[3]^2) / (cons[1])) + U = (p / (equations.gamma - 1) + 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1])) return U end diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl index 737e6e5..863b3e7 100644 --- a/test/test_2d_euler_potential_temperature.jl +++ b/test/test_2d_euler_potential_temperature.jl @@ -34,6 +34,33 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory end end +@trixiatmo_testset "elixir_euler_potential_temperature_gravity_wave" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_potential_temperature_gravity_wave.jl"), + l2=[ + 8.92434879991671e-9, + 1.8131676850378482e-7, + 2.6374650049135436e-5, + 1.4620631924879524e-6, + ], + linf=[ + 7.544232794032268e-8, + 1.654911628179434e-6, + 0.00023724858495745153, + 1.8884994915424613e-5, + ], + polydeg=3, + tspan=(0.0, 1.0)) + # 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 TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + @trixiatmo_testset "elixir_euler_potential_temperature_ec" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_potential_temperature_ec.jl"), From 6c91b292880c1d07b6f2f6d5001897c9469f3a97 Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Wed, 25 Sep 2024 16:39:37 +0200 Subject: [PATCH 27/27] minor changes --- examples/elixir_euler_potential_temperature_dry_bubble.jl | 4 ++-- examples/elixir_euler_potential_temperature_gravity_wave.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl index 1491ee8..33216b2 100644 --- a/examples/elixir_euler_potential_temperature_dry_bubble.jl +++ b/examples/elixir_euler_potential_temperature_dry_bubble.jl @@ -1,6 +1,6 @@ using OrdinaryDiffEq -using Trixi -using TrixiAtmo +using Trixi, TrixiAtmo +using TrixiAtmo: source_terms_gravity function initial_condition_warm_bubble(x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) diff --git a/examples/elixir_euler_potential_temperature_gravity_wave.jl b/examples/elixir_euler_potential_temperature_gravity_wave.jl index 64900ef..880e41c 100644 --- a/examples/elixir_euler_potential_temperature_gravity_wave.jl +++ b/examples/elixir_euler_potential_temperature_gravity_wave.jl @@ -1,6 +1,6 @@ using OrdinaryDiffEq -using Trixi -using TrixiAtmo +using Trixi, TrixiAtmo +using TrixiAtmo: source_terms_gravity function initial_condition_gravity_wave(x, t, equations::CompressibleEulerPotentialTemperatureEquations2D)