Update to Julia 1.6 and Oceananigans.jl v0.58

[ali-ramadhan]

Just upgrading LESbrary.jl to use the latest version of Oceananigans.jl (v0.57.1) but pausing to discuss how we should refactor the diagnostic field structs like TurbulentKineticEnergy and ShearProduction since Oceananigans.Fields has recently changed the AbstractField type.

We could subtype AbstractDataField (I did this for TurbulentKineticEnergy) but maybe a better approach would be to instead use KernelComputedField or even KernelFunctionOperation? Some of these operations are also available through https://github.com/tomchor/Oceanostics.jl.

[glwagner]

I think we should depend on Oceanostics.jl so we don't repeat work here.

If we can't use AbstractOperations directly (we may be able to use this now for turbulent kinetic energy), we should define KernelFunctionOperations rather than subtyping AbstractField.

[ali-ramadhan]

Agreed it would be good to rely on Oceanostics.jl to avoid duplicating work.

Should be pretty easy to just use the already available KernelComputedFields from Oceanostics.jl in this PR, then Oceanostics.jl could be refactored to make use of KernelFunctionOperation (or even just AbstractOperations) if @tomchor thinks it's a good idea as well.

[glwagner]

I think that'll work.

[ali-ramadhan]

This is now a Julia 1.6 + Oceananigans v0.58 PR.

[ali-ramadhan]

@tomchor @glwagner There's a different between how LESbrary.jl and Oceanostics.jl defines (vertical) shear production with U and V.

LESbrary.jl uses: (u - U) w ∂z(U) + (v - V) w ∂z(V) Oceanostics.jl uses: u w ∂z(U) + v w ∂z(V)

Note: I've excluded the w w ∂z(W) term which is not included in LESbrary.jl.

Which one should we be using?

[glwagner]

In Oceanostics you can pass u = model.velocities.u - AveragedField(model.velocities.u, dims=(1, 2)) which I think results in the same formula?

[ali-ramadhan]

That's true. So I guess we want to use (u - U) w ∂z(U) + (v - V) w ∂z(V) in LESbrary.jl to maintain existing behavior and not introduce any breaking changes into Oceanostics.jl.

[glwagner]

I think we might be able to use Oceanostics as is --- we just have to pass in u - U for our LESbrary calculation. @tomchor might be able to explain better.

[ali-ramadhan]

Yes we can use it as is but I wanted to confirm that we indeed want to compute (u - U) w ∂z(U) + (v - V) w ∂z(V) and not u w ∂z(U) + v w ∂z(V) in LESbrary.jl.

[glwagner]

Yes we can use it as is but I wanted to confirm that we indeed want to compute (u - U) w ∂z(U) + (v - V) w ∂z(V) and not u w ∂z(U) + v w ∂z(V) in LESbrary.jl.

Ah! Yes, (u - U) * w * ∂z(U) + (v - V) * w * ∂z(V) describes a conversion of kinetic energy from mean shear (U, V) to turbulent motions represented by a turbulent kinetic energy e = 1/2 * ((u - U)^2 + (v - V)^2 + w^2).

[ali-ramadhan]

Not sure if a Julia 1.6 issue but creating the JLD2 output writers with time averaged output seem to take much longer to compile now. Compiling run! also seems to take forever. All in all, 20+ minutes of compiling. See below.

And so CI here now takes over an hour 💀

Should we still merge?

We can tag LESbrary.jl v0.1.0 (not as a Julia package) as the latest version of LESbrary.jl that supports Julia 1.5 then merge this and continue with Julia 1.6 while looking into the cause of the length compile times.

cc @adelinehillier @xkykai in case you're using the three_layer_constant_fluxes.jl script.

---

Here I'm also timing the creation of the JLD2 output writers with then without halos:

[2021/05/20 16:08:16.979] INFO  Parsing command line arguments...
[2021/05/20 16:08:21.326] INFO  Mapping grid...
[2021/05/20 16:08:23.839] INFO  Enforcing boundary conditions...
[2021/05/20 16:08:24.072] INFO  Forcing and sponging tracers...
[2021/05/20 16:08:24.297] INFO  Building the wall model...
[2021/05/20 16:08:24.315] INFO  Framing the model...
[2021/05/20 16:08:27.359] WARN  Inflating model grid halo size to (3, 3, 3) and recreating grid. The model grid will be different from the input grid. To avoid this warning, pass halo=(3, 3, 3) when constructing the grid. -@-> /home/alir/.julia/packages/Oceananigans/LKGjq/src/Grids/automatic_halo_sizing.jl:41
[2021/05/20 16:08:37.918] INFO  Setting initial conditions...
[2021/05/20 16:10:05.150] INFO  Conjuring the simulation...
[2021/05/20 16:10:05.476] INFO  Strapping on checkpointer...
[2021/05/20 16:10:06.419] INFO  Squeezing out statistics...
[2021/05/20 16:10:35.014] INFO  Garnishing output writers...
[2021/05/20 16:10:35.337] INFO  xy
[2021/05/20 16:11:08.880] INFO  xz
[2021/05/20 16:11:27.394] INFO  yz
[2021/05/20 16:11:45.950] INFO  instant
[2021/05/20 16:11:46.419] INFO  time averaged
[2021/05/20 16:17:06.459] INFO  xy
[2021/05/20 16:17:06.463] INFO  xz
[2021/05/20 16:17:06.467] INFO  yz
[2021/05/20 16:17:06.471] INFO  instant
[2021/05/20 16:17:06.474] INFO  time averaged
[2021/05/20 16:21:22.512] INFO  Teaching simulation to run!...
[2021/05/20 16:29:20.512] INFO  [002.78%] iteration:     10, time: 10 seconds, Δt:   1 second, wall time: 19.244 minutes (1.924 minutes / time step)
[2021/05/20 16:29:25.302] INFO            └── u⃗_max: (1.25e-04, 6.29e-08, 4.75e-10) m/s, CFL: 7.81e-06, ν_max: 1.65e-09 m²/s, νCFL: 1.29e-11

[ali-ramadhan]

We can also disable time averaged outputs temporarily (or make it into a command line argument) since I don't think any of us use that data (we train on the instantaneous statistics).

[glwagner]

I'm fine to leave it out by default. Why would creating the time-averaged output writers take so much time? The outputs themselves don't get compiled until run!, right?

[ali-ramadhan]

Yeah I'm not sure either... seems like we should be able to reproduce in a minimal example but I don't have much experience debugging long compile times. Will ask on Julia Slack.