Add ConformalCubedSphere grid via MultiRegion module

.buildkite/pipeline.yml

@@ -494,7 +494,7 @@ steps:
       architecture: CPU
     depends_on: "init_cpu"
 
-  - label: "🦾 cpu distributed sovers tests"
+  - label: "🦾 cpu distributed solvers tests"
     env:
       JULIA_DEPOT_PATH: "$TARTARUS_HOME/.julia-$BUILDKITE_BUILD_NUMBER"
       TEST_GROUP: "distributed_solvers"

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.2"
 manifest_format = "2.0"
-project_hash = "9044a5920a57a78e707ac7490810d9c46b9bc546"
+project_hash = "ea620612cb5f84093b962d4345aa7d1b1271739c"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -157,6 +157,12 @@ version = "1.0.0"
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[deps.Distances]]
+deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
+git-tree-sha1 = "49eba9ad9f7ead780bfb7ee319f962c811c6d3b2"
+uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+version = "0.10.8"
+
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -676,6 +682,12 @@ deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 version = "1.9.0"
 
+[[deps.StatsAPI]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7"
+uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+version = "1.6.0"
+
 [[deps.Strided]]
 deps = ["LinearAlgebra", "StridedViews", "TupleTools"]
 git-tree-sha1 = "b32eadf6ac726a790567fdc872b63117712e16a8"

Project.toml

@@ -9,6 +9,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
 CubedSphere = "7445602f-e544-4518-8976-18f8e8ae6cdb"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
@@ -39,6 +40,7 @@ Adapt = "3"
 CUDA = "4"
 Crayons = "4"
 CubedSphere = "0.1, 0.2"
+Distances = "0.10"
 DocStringExtensions = "0.8, 0.9"
 FFTW = "1"
 Glob = "1.3"

docs/Manifest.toml

@@ -271,9 +271,9 @@ version = "0.9.3"
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "58fea7c536acd71f3eef6be3b21c0df5f3df88fd"
+git-tree-sha1 = "39fd748a73dce4c05a9655475e437170d8fb1b67"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.27.24"
+version = "0.27.25"
 
 [[deps.DocumenterCitations]]
 deps = ["Bibliography", "DataStructures", "Documenter", "Markdown", "Unicode"]
@@ -355,9 +355,9 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.FillArrays]]
 deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
-git-tree-sha1 = "0b3b52afd0f87b0a3f5ada0466352d125c9db458"
+git-tree-sha1 = "2250347838b28a108d1967663cba57bfb3c02a58"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.2.1"
+version = "1.3.0"
 
 [[deps.FixedPointNumbers]]
 deps = ["Statistics"]
@@ -876,9 +876,9 @@ version = "0.5.4"
 
 [[deps.OffsetArrays]]
 deps = ["Adapt"]
