"Best" advection scheme for hydrostatic free-surface models?

[ali-ramadhan]

I'm curious to hear if there's a community consensus around what the "best" advection scheme is for hydrostatic models.

By "best" I am considering both accuracy and performance as I know the answer probably depends on the architecture with WENO being relatively quite a bit more expensive on CPUs, although I'm mostly thinking about GPUs.

For non-hydrostatic models, I feel like the consensus is that the best advection scheme on GPUs is regular WENO with order between 5 and 9, depending on the scenario. I think this still holds true for hydrostatic models on rectilinear grids?

Is there a consensus regarding advection schemes for hydrostatic models on curvilinear grids?

Based on @simone-silvestri's WENO paper, mainly figures 3 and 7, I also get the impression that the answer is somewhere between 5th and 9th order WENO (inclusive).

So probably the answer is some form of WENOVectorInvariant?

Maybe I'm also a bit confused, e.g. CliMA/ClimaOcean.jl#129 and looking at the constructor I'm wondering why vorticity gets WENO-9 by default with the others getting WENO-5. Should I be upwinding? Which vorticity stencil should I pick? Should I use the multi-dimensional stencil?

Oceananigans.jl/src/Advection/vector_invariant_advection.jl

Lines 183 to 224 in 82099ac

	function WENOVectorInvariant(FT::DataType = Float64;
	upwinding = nothing,
	vorticity_stencil = VelocityStencil(),
	order = nothing,
	vorticity_order = nothing,
	vertical_order = nothing,
	divergence_order = nothing,
	kinetic_energy_gradient_order = nothing,
	multi_dimensional_stencil = false,
	weno_kw...)
	if isnothing(order) # apply global defaults
	vorticity_order = nothing_to_default(vorticity_order, default = 9)
	vertical_order = nothing_to_default(vertical_order, default = 5)
	divergence_order = nothing_to_default(divergence_order, default = 5)
	kinetic_energy_gradient_order = nothing_to_default(kinetic_energy_gradient_order, default = 5)
	else # apply user supplied `order` unless overridden by more specific value
	vorticity_order = nothing_to_default(vorticity_order, default = order)
	vertical_order = nothing_to_default(vertical_order, default = order)
	divergence_order = nothing_to_default(divergence_order, default = order)
	kinetic_energy_gradient_order = nothing_to_default(kinetic_energy_gradient_order, default = order)
	end
	vorticity_scheme = WENO(FT; order=vorticity_order, weno_kw...)
	vertical_scheme = WENO(FT; order=vertical_order, weno_kw...)
	kinetic_energy_gradient_scheme = WENO(FT; order=kinetic_energy_gradient_order, weno_kw...)
	divergence_scheme = WENO(FT; order=divergence_order, weno_kw...)
	default_upwinding = OnlySelfUpwinding(cross_scheme = divergence_scheme)
	upwinding = nothing_to_default(upwinding; default = default_upwinding)
	schemes = (vorticity_scheme, vertical_scheme, kinetic_energy_gradient_scheme, divergence_scheme)
	N = maximum(required_halo_size(s) for s in schemes)
	FT = eltype(vorticity_scheme) # assumption
	return VectorInvariant{N, FT, multi_dimensional_stencil}(vorticity_scheme,
	vorticity_stencil,
	vertical_scheme,
	kinetic_energy_gradient_scheme,
	divergence_scheme,
	upwinding)
	end

I'm also curious why Centered advection in the vertical is used in ClimaOcean.jl: https://github.com/CliMA/ClimaOcean.jl/blob/887bc467eea1852a98b212b87b463a980d999323/src/OceanSimulations/OceanSimulations.jl#L31-L37

[simone-silvestri]

I guess you are referring to momentum advection and not tracer.
For tracer advection, the same considerations from the non-hydrostatic model (more or less) carry over to the hydrostatic model.

The cop-out answer for the momentum advection is that it depends.
Hydrostatic simulations can span many different regimes, from very coarse and non-eddying to eddy-resolving, so the advective term might be significant or not matter at all because it is completely dominated by Coriolis and viscous terms.

Performance considerations suggest using simple advection schemes like centered or upwind third order in the latter regime. Also, if the grid is very coarse, more than the implicit dissipation in WENO might be needed to give satisfactory results. It might be tricky to use a high-order WENO for very coarse simulations (where a 9-point stencil might span half the Atlantic Ocean).

In eddying regimes, where some of the turbulence is resolved, we found that a WENO approach allows a higher effective resolution because it is highly selective and removes the need for explicit dissipation. I am not sure which approach is better on a rectilinear grid, whether flux form or vector invariant. Last time I tested it, the vector invariant form was less dissipative, but if you reach very high resolution (near non-hydrostatic) probably the flux form is a better choice.

Regarding the different orders for the WENOVectorInvariant, in general, most of the energy is contained in the rotational motions, so using a higher order for vorticity flux is beneficial to preserve kinetic energy. Take care that 9th order is high! It should definitely not be used in low-resolution simulations.

The other terms carry less energy than the vorticity flux, so a 5th order is a good compromise between accuracy and performance. Note that vertical profiles are typically very smooth in a hydrostatic simulation with not very much grid cells, so an order higher than 5th for vertical advection is a bit overkill, especially because all the gradients are near vertical boundaries where the order is lower anyway (this is also a valid consideration for tracer advection). Not needing a high order in the vertical is also quite positive for performance since vertical advection is by far the most expensive computation on the GPU given the memory layout (k is the slowest moving index)

