Merge branch 'tm/test_spherical_advection' into tm/covariant_form

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0-DEV"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
@@ -17,8 +18,10 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 [compat]
 LinearAlgebra = "1"
 MuladdMacro = "0.2.2"
-Static = "0.8"
-StaticArrayInterface = "1"
+NLsolve = "4.5.1"
+Reexport = "1.0"
+Static = "0.8, 1"
+StaticArrayInterface = "1.5.1"
 StaticArrays = "1"
 StrideArrays = "0.1.28"
 Trixi = "0.7, 0.8"

examples/elixir_euler_spherical_advection_cartesian.jl

@@ -6,6 +6,8 @@ using TrixiAtmo
 ###############################################################################
 # semidiscretization of the linear advection equation
 
+# We use the Euler equations structure but modify the rhs! function to convert it to a
+# variable-coefficient advection equation
 equations = CompressibleEulerEquations3D(1.4)
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux

examples/elixir_moist_euler_EC_bubble.jl

@@ -9,7 +9,6 @@ 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
@@ -46,35 +45,36 @@ function generate_function_of_y(dz, y0, r_t0, theta_e0,
         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}
+struct AtmosphereLayers{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
+    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
-    equivalentpotential_temperature::RealT  # Value for theta_e since we have a constant temperature theta_e0=theta_e.
+    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 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)
+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 = equivalentpotential_temperature
+    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
-    LayerData = zeros(RealT, n + 1, 4)
+    layer_data = zeros(RealT, n + 1, 4)
 
     F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations)
     sol = nlsolve(F, y0)
@@ -85,7 +85,7 @@ function AtmossphereLayers(equations; total_height = 10010.0, preciseness = 10,
     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]
+    layer_data[1, :] = [rho, rho_theta, rho_qv, rho_ql]
     for i in (1:n)
         y0 = deepcopy(sol.zero)
         dz = preciseness
