Internal APIs

return_eltype(::Op, fields...)

Return the element type of the Field generated by an AbstractOperator when applied to the given fields.

return_space(::Op, spaces...)

Return the space (axes) of the Field generated by an AbstractOperator when applied to fields on the given spaces.

stencil_interior_width(::Op, args...)

Return the width of the interior stencil for the operator Op with the given arguments, as a tuple of 2-tuples: each 2-tuple holds the lower and upper bounds of the index offsets of the stencil for the corresponding argument.

Examples

stencil_interior_width(::Op, arg1, arg2) = ((-half, 1 + half), (0, 0))

implies that at index i, the stencil accesses arg1 at i - half, i + half and i + 1 + half, and arg2 at index i.

stencil_interior(::Op, space, idx, hidx, args...)

Return the value of the interior stencil of the operator Op at vertical index idx and horizontal index hidx of space; args are the input arguments.

boundary_width(::Op, ::BC, args...)

Return the width of a boundary condition BC on an operator Op: the number of locations at which a modified stencil is used. Either this function, or left_interior_idx and right_interior_idx, must be defined for a specific Op/BC combination.

stencil_left_boundary(op, bc, space, idx, hidx, args...)

Return the result of the stencil operator op at horizontal index hidx and vertical index idx of space near the left boundary, with boundary condition bc. For operators that cannot be evaluated without a boundary condition, a NullBoundaryCondition generates NaN values here.

Operators that are rewritten into an operator matrix multiply do not reach this method: their boundary rows come from MatrixFields instead, where a NullBoundaryCondition row is filled with NaNs, so the boundary output is NaN there as well.

stencil_right_boundary(op, bc, space, idx, hidx, args...)

Return the result of the stencil operator op at horizontal index hidx and vertical index idx of space near the right boundary, with boundary condition bc. For operators that cannot be evaluated without a boundary condition, a NullBoundaryCondition generates NaN values here.

Operators that are rewritten into an operator matrix multiply do not reach this method: their boundary rows come from MatrixFields instead, where a NullBoundaryCondition row is filled with NaNs, so the boundary output is NaN there as well.

left_interior_idx(space, op::FiniteDifferenceOperator, bc::VerticalBoundaryCondition, args...)

Return the index of the left-most interior point of the operator op with boundary bc when used with arguments args.... By default, this is

left_idx(space) + boundary_width(op, bc)

but it can be overridden for specific stencil types (e.g. if the stencil is asymmetric).

right_interior_idx(space, op::FiniteDifferenceOperator, bc::VerticalBoundaryCondition, args...)

Return the index of the right-most interior point of the operator op with boundary bc when used with arguments args.... By default, this is

right_idx(space) - boundary_width(op, bc)

but it can be overridden for specific stencil types (e.g. if the stencil is asymmetric).

fd_shmem_is_supported(bc::Base.Broadcast.AbstractBroadcasted)

Return whether the broadcasted object bc supports shared memory, so that the GPU extension can dispatch to an optimized kernel.

This function and its dispatch become unnecessary once all operators support shared memory.

AbstractOperator

Supertype for ClimaCore operators. An operator is a pseudo-function: it cannot be called directly, but can be broadcasted over Fields. Subtypes are FiniteDifferenceOperator and SpectralElementOperator.

SpectralElementOperator

Operator applied to the quadrature points in each spectral element of a Field. Each subtype must define return_eltype and apply_operator.

apply_operator(op, args...)

Eagerly evaluate the SpectralElementOperator op over one slab of slab Fields and/or lazy pointwise broadcasts over slabs.

Notes

Buffer reuse invariant. On GPUs, buffers of equal byte size can be assigned to the same shared memory (see DataLayouts.scoped_static_array), so each apply_operator method keeps at most one buffer of each size live at a time, separating every reuse from the previous use with a DataLayouts.synchronize. The destination stays live for the whole application and may alias a different live application's destination (as in @. wdiv(grad(a)) + wdiv(grad(b))), so it lives in per-thread registers (register_similar), as do temporaries materialized inside a fused slice loop; materialize_buffer republishes register-resident values whenever an operator must read them across threads. Two configurations rely on their allocation sites compiling into one unit (equal-size globals merge only across separately compiled functions): interpolation to an equal or lower degree gives a buffered argument and sequential_muladd_slab!'s partial result the same byte size, and SplitDivergence's buffered second argument stays live across its equal-size per-dimension intermediates. Both are pinned by GPU value tests in test/Operators/spectralelement/gpu_rectilinear.jl.

register_similar(arg, T)

Allocate a Field like Base.similar(arg, T), but with its data in each thread's registers; used for every apply_operator destination, which only its own thread reads and writes. This lets two applications in one fused expression be live at once (see the buffer reuse invariant in apply_operator).

buffer_similar(arg, T)

