diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f6322dc5..89b0fd3c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -26,7 +26,7 @@ concurrency: jobs: test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -35,16 +35,14 @@ jobs: - '1.10' os: - ubuntu-latest - - macOS-latest + - macos-latest - windows-latest - arch: - - x64 steps: - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} - arch: ${{ matrix.arch }} + arch: default show-versioninfo: true - uses: julia-actions/cache@v2 diff --git a/examples/elixir_euler_potential_temperature_baroclinic_instability.jl b/examples/elixir_euler_potential_temperature_baroclinic_instability.jl index c04bdc18..466ec646 100644 --- a/examples/elixir_euler_potential_temperature_baroclinic_instability.jl +++ b/examples/elixir_euler_potential_temperature_baroclinic_instability.jl @@ -18,16 +18,16 @@ using LinearAlgebra: norm function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) # Parameters from Table 1 in the paper # Corresponding names in the paper are commented - radius_earth = 6.371229e6 # a - half_width_parameter = 2 # b - gravitational_acceleration = 9.81 # g - k = 3 # k - surface_pressure = 1.0f5 # p₀ - gas_constant = 287 # R - surface_equatorial_temperature = 310 # T₀ᴱ - surface_polar_temperature = 240 # T₀ᴾ - lapse_rate = 0.005 # Γ - angular_velocity = 7.29212e-5 # Ω + radius_earth = EARTH_RADIUS # a + half_width_parameter = 2 # b + gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g + k = 3 # k + surface_pressure = 1.0f5 # p₀ + gas_constant = 287 # R + surface_equatorial_temperature = 310 # T₀ᴱ + surface_polar_temperature = 240 # T₀ᴾ + lapse_rate = 0.005 # Γ + angular_velocity = EARTH_ROTATION_RATE # Ω # Distance to the center of the Earth r = z + radius_earth @@ -138,7 +138,7 @@ end function initial_condition_baroclinic_instability(x, t, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) lon, lat, r = cartesian_to_sphere(x) - radius_earth = 6.371229e6 + radius_earth = EARTH_RADIUS # Make sure that the r is not smaller than radius_earth z = max(r - radius_earth, 0.0) @@ -155,8 +155,7 @@ function initial_condition_baroclinic_instability(x, t, v1 = -sin(lon) * u - sin(lat) * cos(lon) * v v2 = cos(lon) * u - sin(lat) * sin(lon) * v v3 = cos(lat) * v - radius_earth = 6.371229e6 # a - gravitational_acceleration = 9.81 # g + gravitational_acceleration = EARTH_GRAVITATIONAL_ACCELERATION # g r = norm(x) # Make sure that r is not smaller than radius_earth @@ -174,8 +173,8 @@ end @inline function source_terms_baroclinic_instability(u, x, t, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D) - radius_earth = 6.371229e6 # a - angular_velocity = 7.29212e-5 # Ω + radius_earth = EARTH_RADIUS # a + angular_velocity = EARTH_ROTATION_RATE # Ω r = norm(x) # Make sure that r is not smaller than radius_earth @@ -195,7 +194,7 @@ end end equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004, c_v = 717, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) initial_condition = initial_condition_baroclinic_instability @@ -209,7 +208,7 @@ volume_flux = (flux_tec, flux_nonconservative_souza_etal) solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) trees_per_cube_face = (8, 4) -mesh = P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000, +mesh = P4estMeshCubedSphere(trees_per_cube_face..., EARTH_RADIUS, 30000, polydeg = polydeg, initial_refinement_level = 0) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, diff --git a/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl b/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl index aa365ab8..33732d63 100644 --- a/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl +++ b/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl @@ -38,7 +38,7 @@ end equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004, c_v = 717, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) # We have an isothermal background state with T0 = 250 K. # The reference speed of sound can be computed as: diff --git a/examples/elixir_euler_potential_temperature_linear_hydrostatic.jl b/examples/elixir_euler_potential_temperature_linear_hydrostatic.jl index 5859aadd..acea1d3b 100644 --- a/examples/elixir_euler_potential_temperature_linear_hydrostatic.jl +++ b/examples/elixir_euler_potential_temperature_linear_hydrostatic.jl @@ -97,7 +97,7 @@ end equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004, c_v = 717, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) alpha = 0.035 xr_B = 60000.0 linear_hydrostatic_setup = HydrostaticSetup(alpha, xr_B, equations) diff --git a/examples/elixir_euler_potential_temperature_linear_nonhydrostatic.jl b/examples/elixir_euler_potential_temperature_linear_nonhydrostatic.jl index 9ed7e74b..1fb5177d 100644 --- a/examples/elixir_euler_potential_temperature_linear_nonhydrostatic.jl +++ b/examples/elixir_euler_potential_temperature_linear_nonhydrostatic.jl @@ -93,7 +93,7 @@ end equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004, c_v = 717, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) alpha = 0.03 xr_B = 40000.0 diff --git a/examples/elixir_euler_potential_temperature_robert_bubble.jl b/examples/elixir_euler_potential_temperature_robert_bubble.jl index 229c5590..ed295296 100644 --- a/examples/elixir_euler_potential_temperature_robert_bubble.jl +++ b/examples/elixir_euler_potential_temperature_robert_bubble.jl @@ -12,7 +12,7 @@ function initial_condition_robert_bubble(x, t, # center of perturbation center_x = 500.0 center_z = 260.0 - g = 9.81 + g = EARTH_GRAVITATIONAL_ACCELERATION # radius of perturbation radius = 250.0 # distance of current x to center of perturbation @@ -47,7 +47,7 @@ end @inline function source_terms_gravity(u, x, t, equations::CompressibleEulerPotentialTemperatureEquations2D) rho, _, _, _ = u - g = 9.81 + g = EARTH_GRAVITATIONAL_ACCELERATION return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, zero(eltype(u))) end diff --git a/examples/elixir_euler_potential_temperature_schaer_mountain.jl b/examples/elixir_euler_potential_temperature_schaer_mountain.jl index 9c166c70..a7cb9900 100644 --- a/examples/elixir_euler_potential_temperature_schaer_mountain.jl +++ b/examples/elixir_euler_potential_temperature_schaer_mountain.jl @@ -93,7 +93,7 @@ end # semidiscretization of the compressible Euler equations equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004, c_v = 717, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) alpha = 0.03 xr_B = 20000 schär_setup = SchärSetup(alpha, xr_B) diff --git a/examples/elixir_euler_potential_temperature_well_balanced.jl b/examples/elixir_euler_potential_temperature_well_balanced.jl index 3d8cb16a..1ef76f37 100644 --- a/examples/elixir_euler_potential_temperature_well_balanced.jl +++ b/examples/elixir_euler_potential_temperature_well_balanced.jl @@ -34,7 +34,7 @@ end equations = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004, c_v = 717, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) polydeg = 2 basis = LobattoLegendreBasis(polydeg) cs = 340.0 diff --git a/examples/elixir_moist_euler_EC_bubble.jl b/examples/elixir_moist_euler_EC_bubble.jl index cdbea70e..fc0b456d 100644 --- a/examples/elixir_moist_euler_EC_bubble.jl +++ b/examples/elixir_moist_euler_EC_bubble.jl @@ -11,104 +11,8 @@ c_vd = 717 # specific heat at constant volume for dry air c_pv = 1885 # specific heat at constant pressure for moist air c_vv = 1424 # specific heat at constant volume for moist air equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv, - c_vv = c_vv, gravity = 9.81) - -function moist_state(y, dz, y0, r_t0, theta_e0, - equations::CompressibleMoistEulerEquations2D) - @unpack p_0, g, c_pd, c_pv, c_vd, c_vv, R_d, R_v, c_pl, L_00 = equations - (p, rho, T, r_t, r_v, rho_qv, theta_e) = y - p0 = y0[1] - - F = zeros(7, 1) - rho_d = rho / (1 + r_t) - p_d = R_d * rho_d * T - T_C = T - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) - L = L_00 - (c_pl - c_pv) * T - - F[1] = (p - p0) / dz + g * rho - F[2] = p - (R_d * rho_d + R_v * rho_qv) * T - # H = 1 is assumed - F[3] = (theta_e - - T * (p_d / p_0)^(-R_d / (c_pd + c_pl * r_t)) * - exp(L * r_v / ((c_pd + c_pl * r_t) * T))) - F[4] = r_t - r_t0 - F[5] = rho_qv - rho_d * r_v - F[6] = theta_e - theta_e0 - a = p_vs / (R_v * T) - rho_qv - b = rho - rho_qv - rho_d - # H=1 => phi=0 - F[7] = a + b - sqrt(a * a + b * b) - - return F -end - -function generate_function_of_y(dz, y0, r_t0, theta_e0, - equations::CompressibleMoistEulerEquations2D) - function function_of_y(y) - return moist_state(y, dz, y0, r_t0, theta_e0, equations) - end -end - -#Create Initial atmosphere by generating a layer data set -struct AtmosphereLayers{RealT <: Real} - equations::CompressibleMoistEulerEquations2D - # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql - layer_data::Matrix{RealT} # Contains the layer data for each height - total_height::RealT # Total height of the atmosphere - preciseness::Int # Space between each layer data (dz) - layers::Int # Amount of layers (total height / dz) - ground_state::NTuple{2, RealT} # (rho_0, p_tilde_0) to define the initial values at height z=0 - equivalent_potential_temperature::RealT # Value for theta_e since we have a constant temperature theta_e0=theta_e. - mixing_ratios::NTuple{2, RealT} # Constant mixing ratio. Also defines initial guess for rho_qv0 = r_v0 * rho_0. -end - -function AtmosphereLayers(equations; total_height = 10010.0, preciseness = 10, - ground_state = (1.4, 100000.0), - equivalent_potential_temperature = 320, - mixing_ratios = (0.02, 0.02), RealT = Float64) - @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl = equations - rho0, p0 = ground_state - r_t0, r_v0 = mixing_ratios - theta_e0 = equivalent_potential_temperature - - rho_qv0 = rho0 * r_v0 - T0 = theta_e0 - y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] - - n = convert(Int, total_height / preciseness) - dz = 0.01 - layer_data = zeros(RealT, n + 1, 4) - - F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations) - sol = nlsolve(F, y0) - p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero - - rho_d = rho / (1 + r_t) - rho_ql = rho - rho_d - rho_qv - kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) - rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) - - layer_data[1, :] = [rho, rho_theta, rho_qv, rho_ql] - for i in (1:n) - y0 = deepcopy(sol.zero) - dz = preciseness - F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations) - sol = nlsolve(F, y0) - p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero - - rho_d = rho / (1 + r_t) - rho_ql = rho - rho_d - rho_qv - kappa_M = (R_d * rho_d + R_v * rho_qv) / - (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) - rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) - - layer_data[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] - end - - return AtmosphereLayers{RealT}(equations, layer_data, total_height, dz, n, ground_state, - theta_e0, mixing_ratios) -end + c_vv = c_vv, + gravity = EARTH_GRAVITATIONAL_ACCELERATION) # Generate background state from the Layer data set by linearely interpolating the layers function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D, @@ -137,8 +41,8 @@ function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerE rho, rho_e, rho_qv, rho_ql = perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D) - v1 = 60.0 - v2 = 60.0 + v1 = 60 + v2 = 60 rho_v1 = rho * v1 rho_v2 = rho * v2 rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2) @@ -148,7 +52,7 @@ end # Add perturbation to the profile function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, - equations::CompressibleMoistEulerEquations2D) + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl, L_00 = equations xc = 2000 zc = 2000 @@ -166,16 +70,17 @@ function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, # Assume pressure stays constant if (r < rc && Δθ > 0) θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) - θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5 * r / rc)^2 / 300) + θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5f0 * r / rc)^2 / 300) rt = (rho_qv + rho_ql) / rho_d rv = rho_qv / rho_d θ_loc = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rv) if rt > 0 while true T_loc = θ_loc * (p_loc / p_0)^kappa - T_C = T_loc - 273.15 + T_C = T_loc - convert(RealT, 273.15) # SaturVapor - pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + pvs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) rho_d_new = (p_loc - pvs) / (R_d * T_loc) rvs = pvs / (R_v * rho_d_new * T_loc) θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) @@ -223,14 +128,13 @@ source_term = source_terms_geopotential # Get the DG approximation space polydeg = 4 -basis = LobattoLegendreBasis(polydeg) surface_flux = flux_chandrashekar volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) -solver = DGSEM(basis, surface_flux, volume_integral) +solver = DGSEM(polydeg, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (4000.0, 4000.0) diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index 4d61c1b4..a89b15ca 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -10,7 +10,8 @@ c_vd = 717 # specific heat at constant volume for dry air c_pv = 1885 # specific heat at constant pressure for moist air c_vv = 1424 # specific heat at constant volume for moist air equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv, - c_vv = c_vv, gravity = 9.81) + c_vv = c_vv, + gravity = EARTH_GRAVITATIONAL_ACCELERATION) # Warm bubble test from paper: # Wicker, L. J., and W. C. Skamarock, 1998: A time-splitting scheme @@ -18,15 +19,15 @@ equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c # time differencing. Mon. Wea. Rev., 126, 1992–1999. function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEquations2D) @unpack p_0, kappa, g, c_pd, c_vd, R_d, R_v = equations - xc = 10000.0 - zc = 2000.0 + xc = 10000 + zc = 2000 r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) - rc = 2000.0 - θ_ref = 300.0 - Δθ = 0.0 + rc = 2000 + θ_ref = 300 + Δθ = 0 if r <= rc - Δθ = 2 * cospi(0.5 * r / rc)^2 + Δθ = 2 * cospi(0.5f0 * r / rc)^2 end # Perturbed state: @@ -41,9 +42,9 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEq rho = p / ((p / p_0)^kappa * R_d * θ) T = p / (R_d * rho) - v1 = 20.0 - # v1 = 0.0 - v2 = 0.0 + v1 = 20 + # v1 = 0 + v2 = 0 rho_v1 = rho * v1 rho_v2 = rho * v2 rho_E = rho * c_vd * T + 1 / 2 * rho * (v1^2 + v2^2) @@ -61,7 +62,6 @@ source_term = source_terms_geopotential ############################################################################### # Get the DG approximation space polydeg = 4 -basis = LobattoLegendreBasis(polydeg) surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar @@ -70,7 +70,7 @@ volume_integral = VolumeIntegralFluxDifferencing(volume_flux) # Create DG solver with polynomial degree = 4 and LMARS flux as surface flux # and the EC flux (chandrashekar) as volume flux -solver = DGSEM(basis, surface_flux, volume_integral) +solver = DGSEM(polydeg, surface_flux, volume_integral) coordinates_min = (0.0, -5000.0) coordinates_max = (20000.0, 15000.0) diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 30b566d0..d5a17088 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -11,103 +11,8 @@ c_vd = 717 # specific heat at constant volume for dry air c_pv = 1885 # specific heat at constant pressure for moist air c_vv = 1424 # specific heat at constant volume for moist air equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv, - c_vv = c_vv, gravity = 9.81) - -function moist_state(y, dz, y0, r_t0, theta_e0, - equations::CompressibleMoistEulerEquations2D) - @unpack p_0, g, c_pd, c_pv, c_vd, c_vv, R_d, R_v, c_pl, L_00 = equations - (p, rho, T, r_t, r_v, rho_qv, theta_e) = y - p0 = y0[1] - - F = zeros(7, 1) - rho_d = rho / (1 + r_t) - p_d = R_d * rho_d * T - T_C = T - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) - L = L_00 - (c_pl - c_pv) * T - - F[1] = (p - p0) / dz + g * rho - F[2] = p - (R_d * rho_d + R_v * rho_qv) * T - # H = 1 is assumed - F[3] = (theta_e - - T * (p_d / p_0)^(-R_d / (c_pd + c_pl * r_t)) * - exp(L * r_v / ((c_pd + c_pl * r_t) * T))) - F[4] = r_t - r_t0 - F[5] = rho_qv - rho_d * r_v - F[6] = theta_e - theta_e0 - a = p_vs / (R_v * T) - rho_qv - b = rho - rho_qv - rho_d - # H=1 => phi=0 - F[7] = a + b - sqrt(a * a + b * b) - - return F -end - -function generate_function_of_y(dz, y0, r_t0, theta_e0, - equations::CompressibleMoistEulerEquations2D) - function function_of_y(y) - return moist_state(y, dz, y0, r_t0, theta_e0, equations) - end -end - -struct AtmosphereLayers{RealT <: Real} - equations::CompressibleMoistEulerEquations2D - # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql - layer_data::Matrix{RealT} - total_height::RealT - preciseness::Int - layers::Int - ground_state::NTuple{2, RealT} - equivalent_potential_temperature::RealT - mixing_ratios::NTuple{2, RealT} -end - -function AtmosphereLayers(equations; total_height = 10010.0, preciseness = 10, - ground_state = (1.4, 100000.0), - equivalent_potential_temperature = 320, - mixing_ratios = (0.02, 0.02), RealT = Float64) - @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl = equations - rho0, p0 = ground_state - r_t0, r_v0 = mixing_ratios - theta_e0 = equivalent_potential_temperature - - rho_qv0 = rho0 * r_v0 - T0 = theta_e0 - y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] - - n = convert(Int, total_height / preciseness) - dz = 0.01 - layer_data = zeros(RealT, n + 1, 4) - - F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations) - sol = nlsolve(F, y0) - p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero - - rho_d = rho / (1 + r_t) - rho_ql = rho - rho_d - rho_qv - kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) - rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) - - layer_data[1, :] = [rho, rho_theta, rho_qv, rho_ql] - for i in (1:n) - y0 = deepcopy(sol.zero) - dz = preciseness - F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations) - sol = nlsolve(F, y0) - p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero - - rho_d = rho / (1 + r_t) - rho_ql = rho - rho_d - rho_qv - kappa_M = (R_d * rho_d + R_v * rho_qv) / - (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) - rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) - - layer_data[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] - end - - return AtmosphereLayers{RealT}(equations, layer_data, total_height, dz, n, ground_state, - theta_e0, mixing_ratios) -end + c_vv = c_vv, + gravity = EARTH_GRAVITATIONAL_ACCELERATION) # Moist bubble test case from paper: # G.H. Bryan, J.M. Fritsch, A Benchmark Simulation for Moist Nonhydrostatic Numerical @@ -139,8 +44,8 @@ function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerE rho, rho_e, rho_qv, rho_ql = perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D) - v1 = 0.0 - v2 = 0.0 + v1 = 0 + v2 = 0 rho_v1 = rho * v1 rho_v2 = rho * v2 rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2) @@ -149,12 +54,12 @@ function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerE end function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, - equations::CompressibleMoistEulerEquations2D) + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl, L_00 = equations - xc = 10000.0 - zc = 2000.0 - rc = 2000.0 - Δθ = 2.0 + xc = 10000 + zc = 2000 + rc = 2000 + Δθ = 2 r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) rho_d = rho - rho_qv - rho_ql @@ -165,8 +70,9 @@ function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, p_v = rho_qv * R_v * T_loc p_d = p_loc - p_v - T_C = T_loc - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + T_C = T_loc - convert(RealT, 273.15) + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) H = p_v / p_vs r_v = rho_qv / rho_d r_l = rho_ql / rho_d @@ -184,7 +90,7 @@ function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, # Calculate background density potential temperature θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) # Calculate perturbed density potential temperature - θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5 * r / rc)^2 / 300) + θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5f0 * r / rc)^2 / 300) rt = (rho_qv + rho_ql) / rho_d rv = rho_qv / rho_d # Calculate moist potential temperature @@ -193,9 +99,10 @@ function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, if rt > 0 while true T_loc = θ_loc * (p_loc / p_0)^kappa - T_C = T_loc - 273.15 + T_C = T_loc - convert(RealT, 273.15) # SaturVapor - pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + pvs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) rho_d_new = (p_loc - pvs) / (R_d * T_loc) rvs = pvs / (R_v * rho_d_new * T_loc) θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) @@ -244,14 +151,13 @@ source_term = source_terms_moist_bubble # Get the DG approximation space polydeg = 4 -basis = LobattoLegendreBasis(polydeg) surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) -solver = DGSEM(basis, surface_flux, volume_integral) +solver = DGSEM(polydeg, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (20000.0, 10000.0) diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index 90da724d..576cdf28 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -13,12 +13,12 @@ using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential, # Review Vol. 128.4, pages 1153–1164, 2000, # https://doi.org/10.1175/1520-0493(2000)128<1153:BOFOSO>2.0.CO;2. function initial_condition_nonhydrostatic_gravity_wave(x, t, - equations::CompressibleMoistEulerEquations2D) + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack p_0, kappa, gamma, g, c_pd, c_vd, R_d, R_v = equations z = x[2] - T_0 = 280.0 + T_0 = 280 theta_0 = T_0 - N = 0.01 + N = convert(RealT, 0.01) theta = theta_0 * exp(N^2 * z / g) p = p_0 * (1 + g^2 * inv(c_pd * theta_0 * N^2) * (exp(-z * N^2 / g) - 1))^(1 / kappa) @@ -29,7 +29,7 @@ function initial_condition_nonhydrostatic_gravity_wave(x, t, v2 = 0 rho_v1 = rho * v1 rho_v2 = rho * v2 - rho_E = rho * c_vd * T + 0.5 * rho * (v1^2 + v2^2) + rho_E = rho * c_vd * T + 0.5f0 * rho * (v1^2 + v2^2) rho_qv = 0 rho_ql = 0 return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) @@ -40,7 +40,8 @@ c_vd = 717 # specific heat at constant volume for dry air c_pv = 1885 # specific heat at constant pressure for moist air c_vv = 1424 # specific heat at constant volume for moist air equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv, - c_vv = c_vv, gravity = 9.81) + c_vv = c_vv, + gravity = EARTH_GRAVITATIONAL_ACCELERATION) initial_condition = initial_condition_nonhydrostatic_gravity_wave @@ -59,34 +60,33 @@ boundary_conditions = (x_neg = boundary_condition_periodic, y_pos = boundary_condition_slip_wall) polydeg = 4 -basis = LobattoLegendreBasis(polydeg) surface_flux = FluxLMARS(360.0) volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) -solver = DGSEM(basis, surface_flux, volume_integral) +solver = DGSEM(polydeg, surface_flux, volume_integral) # Deformed rectangle that has "witch of agnesi" as bottom function bottom(x) - h = 400.0 - a = 1000.0 - x_length = 40000.0 + h = 400 + a = 1000 + x_length = 40000 # linear function cx for x in [-1,1] c = x_length / 2 # return (cx , f(cx)-f(c)) return SVector(c * x, (h * a^2 * inv((c * x)^2 + a^2)) - (h * a^2 * inv((c)^2 + a^2))) end -f1(s) = SVector(-20000.0, 8000.0 * s + 8000.0) -f2(s) = SVector(20000.0, 8000.0 * s + 8000.0) +f1(s) = SVector(-20000, 8000 * s + 8000) +f2(s) = SVector(20000, 8000 * s + 8000) function f3(s) - SVector(20000.0 * s, - (400.0 * 1000.0^2 * inv((20000.0 * s)^2 + 1000.0^2)) - - (400.0 * 1000.0^2 * inv((20000.0)^2 + 1000.0^2))) + SVector(20000 * s, + (400 * 1000.0^2 * inv((20000 * s)^2 + 1000^2)) - + (400 * 1000.0^2 * inv((20000)^2 + 1000^2))) end -f4(s) = SVector(20000.0 * s, 16000.0) +f4(s) = SVector(20000 * s, 16000) faces = (f1, f2, f3, f4) diff --git a/examples/elixir_moist_euler_source_terms_dry.jl b/examples/elixir_moist_euler_source_terms_dry.jl index 73b68533..19580fa1 100644 --- a/examples/elixir_moist_euler_source_terms_dry.jl +++ b/examples/elixir_moist_euler_source_terms_dry.jl @@ -10,7 +10,7 @@ c_vd = 717 # specific heat at constant volume for dry air c_pv = 1885 # specific heat at constant pressure for moist air c_vv = 1424 # specific heat at constant volume for moist air equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv, - c_vv = c_vv, gravity = 0.0) + c_vv = c_vv, gravity = 0) initial_condition = initial_condition_convergence_test_dry diff --git a/examples/elixir_moist_euler_source_terms_moist.jl b/examples/elixir_moist_euler_source_terms_moist.jl index 7187cce1..778761f5 100644 --- a/examples/elixir_moist_euler_source_terms_moist.jl +++ b/examples/elixir_moist_euler_source_terms_moist.jl @@ -11,16 +11,16 @@ c_vd = 717 # specific heat at constant volume for dry air c_pv = 1885 # specific heat at constant pressure for moist air c_vv = 1424 # specific heat at constant volume for moist air equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv, - c_vv = c_vv, gravity = 9.81) + c_vv = c_vv, + gravity = EARTH_GRAVITATIONAL_ACCELERATION) initial_condition = initial_condition_convergence_test_moist polydeg = 3 -basis = LobattoLegendreBasis(polydeg) surface_flux = flux_lax_friedrichs -solver = DGSEM(basis, surface_flux) +solver = DGSEM(polydeg, surface_flux) coordinates_min = (0.0, 0.0) coordinates_max = (2.0, 2.0) diff --git a/examples/elixir_moist_euler_source_terms_split_moist.jl b/examples/elixir_moist_euler_source_terms_split_moist.jl index 1bc1af6e..8526da66 100644 --- a/examples/elixir_moist_euler_source_terms_split_moist.jl +++ b/examples/elixir_moist_euler_source_terms_split_moist.jl @@ -11,19 +11,19 @@ c_vd = 717 # specific heat at constant volume for dry air c_pv = 1885 # specific heat at constant pressure for moist air c_vv = 1424 # specific heat at constant volume for moist air equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv, - c_vv = c_vv, gravity = 9.81) + c_vv = c_vv, + gravity = EARTH_GRAVITATIONAL_ACCELERATION) initial_condition = initial_condition_convergence_test_moist polydeg = 4 -basis = LobattoLegendreBasis(polydeg) surface_flux = flux_chandrashekar volume_flux = flux_chandrashekar volume_integral = VolumeIntegralFluxDifferencing(volume_flux) -solver = DGSEM(basis, surface_flux, volume_integral) +solver = DGSEM(polydeg, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (2.0, 2.0) diff --git a/examples/elixir_shallowwater_cartesian_advection_quad_icosahedron.jl b/examples/elixir_shallowwater_cartesian_advection_quad_icosahedron.jl index 70b1bfae..a5bcd240 100644 --- a/examples/elixir_shallowwater_cartesian_advection_quad_icosahedron.jl +++ b/examples/elixir_shallowwater_cartesian_advection_quad_icosahedron.jl @@ -20,7 +20,7 @@ cells_per_dimension = (2, 2) # We use the ShallowWaterEquations3D equations structure but modify the rhs! function to # convert it to a variable-coefficient advection equation -equations = ShallowWaterEquations3D(gravity = 0.0) +equations = ShallowWaterEquations3D(gravity = 0) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg = polydeg, @@ -38,7 +38,7 @@ function source_terms_convert_to_linear_advection(u, du, x, t, s3 = du[1] * v2 - du[3] s4 = du[1] * v3 - du[4] - return SVector(0.0, s2, s3, s4, 0.0) + return SVector(0, s2, s3, s4, 0) end # Hack to use the weak form kernel with ShallowWaterEquations3D (a non-conservative equation). diff --git a/examples/elixir_shallowwater_covariant_well_balanced.jl b/examples/elixir_shallowwater_covariant_well_balanced.jl index a9de6da3..a0bb1b02 100644 --- a/examples/elixir_shallowwater_covariant_well_balanced.jl +++ b/examples/elixir_shallowwater_covariant_well_balanced.jl @@ -10,7 +10,7 @@ using OrdinaryDiffEqLowStorageRK, Trixi, TrixiAtmo # Parameters # Initial condition is atmosphere at rest with constant total geopotential height -initial_condition = (x, t, equations) -> SVector(5960.0, 0.0, 0.0, 0.0) +initial_condition = (x, t, equations) -> SVector(5960, 0, 0, 0) polydeg = 3 cells_per_dimension = (5, 5) diff --git a/examples/moist_euler/convergence_test/convergence_test_rainy.jl b/examples/moist_euler/convergence_test/convergence_test_rainy.jl new file mode 100644 index 00000000..ef215c9f --- /dev/null +++ b/examples/moist_euler/convergence_test/convergence_test_rainy.jl @@ -0,0 +1,223 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiAtmo +using TrixiAtmo: saturation_vapour_pressure, saturation_vapour_pressure_derivative, + saturation_residual, saturation_residual_jacobian, + terminal_velocity_rain + +function initial_condition_convergence_test_rainy(x, t, + equations::CompressibleRainyEulerEquations2D{RealT}) where {RealT} + # constants + c_l = equations.c_liquid_water + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + R_v = equations.R_vapour + ref_L = equations.ref_latent_heat_vap_temp + + # define rho like in dry convergence test + c = 2 + A = convert(RealT, 0.1) + L = 2 + f = 1 / L + ω = 2 * pi * f + rho = c + A * sin(ω * (x[1] + x[2] - t)) + + # define variables of rho + temperature = rho + 250.0 + rho_vapour = saturation_vapour_pressure(temperature, equations) / (R_v * temperature) + rho_cloud = rho / c_l * 3000 + rho_rain = rho / c_l * 1000 + rho_moist = rho_vapour + rho_cloud + rho_dry = rho - rho_moist - rho_rain + + # define matching energy density with v1 := 1 and v2 := 1 , initially + energy = (c_vd * rho_dry + c_vv * rho_vapour + c_l * (rho_cloud + rho_rain)) * + temperature + energy += rho_vapour * ref_L + rho + + return SVector(rho_dry, rho_moist, rho_rain, rho, rho, energy, rho_vapour, rho_cloud, + temperature) +end + +function source_terms_convergence_test_rainy(u, x, t, + equations::CompressibleRainyEulerEquations2D{RealT}) where {RealT} + # constants + c_l = equations.c_liquid_water + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + R_d = equations.R_dry_air + R_v = equations.R_vapour + ref_L = equations.ref_latent_heat_vap_temp + N_0 = equations.rain_water_distr + v_0 = equations.v_mean_rain + + # help constant for terminal rain velocity derivative ( \Gamma(4.5) / 6 ~= 1.9386213994279082 ) + c_help = v_0 * convert(RealT, 1.9386213994279082) * (pi * N_0)^(-0.125f0) + + # define rho like initial condition + c = 2 + A = convert(RealT, 0.1) + L = 2 + f = 1 / L + ω = 2 * pi * f + si, co = sincos(ω * (x[1] + x[2] - t)) + rho = c + A * si + rho_x = ω * A * co + rho_t = -rho_x + + # define variables which depend on rho + temperature = rho + 250.0 + sat_vap_p = saturation_vapour_pressure(temperature, equations) + rho_vapour = sat_vap_p / (R_v * temperature) + rho_cloud = rho / c_l * 3000 + rho_rain = rho / c_l * 1000 + rho_moist = rho_vapour + rho_cloud + rho_dry = rho - rho_moist - rho_rain + vr = terminal_velocity_rain(rho_moist, rho_rain, equations) + + # define needed derivatives + sat_vap_p_t = rho_t * saturation_vapour_pressure_derivative(temperature, equations) + sat_vap_p_x = rho_x * saturation_vapour_pressure_derivative(temperature, equations) + + rho_vapour_t = (sat_vap_p_t * temperature - rho_t * sat_vap_p) / (R_v * temperature^2) + rho_vapour_x = (sat_vap_p_x * temperature - rho_x * sat_vap_p) / (R_v * temperature^2) + + rho_cloud_t = rho_t / c_l * 3000 + rho_cloud_x = rho_x / c_l * 3000 + + rho_rain_t = rho_t / c_l * 1000 + rho_rain_x = rho_x / c_l * 1000 + + rho_moist_t = rho_vapour_t + rho_cloud_t + rho_moist_x = rho_vapour_x + rho_cloud_x + + rho_dry_t = rho_t - rho_moist_t - rho_rain_t + rho_dry_x = rho_x - rho_moist_x - rho_rain_x + + energy_t = (c_vd * rho_dry_t + c_vv * rho_vapour_t + c_l * (rho_cloud_t + rho_rain_t)) * + temperature + energy_t += (c_vd * rho_dry + c_vv * rho_vapour + c_l * (rho_cloud + rho_rain)) * rho_t + energy_t += rho_vapour_t * ref_L + rho_t + + energy_x = (c_vd * rho_dry_x + c_vv * rho_vapour_x + c_l * (rho_cloud_x + rho_rain_x)) * + temperature + energy_x += (c_vd * rho_dry + c_vv * rho_vapour + c_l * (rho_cloud + rho_rain)) * rho_x + energy_x += rho_vapour_x * ref_L + rho_x + + pressure_x = (rho_dry_x * R_d + rho_vapour_x * R_v) * temperature + pressure_x += (rho_dry * R_d + rho_vapour * R_v) * rho_x # temperature_x = rho_x + + vr_x = c_help * 0.125f0 * + ((rho_rain_x * rho_moist - rho_rain * rho_moist_x) / (rho_moist + rho_rain)^2) + vr_x *= (rho_rain / (rho_moist + rho_rain))^(-0.875f0) + + rhor_vr__x = rho_rain_x * vr + rho_rain * vr_x + + # calculate source terms for manufactured solution + # density + S_rho_dry = rho_dry_t + 2 * rho_dry_x + S_rho_moist = rho_moist_t + 2 * rho_moist_x + S_rho_rain = rho_rain_t + 2 * rho_rain_x - rhor_vr__x + + # momentum + S_rho_v1 = rho_x + pressure_x - rhor_vr__x + S_rho_v2 = rho_x + pressure_x - rhor_vr__x + + # energy + S_energy = energy_t + 2 * (energy_x + pressure_x) - (c_l * rho_x * rho_rain * vr) + S_energy -= (c_l * temperature + 1) * rhor_vr__x + + return SVector(S_rho_dry, S_rho_moist, S_rho_rain, S_rho_v1, S_rho_v2, S_energy, 0, + 0, 0) +end + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleRainyEulerEquations2D() + +initial_condition = initial_condition_convergence_test_rainy + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (0.0, 0.0) +coordinates_max = (2.0, 2.0) + +cells_per_dimension = (6, 6) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test_rainy) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 10000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + 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) + +stage_limiter! = NonlinearSolveDG(saturation_residual, saturation_residual_jacobian, + SVector(7, 8, 9); tolerance = 1e-8) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, stage_limiter!), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# For copy-paste convenience: +#convergence_test("TrixiAtmo.jl/examples/convergence_test_elixirs/convergence_test_rainy.jl", 5) + +#= +#################################################################################################### +l2 +rho_dry rho_moist rho_rain rho_v1 rho_v2 energy_density rho_vapour rho_cloud temperature + +error EOC error EOC error EOC error EOC error EOC error EOC error EOC error EOC error EOC +2.40e-06 - 3.89e-05 - 1.52e-05 - 7.53e-03 - 2.09e-02 - 5.20e+01 - 1.75e-08 - 3.89e-05 - 2.16e-04 - + +1.43e-07 4.07 2.32e-06 4.07 9.72e-07 3.97 2.90e-04 4.70 1.39e-03 3.91 3.14e+00 4.05 1.15e-09 3.93 2.32e-06 4.07 1.42e-05 3.93 +8.48e-09 4.08 1.37e-07 4.08 1.20e-07 3.02 1.13e-05 4.68 3.10e-05 5.49 1.78e-01 4.14 2.69e-11 5.42 1.37e-07 4.08 3.32e-07 5.41 +5.29e-10 4.00 8.57e-09 4.00 5.27e-09 4.51 6.34e-07 4.16 9.54e-07 5.02 1.39e-02 3.68 1.24e-12 4.43 8.57e-09 4.00 1.54e-08 4.43 +3.31e-11 4.00 5.34e-10 4.00 2.02e-10 4.70 3.85e-08 4.04 3.80e-08 4.65 7.79e-04 4.16 7.09e-14 4.13 5.34e-10 4.00 8.78e-10 4.13 + +mean 4.04 mean 4.04 mean 4.05 mean 4.39 mean 4.77 mean 4.01 mean 4.48 mean 4.04 mean 4.48 +---------------------------------------------------------------------------------------------------- +linf +rho_dry rho_moist rho_rain rho_v1 rho_v2 energy_density rho_vapour rho_cloud temperature + +error EOC error EOC error EOC error EOC error EOC error EOC error EOC error EOC error EOC +1.68e-05 - 2.74e-04 - 8.69e-05 - 2.34e-02 - 6.25e-02 - 3.85e+02 - 6.00e-08 - 2.74e-04 - 7.43e-04 - + +1.07e-06 3.97 1.73e-05 3.98 4.36e-06 4.32 1.13e-03 4.38 3.11e-03 4.33 2.29e+01 4.07 5.72e-09 3.39 1.73e-05 3.98 7.09e-05 3.39 +6.01e-08 4.15 9.73e-07 4.16 3.46e-07 3.66 4.43e-05 4.67 6.53e-05 5.57 1.08e+00 4.40 1.63e-10 5.13 9.73e-07 4.16 2.00e-06 5.15 +3.69e-09 4.02 5.98e-08 4.03 2.91e-08 3.57 2.47e-06 4.16 2.37e-06 4.79 9.44e-02 3.52 9.40e-12 4.12 5.97e-08 4.03 1.15e-07 4.12 +2.30e-10 4.01 3.71e-09 4.01 1.47e-09 4.31 1.50e-07 4.04 1.16e-07 4.35 5.51e-03 4.10 5.70e-13 4.04 3.71e-09 4.01 6.97e-09 4.04 + +mean 4.04 mean 4.04 mean 3.96 mean 4.31 mean 4.76 mean 4.02 mean 4.17 mean 4.04 mean 4.18 +---------------------------------------------------------------------------------------------------- +Dict{Symbol, Any} with 3 entries: + :variables => ("rho_dry", "rho_moist", "rho_rain", "rho_v1", "rho_v2", "energy_density", "rho_vapour", "rho_cloud", "temperature") + :l2 => [4.03641, 4.03834, 4.04901, 4.39454, 4.7665, 4.00677, 4.47752, 4.03837, 4.47639] + :linf => [4.03824, 4.04306, 3.96283, 4.31327, 4.76078, 4.023, 4.17127, 4.04306, 4.17577] +=# diff --git a/examples/moist_euler/moist_bubble/elixir_rainy_euler_moist_bubble.jl b/examples/moist_euler/moist_bubble/elixir_rainy_euler_moist_bubble.jl new file mode 100644 index 00000000..df579415 --- /dev/null +++ b/examples/moist_euler/moist_bubble/elixir_rainy_euler_moist_bubble.jl @@ -0,0 +1,192 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiAtmo +using TrixiAtmo: source_terms_no_phase_change, saturation_residual, + saturation_residual_jacobian, NonlinearSolveDG, + cons2eq_pot_temp, flux_chandrashekar, + boundary_condition_simple_slip_wall +using NLsolve: nlsolve + +# Moist bubble test case from paper: +# G.H. Bryan, J.M. Fritsch, A Benchmark Simulation for Moist Nonhydrostatic Numerical +# Models, MonthlyWeather Review Vol.130, 2917–2928, 2002, +# https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml. +function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D, + atmosphere_layers::AtmosphereLayers) + @unpack layer_data, preciseness, total_height = atmosphere_layers + dz = preciseness + z = x[2] + if (z > total_height && !(isapprox(z, total_height))) + error("The atmosphere does not match the simulation domain") + end + n = convert(Int, floor((z + eps()) / dz)) + 1 + z_l = (n - 1) * dz + (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = layer_data[n, :] + z_r = n * dz + if (z_l == total_height) + z_r = z_l + dz + n = n - 1 + end + (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = layer_data[n + 1, :] + rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz + rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / + dz + rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz + rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz + + rho, rho_e, rho_qv, rho_ql, T_loc = perturb_moist_profile!(x, rho, rho_theta, rho_qv, + rho_ql, + equations::CompressibleMoistEulerEquations2D) + + v1 = 0 + v2 = 0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2) + + return SVector(rho - rho_qv - rho_ql, rho_qv + rho_ql, 0, rho_v1, rho_v2, rho_E, + rho_qv, rho_ql, T_loc) +end + +function perturb_moist_profile!(x, rho, rho_theta, rho_qv, rho_ql, + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} + @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl, L_00 = equations + xc = 10000 + zc = 2000 + rc = 2000 + Δθ = 2 + + r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + p_loc = p_0 * (R_d * rho_theta / p_0)^(1 / (1 - kappa_M)) + T_loc = p_loc / (R_d * rho_d + R_v * rho_qv) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + + # Assume pressure stays constant + if (r < rc && Δθ > 0) + # Calculate background density potential temperature + θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) + # Calculate perturbed density potential temperature + θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5f0 * r / rc)^2 / 300) + rt = (rho_qv + rho_ql) / rho_d + rv = rho_qv / rho_d + # Calculate moist potential temperature + θ_loc = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rv) + # Adjust varuables until the temperature is met + if rt > 0 + while true + T_loc = θ_loc * (p_loc / p_0)^kappa + T_C = T_loc - convert(RealT, 273.15) + # SaturVapor + pvs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) + rho_d_new = (p_loc - pvs) / (R_d * T_loc) + rvs = pvs / (R_v * rho_d_new * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + if abs(θ_new - θ_loc) <= θ_loc * 1.0e-12 + break + else + θ_loc = θ_new + end + end + else + rvs = 0 + T_loc = θ_loc * (p_loc / p_0)^kappa + rho_d_new = p_loc / (R_d * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + end + rho_qv = rvs * rho_d_new + rho_ql = (rt - rvs) * rho_d_new + rho = rho_d_new * (1 + rt) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + end + return SVector(rho, rho_e, rho_qv, rho_ql, T_loc) +end + +equations_moist = CompressibleMoistEulerEquations2D(c_pd = 1004, c_vd = 717, c_pv = 1885, + c_vv = 1424, + gravity = EARTH_GRAVITATIONAL_ACCELERATION) + +# Create background atmosphere data set +atmosphere_data = AtmosphereLayers(equations_moist) + +# Create the initial condition with the initial data set +function initial_condition_moist(x, t, equations::CompressibleRainyEulerEquations2D) + return initial_condition_moist_bubble(x, t, equations_moist, atmosphere_data) +end + +############################################################################### +# semidiscretization of the compressible rainy Euler equations + +equations = CompressibleRainyEulerEquations2D() + +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 + +surface_flux = flux_LMARS +volume_flux = flux_ec_rain + +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(polydeg, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +cells_per_dimension = (128, 64) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) + +semi = SemidiscretizationHyperbolic(mesh, equations, + initial_condition_moist, solver, + source_terms = source_terms_no_phase_change, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = 1000) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2eq_pot_temp) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +stage_limiter! = NonlinearSolveDG(saturation_residual, saturation_residual_jacobian, + SVector(7, 8, 9)) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, stage_limiter!), + maxiters = 1.0e7, + dt = 1.0, + save_everystep = false, callback = callbacks); diff --git a/examples/moist_euler/rainy_bubble/elixir_rainy_euler_rainy_bubble.jl b/examples/moist_euler/rainy_bubble/elixir_rainy_euler_rainy_bubble.jl new file mode 100644 index 00000000..22427026 --- /dev/null +++ b/examples/moist_euler/rainy_bubble/elixir_rainy_euler_rainy_bubble.jl @@ -0,0 +1,83 @@ +using OrdinaryDiffEqSSPRK +using Trixi +using TrixiAtmo +using TrixiAtmo: source_terms_rainy, initial_condition_bubble_rainy, + saturation_residual, saturation_residual_jacobian, + cons2eq_pot_temp, saturation_vapour_pressure, + boundary_condition_simple_slip_wall, + generate_hydrostatic_residual, generate_perturbation_residual +using NLsolve: nlsolve + +# domain +coordinates_min = (0.0, 0.0) +coordinates_max = (2400.0, 2400.0) + +# create layers for initial condition +equations = CompressibleRainyEulerEquations2D() +atmosphere_data = AtmosphereLayersRainyBubble(equations; + total_height = coordinates_max[2] + 1) + +# Create the initial condition with the initial data set +function initial_condition_rainy(x, t, equations::CompressibleRainyEulerEquations2D) + return initial_condition_bubble_rainy(x, t, equations; + atmosphere_layers = atmosphere_data) +end + +############################################################################### +# semidiscretization of the compressible rainy Euler equations + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_simple_slip_wall, + y_pos = boundary_condition_simple_slip_wall) + +polydeg = 3 + +surface_flux = flux_lax_friedrichs +volume_integral = VolumeIntegralFluxDifferencing(flux_ec_rain) + +solver = DGSEM(polydeg, surface_flux, volume_integral) + +cells_per_dimension = (64, 64) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_rainy, solver, + source_terms = source_terms_rainy, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 600.0) + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +# entropy +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2eq_pot_temp) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution) + +stage_limiter! = NonlinearSolveDG(saturation_residual, saturation_residual_jacobian, + SVector(7, 8, 9)) + +############################################################################### +# run the simulation +sol = solve(ode, SSPRK43(stage_limiter!); ode_default_options()..., + maxiters = 1.0e7, save_everystep = false, callback = callbacks); diff --git a/examples/moist_euler/rainy_bubble/elixir_rainy_euler_rainy_bubble_diffusion.jl b/examples/moist_euler/rainy_bubble/elixir_rainy_euler_rainy_bubble_diffusion.jl new file mode 100644 index 00000000..e23facd9 --- /dev/null +++ b/examples/moist_euler/rainy_bubble/elixir_rainy_euler_rainy_bubble_diffusion.jl @@ -0,0 +1,95 @@ +using OrdinaryDiffEqSSPRK +using Trixi +using TrixiAtmo +using TrixiAtmo: source_terms_rainy, initial_condition_bubble_rainy, + saturation_residual, saturation_residual_jacobian, + cons2eq_pot_temp, saturation_vapour_pressure, + source_terms_no_phase_change, + boundary_condition_laplace, + boundary_condition_simple_slip_wall +using NLsolve: nlsolve + +# domain +coordinates_min = (0.0, 0.0) +coordinates_max = (2400.0, 2400.0) + +# create layers for initial condition +equations = CompressibleRainyEulerEquations2D() +atmosphere_data = AtmosphereLayersRainyBubble(equations; + total_height = coordinates_max[2] + 1) + +# Create the initial condition with the initial data set +function initial_condition_rainy(x, t, equations::CompressibleRainyEulerEquations2D) + return initial_condition_bubble_rainy(x, t, equations; + atmosphere_layers = atmosphere_data) +end + +############################################################################### +# semidiscretization of the compressible rainy Euler equations + +diffusivity = 0.5f0 +equations_parabolic = LaplaceDiffusion2D(diffusivity, equations) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_simple_slip_wall, + y_pos = boundary_condition_simple_slip_wall) + +boundary_conditions_parabolic = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_laplace, + y_pos = boundary_condition_laplace) + +polydeg = 3 + +surface_flux = flux_lax_friedrichs +volume_integral = VolumeIntegralFluxDifferencing(flux_ec_rain) + +solver = DGSEM(polydeg, surface_flux, volume_integral) + +initial_refinement_level = 6 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = initial_refinement_level, + periodicity = (true, false), n_cells_max = 1_000_000) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition_rainy, solver; + source_terms = source_terms_rainy, + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 600.0) + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +# entropy +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2eq_pot_temp) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution) + +stage_limiter! = NonlinearSolveDG(saturation_residual, saturation_residual_jacobian, + SVector(7, 8, 9)) + +############################################################################### +# run the simulation +sol = solve(ode, SSPRK43(stage_limiter!); ode_default_options()..., + maxiters = 1.0e7, save_everystep = false, callback = callbacks); diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index ed316d03..24e83831 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -41,10 +41,12 @@ include("meshes/meshes.jl") include("semidiscretization/semidiscretization.jl") include("solvers/solvers.jl") include("callbacks_step/callbacks_step.jl") +include("callbacks_stage/callbacks_stage.jl") -export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, +export CompressibleMoistEulerEquations2D, + CompressibleRainyEulerEquations2D, CovariantLinearAdvectionEquation2D, CovariantShallowWaterEquations2D, - SplitCovariantShallowWaterEquations2D, + ShallowWaterEquations3D, SplitCovariantShallowWaterEquations2D, CompressibleEulerPotentialTemperatureEquations1D, CompressibleEulerPotentialTemperatureEquations2D, CompressibleEulerPotentialTemperatureEquations3D, @@ -56,6 +58,8 @@ export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, export GlobalCartesianCoordinates, GlobalSphericalCoordinates +export NonlinearSolveDG + export flux_chandrashekar, FluxLMARS export flux_nonconservative_zeros, flux_nonconservative_ec, @@ -63,7 +67,8 @@ export flux_nonconservative_zeros, flux_nonconservative_ec, source_terms_coriolis, source_terms_coriolis_lagrange_multiplier, flux_tec, flux_etec, flux_nonconservative_souza_etal, flux_nonconservative_artiano_etal, - flux_nonconservative_waruzewski_etal, flux_zero + flux_nonconservative_waruzewski_etal, flux_zero, + flux_ec_rain, flux_LMARS export source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier! @@ -85,6 +90,8 @@ export initial_condition_gaussian, initial_condition_geostrophic_balance, export bottom_topography_isolated_mountain, bottom_topography_unsteady_solid_body_rotation +export AtmosphereLayers, AtmosphereLayersRainyBubble + export examples_dir end # module TrixiAtmo diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl new file mode 100644 index 00000000..69661f16 --- /dev/null +++ b/src/callbacks_stage/callbacks_stage.jl @@ -0,0 +1 @@ +include("nonlinear_solve_dg.jl") diff --git a/src/callbacks_stage/nonlinear_solve_dg.jl b/src/callbacks_stage/nonlinear_solve_dg.jl new file mode 100644 index 00000000..63ac847a --- /dev/null +++ b/src/callbacks_stage/nonlinear_solve_dg.jl @@ -0,0 +1,47 @@ +using Trixi: wrap_array, AbstractSemidiscretization, TimerOutputs, @trixi_timeit, timer + +@muladd begin +#! format: noindent + +""" + NonlinearSolveDG + +Newton method, which can be called in every stage via callbacks. + +# Parameters +- `residual::Function`: function evaluating the residual +- `jacobian::Function`: function evaluating the Jacobian of the residual +- `variables_index_vector::Vector{Int64}`: vector of indices of entries of the solution vector the Newton method operates on +- `tolerance::Real`: tolerance for termination of the Newton method +- `max_iterations::Int64`: maximal number of iterations of the Newton method +""" +struct NonlinearSolveDG{RealT <: Real, ResidualFunctionT, JacobianFunctionT} + residual :: ResidualFunctionT + jacobian :: JacobianFunctionT + variables_index_vector :: Vector{Int64} + tolerance :: RealT + max_iterations :: Int64 + + function NonlinearSolveDG(residual, jacobian, variables_index_vector; + tolerance = 1e-9, max_iterations = 20) + new{typeof(tolerance), typeof(residual), typeof(jacobian)}(residual, jacobian, + variables_index_vector, + tolerance, + max_iterations) + end +end + +function (limiter!::NonlinearSolveDG)(u_ode, integrator, + semi::AbstractSemidiscretization, t) + u = wrap_array(u_ode, semi) + + @trixi_timeit timer() "nonlinear system solver" begin + nonlinear_solve_dg_2d!(u, limiter!.residual, limiter!.jacobian, + limiter!.variables_index_vector, + limiter!.tolerance, limiter!.max_iterations, + semi.equations, semi.solver, semi.cache, semi.mesh) + end +end +end + +include("nonlinear_solve_dg2d.jl") diff --git a/src/callbacks_stage/nonlinear_solve_dg2d.jl b/src/callbacks_stage/nonlinear_solve_dg2d.jl new file mode 100644 index 00000000..7d441f6f --- /dev/null +++ b/src/callbacks_stage/nonlinear_solve_dg2d.jl @@ -0,0 +1,45 @@ +using Trixi: @threaded, get_inverse_jacobian, set_node_vars!, get_node_coords + +@muladd begin +#! format: noindent + +function nonlinear_solve_dg_2d!(u, residual, jacobian, variables_index_vector, + tolerance, max_iterations, + equations::AbstractCompressibleRainyEulerEquations, + dg::DGSEM, cache, mesh) + # iterate over every DGSEM element + @threaded for element in eachelement(dg, cache) + # iterate over every node + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + guess = SVector(u_node[7], u_node[8], u_node[9]) + + # keep rain positive + if (u_node[3] < 0.0) + u[3, i, j, element] = 0 + end + + # newton method + for iteration in range(1, max_iterations) + res_vector = residual(u_node, guess, equations) + + if (maximum(abs.(res_vector)) < tolerance) + break + end + + jac_matrix = jacobian(u_node, guess, equations) + guess += -jac_matrix \ res_vector + + if iteration == max_iterations + @warn "Newton method: tolerance not met" + end + end + + # similar to set_node_vars! + for index in eachindex(variables_index_vector) + u[variables_index_vector[index], i, j, element] = guess[index] + end + end + end +end +end # muladd end diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl index 32a50a7e..92fc3984 100644 --- a/src/equations/compressible_moist_euler_2d_lucas.jl +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -4,18 +4,19 @@ #! format: noindent struct CompressibleMoistEulerEquations2D{RealT <: Real} <: AbstractCompressibleMoistEulerEquations{2, 6} - p_0::RealT # constant reference pressure 1000 hPa(100000 Pa) + p_0::RealT # constant reference pressure 1000 hPa(100000 Pa) c_pd::RealT # dry air constant c_vd::RealT # dry air constant - R_d::RealT # dry air gas constant + R_d::RealT # dry air gas constant c_pv::RealT # moist air constant c_vv::RealT # moist air constant - R_v::RealT # moist air gas constant - c_pl::RealT # liqid water constant - g::RealT # gravitation constant - 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 + R_v::RealT # moist air gas constant + c_pl::RealT # liqid water constant + g::RealT # gravitation constant + 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 + function CompressibleMoistEulerEquations2D(; c_pd, c_vd, c_pv, c_vv, gravity) c_pd, c_vd, c_pv, c_vv, g = promote(c_pd, c_vd, c_pv, c_vv, gravity) p_0 = 100_000 @@ -350,7 +351,7 @@ end z = x[2] # relaxed background velocity - vr1, vr2 = (10.0, 0.0) + vr1, vr2 = (10, 0) # damping threshold z_s = 9000.0 # boundary top @@ -401,7 +402,8 @@ end # Calculate Q_ph for a state u. # This source term models the phase chance between could water and vapor. -@inline function phase_change_term(u, equations::CompressibleMoistEulerEquations2D) +@inline function phase_change_term(u, + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack R_v = equations rho, _, _, _, rho_qv, rho_ql = u _, T = get_current_condition(u, equations) @@ -409,7 +411,8 @@ end T_C = T - 273.15 # saturation vapor pressure - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) # saturation density of vapor rho_star_qv = p_vs / (R_v * T) @@ -600,7 +603,8 @@ end end # Convert conservative variables to moisture related variables. -@inline function cons2moist(u, equations::CompressibleMoistEulerEquations2D) +@inline function cons2moist(u, + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack R_d, R_v, c_pd, c_pv, c_pl, p_0 = equations rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u @@ -612,7 +616,8 @@ end p_v = rho_qv * R_v * T T_C = T - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) H = p_v * inv(p_vs) rho_d = rho - (rho_qv + rho_ql) @@ -666,13 +671,15 @@ end return T end -@inline function saturation_pressure(u, equations::CompressibleMoistEulerEquations2D) +@inline function saturation_pressure(u, + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack R_v = equations rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u T = get_current_condition(u, equations)[2] p_v = rho_qv * R_v * T T_C = T - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) H = p_v / p_vs return H end @@ -799,15 +806,16 @@ end # Calculate the equivalent potential temperature for a conservative state `cons`. @inline function equivalent_pottemp_thermodynamic(cons, - equations::CompressibleMoistEulerEquations2D) + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack c_pd, c_pv, c_pl, R_d, R_v, p_0, kappa, L_00 = equations rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = cons rho_d = rho - rho_qv - rho_ql p, T = get_current_condition(cons, equations) p_v = rho_qv * R_v * T p_d = p - p_v - T_C = T - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + T_C = T - convert(RealT, 273.15) + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) H = p_v / p_vs r_v = rho_qv / rho_d r_l = rho_ql / rho_d @@ -825,15 +833,17 @@ end # Convert conservative variables to primitive varuables with # equivalent potential temperature instead of pressure # and mixing ratios innstead of specific contend. -@inline function cons2aeqpot(cons, equations::CompressibleMoistEulerEquations2D) +@inline function cons2aeqpot(cons, + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} @unpack c_pd, c_pv, c_pl, R_d, R_v, p_0, L_00 = equations rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = cons rho_d = rho - rho_qv - rho_ql p, T = get_current_condition(cons, equations) p_v = rho_qv * R_v * T p_d = p - p_v - T_C = T - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + T_C = T - convert(RealT, 273.15) + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) H = p_v / p_vs r_v = rho_qv / rho_d r_l = rho_ql / rho_d @@ -992,6 +1002,106 @@ end return SVector(f1, f2, f3, f4, f_1v, f_1l) end +function moist_state(y, dz, y0, r_t0, theta_e0, + equations::CompressibleMoistEulerEquations2D{RealT}) where {RealT} + @unpack p_0, g, c_pd, c_pv, c_vd, c_vv, R_d, R_v, c_pl, L_00 = equations + (p, rho, T, r_t, r_v, rho_qv, theta_e) = y + p0 = y0[1] + + F = zeros(7, 1) + rho_d = rho / (1 + r_t) + p_d = R_d * rho_d * T + T_C = T - convert(RealT, 273.15) + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) + L = L_00 - (c_pl - c_pv) * T + + F[1] = (p - p0) / dz + g * rho + F[2] = p - (R_d * rho_d + R_v * rho_qv) * T + # H = 1 is assumed + F[3] = (theta_e - + T * (p_d / p_0)^(-R_d / (c_pd + c_pl * r_t)) * + exp(L * r_v / ((c_pd + c_pl * r_t) * T))) + F[4] = r_t - r_t0 + F[5] = rho_qv - rho_d * r_v + F[6] = theta_e - theta_e0 + a = p_vs / (R_v * T) - rho_qv + b = rho - rho_qv - rho_d + # H=1 => phi=0 + F[7] = a + b - sqrt(a * a + b * b) + + return F +end + +function generate_function_of_y(dz, y0, r_t0, theta_e0, + equations::CompressibleMoistEulerEquations2D) + function function_of_y(y) + return moist_state(y, dz, y0, r_t0, theta_e0, equations) + end +end + +#Create Initial atmosphere by generating a layer data set +struct AtmosphereLayers{RealT <: Real} + equations::CompressibleMoistEulerEquations2D + # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql + layer_data::Matrix{RealT} # Contains the layer data for each height + total_height::RealT # Total height of the atmosphere + preciseness::Int # Space between each layer data (dz) + layers::Int # Amount of layers (total height / dz) + ground_state::NTuple{2, RealT} # (rho_0, p_tilde_0) to define the initial values at height z=0 + equivalent_potential_temperature::RealT # Value for theta_e since we have a constant temperature theta_e0=theta_e. + mixing_ratios::NTuple{2, RealT} # Constant mixing ratio. Also defines initial guess for rho_qv0 = r_v0 * rho_0. +end + +function AtmosphereLayers(equations; total_height = 10010, preciseness = 10, + ground_state = (1.4, 100000), + equivalent_potential_temperature = 320, + mixing_ratios = (0.02, 0.02), RealT = Float64) + @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl = equations + rho0, p0 = convert.(RealT, ground_state) + r_t0, r_v0 = convert.(RealT, mixing_ratios) + theta_e0 = equivalent_potential_temperature + + rho_qv0 = rho0 * r_v0 + T0 = theta_e0 + y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] + + n = convert(Int, total_height / preciseness) + dz = convert(RealT, 0.01) + layer_data = zeros(RealT, n + 1, 4) + + F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations) + sol = nlsolve(F, y0) + p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero + + rho_d = rho / (1 + r_t) + rho_ql = rho - rho_d - rho_qv + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + layer_data[1, :] = [rho, rho_theta, rho_qv, rho_ql] + for i in (1:n) + y0 = deepcopy(sol.zero) + dz = preciseness + F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations) + sol = nlsolve(F, y0) + p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero + + rho_d = rho / (1 + r_t) + rho_ql = rho - rho_d - rho_qv + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + layer_data[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] + end + + return AtmosphereLayers{RealT}(equations, layer_data, total_height, dz, n, + ground_state, + theta_e0, mixing_ratios) +end + varnames(::typeof(cons2cons), ::CompressibleMoistEulerEquations2D) = ("rho", "rho_v1", "rho_v2", "rho_E", "rho_qv", diff --git a/src/equations/compressible_rainy_euler_2d.jl b/src/equations/compressible_rainy_euler_2d.jl new file mode 100644 index 00000000..f663c19e --- /dev/null +++ b/src/equations/compressible_rainy_euler_2d.jl @@ -0,0 +1,1205 @@ +# Implemented by Fabian Höck following +# Sabine Doppler, Philip L. Lederer, Joachim Schöberl, Henry von Wahl +# A discontinuous Galerkin approach for atmospheric flows with implicit condensation +# Journal of Computational Physics, Volume 499, 2024 +# https://doi.org/10.1016/j.jcp.2023.112713 + +using Trixi +using NLsolve: nlsolve +import Trixi: varnames, cons2prim, cons2entropy, + flux, max_abs_speeds, max_abs_speed_naive, + boundary_condition_slip_wall, LaplaceDiffusion2D + +@muladd begin +#! format: noindent + +### equation, parameters and constants ### + +@doc raw""" + CompressibleRainyEulerEquations2D{RealT <: Real} <: + AbstractCompressibleRainyEulerEquations{2, 9} + +The compressible Euler equations in two dimensions with gravity and separate densities for +dry air ($\rho_d$), moist air ($\rho_m$), and rain ($\rho_r). +```math +\begin{alignat}{4} + \partial_t \, &&(\rho_d) &+ \nabla \cdot (\rho_d\, v) &&= & 0\,,\\ + \partial_t \, &&(\rho_m) &+ \nabla \cdot (\rho_m\, v) &&= & -S_r \,,\\ + \partial_t \, &&(\rho_r) &+ \nabla \cdot (\rho_r\, (v-v_r\, e_z)) &&= & S_r\,,\\ + \partial_t\,&&(\rho v) &+ \nabla \cdot (\rho v \otimes v - \rho_r\, v_r \,v\otimes e_z + p\cdot\textrm{Id}) &&= & -\rho\, g\, e_z\,,\\ + \partial_t\, &&(E) &+ \nabla \cdot ((E + p)\,v - (c_l\, T + 1/2\, \langle v, v\rangle)\, \rho_r\, v_r\, e_z) &&= &\,\,-\rho\, g\, \langle e_z, v\rangle\,. +\end{alignat} +``` +Here, moist air is the sum of vapor ($\rho_v$) and cloud water ($\rho_c$), which are given implicitly by the nonlinear system +```math +\begin{align} + (c_{vd}\, \rho_d + c_{vv}\, \rho_v + c_l\, \rho_l)\,T + \rho_v\, L_{\textrm{ref}} + 1/2\,\rho\,\langle v, v\rangle - E &= 0\,,\\ + \min \left( \frac{e_s(T)}{R_v\,T}, \rho_m\right) - \rho_v &=0\,,\\ + \rho_m - \rho_v - \rho_c &=0\,, +\end{align} +``` +where $e_s(T)$ is the modeled saturation vapor pressure, $v_r$ the modeled terminal rain, +fall velocity, and $S_r$ a modeled source term for rain conversion. + +## Reference +S. Doppler et al. (2024). A discontinuous Galerkin approach for atmospheric flows with +implicit condensation. Journal of Computational Physics. 499:112713. +[DOI: 10.1016/j.jcp.2023.112713](https://doi.org/10.1016/j.jcp.2023.112713) +""" +struct CompressibleRainyEulerEquations2D{RealT <: Real} <: + AbstractCompressibleRainyEulerEquations{2, 9} + # Specific heat capacities: + c_liquid_water :: RealT + c_dry_air_const_pressure :: RealT + c_dry_air_const_volume :: RealT + c_vapour_const_pressure :: RealT + c_vapour_const_volume :: RealT + + # Gas constants: + R_dry_air :: RealT + R_vapour :: RealT + eps :: RealT + + # Reference values: + ref_saturation_pressure :: RealT + ref_temperature :: RealT + ref_latent_heat_vap_temp :: RealT + ref_pressure :: RealT + + # Other: + gravity :: RealT + rain_water_distr :: RealT + v_mean_rain :: RealT +end + +function CompressibleRainyEulerEquations2D(RealT = Float64) + # Specific heat capacities: + c_liquid_water = 4186 + c_dry_air_const_pressure = 1004 + c_dry_air_const_volume = 717 + c_vapour_const_pressure = 1885 + c_vapour_const_volume = 1424 + + # Gas constants: + R_dry_air = c_dry_air_const_pressure - c_dry_air_const_volume + R_vapour = c_vapour_const_pressure - c_vapour_const_volume + eps = convert(RealT, R_dry_air / R_vapour) + + # Reference values: + ref_saturation_pressure = convert(RealT, 610.7) # This needs to be adjusted if ref_temperature is changed! + ref_temperature = convert(RealT, 273.15) + ref_latent_heat_vap_temp = 2500000 #3147620.0 + ref_pressure = 100000 + + # Other: + gravity = convert(RealT, EARTH_GRAVITATIONAL_ACCELERATION) + rain_water_distr = 8000000 + v_mean_rain = 130 + + return CompressibleRainyEulerEquations2D{RealT}(c_liquid_water, + c_dry_air_const_pressure, + c_dry_air_const_volume, + c_vapour_const_pressure, + c_vapour_const_volume, R_dry_air, + R_vapour, eps, + ref_saturation_pressure, + ref_temperature, + ref_latent_heat_vap_temp, + ref_pressure, gravity, + rain_water_distr, v_mean_rain) +end + +### conversion ### + +@inline function cons2prim(u, equations::CompressibleRainyEulerEquations2D) + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + + # energy density + energy = energy_density(u, equations) + + # nonlinear system + rho_vapour, rho_cloud, temperature = cons2nonlinearsystemsol(u, equations) + + return SVector(rho_dry, rho_moist, rho_rain, v1, v2, energy, rho_vapour, rho_cloud, + temperature) +end + +# converts consverved to entropy variables +@inline function cons2entropy(u, equations::CompressibleRainyEulerEquations2D) + # constants + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + c_l = equations.c_liquid_water + R_d = equations.R_dry_air + R_v = equations.R_vapour + L_ref = equations.ref_latent_heat_vap_temp + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # nonlinear system + rho_vapour, rho_cloud, temperature = cons2nonlinearsystemsol(u, equations) + ln_temperature = log(temperature) + inv_temperature = inv(temperature) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + v_squared_temp = 0.5f0 * (v1^2 + v2^2) * inv_temperature + + # check for zero density + rho_vapour_log = 0 + + if (rho_vapour > 0.0) + rho_vapour_log = log(rho_vapour) + end + + omega_dry = c_vd * ln_temperature - R_d * log(rho_dry) + v_squared_temp - c_vd - R_d + omega_vapour = c_vv * ln_temperature - R_v * rho_vapour_log + v_squared_temp - + c_vv - R_v - L_ref * inv_temperature + omega_liquid = c_l * ln_temperature + v_squared_temp - c_l + omega_momentum_1 = -v1 * inv_temperature + omega_momentum_2 = -v2 * inv_temperature + omega_energy = inv_temperature + + return SVector(omega_dry, omega_vapour + omega_liquid, omega_liquid, + omega_momentum_1, omega_momentum_2, omega_energy, 0, 0, 0) +end + +# adapted from compressible_moist_euler_2d.jl +@inline function cons2eq_pot_temp(u, + equations::CompressibleRainyEulerEquations2D{RealT}) where {RealT} + # constants + c_l = equations.c_liquid_water + c_pd = equations.c_dry_air_const_pressure + c_pv = equations.c_vapour_const_pressure + R_d = equations.R_dry_air + R_v = equations.R_vapour + ref_p = equations.ref_pressure + ref_temp = equations.ref_temperature + ref_L = equations.ref_latent_heat_vap_temp + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + + # energy density + energy = energy_density(u, equations) + + # nonlinear system + rho_vapour, rho_cloud, temperature = cons2nonlinearsystemsol(u, equations) + + # pressure + p = pressure(u, equations) + + p_v = rho_vapour * R_v * temperature + p_d = p - p_v + T_C = temperature - ref_temp + p_vs = convert(RealT, 611.2) * + exp(convert(RealT, 17.62) * T_C / (convert(RealT, 243.12) + T_C)) + H = p_v / p_vs + r_v = rho_vapour / rho_dry + r_c = rho_cloud / rho_dry + r_r = rho_rain / rho_dry + L_v = ref_L + (c_pv - c_l) * temperature + c_p = c_pd + (r_v + r_c + r_r) * c_l + + # equivalent potential temperature + eq_pot = (temperature * (ref_p / p_d)^(R_d / c_p) * H^(-r_v * R_v / c_p) * + exp(L_v * r_v * inv(c_p * temperature))) + + return SVector(rho_dry, rho_vapour, rho_cloud, rho_rain, v1, v2, eq_pot, p) +end + +@inline function cons2nonlinearsystemsol(u, + equations::CompressibleRainyEulerEquations2D) + # constants + c_vd = equations.c_dry_air_const_volume + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + + # energy density + energy = energy_density(u, equations) + + # recover temperature explicitly from energy when other variables are zero + if (rho_moist == 0.0 && rho_rain == 0.0) + # energy density definition without ref_temp for dry case + energy_kinetic = 0.5f0 * (v1^2 + v2^2) * rho + temperature = (energy - energy_kinetic) / (c_vd * rho) #+ ref_temp + + if (temperature < 0.0) + error("temp negative") + end + + return SVector(0, 0, temperature) + else + # experimental and overly simple positivity check + rho_vapour = u[7] + rho_cloud = u[8] + temperature = u[9] + + if (rho_vapour < 0.0 && isapprox(rho_vapour, 0.0, atol = 1e-15)) + rho_vapour = 0 + end + + if (rho_cloud < 0.0 && isapprox(rho_cloud, 0.0, atol = 1e-15)) + rho_cloud = 0 + end + + return SVector(rho_vapour, rho_cloud, temperature) + end +end + +# for convenience TODO rename +@inline function cons2speeds(u, equations::CompressibleRainyEulerEquations2D) + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + + # get speed of sound + v_sound = speed_of_sound(u, equations)[1] + + # get terminal velocity rain + v_r = terminal_velocity_rain(rho_moist, rho_rain, equations) + + return SVector(v1, v2, v_sound, v_r) +end + +### varnames ### + +varnames(::typeof(cons2cons), ::CompressibleRainyEulerEquations2D) = ("rho_dry", + "rho_moist", + "rho_rain", + "rho_v1", + "rho_v2", + "energy_density", + "rho_vapour", + "rho_cloud", + "temperature") + +varnames(::typeof(cons2prim), ::CompressibleRainyEulerEquations2D) = ("rho_dry", + "rho_moist", + "rho_rain", + "v1", "v2", + "energy_density", + "rho_vapour", + "rho_cloud", + "temperature") + +varnames(::typeof(cons2eq_pot_temp), ::CompressibleRainyEulerEquations2D) = ("rho_dry", + "rho_vapour", + "rho_cloud", + "rho_rain", + "v1", "v2", + "eq_pot_temp", + "pressure") + +### physics variables ### + +@inline function densities(u, equations::CompressibleRainyEulerEquations2D) + rho_dry = u[1] + rho_moist = u[2] + rho_rain = u[3] + rho = rho_dry + rho_moist + rho_rain + rho_inv = inv(rho) + + return SVector(rho_dry, rho_moist, rho_rain, rho, rho_inv) +end + +@inline function velocities(u, rho_inv, equations::CompressibleRainyEulerEquations2D) + return SVector(u[4] * rho_inv, u[5] * rho_inv) +end + +@inline function energy_density(u, equations::CompressibleRainyEulerEquations2D) + return u[6] +end + +@inline function pressure(u, equations::CompressibleRainyEulerEquations2D) + # constants + R_d = equations.R_dry_air + R_v = equations.R_vapour + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + rho_vapour, _, temperature = cons2nonlinearsystemsol(u, equations) + + p = (R_d * rho_dry + R_v * rho_vapour) * temperature + + return p +end + +@inline function speed_of_sound(u, equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + R_d = equations.R_dry_air + R_v = equations.R_vapour + + # densities + rho_dry, _, rho_rain, _, rho_inv = densities(u, equations) + + # recover rho_vapour, rho_cloud, temperature from nonlinear system + rho_vapour, rho_cloud, temperature = cons2nonlinearsystemsol(u, equations) + if (rho_vapour < 0.0) + error("rho vapour less than zero") + end + if (rho_cloud < 0.0) + error("rho cloud less than zero") + end + + # formula + p = pressure(u, equations) + q_v = rho_vapour / rho_dry + q_l = (rho_cloud + rho_rain) / rho_dry + gamma_m = 1 + (R_d + R_v * q_v) / (c_vd + c_vv * q_v + c_l * q_l) + + if (rho_inv < 0.0) + error("rho less than zero") + elseif (p < 0.0) + error("pressure less than zero") + end + + v_sound = sqrt(gamma_m * p * rho_inv) + + return SVector(v_sound, gamma_m) +end + +@inline function terminal_velocity_rain(rho_moist, rho_rain, + equations::CompressibleRainyEulerEquations2D{RealT}) where {RealT} + # constants + N_0 = equations.rain_water_distr + v_0 = equations.v_mean_rain + + # formula ( \Gamma(4.5) / 6 ~= 1.9386213994279082 ) + if (rho_rain > 0.0) + v_terminal_rain = v_0 * convert(RealT, 1.9386213994279082) * + (rho_rain / (pi * (rho_moist + rho_rain) * N_0))^(0.125f0) + else + v_terminal_rain = RealT(0) + end + + return v_terminal_rain +end + +@inline function saturation_vapour_pressure(temperature, + equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + c_pv = equations.c_vapour_const_pressure + R_v = equations.R_vapour + ref_s_p = equations.ref_saturation_pressure + ref_temp = equations.ref_temperature + ref_L = equations.ref_latent_heat_vap_temp + + # testing + if (temperature < 0.0) + display(temperature) + error("temp less than zero") + end + + # Clausius Clapeyron formula + p_vapour_saturation = ref_s_p * (temperature / ref_temp)^((c_pv - c_l) / R_v) + p_vapour_saturation *= exp(((ref_L - (c_pv - c_l) * ref_temp) / R_v) * + (1 / ref_temp - 1 / temperature)) + + return p_vapour_saturation +end + +@inline function saturation_vapour_pressure_derivative(temperature, + equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + c_pv = equations.c_vapour_const_pressure + R_v = equations.R_vapour + ref_s_p = equations.ref_saturation_pressure + ref_temp = equations.ref_temperature + ref_L = equations.ref_latent_heat_vap_temp + + # testing + if (temperature < 0.0) + display(temperature) + error("temp less than zero") + end + + const_1 = (c_pv - c_l) / R_v + const_2 = (ref_L - (c_pv - c_l) * ref_temp) / R_v + + p_vapour_saturation_derivative = ref_s_p / (ref_temp^const_1) + p_vapour_saturation_derivative *= (const_1 * temperature^(const_1 - 1) + + const_2 * temperature^(const_1 - 2)) + p_vapour_saturation_derivative *= exp(const_2 * (1 / ref_temp - 1 / temperature)) + + return p_vapour_saturation_derivative +end + +### pde discretization ### + +@inline function flux(u, orientation::Integer, + equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + + #densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # recover rho_vapour, rho_cloud, temperature from nonlinear system + _, _, temperature = cons2nonlinearsystemsol(u, equations) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + v_r = terminal_velocity_rain(rho_moist, rho_rain, equations) + + # pressure + p = pressure(u, equations) + + # energy density + energy = energy_density(u, equations) + + # flux for orientation cases + if (orientation == 1) + # mass + f1 = rho_dry * v1 + f2 = rho_moist * v1 + f3 = rho_rain * v1 + + # momentum + f4 = rho * v1 * v1 + p + f5 = rho * v1 * v2 + + # energy + f6 = (energy + p) * v1 + + else + # mass + f1 = rho_dry * v2 + f2 = rho_moist * v2 + f3 = rho_rain * (v2 - v_r) + + # momentum + f4 = rho * v1 * v2 - rho_rain * v_r * v1 + f5 = rho * v2 * v2 + p - rho_rain * v_r * v2 + + # energy + f6 = (energy + p) * v2 - + (c_l * temperature + 0.5f0 * (v1^2 + v2^2)) * rho_rain * v_r + end + + return SVector(f1, f2, f3, f4, f5, f6, + 0, 0, 0) +end + +@inline function flux(u, normal_direction::AbstractVector, + equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # recover rho_vapour, rho_cloud, temperature from nonlinear system + _, _, temperature = cons2nonlinearsystemsol(u, equations) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + v_r = terminal_velocity_rain(rho_moist, rho_rain, equations) + + # normal velocities + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + v_r_normal = v_r * normal_direction[2] + + # pressure + p = pressure(u, equations) + + # energy density + energy = energy_density(u, equations) + + # flux + # mass + f1 = rho_dry * v_normal + f2 = rho_moist * v_normal + f3 = rho_rain * (v_normal - v_r_normal) + + # momentum + f4 = rho * v_normal * v1 + p * normal_direction[1] - rho_rain * v_r_normal * v1 + f5 = rho * v_normal * v2 + p * normal_direction[2] - rho_rain * v_r_normal * v2 + + # energy + f6 = (energy + p) * v_normal - + (c_l * temperature + 0.5f0 * (v1^2 + v2^2)) * rho_rain * v_r_normal + + return SVector(f1, f2, f3, f4, f5, f6, + 0, 0, 0) +end + +# no Coriolis term +@inline function source_terms_rainy(u, x, t, + equations::CompressibleRainyEulerEquations2D{RealT}) where {RealT} + # constants + R_v = equations.R_vapour + ref_temp = equations.ref_temperature + g = equations.gravity + + # name needed variables + rho_v2 = u[5] + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # recover rho_vapour, rho_cloud, temperature from nonlinear system + rho_vapour, rho_cloud, temperature = cons2nonlinearsystemsol(u, equations) + + rho_vs = saturation_vapour_pressure(temperature, equations) / (R_v * temperature) + + # source terms phase change + S_evaporation = (convert(RealT, 3.86e-3) - + convert(RealT, 9.41e-5) * (temperature - ref_temp)) * + (1 + convert(RealT, 9.1) * rho_rain^(0.1875f0)) + S_evaporation *= (rho_vs - rho_vapour) * rho_rain^(0.5f0) + S_auto_conversion = convert(RealT, 0.001) * rho_cloud + S_accretion = convert(RealT, 1.72) * rho_cloud * rho_rain^(0.875f0) + S_rain = S_auto_conversion + S_accretion - S_evaporation + + return SVector(0, -S_rain, S_rain, 0, + -rho * g, -rho_v2 * g, 0, 0, 0) +end + +# no phase changes and no Coriolis term +@inline function source_terms_no_phase_change(u, x, t, + equations::CompressibleRainyEulerEquations2D) + # constants + g = equations.gravity + + # name needed variables + rho_v2 = u[5] + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + return SVector(0, 0, 0, 0, + -g * rho, -g * rho_v2, 0, 0, 0) +end + +@inline function max_abs_speeds(u, equations::CompressibleRainyEulerEquations2D) + # name needed variables + v1, v2, v_sound, v_r = cons2speeds(u, equations) + + return SVector((abs(v1) + v_sound), (abs(v2) + v_sound + abs(v_r))) +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleRainyEulerEquations2D) + # name needed variables + v1_ll, v2_ll, v_sound_ll, v_r_ll = cons2speeds(u_ll, equations) + v1_rr, v2_rr, v_sound_rr, v_r_rr = cons2speeds(u_rr, equations) + + # calculate upper bounds for left and right speed + v_ll_max = abs(v1_ll * normal_direction[1] + v2_ll * normal_direction[2]) + v_ll_max += abs(v_r_ll * normal_direction[2]) + + v_rr_max = abs(v1_rr * normal_direction[1] + v2_rr * normal_direction[2]) + v_rr_max += abs(v_r_rr * normal_direction[2]) + + return max(v_ll_max, v_rr_max) + + max(v_sound_ll, v_sound_rr) * norm(normal_direction) +end + +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleRainyEulerEquations2D) + # name needed variables + v1_ll, v2_ll, v_sound_ll, v_r_ll = cons2speeds(u_ll, equations) + v1_rr, v2_rr, v_sound_rr, v_r_rr = cons2speeds(u_rr, equations) + + if (orientation == 1) + v_ll = abs(v1_ll) + v_rr = abs(v1_rr) + else + v_ll = abs(v2_ll) + v_rr = abs(v2_rr) + end + # experimental + return max(v_ll, v_rr) + max(v_sound_ll, v_sound_rr) + max(abs(v_r_ll), abs(v_r_rr)) +end + +### boundary conditions ### + +# adapted from compressible_moist_euler_2d.jl which was probably adapted from compressible_euler_2d.jl +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressibleRainyEulerEquations2D) + 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) + + # name needed variables + rho_v1 = u_local[4] + + # densities + rho_dry_local, rho_moist_local, rho_rain_local, rho_local, rho_inv_local = densities(u_local, + equations) + + # velocities + v_normal = rho_v1 * rho_inv_local + v_sound, gamma = speed_of_sound(u_local, equations) + + # pressure + p_local = pressure(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 + p_star = p_local * + (1 + 0.5f0 * (gamma - 1) * v_normal / v_sound)^(2 * gamma * + inv(gamma - 1)) + else # v_normal > 0.0 + A = 2 / ((gamma + 1) * rho_local) + B = p_local * (gamma - 1) / (gamma + 1) + p_star = p_local + + 0.5f0 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return SVector(0, 0, 0, + p_star * normal[1] * norm_, + p_star * normal[2] * norm_, + 0, 0, 0, 0) +end + +# same as in compressible_euler_2d.jl +@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressibleRainyEulerEquations2D) + # 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 + +# same as in compressible_euler_2d.jl +@inline function rotate_to_x(u, normal_vector, + equations::CompressibleRainyEulerEquations2D) + # 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] + + return SVector(u[1], u[2], u[3], + c * u[4] + s * u[5], + -s * u[4] + c * u[5], + u[6], u[7], u[8], u[9]) +end + +# for parabolic terms (LaplaceDiffusion2D) +@inline function boundary_condition_laplace(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Trixi.Gradient, + equations_parabolic::LaplaceDiffusion2D) + return u_inner +end + +@inline function boundary_condition_laplace(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Trixi.Divergence, + equations_parabolic::LaplaceDiffusion2D) + return flux_inner +end + +### Nonlinear System Residual ### + +# in preparation for a callback to solve the nonlinear system +@inline function saturation_residual(u, guess, + equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + R_v = equations.R_vapour + L_ref = equations.ref_latent_heat_vap_temp + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # velocity + v1, v2 = velocities(u, rho_inv, equations) + + # energy density + energy = energy_density(u, equations) + + # define residual + residual1 = (c_vd * rho_dry + c_vv * guess[1] + c_l * (guess[2] + rho_rain)) * + guess[3] + residual1 += guess[1] * L_ref + residual1 -= (energy - rho * 0.5f0 * (v1^2 + v2^2)) + + residual2 = min(saturation_vapour_pressure(guess[3], equations) / (R_v * guess[3]), + rho_moist) + residual2 -= guess[1] + residual2 *= 1.0f7 + + residual3 = rho_moist + residual3 -= guess[1] + guess[2] + residual3 *= 1.0f7 + + return SVector(residual1, residual2, residual3) +end + +@inline function saturation_residual_jacobian(u, guess, + equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + R_v = equations.R_vapour + L_ref = equations.ref_latent_heat_vap_temp + + # densities + rho_dry, rho_moist, rho_rain, rho, rho_inv = densities(u, equations) + + # saturation + svp = saturation_vapour_pressure(guess[3], equations) + svp_t = saturation_vapour_pressure_derivative(guess[3], equations) + + # define jacobian + J_11 = c_vv * guess[3] + L_ref + J_12 = c_l * guess[3] + J_13 = c_vd * rho_dry + c_vv * guess[1] + c_l * (guess[2] + rho_rain) + + J_21 = -10000000 + J_22 = 0 + + if (svp / (R_v * guess[3]) < rho_moist) + J_23 = (svp_t * guess[3] - svp) / (R_v * guess[3]^2) * 1.0f7 + else + J_23 = 0 + end + + J_31 = -10000000 + J_32 = -10000000 + J_33 = 0 + + return SMatrix{3, 3}(J_11, J_21, J_31, J_12, J_22, J_32, J_13, J_23, J_33) +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. +# adapted from compressible_moist_euler_2d.jl, does NOT work with rain! +@inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleRainyEulerEquations2D) + # constants + a = 360 + + # densities + rho_dry_ll, rho_moist_ll, rho_rain_ll, rho_ll, rho_inv_ll = densities(u_ll, + equations) + rho_dry_rr, rho_moist_rr, rho_rain_rr, rho_rr, rho_inv_rr = densities(u_rr, + equations) + + # pressure + p_ll = pressure(u_ll, equations) + p_rr = pressure(u_rr, equations) + + # velocities + v1_ll, v2_ll = velocities(u_ll, rho_inv_ll, equations) + v1_rr, v2_rr = velocities(u_rr, rho_inv_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] + norm_ = norm(normal_direction) + + # diffusion parameter 0.0 < beta <= 1.0 + beta = 1 + + # interface flux components + 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, _, f4, f5, f6, _, _, _ = u_ll * v_interface + f6 += p_ll * v_interface + else + f1, f2, _, f4, f5, f6, _, _, _ = u_rr * v_interface + f6 += p_rr * v_interface + end + + return SVector(f1, f2, 0, + f4 + p_interface * normal_direction[1], + f5 + p_interface * normal_direction[2], + f6, 0, 0, 0) +end + +""" + flux_ec_rain(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleRainyEulerEquations2D) + +Entropy-conserving flux including rain, derived in: +Fabian Höck +A Discontinuous Galerkin Method for Moist Atmospheric Dynamics with Rain +Master's thesis, University of Cologne, 2025 +""" +@inline function flux_ec_rain(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleRainyEulerEquations2D) + # constants + c_l = equations.c_liquid_water + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + R_d = equations.R_dry_air + R_v = equations.R_vapour + L_ref = equations.ref_latent_heat_vap_temp + + # densities and temperatures + rho_dry_ll, rho_moist_ll, rho_rain_ll, rho_ll, rho_inv_ll = densities(u_ll, + equations) + rho_dry_rr, rho_moist_rr, rho_rain_rr, rho_rr, rho_inv_rr = densities(u_rr, + equations) + rho_vapour_ll, rho_cloud_ll, temperature_ll = cons2nonlinearsystemsol(u_ll, + equations) + rho_vapour_rr, rho_cloud_rr, temperature_rr = cons2nonlinearsystemsol(u_rr, + equations) + inv_temperature_ll = inv(temperature_ll) + inv_temperature_rr = inv(temperature_rr) + + # velocities + v1_ll, v2_ll = velocities(u_ll, rho_inv_ll, equations) + v1_rr, v2_rr = velocities(u_rr, rho_inv_rr, equations) + vr_ll = terminal_velocity_rain(rho_vapour_ll + rho_cloud_ll, rho_rain_ll, equations) + vr_rr = terminal_velocity_rain(rho_vapour_rr + rho_cloud_rr, rho_rain_rr, equations) + + # velocity averages + 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) + K_avg = 0.5f0 * (v1_square_avg + v2_square_avg) + v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg + vr_normal_avg = 0.5f0 * normal_direction[2] * (vr_ll + vr_rr) + + # density averages + rho_dry_avg = 0.5f0 * (rho_dry_ll + rho_dry_rr) + rho_vapour_avg = 0.5f0 * (rho_vapour_ll + rho_vapour_rr) + rho_cloud_avg = 0.5f0 * (rho_cloud_ll + rho_cloud_rr) + rho_rain_avg = 0.5f0 * (rho_rain_ll + rho_rain_rr) + + # density log means + rho_dry_log = 0 + rho_vapour_log = 0 + rho_cloud_log = 0 + rho_rain_log = 0 + + if (!(rho_dry_ll == 0.0) && !(rho_dry_rr == 0.0)) + rho_dry_log = ln_mean(rho_dry_ll, rho_dry_rr) + end + + if (!(rho_vapour_ll == 0.0) && !(rho_vapour_rr == 0.0)) + rho_vapour_log = ln_mean(rho_vapour_ll, rho_vapour_rr) + end + + if (!(rho_cloud_ll == 0.0) && !(rho_cloud_rr == 0.0)) + rho_cloud_log = ln_mean(rho_cloud_ll, rho_cloud_rr) + end + + if (!(rho_rain_ll == 0.0) && !(rho_rain_rr == 0.0)) + rho_rain_log = ln_mean(rho_rain_ll, rho_rain_rr) + end + + # other averages + inv_temperature_avg = 0.5f0 * (inv_temperature_ll + inv_temperature_rr) + inv_temperature_log = inv_ln_mean(inv_temperature_ll, inv_temperature_rr) + p_int = inv(inv_temperature_avg) * (R_d * rho_dry_avg + R_v * rho_vapour_avg) + + # density flux + f_vapour = rho_vapour_log * v_dot_n_avg + f_cloud = rho_cloud_avg * v_dot_n_avg + f_dry = rho_dry_log * v_dot_n_avg + f_rain = rho_rain_avg * (v_dot_n_avg - vr_normal_avg) + f_moist = f_vapour + f_cloud + f_rho = f_dry + f_moist + f_rain + + # momentum flux + f_rhov1 = f_rho * v1_avg + p_int * normal_direction[1] + f_rhov2 = f_rho * v2_avg + p_int * normal_direction[2] + + # energy flux + f_energy = (c_vd * inv_temperature_log - K_avg) * f_dry + + (c_vv * inv_temperature_log - K_avg + L_ref) * f_vapour + + (c_l * inv_temperature_log - K_avg) * (f_cloud + f_rain) + + (v1_avg * f_rhov1 + v2_avg * f_rhov2) + + return SVector(f_dry, f_moist, f_rain, f_rhov1, f_rhov2, f_energy, 0, 0, 0) +end + +@inline function flux_ec_rain(u_ll, u_rr, orientation::Int, + equations::CompressibleRainyEulerEquations2D) + if (orientation == 1) + return flux_ec_rain(u_ll, u_rr, SVector(1, 0), equations) + else + return flux_ec_rain(u_ll, u_rr, SVector(0, 1), equations) + end +end + +# adapted from ShallowWaterEquations2D (Recommended with rain!) +@inline function boundary_condition_simple_slip_wall(u_inner, orientation, direction, x, + t, surface_flux_function, + equations::CompressibleRainyEulerEquations2D) + ## get the appropriate normal vector from the orientation + if orientation == 1 + u_boundary = SVector(u_inner[1], u_inner[2], u_inner[3], -u_inner[4], + u_inner[5], u_inner[6], u_inner[7], u_inner[8], u_inner[9]) + else # orientation == 2 + u_boundary = SVector(u_inner[1], u_inner[2], u_inner[3], u_inner[4], + -u_inner[5], u_inner[6], u_inner[7], u_inner[8], + u_inner[9]) + end + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + + return flux +end + +# hydrostatic base state residual +function generate_hydrostatic_residual(pressure_lower, humidity_rel0, z, dz, + equations::CompressibleRainyEulerEquations2D) + # equations constants + c_pd = equations.c_dry_air_const_pressure + R_d = equations.R_dry_air + R_v = equations.R_vapour + eps = equations.eps + ref_pressure = equations.ref_pressure + g = equations.gravity + + function hydrostatic_residual!(residual, guess) + # variables + pressure, rho_dry, rho_vapour, temperature = guess + + rho_vs = saturation_vapour_pressure(temperature, equations) / + (R_v * temperature) + + # pressure derivative residual approximation + residual[1] = (pressure - pressure_lower) / dz + (rho_dry + rho_vapour) * g + + # pressure residual + residual[2] = pressure - temperature * (rho_dry * R_d + rho_vapour * R_v) + + # hydrostatic dry potential temperature residual + residual[3] = theta_d(z, equations) - + temperature * (ref_pressure / pressure)^(R_d / c_pd) + + # humidity residual + residual[4] = rho_vs * (rho_dry + rho_vapour / eps) * humidity_rel0 + residual[4] -= rho_vapour * (rho_dry + rho_vs / eps) + residual[4] *= 1000 + end + + return hydrostatic_residual! +end + +function generate_perturbation_residual(pressure_hydrostatic, H_init, z, + equations::CompressibleRainyEulerEquations2D) + # equations constants + c_pd = equations.c_dry_air_const_pressure + R_d = equations.R_dry_air + R_v = equations.R_vapour + eps = equations.eps + ref_pressure = equations.ref_pressure + + function perturbation_residual!(residual, guess) + # variables + rho_dry, rho_vapour, temperature = guess + + rho_vs = saturation_vapour_pressure(temperature, equations) / + (R_v * temperature) + pressure = (rho_dry * R_d + rho_vapour * R_v) * temperature + + # humidity residual + residual[1] = rho_vs * (rho_dry + rho_vapour / eps) * H_init + residual[1] -= rho_vapour * (rho_dry + rho_vs / eps) + residual[1] *= 30 + + # hydrostatic dry potential temperature residual + residual[2] = theta_d(z, equations) - + temperature * (ref_pressure / pressure_hydrostatic)^(R_d / c_pd) + + # pressure residual + residual[3] = pressure_hydrostatic - pressure + end + + return perturbation_residual! +end + +# hydrostatic dry potential temperature +function theta_d(z, equations::CompressibleRainyEulerEquations2D{RealT}) where {RealT} + # constants + c_pd = equations.c_dry_air_const_pressure + R_d = equations.R_dry_air + ref_pressure = equations.ref_pressure + + # problem specific constants + surface_temperature = 283 + surface_pressure = 85000 + stratification = convert(RealT, 1.3e-5) + + # dry potential temperature at surface + Theta0 = surface_temperature * (ref_pressure / surface_pressure)^(R_d / c_pd) + # at height z + theta_d = Theta0 * exp(stratification * z) + + return theta_d +end + +# for approximating the dz pressure gradient +struct AtmosphereLayersRainyBubble{RealT <: Real} + layer_data :: Matrix{RealT} + total_height :: RealT + precision :: RealT +end + +function AtmosphereLayersRainyBubble(equations::CompressibleRainyEulerEquations2D; + total_height, precision = 1, RealT = Float64) + # constants + humidity_rel0 = convert(RealT, 0.2) # hydrostatic relative humidity + surface_pressure = 85000 + + # surface layer with initial guesses for rho_dry, rho_vapour and temperature + surface_layer = [surface_pressure, convert(RealT, 1.4), convert(RealT, 0.04), 300] + + # allocate layer_data + n = convert(Int, total_height / precision) + layer_data = zeros(RealT, n + 1, 4) + + # solve (slightly above) surface layer + dz = convert(RealT, 0.01) + z = convert(RealT, 0.01) + residual_function! = generate_hydrostatic_residual(surface_pressure, humidity_rel0, + z, + dz, equations) + layer_data[1, :] .= nlsolve(residual_function!, surface_layer).zero + + # adjust to chosen precision + dz = precision + + # iterate up the atmosphere + for i in (1:n) + z += dz + residual_function! = generate_hydrostatic_residual(layer_data[i, 1], + humidity_rel0, + z, dz, equations) + guess = deepcopy(layer_data[i, :]) + layer_data[i + 1, :] .= nlsolve(residual_function!, guess, ftol = 1e-10, + iterations = 20).zero + end + + return AtmosphereLayersRainyBubble{RealT}(layer_data, total_height, precision) +end + +function initial_condition_bubble_rainy(x, t, + equations::CompressibleRainyEulerEquations2D{RealT}; + atmosphere_layers) where {RealT} + # equations constants + c_vd = equations.c_dry_air_const_volume + c_vv = equations.c_vapour_const_volume + ref_L = equations.ref_latent_heat_vap_temp + + # problem specific constants + humidity_rel_bar = convert(RealT, 0.2) # background relative humidity field + humidity_max = 1 + + # bubble parameters + radius_outer, radius_inner = 300, 200 # radii of humidity bubble + x_center, z_center = 1200, 800 # center of humidity bubble + + # radius relative to bubble center + r = sqrt((x[1] - x_center)^2 + (x[2] - z_center)^2) + + # humidity definition + if (r > radius_outer) + # outside the bubble + humidity = humidity_rel_bar + elseif (r > radius_inner) + # outer layers of the bubble + humidity = humidity_rel_bar + + (humidity_max - humidity_rel_bar) * + cos(pi * (r - radius_inner) / (2 * (radius_outer - radius_inner)))^2 + else + # inner layer + humidity = humidity_max + end + + # get atmosphere layer and height information + @unpack layer_data, total_height, precision = atmosphere_layers + dz = precision + z = x[2] + n = convert(Int, floor((z + eps()) / dz)) + 1 + z_lower = (n - 1) * dz + z_upper = n * dz + + if (z_lower == total_height) + z_upper = z_lower + dz + n = n - 1 + end + + if (n == 0) + n = 1 + end + + # check height consistency + if (z > total_height && !(isapprox(z, total_height))) + error("The atmosphere does not match the simulation domain") + end + + # get hydrostatic pressures and approximate between lower and upper data point + pressure_hydrostatic_lower = layer_data[n, 1] + pressure_hydrostatic_upper = layer_data[n + 1, 1] + pressure_hydrostatic = (pressure_hydrostatic_upper * (z - z_lower) + + pressure_hydrostatic_lower * (z_upper - z)) / dz + + # solve perturbation + residual_function! = generate_perturbation_residual(pressure_hydrostatic, humidity, + z, + equations) + rho_dry, rho_vapour, temperature = nlsolve(residual_function!, layer_data[n, 2:4], + ftol = 1e-9, iterations = 20).zero + + energy_density = (c_vd * rho_dry + c_vv * rho_vapour) * temperature + + rho_vapour * ref_L + + return SVector(rho_dry, rho_vapour, 0, 0, 0, energy_density, rho_vapour, 0, + temperature) +end +end # muladd end diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 17e12bb9..9be009c0 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -330,6 +330,9 @@ abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: abstract type AbstractCovariantShallowWaterEquations2D{GlobalCoordinateSystem} <: AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} end +abstract type AbstractCompressibleRainyEulerEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end + @inline function flux_zero(u_ll, u_rr, normal_direction, equations) return zero(u_ll) end @@ -338,6 +341,7 @@ include("covariant_advection.jl") include("covariant_shallow_water.jl") include("covariant_shallow_water_split.jl") include("compressible_moist_euler_2d_lucas.jl") +include("compressible_rainy_euler_2d.jl") include("compressible_euler_potential_temperature_1d.jl") include("compressible_euler_potential_temperature_2d.jl") include("compressible_euler_potential_temperature_3d.jl") diff --git a/src/meshes/p4est_icosahedron_quads.jl b/src/meshes/p4est_icosahedron_quads.jl index 8cb96118..b72c6c8d 100644 --- a/src/meshes/p4est_icosahedron_quads.jl +++ b/src/meshes/p4est_icosahedron_quads.jl @@ -732,17 +732,17 @@ function calc_node_coordinates_triangle_vertices!(triangle_vertices, triangle_vertices[:, 1] = corners_triangle[:, 1] - v2_bilinear = 0.5 * (corners_triangle[:, 1] + corners_triangle[:, 2]) + v2_bilinear = 0.5f0 * (corners_triangle[:, 1] + corners_triangle[:, 2]) triangle_vertices[:, 2] = radius * v2_bilinear / norm(v2_bilinear) triangle_vertices[:, 3] = corners_triangle[:, 2] - v4_bilinear = 0.5 * (corners_triangle[:, 2] + corners_triangle[:, 3]) + v4_bilinear = 0.5f0 * (corners_triangle[:, 2] + corners_triangle[:, 3]) triangle_vertices[:, 4] = radius * v4_bilinear / norm(v4_bilinear) triangle_vertices[:, 5] = corners_triangle[:, 3] - v6_bilinear = 0.5 * (corners_triangle[:, 1] + corners_triangle[:, 3]) + v6_bilinear = 0.5f0 * (corners_triangle[:, 1] + corners_triangle[:, 3]) triangle_vertices[:, 6] = radius * v6_bilinear / norm(v6_bilinear) v7_bilinear = (corners_triangle[:, 1] + corners_triangle[:, 2] + diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl index 8b0ca624..3190bef8 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl @@ -274,7 +274,7 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr for ii in eachnode(basis) # Multiply derivative_matrix to j-dimension to differentiate wrt η - result += 0.5 * derivative_matrix[j, ii] * + result += 0.5f0 * derivative_matrix[j, ii] * (node_coordinates[m, i, ii, element] * node_coordinates[l, i, ii, element] - node_coordinates[l, i, ii, element] * @@ -318,7 +318,7 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr for ii in eachnode(basis) # Multiply derivative_matrix to i-dimension to differentiate wrt ξ - result += 0.5 * derivative_matrix[i, ii] * + result += 0.5f0 * derivative_matrix[i, ii] * (node_coordinates[m, ii, j, element] * node_coordinates[l, ii, j, element] - node_coordinates[l, ii, j, element] * @@ -336,7 +336,7 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr for ii in eachnode(basis) # Multiply derivative_matrix to i-dimension to differentiate wrt ξ - result += 0.5 * derivative_matrix[i, ii] * + result += 0.5f0 * derivative_matrix[i, ii] * (node_coordinates[m, ii, j, element] * jacobian_matrix[l, 2, ii, j, element] - node_coordinates[l, ii, j, element] * @@ -352,7 +352,7 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr for ii in eachnode(basis) # Multiply derivative_matrix to j-dimension to differentiate wrt η - result += 0.5 * derivative_matrix[j, ii] * + result += 0.5f0 * derivative_matrix[j, ii] * (node_coordinates[m, i, ii, element] * jacobian_matrix[l, 1, i, ii, element] - node_coordinates[l, i, ii, element] * diff --git a/test/runtests.jl b/test/runtests.jl index 033414d0..4d089ed1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,7 +18,8 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 @test true - run(`$(mpiexec()) -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --threads=1 --check-bounds=yes $(abspath("test_mpi.jl"))`) + status = run(ignorestatus(`$(mpiexec()) -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --threads=1 --check-bounds=yes $(abspath("test_mpi.jl"))`)) + @test success(status) end @time if TRIXI_TEST == "all" || TRIXI_TEST == "threaded" @@ -27,7 +28,8 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 @test true - run(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none $(abspath("test_threaded.jl"))`) + status = run(ignorestatus(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none $(abspath("test_threaded.jl"))`)) + @test success(status) end @time if TRIXI_TEST == "all" || TRIXI_TEST == "trixi_consistency" @@ -44,6 +46,7 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXI_TEST == "all" || TRIXI_TEST == "moist_euler" include("test_2d_moist_euler.jl") + include("test_2d_rainy_euler.jl") end @time if TRIXI_TEST == "all" || TRIXI_TEST == "spherical_advection" diff --git a/test/test_2d_rainy_euler.jl b/test/test_2d_rainy_euler.jl new file mode 100644 index 00000000..5639761c --- /dev/null +++ b/test/test_2d_rainy_euler.jl @@ -0,0 +1,159 @@ +module TestExamples2DRainyEuler + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = joinpath(EXAMPLES_DIR, "moist_euler") + +@trixi_testset "convergence_test_rainy" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "convergence_test", + "convergence_test_rainy.jl"), + l2=[ + 2.39895785368954e-6, + 3.892427386570826e-5, + 1.5186141589998252e-5, + 0.007532307120716242, + 0.02085654062623111, + 51.988569260688045, + 1.7470730680079724e-8, + 3.892376665583286e-5, + 0.0002156633044466944 + ], + linf=[ + 1.675032502916618e-5, + 0.0002740479581908595, + 8.691475272920579e-5, + 0.023414703432404593, + 0.06248365429859781, + 385.10196906188503, + 6.003313225070965e-8, + 0.000274026861750043, + 0.0007433951467703537 + ]) + # 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 + +@trixi_testset "elixir_rainy_euler_moist_bubble" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "moist_bubble", + "elixir_rainy_euler_moist_bubble.jl"), + l2=[ + 0.003147716446423741, + 6.295432888771655e-5, + 0.0, + 0.030093419727259596, + 1.2721952635118658, + 1040.525033262258, + 0.0011553981102543557, + 0.0011988489480496496, + 2.079644612712912 + ], + linf=[ + 0.01077582487256823, + 0.00021551948971416587, + 0.0, + 0.17975465277331062, + 2.4701908567341215, + 3726.9826530935243, + 0.003030525617948497, + 0.003225736615268841, + 3.2422474910159167 + ], + cells_per_dimension=(64, 32), + tspan=(0.0, 10.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 + +@trixi_testset "elixir_rainy_euler_rainy_bubble" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "rainy_bubble", + "elixir_rainy_euler_rainy_bubble.jl"), + l2=[ + 7.958084854135629e-5, + 4.8974003075933214e-5, + 1.2699485312744632e-9, + 0.01995521719771058, + 0.03383695227022311, + 126.22998495250405, + 4.8983831176389294e-5, + 3.4524261261367737e-7, + 0.0016458343657569546 + ], + linf=[ + 0.001181534085815561, + 0.0007152546342687649, + 2.065750070675263e-8, + 0.12231969676871683, + 0.176070286893087, + 1831.1175623952004, + 0.0007152546342687649, + 6.114801650025656e-6, + 0.008508910564898997 + ], + cells_per_dimension=(16, 16), + tspan=(0.0, 10.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 + +# For unknown reasons, github's macos runners produce results exceeding the default +# tolerance +@trixi_testset "elixir_rainy_euler_rainy_bubble_diffusion" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "rainy_bubble", + "elixir_rainy_euler_rainy_bubble_diffusion.jl"), + l2=[ + 8.024161547037353e-5, + 4.927749867325615e-5, + 1.1775025026817308e-9, + 0.019941027438108397, + 0.03381791952397996, + 126.96661069448412, + 4.925578049099529e-5, + 3.0245852299389407e-7, + 0.0016287477694402524 + ], + linf=[ + 0.001187049337148749, + 0.0007187571931339255, + 2.67415192177334e-8, + 0.12199924771843951, + 0.1760655719865429, + 1840.1706161463226, + 0.0007187571931339255, + 7.049541733574276e-6, + 0.008742438145986853 + ], + rtol=6e-6, + initial_refinement_level=4, + tspan=(0.0, 10.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 +end # module diff --git a/test/test_mpi.jl b/test/test_mpi.jl index 2cc153aa..c25b9be5 100644 --- a/test/test_mpi.jl +++ b/test/test_mpi.jl @@ -9,13 +9,7 @@ outdir = "out" Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive = true) Trixi.MPI.Barrier(Trixi.mpi_comm()) -# CI with MPI and some tests fails often on Windows. Thus, we check whether this -# is the case here. We use GitHub Actions, so we can check whether we run CI -# in the cloud with Windows as follows, see also -# https://docs.github.com/en/actions/learn-github-actions/environment-variables -CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows() - -@testset "MPI tests" begin +@testset verbose=true showtiming=true "MPI tests" begin #! format: noindent @trixi_testset "elixir_moist_euler_moist_bubble" begin diff --git a/test/test_threaded.jl b/test/test_threaded.jl index 628ab005..4627fd50 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -2,7 +2,7 @@ module TestThreaded include("test_trixiatmo.jl") -@testset "Threaded tests" begin +@testset verbose=true showtiming=true "Threaded tests" begin #! format: noindent @trixi_testset "elixir_moist_euler_moist_bubble" begin diff --git a/test/test_trixiatmo.jl b/test/test_trixiatmo.jl index a3edeeda..e9e06e48 100644 --- a/test/test_trixiatmo.jl +++ b/test/test_trixiatmo.jl @@ -7,6 +7,11 @@ using TrixiAtmo: examples_dir EXAMPLES_DIR = examples_dir() +# Check whether we run CI in the cloud with Windows or Mac, see also +# https://docs.github.com/en/actions/learn-github-actions/environment-variables +CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows() +CI_ON_MAC = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.isapple() + macro test_trixi_include(expr, args...) local add_to_additional_ignore_content = [ # This is needed because we overwrite `Trixi.weak_form_kernel!`, e.g., here: diff --git a/test/test_type.jl b/test/test_type.jl index 52aa644c..e919f9b3 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -10,7 +10,7 @@ include("test_trixiatmo.jl") outdir = "out" isdir(outdir) && rm(outdir, recursive = true) -@testset "Test Type Stability" begin +@testset verbose=true showtiming=true "Test Type Stability" begin @timed_testset "Compressible Euler Potential Temperature 1D" begin for RealT in (Float32, Float64) equations = @inferred CompressibleEulerPotentialTemperatureEquations1D(c_p = RealT(1004), @@ -45,7 +45,7 @@ isdir(outdir) && rm(outdir, recursive = true) for RealT in (Float32, Float64) equations = @inferred CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = RealT(1004), c_v = RealT(717), - gravity = RealT(9.81)) + gravity = RealT(EARTH_GRAVITATIONAL_ACCELERATION)) x = SVector(zero(RealT)) t = zero(RealT) @@ -152,7 +152,7 @@ isdir(outdir) && rm(outdir, recursive = true) for RealT in (Float32, Float64) equations = @inferred CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = RealT(1004), c_v = RealT(717), - gravity = RealT(9.81)) + gravity = RealT(EARTH_GRAVITATIONAL_ACCELERATION)) x = SVector(zero(RealT)) t = zero(RealT) @@ -245,7 +245,7 @@ isdir(outdir) && rm(outdir, recursive = true) for RealT in (Float32, Float64) equations = @inferred CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = RealT(1004), c_v = RealT(717), - gravity = RealT(9.81)) + gravity = RealT(EARTH_GRAVITATIONAL_ACCELERATION)) x = SVector(zero(RealT)) t = zero(RealT) @@ -340,6 +340,84 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Compressible Rainy Euler" begin + using TrixiAtmo: boundary_condition_simple_slip_wall, cons2eq_pot_temp, + cons2nonlinearsystemsol, cons2speeds, densities, velocities, + energy_density, speed_of_sound, terminal_velocity_rain, + saturation_vapour_pressure, saturation_vapour_pressure_derivative, + saturation_residual, saturation_residual_jacobian, theta_d, + AtmosphereLayersRainyBubble + for RealT in (Float32, Float64) + equations = @inferred CompressibleRainyEulerEquations2D(RealT) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(ntuple(i -> one(RealT), 9)) + + normal_direction = SVector(one(RealT), one(RealT)) + surface_flux_function = flux_lax_friedrichs + orientations = [1, 2] + directions = [1, 2] + + for direction in directions, orientation in orientations + @test eltype(@inferred boundary_condition_simple_slip_wall(u, orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end + for direction in directions + @test eltype(@inferred boundary_condition_slip_wall(u, normal_direction, + direction, + x, t, + surface_flux_function, + equations)) == RealT + end + for orientation in orientations + @test eltype(@inferred max_abs_speed(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_ec_rain(u_ll, u_rr, orientation, equations)) == + RealT + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_ec_rain(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred flux_LMARS(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test eltype(@inferred cons2eq_pot_temp(u, equations)) == RealT + @test eltype(@inferred cons2nonlinearsystemsol(u, equations)) == RealT + @test eltype(@inferred cons2speeds(u, equations)) == RealT + @test eltype(@inferred densities(u, equations)) == RealT + @test eltype(@inferred velocities(u, u[1], equations)) == RealT + @test eltype(@inferred energy_density(u, equations)) == RealT + @test eltype(@inferred pressure(u, equations)) == RealT + @test eltype(@inferred speed_of_sound(u, equations)) == RealT + @test eltype(@inferred terminal_velocity_rain(u[2], u[3], equations)) == RealT + @test eltype(@inferred saturation_vapour_pressure(u[9], equations)) == RealT + @test eltype(@inferred saturation_vapour_pressure_derivative(u[9], equations)) == + RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred boundary_condition_slip_wall(u, normal_direction, + x, t, + surface_flux_function, + equations)) == RealT + @test eltype(@inferred saturation_residual(u, u, equations)) == RealT + @test eltype(@inferred saturation_residual_jacobian(u, u, equations)) == RealT + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred theta_d(x[1], equations)) == RealT + + @inferred AtmosphereLayersRainyBubble(equations; total_height = x[1]) + end + end + @timed_testset "Shallow Water 3D" begin for RealT in (Float32, Float64) equations = @inferred ShallowWaterEquations3D(gravity = RealT(1), diff --git a/test/test_unit.jl b/test/test_unit.jl index 8cc426ad..24c2a7a3 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -244,7 +244,7 @@ end # Set up equations and dummy conservative variables state equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) u = SVector(1.1, -0.5, 2.34, 330.0, 1500) normal_directions = [SVector(1.0, 0.0), @@ -264,7 +264,7 @@ end normal_2d = SVector(normal_1d[1], 0.0) equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) equations_2d = equations flux_1d = normal_1d[1] * flux_ec(u_1d, u_1d, 1, equations_1d) flux_2d = flux_ec(u_2d, u_2d, normal_2d, equations_2d) @@ -281,7 +281,7 @@ end # check consistency for 3D EC flux equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1500) normal_directions = [SVector(1.0, 0.0, 0.0), @@ -308,7 +308,7 @@ end # Set up equations and dummy conservative variables state equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) u = SVector(1.1, -0.5, 2.34, 330.0, 1500) normal_directions = [SVector(1.0, 0.0), @@ -328,7 +328,7 @@ end normal_2d = SVector(normal_1d[1], 0.0) equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) equations_2d = equations flux_1d = normal_1d[1] * flux_tec(u_1d, u_1d, 1, equations_1d) flux_2d = flux_tec(u_2d, u_2d, normal_2d, equations_2d) @@ -345,7 +345,7 @@ end # check consistency for 3D EC flux equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1500) normal_directions = [SVector(1.0, 0.0, 0.0), @@ -372,7 +372,7 @@ end # Set up equations and dummy conservative variables state equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) u = SVector(1.1, -0.5, 2.34, 330.0, 1500) normal_directions = [SVector(1.0, 0.0), @@ -392,7 +392,7 @@ end normal_2d = SVector(normal_1d[1], 0.0) equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) equations_2d = equations flux_1d = normal_1d[1] * flux_etec(u_1d, u_1d, 1, equations_1d) flux_2d = flux_etec(u_2d, u_2d, normal_2d, equations_2d) @@ -409,7 +409,7 @@ end # check consistency for 3D EC flux equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1500) normal_directions = [SVector(1.0, 0.0, 0.0), @@ -436,7 +436,7 @@ end # Set up equations and dummy conservative variables state equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) flux_lmars = FluxLMARS(340) u = SVector(1.1, -0.5, 2.34, 330.0, 1700) @@ -457,7 +457,7 @@ end normal_2d = SVector(normal_1d[1], 0.0) equations_1d = CompressibleEulerPotentialTemperatureEquationsWithGravity1D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) equations_2d = equations flux_1d = normal_1d[1] * flux_lmars(u_1d, u_1d, 1, equations_1d) flux_2d = flux_lmars(u_2d, u_2d, normal_2d, equations_2d) @@ -467,7 +467,7 @@ end # check consistency for 3D EC flux equations = CompressibleEulerPotentialTemperatureEquationsWithGravity3D(c_p = 1004.0, c_v = 717.0, - gravity = 9.81) + gravity = EARTH_GRAVITATIONAL_ACCELERATION) u = SVector(1.1, -0.5, 2.34, 2.4, 330.0, 1700) normal_directions = [SVector(1.0, 0.0, 0.0), @@ -523,4 +523,23 @@ end flux(u, aux_vars, orientation, equations) end +@testset "Consistency check for EC flux with Rainy Euler" begin + # Set up equations and dummy conservative variables state + equations = CompressibleRainyEulerEquations2D() + # Example state vector (ρ_d, ρ_m, ρ_r, ρu, ρv, ρe, ρq_v, ρq_c, T) + u = SVector(1.0, 0.2, 0.1, 0.5, -0.4, 2.2, 0.1, 0.1, 300) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + equal = flux_ec_rain(u, u, normal_direction, equations) .≈ + flux(u, normal_direction, equations) + # TODO + expected = [true, true, true, true, true, false, true, true, true] + @test equal == expected + end +end end