Merge branch 'main' into compathelper/new_version/2024-11-10-00-32-34…

…-766-00123223224

Project.toml

@@ -33,17 +33,22 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BenchmarkTools = "1"
+Colors = "0.13"
 FLOWMath = "0.4"
 FastGaussQuadrature = "1"
 FiniteDiff = "2"
 FixedPoint = "1"
 ForwardDiff = "0.10"
 ImplicitAD = "0.3"
+LaTeXStrings = "1"
 LineSearches = "7"
 MINPACK = "1"
 NLsolve = "4"
 PreallocationTools = "0.4"
+ProgressMeter = "1"
 QuadGK = "2"
+RecipesBase = "1"
+Roots = "2"
 SIAMFANLEquations = "1"
 SimpleNonlinearSolve = "1, 2"
 SpecialFunctions = "2"

src/analysis/setup.jl

@@ -118,10 +118,10 @@ function setup_analysis(
 
     # copy over operating point
     for f in fieldnames(typeof(operating_point))
-        if f!=:units
-        solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f)
+        if f != :units
+            solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f)
+        end
     end
-end
 
     ##### ----- PERFORM PREPROCESSING COMPUTATIONS ----- #####
     if options.verbose

src/postprocess/boundary_layer_utils.jl

@@ -57,21 +57,32 @@ function split_at_stagnation_point(duct_panel_lengths, duct_panel_tangents, Vtot
         dp[1] = dp[2]
     end
 
-    # interpolate the lengths between control points
-    sum_length = 0.5 * sum(duct_panel_lengths[stag_ids])
-    stag_interp = FLOWMath.linear(dp, [0.0, sum_length], 0.0)
+    if dp[1] == dp[2]
+        # we're likely in a hover-ish case and there is no stagnation point.
 
-    partial_panel_lengths = [stag_interp, sum_length - stag_interp]
+        s_upper = nothing
+        s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[end:-1:1])
+        split_ratio = 1.0
+        stag_ids .= length(duct_panel_lengths)
+    else
 
-    s_upper = arc_lengths_from_panel_lengths(
-        [abs(partial_panel_lengths[2]); duct_panel_lengths[stag_ids[2]:end]]
-    )
+        # interpolate the lengths between control points
+        sum_length = 0.5 * sum(duct_panel_lengths[stag_ids])
+        stag_interp = FLOWMath.linear(dp, [0.0, sum_length], 0.0)
 
-    s_lower = arc_lengths_from_panel_lengths(
-        [abs(partial_panel_lengths[1]); duct_panel_lengths[stag_ids[1]:-1:1]]
-    )
+        partial_panel_lengths = [stag_interp, sum_length - stag_interp]
+
+        s_upper = arc_lengths_from_panel_lengths(
+            [abs(partial_panel_lengths[2]); duct_panel_lengths[stag_ids[2]:end]]
+        )
+
+        s_lower = arc_lengths_from_panel_lengths(
+            [abs(partial_panel_lengths[1]); duct_panel_lengths[stag_ids[1]:-1:1]]
+        )
+        split_ratio = stag_interp / sum_length
+    end
 
-    return s_upper, s_lower, stag_ids, stag_interp / sum_length
+    return s_upper, s_lower, stag_ids, split_ratio
 end
 
 """

src/postprocess/postprocess.jl

@@ -517,7 +517,7 @@ function post_process(
             body_vortex_panels.tangent[:, 1:Int(body_vortex_panels.npanel[1])],
             body_vortex_panels.node[2, Int(body_vortex_panels.nnode[1])],
             operating_point;
-            verbose=verbose,
+            verbose=false,
         )
 
         body_thrust[1] -= duct_viscous_drag

src/postprocess/viscous_drag.jl

@@ -20,7 +20,7 @@ Note: if no separation occurs, the inputs are simply the final values for the bo
 """
 function squire_young(d2, Ue, Uinf, H12)
     # note: formula also divides by chord, but we're going to multiply by chord later.
