Add mountain wave validation case

[glwagner]

cc @kaiyuan-cheng

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

[giordano]

Should the new validation script be run in CI? Maybe after #67, because that touches the Validation workflow

[glwagner]

Should the new validation script be run in CI? Maybe after #67, because that touches the Validation workflow

One issue with always running validation in CI is that there will probably need to be some heavy validation cases. So we may not be able to run them all (or at least not all the time... maybe infrequently, eg nightly or weekly)

Another thought is that I think its nice to implement validation cases as "examples" which get built into the documentation (Oceananigans does this to some extent, though we may be able to do an even better job here).

However one may want to have additional validation cases that are not examples. Some of the non-example validation cases could be lightweight enough to run regularly in CI.

Maybe we can flesh out a strategy that incorporates these considerations?

[navidcy]

Yes. Agree.

It's bit tricky to run the validation cases in a CI. Especially if validation cases also start plotting things and one needs to visually inspect to ensure that nothing changed. An idea would be that validation scripts finish at

run!(simulation)

plus perhaps a very simple test after (eg something to ensure the flow is not nonsense?)

simulation.model.tracers.temperature < 300

And then different scripts for plotting? It does make the validation directory a bit more complex...

[glwagner]

Possibly we want to implement validation tests as regression tests.

[kaiyuan-cheng]

This mountain wave case is not working too well in terms of speed and fidelity for the following reasons

1. Immersed boundary requires very high resolution to resolve the short mountain.
2. Open boundary condition is needed to improve the solution.

At least, the phase of the most significant wave looks not too bad!

[glwagner]

This mountain wave case is not working too well in terms of speed and fidelity for the following reasons

1. Immersed boundary requires very high resolution to resolve the short mountain.
2. Open boundary condition is needed to improve the solution.

At least, the phase of the most significant wave looks not too bad!

Great work @kaiyuan-cheng ! This is indeed a challenging test for z-coordinate, so it is quite useful. My two big comments are 1) can we see what the solution looks like at high resolution (I understand that we need to get GPU working for that --- I will look at that in parallel) and 2) can we investigate whether vertical grid stretching helps?

Some more detail below.

Here are my takeaways:

• for simulations that use a z-coordinate wherein the terrain is very low-lying, we may need very high resolution to resolve the physics correctly
• for low resolution simulations with low-lying terrain and/or deep domains, we will probably want to implement a terrain following coordinate
• (we should work on validation cases that don't have topography too)

here are my existing questions:

• does grid stretching (packing resolution near the bottom of the domain) improve accuracy?
• what resolution is required for accuracy (perhaps expressed as "how well is the terrain resolved")?
• what are the characteristics of simulations where we deem that we should move to a terrain following coordinate (horizontal resolution, vertical resolution, domain depth, topography characteristics)?
• "how much" can we stretch the grid before accuracy breaks down? The usefulness of stretching has a limit.
• How much fidelity do we buy by using a better representation of the terrain in a z-coordinate --- eg cut cells versus partial cells?

a few comments:

• if we are using a terrain-following coordinate, I think we should also probably move to a fully compressible formulation. The reason is that an anelastic solver for terrain-following coordinates may incur a similar expense as acoustic substepping. (I do think it should be theoretically possible to develop a fast enough solver, but that work is difficult and the gains will probably be marginal)
• I think a cut cell implementation is similar or harder than terrain following coordinate, so we shouldn't think of cut cells as a "silver bullet". We should pursue cut cells if partial cells work "ok" and cut cells are an optimization. (the best is of course to pursue both...)

At some point we should open an issue to discuss the pros and cons.

Tagging @OsKnoth, @bischtob because he basically told me all of the points that I made above 😂

[glwagner]

Also tagging @simone-silvestri since I suspect he will be intrigued by acoustic substepping AND terrain following coordinate implementations (whose work could possibly inform the ocean model as well) --- although my hunch is that maintaining restful hydrostatic balance is harder in the ocean due to (1) bigger topography relative to depth (it goes to the top) and (2) possibly, equation of state difficulties.

[glwagner]

@kaiyuan-cheng can you also push your changes to this PR so we can see your code?

[kaiyuan-cheng]

@kaiyuan-cheng can you also push your changes to this PR so we can see your code?

I pushed it last night. Here is the link. Let me know if you have trouble accessing it.

[kaiyuan-cheng]

This mountain wave case is not working too well in terms of speed and fidelity for the following reasons

1. Immersed boundary requires very high resolution to resolve the short mountain.
2. Open boundary condition is needed to improve the solution.

At least, the phase of the most significant wave looks not too bad!

Great work @kaiyuan-cheng ! This is indeed a challenging test for z-coordinate, so it is quite useful. My two big comments are 1) can we see what the solution looks like at high resolution (I understand that we need to get GPU working for that --- I will look at that in parallel) and 2) can we investigate whether vertical grid stretching helps?

Some more detail below.

Here are my takeaways:

• for simulations that use a z-coordinate wherein the terrain is very low-lying, we may need very high resolution to resolve the physics correctly

Yes, but this is not an ideal solution as it will increase the computational cost significantly (more vertical grids, together with smaller dt).

• for low resolution simulations with low-lying terrain and/or deep domains, we will probably want to implement a terrain following coordinate

I think we probably will need a terrain-following coordinate.

• (we should work on validation cases that don't have topography too)

Totally.

here are my existing questions:

• does grid stretching (packing resolution near the bottom of the domain) improve accuracy?

I guess stretching grid will help, but I expect it will still require a finer grid (and a smaller dt) to resolve the topography compared to a terrain-following coordinate, which will increase the computational cost.

• what resolution is required for accuracy (perhaps expressed as "how well is the terrain resolved")?

I think it depends on the slope of the terrain.

• what are the characteristics of simulations where we deem that we should move to a terrain following coordinate (horizontal resolution, vertical resolution, domain depth, topography characteristics)?

I guess a terrain-following coordinate would better capture weather phenomena related to topography, such as orographic precipitation and foehn winds. By "better", I mean more computationally efficient.

• "how much" can we stretch the grid before accuracy breaks down? The usefulness of stretching has a limit.
• How much fidelity do we buy by using a better representation of the terrain in a z-coordinate --- eg cut cells versus partial cells?

a few comments:

• if we are using a terrain-following coordinate, I think we should also probably move to a fully compressible formulation. The reason is that an anelastic solver for terrain-following coordinates may incur a similar expense as acoustic substepping. (I do think it should be theoretically possible to develop a fast enough solver, but that work is difficult and the gains will probably be marginal)

Agree.

• I think a cut cell implementation is similar or harder than terrain following coordinate, so we shouldn't think of cut cells as a "silver bullet". We should pursue cut cells if partial cells work "ok" and cut cells are an optimization. (the best is of course to pursue both...)

Agree.

At some point we should open an issue to discuss the pros and cons.

Tagging @OsKnoth, @bischtob because he basically told me all of the points that I made above 😂

[glwagner]

I guess a terrain-following coordinate would better capture weather phenomena related to topography, such as orographic precipitation and foehn winds. By "better", I mean more computationally efficient.

The key question is whether we require high resolution for other reasons (ie to resolve microphysical / turbulence interactions). In the case that we already would like to use high resolution for fidelity, then we do not accept any additional cost for resolving the topography.

However, if we decide that we would like to use a low resolution for a particular case, then it may not be worth it to increase the resolution massively merely to resolve the topography, and those cases are better served by a terrain following coordinate.

Note that the terrain following coordinate will also increase the computational cost somewhat (we might expect something like a 2-10x slowdown per time-step --- hopefully closer to 2x). So this increased cost has to be balanced against the cost savings of lower resolution. High resolution yields benefits everywhere (not only for resolving topography), so it would be better to use a higher resolution in the case that the two options are 1) high resolution with immersed boundary, which is appropriate for the terrain or 2) terrain following coordinate with lower resolution.

@OsKnoth has written papers about the benefits of cut cell immersed boundary methods so he might have more to add about the pros and cons.

Generally I would expect that immersed boundary starts to make more sense if one is considering steeper, larger terrain/obstacles (relative to domain height); for example simulations of flow around buildings, very large and steep mountains, etc. This Schar case is rather extreme in that we have a domain 20 km deep but the topography is only 250 m high.

[glwagner]

@kaiyuan-cheng can you add a sponge layer at the top, and sponge layers on the side of the domain (making the domain bigger)? If it is too expensive, you can try running on GPU.

(The sponge layer may only need to be on the vertical velocity component).

Also let us know what resolution you are using. I think the typical resolution is dx = 500, dz = 300.

[kaiyuan-cheng]

@kaiyuan-cheng can you add a sponge layer at the top, and sponge layers on the side of the domain (making the domain bigger)? If it is too expensive, you can try running on GPU.

(The sponge layer may only need to be on the vertical velocity component).

Also let us know what resolution you are using. I think the typical resolution is dx = 500, dz = 300.

I believe I have a sponge layer set up at the top of the domain (quite thick; 1/2 of the height, as suggested by a couple of studies running a similar simulation). As for the lateral boundary, we probably don't need a sponge layer, as the information of interest has not yet propagated across the entire domain. I also tested a setup with the horizontal domain size doubled, but the solution didn’t change much. For the vertical domain size, however, a larger model height gives a solution closer to the analytical solution, which suggests that the current sponge layer configuration does not eliminate wave reflection effectively.

The resolution I am using is dx = 1 km and dz = 25 m. I also tested a stretching vertical coordinate, which reduces nz by more than 50% and still runs pretty well.

[glwagner]

There's a helper in Oceananigans called ExponentialDiscretization (for making the z-coordinate) that might useful too!

[OsKnoth]

This means for a coarse equidistant vertical resolution (dz = 300m) the "immersed scheme" does not work. But see my simulation with a two-dimensional cut cell model.
200km x 19.5km, dx=500, dz=300, dt=20s, EndTime=3600s, figure -30km:30km x 0km:1km

[glwagner]

@kaiyuan-cheng I updated the validation case to use a discrete form Forcing rather than Relaxation. It seems that Relaxation has issues with AtmosphereModel on GPU --- we need tests for it.

[kaiyuan-cheng]

@navidcy Thank you for the updates!