Allocate a Field like Base.similar(arg, T), but always through the argument's DataLayouts.DataScope (shared memory on GPUs), never in per-thread registers; used for every buffer whose values cross a thread boundary.

materialize_buffer(arg)

Materialize arg like Base.materialize(arg), but store the result of a lazy Broadcasted expression in a buffer from buffer_similar that every thread in the argument's scope can read; register-resident Fields (see register_similar) are also copied into such a buffer, and other arguments are returned unchanged. On GPUs two buffers of equal byte size can share memory, so callers must keep their lifetimes disjoint; see the buffer reuse invariant in apply_operator.

DGConnectivity

Cached, device-resident connectivity and face geometry for the DG interior-face operators (the DSS-buffer analog for DG).

Fields

• nfaces: Number of faces.
• nbnodes: Number of element boundary nodes in the gather map.
• faces: 5 × nfaces Int32 matrix of interior faces (elem⁻, face⁻, elem⁺, face⁺, reversed).
• sgeom: Precomputed Geometry.SurfaceGeometry per (level, q, face) (level = 1 for pure 2D spaces), evaluated from the minus side exactly as the CPU loops do.
• node_elem, node_i, node_j, node_offset, contrib_face, contrib_side, contrib_q: A deterministic gather map from element boundary nodes to their face contributions, in ragged-array form: each boundary node (elem, i, j) lists the (face, side, q) face-node slots that accumulate into it (2 entries at element corners, 1 elsewhere), sorted at construction so the GPU gather is bitwise deterministic.

Built once per space by dg_connectivity and stored with the array type of the space's device (ClimaComms.array_type). dg_ghost_connectivity builds the same structure over the ghost (inter-rank) faces, with the meanings documented there.

dg_connectivity(space)

Return the memoized DGConnectivity for space, keyed on the underlying grid and the space type (so center/face extruded spaces get separate buffers). Stored in Utilities.Cache.OBJECT_CACHE alongside the grid objects, so Utilities.Cache.clean_cache! releases it (the buffer holds device arrays).

dg_ghost_connectivity(space)

Return the memoized DGConnectivity over the topology's ghost (inter-rank) faces, or nothing when there are none. The entries differ from dg_connectivity in two ways: the third faces row holds the strip slot into the recv buffer of a Topologies.GhostFaceExchange (the plus-side value at loop node q is strip node reversed ? Nq - q + 1 : q) rather than a local element, and the gather map contains only minus-side contributions (side == 1) — the mirror ghost face on the neighbouring rank accumulates the other side (see _add_dg_ghost_faces! for the operator contract that makes the two sides consistent).

RegisterArray{Sz, F, Stride}(array)

Array that presents the full size Sz of a DataScope's data while only storing the entries that belong to the thread reading it, in a StaticArrays.MArray or SArray that a compiler can promote to registers. F is the position of the field dimension (see f_dim), Stride is the number of threads in the scope it was allocated for (also the stride between consecutive points of the reading thread), and register_array_params derives the remaining size parameters from Sz and F.

Every thread is assigned the strided subset rank:Stride:Np of the scope's points (see subscope_indices), so the stored array has the compile-time size cld(Np, Stride) * Nf. When a thread holds a single point (the usual case; see slice_subscope), a component's index does not depend on the point, so every component can live in a register.

Reading or writing a point that belongs to a different thread silently accesses the reading thread's own point at the same subset position, so a RegisterArray may only hold data that no other thread reads; values that cross threads must be published through scoped_static_array buffers. Unlike shared memory globals, identified by byte size alone (see the equal sizes invariant in ext/cuda/scopes.jl), register arrays are distinct values, which keeps two live operator destinations from aliasing.

register_similar(data, T)
register_similar(bc, T)

Allocate a DataLayout like Base.similar(data, T), but backed by a RegisterArray that only stores each thread's own points; the result must not be read by any other thread (see RegisterArray). This falls back to buffer_similar when registers cannot be used (a non-inferrable size, a non-constant thread count, or a single-thread scope, which is already an MArray); lazy layouts are first converted into the layout type that buffer_similar would allocate.

buffer_similar(data, T)
buffer_similar(bc, T)

Allocate a new DataLayout through its DataScope's memory (see scoped_array and scoped_static_array), so that every thread in the scope can read the result. Base.similar on a LazyDataLayout instead routes through register_similar, whose result may only be read by the writing thread, so any buffer whose values cross a thread boundary has to be allocated with this function.

static_num_threads(scope)

Return the num_threads of a DataScope when it is a compile-time constant, and nothing when it is only known at run time. Used by register_similar to determine how many points of a slice each thread can be assigned.

bilinear_interpolate(coords::NTuple{4}, ξ1, ξ2)

Interpolate bilinearly between coords with the reference parameters ξ1, ξ2, each in the interval [-1, 1].

