Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Begins the implementation of Oceananigans.HydrostaticFreeSurfaceModel for general circulation modeling #1349

Merged
merged 7 commits into from
Feb 10, 2021

Conversation

glwagner
Copy link
Member

@glwagner glwagner commented Feb 9, 2021

This PR adds a new AbstractModel called HydrostaticFreeSurfaceModel to IncompressibleModel and ShallowWaterModel in the Oceananigans model suite.

HydrostaticFreeSurfaceModel is intended to solve the hydrostatic Boussinesq equations beneath a free surface. This PR adds the basic infrastructure needed to time step models on a RegularCartesianGrid with explicit free surface dynamics.

I think we should save dynamics tests and an example for a future PR, and just work on getting some of the basic ingredients into this model. The time-stepping tests are currently in place, but most of time_step!(model::HydrostaticFreeSurface, dt) needs to be fleshed out.

Todo for this PR:

  • correct tendency calculation and time-stepping for an Adams-Bashforth timestepping method
  • calculation of vertical velocity from the continuity equation

function HydrostaticFreeSurfaceModel(; grid,
architecture::AbstractArchitecture = CPU(),
clock = Clock{eltype(grid)}(0, 0, 1),
advection = (momentum=VectorInvariant(), tracers=CenteredSecondOrder()),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💪

Out of curiousity I thought VectorInvariant changes the form of the equations being solved, but can we think of it as a different advection scheme?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm generalizing the kwarg advection to mean something like "advection term" (both continuous equation aspects and numerical aspects).

We could try to distinguish between continuous aspects (whether we use a vector invariance form, an "interior" form, or a conservation form) and numerical aspects (numerical approximations / stencils used to implement that form).

We could also have a more hierarchical interface; something like

advection = (momentum=VectorInvariantForm(), tracers=ConservationForm(scheme=CenteredSecondOrder()))

I'm open to suggestion; I just did something lazy and simple here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh yeah no I think what you did is great. Just wanted to understand. Makes sense that it encapsulates both the continuous and numerical aspects.

function hydrostatic_free_surface_model_tracers_and_forcings_work(arch)
grid = RegularCartesianGrid(size=(1, 1, 1), extent=(2π, 2π, 2π))
model = HydrostaticFreeSurfaceModel(grid=grid, architecture=arch, tracers=(:T, :S, :c, :d))
set!(model, η=1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah I guess we'll have to do this for a lot of tests, similar to ShallowWaterModel, to avoid blowup when timestepping due to division by zero?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You mean setting the displacement to some constant value?

Since we're modeling displacement rather than the total thickness of the layer, I think we should get away with an initial value of 0. I can remove the call to set! since it's tested in another routine.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah right it's testing set!.

test/test_hydrostatic_free_surface_models.jl Outdated Show resolved Hide resolved
Copy link
Collaborator

@francispoulin francispoulin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have not tried to run it but things look good to me.

One thing that I noticed is that advection seemed to be commented out for tracers. Did I read that right?

Also, time stepping is done differently, and not exactly sure why that is but any time stepping that works is good time stepping!

@glwagner
Copy link
Member Author

One thing that I noticed is that advection seemed to be commented out for tracers. Did I read that right?

Yes --- because we are still finalizing the user interface that controls the advection term + advection scheme for momentum and tracers.

I think ideally we want to preserve the ability to use high-order advection schemes for momentum when we are on a RegularCartesianGrid. However, we also have to use a "vector-invariant form" of the advection term when running on curvilinear grids. So we need a way to specify both continuous and numerical aspects of the advection term implementation.

@glwagner
Copy link
Member Author

@ali-ramadhan @francispoulin I propose that for the time being we use an abstraction where advection has a field momentum and a field for every tracer.

The value of this field controls the implementation of the advection term, with continuous and numerical aspects "flattened" into a single spec. If one writes VectorInvariant, then we implement a second-order VectorInvariant form of the momentum equations. If one provides an advection scheme, when we implement a conservation form of the momentum / tracer equations with the associated stencil.

When we have curvilinear grids we will need a function validate_advection(advection, grid) to ensure that a conservation form scheme isn't used on a curvilinear grid (for example). There's a similar consideration for stretched grids and our current WENO implementation.

@glwagner
Copy link
Member Author

I also wonder if we really need to have two kwargs (momentum_advection and tracer_advection). It's a little messy and complicated to parse user-input if we try to smash everything into one kwarg.

@ali-ramadhan
Copy link
Member

I think a single advection kwarg would make for a cleaner UI, since there are already lots of model kwargs.

I agree it's a bit messier to parse and validate but the parsing and validation can be done by small functions in a separate file.

@francispoulin
Copy link
Collaborator

Hmm, as a user I think there are reasons why you might wnat to have the choice to use different advection schemes for momentum and tracers (and buoyancy). But could there be a default where if there is only momentum advection specified, then tracers are the same?

Copy link
Collaborator

@francispoulin francispoulin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changes look reasonable to me.

It seems that nonlinear advection was added into the tendencies, but were the Coriolis terms also removed?

@glwagner
Copy link
Member Author

It seems that nonlinear advection was added into the tendencies, but were the Coriolis terms also removed?

No they are there.

clock)

return ( - U_dot_∇v(i, j, k, grid, advection, velocities)
- y_f_cross_U(i, j, k, grid, coriolis, velocities)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@francispoulin here's a Coriolis term

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I clearly misread something.

@glwagner glwagner merged commit 0057819 into master Feb 10, 2021
@glwagner glwagner deleted the glw/hydrostatic-free-surface-model branch February 10, 2021 21:37
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants