cleanup and increse boundary layer robustness, also move back to zero…

… efficiency if power or thrust go negative
byuflowlab · Nov 20, 2024 · d7cbf9e · d7cbf9e
1 parent 1b0b2d6
commit d7cbf9e
Show file tree

Hide file tree

Showing 10 changed files with 100 additions and 88 deletions.
diff --git a/src/postprocess/boundary_layer_utils.jl b/src/postprocess/boundary_layer_utils.jl
@@ -47,7 +47,8 @@ function split_at_stagnation_point(duct_panel_lengths, duct_panel_tangents, Vtot
     dp[1] = dot(duct_panel_tangents[:, 1], Vtot_duct[:, 1])
 
     # loop through panels and stop when dot product of panel vector and velocity vector changes sign
-    for i in 2:length(duct_panel_lengths)
+    # TODO: have to start further than the size of the first step size for the ODE solver
+    for i in 3:length(duct_panel_lengths)
         dp[2] = dot(duct_panel_tangents[:, i], Vtot_duct[:, i])
 
         if sign(dp[1]) != sign(dp[2])
@@ -64,6 +65,14 @@ function split_at_stagnation_point(duct_panel_lengths, duct_panel_tangents, Vtot
         s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[end:-1:1])
         split_ratio = 1.0
         stag_ids .= length(duct_panel_lengths)
+
+        elseif stag_ids[2] == length(duct_panel_lengths)
+        # things flipped on the last panel
+
+        s_upper = nothing
+        s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[end-1:-1:1])
+        split_ratio = 1.0
+        stag_ids .= length(duct_panel_lengths)-1
     else
 
         # interpolate the lengths between control points
