Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
44 changes: 22 additions & 22 deletions mwes/mwe_01.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# Replace va_norm_array = norm.(eachrow(solver.sol.va_array)) with a for loop
# Replace va_norm_dist = norm.(eachrow(solver.sol._va_dist)) with a for loop

# Testcase that shows that the new function is equivalent to the old, allocating line of code.
using Test
Expand All @@ -9,31 +9,31 @@ struct MockSolver
sol::NamedTuple
end

function calc_norm_array!(va_norm_array, va_array)
for i in 1:size(va_array, 1)
va_norm_array[i] = norm(view(va_array, i, :))
function calc_norm_dist!(va_norm_dist, va_dist)
for i in 1:size(va_dist, 1)
va_norm_dist[i] = norm(view(va_dist, i, :))
end
end

@testset "va_norm_array calculation" begin
global va_norm_array
@testset "va_norm_dist calculation" begin
global va_norm_dist

# Create a sample 2D array
sample_va_array = [
sample_va_dist = [
1.0 2.0 3.0;
4.0 5.0 6.0;
7.0 8.0 9.0
]

# Create a mock solver with the sample array
mock_solver = MockSolver((va_array = sample_va_array,))
mock_solver = MockSolver((va_dist = sample_va_dist,))

# Calculate va_norm_array
n = @allocated va_norm_array = norm.(eachrow(mock_solver.sol.va_array))
# Calculate va_norm_dist
n = @allocated va_norm_dist = norm.(eachrow(mock_solver.sol.va_dist))
println(n)

va_norm_array2 = zeros(3)
m = @allocated calc_norm_array!(va_norm_array2, sample_va_array)
va_norm_dist2 = zeros(3)
m = @allocated calc_norm_dist!(va_norm_dist2, sample_va_dist)
println(m)

# Expected results (calculated manually)
Expand All @@ -44,20 +44,20 @@ end
]

# Test the results
@test length(va_norm_array) == size(sample_va_array, 1)
@test va_norm_array ≈ expected_norms atol=1e-10
@test length(va_norm_dist) == size(sample_va_dist, 1)
@test va_norm_dist ≈ expected_norms atol=1e-10

# Test individual values
@test va_norm_array[1] ≈ norm(sample_va_array[1, :]) atol=1e-10
@test va_norm_array[2] ≈ norm(sample_va_array[2, :]) atol=1e-10
@test va_norm_array[3] ≈ norm(sample_va_array[3, :]) atol=1e-10
@test va_norm_dist[1] ≈ norm(sample_va_dist[1, :]) atol=1e-10
@test va_norm_dist[2] ≈ norm(sample_va_dist[2, :]) atol=1e-10
@test va_norm_dist[3] ≈ norm(sample_va_dist[3, :]) atol=1e-10

@test length(va_norm_array2) == size(sample_va_array, 1)
@test va_norm_array2 ≈ expected_norms atol=1e-10
@test length(va_norm_dist2) == size(sample_va_dist, 1)
@test va_norm_dist2 ≈ expected_norms atol=1e-10

# Test individual values
@test va_norm_array2[1] ≈ norm(sample_va_array[1, :]) atol=1e-10
@test va_norm_array2[2] ≈ norm(sample_va_array[2, :]) atol=1e-10
@test va_norm_array2[3] ≈ norm(sample_va_array[3, :]) atol=1e-10
@test va_norm_dist2[1] ≈ norm(sample_va_dist[1, :]) atol=1e-10
@test va_norm_dist2[2] ≈ norm(sample_va_dist[2, :]) atol=1e-10
@test va_norm_dist2[3] ≈ norm(sample_va_dist[3, :]) atol=1e-10
end
nothing
12 changes: 6 additions & 6 deletions mwes/mwe_warntype.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ using VortexStepMethod: calculate_AIC_matrices!, gamma_loop!, calculate_results,
velocity_3D_bound_vortex!,
velocity_3D_trailing_vortex!,
velocity_3D_trailing_vortex_semiinfinite!,
cross3!, calc_norm_array!,
cross3!, calc_norm_dist!,
Panel, reinit!, solve_base!
using LinearAlgebra
using StaticArrays
Expand All @@ -29,8 +29,8 @@ vel_app = [cos(alpha), 0.0, sin(alpha)] .* 20.0
set_va!(body_aero, vel_app)
solver = Solver(body_aero)

