|
| 1 | +using Pkg |
| 2 | +if Base.active_project() != joinpath(@__DIR__, "Project.toml") |
| 3 | + Pkg.activate(@__DIR__) |
| 4 | +end |
| 5 | + |
| 6 | +using GLMakie |
| 7 | +using DifferentiationInterface |
| 8 | +using LinearAlgebra |
| 9 | +using VortexStepMethod |
| 10 | +using VortexStepMethod: linearize, unrefined_deform!, reinit! |
| 11 | + |
| 12 | +# Sweep each linearize input around the operating point and overlay the |
| 13 | +# FiniteDiff and ForwardDiff tangents on the sweep curve. |
| 14 | + |
| 15 | +n_unrefined = 4 |
| 16 | + |
| 17 | +wing = ObjWing( |
| 18 | + joinpath("data", "ram_air_kite", "ram_air_kite_body.obj"), |
| 19 | + joinpath("data", "ram_air_kite", "ram_air_kite_foil.dat"); |
| 20 | + n_unrefined_sections=n_unrefined, |
| 21 | + prn=false, |
| 22 | +) |
| 23 | +body_aero = BodyAerodynamics([wing]) |
| 24 | + |
| 25 | +solver = Solver(body_aero; |
| 26 | + aerodynamic_model_type=VSM, |
| 27 | + is_with_artificial_damping=false, |
| 28 | + rtol=1e-7, |
| 29 | + solver_type=LOOP, |
| 30 | + use_gamma_prev=false, |
| 31 | +) |
| 32 | + |
| 33 | +v_a = 15.0 |
| 34 | +aoa_deg = 10.0 |
| 35 | +aoa_rad = deg2rad(aoa_deg) |
| 36 | +side_slip = 0.0 |
| 37 | +va_b_0 = [ |
| 38 | + cos(aoa_rad) * cos(side_slip), |
| 39 | + sin(side_slip), |
| 40 | + sin(aoa_rad), |
| 41 | +] * v_a |
| 42 | +omega_b_0 = zeros(3) |
| 43 | +theta_0 = zeros(n_unrefined) |
| 44 | + |
| 45 | +theta_idxs = 1:n_unrefined |
| 46 | +va_idxs = (n_unrefined + 1):(n_unrefined + 3) |
| 47 | +omega_idxs = (n_unrefined + 4):(n_unrefined + 6) |
| 48 | +y0 = [theta_0; va_b_0; omega_b_0] |
| 49 | + |
| 50 | +@info "Computing FiniteDiff Jacobian …" |
| 51 | +t_fd = @elapsed begin |
| 52 | + jac_fd, x0_fd, conv_fd = linearize( |
| 53 | + solver, body_aero, y0; |
| 54 | + theta_idxs, va_idxs, omega_idxs, |
| 55 | + aero_coeffs=true, |
| 56 | + backend=AutoFiniteDiff(absstep=1e-5, relstep=1e-5), |
| 57 | + ) |
| 58 | +end |
| 59 | +conv_fd || @warn "FiniteDiff linearize did not converge at operating point" |
| 60 | + |
| 61 | +@info "Computing ForwardDiff Jacobian …" |
| 62 | +t_fwd = @elapsed begin |
| 63 | + jac_fwd, x0_fwd, conv_fwd = linearize( |
| 64 | + solver, body_aero, y0; |
| 65 | + theta_idxs, va_idxs, omega_idxs, |
| 66 | + aero_coeffs=true, |
| 67 | + backend=AutoForwardDiff(), |
| 68 | + ) |
| 69 | +end |
| 70 | +conv_fwd || @warn "ForwardDiff linearize did not converge at operating point" |
| 71 | + |
| 72 | +x0 = x0_fd |
| 73 | + |
| 74 | +@info "Operating point" CFx=x0[1] CFy=x0[2] CFz=x0[3] CM=x0[4:6] |
| 75 | +@info "Timing" t_finitediff_s=t_fd t_forwarddiff_s=t_fwd speedup=t_fd / t_fwd |
| 76 | + |
| 77 | +denom = max(maximum(abs, jac_fwd), eps()) |
| 78 | +rel_err_fwd_fd = maximum(abs.(jac_fwd .- jac_fd)) / denom |
| 79 | +@info "AutoForwardDiff vs AutoFiniteDiff" rel_err=rel_err_fwd_fd |
| 80 | + |
| 81 | +input_labels = [ |
| 82 | + ["θ_$i" for i in 1:n_unrefined]..., |
| 83 | + "va_x", "va_y", "va_z", |
| 84 | + "ω_x", "ω_y", "ω_z", |
| 85 | +] |
| 86 | +output_labels = [ |
| 87 | + "CFx", "CFy", "CFz", |
| 88 | + "CMx", "CMy", "CMz", |
| 89 | + ["cm_$i" for i in 1:n_unrefined]..., |
| 90 | +] |
| 91 | +n_inputs = length(input_labels) |
| 92 | +n_outputs = length(output_labels) |
| 93 | + |
| 94 | +@assert size(jac_fd) == (n_outputs, n_inputs) |
| 95 | +@assert size(jac_fwd) == (n_outputs, n_inputs) |
| 96 | + |
| 97 | +input_scales = [ |
| 98 | + fill(0.05, n_unrefined)..., # θ [rad] : ±0.05 rad ≈ ±2.9° |
| 99 | + fill(1.0, 3)..., # va: ±1 m/s |
| 100 | + fill(0.05, 3)..., # ω : ±0.05 rad/s |
| 101 | +] |
| 102 | +n_sweep = 11 |
| 103 | +sweep_frac = range(-1.0, 1.0; length=n_sweep) |
| 104 | + |
| 105 | +results_sw = [zeros(n_sweep, n_outputs) for _ in 1:n_inputs] |
| 106 | + |
| 107 | +last_theta = fill(NaN, n_unrefined) |
| 108 | + |
| 109 | +function solve_at!(y) |
| 110 | + theta = y[theta_idxs] |
| 111 | + va = y[va_idxs] |
| 112 | + omega = y[omega_idxs] |
| 113 | + if !all(theta .== last_theta) |
| 114 | + unrefined_deform!(wing, theta, nothing; smooth=false) |
| 115 | + reinit!(body_aero; init_aero=false) |
| 116 | + last_theta .= theta |
| 117 | + end |
| 118 | + set_va!(body_aero, va, omega) |
| 119 | + solve!(solver, body_aero; log=false) |
| 120 | + return [ |
| 121 | + solver.sol.force_coeffs..., |
| 122 | + solver.sol.moment_coeffs..., |
| 123 | + solver.sol.cm_unrefined_dist..., |
| 124 | + ] |
| 125 | +end |
| 126 | + |
| 127 | +@info "Sweeping each input …" |
| 128 | +for ci in 1:n_inputs |
| 129 | + @info " $(input_labels[ci])" |
| 130 | + for (si, frac) in enumerate(sweep_frac) |
| 131 | + y = copy(y0) |
| 132 | + y[ci] += frac * input_scales[ci] |
| 133 | + results_sw[ci][si, :] .= solve_at!(y) |
| 134 | + end |
| 135 | +end |
| 136 | + |
| 137 | +solve_at!(y0) |
| 138 | + |
| 139 | +fig = Figure(size=(180 * n_inputs + 80, 90 * n_outputs + 100)) |
| 140 | +Label(fig[0, 1:n_inputs], |
| 141 | + "linearize check: sweeps (blue) vs FiniteDiff tangent (red dashed) " * |
| 142 | + "vs ForwardDiff tangent (green dashed) — " * |
| 143 | + "FD: $(round(t_fd; digits=2)) s, FwdDiff: $(round(t_fwd; digits=2)) s, " * |
| 144 | + "speedup: $(round(t_fd / t_fwd; digits=2))×"; |
| 145 | + fontsize=16, font=:bold, tellwidth=false) |
| 146 | + |
| 147 | +for ri in 1:n_outputs |
| 148 | + for ci in 1:n_inputs |
| 149 | + ax = Axis(fig[ri, ci]; |
| 150 | + xlabel = ri == n_outputs ? input_labels[ci] : "", |
| 151 | + ylabel = ci == 1 ? output_labels[ri] : "", |
| 152 | + xticklabelsvisible = ri == n_outputs, |
| 153 | + yticklabelsvisible = ci == 1, |
| 154 | + xticksvisible = ri == n_outputs, |
| 155 | + yticksvisible = ci == 1, |
| 156 | + ) |
| 157 | + delta_input = sweep_frac .* input_scales[ci] |
| 158 | + ys_sweep = results_sw[ci][:, ri] |
| 159 | + ys_linear_fd = x0[ri] .+ jac_fd[ri, ci] .* delta_input |
| 160 | + ys_linear_fwd = x0[ri] .+ jac_fwd[ri, ci] .* delta_input |
| 161 | + |
| 162 | + lines!(ax, delta_input, ys_sweep; |
| 163 | + color=:steelblue, linewidth=1.5) |
| 164 | + lines!(ax, delta_input, ys_linear_fd; |
| 165 | + color=:crimson, linestyle=:dash, linewidth=1.2) |
| 166 | + lines!(ax, delta_input, ys_linear_fwd; |
| 167 | + color=:seagreen, linestyle=:dash, linewidth=1.2) |
| 168 | + scatter!(ax, [0.0], [x0[ri]]; |
| 169 | + color=:black, markersize=6) |
| 170 | + end |
| 171 | +end |
| 172 | + |
| 173 | +colgap!(fig.layout, 6) |
| 174 | +rowgap!(fig.layout, 4) |
| 175 | + |
| 176 | +display(fig) |
| 177 | + |
| 178 | +# Skip near-zero entries — relative error there is sweep rounding noise. |
| 179 | +function worst_jac_vs_sweep(jac, label) |
| 180 | + sig_threshold = 1e-3 * maximum(abs, jac) |
| 181 | + max_rel = 0.0 |
| 182 | + worst = (0, 0, 0.0, 0.0) |
| 183 | + for ri in 1:n_outputs, ci in 1:n_inputs |
| 184 | + delta = sweep_frac .* input_scales[ci] |
| 185 | + mid = div(n_sweep, 2) + 1 |
| 186 | + slope = (results_sw[ci][mid + 1, ri] - |
| 187 | + results_sw[ci][mid - 1, ri]) / |
| 188 | + (delta[mid + 1] - delta[mid - 1]) |
| 189 | + jc = jac[ri, ci] |
| 190 | + max(abs(jc), abs(slope)) < sig_threshold && continue |
| 191 | + rel = abs(slope - jc) / max(abs(jc), abs(slope)) |
| 192 | + if rel > max_rel |
| 193 | + max_rel = rel |
| 194 | + worst = (ri, ci, slope, jc) |
| 195 | + end |
| 196 | + end |
| 197 | + @info "Worst $label vs sweep mismatch (significant entries)" rel=max_rel output=output_labels[worst[1]] input=input_labels[worst[2]] sweep_slope=worst[3] jac_entry=worst[4] sig_threshold |
| 198 | +end |
| 199 | + |
| 200 | +worst_jac_vs_sweep(jac_fd, "FiniteDiff") |
| 201 | +worst_jac_vs_sweep(jac_fwd, "ForwardDiff") |
| 202 | + |
| 203 | +nothing |
0 commit comments