@@ -2,6 +2,7 @@ using VortexStepMethod
22using DifferentiationInterface
33using LinearAlgebra
44using Test
5+ using ForwardDiff
56
67relative_error (jac, reference) = maximum (abs .(jac .- reference)) / maximum (abs, reference)
78
@@ -97,6 +98,45 @@ relative_error(jac, reference) = maximum(abs.(jac .- reference)) / maximum(abs,
9798
9899 @info " POLAR_MATRICES jacobian norms" norm_fwd= norm (jac_fwd) norm_fd= norm (jac_fd)
99100 @test relative_error (jac_fd, jac_fwd) < 1e-4
100- include (" diag_forwarddiff.jl" )
101+ println (" DIAG cpu=" , Sys. cpu_info ()[1 ]. model, " threads=" , Threads. nthreads (),
102+ " blas=" , LinearAlgebra. BLAS. get_num_threads (), " os=" , Sys. KERNEL)
103+ println (" DIAG table_hash=" , hash ([s. aero_data for s in ram_wing. unrefined_sections]))
104+ colerr = [maximum (abs .(jac_fd[:, j] .- jac_fwd[:, j])) for j in axes (jac_fd, 2 )]
105+ println (" DIAG colerr=" , colerr)
106+ println (" DIAG argmax=" , Tuple (argmax (abs .(jac_fd .- jac_fwd))))
107+ println (" DIAG fwd_minus_fd=" , repr (jac_fwd .- jac_fd))
108+ for rep in 1 : 3
109+ jac_again, _, _ = VortexStepMethod. linearize (ram_solver, ram_body, y_op;
110+ theta_idxs= 1 : 4 , va_vec_idxs= 5 : 7 , omega_idxs= 8 : 10 ,
111+ aero_coeffs= true , backend= AutoForwardDiff ())
112+ println (" DIAG fwd_repeat$rep identical=" , jac_again == jac_fwd,
113+ " maxdiff=" , maximum (abs .(jac_again .- jac_fwd)))
114+ end
115+ N = length (y_op)
116+ TD = ForwardDiff. Dual{Nothing, Float64, N}
117+ y_d = [TD (y_op[i], ForwardDiff. Partials (ntuple (j -> Float64 (i == j), N))) for i in 1 : N]
118+ body_d, solver_d = VortexStepMethod. make_dual_shadow (ram_solver, ram_body, TD)
119+ VortexStepMethod. unrefined_deform! (body_d, y_d[1 : 4 ], nothing )
120+ VortexStepMethod. reinit! (body_d; init_aero= false )
121+ set_va! (body_d, y_d[5 : 7 ], y_d[8 : 10 ])
122+ solve! (solver_d, body_d; log= true )
123+ body_f = BodyAerodynamics ([ram_wing])
124+ VortexStepMethod. unrefined_deform! (body_f, y_op[1 : 4 ], nothing )
125+ VortexStepMethod. reinit! (body_f; init_aero= false )
126+ set_va! (body_f, y_op[5 : 7 ], y_op[8 : 10 ])
127+ solve! (ram_solver, body_f; log= true )
128+ alpha_f = collect (ram_solver. lr. alpha_dist)
129+ alpha_d = ForwardDiff. value .(collect (solver_d. lr. alpha_dist))
130+ println (" DIAG alpha_f_deg=" , repr (rad2deg .(alpha_f)))
131+ println (" DIAG alpha_dual_minus_f=" , repr (alpha_d .- alpha_f))
132+ println (" DIAG gamma_dual_minus_f=" , repr (ForwardDiff. value .(collect (solver_d. sol. gamma_distribution)) .- collect (ram_solver. sol. gamma_distribution)))
133+ println (" DIAG dalpha_dy=" , repr ([collect (ForwardDiff. partials (a)) for a in solver_d. lr. alpha_dist]))
134+ println (" DIAG dual_coeffs_partials=" , repr (collect (ForwardDiff. partials (solver_d. sol. force_coeffs[1 ]))))
135+ body_d2, solver_d2 = VortexStepMethod. make_dual_shadow (ram_solver, ram_body, TD)
136+ VortexStepMethod. unrefined_deform! (body_d2, y_d[1 : 4 ], nothing )
137+ VortexStepMethod. reinit! (body_d2; init_aero= false )
138+ set_va! (body_d2, y_d[5 : 7 ], y_d[8 : 10 ])
139+ solve! (solver_d2, body_d2)
140+ println (" DIAG second_dual_identical=" , collect (solver_d2. sol. force_coeffs) == collect (solver_d. sol. force_coeffs))
101141 end
102142end
0 commit comments