Grid Topology concept

[glwagner]

I think the concept of grid Topology will be useful.

The Topology codifies information about the domain in which the grid is defined. We have three (currently implicit) topologies for each orthogonal direction x, y, z: Flat, Bounded, and Periodic.

A Flat topology has N=1. A Bounded topology has boundaries (either closed or open, though we only supported closed boundaries right now) on left and right. A Periodic topology is periodic in the specified direction.

I propose that we create 3 new type parameters of the grid that codify whether the grid is Flat, Bounded, or Periodic in each of x, y, z.

We can add a keyword to the constructors for grids that allow users to specify the grid topology. For example, the topology that we call a 'Channel Model' would be specified via

grid = RegularCartesianGrid(size=(64, 64, 64), length=(1, 1, 1), topology=(Periodic, Bounded, Bounded)).

A two-dimensional grid can also be defined.

Topologies imply default boundary conditions: a Bounded direction has no-flux boundary conditions; a Periodic direction has Periodic boundary conditions. For a Flat direction we can also use no flux conditions, though it should be irrelevant.

This change will eliminate the need for a special ChannelModel constructor, and also eliminates the need for constructors with names like HorizontallyPeriodicBCs. Instead, we have a generic constructor for FieldBoundaryConditions that lets a user specify boundary conditions in any direction. User-supplied boundary conditions are then checked for consistency with the specified topology, and missing default boundary conditions are added where none exist.

A Flat topology will enable us to set halo sizes to 0 for Flat directions, and elide the unnecessary interpolation and differentiation operations in Flat directions that currently plagues two and one-dimensional models.

A related issue is #330, but this proposal has wider-ranging consequences so I decided to open a new issue.

[glwagner]

@ali-ramadhan do you have any thoughts on the design of x, y, z Topology traits for the grid? I think this is becoming more important, especially because we now have two grids. This may be an essential abstraction, so we may want to make it a "medium" priority to implement it (and also eliminate a lot of boiler plate associated with ChannelModel constructor. It also will break the user API because it will enable a much simpler and straightforward interface for specifying boundary conditions.

[ali-ramadhan]

My biggest concern is that introducing a grid topology and keeping the current implementation of boundary conditions will cause confusion as they perform overlapping functions. But I think with a bit of thinking there can be lots of benefits to adding a grid topology as you say.

This change will eliminate the need for a special ChannelModel constructor, and also eliminates the need for constructors with names like HorizontallyPeriodicBCs. Instead, we have a generic constructor for FieldBoundaryConditions that lets a user specify boundary conditions in any direction.

That would be great I think. If each field now carries around its own boundary conditions and horizontally periodic and channel BCs are treated a grid topology then we can get rid of HorizontallyPeriodicBCs, HorizontallyPeriodicSolutionBCs, ChannelBCs, and ChannelSolutionBCs (+ same for triply periodic and boxes) which would greatly simplify model setup and eliminate most overlap between boundary conditions and grid topology.

A Flat topology will enable us to set halo sizes to 0 for Flat directions, and elide the unnecessary interpolation and differentiation operations in Flat directions that currently plagues two and one-dimensional models.

True. We currently don't have a clean way of eliding unnecessary operations for 1D and 2D models. Dispatching on the grid topology would provide a way of doing this very cleanly.

[ali-ramadhan]

Following up from our discussion:

We can just add three parametric types to AbstractGrid{FT, TX, TY, TZ} to specify three topologies, one for each dimension. Choices are Periodic, Bounded, and Uniform. Users will have to specify the grid topology, e.g.

grid = RegularCartesianGrid(size=(256, 256, 512), length=(500, 500, 100), topology=(Periodic, Periodic, Bounded)

This will be a breaking change. Initially the grid topology won't actually be used, that will come in later PRs.

The grid topology can then be used to construct boundary conditions. We can get rid of model boundary conditions and solution boundary conditions. Here's the API I'm thinking of:

Passing (the default)

boundary_conditions = nothing  # Not sure what else to put?

to the Model constructor will assign the "default boundary conditions" to all fields, e.g. a grid with {Periodic, Periodic, Bounded} will, barring exceptions like pressure maybe, assign horizontally periodic field boundary conditions to all fields.

To assign explicit boundary conditions:

u_bcs = FieldBoundaryConditions(top = BoundaryCondition(Flux, 1e-5))
T_bcs = FieldBoundaryConditions(top    = BoundaryCondition(Flux, -5e-4),
                                bottom = BoundaryCondition(Gradient, 0.01))
...
model = Model(
...
    boundary_conditions = (u=u_bcs, T=T_bcs)
)

will assign specific boundary conditions on u and T while all other fields get the default. This should also work with e.g. (ρu=U_bcs, ρs=S_bcs) for other models.

So users will only have to specify a named tuple of FieldBoundaryConditions and the Model constructor figures it out.

FieldBoundaryConditions will now accept all six kwargs north, south, east, etc. and an error will be thrown if you try to explicitly assign a boundary condition along a Periodic or Uniform dimension.

Thinking a bit ahead: If we want boundary conditions to be a property of each field, we need to create the boundary conditions first then create the fields. Currently fields are instantiated in the Model constructor though.

I think we can get around this (and make #416 easier to resolve) by adding something like a fields kwarg that replaces the velocities and diffusivities kwarg. Passing something like fields = (u=u_array, v=v_array, ...) would initialize the fields with specific values from u_array, v_array, etc.

X-Ref: #606

[glwagner]

This will be a breaking change.

It will break the boundary condition constructors, but need not break the grid constructor if we use a doubly-periodic default.

Passing (the default) to the Model constructor will assign the "default boundary conditions" to all fields

I think we want to write a constructor to create default boundary conditions. This can be the default in the Model constructor.

FieldBoundaryConditions will now accept all six kwargs north, south, east, etc. and an error will be thrown if you try to explicitly assign a boundary condition along a Periodic or Uniform dimension.

It does seem that we unfortunately need to manually check these cases. Ideas welcome for alternative solutions...

To assign explicit boundary conditions:

u_bcs = FieldBoundaryConditions(top = BoundaryCondition(Flux, 1e-5))
T_bcs = FieldBoundaryConditions(top    = BoundaryCondition(Flux, -5e-4),
                                bottom = BoundaryCondition(Gradient, 0.01))

Your FieldBoundaryConditions constructor should take a positional argument grid.

Thinking a bit ahead: If we want boundary conditions to be a property of each field, we need to create the boundary conditions first then create the fields. Currently fields are instantiated in the Model constructor though.

Rather than requiring users to constructor field arrays, which I think will introduce a lot of boilerplate to scripts and make the code harder to use, we can move the construction of fields from the constructor function signature into the constructor function body. This way we can ensure that boundary conditions are set after the fields are created in the case that users do not create the fields themselves, and also allow users to directly pass in the velocity fields with pre-set boundary conditions to velocities, if they desire.