Boussinesq model vs incompressible model

[tomchor]

Should we rename IncompressibleModel to BoussinesqModel or something along those lines? Those two are sometimes used interchangeably but really what we have is a Boussinesq model, no?

[ali-ramadhan]

If you use the current IncompressibleModel without any tracers so the Boussinesq approximation doesn't come into play then isn't it truly an incompressible model?

I'm not necessarily opposed to the change but one potential benefit of keeping the IncompressibleModel name is that it might be more familiar to most new users, some who may not be using Oceananigans.jl to model geophysical flows.

[tomchor]

If you use the current IncompressibleModel without any tracers so the Boussinesq approximation doesn't come into play then isn't it truly an incompressible model?

It is, I agree with that. But then again, if you include a buoyancy model (which is the general case) then it's a Boussinesq model.

Also the first google result for "boussinesq model" gives the inexperienced user a wikipedia article that's quite clear (imo) in explaining explaining the model, so I think confusing new users wouldn't be much of a problem.

[ali-ramadhan]

I agree that BoussinesqModel is the more general case.

I also agree that Googling for "Boussinesq model" returns the relevant result but I wouldn't say it's super clear, especially to a newcomer (perhaps an undergrad), how it might relate to what they want to simulate. To someone wanting to simulate a simple CFD setup like a lid-driven cavity, having to read about approximations for buoyancy-driven flows on Wikipedia might obfuscate the fact that they indeed want to use the BoussinesqModel even though they're not using any buoyancy model.

We could use both names if IncompressibleModel becomes an alias for BoussinesqModel but with tracers = nothing, buoyancy = nothing, etc. Not sure we want to use aliases though as we've been on a trend of using fewer aliases (e.g. #613).

As a side note: does it make sense to also consider renaming HydrostaticFreeSurfaceModel to something less wordy at the same time as we consider renaming IncompressibleModel?

[tomchor]

I think we are never going to be both 100% accurate and 100% clear for everyone at the same time, so I guess we need to pick our poison. In this matter my personal opinion is that we should choose the more accurate name, but I'm not 100% sure either.

And I also think this is a pretty big breaking change (given that atm most Oceananigans scripts probably use IncompressibleModel), so maybe it would be good to get input from other people as well?

CC: @glwagner @francispoulin @navidcy

[francispoulin]

I would argue that all the models in Oceananigans to date are Boussinesq since they make the Boussinesq approximation. Also, all models are incompressible in 3D so that is not a great name.

IncompressibleModel is really a non-hydrostatic model with a rigid-lid. If we wanted to be precise, and follow the naming convention of HydrostaticLinearFreeSurface, then I would suggest NonhydrostaticRigidLId. That doesn't roll of the tongue though now does it?

I think this is a good discussion to have and agree that lots of input is important as any change could be a big one and good to make it one that we think will last, a long time.

[tomchor]

I actually like this nomenclature. It's much more correct and descriptive. I'd slighly change it to something like HydrostaticLinearFreeSurfaceModel and NonhydrostaticRigidLidModel.

[francispoulin]

Thanks @tomchor for both correcting my typo and adding in Model. I prefer yours as well.

[glwagner]

Just to frame the discussion: this is probably most important in the context of new users or scientists trying to read examples, docs, and scripts that implement numerical experiments for scientific purposes.

In that context I agree with @ali-ramadhan that BoussinesqModel is misleading for models with tracers = nothing and buoyancy = nothing. I think we've also discussed making this default, since its simple...

As for the distinction between hydrostatic and non-hydrostatic models I agree that with currently supported configurations, NonhydrostaticRigidLid is better than HydrostaticFreeSurfaceModel.

I'm only hesitant about this change because I wonder if we may want to combine these two models in the future. The motivation to develop HydrostaticFreeSurfaceModel separately is largely practical and caves to short term timeline pressures (it's easier to break things when nobody's using it, and we throw a warning right in the constructor so people know its experimental).

In the long run it might make sense to have just a single interface to all incompressible models and support all combinations of hydrostatic, non-hydrostatic, free surface / rigid lid, etc through a single interface, similar to what MITgcm offers. The interfaces / model specific code is expensive to maintain, which motivates having just one interface. If this is our goal, it'd be premature to change the name of IncompressibleModel.

If people feel they'd like to commit to maintaining two interfaces / model constructors / time-stepping code then I think changing the name to be clearer but less general is an option to consider. I personally would rather maintain as few AbstractModel as possible because I feel they are fairly expensive in terms of person-time to maintain. And we probably have a lot of improvements that need to be made (features like model.auxiliary_fields, VectorField for velocity field rather than baking in the coordinate system which fails on the cubed sphere), proper support for systems of coupled nonlinear boundary conditions, etc). It's a good discussion to have regardless.

[francispoulin]

Thanks @glwagner for sharing perspective, that is greatly helpful.

I agree that maintaining different models is a lot of work (as you know better than any one else) and I really like the idea of converging to one model, which has different flags you flip, if you want non-hydrostatic vs hydrostatic, rigid lid vs linear free-surface, and the like. This is a very exciting goal and look forward to working towards making that happen.

