Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove some type instabilities #3271

Merged
merged 4 commits into from
Sep 19, 2023
Merged

Remove some type instabilities #3271

merged 4 commits into from
Sep 19, 2023

Conversation

navidcy
Copy link
Collaborator

@navidcy navidcy commented Sep 18, 2023

This PR fixes a few type instabilities.

@navidcy
Copy link
Collaborator Author

navidcy commented Sep 18, 2023

There are a few more 0.0 in the code.

Question:

If there is

fill!(var, 0.0)

where var::Array and if we don't have access to grid, then is

fill!(var, 0)

safe to do?

@navidcy navidcy added the numerics 🧮 So things don't blow up and boil the lobsters alive label Sep 18, 2023
@navidcy navidcy changed the title Removes some more type instabilities Remove some type instabilities Sep 18, 2023
@glwagner
Copy link
Member

glwagner commented Sep 18, 2023

0 is almost always safe.

But neither case is really an issue.

The main issues are when the entirety of a heavy kernel (like one that calculates a tendency) may be promoted to higher precision. fill(var, 0) fill!(var, 0) is cheap and unlikely to affect performance.

That said it's just more precise to write fill(var, zero(eltype(var))) (this is exactly what you are trying to do) and therefore the preferred way to write it.

src/Models/Models.jl Outdated Show resolved Hide resolved
@@ -25,15 +25,15 @@ import Oceananigans.TimeSteppers: reset!
function reset!(model::AbstractModel)

for field in fields(model)
fill!(field, 0.0)
fill!(field, zero(model.grid))
end

for field in model.timestepper.G⁻
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Though assuming model.timestepper is obviously more problematic... (but beyond the scope of this PR I guess)

@milankl
Copy link
Collaborator

milankl commented Sep 18, 2023

The main issues are when the entirety of a heavy kernel (like one that calculates a tendency) may be promoted to higher precision. fill(var, 0) is cheap and unlikely to affect performance.

That said it's just more precise to write fill(var, zero(eltype(var))) (this is exactly what you are trying to do) and therefore the preferred way to write it.

You were meant to write fill! not fill right? For fill!(::AbstractArray,0) this is what happens

julia> a = Float32[1,2,3]
3-element Vector{Float32}:
 1.0
 2.0
 3.0

julia> @code_llvm fill!(a,0)
;  @ array.jl:346 within `fill!`
define nonnull {}* @"julia_fill!_127"({}* noundef nonnull align 16 dereferenceable(40) %0, i64 signext %1) #0 {
top:
;  @ array.jl:347 within `fill!`
; ┌ @ number.jl:7 within `convert`
; │┌ @ float.jl:159 within `Float32`
    %2 = sitofp i64 %1 to float
; └└

So the very first thing is that if eltype of the array and type of second argument aren't the same then it's converted (the %2 ... line)
Using fill!(var, zero(eltype(var))) then can skip this conversion (it's compiled away) but the result is the same. You can make and educated guess of the type of the zero, but honestly, I wouldn't even bother. It has probably zero impact on performance for any larger than a few elements and I find fill!(A,0) very clear to read too, so that's what I now always try to write.

@navidcy
Copy link
Collaborator Author

navidcy commented Sep 18, 2023

okie! how does it look now?

Comment on lines +52 to +55
@inline δxᶠᵃᵃ_η(i, j, k, grid, ::Type{Bounded}, η★::Function, args...) = ifelse(i == 1, zero(grid), δxᶠᵃᵃ(i, j, k, grid, η★, args...))
@inline δyᵃᶠᵃ_η(i, j, k, grid, ::Type{Bounded}, η★::Function, args...) = ifelse(j == 1, zero(grid), δyᵃᶠᵃ(i, j, k, grid, η★, args...))
@inline δxᶠᵃᵃ_η(i, j, k, grid, ::Type{RightConnected}, η★::Function, args...) = ifelse(i == 1, zero(grid), δxᶠᵃᵃ(i, j, k, grid, η★, args...))
@inline δyᵃᶠᵃ_η(i, j, k, grid, ::Type{RightConnected}, η★::Function, args...) = ifelse(j == 1, zero(grid), δyᵃᶠᵃ(i, j, k, grid, η★, args...))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you defined zero(grid) as the zero element of the eltype on the grid? Because I assume this to return a whole grid with zeros in it, similar as to how

julia> zero(rand(2,2))
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

As far as I understand these functions (because of @inline and provided indices they are applied to all elements in an array), it's probably wise to directly pass on the zero in the correct type to avoid many convert(T,0) though!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

zero(grid) returns zero(eltype(grid)). But maybe this is poor design. We use the pattern a lot so its nice to make it a little more readable...

I don't think there's a purpose to a grid full of zeros, or at least it's unclear exactly what that means, since we can't use a grid in arithmetic

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's also the zeros function which could be used as zeros(::Type{<:AbstractGrid},n,m,...) in order to create a grid full of zeros. I think the zero function is technically meant to give for any type the zero element so that A + zero(A) == A and A * zero(A) == zero(A)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hm... seems you are right

help?> zero
search: zero zeros iszero set_zero_subnormals get_zero_subnormals count_zeros RoundToZero leading_zeros trailing_zeros RoundFromZero finalizer UndefInitializer

  zero(x)
  zero(::Type)


  Get the additive identity element for the type of x (x can also specify the type itself).

  See also iszero, one, oneunit, oftype.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> zero(1)
  0

  julia> zero(big"2.0")
  0.0

  julia> zero(rand(2,2))
  2×2 Matrix{Float64}:
   0.0  0.0
   0.0  0.0

so perhaps we need to re-think about zero(grid) behavior (in a different PR tho). @glwagner what do you reckon?

@navidcy
Copy link
Collaborator Author

navidcy commented Sep 19, 2023

@milankl are you happy with the PR given that zero(grid) is indeed what you thought, i.e., zero(eltype(grid))?

@milankl
Copy link
Collaborator

milankl commented Sep 19, 2023

Yes, absolutely! Also in SpeedyWeather we use zeros(::Type{<:AbstractGrid},n::Integer) (with s!) to create a grid full with zeros.

@navidcy navidcy merged commit 7f461bd into main Sep 19, 2023
@navidcy navidcy deleted the ncc/type-instabilities branch September 19, 2023 15:52
navidcy added a commit that referenced this pull request Sep 20, 2023
* remove type instabilities

* some 0.0 -> 0 to avoid type instability

* 0 -> zero(grid) + some missing @inbounds + some code alignment

* no need to assume model has a grid
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
numerics 🧮 So things don't blow up and boil the lobsters alive
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants