Try using 2 elements

config/model_configs/box_analytic_cosine_mountain_test.yml

@@ -11,7 +11,6 @@ z_elem: 100
 z_stretch: false
 dt: "6secs"
 t_end: "24hours"
-use_reference_state: false
 rayleigh_sponge: true
 toml: [toml/analytic_mountain_test.toml]
 analytic_check: true

config/model_configs/plane_analytic_cosine_mountain_test.yml

@@ -1,17 +1,14 @@
-config: "box"
+config: "plane"
 FLOAT_TYPE: "Float64"
 initial_condition: "ConstantBuoyancyFrequencyProfile"
 topography: "Cosine2D"
 x_max: 300e3
-y_max: 300e3
 z_max: 21e3
 x_elem: 30
-y_elem: 1
 z_elem: 100
 z_stretch: false
 dt: "6secs"
 t_end: "24hours"
-use_reference_state: false
 rayleigh_sponge: true
 toml: [toml/analytic_mountain_test.toml]
 analytic_check: true

config/model_configs/plane_analytic_schar_mountain_test.yml

@@ -3,15 +3,14 @@ FLOAT_TYPE: "Float64"
 initial_condition: "ConstantBuoyancyFrequencyProfile"
 topography: "Schar"
 x_max: 300e3
-y_max: 300e3
+y_max: 2e3
 z_max: 21e3
 x_elem: 300
-y_elem: 1
+y_elem: 2
 z_elem: 100
 z_stretch: false
 dt: "0.5secs"
 t_end: "24hours"
-use_reference_state: false
 rayleigh_sponge: true
 toml: [toml/analytic_mountain_test.toml]
 analytic_check: true

post_processing/ci_plots.jl

@@ -699,17 +699,18 @@ function column_rms(var)
 end
 
 const AnalyticMountainTest = Union{
-    Val{:plane_analytic_schar_mountain_test},
     Val{:plane_analytic_cosine_mountain_test},
-    Val{:box_analytic_cosine_mountain_test},
+    Val{:gpu_plane_analytic_schar_mountain_test},
+    Val{:gpu_plane_analytic_cosine_mountain_test},
+    Val{:gpu_box_analytic_cosine_mountain_test},
 }
 
 function make_plots(
     val::AnalyticMountainTest,
     output_paths::Vector{<:AbstractString},
 )
     simdirs = SimDir.(output_paths)
-    is_3d = val == Val(:box_analytic_cosine_mountain_test)
+    is_3d = val != Val(:plane_analytic_cosine_mountain_test)
 
     short_names =
         ["ua", "wa", "uapredicted", "wapredicted", "uaerror", "waerror"]
@@ -731,7 +732,7 @@ function make_plots(
         output_name = "slice_summary",
     )
 
-    if val == Val(:plane_analytic_schar_mountain_test)
+    if val == Val(:gpu_plane_analytic_schar_mountain_test)
         zoomed_slice_vars =
             Iterators.flatmap(["wa", "wapredicted", "waerror"]) do short_name
                 map(simdirs) do simdir
@@ -767,60 +768,60 @@ function make_plots(
         output_name = "rms_errors_24h",
     )
 
