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

[WIP] Add ShallowWater module #163

Draft
wants to merge 22 commits into
base: main
Choose a base branch
from
Draft

[WIP] Add ShallowWater module #163

wants to merge 22 commits into from

Conversation

navidcy
Copy link
Member

@navidcy navidcy commented Dec 9, 2020

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 navidcy marked this pull request as draft December 9, 2020 22:08
@navidcy
Copy link
Member Author

navidcy commented Dec 9, 2020

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

\begin{aligned}
\frac{\partial u}{\partial t} + \boldsymbol{u \cdot \nabla} u - f v & = - g \frac{\partial \eta}{\partial x} - \mathrm{D} u, \\
\frac{\partial v}{\partial t} + \boldsymbol{u \cdot \nabla} v + f u & = - \mathrm{D} v, \\
\frac{\partial \eta}{\partial t} + \boldsymbol{\nabla \cdot} ( \boldsymbol{u} h ) & = - \mathrm{D} \eta.
Copy link
Member

Choose a reason for hiding this comment

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

Should we write in conservative form?

Copy link
Member Author

@navidcy navidcy Dec 10, 2020

Choose a reason for hiding this comment

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

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.)

Copy link
Member Author

@navidcy navidcy Jan 11, 2021

Choose a reason for hiding this comment

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

Here they are in conservative form:

Screen Shot 2021-01-11 at 2 07 49 pm

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?

Copy link
Member

Choose a reason for hiding this comment

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

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

@glwagner
Copy link
Member

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

@navidcy navidcy self-assigned this Jan 11, 2021
@navidcy navidcy added the 🤥 enhancement New feature or request label Jan 11, 2021
@aramirezreyes
Copy link

aramirezreyes commented Mar 1, 2021

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
Copy link
Member Author

navidcy commented Mar 1, 2021

@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
Copy link
Member Author

navidcy commented Mar 4, 2021

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
Copy link
Member Author

navidcy commented Mar 4, 2021

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
Copy link
Member

glwagner commented Mar 4, 2021

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
🤥 enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2D Shallow Water module
3 participants