From da260cc43c03a8ec79d4ee90d96b5fd5729aa7e9 Mon Sep 17 00:00:00 2001 From: Andrew Ning Date: Mon, 27 Jan 2025 22:13:50 -0700 Subject: [PATCH 1/5] fix on julia 1.11 for extracting data from array of structs --- docs/src/tutorial.md | 8 ++++-- src/CCBlade.jl | 64 ++++++++++++++++++++++++++------------------ test/runtests.jl | 36 ++++++++++++------------- 3 files changed, 62 insertions(+), 46 deletions(-) diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 3c9b0fe..20d3914 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -2,12 +2,16 @@ This starter tutorial walks through the mechanics of running a analysis. This is designing as a starting point and does not explain every option or consideration. More specific and advanced usage are described in the [how to guide](howto.md). -We will simulate the APC thin electric 10 x 5 propeller. The geometry, and wind tunnel data for this propeller is available from [UIUC](https://m-selig.ae.illinois.edu/props/volume-1/propDB-volume-1.html#APC). Let's load CCBlade and a plotting package (I chose to use PyPlot in this example). +We will simulate the APC thin electric 10 x 5 propeller. The geometry, and wind tunnel data for this propeller is available from [UIUC](https://m-selig.ae.illinois.edu/props/volume-1/propDB-volume-1.html#APC). Let's load CCBlade and a plotting package (I chose to use PyPlot in this example). ```@setup prop using CCBlade using PyPlot ``` +```julia +using CCBlade +using PyPlot +``` There are two parts of the geometry we need to define: the first are quantities for the whole rotor, and the second part defines properties that vary along the blade radius. For the rotor, we only need to define the hub radius, tip radius, and the number of blades. This is a two-bladed propeller, and the 10 in the name (10 x 5) means it has a diameter of 10 inches (and so the tip radius is half of that). I prefer to do all calculations in metric units so I'll convert it. From the [geometry table](https://m-selig.ae.illinois.edu/props/volume-1/data/apce_10x5_geom.txt) for this propeller we see that the hub radius is less than ``0.15 R_{tip}``. It isn't defined exactly, and is less critical, we'll assume ``0.10 R_{tip}`` for this example. @@ -151,7 +155,7 @@ J &= \frac{V}{n D} ``` !!! note - Efficiency is set to zero if the thrust is negative (producing drag). + Efficiency is set to zero if the thrust is negative (producing drag). The code below performs this analysis then plots thrust coefficient, power coefficient, and efficiency as a function of advance ratio as compared to the experimental data. diff --git a/src/CCBlade.jl b/src/CCBlade.jl index 5a3f21f..f9bd4e5 100644 --- a/src/CCBlade.jl +++ b/src/CCBlade.jl @@ -26,10 +26,10 @@ include("airfoils.jl") # all the code related to airfoil data # --------- structs ------------- """ - Rotor(Rhub, Rtip, B; precone=0.0, turbine=false, + Rotor(Rhub, Rtip, B; precone=0.0, turbine=false, mach=nothing, re=nothing, rotation=nothing, tip=PrandtlTipHub()) -Parameters defining the rotor (apply to all sections). +Parameters defining the rotor (apply to all sections). **Arguments** - `Rhub::Float64`: hub radius (along blade length) @@ -42,8 +42,8 @@ Parameters defining the rotor (apply to all sections). - `rotation::RotationCorrection`: correction method for blade rotation - `tip::TipCorrection`: correction method for hub/tip loss """ -struct Rotor{TF, TI, TB, - T1 <: Union{Nothing, MachCorrection}, T2 <: Union{Nothing, ReCorrection}, +struct Rotor{TF, TI, TB, + T1 <: Union{Nothing, MachCorrection}, T2 <: Union{Nothing, ReCorrection}, T3 <: Union{Nothing, RotationCorrection}, T4 <: Union{Nothing, TipCorrection}} Rhub::TF Rtip::TF @@ -65,7 +65,7 @@ Rotor(Rhub, Rtip, B; precone=0.0, turbine=false, mach=nothing, re=nothing, rotat Section(r, chord, theta, af) Define sectional properties for one station along rotor - + **Arguments** - `r::Float64`: radial location along blade - `chord::Float64`: corresponding local chord length @@ -77,7 +77,7 @@ struct Section{TF, TAF} chord::TF theta::TF af::TAF -end +end # promote to same type, e.g., duals function Section(r, chord, theta, af) @@ -88,14 +88,18 @@ end # convenience function to access fields within an array of structs function Base.getproperty(obj::AbstractVector{<:Section}, sym::Symbol) - return getfield.(obj, sym) + if sym == :ref + return getfield(obj, sym) + else + return getfield.(obj, sym) + end end # This is not always type stable b/c we don't know if the return type will be float or af function. """ OperatingPoint(Vx, Vy, rho; pitch=0.0, mu=1.0, asound=1.0) -Operation point for a rotor. -The x direction is the axial direction, and y direction is the tangential direction in the rotor plane. +Operation point for a rotor. +The x direction is the axial direction, and y direction is the tangential direction in the rotor plane. See Documentation for more detail on coordinate systems. `Vx` and `Vy` vary radially at same locations as `r` in the rotor definition. @@ -111,7 +115,7 @@ struct OperatingPoint{TF} Vx::TF Vy::TF rho::TF - pitch::TF + pitch::TF mu::TF asound::TF end @@ -125,7 +129,11 @@ OperatingPoint(Vx, Vy, rho; pitch=zero(rho), mu=one(rho), asound=one(rho)) = Ope # convenience function to access fields within an array of structs function Base.getproperty(obj::AbstractVector{<:OperatingPoint}, sym::Symbol) - return getfield.(obj, sym) + if sym == :ref + return getfield(obj, sym) + else + return getfield.(obj, sym) + end end @@ -177,7 +185,11 @@ Outputs() = Outputs(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, # convenience function to access fields within an array of structs function Base.getproperty(obj::AbstractVector{<:Outputs}, sym::Symbol) - return getfield.(obj, sym) + if sym == :ref + return getfield(obj, sym) + else + return getfield.(obj, sym) + end end # ------------------------------- @@ -195,7 +207,7 @@ function residual_and_outputs(phi, x, p) #rotor, section, op) # unpack inputs r, chord, theta, Rhub, Rtip, Vx, Vy, rho, pitch, mu, asound = x # variables af, B, turbine, re_corr, mach_corr, rotation_corr, tip_corr = p # parameters - + # constants sigma_p = B*chord/(2.0*pi*r) sphi = sin(phi) @@ -235,7 +247,7 @@ function residual_and_outputs(phi, x, p) #rotor, section, op) # hub/tip loss F = 1.0 if !isnothing(tip_corr) - F = tip_correction(tip_corr, r, Rhub, Rtip, phi, B) + F = tip_correction(tip_corr, r, Rhub, Rtip, phi, B) end # sec parameters @@ -252,13 +264,13 @@ function residual_and_outputs(phi, x, p) #rotor, section, op) R = sign(phi) - k elseif isapprox(Vy, 0.0, atol=1e-6) - + u = zero(phi) v = k*ct/cn*abs(Vx) a = zero(phi) ap = zero(phi) R = sign(Vx) + kp - + else if phi < 0 @@ -272,8 +284,8 @@ function residual_and_outputs(phi, x, p) #rotor, section, op) if k >= -2.0/3 # momentum region a = k/(1 - k) - else # empirical region. Not Buhl's correction but instead uses Buhl with F = 1 then multiplied by F. - # (Buhl(F = 1)*F). The original method does not force CT -> 0 as f->0. This can be problematic if + else # empirical region. Not Buhl's correction but instead uses Buhl with F = 1 then multiplied by F. + # (Buhl(F = 1)*F). The original method does not force CT -> 0 as f->0. This can be problematic if # using CT/CP directly as a design variables. Suggestion courtesy of Kenneth Lønbæk. g1 = 2*k + 1.0/9 g2 = -2*k - 1.0/3 @@ -305,14 +317,14 @@ function residual_and_outputs(phi, x, p) #rotor, section, op) Np = cn*0.5*rho*W^2*chord Tp = ct*0.5*rho*W^2*chord - # The BEM methodology applies hub/tip losses to the loads rather than to the velocities. - # This is the most common way to implement a BEM, but it means that the raw velocities are misleading + # The BEM methodology applies hub/tip losses to the loads rather than to the velocities. + # This is the most common way to implement a BEM, but it means that the raw velocities are misleading # as they do not contain any hub/tip loss corrections. # To fix this we compute the effective hub/tip losses that would produce the same thrust/torque. # In other words: # CT = 4 a (1 + a) F = 4 a G (1 + a G)\n # This is solved for G, then multiplied against the wake velocities. - + if isapprox(Vx, 0.0, atol=1e-6) G = sqrt(F) elseif isapprox(Vy, 0.0, atol=1e-6) @@ -512,13 +524,13 @@ function solve(rotor, section, op; npts=10, forcebackwardsearch=false, epsilon_e end _, outputs = residual_and_outputs(phistar, xv, pv) return outputs - end - end + end + end # it shouldn't get to this point. if it does it means no solution was found # it will return empty outputs # alternatively, one could increase npts and try again - + @warn "Invalid data (likely) for this section. Zero loading assumed." return Outputs() end @@ -551,7 +563,7 @@ function simple_op(Vinf, Omega, r, rho; pitch=zero(rho), mu=one(rho), asound=one error("You passed in an vector for r, but this function does not accept an vector.\nProbably you intended to use broadcasting") end - Vx = Vinf * cos(precone) + Vx = Vinf * cos(precone) Vy = Omega * r * cos(precone) return OperatingPoint(Vx, Vy, rho, pitch, mu, asound) @@ -626,7 +638,7 @@ end """ thrusttorque(rotor, sections, outputs::AbstractVector{TO}) where TO -integrate the thrust/torque across the blade, +integrate the thrust/torque across the blade, including 0 loads at hub/tip, using a trapezoidal rule. **Arguments** diff --git a/test/runtests.jl b/test/runtests.jl index 46f7d29..67eea9a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,7 +37,7 @@ function affunc(alpha, Re, M) cl = 0.084*alpha*180/pi return cl, 0.0 -end +end sections = Section.(r, chord, theta, Ref(affunc)) @@ -79,8 +79,8 @@ betavec = 90 .- out.phi*180/pi @test isapprox(betavec[6], 84.7113, atol=1e-3) # using my more converged solution -# -# +# +# # idx = 6 @@ -143,7 +143,7 @@ function affunc2(alpha, Re, M) cd = 0.008 - 0.003*cl + 0.01*cl*cl return cl, cd -end +end sections = Section.(r, chord, theta, Ref(affunc2)) @@ -159,7 +159,7 @@ qsim = 1e2*[0.803638686218187, 0.806984572453978, 0.809709290183008, 0.811743686 for i = 1:60 Vinf = float(i) - Omega = RPM * pi/30 + Omega = RPM * pi/30 ops = simple_op.(Vinf, Omega, r, rho) @@ -209,13 +209,13 @@ function affunc3(alpha, Re, M) cd = 0.008 - 0.003*cl + 0.01*cl*cl return cl, cd -end +end sections = Section.(r, chord, theta, Ref(affunc3)) Vinf = 5.0 -Omega = RPM * pi/30 +Omega = RPM * pi/30 ops = simple_op.(Vinf, Omega, r, rho) out = solve.(Ref(rotor_no_F), sections, ops) @@ -266,7 +266,7 @@ end # empirical data. Rather they are from the figures used in the documentaiton. Qualitatively # the output is about right. Minor changes in, for example, the airfoil interpolation method # coudl slightly change the outputs. The main purpose of these tests is to alert us if something -# significant changes. +# significant changes. Rhub = 1.5 Rtip = 63.0 @@ -298,7 +298,7 @@ aftypes[8] = AlphaAF("airfoils/NACA64_A17.dat", radians=false) # indices correspond to which airfoil is used at which station af_idx = [1, 1, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 8, 8] -# create airfoil array +# create airfoil array airfoils = aftypes[af_idx] sections = Section.(r, chord, theta, airfoils) @@ -491,7 +491,7 @@ CQ = zeros(nP) for i = 1:nP - + op = simple_op.(Vinf, Omega, r, rho, pitch=pitch[i]) outputs = solve.(Ref(rotor), sections, op) T, Q = thrusttorque(rotor, sections, outputs) @@ -674,7 +674,7 @@ function affunc(alpha, Re, M) cd = 0.008 - 0.003*cl + 0.01*cl*cl return cl, cd -end +end n = length(r) airfoils = fill(affunc, n) @@ -687,10 +687,10 @@ precone = 0.0 rho = 1.225 Vinf = 30.0 RPM = 2100 -Omega = RPM * pi/30 +Omega = RPM * pi/30 function ccbladewrapper(x) - + # unpack nall = length(x) nvec = nall - 7 @@ -726,7 +726,7 @@ J = ForwardDiff.jacobian(ccbladewrapper, x) # using BenchmarkTools # @btime ForwardDiff.jacobian($ccbladewrapper, $x) -# original: 584.041 μs (9910 allocations: 1.15 MiB) +# original: 584.041 μs (9910 allocations: 1.15 MiB) # with ImplicitAD: 323.208 μs (11862 allocations: 747.50 KiB) import FiniteDiff @@ -749,10 +749,10 @@ function ccbladewrapper(x) cl = 6.2*alpha cd = 0.008 - 0.003*cl + 0.01*cl*cl - + return cl, cd - end - + end + # unpack nall = length(x) nvec = nall - 7 @@ -797,7 +797,7 @@ precone = 0.0 rho = 1.225 Vinf = 30.0 RPM = 2100 -Omega = RPM * pi/30 +Omega = RPM * pi/30 xvec = [r; chord; theta; Rhub; Rtip; pitch; precone; Vinf; Omega; rho] From a7723f5827d6fdb80e4338d0fe8252d85019c1c3 Mon Sep 17 00:00:00 2001 From: Andrew Ning Date: Mon, 27 Jan 2025 22:45:58 -0700 Subject: [PATCH 2/5] doc update --- docs/Project.toml | 2 +- docs/make.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 00559b3..383cee9 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -11,4 +11,4 @@ SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] -Documenter = "0.25" \ No newline at end of file +Documenter = "1.0" \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index da21ae5..493086a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,6 +13,7 @@ makedocs( repo="https://github.com/byuflowlab/CCBlade.jl/blob/{commit}{path}#L{line}", sitename="CCBlade.jl", authors="Andrew Ning ", + warnonly = Documenter.except(:linkcheck, :footnote), ) deploydocs( From 9b5ab8a8bc84a481905b81017bd37bdc965726ae Mon Sep 17 00:00:00 2001 From: Andrew Ning Date: Mon, 27 Jan 2025 22:53:16 -0700 Subject: [PATCH 3/5] workflow update --- .github/workflows/docs.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index b363031..0ac2884 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -10,7 +10,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: [1.6.0] + julia-version: [1.11.0] julia-arch: [x86] os: [ubuntu-latest] steps: From 0d4ab0c28bd1aa3059c53b5505c6ab2dfb4d0df5 Mon Sep 17 00:00:00 2001 From: Daniel Ingraham Date: Tue, 28 Jan 2025 10:55:25 -0500 Subject: [PATCH 4/5] Tweaks to `getproperty` convenience functions --- src/CCBlade.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/CCBlade.jl b/src/CCBlade.jl index f9bd4e5..09b1ca5 100644 --- a/src/CCBlade.jl +++ b/src/CCBlade.jl @@ -87,8 +87,8 @@ end # convenience function to access fields within an array of structs -function Base.getproperty(obj::AbstractVector{<:Section}, sym::Symbol) - if sym == :ref +function Base.getproperty(obj::Vector{<:Section}, sym::Symbol) + if sym in (:ref, :size) return getfield(obj, sym) else return getfield.(obj, sym) @@ -128,8 +128,8 @@ OperatingPoint(Vx, Vy, rho, pitch, mu, asound) = OperatingPoint(promote(Vx, Vy, OperatingPoint(Vx, Vy, rho; pitch=zero(rho), mu=one(rho), asound=one(rho)) = OperatingPoint(Vx, Vy, rho, pitch, mu, asound) # convenience function to access fields within an array of structs -function Base.getproperty(obj::AbstractVector{<:OperatingPoint}, sym::Symbol) - if sym == :ref +function Base.getproperty(obj::Vector{<:OperatingPoint}, sym::Symbol) + if sym in (:ref, :size) return getfield(obj, sym) else return getfield.(obj, sym) @@ -184,8 +184,8 @@ Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G) = Outputs(prom Outputs() = Outputs(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) # convenience function to access fields within an array of structs -function Base.getproperty(obj::AbstractVector{<:Outputs}, sym::Symbol) - if sym == :ref +function Base.getproperty(obj::Vector{<:Outputs}, sym::Symbol) + if sym in (:ref, :size) return getfield(obj, sym) else return getfield.(obj, sym) From 10f7ad99a3381862140481645f3908ac9b4f4558 Mon Sep 17 00:00:00 2001 From: Daniel Ingraham Date: Tue, 28 Jan 2025 14:40:10 -0500 Subject: [PATCH 5/5] Add `esoteric arrays` tests --- test/Project.toml | 2 ++ test/runtests.jl | 63 ++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 64 insertions(+), 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index b3d7195..621909d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,6 @@ [deps] +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/runtests.jl b/test/runtests.jl index 67eea9a..7c03018 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -816,4 +816,65 @@ end @test checkstability() -end \ No newline at end of file +end + +using OffsetArrays +using FillArrays + +@testset "esoteric arrays" begin + radii = rand(10) + chord = rand(length(radii)) + theta = rand(length(radii)) + sections = Section.(radii, chord, theta, alpha->(1.0, 1.0)) + + Vx = rand(length(radii)) + Vy = rand(length(radii)) + rho = rand(length(radii)) + pitch = rand(length(radii)) + mu = rand(length(radii)) + asound = rand(length(radii)) + ops = OperatingPoint.(Vx, Vy, rho, pitch, mu, asound) + + Np = rand(length(radii)) + Tp = rand(length(radii)) + a = rand(length(radii)) + ap = rand(length(radii)) + u = rand(length(radii)) + v = rand(length(radii)) + phi = rand(length(radii)) + alpha = rand(length(radii)) + W = rand(length(radii)) + cl = rand(length(radii)) + cd = rand(length(radii)) + cn = rand(length(radii)) + ct = rand(length(radii)) + F = rand(length(radii)) + G = rand(length(radii)) + outs = Outputs.(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G) + + @testset "OffsetArrays" begin + sections_oa = OffsetArray(sections, 0:length(sections)-1) + @test sections[1].r ≈ sections_oa[0].r + + ops_oa = OffsetArray(ops, 0:length(ops)-1) + @test ops[1].Vx ≈ ops_oa[0].Vx + + outs_oa = OffsetArray(outs, 0:length(outs)-1) + @test outs[1].Np ≈ outs_oa[0].Np + end + + @testset "FillArrays" begin + sections_fill = Fill(sections[1], 3) + @test length(sections_fill) == 3 + sections_fill[3].r ≈ sections[1].r + + ops_fill = Fill(ops[1], 3) + @test length(ops_fill) == 3 + ops_fill[3].Vx ≈ ops[1].Vx + + outs_fill = Fill(outs[1], 3) + @test length(outs_fill) == 3 + outs_fill[3].Np ≈ outs[1].Np + end + +end