Fix branch 3 direction (#241)

* Fix branch 3 direction

* Fix tests

Author: 1-Bart-1

src/filament.jl

@@ -84,17 +84,19 @@ function velocity_3D_bound_vortex!(
         @debug "inside core radius"
         @debug "distance from control point to filament: $(nr1Xr0 / nr0)"
 
-        # Project onto core radius
-        cross3!(r2Xr0, r2, r0)
         nr0sq = nr0 * nr0
-        nr2Xr0 = norm3(r2Xr0)
         d_r1_r0 = dot3(r1, r0)
         d_r2_r0 = dot3(r2, r0)
+        r_rad = r1Xr0
+        @inbounds for k in 1:3
+            r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq
+        end
+        nr_rad = norm3(r_rad)
         @inbounds for k in 1:3
             r1_proj[k] = d_r1_r0 * r0[k] / nr0sq +
-                         epsilon * r1Xr0[k] / nr1Xr0
+                         epsilon * r_rad[k] / nr_rad
             r2_proj[k] = d_r2_r0 * r0[k] / nr0sq +
-                         epsilon * r2Xr0[k] / nr2Xr0
+                         epsilon * r_rad[k] / nr_rad
         end
         cross3!(r1_projXr2_proj, r1_proj, r2_proj)
 
@@ -289,13 +291,13 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     else
         r1_proj = work_vectors[4]
         cross_tmp = work_vectors[5]
-        # temp = r1/nr1 - Vf
+        nVfsq = nVf * nVf
         @inbounds for k in 1:3
-            cross_tmp[k] = r1[k]/nr1 - Vf[k]
+            cross_tmp[k] = r1[k] - d_r1_Vf * Vf[k] / nVfsq
         end
         n_tmp = norm3(cross_tmp)
         @inbounds for k in 1:3
-            r1_proj[k] = d_r1_Vf * Vf[k] +
+            r1_proj[k] = d_r1_Vf * Vf[k] / nVfsq +
                          epsilon * cross_tmp[k] / n_tmp
         end
         cross3!(cross_tmp, r1_proj, Vf)

test/filament/test_bound_filament.jl

@@ -74,7 +74,7 @@ end
         filament = BoundFilament{Float64}()
         reinit!(filament, [0.0, 0.0, 0.0], [1e6, 0.0, 0.0])
         control_point = [5e5, 1.0, 0.0]
-        
+
         velocity_3D_bound_vortex!(
             induced_velocity,
             filament,
@@ -85,8 +85,8 @@ end
         )
         @test !any(isnan.(induced_velocity))
         @test isapprox(induced_velocity[1], 0.0, atol=1e-8)
-        @test abs(induced_velocity[2]) < 1e-8
-        @test isapprox(induced_velocity[3], 0.0)
+        @test isapprox(induced_velocity[2], 0.0, atol=1e-8)
+        @test abs(induced_velocity[3]) > 1e-10
     end
 
     @testset "Point Far from Filament" begin
@@ -133,35 +133,35 @@ end
     @testset "Around Core Radius" begin
         filament = create_test_filament()
         delta = 1e-5
-        
+
         points = [
             [0.5, core_radius_fraction - delta, 0.0],
             [0.5, core_radius_fraction, 0.0],
             [0.5, core_radius_fraction + delta, 0.0]
         ]
-        
+
         velocities = [zeros(3) for p in points]
         [
             velocity_3D_bound_vortex!(velocities[i], filament, p, gamma, core_radius_fraction, work_vectors)
             for (i, p) in enumerate(points)
         ]
-        
+
         # Check for NaN and finite values
         @test all(!any(isnan.(v)) for v in velocities)
         @test all(all(isfinite.(v)) for v in velocities)
-        
+
         # Check magnitude is maximum at core radius
         @test norm(velocities[2]) > norm(velocities[1])
         @test norm(velocities[2]) > norm(velocities[3])
-        
+
         # Check continuity around core radius
         @test isapprox(velocities[1], velocities[2], rtol=1e-2)
-        
+
         # Check non-zero velocities
         @test !all(isapprox.(velocities[1], zeros(3), atol=1e-10))
         @test !all(isapprox.(velocities[2], zeros(3), atol=1e-10))
         @test !all(isapprox.(velocities[3], zeros(3), atol=1e-10))
-        
+
         # Check symmetry
         v_neg = zeros(3)
         velocity_3D_bound_vortex!(
@@ -174,4 +174,80 @@ end
         )
         @test isapprox(velocities[2], -v_neg)
     end
+
+    @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin
+        filament = create_test_filament()
+        r0 = [1.0, 0.0, 0.0]
+
+        for d in (0.25, 0.5, 0.99, 1.0, 1.01, 2.0) .* core_radius_fraction
+            for phi in (0.0, π/4, π/2, π, -π/3)
+                p = [0.5, d * cos(phi), d * sin(phi)]
+                v = zeros(3)
+                velocity_3D_bound_vortex!(
+                    v, filament, p, gamma,
+                    core_radius_fraction, work_vectors)
+
+                r_radial = [0.0, p[2], p[3]]
+                @test isapprox(dot(v, r0), 0.0; atol=1e-10)
+                @test isapprox(dot(v, r_radial), 0.0; atol=1e-8)
+            end
+        end
+    end
+
+    @testset "Solid-body rotation inside core" begin
+        filament = create_test_filament()
+
+        v_half = zeros(3)
+        v_quarter = zeros(3)
+        velocity_3D_bound_vortex!(
+            v_half, filament,
+            [0.5, 0.5 * core_radius_fraction, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+        velocity_3D_bound_vortex!(
+            v_quarter, filament,
+            [0.5, 0.25 * core_radius_fraction, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+
+        @test isapprox(norm(v_quarter), 0.5 * norm(v_half); rtol=1e-3)
+        @test isapprox(normalize(v_quarter), normalize(v_half); atol=1e-8)
+    end
+
+    @testset "Branch 1 / Branch 3 agree at core boundary" begin
+        filament = create_test_filament()
+        delta = 0.05 * core_radius_fraction
+        v_inside = zeros(3)
+        v_outside = zeros(3)
+        velocity_3D_bound_vortex!(
+            v_inside, filament,
+            [0.5, core_radius_fraction - delta, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+        velocity_3D_bound_vortex!(
+            v_outside, filament,
+            [0.5, core_radius_fraction + delta, 0.0],
+            gamma, core_radius_fraction, work_vectors)
+
+        @test isapprox(norm(v_inside), norm(v_outside); rtol=0.1)
+        @test isapprox(normalize(v_inside), normalize(v_outside); atol=1e-2)
+    end
+
+    @testset "Matches Rankine vortex profile inside core" begin
+        L = 1e6
+        r_c = 0.01 * L
+        filament_long = BoundFilament{Float64}()
+        reinit!(filament_long, [0.0, 0.0, 0.0], [L, 0.0, 0.0])
+
+        for r in (1.0, 10.0, 100.0, 1000.0)
+            v = zeros(3)
+            velocity_3D_bound_vortex!(
+                v, filament_long, [L/2, r, 0.0],
+                gamma, 0.01, work_vectors)
+
+            expected_mag = gamma * r / (2π * r_c^2)
+
+            @test isapprox(v[1], 0.0; atol=expected_mag * 1e-4)
+            @test isapprox(v[2], 0.0; atol=expected_mag * 1e-4)
+            @test isapprox(abs(v[3]), expected_mag; rtol=1e-3)
+            @test v[3] > 0
+        end
+    end
 end

test/filament/test_semi_infinite_filament.jl

@@ -118,7 +118,7 @@ end
         filament = create_test_filament2()
         vel_pos = zeros(3)
         vel_neg = zeros(3)
-        
+
         velocity_3D_trailing_vortex_semiinfinite!(
             vel_pos,
             filament,
@@ -137,7 +137,42 @@ end
             filament.vel_mag,
             work_vectors
         )
-        
+
         @test isapprox(vel_pos, -vel_neg)
     end
+
+    @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin
+        filament = create_test_filament2()
+        direction = filament.direction
+
+        for d in (1e-4, 1e-3, 5e-3, 1e-2, 1e-1)
+            for phi in (0.0, π/4, π/2, π, -π/3)
+                p = [0.5, d * cos(phi), d * sin(phi)]
+                v = zeros(3)
+                velocity_3D_trailing_vortex_semiinfinite!(
+                    v, filament, direction, p, gamma,
+                    filament.vel_mag, work_vectors)
+
+                r_radial = [0.0, p[2], p[3]]
+                @test isapprox(dot(v, direction), 0.0; atol=1e-10)
+                @test isapprox(dot(v, r_radial) / norm(v), 0.0; atol=1e-6)
+            end
+        end
+    end
+
+    @testset "Constant azimuthal direction inside core" begin
+        filament = create_test_filament2()
+        v_a = filament.vel_mag
+
+        d_inside = 1e-4
+        v1 = zeros(3); v2 = zeros(3)
+        velocity_3D_trailing_vortex_semiinfinite!(
+            v1, filament, filament.direction,
+            [0.5, d_inside, 0.0], gamma, v_a, work_vectors)
+        velocity_3D_trailing_vortex_semiinfinite!(
+            v2, filament, filament.direction,
+            [0.5, 2 * d_inside, 0.0], gamma, v_a, work_vectors)
+
+        @test isapprox(normalize(v2), normalize(v1); atol=1e-8)
+    end
 end