Redefines equality (==) for Grids

[navidcy]

This PR supersedes #2019 and closes #2018.

[navidcy]

@tomchor what do you think about this?

(I didn't wanna make changes to your branch in case you are using it for work, so I branched of there...)

[navidcy]

If this seem to be doing the job @tomchor and if others are happy I can write up few tests.

[tomchor]

@tomchor what do you think about this?

(I didn't wanna make changes to your branch in case you are using it for work, so I branched of there...)

Thanks, that very thoughtful. I am indeed using that branch for research :)

The PR also looks very good! This is exactly what I had in mind 👍

Sorry I haven't done this yet in my original PR. I was actually planning on doing it this weekend but I got busy.

Should we add an equality test for grids? I think we can get away with something as simple as @test grid == deepcopy(grid) no? (This used to return false for vertically stretched grids.)

[navidcy]

Should we add an equality test for grids? I think we can get away with something as simple as @test grid == deepcopy(grid) no? (This used to return false for vertically stretched grids.)

We should! But I'll try to test a few cases... different type of grids and on different architectures.

[tomchor]

Assuming this was already tested for GPUs after the removal of @allowscalar I think this is good! Thanks so much for the help, @navidcy.

I do think it's good to add tests before merging (at least for the grids that have arrays in them for both GPU and CPUs), but I'll leave that up to you.

[glwagner]

This is good after changes.

My only question is whether another implementation might be something like

function Base.:(==)(grid1::AbstractGrid, grid2::AbstractGrid)
    !isa(grid2, typeof(grid1).name.wrapper) && return false
    topology(grid1) == topology(grid2) && return false

    x1, y1, z1 = nodes((Face, Face, Face), grid1)
    x2, y2, z2 = nodes((Face, Face, Face), grid2)

    return x1 == x2 && y1 == y2 && z1 == z2
end

This has different behavior than what's implemented here. In particular, two grids will be identified as equal even if one has "number spacing" (eg a constant) whereas the other has array-spacing with all the values equal to that number. In other words, while

julia> 1.0 == ones(2)
false

we may want to identify this case with numeric equality for a grid.

The key here is that two grids are "equal" if they have the same nodes. The other properties (like spacings) are really implied by the placement of the nodes.

[navidcy]

This is good after changes.

My only question is whether another implementation might be something like

function Base.:(==)(grid1::AbstractGrid, grid2::AbstractGrid)
    !isa(grid2, typeof(grid1).name.wrapper) && return false
    topology(grid1) == topology(grid2) && return false

    x1, y1, z1 = nodes((Face, Face, Face), grid1)
    x2, y2, z2 = nodes((Face, Face, Face), grid2)

    return x1 == x2 && y1 == y2 && z1 == z2
end

This has different behavior than what's implemented here. In particular, two grids will be identified as equal even if one has "number spacing" (eg a constant) whereas the other has array-spacing with all the values equal to that number. In other words, while

julia> 1.0 == ones(2)
false

we may want to identify this case with numeric equality for a grid.

The key here is that two grids are "equal" if they have the same nodes. The other properties (like spacings) are really implied by the placement of the nodes.

OK, perhaps then your suggestion is better.
What about when two grids are otherwise the same but one lives on GPU and one on CPU. Will the above render them equal? I guess no since, e.g., z1 == z2 will return false if one is an Array and other a CuArray.

[glwagner]

This is good after changes.
My only question is whether another implementation might be something like

function Base.:(==)(grid1::AbstractGrid, grid2::AbstractGrid)
    !isa(grid2, typeof(grid1).name.wrapper) && return false
    topology(grid1) == topology(grid2) && return false

    x1, y1, z1 = nodes((Face, Face, Face), grid1)
    x2, y2, z2 = nodes((Face, Face, Face), grid2)

    return x1 == x2 && y1 == y2 && z1 == z2
end

This has different behavior than what's implemented here. In particular, two grids will be identified as equal even if one has "number spacing" (eg a constant) whereas the other has array-spacing with all the values equal to that number. In other words, while

julia> 1.0 == ones(2)
false

we may want to identify this case with numeric equality for a grid.
The key here is that two grids are "equal" if they have the same nodes. The other properties (like spacings) are really implied by the placement of the nodes.

OK, perhaps then your suggestion is better. What about when two grids are otherwise the same but one lives on GPU and one on CPU. Will the above render them equal? I guess no since, e.g., z1 == z2 will return false if one is an Array and other a CuArray.

Actually, the julia gods have decided that CPU and GPU arrays can be numerically equal:

julia> using CUDA

julia> a = rand(3)
3-element Vector{Float64}:
 0.3492331921297629
 0.4354818891776633
 0.684771954777823

julia> a == CuArray(a)
true

Since we use == on the nodes (which are always arrays), we retain the behavior that a CPU grid and a GPU grid can be numerically identical via ==.

[navidcy]

@tomchor, can you confirm that the PR at its current stage actually solves the problem you were facing in #2018?

[tomchor]

@tomchor, can you confirm that the PR at its current stage actually solves the problem you were facing in #2018?

Seems to do the trick! Tested for CPUs and GPUs. 🚀

[navidcy]

All right!