From 047bf7a9bfc2c4d422bb03f137dc0c2d8b2c135c Mon Sep 17 00:00:00 2001 From: Bart Date: Sun, 24 May 2026 23:47:50 +0200 Subject: [PATCH 1/2] Fix branch 3 direction --- src/filament.jl | 22 +++-- test/filament/test_bound_filament.jl | 92 ++++++++++++++++++-- test/filament/test_semi_infinite_filament.jl | 52 ++++++++++- 3 files changed, 149 insertions(+), 17 deletions(-) diff --git a/src/filament.jl b/src/filament.jl index e2c5a772..4801cb9c 100644 --- a/src/filament.jl +++ b/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) diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl index 511af4c3..c5022df4 100644 --- a/test/filament/test_bound_filament.jl +++ b/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 \ No newline at end of file diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl index ab870fb5..189c4d23 100644 --- a/test/filament/test_semi_infinite_filament.jl +++ b/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 From 2c37cbdfedd922202f700f14110796e2f03e9617 Mon Sep 17 00:00:00 2001 From: Bart Date: Tue, 26 May 2026 00:24:03 +0200 Subject: [PATCH 2/2] Fix tests --- src/filament.jl | 6 +-- test/filament/test_bound_filament.jl | 52 ++++++++++---------- test/filament/test_semi_infinite_filament.jl | 17 +------ 3 files changed, 28 insertions(+), 47 deletions(-) diff --git a/src/filament.jl b/src/filament.jl index 4801cb9c..c1a5d5d5 100644 --- a/src/filament.jl +++ b/src/filament.jl @@ -84,13 +84,10 @@ function velocity_3D_bound_vortex!( @debug "inside core radius" @debug "distance from control point to filament: $(nr1Xr0 / nr0)" - # 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 d_r1_r0 = dot3(r1, r0) d_r2_r0 = dot3(r2, r0) - r_rad = r1Xr0 # reuse buffer; no longer needed below + r_rad = r1Xr0 @inbounds for k in 1:3 r_rad[k] = r1[k] - d_r1_r0 * r0[k] / nr0sq end @@ -294,7 +291,6 @@ function velocity_3D_trailing_vortex_semiinfinite!( else r1_proj = work_vectors[4] cross_tmp = work_vectors[5] - # Radial direction: perpendicular component of r1 to Vf. nVfsq = nVf * nVf @inbounds for k in 1:3 cross_tmp[k] = r1[k] - d_r1_Vf * Vf[k] / nVfsq diff --git a/test/filament/test_bound_filament.jl b/test/filament/test_bound_filament.jl index c5022df4..11996c51 100644 --- a/test/filament/test_bound_filament.jl +++ b/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 @@ -176,17 +176,9 @@ end 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)] @@ -195,10 +187,7 @@ end 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 @@ -206,9 +195,6 @@ 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) @@ -222,20 +208,13 @@ end [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 + delta = 0.05 * core_radius_fraction v_inside = zeros(3) v_outside = zeros(3) velocity_3D_bound_vortex!( @@ -247,9 +226,28 @@ end [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 + + @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 \ No newline at end of file diff --git a/test/filament/test_semi_infinite_filament.jl b/test/filament/test_semi_infinite_filament.jl index 189c4d23..43b2149f 100644 --- a/test/filament/test_semi_infinite_filament.jl +++ b/test/filament/test_semi_infinite_filament.jl @@ -142,15 +142,9 @@ end 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)] @@ -161,22 +155,16 @@ end 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. + @testset "Constant azimuthal direction inside core" begin 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 + d_inside = 1e-4 v1 = zeros(3); v2 = zeros(3) velocity_3D_trailing_vortex_semiinfinite!( v1, filament, filament.direction, @@ -185,7 +173,6 @@ end 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