In the mean time, while we have different models, it does not matter to me too much what they are called. Changing the names now, which would take a fair bit of work, just to merge the models later on, does not seem necessary to me.

I am happy to keep things as they are and instead spend time thinking about how to generealize the models to cover the different scenarios.

[tomchor]

Yeah I agree with both of you. Apparently unification of models seems advantageous.

But I guess the natural question is: how close does the unification need to be to make sense for us to wait until renaming? I'd say that if the ETA for the unification of models is around 6 months or more, it might be still good to change their names if we can think of more accurate ones.

I have no idea how far we are from combining both models, though!

[glwagner]

There's been some developments since we discussed this last. I think these developments argue for two separately maintained models for the near future and possibly forever. Mostly this regards #1679 and the implementation of general vertical coordinates. While its true that we could implement nonhydrostatic physics in a sigma-coordinate model, I don't personally plan on attempting this any time soon. And I think such an attempt would be a major effort.

It's also just turning out to be not all that difficult to maintain two model constructors. The two models share tons of code under the hood.

So I propose maintaining a split between hydrostatic and nonhydrostatic models for the near future. With that I think it does make sense to rename IncompressibleModel to NonhydrostaticModel. IncompressibleModel is downright confusing given that both the HydrostaticFreeSurfaceModel and, technically, the ShallowWaterModel are all incompressible.

At first I was hesitant to propose this because a blind search and replace might mutilate formatting everywhere in the code. But then I realized that Incompressible and Nonhydrostatic have the same number of characters:

IncompressibleModel
NonhydrostaticModel

wow! :-O So we can actually just blindly search-replace to eternal glory.

Huge API change of course, but there are other huge API changes in our future so I don't think there's any reason to hold back now.

[navidcy]

Unbelievable that characters match. Let’s do it.

[tomchor]

I still like @francispoulin's suggestion of doing NonhydrostaticRigidLidModel, which better matches the name pattern of HydrostaticFreeSurfaceModel. That said, I do agree that it's way easier if the number of characters match and I guess having a rigid lid can be an "implicit characteristic". And I also do think that this is a positive change regardless.

This is a big API change though, and it would be good to hear from other people before merging #1870.

