Results are wonky when running `examples/three_layer_constant_fluxes.jl` with free convection

[xkykai]

I ran examples/three_layer_constant_fluxes.jl with the command as follows: unbuffer ~/julia-1.6.1/bin/julia --project examples/three_layer_constant_fluxes.jl --size 256 256 128 --extent 512 512 256 --hours 192 --buoyancy-flux -4e-8 --momentum-flux 0 --name free_convection_8days_2 --animation, and the results of the simulation is weird, with the temperature profile blowing up as seen below:

https://user-images.githubusercontent.com/50624521/119716832-79d26680-be5d-11eb-98d2-7387d27278bd.mp4

I am using Julia v1.6.1, CUDA v3.2.1, Oceananigans v0.58.0, with the latest version of LESBRARY.jl

[glwagner]

With a negative buoyancy flux I believe you are heating the surface. I'm not sure why it would blow up though.

[ali-ramadhan]

Notice that the temperature is on the order of 1.3e6 in a lot of grid points (roughly the sum of 20 [the temperature] * 256^2 [the number of horizontal grid points]) but the simulations seem to continue fine for 8 days.

I think this is an upstream issue in CUDA.jl: https://github.com/JuliaGPU/CUDA.jl/issues/929

[xkykai]

There is an example here where I used --buoyancy-flux 3e-8 --momentum-flux 0, something odd is going on with a positive buoyancy flux as well.

https://user-images.githubusercontent.com/50624521/119718892-0da53200-be60-11eb-9e38-ffb4bfdf8346.mp4

[ali-ramadhan]

I can confirm this is still an issue with the latest version of packages 😕

[ali-ramadhan]

I tried to reproduce the issue in CUDA.jl by trying to average lots of views but no luck.

using Statistics
using CUDA

function average_views(N; arrays=1, val=1)
    as = [val .* ones(N+2, N+2, N+2) |> CuArray for _ in 1:arrays]
    vs = [view(as[i], 2:N+1, 2:N+1, 2:N+1) for i in 1:arrays]
    As = [zeros(1, 1, N+2) |> CuArray for _ in 1:arrays]
    Vs = [view(As[i], 1, 1, 2:N+1) for i in 1:arrays]
    [mean!(Vs[i], vs[i]) for i in 1:arrays]
    return [Array(Vs[i]) for i in 1:arrays]
end