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

Vector Invariant Formulation of the ShallowWaterModel #2522

Merged
merged 117 commits into from
Jun 3, 2022

Conversation

simone-silvestri
Copy link
Collaborator

@simone-silvestri simone-silvestri commented May 5, 2022

As for now the ShallowWaterModel is implemented in a conservative form.

This PR attempts to introduce a non-conservative form of the equation with a vector invariant formulation of advection

such that we can treat the advection with high order WENO reconstruction.
We should also think about introducing bathymetric terms ( and ) and viscous term (#2403)

with @francispoulin

@francispoulin
Copy link
Collaborator

Wonderful to see @simone-silvestri ! Please let me know how I can help.

One question, could we change the name of this PR to Vector Invariant Formuation of ShallowWaterModel? It seems like a good quality to highlight.

Also, as a test problem, could we use the same initial conditions as we have in the Bickley Jet example? Might make for a good comparison.

@simone-silvestri simone-silvestri changed the title [WIP] Non conservative form of the ShallowWaterModel [WIP] Vector Invariant Formulation of the ShallowWaterModel May 5, 2022
@simone-silvestri
Copy link
Collaborator Author

sounds like a good idea!

@simone-silvestri simone-silvestri added the numerics 🧮 So things don't blow up and boil the lobsters alive label May 5, 2022
@glwagner
Copy link
Member

glwagner commented May 5, 2022

You probably also want to fix #1866

I think what makes the most sense is to store velocities (AbstractOperations) inside the ShallowWaterModel so

velocities = (u = uh / h, v = vh / h)

(these have the correct location as well, which can be checked)

then use velocities.u to advect the total transport solution.uh in the advection operator.

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.

Looks great so far!

@glwagner
Copy link
Member

glwagner commented May 5, 2022

Perhaps

if formulation isa ConservativeFormulation
        solution = (uh = XFaceField(grid),
                    vh = YFaceField(grid),
                     h = CenterField(grid))
     else
        solution = (u = XFaceField(grid),
                    v = YFaceField(grid),
                    h = CenterField(grid))
    end

@simone-silvestri
Copy link
Collaborator Author

@francispoulin I have implemented the vector invariant / non conservative form of the equations. I have also separated the advection between momentum / velocity, height and tracers.

Now we have to think at some nice test cases, I ll start with the bickley jet experiment

@francispoulin
Copy link
Collaborator

That is great news @simone-silvestri, and very fast!

Yes, I think the Bickley jet is a good place to start. I'll think of others that would be fun.

When we put in closure we could do a gyre with a Munk layer. For now, we could do bottom drag and get a Stommel layer. I can look for a reference to compare this too if you like the idea?

@@ -32,6 +32,8 @@ function calculate_tendencies!(model::ShallowWaterModel)
model.grid,
model.gravitational_acceleration,
model.advection,
model.tracer_advection,
model.height_advection,
Copy link
Member

Choose a reason for hiding this comment

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

It might be better to use advection.tracers and advection.mass, so that we don't have to refactor this model too when we change the advection scheme implementations for nonhydrostatic and hydrostatic models in #2454

Copy link
Member

Choose a reason for hiding this comment

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

its ok if you don't want to, just trying to be lazy

@@ -115,7 +125,7 @@ compute!(ω)

# and finally set the "true" initial condition with noise,

set!(model, uh = uhⁱ)
formulation isa VectorInvariantFormulation ? set!(model, u = uhⁱ) : set!(model, uh = uhⁱ)
Copy link
Collaborator

Choose a reason for hiding this comment

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

In the case of vector invariant, should we set u = u^i and not the mass flux?

@@ -77,19 +83,19 @@ function calculate_interior_tendency_contributions!(tendencies,

barrier = Event(device(arch))

Guh_event = calculate_Guh_kernel!(tendencies.uh,
Guh_event = calculate_Guh_kernel!(tendencies[1],
Copy link
Collaborator

Choose a reason for hiding this comment

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

Given that the only difference between the next 3 steps is the index, would it be more compact to put these in a for loop over 1,2,3?

Copy link
Member

Choose a reason for hiding this comment

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

I like it written out personally!

Copy link
Collaborator

Choose a reason for hiding this comment

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

I am happy either way. Just a thought.

U.uh[i, j, k] += Δt * (γⁿ * Gⁿ.uh[i, j, k] + ζⁿ * G⁻.uh[i, j, k])
U.vh[i, j, k] += Δt * (γⁿ * Gⁿ.vh[i, j, k] + ζⁿ * G⁻.vh[i, j, k])
U.h[i, j, k] += Δt * (γⁿ * Gⁿ.h[i, j, k] + ζⁿ * G⁻.h[i, j, k])
U[1][i, j, k] += Δt * (γⁿ * Gⁿ[1][i, j, k] + ζⁿ * G⁻[1][i, j, k])
Copy link
Collaborator

Choose a reason for hiding this comment

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

Could loop over these three too, right?

U.uh[i, j, k] += Δt * γ¹ * G¹.uh[i, j, k]
U.vh[i, j, k] += Δt * γ¹ * G¹.vh[i, j, k]
U.h[i, j, k] += Δt * γ¹ * G¹.h[i, j, k]
U[1][i, j, k] += Δt * γ¹ * G¹[1][i, j, k]
Copy link
Collaborator

Choose a reason for hiding this comment

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

Loop?

Copy link
Member

Choose a reason for hiding this comment

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

I like it written out.

The main thing to fix here are the errors in units in advection terms etc

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.

Look great! This will be very easy to use and build on to add topography and closure schemes.

@francispoulin
Copy link
Collaborator

Should we merge this PR if all tests pass? Is there any other outstanding issues?

Regarding bathymetry properly included in ConservativeFormulation(), if there is more work needed there perhaps it belongs to a different PR (given the title of this PR)?

I agree that we should not wait to get bathymetry working with ConservativeFormulation().

@simone-silvestri and I were playing with a realistic global ocean model using the VectorInvariantForm and had some success. There were some issues when we included realistic topography, but I don't think that should hold us back. That is a hard problem and will take more time. This is a great addition and should allow for some very interesting simulations!

And we have 208 comments. That must be a record?

@navidcy
Copy link
Collaborator

navidcy commented Jun 1, 2022

OK, when test pass let's merge. And let's open an issue or PR for proper treatment of bathymetry in conservative formulation.

@navidcy
Copy link
Collaborator

navidcy commented Jun 1, 2022

Should we export

struct ConservativeFormulation end
struct VectorInvariantFormulation end

?

@francispoulin
Copy link
Collaborator

Should we export

struct ConservativeFormulation end
struct VectorInvariantFormulation end

?

If we don't export it then the user needs to import them manually?

I am fine either way but have a slight preference for exporting it, so it makes things a little easier for the user.

@navidcy
Copy link
Collaborator

navidcy commented Jun 1, 2022

I agree. So I pushed 8eaf216. If there are opposing opinions we can undo this.

With this export we can do:

julia> ShallowWaterModel(; gravitational_acceleration=1, grid, formulation=VectorInvariantFormulation(), momentum_advection = VectorInvariant())

instead of

julia> ShallowWaterModel(; gravitational_acceleration=1, grid, formulation=Oceananigans.Models.ShallowWaterModels.VectorInvariantFormulation(), momentum_advection = VectorInvariant())

@navidcy
Copy link
Collaborator

navidcy commented Jun 2, 2022

is there anything wrong with the regression?

Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>
@simone-silvestri
Copy link
Collaborator Author

I am not sure why it's not passing. I ll fix it, when it passes we can merge!

@navidcy
Copy link
Collaborator

navidcy commented Jun 3, 2022

omg, there is only one test pending...

smells like merging is close by...

visualize.jl Outdated Show resolved Hide resolved
@glwagner
Copy link
Member

glwagner commented Jun 3, 2022

Vorticity in ShallowWaterModel:

near_global_lat_lon_1440_600__fine_surface.6.mp4

and a single-layer HydrostaticFreeSurfaceModel:

near_global_lat_lon_1440_600_1_fine_surface.2.mp4

both with flat bathymetry. The free surface model is using the implicit free surface solver with a long time-step. Note the difference in the color scale.

(via @simone-silvestri and @francispoulin)

@simone-silvestri
Copy link
Collaborator Author

The hydrostatic simulation was done by @francispoulin

@simone-silvestri
Copy link
Collaborator Author

I 'll merge when the last test passes

@simone-silvestri simone-silvestri merged commit 9ebe5d6 into main Jun 3, 2022
@simone-silvestri simone-silvestri deleted the ss-fjp/non-conservative-shallow-water branch June 3, 2022 19:14
@navidcy
Copy link
Collaborator

navidcy commented Jun 3, 2022

Omg

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
numerics 🧮 So things don't blow up and boil the lobsters alive
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Can't we have a 2D simulation on the sphere?
6 participants