Merge remote-tracking branch 'origin/upstream-master' into cs_safe

.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:

docs/Project.toml

@@ -11,4 +11,4 @@ SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-Documenter = "0.25"
+Documenter = "1.0"

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 <aning@byu.edu>",
+    warnonly = Documenter.except(:linkcheck, :footnote),
 )
 
 deploydocs(

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.
 

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)
@@ -87,15 +87,19 @@ end
 
 
 # convenience function to access fields within an array of structs
-function Base.getproperty(obj::AbstractVector{<:Section}, sym::Symbol)
-    return getfield.(obj, sym)
+function Base.getproperty(obj::Vector{<:Section}, sym::Symbol)
+    if sym in (:ref, :size)
+        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
@@ -124,8 +128,12 @@ 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)
-    return getfield.(obj, sym)
+function Base.getproperty(obj::Vector{<:OperatingPoint}, sym::Symbol)
+    if sym in (:ref, :size)
+        return getfield(obj, sym)
+    else
+        return getfield.(obj, sym)
+    end
 end
 
 
@@ -176,8 +184,12 @@ 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)
-    return getfield.(obj, sym)
+function Base.getproperty(obj::Vector{<:Outputs}, sym::Symbol)
+    if sym in (:ref, :size)
+        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
@@ -258,7 +270,7 @@ function residual_and_outputs(phi, x, p)  #rotor, section, op)
         a = zero(phi)
         ap = zero(phi)
         R = sign(Vx) + kp
-    
+
     else
 
         if real(phi) < 0
@@ -272,8 +284,8 @@ function residual_and_outputs(phi, x, p)  #rotor, section, op)
         if real(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,8 +317,8 @@ 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:
@@ -519,13 +531,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
@@ -558,7 +570,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)
@@ -633,7 +645,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**

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"