GPU compilation issue for body using multiple HashedLocators

[JulienL29]

Hi, I am trying to do a simulation with multiple moving parametric bodies.
Here is the script that I use :

using WaterLily,StaticArrays,Plots
using ParametricBodies

function make_sim(;nb_blade=1,L=48,Re=50000,U=1,λ=0.1f0,e=0.6f0,T=Float32,mem=Array)
    X,Y=2^4*floor(Int,L),2^3*floor(Int,L)
    O = SA[0.3f0*X,0.5f0*Y]           #center of rotor
    r=T(2L)
    ω=T(U/(λ*r))   

    function NACA(s)
        return 0.6f0*(0.2969f0(s)-0.126f0(s)^2-0.3516f0(s)^4+0.2843f0(s)^6-0.1036f0(s)^8)
    end
    function foil(s,t)
        return L*SA[(1-s)^2-0.33f0,NACA(1-s)]
    end

    body = sum(1:nb_blade) do blade
        ϕ = T(2π*blade/nb_blade)
        function map(x,t)
            θ = ω*(t*U/L)
            OP = SA[r*sin(θ+ϕ), -r*cos(θ+ϕ)] 
            P=O+OP
            ON=SA[0,-e*r]
            N=O+ON
            NP=P-N
            β=atan(NP[2],NP[1])-atan(OP[2],OP[1])
            R = SA[cos(β+θ+ϕ) sin(β+θ+ϕ); -sin(β+θ+ϕ) cos(β+θ+ϕ)] #rotation matrix for the angle of attack
            ξ = -R*(x-O-OP) 
            return SA[ξ[1],abs(ξ[2])] 
        end
        ParametricBody(foil,HashedLocator(foil,(0,1);T,mem);map)
    end
    Simulation((X,Y),(U,0),L;ν=U*L/Re,body,T,mem)
end

using CUDA
sim = make_sim(nb_blade=2,mem=CuArray)

When I use make_sim with nb_blade=1, the simulation runs correctly on GPU but when nb_blade=2 or more, I have this error message :

ERROR: GPU compilation of MethodInstance for (::WaterLily.var"#gpu_##kern_#567#252"{…})(::KernelAbstractions.CompilerMetadata{…}, ::CuDeviceArray{…}, ::CuDeviceArray{…}, ::CuDeviceArray{…}, ::CuDeviceMatrix{…}, ::CartesianIndex{…}) failed
KernelError: passing non-bitstype argument

I think that this is coming from using a sum of moving parametrics bodies but I don't understand of why.

[weymouth]

This example doesn't have multiple bodies. It has a map which is trying to collapse multiple parts of the domain down onto the same place. It must have a bug in it, since there are lots of examples which do this without error on the GPU.

You should either:

1. Start from a map on the repo that works and slowly add changes until it breaks. Then fix that mistake or post an issue specifically with that change. "This extremely simple map works, but when I add this single operation, it breaks (on the GPU). What gives?"
2. Actually use multiple bodies. Give each one it's own map so you don't need to do all this fairly complicate map logic, and then sum them together. Again, there are examples of this approach in the examples repo.

[JulienL29]

I created one map for each bodies and summed the bodies together like that :

function make_sim_1(;nb_blade=2,L=48,Re=50000,U=1,λ=0.1f0,e=0.6f0,T=Float32,mem=Array)
    # Map from simulation coordinate x to surface coordinate ξ
    X,Y=2^4*floor(Int,L/2),2^3*floor(Int,L/2)
    O = SA[0.3f0*X,0.5f0*Y]           #center of rotor
    r=T(2L)
    ω=T(U/(λ*r))
    
    function map_0(x,t)
        θ = ω*(t*U/L)
        OP = SA[r*sin(θ), -r*cos(θ)] 
        P=O+OP
        ON=SA[0,-e*r]
        N=O+ON
        NP=P-N
        β=atan(NP[2],NP[1])-atan(OP[2],OP[1])
        R = SA[cos(β+θ) sin(β+θ); -sin(β+θ) cos(β+θ)] #rotation matrix for the angle of attack
        ξ = -R*(x-O-OP) 
        return SA[ξ[1],abs(ξ[2])]    # reflect to positive y
    end
    #Define foil using NACA0012 profile equation: https://tinyurl.com/NACA00xx
    NACA_0(s) = 0.6f0*(0.2969f0(s)-0.126f0(s)^2-0.3516f0(s)^4+0.2843f0(s)^6-0.1036f0(s)^8)
    foil_0(s,t) = L*SA[(1-s)^2-0.33f0,NACA_0(1-s)]
    body0 = ParametricBody(foil_0,HashedLocator(foil_0,(0,1);T,mem);map=map_0)

    ϕ=T(2π/nb_blade)
    function map_1(x,t)        
        θ = ω*(t*U/L)
        OP = SA[r*sin(θ+ϕ), -r*cos(θ+ϕ)] 
        P=O+OP
        ON=SA[0,-e*r]
        N=O+ON
        NP=P-N
        β=atan(NP[2],NP[1])-atan(OP[2],OP[1])
        R = SA[cos(β+θ+ϕ) sin(β+θ+ϕ); -sin(β+θ+ϕ) cos(β+θ+ϕ)] #rotation matrix for the angle of attack
        ξ = -R*(x-O-OP) 
        return SA[ξ[1],abs(ξ[2])]    # reflect to positive y
    end
    # Define foil using NACA0012 profile equation: https://tinyurl.com/NACA00xx
    NACA_1(s) = 0.6f0*(0.2969f0(s)-0.126f0(s)^2-0.3516f0(s)^4+0.2843f0(s)^6-0.1036f0(s)^8)
    foil_1(s,t) = L*SA[(1-s)^2-0.33f0,NACA_1(1-s)]
    body1 = ParametricBody(foil_1,HashedLocator(foil_1,(0,1);T,mem);map=map_1)
    
    Simulation((X,Y),(U,0),L;ν=U*L/Re,body=body0+body1 ,T,mem)
