[WIP] Add ShallowWater module

[navidcy]

This PR adds a ShallowWater module that solves the shallow water equations in two dimensions.

The docs preview should be viewable at http://fourierflows.github.io/GeophysicalFlowsDocumentation/previews/PR163.

Closes #127.

cc @glwagner, @liasiegelman

[navidcy]

Ideally we should finish off FourierFlows/FourierFlows.jl#225 and then use the named tuple functionality for the ShallowWater module.

[glwagner]

Should we write in conservative form?

[navidcy]

Yes.

Also there are some typos, e.g., the -g ∂η/∂y term is missing from the v equation.

Also I'd like to include non-trivial bathymetry (which appears as a forcing term on the right-hand side in the momentum equations.)

[navidcy]

Here they are in conservative form:

I made an attempt to write the calcN! function. Bit tricky...

Also, what's the form of dissipation? We want it to conserve the total hu, hv and h... Right?

[glwagner]

We're having a related discussion at CliMA/Oceananigans.jl#1712

[glwagner]

We can avoid using NamedTuple for sol by expanding problem.vars to include a property, say, "state", that provides a view into the different components of sol.

For example, if sol has 3 components (on an nx, ny = (4, 4)) grid:

julia> sol = rand(3, 4, 3)
3×4×3 Array{Float64,3}:
[:, :, 1] =
 0.611945  0.0972204  0.313608   0.279273
 0.709269  0.484545   0.928009   0.988215
 0.855321  0.597252   0.0225207  0.633623

[:, :, 2] =
 0.724406  0.287053  0.36452   0.594918
 0.473921  0.589545  0.457694  0.61121
 0.880811  0.475327  0.067284  0.100218

[:, :, 3] =
 0.860079  0.102246  0.606423  0.532186
 0.529818  0.142198  0.97109   0.0831012
 0.814437  0.633122  0.352723  0.880207

state is defined

julia> state = (hu = view(sol, :, :, 1), hv = view(sol, :, :, 2), h = view(sol, :, :, 3))
(hu = [0.6119450557692874 0.09722035882999713 0.313607543642898 0.27927324296345324; 0.7092687711628598 0.484545143696397 0.9280087864249529 0.9882153380173344; 0.8553207027735736 0.597252372372197 0.022520721078303607 0.6336230102086127], hv = [0.7244063615451659 0.28705291537063293 0.3645199447746996 0.594917539715488; 0.4739208711881049 0.589544872098053 0.45769444095571066 0.6112098564777788; 0.8808107163282806 0.47532706964313487 0.06728398172249617 0.10021829190676201], h = [0.8600794223299935 0.10224636138594412 0.6064226221047209 0.532185818598262; 0.5298180933763801 0.14219800751413225 0.9710902288597041 0.0831011581631147; 0.8144371452255259 0.6331216842925897 0.35272263622882627 0.880206590278507])

Then we can do operations like

julia> state.hu
3×4 view(::Array{Float64,3}, :, :, 1) with eltype Float64:
 0.611945  0.0972204  0.313608   0.279273
 0.709269  0.484545   0.928009   0.988215
 0.855321  0.597252   0.0225207  0.633623

julia> @. state.hu = 2 * state.hv
3×4 view(::Array{Float64,3}, :, :, 1) with eltype Float64:
 1.44881   0.574106  0.72904   1.18984
 0.947842  1.17909   0.915389  1.22242
 1.76162   0.950654  0.134568  0.200437

julia> state.hu
3×4 view(::Array{Float64,3}, :, :, 1) with eltype Float64:
 1.44881   0.574106  0.72904   1.18984
 0.947842  1.17909   0.915389  1.22242
 1.76162   0.950654  0.134568  0.200437

julia> sol
3×4×3 Array{Float64,3}:
[:, :, 1] =
 1.44881   0.574106  0.72904   1.18984
 0.947842  1.17909   0.915389  1.22242
 1.76162   0.950654  0.134568  0.200437

[:, :, 2] =
 0.724406  0.287053  0.36452   0.594918
 0.473921  0.589545  0.457694  0.61121
 0.880811  0.475327  0.067284  0.100218

[:, :, 3] =
 0.860079  0.102246  0.606423  0.532186
 0.529818  0.142198  0.97109   0.0831012
 0.814437  0.633122  0.352723  0.880207

[aramirezreyes]

Hi! I'm interested in this. Is it actively been worked? Is there a checklist to follow to know if I can contribute to?
Thanks!

[navidcy]

@aramirezreyes that would be so great if you can help out with this!

It's been in my list but fell far down... I've made some progress but then left it. Give me few days and I can draft a checklist of items we need to tick and we can move on from there. :)

[navidcy]

The docs include a first write-up of the shallow water equations and can be viewed locally by calling:

julia --project=docs/ -e 'using Pkg; Pkg.instantiate(); Pkg.develop(PackageSpec(path=pwd()))'; julia --project=docs/ docs/make.jl; open docs/build/index.html

from the local repo's main directory. One question is about notation for the equations in conservative form. The docs use qᵤ = u*h and qᵥ = v*h. Are there any better alternatives to qᵤ, qᵥ? I didn't wanna use U=u*h and V=v*h since they usually these simples imply dimensions of velocities and this is not the case here.

hu and hv could be the variables for h*u and h*v... But then, our convention of using a trailing h to denote the Fourier transform becomes a bit cumbersome, e.g.

huh = rfft(hu)

:)

[navidcy]

Some things we should discuss/figure out.

• Decide on notation.
• Do we code up both conservative and non-conservative form of equations? Is there an advantage for the non-conservative form of the equations?
• Hyper-viscosity terms? Where do they appear? Usually one puts them in the primitive form of the equations (the non-conservative form) in such manner so that integrals of u²+v² or h² are dissipated.

When we settle on the above we should:

• Code up the equations (partially done)
• Add basic diagnostics?
• Write up tests for each term of the equations, bathymetry, diagnostics... (follow structure of other module's tests)
• Add at least one example that'll be literated via the Docs. Potentially two examples... one example could be demonstrating something that people find in textbooks, e.g., a near-linear surface wave propagation.

@glwagner, ideas/feedback?

[glwagner]

I don't really know what the advantages of a conservative vs. primitive formulation are for a spectral method. @kburns might have some thoughts? A major disadvantage of the conservative form is that viscous/hyperviscous terms are nonlinear? Probably the most thorough way to make a decision would be to code up both and test.

I wouldn't worry about diagnostics for this PR. It's nice if PRs are not too ambitious, so they can be merged more easily. Diagnostics can be added by users that need them.

I think it's also ok to write up the difficult physics tests after merging the PR (with unit tests), at least if there's time to keep working on this in the future. Perhaps @aramirezreyes would be willing to write physics tests once this module is merged?