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.

Sign up for GitHub

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

Jump to bottom

Fix measure type error on AMD GPU by passing type on loc function, and also fix broadcasting #150

Merged
merged 3 commits into from
Jul 29, 2024
Merged

Fix measure type error on AMD GPU by passing type on loc function, and also fix broadcasting #150

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
890cedd
Measure fixed by passing type on loc. Also rewrote broadcasting opera…
b-fg Jul 23, 2024
852dfdd
Fixed silly bug -.-
b-fg Jul 23, 2024
27b3f51
Added a ForwardDiff test with body
weymouth Jul 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 11 additions & 7 deletions src/Body.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,18 +28,22 @@ at time `t` using an immersion kernel of size `ϵ`.
See Maertens & Weymouth, doi:[10.1016/j.cma.2014.09.007](https://doi.org/10.1016/j.cma.2014.09.007).
"""
function measure!(a::Flow{N,T},body::AbstractBody;t=zero(T),ϵ=1) where {N,T}
a.V .= 0; a.μ₀ .= 1; a.μ₁ .= 0
a.V .= zero(T); a.μ₀ .= one(T); a.μ₁ .= zero(T)
Copy link
Collaborator

Choose a reason for hiding this comment

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

This is the same as fill, which must maintain the type. So do these actually need to be specified?

Copy link
Member Author

@b-fg b-fg Jul 23, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

I don't think they need but this is more correct imho.

Copy link
Collaborator

Choose a reason for hiding this comment

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

how can it be more correct if they compile to the same thing and waste characters? ;-)

Copy link
Member Author

Choose a reason for hiding this comment

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

x)

@fastmath @inline function fill!(μ₀,μ₁,V,d,I)
d[I] = sdf(body,loc(0,I),t)
d[I] = sdf(body,loc(0,I,T),t)
if abs(d[I])<2+ϵ
for i ∈ 1:N
dᵢ,nᵢ,Vᵢ = measure(body,loc(i,I),t)
dᵢ,nᵢ,Vᵢ = measure(body,loc(i,I,T),t)
V[I,i] = Vᵢ[i]
μ₀[I,i] = WaterLily.μ₀(dᵢ,ϵ)
μ₁[I,i,:] .= WaterLily.μ₁(dᵢ,ϵ)*nᵢ
for j ∈ 1:N
μ₁[I,i,j] = WaterLily.μ₁(dᵢ,ϵ)*nᵢ[j]
end
end
elseif d[I]<zero(T)
for i ∈ 1:N
μ₀[I,i] = zero(T)
end
elseif d[I]<0
μ₀[I,:] .= 0
end
end
@loop fill!(a.μ₀,a.μ₁,a.V,a.σ,I) over I ∈ inside(a.p)
Expand All @@ -60,7 +64,7 @@ end

Uses `sdf(body,x,t)` to fill `a`.
"""
measure_sdf!(a::AbstractArray,body::AbstractBody,t=0) = @inside a[I] = sdf(body,loc(0,I),t)
measure_sdf!(a::AbstractArray,body::AbstractBody,t=0) = @inside a[I] = sdf(body,loc(0,I,eltype(a)),t)

"""
NoBody
Expand Down
8 changes: 4 additions & 4 deletions src/Metrics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ Compute ``𝛚⋅𝛉`` at the center of cell `I` where ``𝛉`` is the azimuth
direction around vector `z` passing through `center`.
"""
function ω_θ(I::CartesianIndex{3},z,center,u)
θ = z × (loc(0,I)-SVector{3}(center))
θ = z × (loc(0,I,eltype(u))-SVector{3}(center))
n = norm2(θ)
n<=eps(n) ? 0. : θ'*ω(I,u) / n
end
Expand All @@ -95,7 +95,7 @@ pressure_force(sim) = pressure_force(sim.flow,sim.body)
pressure_force(flow,body) = pressure_force(flow.p,flow.f,body,time(flow))
function pressure_force(p,df,body,t=0,T=promote_type(Float64,eltype(p)))
df .= zero(eltype(p))
@loop df[I,:] .= p[I]*nds(body,loc(0,I),t) over I ∈ inside(p)
@loop df[I,:] .= p[I]*nds(body,loc(0,I,T),t) over I ∈ inside(p)
sum(T,df,dims=ntuple(i->i,ndims(p)))[:] |> Array
end

Expand All @@ -115,7 +115,7 @@ viscous_force(sim) = viscous_force(sim.flow,sim.body)
viscous_force(flow,body) = viscous_force(flow.u,flow.ν,flow.f,body,time(flow))
function viscous_force(u,ν,df,body,t=0,T=promote_type(Float64,eltype(u)))
df .= zero(eltype(u))
@loop df[I,:] .= -ν*∇²u(I,u)*nds(body,loc(0,I),t) over I ∈ inside_u(u)
@loop df[I,:] .= -ν*∇²u(I,u)*nds(body,loc(0,I,T),t) over I ∈ inside_u(u)
sum(T,df,dims=ntuple(i->i,ndims(u)-1))[:] |> Array
end

Expand All @@ -136,6 +136,6 @@ pressure_moment(x₀,sim) = pressure_moment(x₀,sim.flow,sim.body)
pressure_moment(x₀,flow,body) = pressure_moment(x₀,flow.p,flow.f,body,time(flow))
function pressure_moment(x₀,p,df,body,t=0,T=promote_type(Float64,eltype(p)))
df .= zero(eltype(p))
@loop df[I,:] .= p[I]*cross(loc(0,I)-x₀,nds(body,loc(0,I),t)) over I ∈ inside(p)
@loop df[I,:] .= p[I]*cross(loc(0,I,T)-x₀,nds(body,loc(0,I,T),t)) over I ∈ inside(p)
sum(T,df,dims=ntuple(i->i,ndims(p)))[:] |> Array
end
6 changes: 3 additions & 3 deletions src/util.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ Location in space of the cell at CartesianIndex `I` at face `i`.
Using `i=0` returns the cell center s.t. `loc = I`.
"""
@inline loc(i,I::CartesianIndex{N},T=Float32) where N = SVector{N,T}(I.I .- 1.5 .- 0.5 .* δ(i,I).I)
@inline loc(Ii::CartesianIndex) = loc(last(Ii),Base.front(Ii))
@inline loc(Ii::CartesianIndex,T=Float32) = loc(last(Ii),Base.front(Ii),T)
Base.last(I::CartesianIndex) = last(I.I)
Base.front(I::CartesianIndex) = CI(Base.front(I.I))
"""
Expand All @@ -145,8 +145,8 @@ Apply a vector function `f(i,x)` to the faces of a uniform staggered array `c` o
a function `f(x)` to the center of a uniform array `c`.
"""
apply!(f,c) = hasmethod(f,Tuple{Int,CartesianIndex}) ? applyV!(f,c) : applyS!(f,c)
applyV!(f,c) = @loop c[Ii] = f(last(Ii),loc(Ii)) over Ii ∈ CartesianIndices(c)
applyS!(f,c) = @loop c[I] = f(loc(0,I)) over I ∈ CartesianIndices(c)
applyV!(f,c) = @loop c[Ii] = f(last(Ii),loc(Ii,eltype(c))) over Ii ∈ CartesianIndices(c)
applyS!(f,c) = @loop c[I] = f(loc(0,I,eltype(c))) over I ∈ CartesianIndices(c)
"""
slice(dims,i,j,low=1)

Expand Down
22 changes: 19 additions & 3 deletions test/maintests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -221,7 +221,7 @@ end
for f ∈ arrays
p = zeros(4,5) |> f; measure_sdf!(p,body1)
I = CartesianIndex(2,3)
@test GPUArrays.@allowscalar p[I]≈body1.sdf(loc(0,I),0.0)
@test GPUArrays.@allowscalar p[I]≈body1.sdf(loc(0,I,eltype(p)),0.0)
end
end

Expand All @@ -247,14 +247,30 @@ end
WaterLily.L₂(u[:,:,2].-ue[:,:,2]) < 1e-4
end
end
@testset "ForwardDiff of TGV" begin
@testset "ForwardDiff" begin
function TGV_ke(Re)
sim,_ = TGVsim(Array;Re)
sim_step!(sim,π/100)
sum(I->WaterLily.ke(I,sim.flow.u),inside(sim.flow.p))
end
using ForwardDiff:derivative
@test derivative(TGV_ke,1e3) ≈ (TGV_ke(1e3+1)-TGV_ke(1e3-1))/2 rtol=1e-6

# Spinning cylinder lift generation
rot(θ) = SA[cos(θ) -sin(θ); sin(θ) cos(θ)] # rotation matrix
function spinning(ξ;D=16,Re=500)
C,R,U = SA[D,D],D÷2,1
body = AutoBody((x,t)->√(x'*x)-R, # circle sdf
(x,t)->rot(ξ*U*t/R)*(x-C)) # center & spin!
Simulation((2D,2D),(U,0),D;ν=U*D/Re,body,T=typeof(ξ))
end
function lift(ξ,t_end=1)
sim = spinning(ξ)
sim_step!(sim,t_end;remeasure=false)
WaterLily.total_force(sim)[2]/(ξ^2*sim.U^2*sim.L)
end
h = 1e-6
@test derivative(lift,2.0) ≈ (lift(2+h)-lift(2-h))/2h rtol=√h
end

function acceleratingFlow(N;T=Float64,perdir=(1,),jerk=4,mem=Array)
Expand Down Expand Up @@ -305,7 +321,7 @@ import WaterLily: ×
# test force routines
N = 32
p = zeros(N,N) |> f; df₂ = zeros(N,N,2) |> f; df₃ = zeros(N,N,N,3) |> f
@inside p[I] = loc(0, I)[2]
@inside p[I] = loc(0, I, eltype(p))[2]
body = AutoBody((x,t)->√sum(abs2,x.-(N/2))-N÷4,(x,t)->x)
force = WaterLily.pressure_force(p,df₂,body)
@test sum(abs,force/(π*(N/4)^2) - [0,1]) < 2e-3
Expand Down
Loading