glwagner commented 3 years ago

As discussed on #1153 we may start developing a ShallowWaterModel soon.

To reuse as much code as possible, we may want to refactor the TimeSteppers and Models modules. Right now, TimeSteppers is included after Models:

https://github.com/CliMA/Oceananigans.jl/blob/4b972b580e7db9f49c314ab09280a2f99e5c1e32/src/Oceananigans.jl#L136-L137

However, with more than one model, we probably want to use dispatch to control the behavior of key functions in the time-steppers. For example, we'd have

calculate_pressure_correction!(::ShallowWaterModel) = nothing
pressure_correct_velocities!(::ShallowWaterModel) = nothing

Thus I think what we want to do is to include the TimeSteppers module first, with the definitions

function calculate_pressure_correction! end
function pressure_correct_velocities! end

and subsequently add methods in the Models submodule appropriate to the physics / formulation of each model.

glwagner commented 3 years ago

I guess we will also need to overload rk3_substep!, update_state!, and calculate_tendencies! as well for ShallowWaterModel. We should be able to reuse almost everything else I think, including many of the physics / tendency implementations.

Following this pattern hopefully will mean that if we implement new time-steppers in the future we can reuse them for both models. I could see a QuasigeostrophicModel fitting into this framework as well.

glwagner commented 3 years ago

cc @ali-ramadhan @navidcy @francispoulin

francispoulin commented 3 years ago

This is going to be a silly question but I have to ask. I have been browsing the code and getting to know where things are but I'm having difficulties doing something that I thought would be obvious. Where are the modle equations specified?

I looked in Models and that seems to have a lot but don't actually see the governing equations, or what I recognize as a governing equation.

I think my first step is to find these equations and modify them for shallow water. I see that the Incompressible model is defined in terms of a basic state, and I think that would be easy enough to do for this case, and useful for somethings.

glwagner commented 3 years ago

The tendency terms are calculated in this file:

https://github.com/CliMA/Oceananigans.jl/blob/master/src/TimeSteppers/velocity_and_tracer_tendencies.jl

for example:

https://github.com/CliMA/Oceananigans.jl/blob/bd6cb932cc7d6c372709d4b5af59d325b494b30c/src/TimeSteppers/velocity_and_tracer_tendencies.jl#L39-L61

The functions that appears in the tendency for a velocity field are distributed throughout the code. For example, advection terms are defined in Advection, coriolis terms are defined in Coriolis, and turbulence closure terms are defined in TurbulenceClosures.

I believe what we want is to provide an independent implementation of the function calculate_tendencies! for ShallowWaterModel. Right now the function calculate_tendencies! is defined in TimeSteppers:

https://github.com/CliMA/Oceananigans.jl/blob/master/src/TimeSteppers/calculate_tendencies.jl

There's some annoying software stuff (eg the items raised in this issue) that we'd have to solve before being able to cleanly implement a ShallowWaterModel.

francispoulin commented 3 years ago

Ah, sorry that I missed that. I see that the four tendances (what I would call fluxes) are defined there as inline functions. Very nice and readable!

I agree that we need to define different calculate_tendencies. Two thoughts come to mind.

Can define a folder TimeSteppers_SW in the same folder as TimeSteppers and then define the corresponding tendendcies and the like for SW. This seems doable but maybe not incredibly efficient as we would need to redefine AB2 and RK3.
Define ShallowWater_tendencies.jl and calculate_tendencies_ShallowWater.jl (or SW) in the src folder. That way we can use the time steppers but just have to make sure they can deal wtih both types of tendencies. I think this is doable, and maybe slightly more desirable to me. What do you think?

glwagner commented 3 years ago

I'm proposing in this isseu that we reorganize the code a bit, moving calculate_tendencies! out of the TimeSteppers module, and putting it into Models.

The file structure would be

TimeSteppers/
    runge_kutta_3.jl
    ...

