Fix branch 3 direction

Author: 1-Bart-1

src/filament.jl

@@ -84,17 +84,22 @@ 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)
+        # Project onto core cylinder along the radial direction
+        # (perpendicular component of r1 to r0). r1 and r2 share the
+        # same radial component, so one radial unit vector serves both.
         nr0sq = nr0 * nr0
-        nr2Xr0 = norm3(r2Xr0)
         d_r1_r0 = dot3(r1, r0)
         d_r2_r0 = dot3(r2, r0)
+        r_rad = r1Xr0  # reuse buffer; no longer needed below
+        @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 +294,14 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     else
         r1_proj = work_vectors[4]
         cross_tmp = work_vectors[5]
-        # temp = r1/nr1 - Vf
+        # Radial direction: perpendicular component of r1 to 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

@@ -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,82 @@ end
         )
         @test isapprox(velocities[2], -v_neg)
     end
+
+    @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin
+        # A real vortex's induced velocity is purely azimuthal:
+        # perpendicular to BOTH the filament axis r0 AND the radial
+        # direction from the axis to the field point. The old Branch 3
+        # projected the field point in the azimuthal direction instead
+        # of the radial direction, producing a velocity with a nonzero
+        # radial component. This test catches that.
+        filament = create_test_filament()
+        r0 = [1.0, 0.0, 0.0]
+
+        # Probe at multiple distances spanning both branches (in-core
+        # and outside-core) and azimuthal positions.
+        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)
+
+                # Radial direction at p (perpendicular to axis)
+                r_radial = [0.0, p[2], p[3]]
+
+                # v must be perpendicular to both axis and radial
+                @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
+        # Inside the core, velocity magnitude scales linearly with
+        # d_perp (Branch 3 = linear ramp from 0 on axis to peak at
+        # boundary). Direction stays constant along a radial line.
+        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)
+
+        # Magnitudes scale linearly with d_perp
+        @test isapprox(norm(v_quarter), 0.5 * norm(v_half); rtol=1e-3)
+        # Directions are parallel
+        @test isapprox(normalize(v_quarter), normalize(v_half); atol=1e-8)
+    end
+
+    @testset "Branch 1 / Branch 3 agree at core boundary" begin
+        # The two branches are designed to be continuous at d_perp =
+        # epsilon. Probe just inside and just outside the boundary with
+        # enough delta to clearly land in different branches; results
+        # must agree in both magnitude AND direction. Old code returned
+        # orthogonal directions here.
+        filament = create_test_filament()
+        delta = 0.05 * core_radius_fraction  # 5% in/out
+        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)
+
+        # Magnitudes within a few % of each other (peak is at boundary)
+        @test isapprox(norm(v_inside), norm(v_outside); rtol=0.1)
+        # Directions match — would fail with old azimuthal-projection bug
+        @test isapprox(normalize(v_inside), normalize(v_outside); atol=1e-2)
+    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,55 @@ end
             filament.vel_mag,
             work_vectors
         )
-        
+
         @test isapprox(vel_pos, -vel_neg)
     end
+
+    @testset "Velocity is azimuthal (perpendicular to axis and radius)" begin
+        # Vortex induced velocity is purely azimuthal: perpendicular
+        # to BOTH the filament axis and the radial direction. Old
+        # Branch 3 projected along the azimuthal direction instead of
+        # the radial direction, giving a non-azimuthal velocity.
+        filament = create_test_filament2()
+        direction = filament.direction
+
+        # Lamb-Oseen epsilon at x=0.5 is ~sqrt(4·α₀·ν·0.5) ≈ 6e-3.
+        # Probe across that range.
+        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)
+                # Looser tol — semi-infinite velocity has an axial
+                # contribution from the (1 + r1·Vf/|r1|) factor.
+                @test isapprox(dot(v, r_radial) / norm(v), 0.0; atol=1e-6)
+            end
+        end
+    end
+
+    @testset "Solid-body rotation inside core" begin
+        # At small d_perp (inside the Lamb-Oseen core), Branch 3 gives
+        # a linear ramp in magnitude with constant direction.
+        filament = create_test_filament2()
+        v_a = filament.vel_mag
+
+        # epsilon at x=0.5 ≈ sqrt(4·α₀·ν·0.5/v_a) ≈ 6.1e-3
+        # Probe well inside the core.
+        d_inside = 1e-4   # ~60× smaller than epsilon — definitely inside
+        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(norm(v2), 2 * norm(v1); rtol=1e-3)
+        @test isapprox(normalize(v2), normalize(v1); atol=1e-8)
+    end
 end