@@ -79,7 +88,7 @@ function split_at_stagnation_point(duct_panel_lengths, duct_panel_tangents, Vtot
         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
+        split_ratio = s_lower[end]/(s_lower[end]+s_upper[end])
     end
 
     return s_upper, s_lower, stag_ids, split_ratio

diff --git a/src/postprocess/postprocess.jl b/src/postprocess/postprocess.jl
@@ -525,10 +525,10 @@ function post_process(
         boundary_layer_outputs = nothing
     end
 
+
     ### --- TOTAL OUTPUTS --- ###
 
     # - Total Thrust - #
-    total_thrust[] = sum([rotor_inviscid_thrust'; rotor_viscous_thrust'])
     total_thrust[] = sum([rotor_inviscid_thrust'; rotor_viscous_thrust'; body_thrust])
 
     # - Total Torque - #

diff --git a/src/postprocess/rotor_performance.jl b/src/postprocess/rotor_performance.jl
@@ -348,10 +348,8 @@ In-place version of `get_total_efficiency`.
 """
 function get_total_efficiency!(eta, total_thrust, total_power, Vinf)
     for i in 1:length(total_thrust)
-        # if Vinf <= 0.0 || total_power[i] < eps() || total_thrust[i] <= 0.0
-        #do nothing, efficiency can't physically be negative or infinite.
-        if abs(total_power[i]) < eps()
-            # allow negative efficiency
+        if Vinf <= 0.0 || total_power[i] < eps() || total_thrust[i] <= 0.0
+            #do nothing, efficiency can't physically be negative or infinite.
             eta[i] = 0.0
         else
             eta[i] = total_thrust[i] * Vinf / total_power[i]
@@ -388,8 +386,7 @@ In-place version of `get_induced_efficiency`.
 """
 function get_induced_efficiency!(eta_inv, Tinv, Tduct, Pinv, Vinf)
     for (e, ti, p) in zip(eachrow(eta_inv), Tinv, Pinv)
-        # if Vinf <= 0.0 || p <= 0.0
-        if abs(p) <= eps()
+        if Vinf <= 0.0 || p <= eps() || (ti + Tduct) <= 0.0
             e[1] = 0.0
         else
             e[1] = Vinf * (ti + Tduct) / p

diff --git a/src/postprocess/viscous_drag.jl b/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 * FLOWMath.abs_smooth(Ue / Uinf, 2.0*eps())^((5.0 + H12) / 2.0)
+    return 2.0 * d2 * FLOWMath.abs_smooth(Ue / Uinf, 2.0 * eps())^((5.0 + H12) / 2.0)
 end
 
 """

diff --git a/src/process/process.jl b/src/process/process.jl
@@ -127,24 +127,6 @@ function process(
     )
 
     # - Initialize States - #
-    # DFDC-like Initialization: In optimization, tends to break due to sqrt of negative numbers every now and then.
-    # initialize_strengths!(
-    #     solver_options,
-    #     Gamr,
-    #     sigr,
-    #     gamw,
-    #     solve_containers,
-    #     operating_point,
-    #     (; blade_elements..., airfoils...),
-    #     wakeK,
-    #     idmaps.wake_nodemap,
-    #     idmaps.wake_endnodeidxs,
-    #     idmaps.wake_panel_sheet_be_map,
-    #     idmaps.wake_node_sheet_be_map,
-    #     idmaps.wake_node_ids_along_casing_wake_interface,
-    #     idmaps.wake_node_ids_along_centerbody_wake_interface;
-    # )
-
     initialize_strengths!(
         solver_options,
         Gamr,
@@ -187,9 +169,7 @@ function process(
     if options.verbose
         println("\nSolving Nonlinear System using CSOR Method")
     end
-    # return ImplicitAD.implicit(
-    #     solve, CSOR_residual_iad!, solve_parameter_cache_vector, const_cache
-    # )
+
     return solve(solver_options, solve_parameter_cache_vector, const_cache)
 end
 

diff --git a/src/process/solve.jl b/src/process/solve.jl
@@ -162,6 +162,8 @@ function solve(
         solver_options,
         # Constant Parameters
         solve_parameter_cache_dims,
+        # Cache(s)
+        solve_container_cache,
     ) = const_cache
 
     # - Extract Initial Guess Vector for State Variables - #
@@ -175,6 +177,7 @@ function solve(
     if verbose
         println("  " * "Wrapping Residual")
     end
+
     residual_wrapper(r, states) = mod_CSOR_residual!(r, states, inputs, const_cache)
 
     # - Get number of blades for use in relaxation - #

diff --git a/src/utilities/misc.jl b/src/utilities/misc.jl
@@ -38,12 +38,12 @@ Formatted printing when you have lots of debugging to do and want things to line
 """
 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]
+        s = @sprintf "%*s %-*s %2.18f" 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
+            s = @sprintf "%*s %-*s %2.18f" (nspaces) " " colw "" v
             println(rstrip(s, ['0', ' ']))
         end
     end

diff --git a/src/utilities/options.jl b/src/utilities/options.jl
@@ -250,7 +250,7 @@ Type containing all the options for the modified CSOR solver.
 
 # Fields
 - `verbose::Bool = false` : flag to print verbose statements
-- `iteration_limit::Float = 1e3` : maximum number of iterations
+- `iteration_limit::Float = 250` : maximum number of iterations
 - `relaxation_parameters::NamedTuple` = (;
     - `nrf::Float = 0.4` : nominal relaxation factor
     - `bt1::Float = 0.2` : backtracking factor 1
@@ -266,7 +266,7 @@ Type containing all the options for the modified CSOR solver.
 @kwdef struct ModCSORSolverOptions{TB,TF,TI,TT} <: InternalSolverOptions
     # Defaults are DFDC hard-coded values
     verbose::TB = false
-    iteration_limit::TI = 1000
+    iteration_limit::TI = 250
     relaxation_parameters::TT = (;
         nrf=0.4, bt1=0.2, bt2=0.6, pf1=0.4, pf2=0.5, btw=0.6, pfw=1.2
     )

diff --git a/src/visualization/convenience_plots.jl b/src/visualization/convenience_plots.jl
@@ -217,12 +217,18 @@ function generate_plots(
     plot_boundary_layer=false,
     plot_streamlines=false,
     verbose=false,
+    titles=nothing,
+    fps=30,
     kwargs...,
 )
 
     # Extract useful items
     (; body_vortex_panels, wake_vortex_panels, rotor_source_panels) = ins.panels
 
+    if isnothing(titles)
+        titles = fill(nothing, length(outs))
+    end
+
     # - Geometry - #
     if plot_geometry
         verbose && println("Plotting Geometry")
@@ -244,11 +250,15 @@ function generate_plots(
         if verbose
             pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true)
         end
-        for out in outs
+        for (o, out) in enumerate(outs)
 
             # underlay geometry first
             plt = Plots.plot(
-                underlayGeometry(), body_vortex_panels, rotor_source_panels; kwargs...
+                underlayGeometry(),
+                body_vortex_panels,
+                rotor_source_panels;
+                title=titles[o],
+                kwargs...,
             )
 
             # then Plots.plot pressure distributions
@@ -267,7 +277,7 @@ function generate_plots(
             verbose && next!(pbar)
         end
         verbose && println("Saving $(save_path)surface_pressure.gif")
-        Plots.gif(anim, save_path * "surface_pressure.gif")
+        Plots.gif(anim, save_path * "surface_pressure.gif"; fps=fps)
     end
 
     # - Tangential Velocity - #
@@ -277,10 +287,14 @@ function generate_plots(
         if verbose
             pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true)
         end
-        for out in outs
+        for (o, out) in enumerate(outs)
             # underlay geometry first
             plt = Plots.plot(
-                underlayGeometry(), body_vortex_panels, rotor_source_panels; kwargs...
+                underlayGeometry(),
+                body_vortex_panels,
+                rotor_source_panels;
+                title=titles[o],
+                kwargs...,
             )
 
             # then Plots.plot velocity distributions
@@ -299,7 +313,7 @@ function generate_plots(
             verbose && next!(pbar)
         end
         verbose && println("Saving $(save_path)surface_velocity.gif")
-        Plots.gif(anim, save_path * "surface_velocity.gif")
+        Plots.gif(anim, save_path * "surface_velocity.gif"; fps=fps)
     end
 
     # - Boundary Layer Stuff - #
@@ -309,10 +323,15 @@ function generate_plots(
         if verbose
             pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true)
         end
-        for out in outs
+        for (o, out) in enumerate(outs)
             # Plots.plot momentum thicknesses on top
             plt = Plots.plot(
-                plotDuctGeometry(), body_vortex_panels; color=6, linewidth=0.5, kwargs...
+                plotDuctGeometry(),
+                body_vortex_panels;
+                color=6,
+                linewidth=0.5,
+                title=titles[o],
+                kwargs...,
             )
 
             Plots.plot!(
@@ -336,7 +355,7 @@ function generate_plots(
             verbose && next!(pbar)
         end
         verbose && println("Saving $(save_path)boundary_layer.gif")
-        Plots.gif(anim, save_path * "boundary_layer.gif")
+        Plots.gif(anim, save_path * "boundary_layer.gif"; fps=fps)
     end
 
     # - Streamlines - #
@@ -346,8 +365,10 @@ function generate_plots(
             pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true)
         end
         anim = Plots.Animation()
-        for out in outs
-            plt = Plots.plot(plotBodyGeometry(), body_vortex_panels; kwargs...)
+        for (o, out) in enumerate(outs)
+            plt = Plots.plot(
+                plotBodyGeometry(), body_vortex_panels; title=titles[o], kwargs...
+            )
 
             Plots.plot!(
                 plotStreamlines(),
@@ -382,7 +403,7 @@ function generate_plots(
             verbose && next!(pbar)
         end
         verbose && println("Saving $(save_path)streamlines.gif")
-        Plots.gif(anim, save_path * "streamlines.gif")
+        Plots.gif(anim, save_path * "streamlines.gif"; fps=fps)
     end
 
     return nothing

diff --git a/src/visualization/plot_recipes.jl b/src/visualization/plot_recipes.jl
@@ -591,47 +591,6 @@ end
         ylim --> bl_ylim
     end
 
-    @series begin
-        color --> 1
-        label --> false
-        # linewidth --> 0.5
-        if isapprox(blo.separation_point_ratio_upper, 1.0)
-            arrow --> :closed
-        end
-
-        # spline the states vs steps
-        d2_upper = blo.upper_solved_states[1, :]
-
-        # spline the control points vs surface length
-        cpzu = akima_smooth(
-            blo.surface_length_upper,
-            [stagz; bvp.controlpoint[1, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
-            blo.upper_solved_steps,
-        )
-        cpru = akima_smooth(
-            blo.surface_length_upper,
-            [stagr; bvp.controlpoint[2, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
-            blo.upper_solved_steps,
-        )
-
-        nhatzu = akima_smooth(
-            blo.surface_length_upper,
-            [stagnz; bvp.normal[1, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
-            blo.upper_solved_steps,
-        )
-        nhatru = akima_smooth(
-            blo.surface_length_upper,
-            [stagnr; bvp.normal[2, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
-            blo.upper_solved_steps,
-        )
-
-        # thickness point is thickness normal to control points
-        d2zu = cpzu .+ (d2_upper .* nhatzu) .* scale_thickness
-        d2ru = cpru .+ (d2_upper .* nhatru) .* scale_thickness
-
-        return d2zu, d2ru
-    end
-
     @series begin
         color --> 2
         label --> false
@@ -673,6 +632,49 @@ end
         return d2zl, d2rl
     end
 
+    if blo.split_ratio < 1.0
+        @series begin
+            color --> 1
+            label --> false
+            # linewidth --> 0.5
+            if isapprox(blo.separation_point_ratio_upper, 1.0)
+                arrow --> :closed
+            end
+
+            # spline the states vs steps
+            d2_upper = blo.upper_solved_states[1, :]
+
+            # spline the control points vs surface length
+            cpzu = akima_smooth(
+                blo.surface_length_upper,
+                [stagz; bvp.controlpoint[1, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
+                blo.upper_solved_steps,
+            )
+            cpru = akima_smooth(
+                blo.surface_length_upper,
+                [stagr; bvp.controlpoint[2, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
+                blo.upper_solved_steps,
+            )
+
+            nhatzu = akima_smooth(
+                blo.surface_length_upper,
+                [stagnz; bvp.normal[1, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
+                blo.upper_solved_steps,
+            )
+            nhatru = akima_smooth(
+                blo.surface_length_upper,
+                [stagnr; bvp.normal[2, blo.stagnation_indices[2]:Int(bvp.npanel[1])]],
+                blo.upper_solved_steps,
+            )
+
+            # thickness point is thickness normal to control points
+            d2zu = cpzu .+ (d2_upper .* nhatzu) .* scale_thickness
+            d2ru = cpru .+ (d2_upper .* nhatru) .* scale_thickness
+
+            return d2zu, d2ru
+        end
+    end
+
     return nothing
 end