-    return 2.0 * d2 * (Ue / Uinf)^((5.0 + H12) / 2.0)
+    return 2.0 * d2 * FLOWMath.abs_smooth(Ue / Uinf, 2.0*eps())^((5.0 + H12) / 2.0)
 end
 
 """
@@ -183,13 +183,15 @@ function compute_viscous_drag_duct(
     )
 
     # set up boundary layer solve parameters
-    boundary_layer_functions_upper = setup_boundary_layer_functions_head(
-        s_upper,
-        [0.0; Vtan_duct[stag_ids[2]:end]],
-        operating_point,
-        boundary_layer_options;
-        verbose=verbose,
-    )
+    if !isnothing(s_upper)
+        boundary_layer_functions_upper = setup_boundary_layer_functions_head(
+            s_upper,
+            [0.0; Vtan_duct[stag_ids[2]:end]],
+            operating_point,
+            boundary_layer_options;
+            verbose=verbose,
+        )
+    end
 
     # set up boundary layer solve parameters
     boundary_layer_functions_lower = setup_boundary_layer_functions_head(
@@ -201,13 +203,15 @@ function compute_viscous_drag_duct(
     )
 
     # - Set integration steps - #
-    # upper side
-    upper_steps =
-        set_boundary_layer_steps(
-            boundary_layer_options.n_steps,
-            boundary_layer_options.first_step_size,
-            s_upper[end] - boundary_layer_options.offset,
-        ) .+ boundary_layer_options.offset
+    if !isnothing(s_upper)
+        # upper side
+        upper_steps =
+            set_boundary_layer_steps(
+                boundary_layer_options.n_steps,
+                boundary_layer_options.first_step_size,
+                s_upper[end] - boundary_layer_options.offset,
+            ) .+ boundary_layer_options.offset
+    end
 
     # lower side
     lower_steps =
@@ -219,19 +223,26 @@ function compute_viscous_drag_duct(
 
     # - Get drag coeffients - #
 
-    cdc_upper, usol_upper, stepsol_upper, s_sep_upper = compute_single_side_drag_coefficient_head(
-        upper_steps,
-        exit_radius,
-        operating_point,
-        boundary_layer_functions_upper,
-        (;
-            boundary_layer_options.rk,
-            boundary_layer_options.separation_criteria,
-            separation_allowance=boundary_layer_options.separation_allowance_upper,
-            separation_penalty=boundary_layer_options.separation_penalty_upper,
-        );
-        verbose=verbose,
-    )
+    if !isnothing(s_upper)
+        cdc_upper, usol_upper, stepsol_upper, s_sep_upper = compute_single_side_drag_coefficient_head(
+            upper_steps,
+            exit_radius,
+            operating_point,
+            boundary_layer_functions_upper,
+            (;
+                boundary_layer_options.rk,
+                boundary_layer_options.separation_criteria,
+                separation_allowance=boundary_layer_options.separation_allowance_upper,
+                separation_penalty=boundary_layer_options.separation_penalty_upper,
+            );
+            verbose=verbose,
+        )
+    else
+        cdc_upper = 0.0
+        usol_upper = -ones(eltype(Vtan_duct), 3)
+        stepsol_upper = -1
+        s_sep_upper = -1.0
+    end
 
     cdc_lower, usol_lower, stepsol_lower, s_sep_lower = compute_single_side_drag_coefficient_head(
         lower_steps,

src/preprocess/geometry/rotor_geometry.jl

@@ -242,10 +242,10 @@ function get_blade_ends_from_body_geometry!(
     if !silence_warnings
         for (irotor, (R, r)) in enumerate(zip(Rtip, Rhub))
             if r !== rp_centerbody_coordinates[2, ihub[irotor]]
-                @warn "Overwriting Rhub for rotor $(irotor) to place it at the centerbody wall.  Moving from $(r) to $(rp_centerbody_coordinates[2, ihub[irotor]])"
+                @warn "Overwriting Rhub for rotor $(ForwardDiff.value(irotor)) to place it at the centerbody wall.  Moving from $(ForwardDiff.value(r)) to $(ForwardDiff.value(rp_centerbody_coordinates[2, ihub[irotor]]))"
             end
             if R !== rp_duct_coordinates[2, iduct[irotor]] .- tip_gaps[irotor]
-                @warn "Overwriting Rtip for rotor $(irotor) to place it at the correct tip gap relative to the casing wall. Moving from $(R) to $(rp_duct_coordinates[2, iduct[irotor]] .- tip_gaps[irotor])"
+                @warn "Overwriting Rtip for rotor $(ForwardDiff.value(irotor)) to place it at the correct tip gap relative to the casing wall. Moving from $(ForwardDiff.value(R)) to $(ForwardDiff.value(rp_duct_coordinates[2, iduct[irotor]] .- tip_gaps[irotor]))"
             end
         end
     end

src/utilities/bookkeeping.jl

@@ -65,7 +65,6 @@ function get_problem_dimensions(paneling_constants::PanelingConstants)
     nbe = nws - 1
 
     # number of body panels
-    # TODO: need number of casing and centerbody nodes separately as well
     ncp = nduct_inlet
     ncbp = ncenterbody_inlet
     # add rest of panels mutual between centerbody and duct

src/utilities/inputs.jl

@@ -277,7 +277,7 @@ function Rotor(
     fliplift;
     i_know_what_im_doing=false,
 )
-    if length(findall(t -> t > pi / 2, twists)) > 2
+    if length(findall(t -> t > 1.75, twists)) > 2
         @warn "It looks like your input twist angles may be in degrees. Note that the required units for twist are radians. Converting to radians for you (set the `i_know_what_im_doing` keyword argument to true to disable automatic conversion)."
         twists .*= pi / 180.0
     end

src/utilities/misc.jl

@@ -36,11 +36,21 @@ end
 
 Formatted printing when you have lots of debugging to do and want things to line up nicely.
 """
