Improve jacobian calculation

Author: 1-Bart-1

src/kite_geometry.jl

@@ -522,14 +522,15 @@ function RamAirWing(obj_path, dat_path; alpha=0.0, crease_frac=0.75, wind_vel=10
 end
 
 """
-    smooth_deform!(wing::RamAirWing, theta_angles::AbstractVector, delta_angles::AbstractVector)
+    group_deform!(wing::RamAirWing, theta_angles::AbstractVector, delta_angles::AbstractVector)
 
 Distribute control angles across wing panels and apply smoothing using a moving average filter.
 
 # Arguments
 - `wing::RamAirWing`: The wing to deform
 - `theta_angles::AbstractVector`: Twist angles in radians for each control section
 - `delta_angles::AbstractVector`: Trailing edge deflection angles in radians for each control section
+- `smooth::Bool`: Wether to apply smoothing or not
 
 # Algorithm
 1. Distributes each control input to its corresponding group of panels
@@ -541,7 +542,7 @@ Distribute control angles across wing panels and apply smoothing using a moving
 # Returns
 - `nothing` (modifies wing in-place)
 """
-function smooth_deform!(wing::RamAirWing, theta_angles::AbstractVector, delta_angles::AbstractVector)
+function group_deform!(wing::RamAirWing, theta_angles::AbstractVector, delta_angles=nothing; smooth=false)
     !(wing.n_panels % length(theta_angles) == 0) && throw(ArgumentError("Number of angles has to be a multiple of number of panels"))
     n_panels = wing.n_panels
     theta_dist = wing.theta_dist
@@ -552,22 +553,29 @@ function smooth_deform!(wing::RamAirWing, theta_angles::AbstractVector, delta_an
         for _ in 1:(wing.n_panels ÷ length(theta_angles))
             dist_idx += 1
             theta_dist[dist_idx] = theta_angles[angle_idx]
-            delta_dist[dist_idx] = delta_angles[angle_idx]
+            !isnothing(delta_angles) && (delta_dist[dist_idx] = delta_angles[angle_idx])
         end
     end
     @assert (dist_idx == wing.n_panels)
 
-    window_size = wing.n_panels ÷ length(theta_angles)
-    if n_panels > window_size
-        smoothed = wing.cache[1][theta_dist]
-        for i in (window_size÷2 + 1):(n_panels - window_size÷2)
-            @views smoothed[i] = mean(theta_dist[(i - window_size÷2):(i + window_size÷2)])
-        end
-        theta_dist .= smoothed
-        for i in (window_size÷2 + 1):(n_panels - window_size÷2)
-            @views smoothed[i] = mean(delta_dist[(i - window_size÷2):(i + window_size÷2)])
+    if smooth
+        window_size = wing.n_panels ÷ length(theta_angles)
+        if n_panels > window_size
+            smoothed = wing.cache[1][theta_dist]
+            smoothed .= theta_dist
+            for i in (window_size÷2 + 1):(n_panels - window_size÷2)
+                @views smoothed[i] = mean(theta_dist[(i - window_size÷2):(i + window_size÷2)])
+            end
+            theta_dist .= smoothed
+            
+            if !isnothing(delta_angles)
+                smoothed .= delta_dist
+                for i in (window_size÷2 + 1):(n_panels - window_size÷2)
+                    @views smoothed[i] = mean(delta_dist[(i - window_size÷2):(i + window_size÷2)])
+                end
+                delta_dist .= smoothed
+            end
         end
-        delta_dist .= smoothed
     end
     deform!(wing)
     return nothing

src/panel.jl

@@ -113,7 +113,7 @@ function init_pos!(
     panel.x_airf .= x_airf
     panel.y_airf .= y_airf
     panel.z_airf .= z_airf
-    panel.delta = Float64(delta)
+    panel.delta = delta
     return nothing
 end
 
@@ -361,7 +361,7 @@ function calculate_cd_cm(panel::Panel, alpha::Float64)
     elseif panel.aero_model == POLAR_VECTORS
         cd = panel.cd_interp(alpha)::Float64
         cm = panel.cm_interp(alpha)::Float64
-    elseif panel.aero_model == POLAR_MATRICES    
+    elseif panel.aero_model == POLAR_MATRICES
         cd = panel.cd_interp(alpha, panel.delta)::Float64
         cm = panel.cm_interp(alpha, panel.delta)::Float64
     elseif !(panel.aero_model == INVISCID)

src/solver.jl

@@ -110,7 +110,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [solve!](@ref)
     artificial_damping::NamedTuple{(:k2, :k4), Tuple{Float64, Float64}} =(k2=0.1, k4=0.0)
 
     # Additional settings
-    type_initial_gamma_distribution::InitialGammaDistribution = ELLIPTIC
+    type_initial_gamma_distribution::InitialGammaDistribution = ZEROS
     core_radius_fraction::Float64 = 1e-20
     mu::Float64 = 1.81e-5
     is_only_f_and_gamma_output::Bool = false
@@ -121,6 +121,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [solve!](@ref)
     cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:11]
     cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}}  = [LazyBufferCache()]
     cache_solve::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache()]
