Add regional bucket

.buildkite/pipeline.yml

@@ -154,6 +154,13 @@ steps:
         command: "julia --color=yes --project=.buildkite experiments/standalone/Bucket/global_bucket_staticmap.jl"
         artifact_paths: "experiments/standalone/Bucket/artifacts_staticmap/*cpu*"
 
+      - label: "Regional Bucket on CPU (static map albedo)"
+        key: "regional_bucket_staticmap_cpu"
+        command: "julia --color=yes --project=.buildkite experiments/standalone/Bucket/global_bucket_staticmap.jl"
+        artifact_paths: "experiments/standalone/Bucket/artifacts_staticmap_regional/*cpu*"
+        env:
+          CLIMALAND_CI_REGIONAL_BUCKET: true
+
       - label: "Global Bucket on CPU (temporal map albedo)"
         key: "global_bucket_temporalmap_cpu"
         command: "julia --color=yes --project=.buildkite experiments/standalone/Bucket/global_bucket_temporalmap.jl"

experiments/standalone/Bucket/global_bucket_staticmap.jl

@@ -1,7 +1,7 @@
 # # Global bucket run using spatial map albedo
 
 # The code sets up and runs the bucket for 7 days using albedo read in from
-# a file containing static data over the globe, and analytic atmospheric and
+# a file containing static data over the globe or on a region, and analytic atmospheric and
 # radiative forcings.
 # Moving forward, this driver will serve as our more complex global bucket run,
 # eventually running for a longer time period (1+ year) and using temporally
@@ -63,26 +63,45 @@ function compute_clims(v)
     return (minimum(means) - maximum(sigmas), maximum(means) + maximum(sigmas))
 end
 
-anim_plots = true
+anim_plots = false
+# Set to true if you want to run a regional simulation. By default, it is false,
+# unless the `CLIMALAND_CI_REGIONAL_BUCKET` environment variable is defined.
+regional_simulation = haskey(ENV, "CLIMALAND_CI_REGIONAL_BUCKET")
+regional_str = regional_simulation ? "_regional" : ""
 regridder_type = :InterpolationsRegridder
 FT = Float64;
 context = ClimaComms.context()
 earth_param_set = LP.LandParameters(FT);
 outdir = joinpath(
     pkgdir(ClimaLand),
-    "experiments/standalone/Bucket/artifacts_staticmap",
+    "experiments/standalone/Bucket/artifacts_staticmap$(regional_str)",
 )
 !ispath(outdir) && mkpath(outdir)
 
 # Set up simulation domain
 soil_depth = FT(3.5);
-bucket_domain = ClimaLand.Domains.SphericalShell(;
-    radius = FT(6.3781e6),
-    depth = soil_depth,
-    nelements = (30, 10),
-    npolynomial = 1,
-    dz_tuple = FT.((1.0, 0.05)),
-);
+if regional_simulation
+    @info "Running regional simulation"
+    center_long, center_lat = FT(-118.14452), FT(34.14778)
+    # Half size of the grid in meters
+    delta_m = FT(50_000)
+    bucket_domain = ClimaLand.Domains.HybridBox(;
+        xlim = (delta_m, delta_m),
+        ylim = (delta_m, delta_m),
+        zlim = (-soil_depth, FT(0.0)),
+        longlat = (center_long, center_lat),
+        nelements = (30, 30, 10),
+        npolynomial = 1,
+    )
+else
+    bucket_domain = ClimaLand.Domains.SphericalShell(;
+        radius = FT(6.3781e6),
+        depth = soil_depth,
+        nelements = (30, 10),
+        npolynomial = 1,
+        dz_tuple = FT.((1.0, 0.05)),
+    )
+end
 ref_time = DateTime(2021);
 
 # Set up parameters
