From 92001fa7720b88cfdff456f448b82a023ffbdc86 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 17 Jul 2024 11:23:11 +0530 Subject: [PATCH] Formatter --- examples/elixir_moist_euler_EC_bubble.jl | 355 ++++++++-------- examples/elixir_moist_euler_dry_bubble.jl | 107 +++-- examples/elixir_moist_euler_moist_bubble.jl | 392 +++++++++--------- ...oist_euler_nonhydrostatic_gravity_waves.jl | 121 +++--- .../elixir_moist_euler_source_terms_dry.jl | 28 +- .../elixir_moist_euler_source_terms_moist.jl | 29 +- ...ir_moist_euler_source_terms_split_moist.jl | 35 +- 7 files changed, 537 insertions(+), 530 deletions(-) diff --git a/examples/elixir_moist_euler_EC_bubble.jl b/examples/elixir_moist_euler_EC_bubble.jl index e1657ba..0c5b435 100644 --- a/examples/elixir_moist_euler_EC_bubble.jl +++ b/examples/elixir_moist_euler_EC_bubble.jl @@ -10,198 +10,208 @@ using NLsolve: nlsolve equations = CompressibleMoistEulerEquations2D() # TODO - Should the IC functions and struct be in the equation file? -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 +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 +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 AtmossphereLayers{RealT<:Real} - equations::CompressibleMoistEulerEquations2D - # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql - LayerData::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 - equivalentpotential_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. +struct AtmossphereLayers{RealT <: Real} + equations::CompressibleMoistEulerEquations2D + # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql + LayerData::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 + equivalentpotential_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 AtmossphereLayers(equations; total_height = 10010.0, preciseness = 10, + ground_state = (1.4, 100000.0), + equivalentpotential_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 = equivalentpotential_temperature -function AtmossphereLayers(equations ; total_height=10010.0, preciseness=10, ground_state=(1.4, 100000.0), equivalentpotential_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 = equivalentpotential_temperature - - rho_qv0 = rho0 * r_v0 - T0 = theta_e0 - y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] - - n = convert(Int, total_height / preciseness) - dz = 0.01 - LayerData = 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_qv0 = rho0 * r_v0 + T0 = theta_e0 + y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] - 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) + n = convert(Int, total_height / preciseness) + dz = 0.01 + LayerData = zeros(RealT, n + 1, 4) - LayerData[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) - - LayerData[i+1, :] = [rho, rho_theta, rho_qv, rho_ql] - end + kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + LayerData[1, :] = [rho, rho_theta, rho_qv, rho_ql] + for i in (1:n) + 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) + + LayerData[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] + end - return AtmossphereLayers{RealT}(equations, LayerData, total_height, dz, n, ground_state, theta_e0, mixing_ratios) + return AtmossphereLayers{RealT}(equations, LayerData, total_height, dz, n, ground_state, + theta_e0, mixing_ratios) end # Generate background state from the Layer data set by linearely interpolating the layers -function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D, AtmosphereLayers::AtmossphereLayers) - @unpack LayerData, preciseness, total_height = AtmosphereLayers - dz = preciseness - z = x[2] - if (z > total_height && !(isapprox(z, total_height))) - error("The atmossphere does not match the simulation domain") - end - n = convert(Int, floor(z/dz)) + 1 - z_l = (n-1) * dz - (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :] - z_r = n * dz - if (z_l == total_height) - z_r = z_l + dz - n = n-1 - end - (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = LayerData[n+1, :] - rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz - rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / dz - rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz - rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz - - rho, rho_e, rho_qv, rho_ql = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D) - - v1 = 60.0 - v2 = 60.0 - rho_v1 = rho * v1 - rho_v2 = rho * v2 - rho_E = rho_e + 1/2 * rho *(v1^2 + v2^2) - - - return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) +function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D, + AtmosphereLayers::AtmossphereLayers) + @unpack LayerData, preciseness, total_height = AtmosphereLayers + dz = preciseness + z = x[2] + if (z > total_height && !(isapprox(z, total_height))) + error("The atmossphere does not match the simulation domain") + end + n = convert(Int, floor(z / dz)) + 1 + z_l = (n - 1) * dz + (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :] + z_r = n * dz + if (z_l == total_height) + z_r = z_l + dz + n = n - 1 + end + (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = LayerData[n + 1, :] + rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz + rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / + dz + rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz + rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz + + rho, rho_e, rho_qv, rho_ql = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, + equations::CompressibleMoistEulerEquations2D) + + v1 = 60.0 + v2 = 60.0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2) + + return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) end # Add perturbation to the profile -function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D) - @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 - rc = 2000 - # Peak perturbation at center - Δθ = 10 - - 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) - θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) - θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5*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 - # SaturVapor - pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) - rho_d_new = (p_loc - pvs) / (R_d * T_loc) - rvs = pvs / (R_v * rho_d_new * T_loc) - θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) - if abs(θ_new-θ_loc) <= θ_loc * 1.0e-12 - break - else - θ_loc=θ_new - end - end - else - rvs = 0 - T_loc = θ_loc * (p_loc / p_0)^kappa - rho_d_new = p_loc / (R_d * T_loc) - θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) - end - rho_qv = rvs * rho_d_new - rho_ql = (rt - rvs) * rho_d_new - rho = rho_d_new * (1 + rt) +function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, + equations::CompressibleMoistEulerEquations2D) + @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 + rc = 2000 + # Peak perturbation at center + Δθ = 10 + + 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) - rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M) + 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 - end - return SVector(rho, rho_e, 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) + 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 + # SaturVapor + pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + rho_d_new = (p_loc - pvs) / (R_d * T_loc) + rvs = pvs / (R_v * rho_d_new * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + if abs(θ_new - θ_loc) <= θ_loc * 1.0e-12 + break + else + θ_loc = θ_new + end + end + else + rvs = 0 + T_loc = θ_loc * (p_loc / p_0)^kappa + rho_d_new = p_loc / (R_d * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + end + rho_qv = rvs * rho_d_new + rho_ql = (rt - rvs) * rho_d_new + rho = rho_d_new * (1 + rt) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + end + + return SVector(rho, rho_e, rho_qv, rho_ql) end AtmossphereData = AtmossphereLayers(equations) function initial_condition_moist(x, t, equations) - return initial_condition_moist_bubble(x, t, equations, AtmossphereData) + return initial_condition_moist_bubble(x, t, equations, AtmossphereData) end initial_condition = initial_condition_moist -boundary_condition = (x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, - y_neg=boundary_condition_periodic, - y_pos=boundary_condition_periodic) +boundary_condition = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) source_term = source_terms_geopotential @@ -211,16 +221,13 @@ source_term = source_terms_geopotential polydeg = 4 basis = LobattoLegendreBasis(polydeg) - surface_flux = flux_chandrashekar volume_flux = flux_chandrashekar - -volume_integral=VolumeIntegralFluxDifferencing(volume_flux) +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) - coordinates_min = (0.0, 0.0) coordinates_max = (4000.0, 4000.0) @@ -233,8 +240,8 @@ mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) # create the semi discretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - boundary_conditions=boundary_condition, - source_terms=source_term) + boundary_conditions = boundary_condition, + source_terms = source_term) ############################################################################### # ODE solvers, callbacks etc. @@ -248,16 +255,17 @@ summary_callback = SummaryCallback() analysis_interval = 1000 solution_variables = cons2aeqpot -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, extra_analysis_errors=(:entropy_conservation_error, )) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) -alive_callback = AliveCallback(analysis_interval=analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval=1000, - save_initial_solution=true, - save_final_solution=true, - solution_variables=solution_variables) +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = solution_variables) -stepsize_callback = StepsizeCallback(cfl=0.8) +stepsize_callback = StepsizeCallback(cfl = 0.8) callbacks = CallbackSet(summary_callback, analysis_callback, @@ -268,9 +276,8 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation - -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl index 7c10032..3d68bd4 100644 --- a/examples/elixir_moist_euler_dry_bubble.jl +++ b/examples/elixir_moist_euler_dry_bubble.jl @@ -13,45 +13,45 @@ equations = CompressibleMoistEulerEquations2D() # for the elastic equations incorporating second-order Runge–Kutta # time differencing. Mon. Wea. Rev., 126, 1992–1999. function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEquations2D) - @unpack p_0, kappa, g, c_pd, c_vd, R_d, R_v = equations - xc = 10000.0 - zc = 2000.0 - r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) - rc = 2000.0 - θ_ref = 300.0 - Δθ = 0.0 - - if r <= rc - Δθ = 2 * cospi(0.5*r/rc)^2 - end - - #Perturbed state: - θ = θ_ref + Δθ # potential temperature - # π_exner = 1 - g / (c_pd * θ) * x[2] # exner pressure - # rho = p_0 / (R_d * θ) * (π_exner)^(c_vd / R_d) # density - - # calculate background pressure with assumption hydrostatic and neutral - p = p_0 * (1-kappa * g * x[2] / (R_d * θ_ref))^(c_pd / R_d) - - #calculate rho and T with p and theta (now perturbed) rho = p / R_d T, T = θ / π - rho = p / ((p / p_0)^kappa*R_d*θ) - T = p / (R_d * rho) - - v1 = 20.0 - #v1 = 0.0 - v2 = 0.0 - rho_v1 = rho * v1 - rho_v2 = rho * v2 - rho_E = rho * c_vd * T + 1/2 * rho * (v1^2 + v2^2) - return SVector(rho, rho_v1, rho_v2, rho_E, zero(eltype(g)) ,zero(eltype(g))) + @unpack p_0, kappa, g, c_pd, c_vd, R_d, R_v = equations + xc = 10000.0 + zc = 2000.0 + r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) + rc = 2000.0 + θ_ref = 300.0 + Δθ = 0.0 + + if r <= rc + Δθ = 2 * cospi(0.5 * r / rc)^2 + end + + #Perturbed state: + θ = θ_ref + Δθ # potential temperature + # π_exner = 1 - g / (c_pd * θ) * x[2] # exner pressure + # rho = p_0 / (R_d * θ) * (π_exner)^(c_vd / R_d) # density + + # calculate background pressure with assumption hydrostatic and neutral + p = p_0 * (1 - kappa * g * x[2] / (R_d * θ_ref))^(c_pd / R_d) + + #calculate rho and T with p and theta (now perturbed) rho = p / R_d T, T = θ / π + rho = p / ((p / p_0)^kappa * R_d * θ) + T = p / (R_d * rho) + + v1 = 20.0 + #v1 = 0.0 + v2 = 0.0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_E = rho * c_vd * T + 1 / 2 * rho * (v1^2 + v2^2) + return SVector(rho, rho_v1, rho_v2, rho_E, zero(eltype(g)), zero(eltype(g))) end initial_condition = initial_condition_warm_bubble -boundary_condition = (x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, - y_neg=boundary_condition_slip_wall, - y_pos=boundary_condition_slip_wall) +boundary_condition = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) # Gravity source since Q_ph=0 source_term = source_terms_geopotential @@ -64,8 +64,7 @@ basis = LobattoLegendreBasis(polydeg) surface_flux = flux_LMARS volume_flux = flux_chandrashekar -volume_integral=VolumeIntegralFluxDifferencing(volume_flux) - +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 @@ -74,23 +73,22 @@ solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (20000.0, 10000.0) - cells_per_dimension = (64, 32) # Create curved mesh with 64 x 32 elements -mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, periodicity = (true, false)) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - boundary_conditions=boundary_condition, - source_terms=source_term) + boundary_conditions = boundary_condition, + source_terms = source_term) ############################################################################### # ODE solvers, callbacks etc. tspan = (0.0, 1000.0) - # Create ODE problem with time span from 0.0 to 1000.0 ode = semidiscretize(semi, tspan) @@ -99,17 +97,17 @@ summary_callback = SummaryCallback() analysis_interval = 1000 solution_variables = cons2drypot -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, extra_analysis_errors=(:entropy_conservation_error, )) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) -alive_callback = AliveCallback(analysis_interval=analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval=1000, - save_initial_solution=true, - save_final_solution=true, - solution_variables=solution_variables) +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = solution_variables) - -stepsize_callback = StepsizeCallback(cfl=0.2) +stepsize_callback = StepsizeCallback(cfl = 0.2) callbacks = CallbackSet(summary_callback, analysis_callback, @@ -117,11 +115,10 @@ callbacks = CallbackSet(summary_callback, save_solution, stepsize_callback) - ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - maxiters=1.0e7, - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl index 281811e..4a34160 100644 --- a/examples/elixir_moist_euler_moist_bubble.jl +++ b/examples/elixir_moist_euler_moist_bubble.jl @@ -9,203 +9,214 @@ using NLsolve: nlsolve equations = CompressibleMoistEulerEquations2D() -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 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] -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 + 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 -struct AtmossphereLayers{RealT<:Real} - equations::CompressibleMoistEulerEquations2D - # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql - LayerData::Matrix{RealT} - total_height::RealT - preciseness::Int - layers::Int - ground_state::NTuple{2, RealT} - equivalentpotential_temperature::RealT - mixing_ratios::NTuple{2, RealT} +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 AtmossphereLayers{RealT <: Real} + equations::CompressibleMoistEulerEquations2D + # structure: 1--> i-layer (z = total_height/precision *(i-1)), 2--> rho, rho_theta, rho_qv, rho_ql + LayerData::Matrix{RealT} + total_height::RealT + preciseness::Int + layers::Int + ground_state::NTuple{2, RealT} + equivalentpotential_temperature::RealT + mixing_ratios::NTuple{2, RealT} +end -function AtmossphereLayers(equations ; total_height=10010.0, preciseness=10, ground_state=(1.4, 100000.0), equivalentpotential_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 = equivalentpotential_temperature - - rho_qv0 = rho0 * r_v0 - T0 = theta_e0 - y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] - - n = convert(Int, total_height / preciseness) - dz = 0.01 - LayerData = zeros(RealT, n+1, 4) +function AtmossphereLayers(equations; total_height = 10010.0, preciseness = 10, + ground_state = (1.4, 100000.0), + equivalentpotential_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 = equivalentpotential_temperature - 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_qv0 = rho0 * r_v0 + T0 = theta_e0 + y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0] - 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) + n = convert(Int, total_height / preciseness) + dz = 0.01 + LayerData = zeros(RealT, n + 1, 4) - LayerData[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) - - LayerData[i+1, :] = [rho, rho_theta, rho_qv, rho_ql] - end + kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) * r_v) / (1 + r_t) + + LayerData[1, :] = [rho, rho_theta, rho_qv, rho_ql] + for i in (1:n) + 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) + + LayerData[i + 1, :] = [rho, rho_theta, rho_qv, rho_ql] + end - return AtmossphereLayers{RealT}(equations, LayerData, total_height, dz, n, ground_state, theta_e0, mixing_ratios) + return AtmossphereLayers{RealT}(equations, LayerData, total_height, dz, n, ground_state, + theta_e0, mixing_ratios) end # Moist bubble test case from paper: # G.H. Bryan, J.M. Fritsch, A Benchmark Simulation for Moist Nonhydrostatic Numerical # Models, MonthlyWeather Review Vol.130, 2917–2928, 2002, # https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml. -function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D, AtmosphereLayers::AtmossphereLayers) - @unpack LayerData, preciseness, total_height = AtmosphereLayers - dz = preciseness - z = x[2] - if (z > total_height && !(isapprox(z, total_height))) - error("The atmossphere does not match the simulation domain") - end - n = convert(Int, floor((z+eps())/dz)) + 1 - z_l = (n-1) * dz - (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :] - z_r = n * dz - if (z_l == total_height) - z_r = z_l + dz - n = n-1 - end - (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = LayerData[n+1, :] - rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz - rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / dz - rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz - rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz - - rho, rho_e, rho_qv, rho_ql = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D) - - v1 = 0.0 - v2 = 0.0 - rho_v1 = rho * v1 - rho_v2 = rho * v2 - rho_E = rho_e + 1/2 * rho *(v1^2 + v2^2) - - - return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) +function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D, + AtmosphereLayers::AtmossphereLayers) + @unpack LayerData, preciseness, total_height = AtmosphereLayers + dz = preciseness + z = x[2] + if (z > total_height && !(isapprox(z, total_height))) + error("The atmossphere does not match the simulation domain") + end + n = convert(Int, floor((z + eps()) / dz)) + 1 + z_l = (n - 1) * dz + (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :] + z_r = n * dz + if (z_l == total_height) + z_r = z_l + dz + n = n - 1 + end + (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = LayerData[n + 1, :] + rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz + rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / + dz + rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz + rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz + + rho, rho_e, rho_qv, rho_ql = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, + equations::CompressibleMoistEulerEquations2D) + + v1 = 0.0 + v2 = 0.0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_E = rho_e + 1 / 2 * rho * (v1^2 + v2^2) + + return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) end -function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D) - @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 - - r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2) - rho_d = rho - rho_qv - rho_ql - kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) - p_loc = p_0 *(R_d * rho_theta / p_0)^(1/(1-kappa_M)) - T_loc = p_loc / (R_d * rho_d + R_v * rho_qv) - rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv - - p_v = rho_qv * R_v * T_loc - p_d = p_loc - p_v - T_C = T_loc - 273.15 - p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) - H = p_v / p_vs - r_v = rho_qv / rho_d - r_l = rho_ql / rho_d - r_t = r_v + r_l - - # equivalentpotential temperature - a=T_loc * (p_0 / p_d)^(R_d / (c_pd + r_t * c_pl)) - b=H^(- r_v * R_v /c_pd) - L_v = L_00 + (c_pv - c_pl) * T_loc - c=exp(L_v * r_v / ((c_pd + r_t * c_pl) * T_loc)) - aeq_pot = (a * b *c) - - # Assume pressure stays constant - if (r < rc && Δθ > 0) - # Calculate background density potential temperature - θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) - # Calculate perturbed density potential temperature - θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5*r/rc)^2 / 300) - rt =(rho_qv + rho_ql) / rho_d - rv = rho_qv / rho_d - # Calculate moist potential temperature - θ_loc = θ_dens_new * (1 + rt)/(1 + (R_v / R_d) * rv) - # Adjust varuables until the temperature is met - if rt > 0 - while true - T_loc = θ_loc * (p_loc / p_0)^kappa - T_C = T_loc - 273.15 - # SaturVapor - pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) - rho_d_new = (p_loc - pvs) / (R_d * T_loc) - rvs = pvs / (R_v * rho_d_new * T_loc) - θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) - if abs(θ_new-θ_loc) <= θ_loc * 1.0e-12 - break - else - θ_loc=θ_new - end - end - else - rvs = 0 - T_loc = θ_loc * (p_loc / p_0)^kappa - rho_d_new = p_loc / (R_d * T_loc) - θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) - end - rho_qv = rvs * rho_d_new - rho_ql = (rt - rvs) * rho_d_new - rho = rho_d_new * (1 + rt) +function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, + equations::CompressibleMoistEulerEquations2D) + @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 + + 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) - rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M) + 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 - end - return SVector(rho, rho_e, 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)) + H = p_v / p_vs + r_v = rho_qv / rho_d + r_l = rho_ql / rho_d + r_t = r_v + r_l + + # equivalentpotential temperature + a = T_loc * (p_0 / p_d)^(R_d / (c_pd + r_t * c_pl)) + b = H^(-r_v * R_v / c_pd) + L_v = L_00 + (c_pv - c_pl) * T_loc + c = exp(L_v * r_v / ((c_pd + r_t * c_pl) * T_loc)) + aeq_pot = (a * b * c) + + # Assume pressure stays constant + if (r < rc && Δθ > 0) + # Calculate background density potential temperature + θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa) + # Calculate perturbed density potential temperature + θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5 * r / rc)^2 / 300) + rt = (rho_qv + rho_ql) / rho_d + rv = rho_qv / rho_d + # Calculate moist potential temperature + θ_loc = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rv) + # Adjust varuables until the temperature is met + if rt > 0 + while true + T_loc = θ_loc * (p_loc / p_0)^kappa + T_C = T_loc - 273.15 + # SaturVapor + pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + rho_d_new = (p_loc - pvs) / (R_d * T_loc) + rvs = pvs / (R_v * rho_d_new * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + if abs(θ_new - θ_loc) <= θ_loc * 1.0e-12 + break + else + θ_loc = θ_new + end + end + else + rvs = 0 + T_loc = θ_loc * (p_loc / p_0)^kappa + rho_d_new = p_loc / (R_d * T_loc) + θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs) + end + rho_qv = rvs * rho_d_new + rho_ql = (rt - rvs) * rho_d_new + rho = rho_d_new * (1 + rt) + rho_d = rho - rho_qv - rho_ql + kappa_M = (R_d * rho_d + R_v * rho_qv) / + (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql) + rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M) + rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv + end + return SVector(rho, rho_e, rho_qv, rho_ql) end # Create background atmosphere data set @@ -213,15 +224,15 @@ AtmossphereData = AtmossphereLayers(equations) # Create the initial condition with the initial data set function initial_condition_moist(x, t, equations) - return initial_condition_moist_bubble(x, t, equations, AtmossphereData) + return initial_condition_moist_bubble(x, t, equations, AtmossphereData) end initial_condition = initial_condition_moist -boundary_condition = (x_neg=boundary_condition_slip_wall, - x_pos=boundary_condition_slip_wall, - y_neg=boundary_condition_slip_wall, - y_pos=boundary_condition_slip_wall) +boundary_condition = (x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) source_term = source_terms_moist_bubble @@ -234,12 +245,10 @@ basis = LobattoLegendreBasis(polydeg) surface_flux = flux_LMARS volume_flux = flux_chandrashekar - -volume_integral=VolumeIntegralFluxDifferencing(volume_flux) +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) - coordinates_min = (0.0, 0.0) coordinates_max = (20000.0, 10000.0) @@ -253,8 +262,8 @@ mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, # create the semi discretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - boundary_conditions=boundary_condition, - source_terms=source_term) + boundary_conditions = boundary_condition, + source_terms = source_term) ############################################################################### # ODE solvers, callbacks etc. @@ -267,16 +276,19 @@ summary_callback = SummaryCallback() analysis_interval = 1000 solution_variables = cons2aeqpot -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, extra_analysis_errors=(:entropy_conservation_error, ), extra_analysis_integrals=(entropy, energy_total, saturation_pressure)) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,), + extra_analysis_integrals = (entropy, energy_total, + saturation_pressure)) -alive_callback = AliveCallback(analysis_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=solution_variables) +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = solution_variables) -stepsize_callback = StepsizeCallback(cfl=0.2) +stepsize_callback = StepsizeCallback(cfl = 0.2) callbacks = CallbackSet(summary_callback, analysis_callback, @@ -287,7 +299,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl index 8d0fc12..b3b87c3 100644 --- a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl +++ b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl @@ -1,7 +1,7 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential, - source_terms_phase_change, + source_terms_phase_change, source_terms_nonhydrostatic_raylight_sponge, cons2drypot, flux_LMARS @@ -12,26 +12,27 @@ using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential, # W. A. Gallus JR., J. B. Klemp, Behavior of Flow over Step Orography, Monthly Weather # 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) - @unpack p_0, kappa, gamma, g, c_pd, c_vd, R_d, R_v = equations - z = x[2] - T_0 = 280.0 - theta_0 = T_0 - N = 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) - rho = p / ((p / p_0)^kappa*R_d*theta) - T = p / (R_d * rho) - - v1 = 10 - v2 = 0 - rho_v1 = rho * v1 - rho_v2 = rho * v2 - rho_E = rho * c_vd * T + 0.5 * rho * (v1^2 + v2^2) - rho_qv = 0 - rho_ql = 0 - return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) +function initial_condition_nonhydrostatic_gravity_wave(x, t, + equations::CompressibleMoistEulerEquations2D) + @unpack p_0, kappa, gamma, g, c_pd, c_vd, R_d, R_v = equations + z = x[2] + T_0 = 280.0 + theta_0 = T_0 + N = 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) + rho = p / ((p / p_0)^kappa * R_d * theta) + T = p / (R_d * rho) + + v1 = 10 + v2 = 0 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_E = rho * c_vd * T + 0.5 * rho * (v1^2 + v2^2) + rho_qv = 0 + rho_ql = 0 + return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) end equations = CompressibleMoistEulerEquations2D() @@ -39,60 +40,59 @@ equations = CompressibleMoistEulerEquations2D() initial_condition = initial_condition_nonhydrostatic_gravity_wave function source(u, x, t, equations::CompressibleMoistEulerEquations2D) - return (source_terms_geopotential(u, equations) + - source_terms_phase_change(u, equations::CompressibleMoistEulerEquations2D) + - source_terms_nonhydrostatic_raylight_sponge(u, x, t, equations::CompressibleMoistEulerEquations2D)) + return (source_terms_geopotential(u, equations) + + source_terms_phase_change(u, equations::CompressibleMoistEulerEquations2D) + + source_terms_nonhydrostatic_raylight_sponge(u, x, t, + equations::CompressibleMoistEulerEquations2D)) end -source_term=source +source_term = source -boundary_conditions = ( - x_neg=boundary_condition_periodic, - x_pos=boundary_condition_periodic, - y_neg=boundary_condition_slip_wall, - y_pos=boundary_condition_slip_wall, - ) +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) polydeg = 4 basis = LobattoLegendreBasis(polydeg) surface_flux = flux_LMARS volume_flux = flux_chandrashekar - -volume_integral=VolumeIntegralFluxDifferencing(volume_flux) - -solver = DGSEM(basis, surface_flux, volume_integral) +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +solver = DGSEM(basis, 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 - # 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))) + h = 400.0 + a = 1000.0 + x_length = 40000.0 + # 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) -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))) -f4(s) = SVector( 20000.0 * s, 16000.0) - +f2(s) = SVector(20000.0, 8000.0 * s + 8000.0) +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))) +end +f4(s) = SVector(20000.0 * s, 16000.0) faces = (f1, f2, f3, f4) - # dx = 0.2*a dz = 10-200 m für (40,16) km cells_per_dimension = (100, 80) - mesh = StructuredMesh(cells_per_dimension, faces, periodicity = (true, false)) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms=source_term, boundary_conditions=boundary_conditions) + source_terms = source_term, + boundary_conditions = boundary_conditions) ############################################################################### # ODE solvers, callbacks etc. @@ -105,17 +105,16 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() analysis_interval = 10000 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval) - -alive_callback = AliveCallback(analysis_interval=analysis_interval) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -save_solution = SaveSolutionCallback(interval=5000, - save_initial_solution=true, - save_final_solution=true, - solution_variables=cons2drypot) +alive_callback = AliveCallback(analysis_interval = analysis_interval) +save_solution = SaveSolutionCallback(interval = 5000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2drypot) -stepsize_callback = StepsizeCallback(cfl=0.2) +stepsize_callback = StepsizeCallback(cfl = 0.2) callbacks = CallbackSet(summary_callback, analysis_callback, @@ -123,14 +122,12 @@ callbacks = CallbackSet(summary_callback, save_solution, stepsize_callback); - ############################################################################### # run the simulation - sol = solve(ode, SSPRK33(), - maxiters=1.0e7, - dt=1, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); + maxiters = 1.0e7, + dt = 1, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/examples/elixir_moist_euler_source_terms_dry.jl b/examples/elixir_moist_euler_source_terms_dry.jl index 9d88ca8..efe7a52 100644 --- a/examples/elixir_moist_euler_source_terms_dry.jl +++ b/examples/elixir_moist_euler_source_terms_dry.jl @@ -8,11 +8,11 @@ using TrixiAtmo: source_terms_convergence_test_dry, initial_condition_convergenc ############################################################################### # semidiscretization of the compressible Euler equations -equations = CompressibleMoistEulerEquations2D( ;g=0.0) +equations = CompressibleMoistEulerEquations2D(; g = 0.0) initial_condition = initial_condition_convergence_test_dry -solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) coordinates_min = (0.0, 0.0) coordinates_max = (2.0, 2.0) @@ -21,10 +21,8 @@ cells_per_dimension = (16, 16) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms=source_terms_convergence_test_dry) - + source_terms = source_terms_convergence_test_dry) ############################################################################### # ODE solvers, callbacks etc. @@ -35,16 +33,16 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -alive_callback = AliveCallback(analysis_interval=analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval=100, - save_initial_solution=true, - save_final_solution=true, - solution_variables=cons2prim) +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl=1.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, @@ -54,7 +52,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/examples/elixir_moist_euler_source_terms_moist.jl b/examples/elixir_moist_euler_source_terms_moist.jl index 08d9e39..f5550d7 100644 --- a/examples/elixir_moist_euler_source_terms_moist.jl +++ b/examples/elixir_moist_euler_source_terms_moist.jl @@ -3,12 +3,13 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo -using TrixiAtmo: source_terms_convergence_test_moist, initial_condition_convergence_test_moist +using TrixiAtmo: source_terms_convergence_test_moist, + initial_condition_convergence_test_moist ############################################################################### # semidiscretization of the compressible Euler equations -equations = CompressibleMoistEulerEquations2D() +equations = CompressibleMoistEulerEquations2D() initial_condition = initial_condition_convergence_test_moist @@ -26,10 +27,8 @@ cells_per_dimension = (4, 4) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms=source_terms_convergence_test_moist) - + source_terms = source_terms_convergence_test_moist) ############################################################################### # ODE solvers, callbacks etc. @@ -40,16 +39,16 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -alive_callback = AliveCallback(analysis_interval=analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval=100, - save_initial_solution=true, - save_final_solution=true, - solution_variables=cons2prim) +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl=1.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, @@ -59,7 +58,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/examples/elixir_moist_euler_source_terms_split_moist.jl b/examples/elixir_moist_euler_source_terms_split_moist.jl index 73bb0ae..f0e2415 100644 --- a/examples/elixir_moist_euler_source_terms_split_moist.jl +++ b/examples/elixir_moist_euler_source_terms_split_moist.jl @@ -3,24 +3,23 @@ using OrdinaryDiffEq using Trixi, TrixiAtmo -using TrixiAtmo: source_terms_convergence_test_moist, initial_condition_convergence_test_moist +using TrixiAtmo: source_terms_convergence_test_moist, + initial_condition_convergence_test_moist ############################################################################### # semidiscretization of the compressible Euler equations -equations = CompressibleMoistEulerEquations2D() +equations = CompressibleMoistEulerEquations2D() 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) +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) @@ -31,10 +30,8 @@ cells_per_dimension = (4, 4) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms=source_terms_convergence_test_moist) - + source_terms = source_terms_convergence_test_moist) ############################################################################### # ODE solvers, callbacks etc. @@ -45,16 +42,16 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -alive_callback = AliveCallback(analysis_interval=analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval=100, - save_initial_solution=true, - save_final_solution=true, - solution_variables=cons2prim) +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl=1.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, @@ -64,8 +61,8 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - maxiters=1.0e7, - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary