Many of the shallow water advection operators appear to be wrong

[glwagner]

Shallow water advection operators don't look right:

Oceananigans.jl/src/Models/ShallowWaterModels/shallow_water_advection_operators.jl

Lines 30 to 36 in 5fbd8cd

	@inline div_hUu(i, j, k, grid, advection, solution) =
	1 / Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, momentum_flux_huu, advection, solution) +
	δyᵃᶜᵃ(i, j, k, grid, momentum_flux_hvu, advection, solution))
	@inline div_hUv(i, j, k, grid, advection, solution) =
	1 / Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, momentum_flux_huv, advection, solution) +
	δyᵃᶠᵃ(i, j, k, grid, momentum_flux_hvv, advection, solution))

In particular the advective fluxes are not multiplied by cell area, even though the flux differences are divided by cell volume. Surely this can't be right? Perhaps I am missing something.

The mass flux divergence is correct, however:

Oceananigans.jl/src/Models/ShallowWaterModels/shallow_water_advection_operators.jl

Lines 55 to 58 in 5fbd8cd

	@inline function div_Uh(i, j, k, grid, solution)
	1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_uᶠᶜᶜ, solution.uh) +
	δyᵃᶜᵃ(i, j, k, grid, Ay_vᶜᶠᶜ, solution.vh))
	end

But the tracer flux is again incorrect, missing a multiplication by area:

Oceananigans.jl/src/Models/ShallowWaterModels/shallow_water_advection_operators.jl

Lines 80 to 83 in 5fbd8cd

	@inline function div_Uc(i, j, k, grid, advection, solution, c)
	1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, transport_tracer_flux_x, advection, solution.uh, solution.h, c) +
	δyᵃᶜᵃ(i, j, k, grid, transport_tracer_flux_y, advection, solution.vh, solution.h, c))
	end

The functions that calculate velocity are not correct:

Oceananigans.jl/src/Models/ShallowWaterModels/shallow_water_advection_operators.jl

Lines 88 to 89 in 5fbd8cd

	@inline u(i, j, k, grid, solution) = @inbounds solution.uh[i, j, k] / solution.h[i, j, k]
	@inline v(i, j, k, grid, solution) = @inbounds solution.vh[i, j, k] / solution.h[i, j, k]

because they don't take into account the fact that uh is at cell faces in x, while h is at cell centers.

Finally, I think the momentum flux operators aren't quite right because they don't produce the correct high-order advection stencil:

Oceananigans.jl/src/Models/ShallowWaterModels/shallow_water_advection_operators.jl

Lines 14 to 24 in 5fbd8cd

	@inline momentum_flux_huu(i, j, k, grid, advection, solution) =
	@inbounds _advective_momentum_flux_Uu(i, j, k, grid, advection, solution.uh, solution.uh) / solution.h[i, j, k]
	@inline momentum_flux_hvu(i, j, k, grid, advection, solution) =
	@inbounds _advective_momentum_flux_Vu(i, j, k, grid, advection, solution.vh, solution.uh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)
	@inline momentum_flux_huv(i, j, k, grid, advection, solution) =
	@inbounds _advective_momentum_flux_Uv(i, j, k, grid, advection, solution.uh, solution.vh) / ℑxyᶠᶠᵃ(i, j, k, grid, solution.h)
	@inline momentum_flux_hvv(i, j, k, grid, advection, solution) =
	@inbounds _advective_momentum_flux_Vv(i, j, k, grid, advection, solution.vh, solution.vh) / solution.h[i, j, k]

I think what we really want is to express the notion that the velocity u advects the total momentum uh. For this we have to pass uh. / h in the place of the advecting velocity (the first entry) in functions like _advective_momentum_flux_Uu, rather than dividing by h after using second-order interpolation after the advective flux is calculated.

[glwagner]

@francispoulin does it also make sense to write these divergences as two-dimensional (in xy) since that's what we are restricted to for ShallowWaterModel? It might make the code clearer.

[francispoulin]

Very sorry for the problems that you found but I'm glad you found them. I believe when @ali-ramadhan and I put this together we were following other examples but I definitely should have been more careful.

Just so that I understand, instead of having momentum_flux_huu, advection and transport_tracer_flux_x we should have had something involving the area? I'm happy to help fix this where I can.

As for computing the velocity, I hope we can fix that soon as well. I know that ShallowWaterModel is a bit odd as we integrate the mass transports, not the velocities, but we do use the velocity a lot. I wonder if it's worth while computing the velocities (correctly) and then storing those. That should certainly help when we add in closure schemes, since those should be based on the velocities, for the most part.

