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

Fixes error when reducing Fields with a condition on ImmersedBoundaryGrids #3440

Merged
merged 13 commits into from
Feb 1, 2024
3 changes: 3 additions & 0 deletions src/Architectures.jl
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,9 @@ arch_array(::CPU, a::CuArray) = Array(a)
arch_array(::GPU, a::Array) = CuArray(a)
arch_array(::GPU, a::CuArray) = a

arch_array(::CPU, a::BitArray) = a
arch_array(::GPU, a::BitArray) = CuArray(a)

arch_array(::GPU, a::SubArray{<:Any, <:Any, <:CuArray}) = a
arch_array(::CPU, a::SubArray{<:Any, <:Any, <:CuArray}) = Array(a)

Expand Down
7 changes: 7 additions & 0 deletions src/ImmersedBoundaries/immersed_reductions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,13 @@ const IF = AbstractField{<:Any, <:Any, <:Any, <:ImmersedBoundaryGrid}
@inline condition_operand(func::Function, op::IF, ::Nothing, mask) = ConditionalOperation(op; func, condition=NotImmersed(truefunc), mask)
@inline condition_operand(func::typeof(identity), op::IF, ::Nothing, mask) = ConditionalOperation(op; func, condition=NotImmersed(truefunc), mask)

@inline function condition_operand(func::Function, op::IF, cond::AbstractArray, mask)
tomchor marked this conversation as resolved.
Show resolved Hide resolved
arch = architecture(op.grid)
arch_condition = arch_array(arch, cond)
Copy link
Collaborator

Choose a reason for hiding this comment

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

It would be nice if we make this work also for fields, which are AbstractArrays, but do not have an arch_array method defined

Copy link
Collaborator

Choose a reason for hiding this comment

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

something like

Suggested change
arch_condition = arch_array(arch, cond)
arch_condition = arch_array(arch, parent(cond))

might suffice

Copy link
Member

Choose a reason for hiding this comment

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

on_architecture was meant for this

Copy link
Member

Choose a reason for hiding this comment

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

something like

might suffice

I don't like that. Why not keep it as a Field?

ni_condition = NotImmersed(arch_condition)
return ConditionalOperation(op; func, condition=ni_condition, mask)
end

@inline conditional_length(c::IF) = conditional_length(condition_operand(identity, c, nothing, 0))
@inline conditional_length(c::IF, dims) = conditional_length(condition_operand(identity, c, nothing, 0), dims)

Expand Down
24 changes: 22 additions & 2 deletions test/test_field_reductions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -228,13 +228,33 @@ trilinear(x, y, z) = x + y + z
c = Field((Center, Center, Nothing), grid)

set!(c, (x, y) -> y)
@test maximum(c) == CUDA.@allowscalar grid.yᵃᶜᵃ[1]
@test maximum(c) == grid.yᵃᶜᵃ[1]

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom((x, y) -> y < 0.5 ? - 0.6 : -0.4))
c = Field((Center, Center, Nothing), grid)

set!(c, (x, y) -> y)
@test maximum(c) == CUDA.@allowscalar grid.yᵃᶜᵃ[3]
@test maximum(c) == grid.yᵃᶜᵃ[3]

underlying_grid = RectilinearGrid(arch, size = (1, 1, 8), extent=(1, 1, 1))

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom((x, y) -> -3/4))
c = Field((Center, Center, Center), grid)

set!(c, (x, y, z) -> -z)
@test maximum(c) == Array(interior(c))[1, 1, 3]
c_condition = interior(c) .< 0.5
avg_c_smaller_than_½ = Array(interior(compute!(Field(Average(c, condition=c_condition)))))
@test avg_c_smaller_than_½[1, 1, 1] == 0.25

zᶜᶜᶜ = KernelFunctionOperation{Center, Center, Center}(znode, grid, Center(), Center(), Center())
ci = Array(interior(c)) # transfer to CPU
bottom_half_average_manual = (ci[1, 1, 3] + ci[1, 1, 4]) / 2
bottom_half_average = Average(c; condition=(zᶜᶜᶜ .< -1/2))
bottom_half_average_field = Field(bottom_half_average)
compute!(bottom_half_average_field)
bottom_half_average_array = Array(interior(bottom_half_average_field))
@test bottom_half_average_array[1, 1, 1] == bottom_half_average_manual
end
end
end