va_norm_array = ones(n_panels)
va_unit_array = ones(n_panels, 3)
va_norm_dist = ones(n_panels)
va_unit_dist = ones(n_panels, 3)

# Prepare args for individual functions
panel = body_aero.panels[1]
Expand Down Expand Up @@ -73,7 +73,7 @@ printstyled("\n$sep\n calculate_velocity_induced_bound_2D!\n$sep\n"; color=:cyan

printstyled("\n$sep\n calculate_AIC_matrices!\n$sep\n"; color=:cyan)
@code_warntype calculate_AIC_matrices!(
body_aero, VSM, 0.001, va_norm_array, va_unit_array)
body_aero, VSM, 0.001, va_norm_dist, va_unit_dist)

printstyled("\n$sep\n gamma_loop!\n$sep\n"; color=:cyan)
@code_warntype gamma_loop!(
Expand All @@ -82,5 +82,5 @@ printstyled("\n$sep\n gamma_loop!\n$sep\n"; color=:cyan)
printstyled("\n$sep\n solve_base!\n$sep\n"; color=:cyan)
@code_warntype solve_base!(solver, body_aero, nothing)

printstyled("\n$sep\n calc_norm_array!\n$sep\n"; color=:cyan)
@code_warntype calc_norm_array!(solver.br.va_norm_dist, solver.sol._va_dist)
printstyled("\n$sep\n calc_norm_dist!\n$sep\n"; color=:cyan)
@code_warntype calc_norm_dist!(solver.br.va_norm_dist, solver.sol._va_dist)
124 changes: 57 additions & 67 deletions src/body_aerodynamics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -384,10 +384,8 @@ end
end

"""
calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model,
core_radius_fraction,
va_norm_array,
va_unit_array)
calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::Model, core_radius_fraction,
va_norm_dist, va_unit_dist, target=body_aero.AIC)

Calculate Aerodynamic Influence Coefficient matrices.

Expand All @@ -397,8 +395,8 @@ Returns: nothing
"""
@inline function calculate_AIC_matrices!(body_aero::BodyAerodynamics{P, W, T}, model::Model,
core_radius_fraction,
va_norm_array::AbstractVector{T},
va_unit_array::AbstractMatrix{T},
va_norm_dist::AbstractVector{T},
va_unit_dist::AbstractMatrix{T},
target::AbstractArray{T, 3}=body_aero.AIC) where {P, W, T}
# Determine evaluation point based on model
evaluation_point = model == VSM ? :control_point : :aero_center
Expand All @@ -414,7 +412,7 @@ Returns: nothing
panel_areas = [panel.chord * panel.width for panel in body_aero.panels]
va_distribution = zeros(T, length(body_aero.panels), 3)
@inbounds for i in 1:length(body_aero.panels), k in 1:3
va_distribution[i, k] = va_unit_array[i, k] * va_norm_array[i]
va_distribution[i, k] = va_unit_dist[i, k] * va_norm_dist[i]
end
wake_velocity = _compute_reference_velocity_from_distribution(
va_distribution,
Expand Down Expand Up @@ -488,13 +486,9 @@ function calculate_circulation_distribution_elliptical_wing(gamma_i, body_aero::
end

"""
update_effective_angle_of_attack_if_VSM(body_aero::BodyAerodynamics, gamma,
core_radius_fraction,
z_airf_array,
x_airf_array,
va_array,
va_norm_array,
va_unit_array)
update_effective_angle_of_attack!(alpha_corrected, body_aero::BodyAerodynamics, gamma,
core_radius_fraction, z_airf_dist, x_airf_dist,
va_dist, va_norm_dist, va_unit_dist)

Update angle of attack at aerodynamic center for VSM method.

Expand All @@ -505,25 +499,25 @@ function update_effective_angle_of_attack!(alpha_corrected,
body_aero::BodyAerodynamics,
gamma,
core_radius_fraction,
z_airf_array,
x_airf_array,
va_array,
va_norm_array,
va_unit_array)
z_airf_dist,
x_airf_dist,
va_dist,
va_norm_dist,
va_unit_dist)

# Its own buffer: `AIC` holds the control-point matrix the circulation was solved
# against, so overwriting it here would leave post-solve readers on the LLT one.
calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_norm_array,
va_unit_array, body_aero.AIC_aero_center)
calculate_AIC_matrices!(body_aero, LLT, core_radius_fraction, va_norm_dist,
va_unit_dist, body_aero.AIC_aero_center)

induced_velocity = body_aero.cache[1][va_array]
induced_velocity = body_aero.cache[1][va_dist]
for k in 1:3
mul!(view(induced_velocity, :, k), view(body_aero.AIC_aero_center, :, :, k), gamma)
end

# In-place relative velocity calculation
relative_velocity = body_aero.cache[2][va_array]
relative_velocity .= va_array .+ induced_velocity
relative_velocity = body_aero.cache[2][va_dist]
relative_velocity .= va_dist .+ induced_velocity

# Preallocate and compute dot products manually
n = size(relative_velocity, 1)
Expand All @@ -534,8 +528,8 @@ function update_effective_angle_of_attack!(alpha_corrected,
vn = 0.0
vt = 0.0
for j in 1:3
vn += z_airf_array[i, j] * relative_velocity[i, j]
vt += x_airf_array[i, j] * relative_velocity[i, j]
vn += z_airf_dist[i, j] * relative_velocity[i, j]
vt += x_airf_dist[i, j] * relative_velocity[i, j]
end
v_normal[i] = vn
v_tangential[i] = vt
Expand Down Expand Up @@ -614,13 +608,13 @@ end

function find_center_of_pressure(
body_aero::BodyAerodynamics,
force_array,
moment_array,
force,
moment,
reference_point;
force_tol::Float64 = 1e-12
)
F = force_array
M0 = moment_array
F = force

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

MINOR: The parameters are renamed to force/moment, then aliased straight back to F/M0 (and r0), so the new names are never used past the first line. Either use them in the body or name the parameters F/M0; keeping both leaves two names for the same value.

M0 = moment
r0 = reference_point
F_norm_sq = dot3(F, F)
# Treat near-zero forces as "CoP undefined"
Expand Down Expand Up @@ -750,15 +744,11 @@ function set_pitch_rate_dist!(body_aero::BodyAerodynamics, omega)
end

"""
calculate_results(body_aero::BodyAerodynamics, gamma_new,
density,
core_radius_fraction, mu,
alpha_dist, v_a_dist,
chord_array, x_airf_array,
z_airf_array,
va_array, va_norm_array,
va_unit_array, panels::Vector{<:Panel},
is_only_f_and_gamma_output::Bool)
calculate_results(body_aero::BodyAerodynamics, gamma_new, reference_point, density,
core_radius_fraction, mu, alpha_dist, v_a_dist, chord_dist,
x_airf_dist, z_airf_dist, va_dist, va_norm_dist, va_unit_dist,
panels::Vector{<:Panel}, is_only_f_and_gamma_output::Bool;
correct_aoa=false, flow_curvature=false)

Calculate final aerodynamic results. Reference point is in the kite body (KB) frame.

Expand All @@ -777,12 +767,12 @@ function calculate_results(
mu,
alpha_dist,
v_a_dist,
chord_array,
x_airf_array,
z_airf_array,
va_array,
va_norm_array,
va_unit_array,
chord_dist,
x_airf_dist,
z_airf_dist,
va_dist,
va_norm_dist,
va_unit_dist,
panels::Vector{<:Panel},
is_only_f_and_gamma_output::Bool;
correct_aoa::Bool=false,
Expand All @@ -794,10 +784,10 @@ function calculate_results(
append!(body_aero.cache, [LazyBufferCache() for _ in 1:(15 - length(body_aero.cache))])
end

cl_array = body_aero.cache[5][alpha_dist]
cd_array = body_aero.cache[6][alpha_dist]
cm_array = body_aero.cache[7][alpha_dist]
panel_width_array = body_aero.cache[8][alpha_dist]
cl_dist = body_aero.cache[5][alpha_dist]
cd_dist = body_aero.cache[6][alpha_dist]
cm_dist = body_aero.cache[7][alpha_dist]
panel_width_dist = body_aero.cache[8][alpha_dist]
alpha_corrected = body_aero.cache[9][alpha_dist]
cl_prescribed_va = body_aero.cache[10][alpha_dist]
cd_prescribed_va = body_aero.cache[11][alpha_dist]
Expand All @@ -811,15 +801,15 @@ function calculate_results(

# Calculate coefficients and geometric AoA for each panel
for (i, panel) in enumerate(panels)
cl_array[i] = calculate_cl(panel, alpha_dist[i])
cd_array[i], cm_array[i] = calculate_cd_cm(
cl_dist[i] = calculate_cl(panel, alpha_dist[i])
cd_dist[i], cm_dist[i] = calculate_cd_cm(
panel, alpha_dist[i])
if flow_curvature
cm_array[i] += flow_curvature_cm(
body_aero.pitch_rate_dist[i], chord_array[i], v_a_dist[i])
cm_dist[i] += flow_curvature_cm(
body_aero.pitch_rate_dist[i], chord_dist[i], v_a_dist[i])
end
panel_width_array[i] = panel.width
va_norm = va_norm_array[i]
panel_width_dist[i] = panel.width
va_norm = va_norm_dist[i]
x_norm = norm3(panel.x_airf)
z_norm = norm3(panel.z_airf)
if va_norm == 0.0 || x_norm == 0.0 || z_norm == 0.0
Expand All @@ -841,11 +831,11 @@ function calculate_results(
body_aero,
gamma_new,
core_radius_fraction,
z_airf_array,
x_airf_array,
va_array,
va_norm_array,
va_unit_array
z_airf_dist,
x_airf_dist,
va_dist,
va_norm_dist,
va_unit_dist
)
else
alpha_corrected .= alpha_dist
Expand All @@ -862,13 +852,13 @@ function calculate_results(
weighted_speed_sq = 0.0
total_area = 0.0
@inbounds for i in 1:n_panels
area_i = chord_array[i] * panel_width_array[i]
area_i = chord_dist[i] * panel_width_dist[i]
total_area += area_i
speed_i = va_norm_array[i]
speed_i = va_norm_dist[i]
weighted_speed_sq += area_i * speed_i^2
va_ref_vector[1] += area_i * va_array[i, 1]
va_ref_vector[2] += area_i * va_array[i, 2]
va_ref_vector[3] += area_i * va_array[i, 3]
va_ref_vector[1] += area_i * va_dist[i, 1]
va_ref_vector[2] += area_i * va_dist[i, 2]
va_ref_vector[3] += area_i * va_dist[i, 3]
end
total_area > 0.0 || throw(ArgumentError(
"Total panel area must be positive."))
Expand Down Expand Up @@ -918,14 +908,14 @@ function calculate_results(
dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit)
loads = panel_loads(axes, dirs,
dynamic_pressure(density, density, v_a_dist[i]),
cl_array[i], cd_array[i], cm_array[i])
cl_dist[i], cd_dist[i], cm_dist[i])
moment_i = loads.moment
@inbounds for k in 1:3
lift_induced_va[k] = loads.lift * dirs.dir_lift[k]
drag_induced_va[k] = loads.drag * dirs.dir_drag[k]
end

va_panel_mag = va_norm_array[i]
va_panel_mag = va_norm_dist[i]
va_panel_mag > 0.0 || throw(ArgumentError(
"Panel $i has non-positive apparent " *
"velocity magnitude."))
Expand Down
Loading
Loading