end

But I still have the same GPU error and it works on CPU.
Is that because of the "map=map_1" part in ParametricBody(foil_1,HashedLocator(foil_1,(0,1);T,mem);map=map_1) ?
Because when I tried ParametricBody(foil_1,HashedLocator(foil_1,(0,1);T,mem);map) when I had only one map function named map for one body, it was working on GPU.

[weymouth]

I know it's slower, but please try to reduce these functions down the the simplest possible example that replicates the issue. This is how you debug these problems yourself. And if you're still stuck after you've got two single line geometries that won't work together, we'll know exactly how to help.

[JulienL29]

I tried with two circles with simple maps :

function make_sim_2(;L=48,Re=50000,U=1,T=Float64,mem=Array)

    X,Y=2^4*floor(Int,L/2),2^3*floor(Int,L/2)
    O1 = SA[0.3f0*X,0.5f0*Y]
    O2 = SA[0.6f0*X,0.5f0*Y]

    function map_1(x,t) 
        ξ = x-O1
        return SA[ξ[1],ξ[2]]   
    end
    circle1(s,t)=0.5f0L*SA[cos(s),sin(s)]
    body1 = ParametricBody(circle1,HashedLocator(circle1,(0,2π);T,mem);map=map_1)   

    function map_2(x,t) 
        ξ = x-O2
        return SA[ξ[1],ξ[2]]
    end
    circle2(s,t)=0.5f0L*SA[cos(s),sin(s)]
    body2 = ParametricBody(circle2,HashedLocator(circle2,(0,2π);T,mem);map=map_2)

    Simulation((X,Y),(U,0),L;ν=U*L/Re,body=body1+body2 ,T,mem)
end

using CUDA
sim = make_sim_2(;mem=CuArray,T=Float32)

It works when I use body1 and body2 separately but when I sum them together I still have the compilation error from the GPU

[weymouth]

Thanks Julien. I can confirm that this code works on the GPU for body=body1 or body2, but not body1+body2 on the GPU.

using WaterLily,StaticArrays,ParametricBodies,CUDA
function make_sim_2(;L::Int=48,Re=500,U=1,T=Float64,mem=Array)
    N,M = 2^3*L,2^2*L
    circle(s,t) = L*SA[cos(s)/2,sin(s)/2]
    body1 = HashedBody(circle,(0,2π);T,mem,map=(x,t)->x-SA[N÷3,M÷2])   
    body2 = HashedBody(circle,(0,2π);T,mem,map=(x,t)->x-SA[2N÷3,M÷2])
    Simulation((N,M),(U,0),L;ν=U*L/Re,body=body1+body2,T,mem)
end
sim = make_sim_2(;mem=CuArray,T=Float32)

My guess is that it is the HashedLocator messing this up. We have to copy the hash to the GPU and this must be messing up the body arithmetic.

[weymouth]

Yup. Passing the locate function explicitly solves the problem

function make_sim_3(;L::Int=48,Re=500,U=1,T=Float64,mem=Array)
    N,M = 2^3*L,2^2*L
    circle(s,t) = L*SA[cos(s)/2,sin(s)/2]
    locate(ξ,t) = atan(ξ[2],ξ[1])
    body1 = ParametricBody(circle,locate;T,map=(x,t)->x-SA[N÷3,M÷2])
    body2 = ParametricBody(circle,locate;T,map=(x,t)->x-SA[2N÷3,M÷2])
    Simulation((N,M),(U,0),L;ν=U*L/Re,body=body1+body2,T,mem)
end
sim = make_sim_3(;mem=CuArray,T=Float32)

So the problem is the HashedLocator.

I will try to hunt this down eventually, but I suggest using NURBS instead:

function make_sim_NURBS(;L::Int=48,Re=500,U=1,mem=Array)
    N,M = 2^3*L,2^2*L
    cps = SA_F32[1 1 0 -1 -1 -1  0  1 1
                0 1 1  1  0 -1 -1 -1 0]
    weights = SA_F32[1.,√2/2,1.,√2/2,1.,√2/2,1.,√2/2,1.] # A perfect circle requires...
    knots = SA_F32[0,0,0,1/4,1/4,1/2,1/2,3/4,3/4,1,1,1]  # non-uniform knot and weight vectors
    circle = NurbsCurve(cps,knots,weights)  # unit circle
    body1 = ParametricBody(circle;map=(x,t)->(x-SA[N÷3,M÷2])/L)  # add the scaling to the map
    body2 = ParametricBody(circle;map=(x,t)->(x-SA[2N÷3,M÷2])/L) # ditto
    Simulation((N,M),(U,0),L;ν=U*L/Re,body=body1+body2,T=Float32,mem)
end
sim = make_sim_NURBS(;mem=CuArray)

This is the much better option in general.