-    collected_vars = collect(slice_summary_vars)
-    (; dims, dim_attributes, dim2index, index2dim) = collected_vars[1]
-    u_data = collected_vars[4].data
-    w_data = collected_vars[8].data
-    predicted_u_data = collected_vars[9].data
-    predicted_w_data = collected_vars[10].data
-
-    FT = eltype(predicted_u_data)
-    background_velocity = SVector(CA.background_u(FT), FT(0))
-    perturbation_data = @. SVector(u_data, w_data) - (background_velocity,)
-    predicted_perturbation_data =
-        @. SVector(predicted_u_data, predicted_w_data) - (background_velocity,)
-
-    # Rescale the predicted perturbations so that their range no less than the
-    # background air speed. Apply the same scaling to the actual perturbations.
-    background_speed = norm(background_velocity)
-    predicted_perturbation_range =
-        -(reverse(extrema(norm, predicted_perturbation_data))...)
-    s = max(1, round(Int, background_speed / predicted_perturbation_range))
-    uw_data = @. s * perturbation_data + (background_velocity,)
-    predicted_uw_data =
-        @. s * predicted_perturbation_data + (background_velocity,)
-
-    uw_var = ClimaAnalysis.OutputVar(
-        Dict{String, String}(
-            "units" => "m s^-1",
-            "short_name" => "ū + $s u′",
-            "long_name" => "Air Speed with $(s)x Perturbations, Instantaneous",
-        ),
-        dims,
-        dim_attributes,
-        uw_data,
-        dim2index,
-        index2dim,
-    )
-    predicted_uw_var = ClimaAnalysis.OutputVar(
-        Dict{String, String}(
-            "units" => "m s^-1",
-            "short_name" => "predicted ū + $s u′",
-            "long_name" => "Predicted Air Speed with $(s)x Perturbations, Instantaneous",
-        ),
-        dims,
-        dim_attributes,
-        predicted_uw_data,
-        dim2index,
-        index2dim,
-    )
-
-    make_plots_generic(
-        output_paths,
-        [uw_var, predicted_uw_var];
-        plot_fn = plot_streamlines!,
-        output_name = "streamlines_24h",
-    )
+    # collected_vars = collect(slice_summary_vars)
+    # (; dims, dim_attributes, dim2index, index2dim) = collected_vars[1]
+    # u_data = collected_vars[4].data
+    # w_data = collected_vars[8].data
+    # predicted_u_data = collected_vars[9].data
+    # predicted_w_data = collected_vars[10].data
+
+    # FT = eltype(predicted_u_data)
+    # background_velocity = SVector(CA.background_u(FT), FT(0))
+    # perturbation_data = @. SVector(u_data, w_data) - (background_velocity,)
+    # predicted_perturbation_data =
+    #     @. SVector(predicted_u_data, predicted_w_data) - (background_velocity,)
+
+    # # Rescale the predicted perturbations so that their range no less than the
+    # # background air speed. Apply the same scaling to the actual perturbations.
+    # background_speed = norm(background_velocity)
+    # predicted_perturbation_range =
+    #     -(reverse(extrema(norm, predicted_perturbation_data))...)
+    # s = max(1, round(Int, background_speed / predicted_perturbation_range))
+    # uw_data = @. s * perturbation_data + (background_velocity,)
+    # predicted_uw_data =
+    #     @. s * predicted_perturbation_data + (background_velocity,)
+
+    # uw_var = ClimaAnalysis.OutputVar(
+    #     Dict{String, String}(
+    #         "units" => "m s^-1",
+    #         "short_name" => "ū + $s u′",
+    #         "long_name" => "Air Speed with $(s)x Perturbations, Instantaneous",
+    #     ),
+    #     dims,
+    #     dim_attributes,
+    #     uw_data,
+    #     dim2index,
+    #     index2dim,
+    # )
+    # predicted_uw_var = ClimaAnalysis.OutputVar(
+    #     Dict{String, String}(
+    #         "units" => "m s^-1",
+    #         "short_name" => "predicted ū + $s u′",
+    #         "long_name" => "Predicted Air Speed with $(s)x Perturbations, Instantaneous",
+    #     ),
+    #     dims,
+    #     dim_attributes,
+    #     predicted_uw_data,
+    #     dim2index,
+    #     index2dim,
+    # )
+
+    # make_plots_generic(
+    #     output_paths,
+    #     [uw_var, predicted_uw_var];
+    #     plot_fn = plot_streamlines!,
+    #     output_name = "streamlines_24h",
+    # )
 end
 
 MountainPlots = Union{

src/analytic_solutions/analytic_solutions.jl

@@ -11,8 +11,10 @@ import StaticArrays: @SVector
 # of quadrature. TODO: QuadGK does not seem to work on GPUs.
 # definite_integral(f::F, x1, x2, maxevals, order) where {F} =
 #     QuadGK.quadgk(f, x1, x2; maxevals, order, rtol = 0)[1]
-definite_integral(f::F, x1, x2, maxevals, _) where {F} =
-    sum(f, LinRange(x1, x2, maxevals)) * (x2 - x1) / (maxevals - 1)
+function definite_integral(f::F, x1, x2, maxevals, _) where {F}
+    dx = (x2 - x1) / maxevals
+    return sum(index -> f(x1 + index * dx), 2:maxevals; init = f(x1 + dx)) * dx
+end
 
 background_N(::Type{FT}) where {FT} = FT(0.01)
 background_T_sfc(::Type{FT}) where {FT} = FT(288)
@@ -166,7 +168,7 @@ function schar_predicted_velocity(params, coord)
     n_efolds = -log(eps(FT))
     k_x_max = k_peak + 2 * sqrt(n_efolds) / a
     (Δu, Δv, Δw) =
-        2 * definite_integral(FT(0), k_x_max, 10, 20) do k_x
+        2 * definite_integral(FT(0), k_x_max, 100, 20) do k_x
             _FT = typeof(k_x) # Redefine FT to prevent it from getting boxed.
             Fh =
                 h_max * a / (8 * sqrt(_FT(π))) * (