Improve type hygene in wavelike

No allocations for Float32/64/Duals. Also tested nearfield and NKPanelSystem which were fine.

Author: weymouth

src/NKPanelSystem.jl

@@ -103,22 +103,23 @@ Ngk(X::SVector{3}) = Ngk(X...)
 # Wave-like disturbance
 function wavelike(x::T,y::T,z::T)::T where T
     (x≥0 || z≤-10) && return zero(T)                 # trivial case
-    R = √(-10/z-1)                                   # heuristic t limit
-    S = filter(s->-R<s<R,stationary_points(x,y))     # g'=0 points
+    xv,yv,zv = value.((x,y,z))                       # strip any Duals for ranges
+    R,Δg = √(-10/zv-1),typeof(xv).(5.4)              # heuristic angle & phase limits
+    S = filter(s->-R<s<R,stationary_points(xv,yv))   # g'=0 points
     f(t) = exp(z*(1+t^2))*sin(g(x,y,t))              # integrand
-    length(S)==0 && return quadgl(f,-R,R)            # easy case
-    rngs = finite_ranges(S,t->g(x,y,t),5.4,R)        # finite phase ranges
-    4complex_path(t->g(x,y,t)-im*z*(1+t^2),          # complex phase
-                  t->dg(x,y,t)-2im*z*t, rngs; f)     # it's derivative
+    length(S)==0 && return quadgl(f,-R,R)            # real-line case
+    rngs = finite_ranges(S,t->g(xv,yv,t),Δg,R)       # finite phase ranges
+    4complex_path(t->g(x,y,t)-im*z*(1+t^2),          # complex case: complex phase
+                  t->dg(x,y,t)-2im*z*t, rngs; f)     # ...and it's derivative
 end
 g(x,y,t) = (x+y*t)*⎷(1+t^2)               # phase function
 dg(x,y,t) = (x*t+y*(2t^2+1))/⎷(1+t^2)     # it's derivative
 ⎷(z::Complex) = π/2≤angle(z)≤π ? -√z : √z # move √ branch-cut
 ⎷(x) = √x
 
 # Return points where dg=0
-function stationary_points(x,y)
-    abs(y)≤√eps() && return (0.,)
+function stationary_points(x::T,y::T) where T
+    abs(y)≤√eps(T) && return (zero(T),)
     diff = x^2-8y^2
     diff≤√eps() ? (-x/4y,) : @. (-x+(-1,1)*√diff)/4y
 end

src/panels.jl

@@ -9,16 +9,17 @@ The parameter `(hᵤ+hᵥ)*c` sets the max deviation of the panel from a flat pl
 panel size in regions of high curvature. Errors if the adaptive routine gives more than `N_max` panels.
 """
 function panelize(surface,u₀=0.,u₁=1.,v₀=0.,v₁=1.;hᵤ=1.,hᵥ=hᵤ,c=0.05,
-                  transpose=false,flip=false,N_max=1000,verbose=false,kwargs...)
+                  transpose=false,flip=false,N_max=1000,verbose=false,T=Float64,kwargs...)
     # Transpose arguments u,v -> v,u
     transpose && return panelize((v,u)->surface(u,v),v₀,v₁,u₀,u₁,;hᵤ=hᵥ,hᵥ=hᵤ,c,
-                                 transpose=false,flip=!flip,N_max,verbose,kwargs...)
+                                 transpose=false,flip=!flip,N_max,verbose,T,kwargs...)
 
     # Check inputs and get output type
     (u₀≥u₁ || v₀≥v₁) && throw(ArgumentError("Need `u₀<u₁` and `v₀<v₁`. Got [$u₀,$u₁],[$v₀,$v₁]."))
     (hᵤ≤0 || hᵥ≤0 || c≤0) && throw(ArgumentError("Need positive `hᵤ,hᵥ,c`. Got $hᵤ,$hᵥ,$c."))
+    u₀,u₁,v₀,v₁ = T.((u₀,u₁,v₀,v₁))
     !(typeof(surface(u₀,v₀)) <: SVector) && throw(ArgumentError("`surface` function doesn't return an SVector."))
-    init = typeof(measure(surface,0.5u₀+0.5u₁,0.5v₀+0.5v₁,u₁-u₀,v₁-v₀))[]
+    init = typeof(measure(surface,(u₀+u₁)/2,(v₀+v₁)/2,u₁-u₀,v₁-v₀))[]
 
     # Get arcslength and inverse along bottom & top edges
     S₀,s₀⁻¹ = arclength(u->surface(u,v₀),hᵤ,c,u₀,u₁)
@@ -45,7 +46,7 @@ function panelize(surface,u₀=0.,u₁=1.,v₀=0.,v₁=1.;hᵤ=1.,hᵥ=hᵤ,c=0.
         dv = ve[2:end]-ve[1:end-1]             # panel heights
 
         # Measure panels along strip
-        @. measure(surface,u,v,du,dv;flip,kwargs...)
+        @. measure(surface,u,v,du,dv;flip,T,kwargs...)
     end
 
     # Check length and return as a Table

src/quad.jl

@@ -17,15 +17,15 @@ quadgl(f,a,b;x=xg32,w=wg32) = (b-a)/2*quadgl(t->f((b+a)/2+t*(b-a)/2);x,w)
 """
     complex_path(g,dg,rngs;atol=1e-3,γ=one,f=Im(γ*exp(im*g)))
 