-function printdebug(variable_name, variable, nspaces=4)
-    @printf "%*s %-14s %2.12f\n" nspaces " " variable_name variable
+function printdebug(variable_name, variable, nspaces=4; colw=16)
+    if length(variable) == 1
+        s = @sprintf "%*s %-*s %2.12f" nspaces " " colw variable_name variable[1]
+        println(rstrip(s, ['0', ' ']))
+    else
+        @printf "%*s %-*s\n" nspaces " " colw variable_name
+        for v in variable
+            s = @sprintf "%*s %-*s %2.12f" (nspaces) " " colw "" v
+            println(rstrip(s, ['0', ' ']))
+        end
+    end
     # println("\t$(variable_name)" * v)
     return nothing
 end
+
 """
     dot(A, B) = sum(a * b for (a, b) in zip(A, B))
 
@@ -150,32 +160,33 @@ function reset_containers!(c; exception_keys=[])
 end
 
 """
-    promote_ducted_rotor_type(ducted_rotor)
+    promote_ducted_rotor_type(ducted_rotor, operating_point)
 
 Convenience function for promoting types based on any potential elements of the ducted_rotor object dependent on optimization design variables.
 
 # Arguments
 - `ducted_rotor::DuctedRotor` : the ducted_rotor input
+- `operating_point::OperatingPoint` : the operating_point input
 
 # Returns
 - `TP::Type` : the promoted type
 """
-function promote_ducted_rotor_type(p)
+function promote_ducted_rotor_type(d,o)
     return promote_type(
-        eltype(p.duct_coordinates),
-        eltype(p.centerbody_coordinates),
-        eltype(p.operating_point.Vinf),
-        eltype(p.operating_point.rhoinf),
-        eltype(p.operating_point.muinf),
-        eltype(p.operating_point.asound),
-        eltype(p.operating_point.Omega),
-        eltype(p.rotor.B),
-        eltype(p.rotor.rotorzloc),
-        eltype(p.rotor.r),
-        eltype(p.rotor.Rhub),
-        eltype(p.rotor.Rtip),
-        eltype(p.rotor.chords),
-        eltype(p.rotor.twists),
+        eltype(d.duct_coordinates),
+        eltype(d.centerbody_coordinates),
+        eltype(o.Vinf),
+        eltype(o.rhoinf),
+        eltype(o.muinf),
+        eltype(o.asound),
+        eltype(o.Omega),
+        eltype(d.rotor.B),
+        eltype(d.rotor.rotorzloc),
+        eltype(d.rotor.r),
+        eltype(d.rotor.Rhub),
+        eltype(d.rotor.Rtip),
+        eltype(d.rotor.chords),
+        eltype(d.rotor.twists),
     )
 end
 

src/utilities/options.jl

@@ -660,7 +660,7 @@ end
 
 # Fields:
 - `model_drag::Tb=false` : flag to turn on viscous drag approximation
-- `n_steps::Int = Int(2e2)` : number of steps to use in boundary layer integration
+- `n_steps::Int = Int(5e2)` : number of steps to use in boundary layer integration
 - `first_step_size::Float = 1e-6` : size of first step in boundary layer integration
 - `offset::Float = 1e-3` : size of offset for (where to initialize) boundary layer integration
 - `rk::Function = RK4` : solver to use for boundary layer integration (RK4 or RK2 available)
@@ -672,7 +672,7 @@ end
 """
 @kwdef struct HeadsBoundaryLayerOptions{Tb,Tf,Tfun,Ti,To,Tp,Ts} <: BoundaryLayerOptions
     model_drag::Tb=false