+    cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4]
 
     # Solution
     sol::VSMSolution{P} = VSMSolution(P)
@@ -624,26 +625,43 @@ Linearize the ram air wing aero model around an operating point using automatic
 - `forces_op`: Forces at the operating point [Fx, Fy, Fz, Mx, My, Mz]
 """
 function linearize(solver::Solver, body_aero::BodyAerodynamics, wing::RamAirWing, y; 
-        alpha_idxs=1:4, 
+        theta_idxs=1:4, 
         delta_idxs=nothing,
         va_idxs=nothing,
         omega_idxs=nothing,
         kwargs...)
 
-    init_va = copy(body_aero.va)
-    init_alpha_beta = copy(y[alpha_idxs])
-    zero_delta = zeros(alpha_idxs)
+    init_va           = solver.cache_lin[1][body_aero.va]
+    init_va .= body_aero.va
+    if !isnothing(theta_idxs)
+        last_theta = solver.cache_lin[2][y[theta_idxs]]
+        last_theta .= y[theta_idxs]
+    end
+    if !isnothing(delta_idxs)
+        last_delta = solver.cache_lin[3][y[delta_idxs]]
+        last_delta .= y[delta_idxs]
+    end
 
     # Function to compute forces for given control inputs
     function calc_results!(results, y)
-        if !isnothing(alpha_idxs) && isnothing(delta_idxs)
-            @views VortexStepMethod.smooth_deform!(wing, y[alpha_idxs], zero_delta)
-        elseif !isnothing(alpha_idxs) && !isnothing(delta_idxs)
-            @views VortexStepMethod.smooth_deform!(wing, y[alpha_idxs], y[delta_idxs])
+        @views theta_angles = isnothing(theta_idxs) ? nothing : y[theta_idxs]
+        @views delta_angles = isnothing(delta_idxs) ? nothing : y[delta_idxs]
+
+        if !isnothing(theta_angles)
+            if isnothing(delta_angles)
+                if !all(theta_angles .== last_theta)
+                    VortexStepMethod.group_deform!(wing, theta_angles; smooth=false)
+                    VortexStepMethod.init!(body_aero; init_aero=false)
+                    last_theta .= theta_angles
+                end
+            elseif !all(delta_angles .== last_delta) || !all(theta_angles .== last_theta)
+                VortexStepMethod.group_deform!(wing, theta_angles, delta_angles; smooth=false)
+                VortexStepMethod.init!(body_aero; init_aero=false)
+                last_theta .= theta_angles
+                last_delta .= delta_angles
+            end
         end
 
-        VortexStepMethod.init!(body_aero; init_aero=false)
-
         if !isnothing(va_idxs) && isnothing(omega_idxs)
             set_va!(body_aero, y[va_idxs])
         elseif !isnothing(va_idxs) && !isnothing(omega_idxs)
@@ -661,7 +679,7 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, wing::RamAirWing
 
     results = zeros(3+3+length(solver.sol.moment_coeff_dist))
     jac = zeros(length(results), length(y))
-    backend = AutoFiniteDiff(absstep=1e-3, relstep=1e-3)
+    backend = AutoFiniteDiff(absstep=1e2solver.atol, relstep=1e2solver.rtol)
     prep = prepare_jacobian(calc_results!, results, backend, y)
     jacobian!(calc_results!, results, jac, prep, backend, y)