From 23f172aecd036720fb7749459ce7ad945566a4e9 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Tue, 8 Apr 2025 15:27:08 +0200 Subject: [PATCH 1/4] Remove parallel compute --- scripts/polars.jl | 195 ---------------------- src/VortexStepMethod.jl | 5 +- src/polars.jl | 168 +++++++++++++++++++ src/{kite_geometry.jl => ram_geometry.jl} | 6 +- 4 files changed, 173 insertions(+), 201 deletions(-) delete mode 100644 scripts/polars.jl create mode 100644 src/polars.jl rename src/{kite_geometry.jl => ram_geometry.jl} (98%) diff --git a/scripts/polars.jl b/scripts/polars.jl deleted file mode 100644 index 98bf95ac..00000000 --- a/scripts/polars.jl +++ /dev/null @@ -1,195 +0,0 @@ -using Distributed, Timers, Serialization, SharedArrays, StaticArrays -using Interpolations -using Xfoil -using Logging - -const SPEED_OF_SOUND = 343 # [m/s] at 20 °C, see: https://en.wikipedia.org/wiki/Speed_of_sound -const KINEMATIC_VISCOSITY = 1.460e-5 # [m²/s] for the atmosphere at sea level. - # see: https://en.wikipedia.org/wiki/Reynolds_number - -@info "Creating polars. This can take several minutes." -tic() - -procs = addprocs() - -try - function normalize!(x, y) - x_min = minimum(x) - x_max = maximum(x) - for i in eachindex(x) - x[i] = (x[i] - x_min) / (x_max - x_min) - y[i] = (y[i] - x_min) / (x_max - x_min) - end - end - - @eval @everywhere using VortexStepMethod, Xfoil, Statistics, SharedArrays - - # alpha_range = deg2rad.(-5:1:20) - # delta_range = deg2rad.(-5:1:20) - - cl_matrix = SharedArray{Float64}((length(alpha_range), length(delta_range))) - cd_matrix = SharedArray{Float64}((length(alpha_range), length(delta_range))) - cm_matrix = SharedArray{Float64}((length(alpha_range), length(delta_range))) - - @everywhere begin - function turn_trailing_edge!(angle, x, y, lower_turn, upper_turn, x_turn) - turn_distance = upper_turn - lower_turn - smooth_idx = [] - rm_idx = [] - - sign = angle > 0 ? 1 : -1 - y_turn = angle > 0 ? upper_turn : lower_turn - for i in eachindex(x) - if x_turn - turn_distance < x[i] < x_turn + turn_distance && sign * y[i] > 0 - append!(smooth_idx, i) - elseif sign * y[i] < 0 && x_turn > x[i] > x_turn - turn_distance - append!(rm_idx, i) - end - if x[i] > x_turn - x_rel = x[i] - x_turn - y_rel = y[i] - y_turn - x[i] = x_turn + x_rel * cos(angle) + y_rel * sin(angle) - y[i] = y_turn - x_rel * sin(angle) + y_rel * cos(angle) - if angle > 0 && x[i] < x_turn - turn_distance/2 && y[i] > lower_turn - append!(rm_idx, i) - elseif angle < 0 && x[i] < x_turn - turn_distance/2 && y[i] < upper_turn - append!(rm_idx, i) - end - end - end - - #TODO: lower and upper is slightly off because of smoothing - lower_i, upper_i = minimum(smooth_idx), maximum(smooth_idx) - for i in smooth_idx - window = min(i - lower_i + 1, upper_i - i + 1) - x[i] = mean(x[i-window:i+window]) - end - deleteat!(x, rm_idx) - deleteat!(y, rm_idx) - nothing - end - - function solve_alpha!(cls, cds, cms, alpha_range, alpha_idxs, delta, re, x_, y_, lower, upper, kite_speed, speed_of_sound, x_turn) - x = deepcopy(x_) - y = deepcopy(y_) - turn_trailing_edge!(delta, x, y, lower, upper, x_turn) - Xfoil.set_coordinates(x, y) - x, y = Xfoil.pane(npan=140) - reinit = true - for (alpha, alpha_idx) in zip(alpha_range, alpha_idxs) - converged = false - cl = 0.0 - cd = 0.0 - # Solve for the given angle of attack - cl, cd, _, cm, converged = Xfoil.solve_alpha(rad2deg(alpha), re; iter=50, reinit=reinit, mach=kite_speed/speed_of_sound, xtrip=(0.05, 0.05)) - reinit = false - if converged - cls[alpha_idx] = cl - cds[alpha_idx] = cd - cms[alpha_idx] = cm - end - end - return nothing - end - - function run_solve_alpha(alpha_range, delta, re, x_, y_, lower, upper, kite_speed, speed_of_sound, x_turn) - @info "solving alpha with trailing edge angle: $(rad2deg(delta)) degrees" - cls = Float64[NaN for _ in alpha_range] - cds = Float64[NaN for _ in alpha_range] - cms = Float64[NaN for _ in alpha_range] - neg_idxs = sort(findall(alpha_range .< 0.0), rev=true) - neg_alpha_range = alpha_range[neg_idxs] - pos_idxs = sort(findall(alpha_range .>= 0.0)) - pos_alpha_range = alpha_range[pos_idxs] - solve_alpha!(cls, cds, cms, neg_alpha_range, neg_idxs, delta, - re, x_, y_, lower, upper, kite_speed, speed_of_sound, x_turn) - solve_alpha!(cls, cds, cms, pos_alpha_range, pos_idxs, delta, - re, x_, y_, lower, upper, kite_speed, speed_of_sound, x_turn) - return cls, cds, cms - end - end - - function get_lower_upper(x, y) - lower_trailing_edge = 0.0 - upper_trailing_edge = 0.0 - min_lower_distance = Inf - min_upper_distance = Inf - for (xi, yi) in zip(x, y) - if yi < 0 - lower_distance = abs(xi - x_turn) - if lower_distance < min_lower_distance - min_lower_distance = lower_distance - lower_trailing_edge = yi - end - else - upper_distance = abs(xi - x_turn) - if upper_distance < min_upper_distance - min_upper_distance = upper_distance - upper_trailing_edge = yi - end - end - end - return lower_trailing_edge, upper_trailing_edge - end - - kite_speed = v_wind - chord_length = area / width - local reynolds_number = kite_speed * chord_length / KINEMATIC_VISCOSITY # https://en.wikipedia.org/wiki/Reynolds_number - - # Read airfoil coordinates from a file. - local x, y = open(foil_path, "r") do f - x = Float64[] - y = Float64[] - for line in eachline(f) - entries = split(chomp(line)) - try - push!(x, parse(Float64, entries[1])) - push!(y, parse(Float64, entries[2])) - catch ArgumentError - println(entries) - end - end - x, y - end - normalize!(x, y) - Xfoil.set_coordinates(x, y) - x, y = Xfoil.pane(npan=140) - lower, upper = get_lower_upper(x, y) - - @everywhere begin - x = $x - y = $y - x_turn = $x_turn - reynolds_number = $reynolds_number - end - - @sync @distributed for j in eachindex(delta_range) - cl_matrix[:, j], cd_matrix[:, j], cm_matrix[:, j] = run_solve_alpha(alpha_range, delta_range[j], - reynolds_number, x, y, lower, upper, kite_speed, SPEED_OF_SOUND, x_turn) - end - cl_matrix = Matrix{Float64}(cl_matrix) - cd_matrix = Matrix{Float64}(cd_matrix) - cm_matrix = Matrix{Float64}(cm_matrix) - - println("Generated lift matrix:") - display(cl_matrix) - println("Generated drag matrix:") - display(cd_matrix) - println("Generated moment matrix:") - display(cm_matrix) - - @info "Relative trailing_edge height: $(upper - lower)" - @info "Reynolds number for flying speed of $kite_speed is $reynolds_number" - - serialize(polar_path, (alpha_range, delta_range, cl_matrix, cd_matrix, cm_matrix)) - -catch e - @info "Removing processes" - rmprocs(procs) - throw(e) -finally - @info "Removing processes" - rmprocs(procs) -end - -toc() diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index bed88f94..13753c5a 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -15,12 +15,14 @@ using Interpolations: Extrapolation using Parameters using Serialization using SharedArrays +using Timers using PreallocationTools using PrecompileTools using Pkg using DifferentiationInterface import SciMLBase: successful_retcode import YAML +using Xfoil # Export public interface export VSMSettings, WingSettings, SolverSettings, vs @@ -266,7 +268,8 @@ end # Include core functionality include("settings.jl") include("wing_geometry.jl") -include("kite_geometry.jl") +include("polars.jl") +include("ram_geometry.jl") include("filament.jl") include("panel.jl") include("body_aerodynamics.jl") diff --git a/src/polars.jl b/src/polars.jl new file mode 100644 index 00000000..bf87e22f --- /dev/null +++ b/src/polars.jl @@ -0,0 +1,168 @@ + +const SPEED_OF_SOUND = 343 # [m/s] at 20 °C, see: https://en.wikipedia.org/wiki/Speed_of_sound +const KINEMATIC_VISCOSITY = 1.460e-5 # [m²/s] for the atmosphere at sea level. + # see: https://en.wikipedia.org/wiki/Reynolds_number + +function normalize_foil!(x, y) + x_min = minimum(x) + x_max = maximum(x) + for i in eachindex(x) + x[i] = (x[i] - x_min) / (x_max - x_min) + y[i] = (y[i] - x_min) / (x_max - x_min) + end +end + +function turn_trailing_edge!(angle, x, y, lower_turn, upper_turn, crease_frac) + turn_distance = upper_turn - lower_turn + smooth_idx = [] + rm_idx = [] + + sign = angle > 0 ? 1 : -1 + y_turn = angle > 0 ? upper_turn : lower_turn + for i in eachindex(x) + if crease_frac - turn_distance < x[i] < crease_frac + turn_distance && sign * y[i] > 0 + append!(smooth_idx, i) + elseif sign * y[i] < 0 && crease_frac > x[i] > crease_frac - turn_distance + append!(rm_idx, i) + end + if x[i] > crease_frac + x_rel = x[i] - crease_frac + y_rel = y[i] - y_turn + x[i] = crease_frac + x_rel * cos(angle) + y_rel * sin(angle) + y[i] = y_turn - x_rel * sin(angle) + y_rel * cos(angle) + if angle > 0 && x[i] < crease_frac - turn_distance/2 && y[i] > lower_turn + append!(rm_idx, i) + elseif angle < 0 && x[i] < crease_frac - turn_distance/2 && y[i] < upper_turn + append!(rm_idx, i) + end + end + end + + #TODO: lower and upper is slightly off because of smoothing + lower_i, upper_i = minimum(smooth_idx), maximum(smooth_idx) + for i in smooth_idx + window = min(i - lower_i + 1, upper_i - i + 1) + x[i] = mean(x[i-window:i+window]) + end + deleteat!(x, rm_idx) + deleteat!(y, rm_idx) + nothing +end + +function solve_alpha!(cls, cds, cms, alpha_range, alpha_idxs, delta, re, x_, y_, lower, upper, kite_speed, speed_of_sound, crease_frac) + x = deepcopy(x_) + y = deepcopy(y_) + turn_trailing_edge!(delta, x, y, lower, upper, crease_frac) + Xfoil.set_coordinates(x, y) + x, y = Xfoil.pane(npan=140) + reinit = true + for (alpha, alpha_idx) in zip(alpha_range, alpha_idxs) + converged = false + cl = 0.0 + cd = 0.0 + # Solve for the given angle of attack + cl, cd, _, cm, converged = Xfoil.solve_alpha(rad2deg(alpha), re; iter=50, reinit=reinit, mach=kite_speed/speed_of_sound, xtrip=(0.05, 0.05)) + reinit = false + if converged + cls[alpha_idx] = cl + cds[alpha_idx] = cd + cms[alpha_idx] = cm + end + end + return nothing +end + +function run_solve_alpha(alpha_range, delta, re, x_, y_, lower, upper, kite_speed, speed_of_sound, crease_frac) + @info "solving alpha with trailing edge angle: $(rad2deg(delta)) degrees" + cls = Float64[NaN for _ in alpha_range] + cds = Float64[NaN for _ in alpha_range] + cms = Float64[NaN for _ in alpha_range] + neg_idxs = sort(findall(alpha_range .< 0.0), rev=true) + neg_alpha_range = alpha_range[neg_idxs] + pos_idxs = sort(findall(alpha_range .>= 0.0)) + pos_alpha_range = alpha_range[pos_idxs] + + solve_alpha!(cls, cds, cms, neg_alpha_range, neg_idxs, delta, + re, x_, y_, lower, upper, kite_speed, speed_of_sound, crease_frac) + solve_alpha!(cls, cds, cms, pos_alpha_range, pos_idxs, delta, + re, x_, y_, lower, upper, kite_speed, speed_of_sound, crease_frac) + return cls, cds, cms +end + +function get_lower_upper(x, y, crease_frac) + lower_trailing_edge = 0.0 + upper_trailing_edge = 0.0 + min_lower_distance = Inf + min_upper_distance = Inf + for (xi, yi) in zip(x, y) + if yi < 0 + lower_distance = abs(xi - crease_frac) + if lower_distance < min_lower_distance + min_lower_distance = lower_distance + lower_trailing_edge = yi + end + else + upper_distance = abs(xi - crease_frac) + if upper_distance < min_upper_distance + min_upper_distance = upper_distance + upper_trailing_edge = yi + end + end + end + return lower_trailing_edge, upper_trailing_edge +end + +function create_polars(; dat_path, polar_path, wind_vel, area, width, crease_frac, alpha_range, delta_range) + @info "Creating polars. This can take several minutes." + tic() + + cl_matrix = zeros(length(alpha_range), length(delta_range)) + cd_matrix = zeros(length(alpha_range), length(delta_range)) + cm_matrix = zeros(length(alpha_range), length(delta_range)) + + kite_speed = wind_vel + chord_length = area / width + local reynolds_number = kite_speed * chord_length / KINEMATIC_VISCOSITY # https://en.wikipedia.org/wiki/Reynolds_number + + # Read airfoil coordinates from a file. + local x, y = open(dat_path, "r") do f + x = Float64[] + y = Float64[] + for line in eachline(f) + entries = split(chomp(line)) + try + push!(x, parse(Float64, entries[1])) + push!(y, parse(Float64, entries[2])) + catch ArgumentError + println(entries) + end + end + x, y + end + normalize_foil!(x, y) + Xfoil.set_coordinates(x, y) + x, y = Xfoil.pane(npan=140) + lower, upper = get_lower_upper(x, y, crease_frac) + + for j in eachindex(delta_range) + cl_matrix[:, j], cd_matrix[:, j], cm_matrix[:, j] = run_solve_alpha(alpha_range, delta_range[j], + reynolds_number, x, y, lower, upper, kite_speed, SPEED_OF_SOUND, crease_frac) + end + cl_matrix = Matrix{Float64}(cl_matrix) + cd_matrix = Matrix{Float64}(cd_matrix) + cm_matrix = Matrix{Float64}(cm_matrix) + + println("Generated lift matrix:") + display(cl_matrix) + println("Generated drag matrix:") + display(cd_matrix) + println("Generated moment matrix:") + display(cm_matrix) + + @info "Relative trailing_edge height: $(upper - lower)" + @info "Reynolds number for flying speed of $kite_speed is $reynolds_number" + + serialize(polar_path, (alpha_range, delta_range, cl_matrix, cd_matrix, cm_matrix)) + + toc() +end \ No newline at end of file diff --git a/src/kite_geometry.jl b/src/ram_geometry.jl similarity index 98% rename from src/kite_geometry.jl rename to src/ram_geometry.jl index 9f2497c3..f08db6a5 100644 --- a/src/kite_geometry.jl +++ b/src/ram_geometry.jl @@ -506,11 +506,7 @@ function RamAirWing( if !ispath(polar_path) width = 2gamma_tip * radius area = area_interp(gamma_tip) - @eval Main begin - foil_path, polar_path, v_wind, area, width, x_turn, alpha_range, delta_range = - $dat_path, $polar_path, $wind_vel, $gamma_tip, $width, $crease_frac, $alpha_range, $delta_range - include(joinpath(dirname(@__FILE__), "../scripts/polars.jl")) - end + create_polars(; dat_path, polar_path, wind_vel, area, width, crease_frac, alpha_range, delta_range) end (alpha_range, delta_range, cl_matrix::Matrix, cd_matrix::Matrix, cm_matrix::Matrix) = deserialize(polar_path) From 73348d186a51af33608d0d241258d08e7b5fe00b Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Tue, 8 Apr 2025 15:34:32 +0200 Subject: [PATCH 2/4] Remove unused deps --- Project.toml | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 09f37f97..b35d4a43 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,6 @@ PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" -SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Timers = "21f18d07-b854-4dab-86f0-c15a3821819a" @@ -59,7 +58,6 @@ PreallocationTools = "0.4.25" PrecompileTools = "1.2.1" SciMLBase = "2.77.0" Serialization = "1" -SharedArrays = "1" StaticArrays = "1" Statistics = "1" Test = "1" @@ -74,9 +72,8 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "DataFrames", "CSV", "Distributed", "Documenter", "BenchmarkTools", "ControlPlots", "Aqua"] +test = ["Test", "DataFrames", "CSV", "Documenter", "BenchmarkTools", "ControlPlots", "Aqua"] From 809e8a074ac7f243d03bcc84bc2ec8c7b96d4357 Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Tue, 8 Apr 2025 15:40:10 +0200 Subject: [PATCH 3/4] Remove unused deps --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index b35d4a43..7c8b2f90 100644 --- a/Project.toml +++ b/Project.toml @@ -43,7 +43,6 @@ DataFrames = "1.7" DefaultApplication = "1" DelimitedFiles = "1" DifferentiationInterface = "0.6.48" -Distributed = "1.10" Documenter = "1.8" FiniteDiff = "2.27.0" Interpolations = "0.15" From b15d081bc718863cff31570df4b8d24263500c2e Mon Sep 17 00:00:00 2001 From: 1-Bart-1 Date: Tue, 8 Apr 2025 15:46:01 +0200 Subject: [PATCH 4/4] Remove unused import --- src/VortexStepMethod.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 13753c5a..437b4c6d 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -14,7 +14,6 @@ using Interpolations using Interpolations: Extrapolation using Parameters using Serialization -using SharedArrays using Timers using PreallocationTools using PrecompileTools