-Estimate the integral `∫f(t)dt` from `t=[-∞,∞]` using a complex path, see Gibbs 2024. The 
-finite phase ranges `rngs` are integrated along the real line with QuadGK. The range end 
+Estimate the integral `∫f(t)dt` from `t=[-∞,∞]` using a complex path, see Gibbs 2024. The
+finite phase ranges `rngs` are integrated along the real line with QuadGK. The range end
 points where `flag=true` are integrated to ±∞ in the complex-plane using `±nsp(t₀,g,dg,γ)`.
 """
 @inline function complex_path(g,dg,rngs;γ=one,
     f = t->((u,v)=reim(g(t)); @fastmath γ(t)*exp(-v)*sin(u)))
 
     # Sum the flagged endpoints and interval contributions
-    val = zero(rngs[1][1])
+    val = zero(f(rngs[1][1]))
     for i in 1:2:length(rngs)
         (t₁,∞₁),(t₂,∞₂) = rngs[i],rngs[i+1]
         ∞₁ && (val -= nsp(t₁,g,dg,γ))
@@ -62,15 +62,14 @@ such that `|g(a)-g(aᵢ)|≈Δg`. Ranges do no overlap and limited to `±R`. "Un
 `fᵢ=false` if `aᵢ=±R`.
 """
 function finite_ranges(S::NTuple{N}, g, Δg, R; atol=Δg/10) where N
-    Sv, Rv, gv = map(value,S), value(R), t->value(g(t)) # no Duals
     # helper functions to offset the phase and flag if there's no root
-    dg(a) = t->abs(gv(a)-gv(t))-Δg
-    no(a,b) = abs(gv(a)-gv(b)) ≤ Δg+atol
+    dg(a) = t->abs(g(a)-g(t))-Δg
+    no(a,b) = abs(g(a)-g(b)) ≤ Δg+atol
     # find roots of dg using brackets (Order0) or secant method (Order1)
     fz0(a,b) = no(a,b) ? (return b, false) : (find_zero(dg(a), (a,b), Order0()), true)
-    fz1(a,b) = (isfinite(b) && no(a,b)) ? (return b, false) : (find_zero(dg(a), (a,a+copysign(1,b)), Order1(); atol), true)    
+    fz1(a,b) = (isfinite(b) && no(a,b)) ? (return b, false) : (find_zero(dg(a), (a,a+copysign(1,b)), Order1(); atol), true)
     # return flagged sub-range
-    (fz1(first(Sv), -Rv), mid_ranges(Val(N), Sv, fz0)..., fz1(last(Sv), Rv))
+    (fz1(first(S), -R), mid_ranges(Val(N), S, fz0)..., fz1(last(S), R))
 end
 using TupleTools
 mid_ranges(::Val{N}, S, fz) where N = TupleTools.vcat(ntuple(N-1) do i

test/runtests.jl

@@ -226,8 +226,9 @@ end
 @inline bruteW(x,y,z) = 4quadgk(t->exp(z*(1+t^2))*sin((x+y*t)*hypot(1,t)),-Inf,Inf,atol=1e-10)[1]
 @inline bruteN(x,y,z) = -2*(1-z/(hypot(x,y,z)+abs(x)))+NeumannKelvin.Ngk(x,y,z)
 using SpecialFunctions
+using ForwardDiff: Dual
 @testset "kelvin.jl" begin
-    @test NeumannKelvin.stationary_points(-1,1/sqrt(8))[1]≈1/sqrt(2)
+    @test NeumannKelvin.stationary_points(-1.,1/sqrt(8))[1]≈1/sqrt(2)
     @test NeumannKelvin.wavelike(10.,0.,-0.)==NeumannKelvin.wavelike(0.,10.,-0.)==0.
     @test 4π*bessely1(10)≈NeumannKelvin.wavelike(-10.,0.,-0.) atol=1e-5
 
@@ -241,7 +242,11 @@ using SpecialFunctions
     end
 
     @test @ballocations(NeumannKelvin.nearfield(-1.,0.,0.)) ≤ TEST_ALLOCS
+    @test @ballocations(NeumannKelvin.nearfield(Dual(-1.,0),Dual(0.,0),Dual(0.,0))) ≤ TEST_ALLOCS
+    @test @ballocations(NeumannKelvin.nearfield(Dual(-1f0,0),Dual(0f0,0),Dual(0f0,0))) ≤ TEST_ALLOCS
     @test @ballocations(NeumannKelvin.wavelike(-10.,1.,-1.)) ≤ TEST_ALLOCS
+    @test @ballocations(NeumannKelvin.wavelike(Dual(-10.,1),Dual(1.,1),Dual(-1.,1))) ≤ TEST_ALLOCS
+    @test @ballocations(NeumannKelvin.wavelike(Dual(-10f0,1),Dual(1f0,1),Dual(-1f0,1))) ≤ TEST_ALLOCS
 end
 
 function prism(h;q=0.2,Z=1)