@@ -99,39 +99,39 @@ function AtmossphereLayers(equations; total_height = 10010.0, preciseness = 10,
                   (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]
+        layer_data[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 AtmosphereLayers{RealT}(equations, layer_data, 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
+                                        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 atmossphere does not match the simulation domain")
+        error("The atmosphere 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, :]
+    (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) = LayerData[n + 1, :]
+    (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 = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql,
-                                                     equations::CompressibleMoistEulerEquations2D)
+    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
@@ -143,8 +143,8 @@ function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerE
 end
 
 # Add perturbation to the profile
-function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql,
-                             equations::CompressibleMoistEulerEquations2D)
+function perturb_moist_profile!(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
@@ -200,10 +200,10 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql,
     return SVector(rho, rho_e, rho_qv, rho_ql)
 end
 
-AtmossphereData = AtmossphereLayers(equations)
+atmosphere_data = AtmosphereLayers(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, atmosphere_data)
 end
 
 initial_condition = initial_condition_moist

examples/elixir_moist_euler_dry_bubble.jl

@@ -1,5 +1,5 @@
 using OrdinaryDiffEq
-using Trixi # TODO - Decide. This structure requires Trixi.jl to be in Project.toml of `Test`
+using Trixi
 using TrixiAtmo
 using TrixiAtmo: flux_LMARS, source_terms_geopotential, cons2drypot
 

examples/elixir_moist_euler_moist_bubble.jl

@@ -46,34 +46,34 @@ function generate_function_of_y(dz, y0, r_t0, theta_e0,
     end
 end
 
-struct AtmossphereLayers{RealT <: Real}
+struct AtmosphereLayers{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}
+    layer_data::Matrix{RealT}
     total_height::RealT
     preciseness::Int
     layers::Int
     ground_state::NTuple{2, RealT}
-    equivalentpotential_temperature::RealT
+    equivalent_potential_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)
+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 = equivalentpotential_temperature
+    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
-    LayerData = zeros(RealT, n + 1, 4)
+    layer_data = zeros(RealT, n + 1, 4)
 
     F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations)
     sol = nlsolve(F, y0)
@@ -84,7 +84,7 @@ function AtmossphereLayers(equations; total_height = 10010.0, preciseness = 10,
     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]
+    layer_data[1, :] = [rho, rho_theta, rho_qv, rho_ql]
     for i in (1:n)
         y0 = deepcopy(sol.zero)
         dz = preciseness
@@ -98,42 +98,42 @@ function AtmossphereLayers(equations; total_height = 10010.0, preciseness = 10,
                   (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]
+        layer_data[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 AtmosphereLayers{RealT}(equations, layer_data, 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
+                                        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 atmossphere does not match the simulation domain")
+        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) = LayerData[n, :]
+    (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) = LayerData[n + 1, :]
+    (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 = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql,
-                                                     equations::CompressibleMoistEulerEquations2D)
+    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
@@ -144,8 +144,8 @@ function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerE
     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)
+function perturb_moist_profile!(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
@@ -168,12 +168,12 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql,
     r_l = rho_ql / rho_d
     r_t = r_v + r_l
 
-    # equivalentpotential temperature
+    # equivalent potential 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)
+    aeq_pot = (a * b * c) # TODO: this is not used. remove?
 
     # Assume pressure stays constant
     if (r < rc && Δθ > 0)
@@ -220,11 +220,11 @@ function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql,
 end
 
 # Create background atmosphere data set
-AtmossphereData = AtmossphereLayers(equations)
+atmosphere_data = AtmosphereLayers(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, atmosphere_data)
 end
 
 initial_condition = initial_condition_moist

examples/elixir_moist_euler_source_terms_dry.jl

@@ -1,6 +1,3 @@
-# The same setup as tree_2d_dgsem/elixir_euler_source_terms.jl
-# to verify the StructuredMesh implementation against TreeMesh
-
 using OrdinaryDiffEq
 using Trixi, TrixiAtmo
 using TrixiAtmo: source_terms_convergence_test_dry, initial_condition_convergence_test_dry

examples/elixir_moist_euler_source_terms_moist.jl

@@ -1,6 +1,3 @@
-# The same setup as tree_2d_dgsem/elixir_euler_source_terms.jl
-# to verify the StructuredMesh implementation against TreeMesh
-
 using OrdinaryDiffEq
 using Trixi, TrixiAtmo
 using TrixiAtmo: source_terms_convergence_test_moist,

examples/elixir_moist_euler_source_terms_split_moist.jl

@@ -1,6 +1,3 @@
-# The same setup as tree_2d_dgsem/elixir_euler_source_terms.jl
-# to verify the StructuredMesh implementation against TreeMesh
-
 using OrdinaryDiffEq
 using Trixi, TrixiAtmo
 using TrixiAtmo: source_terms_convergence_test_moist,

src/TrixiAtmo.jl

@@ -11,7 +11,6 @@ module TrixiAtmo
 using Reexport: @reexport
 using Trixi
 using MuladdMacro: @muladd
-@reexport using StaticArrays: SVector, SMatrix
 using Static: True, False
 using StaticArrayInterface: static_size
 using StrideArrays: StrideArray, StaticInt, PtrArray

test/Project.toml

@@ -5,6 +5,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 
 [compat]
+NLsolve = "4.5.1"
 OrdinaryDiffEq = "6.53.2"
 Test = "1"
 Trixi = "0.7, 0.8"

test/runtests.jl

@@ -16,4 +16,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3)
     @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "moist_euler"
         include("test_2d_moist_euler.jl")
     end
+
+    @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection"
+        include("test_spherical_advection.jl")
+    end
 end

test/test_2d_moist_euler.jl

@@ -7,7 +7,7 @@ include("test_trixiatmo.jl") # TODO - This is a repetition from Trixi.jl
 
 EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory for examples?
 
-@trixiatmo_testset "elixir_moist_euler_dry_bubble.jl" begin
+@trixiatmo_testset "elixir_moist_euler_dry_bubble" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_moist_euler_dry_bubble.jl"),
                         l2=[