Breaking revision

Cleaning up a few pain points.
- Get rid of PlotsExt
- Wrong sign on steadyforce
- Don't allow directsolve(FSPanelSystem)

Author: weymouth

Project.toml

@@ -1,7 +1,7 @@
 name = "NeumannKelvin"
 uuid = "7f078b06-e5c4-4cf8-bb56-b92882a0ad03"
 authors = ["Gabriel Weymouth"]
-version = "0.8.2"
+version = "0.9.0"
 
 [deps]
 AcceleratedKernels = "6a4ca0a5-0e36-4168-a932-d9be78d558f1"
@@ -25,13 +25,11 @@ TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NURBS = "dde13934-061e-461b-aa91-2c0fad390a0d"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 
 [extensions]
 NeumannKelvinMakieExt = "Makie"
 NeumannKelvinNURBSExt = "NURBS"
-NeumannKelvinPlotsExt = "Plots"
 NeumannKelvinGeometryBasicsExt = "GeometryBasics"
 
 [compat]

ext/NeumannKelvinPlotsExt.jl

@@ -77,7 +77,7 @@ the body panel area. Computation is accelerated when Threads.nthreads()>1.
 
 See also: [`cₚ`](@ref)
 """
-steadyforce(sys;S=bodyarea(sys)) = surface_integral(cₚ,sys)/S
+steadyforce(sys;S=bodyarea(sys)) = -surface_integral(cₚ,sys)/S
 @inline function surface_integral(f,sys)
     init = neutral = zero(eltype(sys.body.n))
     AK.mapreduce(+, sys.body, AK.get_backend(sys.body.q); init, neutral) do p

src/panel_method.jl

@@ -12,7 +12,7 @@ This function throws an error if the adaptive routine gives more than `N_max` pa
 function panelize(surface,u₀=0.,u₁=1.,v₀=0.,v₁=1.;hᵤ=1.,hᵥ=hᵤ,devlimit=0.05,
                   transpose=false,flip=false,N_max=1000,verbose=false,submerge=false,T=Float64,kwargs...)
     # Transpose arguments u,v -> v,u
-    transpose && return panelize((v,u)->surface(u,v),v₀,v₁,u₀,u₁,;hᵤ,hᵥ,devlimit,
+    transpose && return panelize((v,u)->surface(u,v),v₀,v₁,u₀,u₁,;hᵤ=hᵥ,hᵥ=hᵤ,devlimit,
                                  transpose=false,flip=!flip,submerge,N_max,verbose,T,kwargs...)
 
     # Check inputs and get output type

src/panels.jl

@@ -46,10 +46,9 @@ end
 
 Solve a panel system using a direct construction and solve such that the normal
 velocity boundary condition `∂ₙϕ(pᵢ,pⱼ)*qⱼ = -U⋅nᵢ` is satisfied on body panels.
+This function can not be used to solve an `FSPanelSystem`.
 
-**Note**: This function does not apply the FSBC and ignores free surface panels.
-
-*Note*: This function is memory (and therefore time) intensive for large number of
+**Note**: This function is memory (and therefore time) intensive for large number of
 panels N because it constructs the full N² matrix elements. It is *not* accelerated
 with a PanelTree, but *is* accelerated when Threads.nthreads()>1.
 
@@ -68,8 +67,9 @@ extrema(cₚ(sys))               # measure
 ```
 """
 function directsolve!(sys;verbose=true)
+    typeof(sys)<:FSPanelSystem && throw(ArgumentError("Cannot directsolve! an FSPanelSystem"))
     q = if verbose
-        @warn "This routine ignores free surface panels and is memory intensive. See help?>directsolve!."
+        @warn "This routine is memory intensive. See help?>directsolve!."
         @time influence(sys)\rhs(sys.body,sys.U)
     else
         influence(sys)\rhs(sys.body,sys.U)

src/solvers.jl

@@ -209,11 +209,6 @@ end
     badsurf = measure.((u,v)->SA[u,v,0],-4:1:2,(1/2:1:2)',1,1,flip=true)
     @test_throws ArgumentError FSPanelSystem(body,badsurf;ℓ)
 
-    # Direct solve ignores freesurf
-    directsolve!(sys,verbose=false)
-    sys2 = directsolve!(BodyPanelSystem(body,sym_axes=2),verbose=false)
-    @test sys.body.q == sys2.body.q
-
     # Setting ℓ=0 turns freesurf into reflection wall
     gmressolve!(sys)
     sys2 = gmressolve!(BodyPanelSystem(body,sym_axes=(2,3),wrap=PanelTree))
@@ -223,7 +218,7 @@ end
     # Setting ℓ>0 turns on freesurf, but it's slow to converge
     sys = FSPanelSystem(body,freesurf;sym_axes=2,ℓ,θ²=16)
     gmressolve!(sys,itmax=160)                      # should converge...
-    @test 2steadyforce(sys,S=1)[1] ≈ 6.2e-3 rtol=0.066 # Analytic Linear FSBC solution
+    @test -2steadyforce(sys,S=1)[1] ≈ 6.2e-3 rtol=0.066 # Analytic Linear FSBC solution
     mn,mx = extrema(ζ(sys))
     @test -2mn > mx # trough is much bigger than crest
 end
@@ -268,15 +263,15 @@ wigley(hᵤ;B=0.125,D=0.05,hᵥ=0.25) = measure.(
     sys = NKPanelSystem(spheroid(0.04),sym_axes=2,ℓ=0.5^2) |> directsolve!
     @test @ballocations(Φ($sys.body.x[1],$sys)) ≤ TEST_ALLOCS
     @test @ballocations(cₚ($sys.body.x[1],$sys)) ≤ TEST_ALLOCS
-    @test steadyforce(sys,S=1/2)[1] ≈ 6e-3 rtol=0.02
+    @test -steadyforce(sys,S=1/2)[1] ≈ 6e-3 rtol=0.02
 
     # Compare elliptical prism drag to Guevel/Baar (no WL contour)
     sys = NKPanelSystem(prism(0.1),sym_axes=2,ℓ=0.55^2) |> directsolve!
-    @test steadyforce(sys,S=1/2)[1] ≈ 0.062 rtol=0.03
+    @test -steadyforce(sys,S=1/2)[1] ≈ 0.062 rtol=0.03
 
     # Compare thin-ship wigley to Tuck 2008 (no WL contour)
     sys = NKPanelSystem(wigley(0.05),sym_axes=2,ℓ=0.51^2) |> directsolve!
-    @test steadyforce(sys)[1] ≈ 0.0088-0.003 rtol=0.03 # Remove ITTC Cf
+    @test -steadyforce(sys)[1] ≈ 0.0088-0.003 rtol=0.03 # Remove ITTC Cf
 end
 
 using NURBS,FileIO  # or whatever triggers the extension