diff --git a/src/postprocess/boundary_layer_utils.jl b/src/postprocess/boundary_layer_utils.jl index a7cc56f3..67512628 100644 --- 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 index 1762f738..72b615a9 100644 --- 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 index b35db081..66b3d365 100644 --- 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 index c8f644cb..e30564cc 100644 --- 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 index 644cdf16..8d46d122 100644 --- 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 index ce470bdd..249b5d03 100644 --- 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 index 2dfef3ac..ceda89c4 100644 --- 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 index a7891cf7..7048355d 100644 --- 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 index 5e1e5eff..00421f31 100644 --- 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 index 40d129b9..6d0668f9 100644 --- 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