diff --git a/mwes/mwe_01.jl b/mwes/mwe_01.jl index 685696ff..cfc40ca5 100644 --- a/mwes/mwe_01.jl +++ b/mwes/mwe_01.jl @@ -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 @@ -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) @@ -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 \ No newline at end of file diff --git a/mwes/mwe_warntype.jl b/mwes/mwe_warntype.jl index a7b6401c..1eef5db2 100644 --- a/mwes/mwe_warntype.jl +++ b/mwes/mwe_warntype.jl @@ -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 @@ -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] @@ -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!( @@ -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) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 8ecd60ff..d04a966d 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -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. @@ -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 @@ -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, @@ -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. @@ -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) @@ -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 @@ -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 + M0 = moment r0 = reference_point F_norm_sq = dot3(F, F) # Treat near-zero forces as "CoP undefined" @@ -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. @@ -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, @@ -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] @@ -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 @@ -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 @@ -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.")) @@ -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.")) diff --git a/src/solver.jl b/src/solver.jl index 9a706020..c16d3aa4 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -602,11 +602,11 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n return results end -@inline @inbounds function calc_norm_array!(va_norm_dist, va_array) - for i in axes(va_array, 1) +@inline @inbounds function calc_norm_dist!(va_norm_dist, va_dist) + for i in axes(va_dist, 1) va_norm_dist[i] = sqrt( - va_array[i,1]^2 + va_array[i,2]^2 + - va_array[i,3]^2) + va_dist[i,1]^2 + va_dist[i,2]^2 + + va_dist[i,3]^2) end end @@ -651,7 +651,7 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma end # Calculate unit vectors - calc_norm_array!(solver.br.va_norm_dist, solver.sol._va_dist) + calc_norm_dist!(solver.br.va_norm_dist, solver.sol._va_dist) @inbounds for i in 1:n_panels inv_norm = 1.0 / solver.br.va_norm_dist[i] for k in 1:3 @@ -709,70 +709,70 @@ end velocity_view_x, velocity_view_y, velocity_view_z, - va_array, + va_dist, induced_velocity_all, - relative_velocity_array, - y_airf_array, + relative_velocity_dist, + y_airf_dist, relative_velocity_crossz, - v_acrossz_array, - z_airf_array, - x_airf_array, - v_normal_array, - v_tangential_array, - va_magw_array, + v_acrossz_dist, + z_airf_dist, + x_airf_dist, + v_normal_dist, + v_tangential_dist, + va_magw_dist, cl_dist, - chord_array, + chord_dist, ) mul!(velocity_view_x, AIC_x, gamma_in) mul!(velocity_view_y, AIC_y, gamma_in) mul!(velocity_view_z, AIC_z, gamma_in) - relative_velocity_array .= va_array .+ induced_velocity_all + relative_velocity_dist .= va_dist .+ induced_velocity_all @inbounds for i in 1:n_panels - ax = relative_velocity_array[i,1] - ay = relative_velocity_array[i,2] - az = relative_velocity_array[i,3] - bx = y_airf_array[i,1] - by = y_airf_array[i,2] - bz = y_airf_array[i,3] + ax = relative_velocity_dist[i,1] + ay = relative_velocity_dist[i,2] + az = relative_velocity_dist[i,3] + bx = y_airf_dist[i,1] + by = y_airf_dist[i,2] + bz = y_airf_dist[i,3] relative_velocity_crossz[i,1] = ay*bz - az*by relative_velocity_crossz[i,2] = az*bx - ax*bz relative_velocity_crossz[i,3] = ax*by - ay*bx - ax = va_array[i,1] - ay = va_array[i,2] - az = va_array[i,3] - v_acrossz_array[i,1] = ay*bz - az*by - v_acrossz_array[i,2] = az*bx - ax*bz - v_acrossz_array[i,3] = ax*by - ay*bx + ax = va_dist[i,1] + ay = va_dist[i,2] + az = va_dist[i,3] + v_acrossz_dist[i,1] = ay*bz - az*by + v_acrossz_dist[i,2] = az*bx - ax*bz + v_acrossz_dist[i,3] = ax*by - ay*bx end @inbounds for i in 1:n_panels - v_normal_array[i] = - z_airf_array[i,1]*relative_velocity_array[i,1] + - z_airf_array[i,2]*relative_velocity_array[i,2] + - z_airf_array[i,3]*relative_velocity_array[i,3] - v_tangential_array[i] = - x_airf_array[i,1]*relative_velocity_array[i,1] + - x_airf_array[i,2]*relative_velocity_array[i,2] + - x_airf_array[i,3]*relative_velocity_array[i,3] + v_normal_dist[i] = + z_airf_dist[i,1]*relative_velocity_dist[i,1] + + z_airf_dist[i,2]*relative_velocity_dist[i,2] + + z_airf_dist[i,3]*relative_velocity_dist[i,3] + v_tangential_dist[i] = + x_airf_dist[i,1]*relative_velocity_dist[i,1] + + x_airf_dist[i,2]*relative_velocity_dist[i,2] + + x_airf_dist[i,3]*relative_velocity_dist[i,3] end - solver.lr.alpha_dist .= atan.(v_normal_array, v_tangential_array) + solver.lr.alpha_dist .= atan.(v_normal_dist, v_tangential_dist) @inbounds for i in 1:n_panels solver.lr.v_a_dist[i] = smooth_sqrt( relative_velocity_crossz[i,1]^2 + relative_velocity_crossz[i,2]^2 + relative_velocity_crossz[i,3]^2) - va_magw_array[i] = smooth_sqrt( - v_acrossz_array[i,1]^2 + - v_acrossz_array[i,2]^2 + - v_acrossz_array[i,3]^2) + va_magw_dist[i] = smooth_sqrt( + v_acrossz_dist[i,1]^2 + + v_acrossz_dist[i,2]^2 + + v_acrossz_dist[i,3]^2) end for (i, (panel, alpha)) in enumerate(zip(panels, solver.lr.alpha_dist)) cl_dist[i] = calculate_cl(panel, alpha) end - gamma_out .= 0.5 .* solver.lr.v_a_dist.^2 ./ va_magw_array .* cl_dist .* chord_array + gamma_out .= 0.5 .* solver.lr.v_a_dist.^2 ./ va_magw_dist .* cl_dist .* chord_dist return nothing end @@ -822,7 +822,7 @@ end """ apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, viscosity_matrix, - lift_slope, mu_array, gamma_target, planform_area, factor) + lift_slope, mu_dist, gamma_target, planform_area, factor) Apply one implicit Li/Gaunaa artificial-viscosity step to `gamma` in place and return `true` when it fired. The per-panel viscosity is @@ -834,21 +834,21 @@ post-stall (`mu_i > 0`), `gamma` is replaced by the solution of remaining arguments are preallocated work buffers reused across iterations. """ function apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, viscosity_matrix, - lift_slope, mu_array, gamma_target, planform_area, factor) + lift_slope, mu_dist, gamma_target, planform_area, factor) n_panels = length(panels) local_lift_slope!(lift_slope, panels, alpha_dist) any_stalled = false @inbounds for i in 1:n_panels m = -factor * planform_area * lift_slope[i] / panels[i].width^2 - mu_array[i] = max(zero(eltype(mu_array)), m) - mu_array[i] > 0 && (any_stalled = true) + mu_dist[i] = max(zero(eltype(mu_dist)), m) + mu_dist[i] > 0 && (any_stalled = true) end any_stalled || return false gamma_target .= gamma one_t, zero_t = one(eltype(viscosity_matrix)), zero(eltype(viscosity_matrix)) @inbounds for col in 1:n_panels, row in 1:n_panels viscosity_matrix[row, col] = - (row == col ? one_t : zero_t) - mu_array[row] * laplacian[row, col] + (row == col ? one_t : zero_t) - mu_dist[row] * laplacian[row, col] end ldiv!(gamma, lu!(viscosity_matrix), gamma_target) return true @@ -880,28 +880,28 @@ function gamma_loop!( relaxation_factor; log::Bool = true ) where {P, U, T} - va_array = solver.sol._va_dist - chord_array = solver.sol._chord_dist - x_airf_array = solver.sol._x_airf_dist - y_airf_array = solver.sol._y_airf_dist - z_airf_array = solver.sol._z_airf_dist + va_dist = solver.sol._va_dist + chord_dist = solver.sol._chord_dist + x_airf_dist = solver.sol._x_airf_dist + y_airf_dist = solver.sol._y_airf_dist + z_airf_dist = solver.sol._z_airf_dist solver.lr.converged = false n_panels = length(body_aero.panels) solver.lr.alpha_dist .= body_aero.alpha_dist solver.lr.v_a_dist .= body_aero.v_a_dist - va_magw_array = solver.cache[1][solver.lr.v_a_dist] + va_magw_dist = solver.cache[1][solver.lr.v_a_dist] gamma = solver.cache[2][solver.lr.gamma_new] abs_gamma_new = solver.cache[3][solver.lr.gamma_new] - induced_velocity_all = solver.cache[4][va_array] - relative_velocity_array = solver.cache[5][va_array] - relative_velocity_crossz = solver.cache[6][va_array] - v_acrossz_array = solver.cache[7][va_array] - cl_dist = solver.cache[8][solver.lr.gamma_new] + induced_velocity_all = solver.cache[4][va_dist] + relative_velocity_dist = solver.cache[5][va_dist] + relative_velocity_crossz = solver.cache[6][va_dist] + v_acrossz_dist = solver.cache[7][va_dist] + cl_dist = solver.cache[8][solver.lr.gamma_new] damp = solver.cache[9][solver.lr.gamma_new] damp .= zero(T) - v_normal_array = solver.cache[10][solver.lr.gamma_new] - v_tangential_array = solver.cache[11][solver.lr.gamma_new] + v_normal_dist = solver.cache[10][solver.lr.gamma_new] + v_tangential_dist = solver.cache[11][solver.lr.gamma_new] AIC_x = @view body_aero.AIC[:, :, 1] AIC_y = @view body_aero.AIC[:, :, 2] @@ -928,11 +928,11 @@ function gamma_loop!( residual, gamma_iter, solver, panels, n_panels, AIC_x, AIC_y, AIC_z, velocity_view_x, velocity_view_y, velocity_view_z, - va_array, induced_velocity_all, relative_velocity_array, - y_airf_array, relative_velocity_crossz, v_acrossz_array, - z_airf_array, x_airf_array, - v_normal_array, v_tangential_array, - va_magw_array, cl_dist, chord_array, + va_dist, induced_velocity_all, relative_velocity_dist, + y_airf_dist, relative_velocity_crossz, v_acrossz_dist, + z_airf_dist, x_airf_dist, + v_normal_dist, v_tangential_dist, + va_magw_dist, cl_dist, chord_dist, ) @inbounds for i in 1:n_panels residual[i] -= gamma_iter[i] @@ -950,11 +950,11 @@ function gamma_loop!( residual_perturbed, gamma_perturbed, solver, panels, n_panels, AIC_x, AIC_y, AIC_z, velocity_view_x, velocity_view_y, velocity_view_z, - va_array, induced_velocity_all, relative_velocity_array, - y_airf_array, relative_velocity_crossz, v_acrossz_array, - z_airf_array, x_airf_array, - v_normal_array, v_tangential_array, - va_magw_array, cl_dist, chord_array, + va_dist, induced_velocity_all, relative_velocity_dist, + y_airf_dist, relative_velocity_crossz, v_acrossz_dist, + z_airf_dist, x_airf_dist, + v_normal_dist, v_tangential_dist, + va_magw_dist, cl_dist, chord_dist, ) inv_step = 1.0 / step for i in 1:n_panels @@ -979,11 +979,11 @@ function gamma_loop!( residual_perturbed, gamma_perturbed, solver, panels, n_panels, AIC_x, AIC_y, AIC_z, velocity_view_x, velocity_view_y, velocity_view_z, - va_array, induced_velocity_all, relative_velocity_array, - y_airf_array, relative_velocity_crossz, v_acrossz_array, - z_airf_array, x_airf_array, - v_normal_array, v_tangential_array, - va_magw_array, cl_dist, chord_array, + va_dist, induced_velocity_all, relative_velocity_dist, + y_airf_dist, relative_velocity_crossz, v_acrossz_dist, + z_airf_dist, x_airf_dist, + v_normal_dist, v_tangential_dist, + va_magw_dist, cl_dist, chord_dist, ) @inbounds for i in 1:n_panels residual_perturbed[i] -= gamma_perturbed[i] @@ -1016,13 +1016,13 @@ function gamma_loop!( laplacian = use_viscosity ? zeros(T, n_panels, n_panels) : zeros(T, 0, 0) viscosity_matrix = use_viscosity ? zeros(T, n_panels, n_panels) : zeros(T, 0, 0) lift_slope = use_viscosity ? zeros(T, n_panels) : zeros(T, 0) - mu_array = use_viscosity ? zeros(T, n_panels) : zeros(T, 0) + mu_dist = use_viscosity ? zeros(T, n_panels) : zeros(T, 0) gamma_target = use_viscosity ? zeros(T, n_panels) : zeros(T, 0) planform_area = zero(T) if use_viscosity build_spanwise_laplacian!(laplacian, n_panels) @inbounds for i in 1:n_panels - planform_area += panels[i].width * chord_array[i] + planform_area += panels[i].width * chord_dist[i] end end @@ -1040,25 +1040,25 @@ function gamma_loop!( velocity_view_x, velocity_view_y, velocity_view_z, - va_array, + va_dist, induced_velocity_all, - relative_velocity_array, - y_airf_array, + relative_velocity_dist, + y_airf_dist, relative_velocity_crossz, - v_acrossz_array, - z_airf_array, - x_airf_array, - v_normal_array, - v_tangential_array, - va_magw_array, + v_acrossz_dist, + z_airf_dist, + x_airf_dist, + v_normal_dist, + v_tangential_dist, + va_magw_dist, cl_dist, - chord_array, + chord_dist, ) # Gate the linear solve on any(mu > 0): fires exactly in post-stall. if use_viscosity apply_artificial_viscosity!( gamma_new, panels, solver.lr.alpha_dist, laplacian, - viscosity_matrix, lift_slope, mu_array, gamma_target, + viscosity_matrix, lift_slope, mu_dist, gamma_target, planform_area, solver.artificial_viscosity_factor, ) end diff --git a/test/bench.jl b/test/bench.jl index fb0a2c6b..c956add9 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -83,8 +83,8 @@ using LinearAlgebra AIC_z = similar(AIC_x) v_ind = zeros(3) point = rand(3) - 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) models = [VSM, LLT] core_radius_fractions = [0.001, 10.0] @@ -94,7 +94,7 @@ using LinearAlgebra for model in models for frac in core_radius_fractions @testset "Model $model Core Radius Fraction $frac" begin - result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_norm_array, $va_unit_array) samples=1 evals=1 + result = @benchmark calculate_AIC_matrices!($body_aero, $model, $frac, $va_norm_dist, $va_unit_dist) samples=1 evals=1 @test result.allocs ≤ 30 @info "Model: $(model) \t Core radius fraction: $(frac) \t Allocations: $(result.allocs) \t Memory: $(result.memory)" end @@ -106,19 +106,19 @@ using LinearAlgebra @info "Gamma Loop" # Pre-allocate arrays gamma_new = zeros(n_panels) - va_array = zeros(n_panels, 3) - chord_array = zeros(n_panels) - x_airf_array = zeros(n_panels, 3) - y_airf_array = zeros(n_panels, 3) - z_airf_array = zeros(n_panels, 3) + va_dist = zeros(n_panels, 3) + chord_dist = zeros(n_panels) + x_airf_dist = zeros(n_panels, 3) + y_airf_dist = zeros(n_panels, 3) + z_airf_dist = zeros(n_panels, 3) # Fill arrays with data for (i, panel) in enumerate(body_aero.panels) - va_array[i, :] .= panel.va - chord_array[i] = panel.chord - x_airf_array[i, :] .= panel.x_airf - y_airf_array[i, :] .= panel.y_airf - z_airf_array[i, :] .= panel.z_airf + va_dist[i, :] .= panel.va + chord_dist[i] = panel.chord + x_airf_dist[i, :] .= panel.x_airf + y_airf_dist[i, :] .= panel.y_airf + z_airf_dist[i, :] .= panel.z_airf end alphas = collect(-20.0:30.0) @@ -145,11 +145,11 @@ using LinearAlgebra solver = Solver(body_aero; aerodynamic_model_type=model ) - solver.sol._va_dist .= va_array - solver.sol._chord_dist .= chord_array - solver.sol._x_airf_dist .= x_airf_array - solver.sol._y_airf_dist .= y_airf_array - solver.sol._z_airf_dist .= z_airf_array + solver.sol._va_dist .= va_dist + solver.sol._chord_dist .= chord_dist + solver.sol._x_airf_dist .= x_airf_dist + solver.sol._y_airf_dist .= y_airf_dist + solver.sol._z_airf_dist .= z_airf_dist result = @benchmark gamma_loop!( $solver, $body_aero, @@ -165,29 +165,29 @@ using LinearAlgebra @testset "Results Calculation" begin # Pre-allocate arrays - alpha_array = zeros(n_panels) - v_a_array = zeros(n_panels) - chord_array = zeros(n_panels) - x_airf_array = zeros(n_panels, 3) - y_airf_array = zeros(n_panels, 3) - z_airf_array = zeros(n_panels, 3) - va_array = zeros(n_panels, 3) - va_norm_array = zeros(n_panels) - va_unit_array = zeros(n_panels, 3) + alpha_dist = zeros(n_panels) + v_a_dist = zeros(n_panels) + chord_dist = zeros(n_panels) + x_airf_dist = zeros(n_panels, 3) + y_airf_dist = zeros(n_panels, 3) + z_airf_dist = zeros(n_panels, 3) + va_dist = zeros(n_panels, 3) + va_norm_dist = zeros(n_panels) + va_unit_dist = zeros(n_panels, 3) reference_point = zeros(3) set_va!(body_aero, vel_app) # Fill arrays with panel data to satisfy calculate_results preconditions. for (i, panel) in enumerate(body_aero.panels) - chord_array[i] = panel.chord - x_airf_array[i, :] .= panel.x_airf - y_airf_array[i, :] .= panel.y_airf - z_airf_array[i, :] .= panel.z_airf - va_array[i, :] .= panel.va - va_norm_array[i] = norm(panel.va) - va_unit_array[i, :] .= va_norm_array[i] > 0.0 ? panel.va ./ va_norm_array[i] : [1.0, 0.0, 0.0] - v_a_array[i] = va_norm_array[i] + chord_dist[i] = panel.chord + x_airf_dist[i, :] .= panel.x_airf + y_airf_dist[i, :] .= panel.y_airf + z_airf_dist[i, :] .= panel.z_airf + va_dist[i, :] .= panel.va + va_norm_dist[i] = norm(panel.va) + va_unit_dist[i, :] .= va_norm_dist[i] > 0.0 ? panel.va ./ va_norm_dist[i] : [1.0, 0.0, 0.0] + v_a_dist[i] = va_norm_dist[i] end results = @MVector zeros(3) @@ -198,14 +198,14 @@ using LinearAlgebra $density, 1e-20, 0.0, - $alpha_array, - $v_a_array, - $chord_array, - $x_airf_array, - $z_airf_array, - $va_array, - $va_norm_array, - $va_unit_array, + $alpha_dist, + $v_a_dist, + $chord_dist, + $x_airf_dist, + $z_airf_dist, + $va_dist, + $va_norm_dist, + $va_unit_dist, $body_aero.panels, false ) samples=1 evals=1 diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index ab705b2b..581f73ce 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -63,14 +63,14 @@ end ) # Calculate new matrices - va_norm_array = fill(norm(v_a), length(coord)) - va_unit_array = repeat(reshape(v_a ./ norm(v_a), 1, 3), length(coord)) + va_norm_dist = fill(norm(v_a), length(body_aero.panels)) + va_unit_dist = repeat(reshape(v_a ./ norm(v_a), 1, 3), length(body_aero.panels)) calculate_AIC_matrices!( body_aero, LLT, core_radius_fraction, - va_norm_array, - va_unit_array + va_norm_dist, + va_unit_dist ) AIC_x, AIC_y, AIC_z = @views body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3] @@ -98,14 +98,14 @@ end ) # Calculate new matrices - va_norm_array = fill(norm(v_a), length(coord)) - va_unit_array = repeat(reshape(v_a ./ norm(v_a), 1, 3), length(coord)) + va_norm_dist = fill(norm(v_a), length(body_aero.panels)) + va_unit_dist = repeat(reshape(v_a ./ norm(v_a), 1, 3), length(body_aero.panels)) calculate_AIC_matrices!( body_aero, VSM, core_radius_fraction, - va_norm_array, - va_unit_array + va_norm_dist, + va_unit_dist ) AIC_x, AIC_y, AIC_z = body_aero.AIC[:, :, 1], body_aero.AIC[:, :, 2], body_aero.AIC[:, :, 3] diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index b940c69e..56896de1 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -211,7 +211,7 @@ roughness(v) = sum(abs, @views v[1:end-2] .- 2 .* v[2:end-1] .+ v[3:end]) VortexStepMethod.build_spanwise_laplacian!(laplacian, n) viscosity_matrix = zeros(n, n) lift_slope = zeros(n) - mu_array = zeros(n) + mu_dist = zeros(n) gamma_target = zeros(n) planform_area = sum(p.width * p.chord for p in panels) @@ -223,7 +223,7 @@ roughness(v) = sum(abs, @views v[1:end-2] .- 2 .* v[2:end-1] .+ v[3:end]) # the solve is skipped, and gamma is returned untouched. gamma_attached = spiky() fired_attached = VortexStepMethod.apply_artificial_viscosity!(gamma_attached, - panels, attached, laplacian, viscosity_matrix, lift_slope, mu_array, + panels, attached, laplacian, viscosity_matrix, lift_slope, mu_dist, gamma_target, planform_area, 0.035) @test !fired_attached @test gamma_attached == spiky() @@ -232,7 +232,7 @@ roughness(v) = sum(abs, @views v[1:end-2] .- 2 .* v[2:end-1] .+ v[3:end]) gamma_stalled = spiky() rough_before = roughness(gamma_stalled) fired_stalled = VortexStepMethod.apply_artificial_viscosity!(gamma_stalled, - panels, post_stall, laplacian, viscosity_matrix, lift_slope, mu_array, + panels, post_stall, laplacian, viscosity_matrix, lift_slope, mu_dist, gamma_target, planform_area, 0.035) @test fired_stalled @test roughness(gamma_stalled) < rough_before @@ -240,10 +240,10 @@ roughness(v) = sum(abs, @views v[1:end-2] .- 2 .* v[2:end-1] .+ v[3:end]) # The attached (hot) path must not allocate. gamma_alloc = spiky() VortexStepMethod.apply_artificial_viscosity!(gamma_alloc, panels, attached, - laplacian, viscosity_matrix, lift_slope, mu_array, gamma_target, + laplacian, viscosity_matrix, lift_slope, mu_dist, gamma_target, planform_area, 0.035) allocs = @allocated VortexStepMethod.apply_artificial_viscosity!(gamma_alloc, - panels, attached, laplacian, viscosity_matrix, lift_slope, mu_array, + panels, attached, laplacian, viscosity_matrix, lift_slope, mu_dist, gamma_target, planform_area, 0.035) @test allocs == 0 end