Merge #310

310: Domain stretching r=kmdeck a=kmdeck

## Purpose 
Allows for variable domain stretching (finer resolution at the top of the domain compared to the bottom).

Whenever we need dz in source code, we obtain it using the function `get_Δz`, which uses ClimaCore to extract the `dz` values. so we do not make the assumption of uniform spacing in src code. (we do in some tests, but that seems OK if those tests use uniform spacing.)
## To-do
Future PRs: fix piracy and add the aqua test back in:
#277 #311 

## Content

1. Adds domain stretching to all domains with a z coordinate.
2. Renames height to depth for the SphericalShell domain
3. uses the stretched domain in the ozark example
4. adds unit tests.

Note that the user passes the target (dz_bottom, dz_top) desired, as a Tuple, or `nothing`, if no stretching is needed. We then create the ClimaCore.Meshes.GeneralizedExponentialStretching object in the domain constructor appropriately. Instead, the user could pass the stretching object itself, and the default could be `ClimaCore.Meshes.Uniform`. In a way this is nicer (it's not prescribed that they use Uniform or GeneralizedExponentialStretching, gives them all the options that ClimaCore develops in the future, no `if` statement in the constructor), but I felt this might be more confusing than supplying the target layer widths at the top and bottom. Happy to change if anyone feels strongly!


Review checklist

I have:
- followed the codebase contribution guide: https://clima.github.io/ClimateMachine.jl/latest/Contributing/
- followed the style guide: https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/
- followed the documentation policy: https://github.com/CliMA/policies/wiki/Documentation-Policy
- checked that this PR does not duplicate an open PR.

In the Content, I have included 
- relevant unit tests, and integration tests, 
- appropriate docstrings on all functions, structs, and modules, and included relevant documentation.

----
- [X] I have read and checked the items on the review checklist.


Co-authored-by: kmdeck <kdeck@caltech.edu>

docs/tutorials/Usage/domain_tutorial.jl

@@ -16,7 +16,7 @@ system of differential equations, where at each coordinate point a value of
 `ρ`, `T`, `P`, and `u⃗` are solved for at each step. The choice of domain is a question "only"
 of geometry: you may be interested in a large eddy simulation (using a box domain), or
 in a global model (where you would need a spherical shell domain
-representing the atmosphere or ocean from some depth to a given height).
+representing the atmosphere or ocean from some depth to z_sfc = 0).
 
 For land surface
 models, each variable is not defined everywhere in space. For example,
