Conversion of integers to floats in GeoUnit conversion

[albert-de-montserrat]

Is there any specific reason for converting integers to floats for the convert method corresponding to GeoUnits ?

I am asking because for the stress / strain / viscosity calculations the slowest part is by far the powers with float exponents, when in fact many times most of them are actual integers. Example:

using GeoParams; import GeoParams.fastpow
f, d,  A, FT, E, P, V, R, T, FE, TauII = tuple(rand(11)...)

function f1(f, r, d, p,  A, FT, E, P, V, R, T, FE, TauII, n)
    (fastpow(FT * TauII, n) *fastpow(f, r) * fastpow(d, p) * A * FT * exp((-E - P * V) / (R * T)) * inv(FE))
end

julia> @btime f1($(f, 1.0, d, 2.0,  A, FT, E, P, V, R, T, FE, TauII, 3.0)...)
  16.658 ns (0 allocations: 0 bytes)
0.009911001392750043

julia> @btime f1($(f, 1, d, 2,  A, FT, E, P, V, R, T, FE, TauII, 3)...)
  13.388 ns (0 allocations: 0 bytes)
0.009911001392750043

I would like to not force the conversion to floats in GeoUnits if possible, hopefully to scratch some performance gain here since computing the viscosity with diffusion and/or dislocation creep turns out to be one of the most expensive parts of the PT iterations.

[boriskaus]

I’m fine if you undo this automatic conversion

[boriskaus]

If this is the slowest part, should we perhaps precompute many of these factors? Grain size will rarely be set as spatially variable so one could precompute that. Same with many of the parameters in the thermal part

[albert-de-montserrat]

Pre-computing fastpow(f, r) and fastpow(d, p) would certainly help a lot, by doing that f1 goes down to about 6 ns, as long as n is an integer.

[boriskaus]

the issue is to find a way to tell the code when to precompute vs. when to do the full computation.

[albert-de-montserrat]

A somehow less intrusive change, which seems to perform almost as well as precomputing those factors, is adding a couple of checks:

function f2(f, r, d, p,  A, FT, E, P, V, R, T, FE, TauII, n)
    fr = if iszero(r)
        1.0
    elseif isone(r)
        f
    else
        fastpow(f, r)
    end
    dp = if iszero(p)
        1.0
    elseif isone(p)
        d
    else
        fastpow(d, p)
    end
    (fastpow(FT * TauII, n) * fr * dp * A * FT * exp((-E - P * V) / (R * T)) * inv(FE))
end

In [9]: @btime f1($(f, 1.0, d, 0.0,  A, FT, E, P, V, R, T, FE, TauII, 3.0)...)
  19.038 ns (0 allocations: 0 bytes)
0.0007530828398154311

In [10]: @btime f2($(f, 1.0, d, 0.0,  A, FT, E, P, V, R, T, FE, TauII, 3.0)...)        
  14.300 ns (0 allocations: 0 bytes)
0.0007530828398154311

The boost of pre-computing inv(R * T) seems minimal, in particular when n isa AbstractFloat), so I will leave this part as it is now.