-    n_steps::Ti = Int(2e2)
+    n_steps::Ti = Int(5e2)
     first_step_size::Tf = 1e-6
     offset::To = 1e-3
     rk::Tfun = RK2

src/visualization/plot_recipes.jl

@@ -160,9 +160,8 @@ end
     ylim --> (0.0, maximum(bvp.node[2, :]) * 1.05)
 
     if discrete_labels
-        xticks --> DuctAPE.determine_geometry_xlabels(
-            bvp, rsp, wvp; tol=xticktol, tickdigits=tickdigits
-        )
+        xticks -->
+        determine_geometry_xlabels(bvp, rsp, wvp; tol=xticktol, tickdigits=tickdigits)
         yticks --> determine_geometry_ylabels(bvp, rsp; tol=yticktol, tickdigits=tickdigits)
         xgrid --> true
         ygrid --> true
@@ -406,8 +405,10 @@ function determine_geometry_ylabels(bvp, rsp; tol=1e-2, tickdigits=2)
 
     # Rotor
     for r in 1:Int(rsp.nbodies[])
-        push!(yt, rsp.node[2, Int(rsp.endnodeidxs[r, r])])
-        push!(yl, @sprintf "%3.*f" tickdigits rsp.node[2, Int(rsp.endnodeidxs[r, r])])
+        for i in 1:2
+            push!(yt, rsp.node[2, Int(rsp.endnodeidxs[i, r])])
+            push!(yl, @sprintf "%3.*f" tickdigits rsp.node[2, Int(rsp.endnodeidxs[i, r])])
+        end
     end
 
     # Centerbody
@@ -416,18 +417,22 @@ function determine_geometry_ylabels(bvp, rsp; tol=1e-2, tickdigits=2)
         yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits)
     end
 
-    # Duct front/back
-    for f in [findmax, findmin]
-        tmp = f(bvp.node[1, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])])
-        tmpy = bvp.node[2, tmp[2]]
-        yt, yl = add_ticks(yt, yl, tmpy, tol, tickdigits)
-    end
+    # # Duct front/back
+    # for f in [findmax, findmin]
+    #     tmp = f(bvp.node[1, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])])
+    #     tmpy = bvp.node[2, tmp[2]]
+    #     yt, yl = add_ticks(yt, yl, tmpy, tol, tickdigits)
+    # end
 
-    # Duct inner/outer
-    for f in [maximum, minimum]
-        tmp = f(bvp.node[2, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])])
-        yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits)
-    end
+    # Duct Exit Radius
+    tmpy = bvp.node[2, 1]
+    yt, yl = add_ticks(yt, yl, tmpy, tol, tickdigits)
+
+    # # Duct inner/outer
+    # for f in [maximum, minimum]
+    #     tmp = f(bvp.node[2, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])])
+    #     yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits)
+    # end
 
     return (yt, yl)
 end
@@ -712,7 +717,7 @@ end
     # Plot specific values
     aspect_ratio --> 1
 
-    p, _ = DuctAPE.calculate_streamlines(
+    p, _ = calculate_streamlines(
         bvp,
         bvs,
         wvp,