CliMA / ClimateMachine.jl

Climate Machine: an Earth System Model that automatically learns from data
https://clima.github.io/ClimateMachine.jl/latest/
Other
453 stars 78 forks source link

Teaching RectangularDomain, SpectralElementField, and RectangularElement to walk #1736

Open glwagner opened 4 years ago

glwagner commented 4 years ago

From a discussion with @simonbyrne...

PR #1716 and #1723 introduce CartesianField, CartesianDomain, and RectangularElement.

Rather than CartesianDomain, a better name might be RectangularDomain, since the domain really specifies a shape.

There's no need for a "CartesianField" --- instead we can generalize fully to DiscontinuousField, which stores an array of elements, domain, and the map between its elements and the underlying MPIStateArray. The DiscontinuousField constructor dispatches on domain, and the structure of elements depends on the type of domain (eg, domain::RectangularDomain implies that i, j, k correspond to x, y, z).

We can avoid specifying boundary tags in the constructor for RectangularDomain if we divorce grid from domain and build grid inside HydrostaticBoussinesqSuperModel. domain can still be the Coordinating Object where, for example, the floating point type of the experiment should be specified. This requires keeping both domain and grid inside DiscontinuousField (not just domain). This will allow us to design a utility for boundary_conditions for all the models that's separate from domain specification.

(Ultimately, it'd be great if, somehow, the element abstraction were restricted to true discontinuous spectral elements, and that assemble returns a new AssembledField object rather than just a big Element. Getting the abstractions and patterns right for this change will be a bit tricker so we may want to open a new issue for this once the low hanging fruit is plucked.)

glwagner commented 4 years ago

A few notes:

Separating the domain and grid might be a good idea, but could create some convolution largely because we already have a grid, namely the DiscontinuousSpectralElementGrid. In order to go from the unstructured grid to a Cartesian or Spherical representation, we need (want?) to know how many elements there are in each direction.

I think one solution might be to add a domain property to DiscontinuousSpectralElementGrid.

Then we have the information we need to deduce order from the chaos, aka, to extract a Cartesian view into the unstructured data in DiscontinuousSpectralElementGrid and MPIStateArray. We can add an "empty" domain type with the bland name UnstructuredDomain to use as a default within the inner sanctum Mesh.Grids

We can also write constructors for DiscontinuousSpectralElementGrid that dispatch on domain and build the appropriate topologies and whatnot:

function DiscontinuousSpectralElementGrid(domain::Rectangle; Ne, Np)
    # stacked bricks, MPI, etc
end

domain = Rectangle(x = (-1, 1), y = (-1, 1), z = (0, 1), periodic=(true, false, false))
grid = DiscontinuousSpectralElementGrid(domain, Ne = (24, 24, 24), Np = 4)

Important information about the domain (extents, periodicity, origin) is stored in grid.domain.

The only missing piece is how to retrieve the tuple of elements (24, 24, 24) from grid, which we need for Cartesian and Spherical views into MPIStateArray. It seems like this might be inferable from properties of the grid like topology, so if we can't provide a property grid.Ne at least we can provide a function elements(grid) that returns either something like size(grid.vgeo, 3) for unstructured grids, or a 3-tuple for structured grids. Perhaps @simonbyrne or @jkozdon can advise?

Happy to work on this after #1723 goes through.