-git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4"
+git-tree-sha1 = "2ac17d29c523ce1cd38e27785a7d23024853a4bb"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.12.9"
+version = "1.12.10"
 
 [[deps.Ogg_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -1232,10 +1232,14 @@ uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
 version = "0.1.1"
 
 [[deps.StaticArrays]]
-deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
-git-tree-sha1 = "832afbae2a45b4ae7e831f86965469a24d1d8a83"
+deps = ["LinearAlgebra", "Random", "StaticArraysCore"]
+git-tree-sha1 = "0da7e6b70d1bb40b1ace3b576da9ea2992f76318"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.5.26"
+version = "1.6.0"
+weakdeps = ["Statistics"]
+
+    [deps.StaticArrays.extensions]
+    StaticArraysStatisticsExt = "Statistics"
 
 [[deps.StaticArraysCore]]
 git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
@@ -1397,16 +1401,16 @@ uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
 version = "1.1.34+0"
 
 [[deps.Xorg_libX11_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
-git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
+git-tree-sha1 = "afead5aba5aa507ad5a3bf01f58f82c8d1403495"
 uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
-version = "1.6.9+4"
+version = "1.8.6+0"
 
 [[deps.Xorg_libXau_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "6035850dcc70518ca32f012e46015b9beeda49d8"
 uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
-version = "1.0.9+4"
+version = "1.0.11+0"
 
 [[deps.Xorg_libXdmcp_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1433,10 +1437,10 @@ uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
 version = "0.1.1+0"
 
 [[deps.Xorg_libxcb_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
-git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
+git-tree-sha1 = "b4bfde5d5b652e22b9c790ad00af08b6d042b97d"
 uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
-version = "1.13.0+3"
+version = "1.15.0+0"
 
 [[deps.Xorg_xtrans_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]

src/AbstractOperations/at.jl

@@ -58,52 +58,72 @@ indices(f::Function) = default_indices(3)
 indices(f::Number)   = default_indices(3)
 
 """
-    interpolate_indices(operands..., loc_operation = abstract_operation_location)
+    intersect_indices(loc, operands...)
 
-Utility to propagate operands' indices in `AbstractOperations`s with multiple operands 
-(`BinaryOperation`s and `MultiaryOperation`s).
+Utility to compute the intersection of `operands' indices.
 """
-function interpolate_indices(args...; loc_operation = (Center, Center, Center))
-    idxs = Any[:, :, :]
-    for i in 1:3
-        for arg in args
-            idxs[i] = interpolate_index(indices(arg)[i], idxs[i], location(arg)[i], loc_operation[i])
-        end
-    end
+function intersect_indices(loc, operands...)
 
-    return Tuple(idxs)
+    idx1 = compute_index_intersection(Colon(), loc[1], operands...; dim=1)
+    idx2 = compute_index_intersection(Colon(), loc[2], operands...; dim=2)
+    idx3 = compute_index_intersection(Colon(), loc[3], operands...; dim=3)
+
+    return (idx1, idx2, idx3)
 end
 
-interpolate_index(::Colon, ::Colon, args...)       = Colon()
-interpolate_index(::Colon, b::UnitRange, args...)  = b
-
-instantiate(T::Type) = T()
-instantiate(t) = t
+# Fallback for `KernelFunctionOperation`s with no argument 
+compute_index_intersection(::Colon, to_loc; kw...) = Colon()
 
-function interpolate_index(a::UnitRange, ::Colon, loc, new_loc)  
-    a = corrected_index(a, instantiate(loc), instantiate(new_loc))
+compute_index_intersection(to_idx, to_loc, op; dim) =
+    _compute_index_intersection(to_idx, indices(op)[dim],
+                                to_loc, location(op, dim))
 
-    # Abstract operations that require an interpolation of a sliced fields are not supported!
-    first(a) > last(a) && throw(ArgumentError("Cannot interpolate a sliced field from $loc to $(new_loc)!"))
-    return a
+"""Compute index intersection recursively for `dim`ension ∈ (1, 2, 3)."""
+function compute_index_intersection(to_idx, to_loc, op1, op2, more_ops...; dim)
+    new_to_idx = _compute_index_intersection(to_idx, indices(op1)[dim], to_loc, location(op1, dim))
+    return compute_index_intersection(new_to_idx, to_loc, op2, more_ops...; dim)
 end
 
-function interpolate_index(a::UnitRange, b::UnitRange, loc, new_loc)   
-    a = corrected_index(a, instantiate(loc), instantiate(new_loc))
+# Life is pretty simple in this case.
+_compute_index_intersection(to_idx::Colon, from_idx::Colon, args...) = Colon()
+
+# Because `from_idx` imposes no restrictions, we just return `to_idx`.
+_compute_index_intersection(to_idx::UnitRange, from_idx::Colon, args...) = to_idx
 
-    # Abstract operations that require an interpolation of a sliced fields are not supported!
-    first(a) > last(a) && throw(ArgumentError("Cannot interpolate a sliced field from $loc to $(new_loc)!"))
+# This time we account for the possible range-reducing effect of interpolation on `from_idx`.
+function _compute_index_intersection(to_idx::Colon, from_idx::UnitRange, to_loc, from_loc)
+    shifted_idx = restrict_index_for_interpolation(from_idx, from_loc, to_loc)
+    validate_shifted_index(shifted_idx)
+    return shifted_idx
+end
+
+# Compute the intersection of two index ranges
+function _compute_index_intersection(to_idx::UnitRange, from_idx::UnitRange, to_loc, from_loc)
+    shifted_idx = restrict_index_for_interpolation(from_idx, from_loc, to_loc)
+    validate_shifted_index(shifted_idx)
 
-    indices = UnitRange(max(first(a), first(b)), min(last(a), last(b)))
+    range_intersection = UnitRange(max(first(shifted_idx), first(to_idx)), min(last(shifted_idx), last(to_idx)))
 
-    # Abstract operations between parallel non-intersecating windowed fields are not supported
-    first(indices) > last(indices) && throw(ArgumentError("BinaryOperation operand indices $(a) and $(b) do not intersect!"))
-    return indices
+    # Check validity of the intersection index range
+    first(range_intersection) > last(range_intersection) &&
+        throw(ArgumentError("Indices $(from_idx) and $(to_idx) interpolated from $(from_loc) to $(to_loc) do not intersect!"))
+
+    return range_intersection
 end
 
-# Windowed fields interpolated from `Center`s to `Face`s lose the first index.
-# Viceverse, windowed fields interpolated from `Face`s to `Center`s lose the last index
-corrected_index(a, ::Face,   ::Face)   = UnitRange(first(a),   last(a)  )
-corrected_index(a, ::Center, ::Center) = UnitRange(first(a),   last(a)  )
-corrected_index(a, ::Face,   ::Center) = UnitRange(first(a),   last(a)-1)
-corrected_index(a, ::Center, ::Face)   = UnitRange(first(a)+1, last(a)  )
+validate_shifted_index(shifted_idx) = first(shifted_idx) > last(shifted_idx) &&
+    throw(ArgumentError("Cannot compute index intersection for indices $(from_idx) interpolating from $(from_loc) to $(to_loc)!"))
+
+"""
+    restrict_index_for_interpolation(from_idx, from_loc, to_loc)
+
+Return a "restricted" index range for the result of interpolating from
+`from_loc` to `to_loc`, over the index range `from_idx`:
+
+* Windowed fields interpolated from `Center`s to `Face`s lose the first index.
+* Conversely, windowed fields interpolated from `Face`s to `Center`s lose the last index
+"""
+restrict_index_for_interpolation(from_idx, ::Type{Face},   ::Type{Face})   = UnitRange(first(from_idx),   last(from_idx))
+restrict_index_for_interpolation(from_idx, ::Type{Center}, ::Type{Center}) = UnitRange(first(from_idx),   last(from_idx))
+restrict_index_for_interpolation(from_idx, ::Type{Face},   ::Type{Center}) = UnitRange(first(from_idx),   last(from_idx)-1)
+restrict_index_for_interpolation(from_idx, ::Type{Center}, ::Type{Face})   = UnitRange(first(from_idx)+1, last(from_idx))

src/AbstractOperations/binary_operations.jl

@@ -29,7 +29,7 @@ end
 # Recompute location of binary operation
 @inline at(loc, β::BinaryOperation) = β.op(loc, at(loc, β.a), at(loc, β.b))
 
-indices(β::BinaryOperation) = interpolate_indices(β.a, β.b; loc_operation = location(β))
+indices(β::BinaryOperation) = construct_regionally(intersect_indices, location(β), β.a, β.b)
 
 """Create a binary operation for `op` acting on `a` and `b` at `Lc`, where
 `a` and `b` have location `La` and `Lb`."""

src/AbstractOperations/computed_field.jl

@@ -44,7 +44,7 @@ function Field(operand::AbstractOperation;
     loc = location(operand)
     indices = validate_indices(indices, loc, grid)
 
-    boundary_conditions = FieldBoundaryConditions(indices, boundary_conditions)
+    @apply_regionally boundary_conditions = FieldBoundaryConditions(indices, boundary_conditions)
 
     if isnothing(data)
         data = new_data(grid, loc, indices)
@@ -63,7 +63,7 @@ Compute `comp.operand` and store the result in `comp.data`.
 """
 function compute!(comp::ComputedField, time=nothing)
     # First compute `dependencies`:
-    @apply_regionally compute_at!(comp.operand, time)
+    compute_at!(comp.operand, time)
 
     # Now perform the primary computation
     @apply_regionally compute_computed_field!(comp)

src/AbstractOperations/conditional_operations.jl

@@ -11,10 +11,10 @@ struct ConditionalOperation{LX, LY, LZ, O, F, G, C, M, T} <: AbstractOperation{L
     condition :: C
     mask :: M
 
-     function ConditionalOperation{LX, LY, LZ}(operand::O, func::F, grid::G, condition::C, mask::M) where {LX, LY, LZ, O, F, G, C, M}
-         T = eltype(operand)
-         return new{LX, LY, LZ, O, F, G, C, M, T}(operand, func, grid, condition, mask)
-     end
+    function ConditionalOperation{LX, LY, LZ}(operand::O, func::F, grid::G, condition::C, mask::M) where {LX, LY, LZ, O, F, G, C, M}
+        T = eltype(operand)
+        return new{LX, LY, LZ, O, F, G, C, M, T}(operand, func, grid, condition, mask)
+    end
 end
 
 """

src/AbstractOperations/kernel_function_operation.jl

@@ -71,7 +71,7 @@ struct KernelFunctionOperation{LX, LY, LZ, G, T, K, D} <: AbstractOperation{LX,
 end
 
 @inline Base.getindex(κ::KernelFunctionOperation, i, j, k) = κ.kernel_function(i, j, k, κ.grid, κ.arguments...)
-indices(κ::KernelFunctionOperation) = interpolate_indices(κ.arguments...; loc_operation = location(κ))
+indices(κ::KernelFunctionOperation) = construct_regionally(intersect_indices, location(κ), κ.arguments...)
 compute_at!(κ::KernelFunctionOperation, time) = Tuple(compute_at!(d, time) for d in κ.arguments)
 
 "Adapt `KernelFunctionOperation` to work on the GPU via CUDAnative and CUDAdrv."

src/AbstractOperations/multiary_operations.jl

@@ -20,7 +20,7 @@ end
 ##### MultiaryOperation construction
 #####
 
-indices(Π::MultiaryOperation) = interpolate_indices(Π.args...; loc_operation = location(Π))
+indices(Π::MultiaryOperation) = construct_regionally(intersect_indices, location(Π), Π.args...)
 
 function _multiary_operation(L, op, args, Largs, grid)
     ▶ = Tuple(interpolation_operator(La, L) for La in Largs)

src/Advection/weno_interpolants.jl

@@ -328,7 +328,7 @@ julia> calc_weno_stencil(2, :right, :x)
                                     dir == :y ?
                                     :(@inbounds ψ(i, j + $c, k, args...)) :
                                     :(@inbounds ψ(i, j, k + $c, args...))
-            else    
+            else
                 stencil_point[idx] =  dir == :x ? 
                                     :(@inbounds ψ[i + $c, j, k]) :
                                     dir == :y ?

src/BoundaryConditions/boundary_condition.jl

@@ -81,8 +81,8 @@ const OBC  = BoundaryCondition{<:Open}
 const VBC  = BoundaryCondition{<:Value}
 const GBC  = BoundaryCondition{<:Gradient}
 const ZFBC = BoundaryCondition{Flux, Nothing} # "zero" flux
-const MCBC  = BoundaryCondition{<:MultiRegionCommunication}
-const DCBC  = BoundaryCondition{<:DistributedCommunication}
+const MCBC = BoundaryCondition{<:MultiRegionCommunication}
+const DCBC = BoundaryCondition{<:DistributedCommunication}
 
 # More readable BC constructors for the public API.
                 PeriodicBoundaryCondition() = BoundaryCondition(Periodic,                 nothing)