Models/
    Models.jl

    IncompressibleModels/
        IncompressibleModels.jl
        calculate_incompressible_tendencies.jl
        other_incompressible_specific_files.jl

    ShallowWaterModels/
        ShallowWaterModels.jl
        calculate_shallow_water_tendencies.jl
        other_shallow_water_specific_files.jl

The term calculate_incompressible_tendencies.jl would look like:

import Oceananigans.TimeSteppers: calculate_tendencies!

function calculate_tendencies!(model::IncompressibleModel)
    # the code currently in calculate_tendencies.jl
    return nothing
end

The file calculate_shallow_water_tendencies.jl would contain new functionality

import Oceananigans.TimeSteppers: calculate_tendencies!

function calculate_tendencies!(model::ShallowWaterModel)
    # New functionality
    return nothing
end

This implementation will cause minimal disruption to the code, is consistent with the style we're currently using, and will generalize to additional models other than IncompressibleModel and ShallowWaterModel.

What we call tendencies are terms like du/dt, dv/dt, etc. We typically refer to "fluxes" as terms like wu (vertical momentum flux) and wc (tracer flux). The flux divergences contribute to tendencies.

francispoulin commented 3 years ago

That sounds like a very solid plan and I like it.

Just to be clear, the folders should have the following:

TimeSteppers/
    TimeSteppers.jl
    quasi_adams_bashforth_2.jl
    runge_kutta_3.jl

Models/
    Models.jl

    IncompressibleModels/
        incompressiblemodels.jl
        calculate_incompressible_tendencies.jl
        pressure_correction.jl
        store_tendencies.jl
        update_hydrostatic_pressure.jl
        update_state.jl
        velocity_and_tracer_tendencies.jl

    ShallowWaterModels/
        ShallowWaterModels.jl
        calculate_shallow_water_tendencies.jl
        other_shallow_water_specific_files.jl

Even before SW is created, this requires moving things around, as you have shown above. Did you want me to give this a try? It does not seem like a lots needs to change and I will certainly learn a lot in the process.

On a different note, I've noticed that when I view the docs in chrome (on ubuntu) the math equations aren't actually converted. I have a chrome extension but that asks for $ and $$ and so forth. The docs are using `` instead. What should one use to view this? I went down a rabbit hole of finding markdown viewers and sadly did not succeed in figuring this out.

ali-ramadhan commented 3 years ago

On a different note, I've noticed that when I view the docs in chrome (on ubuntu) the math equations aren't actually converted. I have a chrome extension but that asks for $ and $$ and so forth. The docs are using `` instead. What should one use to view this? I went down a rabbit hole of finding markdown viewers and sadly did not succeed in figuring this out.

