only allow H to increase slightly above 3

src/postprocess/boundary_layer_head.jl

@@ -69,7 +69,8 @@ function calculate_H(H1)
 
     # get each side of the piecewise equation
     hgeq = 0.86 * (H1 - 3.3)^(-0.777) + 1.1
-    hlt = 1.1538 * (H1 - 3.3)^(-0.326) + 0.6778
+    # hlt = 1.1538 * (H1 - 3.3)^(-0.326) + 0.6778
+    hlt = FLOWMath.ksmin([3.1; 1.1538 * (H1 - 3.3)^(-0.326) + 0.6778], 5)
 
     # blend the pieces smoothly
     return FLOWMath.sigmoid_blend(hlt, hgeq, H1, 5.3)
@@ -349,7 +350,7 @@ function solve_head_boundary_layer!(
 
         #spline H and steps from solution, get step at H=3
         sepwrap(x) = parameters.separation_criteria - akima_smooth(stepsol, Hsol, x)
-        hid = findfirst(x->x>=parameters.separation_criteria, Hsol)
+        hid = findfirst(x -> x >= parameters.separation_criteria, Hsol)
         s_sep = Roots.find_zero(sepwrap, stepsol[hid])
         # spline states and steps, get states at step for H=3
         usep = [akima_smooth(stepsol, u, s_sep) for u in eachrow(usol)]

src/postprocess/boundary_layer_utils.jl

@@ -64,9 +64,10 @@ function split_at_stagnation_point(
             s_upper = nothing
             s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[end:-1:1])
             split_ratio = 1.0
+            stag_point = s_tot[1]
             stag_ids = ones(Int, 2) .* length(duct_panel_lengths)
 
-            return s_upper, s_lower, stag_ids, split_ratio, dots
+            return s_upper, s_lower, stag_ids, stag_point, split_ratio, dots
         else
             stag_point = Roots.find_zero(
                 x -> FLOWMath.derivative(vtsp, x),

src/postprocess/viscous_drag.jl

@@ -379,6 +379,8 @@ function compute_viscous_drag_duct(
         )
     else
         cdc_upper = 0.0
+        u_init_upper = 0.0
+        boundary_layer_functions_upper = (;)
         usol_upper = -ones(eltype(Vtan_duct), 3)
         stepsol_upper = -1
         s_sep_upper = -1.0