@@ -258,8 +277,31 @@ sol = ClimaComms.@time ClimaComms.device() SciMLBase.solve(
 
 # Interpolate to grid
 space = axes(coords.surface)
-longpts = range(-180.0, 180.0, 21)
-latpts = range(-90.0, 90.0, 21)
+num_pts = 21
+if regional_simulation
+    radius_earth = FT(6.378e6)
+    xlim = (
+        center_long - delta_m / (2radius_earth),
+        center_long + delta_m / (2radius_earth),
+    )
+    ylim = (
+        center_lat - delta_m / (2radius_earth),
+        center_lat + delta_m / (2radius_earth),
+    )
+    longpts = range(xlim[1], xlim[2], num_pts)
+    latpts = range(ylim[1], ylim[2], num_pts)
+
+    # Temporary monkey patch until we use the new diagnostics.
+
+    # This is used somewhere internally by the remapper.
+    ClimaCore.Geometry._coordinate(pt::ClimaCore.Geometry.LatLongPoint, ::Val{1}) = ClimaCore.Geometry.LongPoint(pt.long)
+    ClimaCore.Geometry._coordinate(pt::ClimaCore.Geometry.LatLongPoint, ::Val{2}) = ClimaCore.Geometry.LatPoint(pt.lat)
+
+else
+    longpts = range(-180.0, 180.0, num_pts)
+    latpts = range(-90.0, 90.0, num_pts)
+    hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
+end
 hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
 remapper = Remapping.Remapper(space, hcoords)
 

test/shared_utilities/domains.jl

@@ -128,6 +128,7 @@ FT = Float32
     @test obtain_surface_domain(box) == Plane{FT}(
         xlim,
         ylim,
+        nothing,
         nelements[1:2],
         (true, true),
         0,
@@ -169,24 +170,82 @@ FT = Float32
           ClimaCore.Spaces.SpectralElementSpace2D
 
     # Plane latlong
+    dxlim = (FT(-50_000), FT(80_000))
+    dylim = (FT(-30_000), FT(40_000))
+    longlat = (FT(-118.14452), FT(34.14778))
+    radius_earth = FT(6.378e6)
+    xlim_longlat = (
+        longlat[1] - dxlim[1] / 2radius_earth,
+        longlat[1] + dxlim[2] / 2radius_earth,
+    )
+    ylim_longlat = (
+        longlat[2] - dylim[1] / 2radius_earth,
+        longlat[2] + dylim[2] / (2radius_earth),
+    )
+
     longlat_plane = Plane(;
-        xlim = xlim,
-        ylim = ylim,
-        longlat = (FT(-118.14452), FT(34.14778)),
+        xlim = dxlim,
+        ylim = dylim,
+        longlat,
         nelements = nelements[1:2],
         npolynomial = 0,
     )
     plane_coords = coordinates(longlat_plane).surface
-    @test eltype(plane_coords) == ClimaCore.Geometry.LongLatPoint{FT}
+    @test eltype(plane_coords) == ClimaCore.Geometry.LatLongPoint{FT}
     @test typeof(plane_coords) <: ClimaCore.Fields.Field
-    @test xy_plane.xlim == FT.(xlim)
-    @test xy_plane.ylim == FT.(ylim)
-    @test xy_plane.nelements == nelements[1:2]
-    @test xy_plane.npolynomial == 0
-    @test xy_plane.periodic == (false, false)
-    @test typeof(xy_plane.space.surface) <:
+    @test longlat_plane.xlim == FT.(xlim_longlat)
+    @test longlat_plane.ylim == FT.(ylim_longlat)
+    @test longlat_plane.nelements == nelements[1:2]
+    @test longlat_plane.npolynomial == 0
+    @test longlat_plane.periodic == (false, false)
+    @test typeof(longlat_plane.space.surface) <:
           ClimaCore.Spaces.SpectralElementSpace2D
 
+    # Box latlong
+    longlat_box = HybridBox(;
+        xlim = dxlim,
+        ylim = dylim,
+        zlim = zlim,
+        longlat,
+        nelements = nelements,
+        npolynomial = 0,
+    )
+    @test longlat_box.fields.z ==
+          ClimaCore.Fields.coordinate_field(longlat_box.space.subsurface).z
+    face_space = obtain_face_space(longlat_box.space.subsurface)
+    z_face = ClimaCore.Fields.coordinate_field(face_space).z
+    @test longlat_box.fields.z_sfc ==
+          top_face_to_surface(z_face, longlat_box.space.surface)
+    Δz_top, Δz_bottom = get_Δz(longlat_box.fields.z)
+    @test longlat_box.fields.Δz_top == Δz_top
+    @test longlat_box.fields.Δz_bottom == Δz_bottom
+    longlat_box_coords = coordinates(longlat_box).subsurface
+    @test eltype(longlat_box_coords) == ClimaCore.Geometry.LatLongZPoint{FT}
+    @test typeof(longlat_box_coords) <: ClimaCore.Fields.Field
+    @test longlat_box.xlim == FT.(xlim_longlat)
+    @test longlat_box.ylim == FT.(ylim_longlat)
+    @test longlat_box.zlim == FT.(zlim)
+    @test longlat_box.nelements == nelements
+    @test longlat_box.npolynomial == 0
+    @test longlat_box.periodic == (false, false)
+    @test typeof(
+        ClimaCore.Spaces.horizontal_space(longlat_box.space.subsurface),
+    ) <: ClimaCore.Spaces.SpectralElementSpace2D
+    @test typeof(longlat_box.space.subsurface) <:
+          ClimaCore.Spaces.CenterExtrudedFiniteDifferenceSpace
+    @test typeof(longlat_box.space.surface) <:
+          ClimaCore.Spaces.SpectralElementSpace2D
+    @test obtain_surface_space(longlat_box.space.subsurface) ==
+          longlat_box.space.surface
+    @test obtain_surface_domain(longlat_box) == Plane{FT}(
+        xlim_longlat,
+        ylim_longlat,
+        longlat,
+        nelements[1:2],
+        (false, false),
+        0,
+        (; surface = longlat_box.space.surface),
+    )
 
     # Column