Just to double check, are you viewing the documentation source code (e.g. https://github.com/CliMA/Oceananigans.jl/tree/master/docs/src) or the built/compiled documentation (e.g. https://clima.github.io/OceananigansDocumentation/stable/generated/eady_turbulence/)?

You shouldn't see any Markdown artifacts like $ and `` if you're viewing the built documentation, but maybe a few equations are not converted to LaTeX properly?

I think the Markdown written for Julia documentation (then built with Documenter.jl) used to use $ but has now switched to using `` and

```math

francispoulin commented 3 years ago

Thanks @ali-ramadhan . You are correct. I see the raw latex when looking at the source code but everything looks pretty with the compiled documentation.

I presume you use Documenter.jl to build the docs?

If yes I will learn how that works.

glwagner commented 3 years ago

Yes, we use Documenter.jl, which supports latex (and therefore is awesome)

glwagner commented 3 years ago

Did you want me to give this a try? It does not seem like a lots needs to change and I will certainly learn a lot in the process.

That's fine with me (and I am very happy to help, or pair program over zoom). I think resolving this issue (moving files around) should happen in a separate PR from implementing ShallowWaterModel. Imports and exports will also change, for example IncompressibleModels will become a submodule of Models, and we'll have to pull the name IncompressibleModel into Models:

# in Models.jl

include("IncompressibleModels/IncompressibleModels.jl")

using IncompressibleModels.IncompressibleModel

Hopefully we can reuse as much code as possible. If your prognostic variables are hu,, hv, h, you may benefit from some fancy julia abstraction that will let you reuse momentum advection and tracer advection operators like div_Uu and div_Uc. I'm happy to help with that.

Your fields will all have to be three-dimensional (with a single grid point in the vertical direction). Also until #1024 is fixed we will unfortunately also need halo regions in the vertical direction. Hopefully this won't matter too much and eventually we will be able to use Flat vertical directions.

As for the file structure, I think some functions can be generalized to work for any model rather than split, like store_tendencies! (it perhaps makes sense that store_tendencies! is independent of the model type).

I think it's nice that the physics-specific functions will now be in Models/IncompressibleModels.

francispoulin commented 3 years ago

Thanks for the offer of doing this together. Since I know that it would go much faster doing this together, and also I will learn more from talking to you, how about we do this over zoom? Let me know by email what you prefer.

Yes, I do think that we can use a lot of the existing code, which will make things easier and consistent.

I do have a question of the form of the equations being solved. The docs suggest that the advection terms are written as u \cdot \nabla u for example. Is it solved like this or is incompressability used to write it as \nabla \cdot (u u)?

Since this is a finite volume method I presume it's the latter but not entirely sure.

Having extra ghost cells for a whlle should not be a problem as this is just temporary.

I agree that much of the code can be shared, and should be shared where possible.

francispoulin commented 3 years ago

going back to @ali-ramadhan 's response.

When preparing a markdown document, it is convenient to view it in the pretty form, as one does in latex. Are there any markdown viewers that people might recommend for linux? If not no big deal, I will learn to use Documenter but that sounds like it might be a bit slow as you have to run code to compile it, but maybe I'm wrong.

ali-ramadhan commented 3 years ago

I do have a question of the form of the equations being solved. The docs suggest that the advection terms are written as u \cdot \nabla u for example. Is it solved like this or is incompressability used to write it as \nabla \cdot (u u)?

Yes you're right, we solve/compute it like ∇⋅(uu). For example, for centered second-order advection: https://github.com/CliMA/Oceananigans.jl/blob/master/src/Advection/centered_second_order.jl#L11

For reference, ℑxᶜᵃᵃ is a second-order interpolation operator from cell faces to cell (C)enters (therefore the superscript c): https://github.com/CliMA/Oceananigans.jl/blob/master/src/Operators/interpolation_operators.jl#L7

These could probably be easily reused to advect e.g. hu for shallow water. Here's what I did to advect ∇⋅(ρuu) for a compressible model (basically advect ρuρu then divide by ρ haha): https://github.com/CliMA/Oceananigans.jl/blob/ar/compressible-model/compressible/src/Operators/compressible_operators.jl#L140-L149

ali-ramadhan commented 3 years ago

When preparing a markdown document, it is convenient to view it in the pretty form, as one does in latex. Are there any markdown viewers that people might recommend for linux? If not no big deal, I will learn to use Documenter but that sounds like it might be a bit slow as you have to run code to compile it, but maybe I'm wrong.

Ah yes it's quite painful :(

I guess it's not enough to look at the formatted Markdown as Documenter.jl might have to process it to add in any cross-links, references, citations, results from code blocks, math, etc.

So usually when I work on documentation I have a REPL open to run make.jl so I can look at my changes (first time is slow but subsequent runs are decently fast). But this is still slow and it's annoying because I have to disable building the examples by commenting them out (examples run as part of building the documentation).

There must be a better way... I'll ask on the #documentation channel on the Julia Slack to see if anyone has a better workflow.

glwagner commented 3 years ago

@francispoulin one strategy for building the documentation more quickly is to comment out everything that you don't need to look at from pages. The examples in particular take a while to build since some of them run 3D simulations...

francispoulin commented 3 years ago

Yesterday we discussed solving the SW model in conservative form, which is a perfectly good approach. However, I remembered that there another formulation in terms of the velocity where we use the vorticity and the Bernouli function. This is also very nice and this approach shows there are the different finite difference formulations that ensure either energy or enstrophy conservation. I don't know how Finite Volume would compare, but I'm tempted to start with this for simplicity.

http://www.aos.wisc.edu/~aos718/sadourny.pdf

francispoulin commented 3 years ago

More thoughts.

I'm reading the discussion on immersed boundary methods (IBMs) with interest on #916. It should be pointed out that this would be a great thing to test in a shalllow water model since there is no pressure inversion and it's much easier to resolve the horizontal, since we have no vertical.
Climate Machine also has a shallow water model. I don't know the details of what this does but I should probably learn how the two models will differ.

glwagner commented 3 years ago

2. Climate Machine also has a shallow water model. I don't know the details of what this does but I should probably learn how the two models will differ.

I'm not 100% sure how its implemented but I think at the end it'd be great if both ClimateMachine and Oceananigans shallow water models used the same equation set, since we can then compare the numerics. We don't have the bandwidth to develop the ClimateMachine shallow water model, but possibly at the time the Oceananigans has a nice one we'll be motivated to flesh out ClimateMachine's model.

I'm not sure about immersed boundaries --- does the method generalize in a simple way? I'm not sure whether a special treatment is required for the height field, compared to an ordinary velocity variable, or tracer.

glwagner commented 3 years ago

I don't know how Finite Volume would compare, but I'm tempted to start with this for simplicity.

Is that method implementable on a C-grid?

francispoulin commented 3 years ago

I agree @glwagner that we want an Oceananigans version and happy to help do that. Comparing them would be a lot of fun.

Using the same equations would be nice too. I looked at the julia code but having problems finding the equations, again.

Yes, it is a C-grid. IncompressibleModels also uses the same grid, right?

johncmarshall54 commented 3 years ago

The climate machine shallow water model used DG numerics. I agree with Greg that it would be nice to do a side-by-side comparison of two numerical approaches to the same underlying equations. Not sure it used a Bernoulli formulation. john ps Francis - great that you are involved with oceananigans!

On Fri, Nov 13, 2020 at 11:53 AM Gregory L. Wagner notifications@github.com wrote:

Climate Machine also has a shallow water model https://github.com/CliMA/ClimateMachine.jl/blob/master/src/Ocean/ShallowWater/ShallowWaterModel.jl. I don't know the details of what this does but I should probably learn how the two models will differ.

I'm not 100% sure how its implemented but I think at the end it'd be great if both ClimateMachine and Oceananigans shallow water models used the same equation set, since we can then compare the numerics. We don't have the bandwidth to develop the ClimateMachine shallow water model, but possibly at the time the Oceananigans has a nice one we'll be motivated to flesh out ClimateMachine's model.

I'm not sure about immersed boundaries --- does the method generalize in a simple way? I'm not sure whether a special treatment is required for the height field, compared to an ordinary velocity variable, or tracer.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/CliMA/Oceananigans.jl/issues/1165#issuecomment-726873787, or unsubscribe https://github.com/notifications/unsubscribe-auth/AKXUEQWKUQJEBFEYXSOVQO3SPVQABANCNFSM4TQVNFEA .

--

John Marshall Earth, Atmospheric and Planetary Sciences, MIT http://oceans.mit.edu/JohnMarshall/

francispoulin commented 3 years ago

Thank you @johncmarshall54 for the information. I will look through the code and figure it out.

I am happy to work on this great project and happy to be back in touch. I hope you are doing well.

glwagner commented 3 years ago

@francispoulin I got confused by this image

But now I see that it is indeed equivalent to a C-grid.

The ClimateMachine implementation uses tools associated with the dynamical core and is thus written in the "balance law framework". It might be worth checking out documentation or tutorials first before looking through the source but I'm not sure the docs will help either...

francispoulin commented 3 years ago

Thanks @glwagner , I did figure out that the equations are referred to as balance laws and was looking in that folder but didn't see anything obvious. I will go to the docs and read through them more carefully to understand what's going on.

francispoulin commented 3 years ago

I looked over the documents and learned much better how the equations, what they call balance laws, are put together. Very nice and clean system.

I do see that in the code the model equations are in conservative form, as you will find in equations 7,8, and 9 of this paper.

Given that's the case I will ignord the Sadourny approach that I mentioned before. Thanks for the discussion to help me figure this out.

I will start wtih a one-dimensional version of the linear equations without rotation (essentially the wave equation) and build it up from there.

glwagner commented 3 years ago

I will start wtih a one-dimensional version of the linear equations without rotation (essentially the wave equation) and build it up from there.

That sounds smart!

johncmarshall54 commented 3 years ago

Good plan

On Fri, Nov 13, 2020, 2:27 PM Francis J. Poulin notifications@github.com wrote:

I looked over the documents and learned much better how the equations, what they call balance laws, are put together. Very nice and clean system.

I do see that in the code https://github.com/CliMA/ClimateMachine.jl/blob/master/src/Ocean/ShallowWater/ShallowWaterModel.jl the model equations are in conservative form, as you will find in equations 7,8, and 9 of this paper. https://pdfs.semanticscholar.org/9981/51bcc42a04b281811b3b02217799c352b52f.pdf

Given that's the case I will ignord the Sadourny approach that I mentioned before. Thanks for the discussion to help me figure this out.

I will start wtih a one-dimensional version of the linear equations without rotation (essentially the wave equation) and build it up from there.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/CliMA/Oceananigans.jl/issues/1165#issuecomment-726988741, or unsubscribe https://github.com/notifications/unsubscribe-auth/AKXUEQTMZA3ATENPYU2Q5VLSPWCBJANCNFSM4TQVNFEA .

francispoulin commented 3 years ago

I don't suspect this is the best place to post this but I'll start here and happy to move the conversation over to wherever it should be happening.

I have created a new branch called "ShallowWaterS1waveeqn". S1 is for step 1. It doesn't do much right now but I'm trying to set up the skeleton so that we can actually do something intereting, like evolve the equation.

At the moment I have the following:

Created a ShallowWaterModels folder
Created a ShallowWaterModels.jl module. It doesn't do much except include shallowwater_model.jl
shallowwater_model.jl defines a mutable struct and a function. I couldn't get a lot of things to work but it does regularlize the boundary condtions, set up the velocities and a new variable called layer_depth. The name is not great and can go with depth or height, as people prefer. It will be in in the governing equations.
Created an example, in the same folder, called onedim_shallowwater.jl. Currently, this includes the module mentioned above, defines a grid (yes, a single point!), defines the initial conditions and sets the velocity.

Questions:

This is pretty bare bones so far but does this follow the philosophy of Oceananigans?
I have tried to inititalze the layer_depth by adding an argument into the set! command but it failed. Below is the command that I thought would work and below is the error because layer_depth is not model.velocities or model.tracers. Can someone help me to see what silliness I am doing?


set!(model, u = u, v = v, layer_deth = layer_depth)
fails with "ERROR: LoadError: ArgumentError: name layer_depth not found in model.velocities or model.tracers."

glwagner commented 3 years ago

@francispoulin perhaps the right place for this post is a new PR? (We can continue discussing stuff there.)

francispoulin commented 3 years ago

Thanks, will do that right now

CliMA / Oceananigans.jl

Making room for `ShallowWaterModel` #1165

--

John Marshall Earth, Atmospheric and Planetary Sciences, MIT http://oceans.mit.edu/JohnMarshall/