One option would be to add model.velocities,u and something similar for v, and then access them when we need them. That has the unfortuante effect of storing 5 instead of 3 fields, so it would make things more memory intensive. I don't know if it's better to just compute the velocities everytime we need them?

[glwagner]

No need to apologize, just trying to help.

I think you can form an AbstractOperation for velocities when the model is constructed, eg:

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

These can be stored in ShallowWaterModel.velocities, and used in kernels. They have a getindex method so they'll work with WENO functions too. This doesn't use any extra memory. AbstractOperations perform the calculation uh/h with correct interpolation on the fly.

Just so that I understand, instead of having momentum_flux_huu, advection and transport_tracer_flux_x we should have had something involving the area? I'm happy to help fix this where I can.

We'd write a horizontal divergence either as

div(Q) = 1 / V * (δx(Ax * Qx) + δy(Ay * Qy))

or

div(Q) = 1 / Az * (δx(Δy * Qx) + δy(Δx * Qy))

I think this is written in the docstrings but doesn't appear to be reflected in the code. Correct me if I'm wrong.

[glwagner]

I don't quite understand how the examples can be right without the areas though, so I feel I might be missing something.

[francispoulin]

Maybe part of my confusion is on whether the solution fields, uh,vh,h, are cell averaged quantities or not. If they are then do we need to multipy by the area?

Maybe having docs on the finite volume method, as discussed previously, would help to clear some of this up?

[glwagner]

All quantities are cell averaged. Certainly it would be good to state this in the docs if it is not stated already.

[navidcy]

Well, it's mentioned in the docs.

[navidcy]

#1900 brings the finite-volume discussion further up in the Docs/Numerical Implementations.

[francispoulin]

#1900 brings the finite-volume discussion further up in the Docs/Numerical Implementations.

I think that's a great idea.

Also, it occurs to me that it would be nice to have another validation code for ShallowWaterModel. The Bickley jet example was good and we confirmed the growth rates, but perhaps a propagatoing equatorial Rossby wave would be fun. This is an exact solution to the equations and we can ensure that the phase speed matches that with theory. Maybe this would help us to ensure that all the integrals are done correctly. I don't think this would be interesting enough to become an example but it might interest some people.

[glwagner]

I think more validation is great.

Integrated cases are split into three categories:

1. Tests (eg the stuff in test_dynamics.jl for NonhydrostaticModel). These run during CI.
2. validation/. These are scientific validation cases that often require scientific interpretation or are expensive. These are similar to "Tests" but may lack a quantitative metric of success.
3. examples/. These are intended to showcase the API and library usage to users. They should not be used as tests, because they are very expensive to run (via Documenter) and to maintain (for one because they have a high standard for code quality).

I suggest adding bona fide Tests and validation, rather than examples, if we are interested in determining the correctness of the code.

[francispoulin]

@glwagner :

I am trying to follow the reasoning as to why the momentum flux does not need the area terms and I have an idea.

First, I have tried to follow the dependencies of the flux function and find the following:

div_hUu -> momentum_flux_huu, -> _advective_momentum_flux_Uu -> advective_momentum_flux_Uu

The final function is defined for either centered of upwinding schemes.

centered_advective_fluxes.jl has a defintion that shows it's proportional to Ax, and hence the area:

Oceananigans.jl/src/Advection/centered_advective_fluxes.jl

Line 15 in f1cd43f

@inline advective_momentum_flux_Uu(i, j, k, grid, scheme::Centered, U, u) = Axᶜᶜᶜ(i, j, k, grid) * _symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, U) * _symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, u)

upwind_biased_advective_fluxes.jl has a definition that shows it is proportional to Ax as well:

Oceananigans.jl/src/Advection/upwind_biased_advective_fluxes.jl

Line 24 in f1cd43f

return Axᶜᶜᶜ(i, j, k, grid) * upwind_biased_product(ũ, uᴸ, uᴿ)

Something similar is true for advective_tracer_flux, and those can be found in the same files.

Does this answer the question why there should not be any area terms in the flux?

If this convention is confusing, do we want to do something different?

[francispoulin]

Also, it occurs to me that dividing the difference terms by the volumes in ShallowWaterModel seems silly. What do people think of replacing those with areas instead?

[francispoulin]

@glwagner : I don't know that this approaches achieves the high order that we can achieve and I am happy to try something else.

Also, it would be nice to have a test that does this. The test that we have for advection assumes that h is constant initially, and that migth be too simple to see whether we always achieve the high order that we want to achieve.

[glwagner]

I'm closing this issue because I'm judging that it's not of current, timely relevance to Oceananigans development. If you would like to make it a higher priority or if you think the issue was closed in error please feel free to re-open.