Guarantee uniqueness of topologies and spaces

[simonbyrne]

Is your feature request related to a problem? Please describe.

It's very easy to get the mismatched spaces warning, especially when loading data from multiple files. It will also become an issue

Related issues

• mismatched spaces #639 was the original issue, and was closed by Change same space different instance error to warning. #640 which turned it into a warning
• File format time series support #960 referred to this problem in the context of time series, but it is more general

Describe the solution you'd like

The basic idea would be to keep track of existing objects (spaces and below) in some sort of internal cache: if a constructor is called with the same arguments, we instead use the cached object.

The major challenge would be making sure that objects get GCed even if they're not referenced anywhere outside the cache. Julia provides WeakRefs for this purpose. See also WeakKeyDict (unfortunately we can't use a WeakKeyDict directly, as not all our objects are mutable).

We may want to change some objects to be mutable structs as that may make the table lookup easier.

Describe alternatives you've considered
The alternative is to allow more operations on duplicate objects, but this will be less efficient in terms of memory usage (as we will end up with duplicate objects) and performance (as we will need to spend more time checking objects are actually the same).

Additional context
Add any other context or screenshots about the feature request here.

[simonbyrne]

One place where this could cause problems is when we extract levels or columns of extruded spaces. Will need to think about it.

[simonbyrne]

After a bit more thought, I think the solution is:

• Our standard Topologies and Spaces are made to be mutable structs
  • levels and columns return distinct objects (LevelSpace, ColumnSpace) that are immutale wrappers which providing a view into the relevant object.
  • We will need some trimmed down space object that we can pass geometry info to the GPU (via Adapt.jl)
• For each Topology/Space type, we use a WeakValueDict as a cache to track inputs => created object
  • the default constructor will first check if the object is in the cache before constructing a new one
  • only the created objects need to be kept in a WeakRef

[simonbyrne]

There is an important design choice on how to handle staggered grids.

Currently we have a field staggering which details whether the field values will live on cell centers or cell faces, and we can construct one from the other by just creating a new object and switching that field. This will no longer work, as it will create new objects, and we will have no way to know if fields match up or not.

So some possible ways to do this:

1. We define one field as a wrapper around the other. Since most values are on cell centers, it would make sense to define the center field as the base case, and the face one as a wrapper:

   mutable struct CenterFiniteDifferenceSpace
    # all the information
   end
   
   # wrapper
   struct FaceFiniteDifferenceSpace{C<:CenterFiniteDifferenceSpace}
      center_space::C
   end

   and then we can switch between them by wrapping/unwrapping.

2. Define some intermediate object (say a "grid" for want of a better name), and make both spaces wrappers around this:

   mutable struct FiniteDifferenceGrid
    # all the information
   end
   
   struct CenterFiniteDifferenceSpace{G<:FiniteDifferenceGrid}
       grid::G
   end
   struct FaceFiniteDifferenceSpace{G<:FiniteDifferenceGrid}
       grid::G
   end

3. Move the staggering into the Field itself (this is what Oceananigans does: https://github.com/CliMA/Oceananigans.jl/blob/2270cc00c490c4d717c0b32b404320d6a8cbae75/src/Fields/field.jl#L122. I'm a bit reluctant to do this, as it would be a considerable breaking change.

[simonbyrne]

One advantage of option 2 is that it may make some aspects of the ClimaAtmos setup easier: cc @Sbozzolo

[charleskawczynski]

Maybe this is a silly question, but is it really a significant risk to simply remove

     if space1 !== space2 
         error_mismatched_spaces(typeof(space1), typeof(space2)) 
     end 

and let broadcasting over spaces be "unsafe"?

I would prefer option 2, though.

[glwagner]

I'd note that the motivation for option 3 is that this corresponds to how people typically communicate about staggered grids; ie there is a single primal grid or space (the tracer / fully centered grid), and the "other grids" that are staggered with respect to the tracer grid are (somehow, implicitly) diagnostic to that grid.

How is option 2 functionally different from 3? What else would the distinction between the two spaces be used for other than constructing fields?

[Sbozzolo]

One advantage of option 2 is that it may make some aspects of the ClimaAtmos setup easier: cc @Sbozzolo

I don't follow the nuances, but a comment against 2. is that there's already a tower of types needed to define something that is directly useful to the end user (Domain, Mesh, Topology, Space). In Atmos, I would like to provide something of immediate relevance to the end user (e.g., CubedSphere(radius, nh_poly, ...)), so I would probably end up packaging the FiniteDifferenceGrid in yet another struct.

[simonbyrne]

@charleskawczynski there are two problems we're trying to solve:

1. avoiding creating the same space multiple times: i.e. so that metric information isn't duplicated, and also make checking that spaces are the same is easier (you can just check if the two pointers are the same).
2. blow up in compilation time when using tuples of fields (Make Spaces mutable #1470)

[charleskawczynski]

I'd prefer option 2 over 1, just for symmetry sake. There's nothing special about center vs face grids, and as Greg was saying, one is a sort of dual grid to the other. Also, I'm not a fan of 3-- I like that fields just have values and a space.

[simonbyrne]

How is option 2 functionally different from 3? What else would the distinction between the two spaces be used for other than constructing fields?

They're basically the same, but 3 would require larger interface changes (probably breaking) as they would modify how Fields work. That might be worthwhile though?

[simonbyrne]

I don't follow the nuances, but a comment against 2. is that there's already a tower of types needed to define something that is directly useful to the end user (Domain, Mesh, Topology, Space). In Atmos, I would like to provide something of immediate relevance to the end user (e.g., CubedSphere(radius, nh_poly, ...)), so I would probably end up packaging the FiniteDifferenceGrid in yet another struct.

That's a reasonable point. I was hoping to find an object that would allow us to avoid creating a wrapper object.

[dennisYatunin]

One other option that avoids additional wrapper types/levels in the hierarchy is

abstract type FiniteDifferenceSpace end
mutable struct CenterFiniteDifferenceSpace <: FiniteDifferenceSpace
    # information at cell centers
end
mutable struct FaceFiniteDifferenceSpace <: FiniteDifferenceSpace
    # information at cell faces
end

I just found out while working on #1445 that column_integral_definite! is incorrect for input fields on cell faces because the field returned by Fields.Δz_field does not include information about half cell widths at boundary faces. If we want to fix this (and other related issues that might crop up), we might need to store different metric terms in a CenterFiniteDifferenceSpace and a FaceFiniteDifferenceSpace, in which case there won't be much downside to making them distinct structs.

[simonbyrne]

@dennisYatunin the challenge is how to link the two, which we need for operators so what we know that the gradient of a field on a CenterFiniteDifferenceSpace returns a field on the corresponding FaceFiniteDifferenceSpace. The obvious answer would be each one keeps a reference to the other, but that's difficult to do as Julia doesn't allow mututally recursive types (JuliaLang/julia#269).