Oceananigans.jl: Do TimeSteppers belong to Models, or to Simulations?

[glwagner]

Currently, IncompressibleModel keeps a reference to TimeStepper:

https://github.com/CliMA/Oceananigans.jl/blob/03a6f855f839504d94cb8cee3c2665b17afbc6d5/src/Models/incompressible_model.jl#L33

This means that every model has one time stepper, and we use time_step!(model, dt) to advance a model forward.

Yet this is not a necessary entanglement. A perfectly valid time-stepping loop is

for i = 1:100
    time_step!(model, forward_euler_stepper, dt)
    time_step!(model, rk3_stepper, dt)
end

(obviously we don't have a forward Euler time-stepper, which is a separate issue, but hopefully the point is made.)

We can disentangle the timestepper from model by moving the timestepper reference over to Simulation. This would also clean up the interface since we wouldn't need to use symbols to specify the time-stepping method (#1119, this might be a hint that we're on the right track...)

model = IncompressibleModel(...)
simulation = Simulation(model, RungeKutta3TimeStepper(model), kwargs...)

I think there's a conceptual advantage to this too. Simulations manage the creation of a time-series, while a Model is more the idealization of a discrete physical system at a particular moment in time.

To make this change, we first have to checkpoint Simulations rather than Models to support time-steppers that require history like AB2. I think this is something we need anyways.

This is a major API change, though it could be mitigated if we put a default into the Simulation constructor.

[navidcy]

I agree with this point.

But, to play devil's advocate one could argue that the advection scheme should not be part of the model either, no? One could intuitively think that model=PDE and simulation=PDE+numerical schemes required to solve it.

[glwagner]

I think it's useful to think of the "model" as a spatially discrete approximation to the PDE. The simulation sends it barreling forward through time.

[navidcy]

Framed this way, it makes absolute sense to include time stepper in simulation.

[francispoulin]

I am happy to follow this discussion and even though I am very new to the table but will share a perspective.

Outside of the context of Oceananigans, I would think of the collection of PDEs as the physical model, and the method by which we discretize them in both space and time, as the numerical model. In this context, when we say model we mean the spatial discretization of the physical model, which is most of but not quite all of the numerical model. The fact that we don't have Advection in the Models folder is consistent and a good idea.

It makes sense to me to put the time stepping in simulations since that is the part where we turn the crank and find the solutions to the prognostic variables.

I suspect the user might benefit from a little bit of discusson on this in the docs so that they can better understand the framework you've set up, which works great. This is of course something to be done when you have spare time, if such a thing really exists. ;)

[glwagner]

Well, right now, the Model includes both discrete and continuous aspects of the PDE. Model also stores architecture, for example, which is unrelated to the PDE.

We did consider an abstraction that represented the equation set independent from discretization. This would not be a replacement for model, but would simply be a way to organize some of the properties of model (eg buoyancy, coriolis, diffusion, etc). Yet this is a challenging abstraction to design. One issue is specification. If we try to separate the terms from their numerical implementation, we have to figure how to distinguish between "diffusion with second-order differences" and "diffusion with fourth-order differences". It's clearly possible to do this, but it's going to take some time and careful thought to implement. We've taken a more incremental approach to development the model instead, resulting in a fairly "flat" interface to IncompressibleModel that combines aspects both of the continuous equations, their discretization, and things like architecture, the pressure_solver, etc.

I think flat is simple (think parameter files as the platonic ideal of a flat API); however, we still want to benefit from modularity where it's a simple change that doesn't overcomplicate the interface (hence this issue).

[francispoulin]

That makes sense and I am a firm believer in starting simple, and I'm glad you made that choice.

[navidcy]

"Platonic ideal" :)