Skip to content

Use smooth norm and sqrt #240

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
Open

Use smooth norm and sqrt #240

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
98169f1
Use smooth norm and sqrt
1-Bart-1 May 23, 2026
62891db
Tighter smooth norm
1-Bart-1 May 23, 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
20 changes: 10 additions & 10 deletions src/body_aerodynamics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -625,7 +625,7 @@ function find_center_of_pressure(
normal[1] = v1y*v2z - v1z*v2y
normal[2] = v1z*v2x - v1x*v2z
normal[3] = v1x*v2y - v1y*v2x
normal_norm = norm3(normal)
normal_norm = smooth_norm3(normal)
if normal_norm != 0
normal[1] /= normal_norm
normal[2] /= normal_norm
Expand Down Expand Up @@ -771,8 +771,8 @@ function calculate_results(
panel, alpha_dist[i])
panel_width_array[i] = panel.width
va_norm = va_norm_array[i]
x_norm = norm3(panel.x_airf)
z_norm = norm3(panel.z_airf)
x_norm = smooth_norm3(panel.x_airf)
z_norm = smooth_norm3(panel.z_airf)
if va_norm == 0.0 || x_norm == 0.0 || z_norm == 0.0
alpha_geometric[i] = NaN
else
Expand Down Expand Up @@ -824,7 +824,7 @@ function calculate_results(
total_area > 0.0 || throw(ArgumentError(
"Total panel area must be positive."))
reference_speed = sqrt(weighted_speed_sq / total_area)
direction_norm = norm3(va_ref_vector)
direction_norm = smooth_norm3(va_ref_vector)
if direction_norm <= 0.0
va_ref_vector .= (1.0, 0.0, 0.0)
direction_norm = 1.0
Expand All @@ -833,7 +833,7 @@ function calculate_results(
va_ref_vector[k] = va_ref_vector[k] / direction_norm *
reference_speed
end
va_ref_mag = norm3(va_ref_vector)
va_ref_mag = smooth_norm3(va_ref_vector)
va_ref_mag > 0.0 || throw(ArgumentError(
"Reference freestream magnitude must be positive."))
va_ref_unit = body_aero.work_vectors[1]
Expand All @@ -843,7 +843,7 @@ function calculate_results(
end
dir_lift_ref = body_aero.work_vectors[2]
cross3!(dir_lift_ref, va_ref_vector, spanwise_direction)
dir_lift_ref_norm = norm3(dir_lift_ref)
dir_lift_ref_norm = smooth_norm3(dir_lift_ref)
dir_lift_ref_norm > 0.0 || throw(ArgumentError(
"Reference lift direction is undefined because " *
"reference flow is parallel to spanwise direction."))
Expand Down Expand Up @@ -874,14 +874,14 @@ function calculate_results(
induced_va_airfoil[k] = c_alpha * panel.x_airf[k] +
s_alpha * panel.z_airf[k]
end
normalize3!(induced_va_airfoil)
smooth_normalize3!(induced_va_airfoil)

cross3!(dir_lift_induced_va,
induced_va_airfoil, panel.y_airf)
normalize3!(dir_lift_induced_va)
smooth_normalize3!(dir_lift_induced_va)
cross3!(dir_drag_induced_va,
spanwise_direction, dir_lift_induced_va)
normalize3!(dir_drag_induced_va)
smooth_normalize3!(dir_drag_induced_va)

q_lift = 0.5 * density * v_a_dist[i]^2
lift_i = cl_array[i] * q_lift * chord_array[i]
Expand All @@ -900,7 +900,7 @@ function calculate_results(
q_panel = 0.5 * density * va_panel_mag^2
cross3!(dir_lift_prescribed_va,
panel.va, spanwise_direction)
normalize3!(dir_lift_prescribed_va)
smooth_normalize3!(dir_lift_prescribed_va)

lift_prescribed_va =
dot3(lift_induced_va, dir_lift_prescribed_va) +
Expand Down
64 changes: 30 additions & 34 deletions src/filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,35 +59,33 @@ function velocity_3D_bound_vortex!(
r2 .= XVP .- filament.x2

# Cut-off radius
nr0 = norm3(r0)
nr0 = smooth_norm3(r0)
epsilon = core_radius_fraction * nr0

cross3!(r1Xr2, r1, r2)
cross3!(r1Xr0, r1, r0)
nr1 = norm3(r1)
nr2 = norm3(r2)
nr1 = smooth_norm3(r1)
nr2 = smooth_norm3(r2)
@inbounds for k in 1:3
r1r2norm[k] = r1[k]/nr1 - r2[k]/nr2
end

# Check point location relative to filament
nr1Xr0 = norm3(r1Xr0)
nr1Xr0 = smooth_norm3(r1Xr0)
if nr1Xr0 / nr0 > epsilon
nr1Xr2 = norm3(r1Xr2)
nr1Xr2 = smooth_norm3(r1Xr2)
coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, r1r2norm)
@inbounds for k in 1:3
vel[k] = coeff * r1Xr2[k]
end
elseif nr1Xr0 / nr0 < 1e-12 * epsilon
vel .= 0.0
else
@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)
nr2Xr0 = smooth_norm3(r2Xr0)
d_r1_r0 = dot3(r1, r0)
d_r2_r0 = dot3(r2, r0)
@inbounds for k in 1:3
Expand All @@ -98,9 +96,9 @@ function velocity_3D_bound_vortex!(
end
cross3!(r1_projXr2_proj, r1_proj, r2_proj)

nr1pXr2p = norm3(r1_projXr2_proj)
nr1_proj = norm3(r1_proj)
nr2_proj = norm3(r2_proj)
nr1pXr2p = smooth_norm3(r1_projXr2_proj)
nr1_proj = smooth_norm3(r1_proj)
nr2_proj = smooth_norm3(r2_proj)
d_sum = 0.0
@inbounds for k in 1:3
d_sum += r0[k] * (r1_proj[k]/nr1_proj -
Expand Down Expand Up @@ -156,30 +154,30 @@ as implemented in KiteAeroDyn".
r2 .= XVP .- filament.x2

# Vector perpendicular to core radius
nr0 = norm3(r0)
nr0 = smooth_norm3(r0)
nr0sq = nr0 * nr0
d_r1_r0 = dot3(r1, r0)
@inbounds for k in 1:3
r_perp[k] = d_r1_r0 * r0[k] / nr0sq
end

# Cut-off radius
epsilon = sqrt(4 * ALPHA0 * NU * norm3(r_perp) / v_a)
epsilon = sqrt(4 * ALPHA0 * NU * smooth_norm3(r_perp) / v_a)

cross3!(r1Xr2, r1, r2)
cross3!(r1Xr0, r1, r0)
cross3!(r2Xr0, r2, r0)

nr1 = norm3(r1)
nr2 = norm3(r2)
nr1 = smooth_norm3(r1)
nr2 = smooth_norm3(r2)
@inbounds for k in 1:3
normr1r2[k] = r1[k]/nr1 - r2[k]/nr2
end

# Check point location relative to filament
nr1Xr0 = norm3(r1Xr0)
nr1Xr0 = smooth_norm3(r1Xr0)
if nr1Xr0 / nr0 > epsilon
nr1Xr2 = norm3(r1Xr2)
nr1Xr2 = smooth_norm3(r1Xr2)
coeff = (gamma / (4π)) / (nr1Xr2^2) * dot3(r0, normr1r2)
@inbounds for k in 1:3
vel[k] = coeff * r1Xr2[k]
Expand All @@ -190,7 +188,7 @@ as implemented in KiteAeroDyn".
# Project onto core radius — reuse r_perp, normr1r2
r1_proj = r_perp
r2_proj = normr1r2
nr2Xr0 = norm3(r2Xr0)
nr2Xr0 = smooth_norm3(r2Xr0)
d_r2_r0 = dot3(r2, r0)
@inbounds for k in 1:3
r1_proj[k] = d_r1_r0 * r0[k] / nr0sq +
Expand All @@ -200,9 +198,9 @@ as implemented in KiteAeroDyn".
end

cross3!(r1Xr2, r1_proj, r2_proj)
nr1Xr2_val = norm3(r1Xr2)
nr1_proj = norm3(r1_proj)
nr2_proj = norm3(r2_proj)
nr1Xr2_val = smooth_norm3(r1Xr2)
nr1_proj = smooth_norm3(r1_proj)
nr2_proj = smooth_norm3(r2_proj)
d_sum = 0.0
@inbounds for k in 1:3
d_sum += r0[k] * (r1_proj[k]/nr1_proj -
Expand Down Expand Up @@ -272,35 +270,33 @@ function velocity_3D_trailing_vortex_semiinfinite!(
@inbounds for k in 1:3
r_perp[k] = d_r1_Vf * Vf[k]
end
epsilon = sqrt(4 * ALPHA0 * NU * norm3(r_perp) / v_a)
epsilon = sqrt(4 * ALPHA0 * NU * smooth_norm3(r_perp) / v_a)

cross3!(r1XVf, r1, Vf)

nr1XVf = norm3(r1XVf)
nVf = norm3(Vf)
nr1 = norm3(r1)
nr1XVf = smooth_norm3(r1XVf)
nVf = smooth_norm3(Vf)
nr1 = smooth_norm3(r1)
if nr1XVf / nVf > epsilon
K = GAMMA / (4π) / (nr1XVf^2) * (1 + d_r1_Vf / nr1)
@inbounds for k in 1:3
vel[k] = K * r1XVf[k]
end
elseif nr1XVf / nVf < 1e-12 * epsilon
vel .= 0.0
else
r1_proj = work_vectors[4]
cross_tmp = work_vectors[5]
# temp = r1/nr1 - Vf
@inbounds for k in 1:3
cross_tmp[k] = r1[k]/nr1 - Vf[k]
end
n_tmp = norm3(cross_tmp)
n_tmp = smooth_norm3(cross_tmp)
@inbounds for k in 1:3
r1_proj[k] = d_r1_Vf * Vf[k] +
epsilon * cross_tmp[k] / n_tmp
end
cross3!(cross_tmp, r1_proj, Vf)
K = GAMMA / (4π) / (norm3(cross_tmp)^2) *
(1 + dot3(r1_proj, Vf) / norm3(r1_proj))
K = GAMMA / (4π) / (smooth_norm3(cross_tmp)^2) *
(1 + dot3(r1_proj, Vf) / smooth_norm3(r1_proj))
@inbounds for k in 1:3
vel[k] = K * cross_tmp[k]
end
Expand All @@ -324,10 +320,10 @@ Compute cross product of 3D vectors in-place.
nothing
end

@inline norm3(a) = sqrt(a[1]*a[1] + a[2]*a[2] + a[3]*a[3])
@inline smooth_norm3(a) = sqrt(a[1]*a[1] + a[2]*a[2] + a[3]*a[3] + 1e-18)
@inline dot3(a, b) = a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
@inline function normalize3!(v)
n = norm3(v)
n > 0 && (v[1] /= n; v[2] /= n; v[3] /= n)
@inline function smooth_normalize3!(v)
n = smooth_norm3(v)
v[1] /= n; v[2] /= n; v[3] /= n
nothing
end
2 changes: 1 addition & 1 deletion src/panel.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -596,7 +596,7 @@ function calculate_velocity_induced_bound_2D!(
cross3!(cross_, r0, r3)

# Calculate induced velocity
coeff = norm3(r0) / (2π * norm3(cross_)^2)
coeff = smooth_norm3(r0) / (2π * smooth_norm3(cross_)^2)
@inbounds for k in 1:3
U_2D[k] = cross_[k] * coeff
end
Expand Down
8 changes: 4 additions & 4 deletions src/solver.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -354,13 +354,13 @@ function solve!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_dist
dir_iva[k] = c_alpha * panel.x_airf[k] +
s_alpha * panel.z_airf[k]
end
normalize3!(dir_iva)
smooth_normalize3!(dir_iva)

# Calculate lift and drag directions
cross3!(dir_lift, dir_iva, panel.y_airf)
normalize3!(dir_lift)
smooth_normalize3!(dir_lift)
cross3!(dir_drag, spanwise_direction, dir_lift)
normalize3!(dir_drag)
smooth_normalize3!(dir_drag)

# Calculate force vectors
li = lift[i]
Expand Down Expand Up @@ -656,7 +656,7 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma
nothing
end

@inline smooth_sqrt(x) = sqrt(x + 1e-30)
@inline smooth_sqrt(x) = sqrt(x + 1e-18)

@inline function update_gamma_candidate!(
gamma_out,
Expand Down
14 changes: 7 additions & 7 deletions test/filament/test_bound_filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
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 Down
7 changes: 5 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 @@ -73,7 +73,10 @@ end
]
induced_velocity = zeros(3)

# Filament start point is singular in the current implementation.
# Singular geometry: smooth_norm3 floors the cross-product
# magnitude, so no division-by-zero produces NaN. Result must
# be finite; structural zero cross-products on the axis give
# zero velocity.
velocity_3D_trailing_vortex_semiinfinite!(
induced_velocity,
filament,
Expand All @@ -83,7 +86,7 @@ end
filament.vel_mag,
work_vectors
)
@test all(isnan.(induced_velocity))
@test all(isfinite.(induced_velocity))

for point in test_points
velocity_3D_trailing_vortex_semiinfinite!(
Expand Down
Loading