coords is a 4-tuple of coordinates (SVectors or points) in counter-clockwise order.

      4-------3
 ^    |       |
 |    |       |
ξ2    |       |
      1-------2

        ξ1-->

mul_with_projection(x, y, lg)

Compute x * y, projecting y first so that Tensor operands do not raise DimensionMismatch errors. For example, if x is a covector along the Covariant3Axis (e.g., Covariant3Vector(1)'), then y is projected onto the Contravariant3Axis. In general, the first axis of y is projected onto the dual of the last axis of x, using the local geometry lg. When x is not a Tensor{2}, this is x * y.

mul_return_type(X, Y)

Return the type of mul_with_projection(x, y, lg) for x::X and y::Y. Equivalent to Base._return_type(mul_with_projection, Tuple{X, Y, LG}), but explicit so that the eltype inference in MatrixFields always sees a concrete type; the internal _return_type can widen to Union or Any and is unstable across Julia versions.

The methods cover six result shapes (scalar×scalar, scalar×tensor, covector×vector, vector×covector, 2-tensor×vector, 2-tensor×2-tensor) in which the output ndims goes up, down, or stays the same; no single formula fits all of them. The fallback method throws an error for unsupported type pairs.

Future cleanup: try replacing with Base._return_type and keep only methods that fail; collapse the two Tensor{2}*Tensor{N} cases via Base.tail(axes(Y)); move tensor_type/basis1/basis2 to tensors.jl.

Meshes.SharedVertices(mesh, elem, vert)

Iterator over the (element, vertex) pairs that share the vertex vert of element elem in mesh, starting with (elem, vert) itself.

Meshes.face_connectivity_matrix(mesh, elemorder = elements(mesh))

Construct a Bool-valued SparseMatrixCSC M containing the face connections of mesh. Elements are indexed according to elemorder.

M[i, i] == true only if two distinct faces of element i are connected.

Meshes.vertex_connectivity_matrix(mesh, elemorder = elements(mesh))

Construct a Bool-valued SparseMatrixCSC M containing the vertex connections of mesh. Elements are indexed according to elemorder.

M[i, i] == true only if two distinct vertices of element i are connected.

opelem, opface, reversed = Meshes.opposing_face(mesh::AbstractMesh, elem, face::Int)

Return the element and face (opelem, opface) opposite face face of element elem, and whether the node ordering along the shared face is reversed relative to face.

Meshes.is_boundary_face(mesh::AbstractMesh, elem, face::Int)::Bool

Return true if face face of element elem is on the boundary of mesh.

elem is an element of elements(mesh).

Meshes.boundary_face_name(mesh::AbstractMesh, elem, face::Int)::Union{Symbol,Nothing}

Return the name of the boundary containing face face of element elem, or nothing if the face is not on the boundary.

GhostFaceExchange

Face-strip halo exchange for the DG ghost-face operators: exchanges only the Nq face-node values (per level and field component) of each rank-boundary face, instead of whole neighbour elements — the face analog of the perimeter-only DSSBuffer. Only face-sharing neighbours participate; vertex-only halo neighbours are excluded on both sides of each pair, so participation is symmetric by construction.

send_data/recv_data are (Nv, Nq, 1, nstrips) layouts with one strip per entry of ghost_faces: slot s of send_data holds this rank's own ("minus") side of a shared face, and slot s of recv_data receives the neighbouring rank's side of the same face. Both ranks assign slots by sorting their shared faces by neighbour pid and then by the canonical face key (the ordered pair of (global element, face) of the two sides), which both ranks can compute, so the k-th strip sent within a pair is the k-th strip the peer expects. Strips are packed in the sending face's natural node order (face_node_index(face, Nq, q, false)); a receiver whose face is reversed relative to the sender reads node q at strip index Nq - q + 1.

Fields

• graph_context: The ClimaComms exchange over face-sharing neighbours.
• send_data, recv_data: The strip layouts described above.
• slot_lidx, slot_face: Strip slot → local element and face, stored with the topology's array type for the device-side pack.
• face_slot: Position in ghost_faces order → strip slot (host vector, consumed when building face connectivity on the host).
• in_flight: Latch set at fill and cleared at finish. Exchanges are memoized per (space, data type, argument position), so distinct fields of the same type share this object; an overlapping second start is an error rather than a silent overwrite of the in-flight send strips.

nsendelems(topology)

Return the number of elements this process sends to its neighbors in topology.

nghostelems(topology)

Return the number of ghost elements in topology.

localelemindex(topology, elem)

Return the local index of element elem; used by distributed topologies.

face_node_index(face, Nq, q, reversed = false)

Return the node indices (i, j) of the qth node on face face, where Nq is the number of nodes in each direction. If reversed, count from the other end of the face.

ghost_faces(topology::AbstractTopology)

Return an iterator over the ghost faces of topology. Each item is a 5-tuple of the form

(elem1, face1, elem2, face2, reversed)

where elemX, faceX are the element and face numbers, and reversed indicates whether they have opposing orientations.

vertex_node_index(vertex_num, Nq)

Return the node indices (i, j) of vertex vertex_num, where Nq is the number of nodes in each direction.

ghost_vertices(topology)

Return an iterator over the ghost vertices of topology. Each vertex is an iterator over (isghost, idx, vert) triples, where idx is a local index lidx if isghost is false and a receive buffer index ridx otherwise.

ghost_neighboring_elements(topology::AbstractTopology, lidx::Integer)

Return an iterator over the receive buffer indices (ridx) of the ghost elements that are neighbors of the local element lidx in topology.

dss_transform(arg, local_geometry, weight, I)

Transform arg[I] to a basis for direct stiffness summation (DSS) and multiply it by the DSS weight. Transformations only apply to vector quantities.

• local_geometry[I]: The LocalGeometry at I. If local_geometry is nothing, no transformation is performed.
• weight[I]: The DSS weight at I. If weight is nothing, the weight is 1 and the result is plain summation.

See ClimaCore.Spaces.weighted_dss!.

dss_transform!(device, dss_buffer, data, local_geometry, dss_weights, perimeter, localelems)

Transform vectors in data from covariant/contravariant axes to physical axes, weight the data at perimeter nodes, and store the result in dss_buffer.perimeter_data.

Arguments

• dss_buffer: DSSBuffer created by create_dss_buffer for data.
• data: Field data.
• local_geometry: Local metric information at each node.
• dss_weights: DSS weights of the horizontal space.
• perimeter: Perimeter iterator.
• localelems: Local elements on which to perform the transformation.

Part of ClimaCore.Spaces.weighted_dss!.

dss_untransform(T, targ, local_geometry, I)

Transform targ back to a value of type T using local_geometry[I], after direct stiffness summation (DSS).

See ClimaCore.Spaces.weighted_dss!.

dss_untransform!(device, dss_buffer, data, local_geometry, perimeter, localelems)

Transform physical vectors in dss_buffer.perimeter_data back to their original covariant/contravariant axes, and store the result in data.

Arguments

• dss_buffer: DSSBuffer created by create_dss_buffer for data.
• data: Field data.
• local_geometry: Local metric information at each node.
• perimeter: Perimeter iterator.
• localelems: Local elements on which to perform the transformation.

Part of ClimaCore.Spaces.weighted_dss!.

compute_element_bounds!(limiter::QuasiMonotoneLimiter, ρq, ρ)

Compute the min and max of q = ρq / ρ over the nodes of each element and store them in limiter.q_bounds.

Part of compute_bounds!.

compute_neighbor_bounds_local!(limiter::QuasiMonotoneLimiter, ρ)

Set limiter.q_bounds_nbr in each element to the min and max of limiter.q_bounds over the element and its process-local neighbors in the topology of axes(ρ).

Part of compute_bounds!.

compute_neighbor_bounds_ghost!(limiter::QuasiMonotoneLimiter, topology)

Widen limiter.q_bounds_nbr in each element with the min and max of limiter.q_bounds in its ghost neighbors, read from limiter.ghost_buffer.recv_data. Call it after the ghost exchange has completed; it does nothing when limiter.ghost_buffer is not a GhostBuffer.

Part of compute_bounds!.

apply_limit_slab!(slab_ρq, slab_ρ, slab_WJ, slab_q_bounds, rtol)

Limit the nodal values of one element's tracer density slab_ρq in place so that q = ρq / ρ lies within slab_q_bounds, the min and max for the element, given the density slab_ρ and the quadrature weights times Jacobians slab_WJ.

The bounds are first widened to include the element mean of q, so that a solution exists. Values outside the bounds are clipped, and the resulting tracer mass change is redistributed over the nodes not at a bound, in proportion to slab_ρ * slab_WJ, until the relative mass change is at most rtol or Nq^2 iterations have been done. Each component of a struct-valued slab_ρq is limited independently against the matching component of slab_q_bounds.

Returns

A tuple (converged, max_rel_err, min_tracer_mass): whether all components met the tolerance, the largest relative mass error, and the smallest absolute tracer mass over the components.

column_massborrow!(
    q_data::AbstractArray,
    ρ_data::AbstractArray,
    ΔV_data::AbstractArray,
    q_min::AbstractFloat,
)

Apply the vertical mass-borrowing limiter in place to q_data, the array backing a single column of scalar tracer data, with column density ρ_data, cell volume ΔV_data, and minimum value q_min. Returns nothing.

default_target_hcoords(space::Spaces.AbstractSpace; hresolution = 180)

Return an array of Geometry.Points that cover the horizontal domain of space uniformly with hresolution points per direction.

On the sphere, the result is a hresolution × hresolution array of LatLongPoints, with latitudes from -90 to 90 along the first dimension and longitudes from -180 to 180 along the second (in degrees). On a plane, the points have the coordinate type of the horizontal domain and span its extent. Return nothing for spaces without a horizontal direction to interpolate over (FiniteDifferenceSpace, MultiColumnFiniteDifferenceSpace, MultiPointSpace).

default_target_zcoords(space::Spaces.AbstractSpace; zresolution = nothing)

Return a vector of Geometry.ZPoints covering the vertical extent of space.

When zresolution is nothing, return the cell-center heights of the model levels, so that vertical interpolation of a center field reproduces its values. Otherwise, return zresolution uniformly spaced heights between Domains.z_min(space) and Domains.z_max(space). Return nothing for spaces without a vertical direction (AbstractSpectralElementSpace, MultiPointSpace).

@drop_recursion_limits f₁, f₂, ...

Remove the inference recursion limit from every listed function's methods that this package owns, along with the package's Core.kwcall methods that target a listed function (or all of them, if Core.kwcall is listed) and their hidden body functions, where keyword arguments put the loop bodies. Listing a module applies to every function it defines. Use this for functions that recurse by design: the default limit widens their argument types, turning downstream calls into dynamic dispatch with runtime allocations, including inside GPU kernels. Only owned methods are touched, since a mutated foreign method would not survive precompilation.

stable_view(array, indices...)

Return a view of array like view(array, indices...), with two modifications that avoid expensive operations:

• Every view is a SubArray, even when array is a GPU array. GPUArrays replaces each contiguous view of a CuArray with a new CuArray derived from the same memory buffer, and the derived array's type is not inferrable, which makes all host code that builds slice or property views type-unstable. The SubArrays constructed here have fully inferred types, and they are converted to SubArrays of CuDeviceArrays when passed to kernels.
• A view along the linear indices of a multidimensional array (a single Integer or range of Integers) wraps the array in a 1-dimensional ReshapedArray, instead of using reshape like Base's view does, which allocates a new object whenever it is applied to an Array. If the array is already a ReshapedArray, its parent gets wrapped instead, since a reshape stores the same values in the same linear order as its parent.

julia> array = rand(3, 1, 4);

julia> parent(view(array, 4:6))
12-element Vector{Float64}

julia> parent(stable_view(array, 4:6))
12-element reshape(::Array{Float64, 3}, 12) with eltype Float64

unionall_type(T)

Return the type T with all of its parameters dropped. If the argument is not a Type, its type is used instead.

Examples

julia> unionall_type(typeof([1, 2, 3]))
Array

julia> unionall_type((; a = 1, b = 2))
NamedTuple

replace_type_parameter(T, P, P′)

Return the type T with every parameter that is a subtype of P recursively replaced by P′. This is like constructing a value of type T and converting subfields of type P to type P′, though no constructors are called or compiled.

fieldtype_vals(T)

Return a statically inferrable analogue of Val.(fieldtypes(T)). Functions of Types are specialized upon successful constant propagation, but functions of Vals are always specialized, so fieldtype_vals can be used in place of fieldtypes to ensure that recursive functions over nested types have inferrable outputs.

new(T, [fields])

Allocate a value of type T with the new pseudo-function, which can otherwise only be called from inner constructors.

If provided, the second argument is used to initialize fields of the new value (unlike the lowered pseudo-function, this does not convert the values to the fieldtypes of T). Otherwise, the fields are initialized with arbitrary data, with special handling of DataType fields to avoid errors during compilation.

Examples

julia> new(Int)
4889520192

julia> new(Complex{Int}, (1, 2))
1 + 2im

julia> new(@NamedTuple{a::Type{Int}, b::Int, c::Complex{Int}})
(a = Int64, b = 4889520192, c = 6162822528 + 8036417625im)

julia> new(@NamedTuple{a::DataType, b::Int, c::Complex{Int}}, (Int, 1, 1 + 2im))
(a = Int64, b = 1, c = 1 + 2im)

is_inferred_type(T)

Return whether T either satisfies isconcretetype or is a Type{..} value (or the more generic DataType value).

return_type(f, T)

Return Core.Compiler.return_type(f, T), after checking that the result satisfies is_inferred_type whenever T does; otherwise throw an InferenceError. Used in place of Core.Compiler.return_type to flag deteriorations in type inference before they can lead to behavioral changes.

unsafe_eltype(itr)

Return the element type of itr, including un-materialized broadcast expressions, adapted from Base.Broadcast.combine_eltypes. This performs no safety checks, and it may return non-concrete types (like an empty Union{}).

safe_eltype(itr)

Return the element type of itr, including un-materialized broadcast expressions, adapted from Base.Broadcast.combine_eltypes. Throw an error when the result does not satisfy is_inferred_type, indicating which part of the expression first encounters a type instability or an error during inference.

safe_mapreduce(f, op, itr; [init])

Apply f to every element of the indexable collection itr and reduce the results with op, like Base.mapreduce(f, op, itr; init), with the guarantee that the call can be compiled in GPU kernels. Unlike Base.mapreduce, this never reaches the empty collection error path, which builds a string that cannot be compiled for a GPU; if init is not provided, itr must be nonempty.

When init is provided, the reduction is a sequential left fold seeded by init. Otherwise, the reduction is only sequential for iterator lengths below 1024; for longer iterators, the reduction is pairwise, so its roundoff error grows logarithmically rather than linearly with length. Sequential reductions use @simd loops, so the order in which op is applied is not guaranteed.

ConvertTo{T}()

GPU-compatible callable that converts its argument to type T, equivalent to Base.Fix1(convert, T) but isbitstype. Base.Fix1 stores a Type{T} field, which is not isbits, so it cannot be captured by GPU kernels. ConvertTo{T} is an empty struct and is always isbits, so it can be used in broadcast expressions that run on the GPU.

Examples

julia> isbitstype(typeof(ConvertTo{Float32}()))
true

julia> isbitstype(typeof(Base.Fix1(convert, Float32))) # cannot enter a kernel
false

AutoBroadcaster(itr)

Wrapper for an iterator that forces certain functions to be broadcasted over the iterator's elements. This allows different types of broadcasting to be applied simultaneously; e.g., ClimaCore's Fields and similar types use the standard dot syntax to denote parallelized iteration over spatial locations, and they wrap their values in AutoBroadcasters for unrolled iteration over subfields. All statically-sized iterators for which is_auto_broadcastable is true are compatible with AutoBroadcasters.

In the context of AutoBroadcasters, broadcasting a function applies it with unrolled_map, iterating over all arguments for which is_auto_broadcastable is true (including those not wrapped in AutoBroadcasters), while other arguments are passed to the function directly. This behavior is triggered by using AutoBroadcasters, optionally in conjunction with compatible iterators that are not wrapped in AutoBroadcasters, in the following ways:

• passing them to standard math functions or constructors,
• passing them to ifelse (for iterating over conditional values),
• applying them as function calls (for iterating over functions),
• explicitly calling nested_broadcast.

Nested AutoBroadcasters constructed with add_auto_broadcasters evaluate broadcasts recursively, mapping across every layer of nested iterators so that broadcasted functions are only applied to non-iterators in the innermost layers. Aside from automatic broadcasting, AutoBroadcasters behave like their underlying iterators, with support for common operations like iterate, propertynames, getindex, and reduce.

Examples

julia> import ClimaCore.Utilities, StaticArrays;

julia> x = Utilities.AutoBroadcaster((1, 2.0, StaticArrays.SVector(3, 4)))
(1, 2.0, [3, 4])

julia> zero(typeof(x))
(0, 0.0, [0, 0])

julia> 2 * x - (2, 3, [4, 5])
(0, 1.0, [2, 3])

julia> y = Utilities.add_auto_broadcasters((1, 2, (a = 3, b = 4, c = (5, 6, (7, 8)))))
(1, 2, (a = 3, b = 4, c = (5, 6, (7, 8))))

julia> min(y, abs(5 - y))
(1, 2, (a = 2, b = 1, c = (0, 1, (2, 3))))

julia> x' * y * x ÷ 5
(0, 1.0, (a = 15, b = 20, c = (25, 30, (35, 40))))

is_auto_broadcastable(T)
is_auto_broadcastable(itr)

Return whether an AutoBroadcaster should broadcast over iterators of type T. By default, this is only true for Tuple and NamedTuple types, but it can be extended to any statically sized type compatible with UnrolledUtilities.jl.

For convenience, is_auto_broadcastable also accepts an iterator instead of its type, but this method should not be extended directly.

add_auto_broadcasters(itr)
add_auto_broadcasters(T)

Recursively apply the AutoBroadcaster constructor to iterators for which is_auto_broadcastable is true, as well as their elements for which it is true, while leaving values for which it is false unmodified. When passed an iterator's type T, return the inferred result type for such an iterator.

drop_auto_broadcasters(itr)
drop_auto_broadcasters(T)

Recursively undo add_auto_broadcasters, extracting the iterator from every AutoBroadcaster in itr. When passed an iterator's type T, return the inferred result type for such an iterator.

auto_broadcasted([style], f, args, [axes])

Construct a Base.Broadcast.Broadcasted(style, f, args, axes), passing the arguments of f through either add_auto_broadcasters or drop_auto_broadcasters if doing so avoids an inferred error.

When the unsafe_eltype of Broadcasted(style, f, args, axes) indicates that f throws an error, a new Broadcasted wrapper is constructed with add_auto_broadcasters applied to every argument, and then another is constructed with drop_auto_broadcasters applied to every argument. If one of the new wrappers no longer corresponds to a guaranteed error, it is returned instead of the original wrapper. Otherwise, the default result of Broadcasted(style, f, args, axes) is returned without modifications.

Examples

julia> import ClimaCore.Utilities;

julia> x = (im, (1, 2.0), [3, 4])
(im, (1, 2.0), [3, 4])

julia> y = [x, x, x, x];

julia> bc = Base.Broadcast.Broadcasted(*, (Base.Broadcast.Broadcasted(adjoint, (y,)), y));

julia> sum(Base.materialize(bc))
ERROR: MethodError: no method matching adjoint(::Tuple{...})
[...]

julia> bc = Utilities.auto_broadcasted(*, (Utilities.auto_broadcasted(adjoint, (y,)), y));

julia> sum(Base.materialize(bc))
(4 + 0im, (4, 16.0), 100)

nested_broadcast(f, args...)

Apply f elementwise to args, like broadcast, recursing over nested iterators, as long as at least one argument is an AutoBroadcaster. All loops over iterator elements are unrolled and inlined.

This function is called automatically when an AutoBroadcaster is passed to any standard math function or constructor, but for generic operations it must be called explicitly.

Examples

julia> import ClimaCore.Utilities;

julia> x = Utilities.add_auto_broadcasters(((:a, :b, :c), (:d, :e, :f), :g))
((:a, :b, :c), (:d, :e, :f), :g)

julia> Utilities.nested_broadcast(string, x)
(("a", "b", "c"), ("d", "e", "f"), "g")

julia> y = Utilities.add_auto_broadcasters((1, 11, (111, 1111, 11111)))
(1, 11, (111, 1111, 11111))

julia> Utilities.nested_broadcast(Symbol, x, y * y)
((:a1, :b1, :c1), (:d121, :e121, :f121), (:g12321, :g1234321, :g123454321))

outer_diagonals(::Type{<:BandMatrixRow})

Return the indices of the lower and upper diagonals, (ld, ud), of the given subtype of BandMatrixRow.

band_matrix_row_type(ld, ud, T)

Return the subtype of BandMatrixRow that has entries of type T on the diagonals with indices in the range ld:ud.

matrix_shape(matrix_field, [matrix_space])

Return the shape of a matrix field whose rows are defined on matrix_space, which defaults to axes(matrix_field).

When matrix_space is a finite difference space (extruded or not), the shape is Square(), FaceToCenter(), or CenterToFace(), depending on whether the diagonal indices of matrix_field are Ints or PlusHalfs and whether matrix_space is on cell centers or cell faces.

When matrix_space is a spectral element or point space, only the Square() shape is supported, and matrix_field must have DiagonalMatrixRow entries.

column_axes(matrix_field, [matrix_space])

Return the space that corresponds to the columns of matrix_field, i.e., the axes of the Fields by which matrix_field can be multiplied. The matrix_space is the space that corresponds to the rows of matrix_field, and it defaults to axes(matrix_field).

AbstractLazyOperator

Supertype for "lazy operators", i.e., operators that users can call without any arguments, as long as they appear in broadcast expressions that contain at least one Field. If lazy_op is an AbstractLazyOperator, the expression lazy_op.() is translated to non_lazy_op.(fields...) when it appears in a broadcast expression with at least one Field. This translation is done by the function replace_lazy_operator(space, lazy_op), which every subtype of AbstractLazyOperator must implement.

replace_lazy_operator(space, lazy_op)

Return an instance of Base.AbstractBroadcasted that corresponds to the expression lazy_op.(), where the broadcast in which this expression appears is evaluated on the given space. The staggering (CellCenter or CellFace) of this space depends on the specifics of the broadcast and is not predetermined.

is_lazy(dict)

Return whether the FieldNameDict dict has any AbstractBroadcasted entries, i.e. entries that have not been materialized.

lazy_main_diagonal(matrix)

Construct a lazy FieldMatrix with the diagonal keys of matrix, whose entries are the main diagonals of the corresponding entries of matrix. Diagonal entries are kept as they are; other ColumnwiseBandMatrixField entries become Broadcasted objects of DiagonalMatrixRows.

lazy_mul(A, args...)

Construct a lazy FieldMatrix that represents the product @. *(A, args...). This involves regular broadcasting when A is a FieldMatrix, but it has more complex behavior for other objects like the LazySchurComplement.

LazySchurComplement(A₁₁, A₁₂, A₂₁, A₂₂, [alg₁, cache₁, A₁₂_x₂, invA₁₁_A₁₂_x₂])

Analogue of a FieldMatrix that represents the Schur complement of A₁₁ in A, A₂₂ - A₂₁ * inv(A₁₁) * A₁₂. Since inv(A₁₁) is generally a dense matrix, computing the Schur complement directly is inefficient, so this object only supports the "lazy" functions lazy_mul, which multiplies it by the vector x₂, and lazy_preconditioner, which approximates it with a FieldMatrix.

The values alg₁, cache₁, A₁₂_x₂, and invA₁₁_A₁₂_x₂ must be specified for lazy_mul to compute inv(A₁₁) * A₁₂ * x₂. When a LazySchurComplement is not passed to lazy_mul, these values can be omitted.

field_matrix_solver_cache(alg, A, b)

Allocate and return the cache required by the FieldMatrixSolverAlgorithm alg to solve the equation A * x = b.

check_field_matrix_solver(alg, cache, A, b)

Check that the sparsity structure of A is supported by the FieldMatrixSolverAlgorithm alg, and that A is compatible with b in the equation A * x = b; throw an error otherwise.

run_field_matrix_solver!(alg, cache, x, A, b)

Set x to the value that solves the equation A * x = b using the FieldMatrixSolverAlgorithm alg and its cache.

solver_algorithm(P_alg)

Return a FieldMatrixSolverAlgorithm that can be used to solve P * x = b for x, where P is the preconditioner generated by the PreconditionerAlgorithm P_alg. If P_alg is nothing instead of a PreconditionerAlgorithm, or if P is a diagonal matrix (and no solver is required to invert it), this returns nothing.

lazy_preconditioner(P_alg, A)

Construct a lazy FieldMatrix (or a concrete one when possible) that approximates A according to the PreconditionerAlgorithm P_alg. If P_alg is nothing instead of a PreconditionerAlgorithm, this returns one(A).

preconditioner_cache(P_alg, A, b)

Allocate and return the cache required to solve the equation P * x = b, where P is the preconditioner generated by the PreconditionerAlgorithm P_alg for A. The cache is an empty NamedTuple when P is diagonal.

check_preconditioner(P_alg, P_cache, A, b)

Check that P is compatible with b in the equation P * x = b, where P is the preconditioner generated by the PreconditionerAlgorithm P_alg for A. If P_alg requires a FieldMatrixSolverAlgorithm alg to solve the equation, this also calls check_field_matrix_solver on alg.

lazy_or_concrete_preconditioner(P_alg, P_cache, A)

Return the preconditioner P from lazy_preconditioner, materialized into the concrete FieldMatrix P_cache.P when the PreconditionerAlgorithm P_alg requires a FieldMatrixSolverAlgorithm to invert it. Return nothing if P_alg is nothing.

apply_preconditioner(P_alg, P_cache, P, lazy_b)

Construct a lazy FieldVectorView (or a concrete one when possible) that represents the product @. inv(P) * b. Here, lazy_b is a (possibly lazy) FieldVectorView that represents b. When P_alg requires a solver, the result is stored in P_cache.x.

get_scalar_keys(dict::FieldMatrix)

Return the FieldMatrixKeys whose entries in dict are scalar blocks: keys for which dict[key] is a ScalingFieldMatrixEntry with values of type target_type or a ColumnwiseBandMatrixField whose band elements are of type target_type. target_type is the element type of the parent array of the first Field entry in dict, or Number when dict has no Field entries. Called from scalar_field_matrix.

get_scalar_keys(::Type{T}, ::Val{FT})

Return a tuple of the FieldNamePairs, relative to an entry with element type T, that address the components of T of type <: FT. A BandMatrixRow is indexed through its element type, the components of a Geometry.Tensor{2} are addressed by row and column index (a Geometry.Covector only by column index), and a rank-1 tensor is indexed through its components field.

field_offset_and_type(name_pair::FieldNamePair, ::Type{T}, ::Type{S}, full_key::FieldNamePair)

Return a tuple (offset, type, index_method) for the component of an object of type S at name_pair, where T is the base type of the parent array.

Returns

• offset: The position of the component within S, in multiples of sizeof(T).
• type: The type of the component.
• index_method: A Val naming how a Field with element type S can be indexed with name_pair:
  • Val(:view): a strided view of the parent array.
  • Val(:view_of_blocks): a view with non-uniform stride; not implemented, so it is handled like Val(:broadcasted_fallback).
  • Val(:broadcasted_fallback): a Broadcasted object that indexes each element.
  • Val(:broadcasted_zero): a Broadcasted object of zeros, because name_pair indexes an off-diagonal block under the implicit tensor structure optimization (see the MatrixFields reference page).

When S is a Geometry.Tensor{2} and only one name of name_pair indexes into its components, the pair addresses a slice, which cannot be a view, so the result is (0, S, Val(:broadcasted_fallback)). When neither name of name_pair is @name() and neither names a field of S, the pair is treated with the implicit tensor structure optimization: the first names of both are dropped, and the block is zero when they differ.

full_key is the complete key, used in the KeyError thrown when name_pair does not index into S. Called from get_internal_entry(::ColumnwiseBandMatrixField, ::FieldNamePair).