(Also, it feels like we're nearing v1.0?)

[glwagner]

I think it's ok that "rigid lid" is implicit. In fact all the sides are rigid and there's no special direction at all.

In fact, it might make sense to omit "free surface" from the hydrostatic model description. We could soon have a rigid lid option there, and then it wouldn't really make sense to continue calling it a "HydrostaticFreeSurfaceModel".

There's a few things I think should go in before 1.0:

• Abstraction for vectors. It's only for simple grids that we can really get away with referring to the velocity field component wise with u, v, w. Building an abstraction for vectors changes the API because we'll specify boundary conditions on the whole velocity field, not just each component separately. I don't think we should do 1.0 without this.
• Resolve Should we store architecture in grid? #1825 (specifying architecture when building grid). Another major API change, not to mention a substantial internal refactor.
• Simplify grids. We really only need one RectilinearGrid, one LatitudeLongitudeGrid, and one "arbitrary" (not aligned with a coordinate system, like what's used for the cubed sphere). Another API change because we don't need RegularRectilinearGrid and VerticallyStretchedRectilinearGrid. Also simplifies the code a lot.

It could also make sense to build some experience with immersed boundaries in case there are change motivated by that feature.

[tomchor]

In fact all the sides are rigid and there's no special direction at all.

I was basing myself more on the fact that the upper boundary is especial when simulating the ocean from a physics point of view (therefore might requires some extra explanation), even though from the numerical point of view there doesn't need to be anything special about it.

That said, I agree with you that, based off of context, leaving "rigid lid" as implicit is okay. I definitely think that the change as it is proposed is positive.

[glwagner]

In fact all the sides are rigid and there's no special direction at all.

I was basing myself more on the fact that the upper boundary is especial when simulating the ocean from a physics point of view (therefore might requires some extra explanation), even though from the numerical point of view there doesn't need to be anything special about it.

That said, I agree with you that, based off of context, leaving "rigid lid" as implicit is okay. I definitely think that the change as it is proposed is positive.

It's a fair point that oceananographers are "primed" to look for how the free surface is modeled (and the package is Oceananigans after all).

I reflected a bit on this and I think one tricky thing is that I think it may actually be possible in principle (though perhaps tricky or only possible in some crude way) for users to implement an explicitly-time-stepped free surface in the nonhydrostatic model. A "free surface" is just a boundary condition I suppose, and we support fairly general boundary conditions. The free surface implementation in the hydrostatic model isn't much different than this (tracers and momentum actually leave the domain through the top, because the coordinates currently do not follow the surface motion).

[navidcy]

This is a big API change though, and it would be good to hear from other people before merging #1870.

Of course, we shouldn't merge until more voices are heard.

However, in my opinion it's not that big of a change. It's just a rename. Once can even add the alias IncompressibleModel = NonhydrostaticModel... E.g., #1843 was a much bigger change since it required changing the way we construct boundary conditions.

[navidcy]

I think it's ok that "rigid lid" is implicit. In fact all the sides are rigid and there's no special direction at all.

In fact, it might make sense to omit "free surface" from the hydrostatic model description. We could soon have a rigid lid option there, and then it wouldn't really make sense to continue calling it a "HydrostaticFreeSurfaceModel".

There's a few things I think should go in before 1.0:

• Abstraction for vectors. It's only for simple grids that we can really get away with referring to the velocity field component wise with u, v, w. Building an abstraction for vectors changes the API because we'll specify boundary conditions on the whole velocity field, not just each component separately. I don't think we should do 1.0 without this.
• Resolve Should we store architecture in grid? #1825 (specifying architecture when building grid). Another major API change, not to mention a substantial internal refactor.
• Simplify grids. We really only need one RectilinearGrid, one LatitudeLongitudeGrid, and one "arbitrary" (not aligned with a coordinate system, like what's used for the cubed sphere). Another API change because we don't need RegularRectilinearGrid and VerticallyStretchedRectilinearGrid. Also simplifies the code a lot.

It could also make sense to build some experience with immersed boundaries in case there are change motivated by that feature.

in my opinion, for v1.0 we should finalize spherical implementations and bathymetry also...

[navidcy]

I merged the name change. (After all it’s just a name.) I think it’s for the best.

[francispoulin]

I think this is all great! I do like the idea of having 'NonhydrostaticModelandHydrostaticModel`. Whether there's a free-surface or a rigid lid is an option, that should be specified. It is nice to build models that can do either a rigid lid or a free-surface.

[glwagner]

Rigid lid is actually just a short step from having an implicit free surface (and looks like we are about to have a direct solver for simple domains ala #1869...) --- the math is identical for the implicit solve, with a few terms zero'd out (eg, take the limit of an infinitely long time step).

This implies too that models with an implicit free surface that take very long time steps compared to the gravity wave time-scale (most large scale models...) are effectively quite close mathematically to the rigid lid case (notwithstanding effects of runoff, precipitation, and evaporation).

It seems an explicit rigid lid formulation is not very practically useful, mostly because iterative solves take longer so rigid lid models in complex domains where there's no direct solve end up being slower. I'm not sure if there are any advantages re: tracer conservation. The story might get more interesting with generalized vertical coordinates, too.

If anyone is keen to run problems with rigid lids, we can implement it and test it without too much trouble.

MITgcm docs explain this well: https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#implicit-time-stepping-backward-method

[francispoulin]

I know that the rigid-lid and linear free-surface models are very closely related as the grid doesn't change, and I would think that we could use the same solver on both cases. Is that true?

I guess one important difference is in the latter we have to evole the free-surrface but that is pretty cheap compared to everything else as it's only a two-dimensional field.

Dealing with complex geometries is, well, complex, and I don't pretend to understand the nuances there, but in simple rectilinear geometries, I think it would be fairly easy to have an option to go between rigid-lid and linear free-surface in both NonhydrostaticModel and HydrostaticModel. Is this true? Would people want these options?

[glwagner]

I know that the rigid-lid and linear free-surface models are very closely related as the grid doesn't change, and I would think that we could use the same solver on both cases. Is that true?

I guess one important difference is in the latter we have to evole the free-surrface but that is pretty cheap compared to everything else as it's only a two-dimensional field.

Dealing with complex geometries is, well, complex, and I don't pretend to understand the nuances there, but in simple rectilinear geometries, I think it would be fairly easy to have an option to go between rigid-lid and linear free-surface in both NonhydrostaticModel and HydrostaticModel. Is this true? Would people want these options?

Yes, any solver for an implicit free surface can be generalized with some small code modifications to the rigid lid case.

When using a direct solve, there is negligible or actually zero difference in cost between an implicit free surface and rigid lid. When using an iterative solve, the rigid lid can be more expensive because the solution may not converge as quickly.

I don't think a linear + porous free surface with a fixed grid is very useful in the non-hydrostatic model (in some respects the resulting model might be a worse approximation to any real ocean scenario). It can be done, but I don't think there is much benefit from a modeling perspective.

With a breathing / moving grid that follows the surface the story might be a bit different. Having such an abstraction would enable some neat problems that are difficult to solve with any code. But for that someone would have to refactor the non hydrostatic model to use generalized vertical coordinates, and also develop a new pressure solver. This paper by Sullivan and McWilliams describe such a model for atmospheric simulations with a moving, weakly distorted lower boundary:

https://journals.ametsoc.org/view/journals/atsc/71/11/jas-d-14-0095.1.xml

In that paper the "undistorted" pressure solve is used a preconditioner for a fixed point iteration. We could do the same thing except with a conjugate gradient iteration (possibly slightly better than a fixed point iteration). It's possibly a nice project for someone interested in LES above or below weakly distorted boundaries. Note that this method probably won't work when the bathymetry is "too big", because the iterative solve may converge too slowly. cc @whitleyv .