Make sure Oceananigans is type stable

[ali-ramadhan]

To start testing Oceananigans at low-precision arithmetic and with alternative types with @milankl:

I think the first thing I would need to do is making sure that Oceananigans is type stable. Meaning that you can pass a type T as an argument and this type is used throughout all computations. Even though you may only use Float64 atm, but having an entirely type-stable code that also means things like

function f(A::AbstractMatrix)
    ...
    0.5*A[i,j]
    ...
end

are written as

function f(A::Array{T,2}) where T
    ...
    one_half = T(0.5)
    one_half*A[i,j]

, makes sure that we can use Sherlog64 to understand where the problems are.

Ideally, all type conversion are explicitly written into the code rather than relying on promotion

[ali-ramadhan]

We use FT(0.5) and zero(FT) in a lot of places but never tested it outside of Float32 and Float64 so I'm sure there are places we missed.

For type-stable FFT we'll probably need something like https://github.com/milankl/coolFFT.jl

[milankl]

Hey everyone, as an example I can point you to Juls.jl which can be executed with type T as the first argument

julia> using Juls
julia> RunJuls(Float64);
Starting Juls on Thu, 05 Dec 2019 15:03:14 without output.
100% Integration done in 1.19s.

julia> RunJuls(Float32);
Starting Juls on Thu, 05 Dec 2019 15:03:19 without output.
100% Integration done in 0.72s.

julia> RunJuls(Float16);
Starting Juls on Thu, 05 Dec 2019 15:03:28 without output.
100% Integration done in 1min, 8s.

without relying on type promotion.

[milankl]

Concerning FFTs:
I've started once a conversation, but not sure how far we are now to have generic FFTs for any number format. We should check for what's happening with FourierTransforms.jl

[ali-ramadhan]

Yeah that would be cool. Might also not be a bad idea to just use coolFFT.jl as it's <40 lines and seems to work (if we restrict ourselves to powers of 2 for Nx, Ny, Nz).

[milankl]

Problem with coolFFT.jl is that it's fairly slow as it doesn't allow for plan_fft and plan_rfft, both should be supported by FourierTransforms.jl though, but I don't know how type-flexible they are.

[navidcy]

@ali-ramadhan did we resolve this?

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