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 4 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
4 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
047bf7a
Fix branch 3 direction
1-Bart-1 May 24, 2026
4c312cf
Merge branch 'fix-projection-radial' into smooth
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
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
84 changes: 43 additions & 41 deletions src/filament.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,48 +59,51 @@ 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)
# 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 = smooth_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)

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 +159,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 +193,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 +203,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 +275,34 @@ 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
# 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)
n_tmp = smooth_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)
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 +326,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
Loading
Loading