Skip to content

Fix branch 3 direction #241

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
1-Bart-1 wants to merge 2 commits into
base: main
Choose a base branch
from
+133 −20
Open

Fix branch 3 direction #241

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
047bf7a
Fix branch 3 direction
1-Bart-1 May 24, 2026
2c37cbd
Fix tests
1-Bart-1 May 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 10 additions & 8 deletions src/filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down
96 changes: 86 additions & 10 deletions test/filament/test_bound_filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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!(
Expand All @@ -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
39 changes: 37 additions & 2 deletions test/filament/test_semi_infinite_filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Loading