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:

instantaneous_statistics.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: JuliaGPU/CUDA.jl#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.

instantaneous_statistics.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