In terms of details for the WENOVectorInvariant, aside from the order that might depend on resolution, and requires a bit of experimenting (mainly decrease vorticity order if it the simulation crashes), the defaults are always the best option.
multi-dimensional stencils implement a 2D horizontal stencil instead of a 1D stencil that performs a centered 5th-order WENO reconstruction of vorticity, divergence and kinetic energy in the tangential direction. In theory, it is more accurate, but it is very computationally expensive and doesn't seem to improve the results much, so I haven't found a good use for this reconstruction.
The upwinding keyword argument is a bit of a vestigial argument because it allows two different forms of unwinding, one that ensures energy dissipation locally and one that might lead to anti-dissipative behavior locally (and thus might be unstable). There might be arguments to remove that option since I haven't really found a practical use for it.

[ali-ramadhan]

Thank you so much @simone-silvestri for all the explanations! Definitely appreciate that it depends especially on the dynamics and the significance of the advective terms.

Also thanks for sharing your experiences with the different flags/options.

I'll be playing around with advection schemes so I'll post back with anything new I learn.

[glwagner]

True! I'd be curious to see 9th order. I think 5th order is generally too diffusive. Also I may be reading your test wrong, but it seemed like in one case you used 7th order horizontal tracer advection + Centered vertical advection, while in the second case it was 5th order all around. Maybe try 7th order all around to just test the change from Centered to 7th order in the vertical?

We should change the default for WENOVectorInvariant to return the scheme we recommend rather than some random other scheme. @simone-silvestri do you know why this is not currently the case? Is it because the scheme we recommend has changed but we haven't yet updated the default?

[ali-ramadhan]

True. I re-ran the benchmarks with more consistent advection schemes to look at Centered vs. WENO-7 in the vertical.

The difference is bigger now. Looking at the median time it's ~185 ms -> ~245 ms so an increase of ~32%. Probably not catastrophic but not small either. However, I used 128 vertical levels which may be more than the usual hydrostatic model. But 32% increased time-to-solution may not be worth it if WENO isn't really improving a smooth solution in the vertical.

---

Centered() in the vertical

        momentum_advection = VectorInvariant(;
            vorticity_scheme = WENO(; order = 7),
            divergence_scheme = WENO(; order = 7),
            kinetic_energy_gradient_scheme = WENO(; order = 7),
            vertical_scheme = Centered()
        ),
        tracer_advection = FluxFormAdvection(
            WENO(; order = 7),
            WENO(; order = 7),
            Centered()
        ),

BenchmarkTools.Trial: 100 samples with 1 evaluation.
 Range (min … max):  184.514 ms … 207.852 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     185.723 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   188.241 ms ±   6.232 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▁▅▄ ▄▁                                                        
  ████▇██▇▇▁▅█▅▁▁▅▇▁▁▁▅▁▁▁▁▁▅▁▁▁▁▁▁▁▁▁▁▁▅▁▁▅▁▁▁▅▅▁▁▁▁▅▁█▁▁▁▁▅▅▅ ▅
  185 ms        Histogram: log(frequency) by time        208 ms <

 Memory estimate: 6.19 MiB, allocs estimate: 23135.

WENO(; order = 7) in the vertical

        momentum_advection = VectorInvariant(;
            vorticity_scheme = WENO(; order = 7),
            divergence_scheme = WENO(; order = 7),
            kinetic_energy_gradient_scheme = WENO(; order = 7),
            vertical_scheme = WENO(; order = 7)
        ),
        tracer_advection = FluxFormAdvection(
            WENO(; order = 7),
            WENO(; order = 7),
            WENO(; order = 7)
        ),

BenchmarkTools.Trial: 100 samples with 1 evaluation.
 Range (min … max):  242.289 ms … 277.566 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     245.354 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   251.157 ms ±  10.572 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                                                              
  █▇█▅▅▆▃▄▃▁▃▃▃▁▃▃▁▁▁▁▁▁▁▃▁▃▃▃▃▃▃▁▁▁▁▃▃▁▃▃▅▃▁▃▃▁▁▃▃▁▁▁▁▃▁▁▁▁▃▁▃ ▃
  242 ms           Histogram: frequency by time          277 ms <

 Memory estimate: 6.19 MiB, allocs estimate: 23134.

[glwagner]

Super useful @ali-ramadhan . @simone-silvestri didn't you find more catastrophic slow downs?

However, I used 128 vertical levels which may be more than the usual hydrostatic model.

Oceananigans is good at high resolution, so I'd expect this to be close to typical for a global simulation. For regional simulations people will probably go even higher.

[glwagner]

Another note is that ClimaOcean uses 7th order tracer advection right now, but there are reasons to believe that we should be using higher-order tracer advection than momentum advection.

Suffice to say that @simone-silvestri and @jm-c are actively working on this a bit (in particularly characterizing numerical dissipation for different choices) but I hope that our current recommendation will always be encoded in the default for ClimaOcean.

[simone-silvestri]

Super useful @ali-ramadhan . @simone-silvestri didn't you find more catastrophic slow downs?

However, I used 128 vertical levels which may be more than the usual hydrostatic model.

Oceananigans is good at high resolution, so I'd expect this to be close to typical for a global simulation. For regional simulations people will probably go even higher.

Huh, I remember worse. The new upwinding stencils probably helped in this regard. This might warrant us having WENO5 as the default in the vertical instead of Centered. Centered seemed to work, but it might cause oscillation in particular cases that WENO5 will not cause. I haven't yet found a case where the solution crashes with Centered in the vertical, but if the cost is not too much, it is probably better to have a safer and more accurate advection scheme in the vertical as the default.