Merge pull request #150 from WaterLily-jl/fix_measure

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

src/Body.jl

@@ -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)
     @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)
@@ -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

src/Metrics.jl

@@ -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
@@ -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
 
@@ -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
 
@@ -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

src/util.jl

@@ -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))
 """
@@ -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)
 

test/maintests.jl

@@ -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
 
@@ -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)
@@ -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