@@ -55,7 +55,7 @@ There is a single key method which take a ClimaLSM domain as an argument.
 - [` coordinates(domain)`](https://clima.github.io/ClimaLSM.jl/dev/APIs/shared_utilities/#ClimaLSM.Domains.coordinates): under the hood, this function  uses
 the function space (domain.space) to create the coordinate field. This returns the coordinates as a ClimaCore.Fields.Field object.
 Depending on the domain, the returned coordinate field will have elements of different names and types. For example,
-the SphericalShell domain has coordinates of latitude, longitude, and height, while a Plane domain has coordinates
+the SphericalShell domain has coordinates of latitude, longitude, and depth, while a Plane domain has coordinates
 of x and y, and a Point domain only has a coordinate z_sfc.
 
 

experiments/integrated/ozark/ozark.jl

@@ -425,8 +425,7 @@ Plots.plot!(
 Plots.plot(plt1, plt2, layout = (2, 1))
 Plots.savefig(joinpath(savedir, "stomatal_conductance.png"))
 
-# Current resolution is 3.333 cm per layer, our nodes are at the center of the
-# layers. The second layer is ~ 5cm
+# Current resolution has the first layer at 0.1 cm, the second at 5cm.
 plt1 = Plots.plot(size = (1500, 800))
 Plots.plot!(
     plt1,

experiments/integrated/ozark/ozark_domain.jl

@@ -1,10 +1,16 @@
 # Domain setup
 # For soil column
-nelements = 60
+nelements = 10
 zmin = FT(-2)
 zmax = FT(0)
-land_domain =
-    LSMSingleColumnDomain(; zlim = (zmin, zmax), nelements = nelements)
+dz_bottom = FT(0.5)
+dz_top = FT(0.025)
+
+land_domain = LSMSingleColumnDomain(;
+    zlim = (zmin, zmax),
+    nelements = nelements,
+    dz_tuple = (dz_bottom, dz_top),
+)
 
 # Number of stem and leaf compartments. Leaf compartments are stacked on top of stem compartments
 n_stem = Int64(1)

experiments/integrated/ozark/ozark_simulation.jl

@@ -1,8 +1,9 @@
 t0 = FT(120 * 3600 * 24)# start mid year
 N_days = 120
 tf = t0 + FT(3600 * 24 * N_days)
-dt = FT(60)
-n = 60
+dt = FT(240)
+# Number of timesteps between saving output
+n = 15
 saveat = Array(t0:(n * dt):tf)
 timestepper = CTS.RK4()
 # Set up timestepper

src/shared_utilities/Domains.jl

@@ -63,6 +63,8 @@ struct Column{FT, S} <: AbstractDomain{FT}
     zlim::Tuple{FT, FT}
     "Number of elements used to discretize the interval"
     nelements::Tuple{Int}
+    "Tuple for mesh stretching specifying *target* (dz_bottom, dz_top) (m). If nothing, no stretching is applied."
+    dz_tuple::Union{Tuple{FT, FT}, Nothing}
     "Boundary face identifiers"
     boundary_tags::Tuple{Symbol, Symbol}
     "The associated ClimaCore Space"
@@ -72,24 +74,55 @@ end
 """
     Column(;
         zlim::Tuple{FT, FT},
-        nelements::Int) where {FT}
-Outer constructor for the `Column` type.
+        nelements::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing) where {FT}
+
+Outer constructor for the `Column` type. 
+
+Using `ClimaCore` tools, the coordinate
+mesh can be stretched such that the top of the domain has finer resolution
+than the bottom of the domain. In order to activate this, a tuple with the
+target dz_bottom and dz_top should be passed via keyword argument. The default is
+uniform spacing. Please note that in order to use this feature, ClimaCore requires that
+the elements of zlim be <=0. Additionally, the dz_tuple you supply may not be compatible
+with the domain boundaries in some cases, in which case you may need to choose
+different values.
+
 The `boundary_tags` field values are used to label the boundary faces 
 at the top and bottom of the domain.
 """
-function Column(; zlim::Tuple{FT, FT}, nelements::Int) where {FT}
+function Column(;
+    zlim::Tuple{FT, FT},
+    nelements::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
+) where {FT}
     @assert zlim[1] < zlim[2]
     boundary_tags = (:bottom, :top)
     column = ClimaCore.Domains.IntervalDomain(
         ClimaCore.Geometry.ZPoint{FT}(zlim[1]),
         ClimaCore.Geometry.ZPoint{FT}(zlim[2]);
         boundary_tags = boundary_tags,
     )
-    mesh = ClimaCore.Meshes.IntervalMesh(column; nelems = nelements)
+    if dz_tuple isa Nothing
+        mesh = ClimaCore.Meshes.IntervalMesh(column; nelems = nelements)
+    else
+        @assert zlim[2] <= 0
+        mesh = ClimaCore.Meshes.IntervalMesh(
+            column,
+            ClimaCore.Meshes.GeneralizedExponentialStretching{FT}(
+                dz_tuple[1],
+                dz_tuple[2],
+            );
+            nelems = nelements,
+            reverse_mode = true,
+        )
+    end
+
     center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(mesh)
     return Column{FT, typeof(center_space)}(
         zlim,
         (nelements,),
+        dz_tuple,
         boundary_tags,
         center_space,
     )
@@ -174,6 +207,7 @@ end
         xlim::Tuple{FT, FT}
         ylim::Tuple{FT, FT}
         zlim::Tuple{FT, FT}
+        dz_tuple::Union{Tuple{FT, FT}, Nothing}
         nelements::Tuple{Int, Int, Int}
         npolynomial::Int
         periodic::Tuple{Bool, Bool}
@@ -194,6 +228,8 @@ struct HybridBox{FT, S} <: AbstractDomain{FT}
     ylim::Tuple{FT, FT}
     "Domain interval limits along z axis, in meters"
     zlim::Tuple{FT, FT}
+    "Tuple for mesh stretching specifying *target* (dz_bottom, dz_top) (m). If nothing, no stretching is applied."
+    dz_tuple::Union{Tuple{FT, FT}, Nothing}
     "Number of elements to discretize interval, (nx, ny,nz)"
     nelements::Tuple{Int, Int, Int}
     " Polynomial order for the horizontal directions"
@@ -211,6 +247,7 @@ end
         zlim::Tuple{FT, FT},
         nelements::Tuple{Int, Int, Int},
         npolynomial::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
         periodic = (true, true),
     ) where {FT}
 Constructs the `HybridBox` domain
@@ -222,13 +259,23 @@ to the x-axis, the second corresponding to the y-axis, and the third
 corresponding to the z-axis. The domain is periodic at the (xy) boundaries,
 and the function space is of polynomial order `npolynomial` in the
 horizontal directions.
+
+Using `ClimaCore` tools, the coordinate
+mesh can be stretched such that the top of the domain has finer resolution
+than the bottom of the domain. In order to activate this, a tuple with the
+target dz_bottom and dz_top should be passed via keyword argument. The default is
+uniform spacing. Please note that in order to use this feature, ClimaCore requires that
+the elements of zlim be <=0. Additionally, the dz_tuple you supply may not be compatible
+with the domain boundaries in some cases, in which case you may need to choose
+different values.
 """
 function HybridBox(;
     xlim::Tuple{FT, FT},
     ylim::Tuple{FT, FT},
     zlim::Tuple{FT, FT},
     nelements::Tuple{Int, Int, Int},
     npolynomial::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     periodic = (true, true),
 ) where {FT}
     @assert xlim[1] < xlim[2]
@@ -240,7 +287,21 @@ function HybridBox(;
         ClimaCore.Geometry.ZPoint(zlim[2]);
         boundary_tags = (:bottom, :top),
     )
-    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, nelems = nelements[3])
+    if dz_tuple isa Nothing
+        vertmesh =
+            ClimaCore.Meshes.IntervalMesh(vertdomain; nelems = nelements[3])
+    else
+        @assert zlim[2] <= 0
+        vertmesh = ClimaCore.Meshes.IntervalMesh(
+            vertdomain,
+            ClimaCore.Meshes.GeneralizedExponentialStretching{FT}(
+                dz_tuple[1],
+                dz_tuple[2],
+            );
+            nelems = nelements[3],
+            reverse_mode = true,
+        )
+    end
     vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)
 
     horzdomain = Plane(;
@@ -260,6 +321,7 @@ function HybridBox(;
         xlim,
         ylim,
         zlim,
+        dz_tuple,
         nelements,
         npolynomial,
         periodic,
@@ -270,7 +332,8 @@ end
 """
     struct SphericalShell{FT} <: AbstractDomain{FT}
         radius::FT
-        height::FT
+        depth::FT
+        dz_tuple::Union{Tuple{FT, FT}, Nothing}
         nelements::Tuple{Int, Int}
         npolynomial::Int
     end
@@ -285,7 +348,9 @@ struct SphericalShell{FT, S} <: AbstractDomain{FT}
     "The radius of the shell"
     radius::FT
     "The radial extent of the shell"
-    height::FT
+    depth::FT
+    "Tuple for mesh stretching specifying *target* (dz_bottom, dz_top) (m). If nothing, no stretching is applied."
+    dz_tuple::Union{Tuple{FT, FT}, Nothing}
     "The number of elements to be used in the non-radial and radial directions"
     nelements::Tuple{Int, Int}
     "The polynomial order to be used in the non-radial directions"
@@ -297,31 +362,49 @@ end
 """
     SphericalShell(;
         radius::FT,
-        height::FT,
+        depth::FT,
         nelements::Tuple{Int, Int},
         npolynomial::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     ) where {FT}
 Outer constructor for the `SphericalShell` domain, using keyword arguments.
+
+Using `ClimaCore` tools, the coordinate
+mesh can be stretched such that the top of the domain has finer resolution
+than the bottom of the domain. In order to activate this, a tuple with the
+target dz_bottom and dz_top should be passed via keyword argument. The default is
+uniform spacing. Please note that the dz_tuple you supply may not be compatible
+with the depth/nelements chosen, in which case you may need to choose
+different values.
 """
 function SphericalShell(;
     radius::FT,
-    height::FT,
+    depth::FT,
     nelements::Tuple{Int, Int},
     npolynomial::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
 ) where {FT}
     @assert 0 < radius
-    @assert 0 < height
+    @assert 0 < depth
     vertdomain = ClimaCore.Domains.IntervalDomain(
-        ClimaCore.Geometry.ZPoint(FT(0)),
-        ClimaCore.Geometry.ZPoint(FT(height));
+        ClimaCore.Geometry.ZPoint(FT(-depth)),
+        ClimaCore.Geometry.ZPoint(FT(0));
         boundary_tags = (:bottom, :top),
     )
-
-    vertmesh = ClimaCore.Meshes.IntervalMesh(
-        vertdomain,
-        ClimaCore.Meshes.Uniform(),
-        nelems = nelements[2],
-    )
+    if dz_tuple isa Nothing
+        vertmesh =
+            ClimaCore.Meshes.IntervalMesh(vertdomain; nelems = nelements[2])
+    else
+        vertmesh = ClimaCore.Meshes.IntervalMesh(
+            vertdomain,
+            ClimaCore.Meshes.GeneralizedExponentialStretching{FT}(
+                dz_tuple[1],
+                dz_tuple[2],
+            );
+            nelems = nelements[2],
+            reverse_mode = true,
+        )
+    end
     vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)
 
     horzdomain = ClimaCore.Domains.SphereDomain(radius)
@@ -336,7 +419,8 @@ function SphericalShell(;
     )
     return SphericalShell{FT, typeof(hv_center_space)}(
         radius,
-        height,
+        depth,
+        dz_tuple,
         nelements,
         npolynomial,
         hv_center_space,
@@ -424,12 +508,18 @@ end
     LSMSingleColumnDomain(;
         zlim::Tuple{FT, FT},
         nelements::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     ) where {FT}
 A constructor for the LSMSingleColumnDomain.
 """
-function LSMSingleColumnDomain(; zlim::Tuple{FT, FT}, nelements::Int) where {FT}
+function LSMSingleColumnDomain(;
+    zlim::Tuple{FT, FT},
+    nelements::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
+) where {FT}
     @assert zlim[1] < zlim[2]
-    subsurface_domain = Column(; zlim = FT.(zlim), nelements = nelements)
+    subsurface_domain =
+        Column(; zlim = FT.(zlim), nelements = nelements, dz_tuple = dz_tuple)
     surface_space = obtain_surface_space(subsurface_domain.space)
     surface_domain = Point{FT, typeof(surface_space)}(zlim[2], surface_space)
     return LSMSingleColumnDomain{
@@ -469,7 +559,8 @@ end
         zlim::Tuple{FT, FT}
         nelements::Tuple{Int, Int, Int}
         npolynomial::Int
-        periodic::Tuple{Bool, Bool}
+        periodic::Tuple{Bool, Bool} = (true, true),
+        dz_tuple::Union{Tuple{FT,FT}, Nothing} = nothing,
     ) where {FT}
 A constructor for the LSMMultiColumnDomain.
 """
@@ -479,7 +570,8 @@ function LSMMultiColumnDomain(;
     zlim::Tuple{FT, FT},
     nelements::Tuple{Int, Int, Int},
     npolynomial::Int,
-    periodic::Tuple{Bool, Bool},
+    periodic::Tuple{Bool, Bool} = (true, true),
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
 ) where {FT}
     @assert xlim[1] < xlim[2]
     @assert ylim[1] < ylim[2]
@@ -489,6 +581,7 @@ function LSMMultiColumnDomain(;
         xlim = xlim,
         ylim = ylim,
         zlim = zlim,
+        dz_tuple = dz_tuple,
         nelements = nelements,
         npolynomial = npolynomial,
         periodic = periodic,
@@ -536,23 +629,26 @@ end
 """
     LSMSphericalShellDomain(;
         radius::FT,
-        height::FT,
+        depth::FT,
         nelements::Tuple{Int, Int},
         npolynomial::Int,
+        dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
     ) where {FT}
 A constructor for the LSMSphericalShellDomain.
 """
 function LSMSphericalShellDomain(;
     radius::FT,
-    height::FT,
+    depth::FT,
     nelements::Tuple{Int, Int},
     npolynomial::Int,
+    dz_tuple::Union{Tuple{FT, FT}, Nothing} = nothing,
 ) where {FT}
     @assert radius > FT(0)
-    @assert height > FT(0)
+    @assert depth > FT(0)
     subsurface = SphericalShell(;
         radius = radius,
-        height = height,
+        depth = depth,
+        dz_tuple = dz_tuple,
         nelements = nelements,
         npolynomial = npolynomial,
     )