Improve linearize test (#151)

Author: 1-Bart-1

src/kite_geometry.jl

@@ -570,31 +570,40 @@ Distribute control angles across wing panels and apply smoothing using a moving
 # Returns
 - `nothing` (modifies wing in-place)
 """
-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"))
+function group_deform!(wing::RamAirWing, theta_angles=nothing, delta_angles=nothing; smooth=false)
+    !isnothing(theta_angles) && !(wing.n_panels % length(theta_angles) == 0) && 
+        throw(ArgumentError("Number of angles has to be a multiple of number of panels"))
+    !isnothing(delta_angles) && !(wing.n_panels % length(delta_angles) == 0) && 
+        throw(ArgumentError("Number of angles has to be a multiple of number of panels"))
+    isnothing(theta_angles) && isnothing(delta_angles) && return nothing
+
     n_panels = wing.n_panels
     theta_dist = wing.theta_dist
     delta_dist = wing.delta_dist
+    n_angles = isnothing(theta_angles) ? length(delta_angles) : length(theta_angles)
 
     dist_idx = 0
-    for angle_idx in eachindex(theta_angles)
-        for _ in 1:(wing.n_panels ÷ length(theta_angles))
+    for angle_idx in 1:n_angles
+        for _ in 1:(wing.n_panels ÷ n_angles)
             dist_idx += 1
-            theta_dist[dist_idx] = theta_angles[angle_idx]
+            !isnothing(theta_angles) && (theta_dist[dist_idx] = theta_angles[angle_idx])
             !isnothing(delta_angles) && (delta_dist[dist_idx] = delta_angles[angle_idx])
         end
     end
     @assert (dist_idx == wing.n_panels)
 
     if smooth
-        window_size = wing.n_panels ÷ length(theta_angles)
+        window_size = wing.n_panels ÷ n_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)])
+
+            if !isnothing(theta_angles)
+                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
             end
-            theta_dist .= smoothed
 
             if !isnothing(delta_angles)
                 smoothed .= delta_dist

src/solver.jl

@@ -703,19 +703,25 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
         @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)
+        if !isnothing(theta_angles) && isnothing(delta_angles)
+            if !all(theta_angles .== last_theta)
+                VortexStepMethod.group_deform!(wing, theta_angles, nothing; smooth=false)
+                VortexStepMethod.init!(body_aero; init_aero=false)
+                last_theta .= theta_angles
+            end
+        elseif !isnothing(theta_angles) && !isnothing(delta_angles)
+            if !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
+        elseif isnothing(theta_angles) && !isnothing(delta_angles)
+            if !all(delta_angles .== last_delta)
+                VortexStepMethod.group_deform!(wing, nothing, delta_angles; smooth=false)
+                VortexStepMethod.init!(body_aero; init_aero=false)
+                last_delta .= delta_angles
+            end
         end
 
         if !isnothing(va_idxs) && isnothing(omega_idxs)

test/test_results.jl

@@ -132,49 +132,125 @@ end
 
     # Test combinations of input variations
     @testset "Combined Input Variations" begin
-        # Sample some key combinations
+        # For combination testing, we'll create targeted test cases
+        # that use only the specific indices we want to test together
         combination_tests = [
-            ("Theta + VA", [dtheta_magnitudes[1] * ones(4); dva_magnitudes[1] * ones(3); zeros(3); zeros(4)]),
-            ("Theta + Omega", [dtheta_magnitudes[1] * ones(4); zeros(3); domega_magnitudes[1] * ones(3); zeros(4)]),
-            ("VA + Omega", [zeros(4); dva_magnitudes[1] * ones(3); domega_magnitudes[1] * ones(3); zeros(4)]),
-            ("Delta + VA", [zeros(4); dva_magnitudes[1] * ones(3); zeros(3); ddelta_magnitudes[1] * ones(4)]),
-            ("All Inputs", [dtheta_magnitudes[1] * ones(4); dva_magnitudes[1] * ones(3); 
-                           domega_magnitudes[1] * ones(3); ddelta_magnitudes[1] * ones(4)])
+            # Name, active indices, perturbation vector, indices mappings
+            (
+                "Theta + VA", 
+                # Use only theta and va indices
+                [1:4; 5:7],  
+                # Perturbation values for these indices
+                [dtheta_magnitudes; dva_magnitudes],
+                # Mappings for linearize function
+                (theta_idxs=1:4, va_idxs=5:7, omega_idxs=nothing, delta_idxs=nothing)
+            ),
+            (
+                "Theta + Omega", 
+                [1:4; 8:10],
+                [dtheta_magnitudes; domega_magnitudes],
+                (theta_idxs=1:4, va_idxs=nothing, omega_idxs=5:7, delta_idxs=nothing)
+            ),
+            (
+                "VA + Omega", 
+                [5:7; 8:10],
+                [dva_magnitudes; domega_magnitudes],
+                (theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6, delta_idxs=nothing)
+            ),
+            (
+                "Delta + VA", 
+                [5:7; 11:14],
+                [dva_magnitudes; ddelta_magnitudes],
+                (theta_idxs=nothing, va_idxs=1:3, omega_idxs=nothing, delta_idxs=4:7)
+            ),
+            (
+                "All Inputs", 
+                [1:4; 5:7; 8:10; 11:14],
+                [dtheta_magnitudes; dva_magnitudes; 
+                domega_magnitudes; ddelta_magnitudes],
+                (theta_idxs=1:4, va_idxs=5:7, omega_idxs=8:10, delta_idxs=11:14)
+            )
         ]
 
-        for (combo_name, perturbation) in combination_tests
-            # Reset to baseline
-            VortexStepMethod.group_deform!(ram_wing, theta, delta; smooth=false)
-            init!(body_aero; init_aero=false, va=va, omega=zeros(3))
-            
-            # Extract components
-            perturbed_theta = theta + perturbation[1:4]
-            perturbed_va = va + perturbation[5:7]
-            perturbed_omega = perturbation[8:10]
-            perturbed_delta = delta + perturbation[11:14]
-            
-            # Apply to nonlinear model
-            VortexStepMethod.group_deform!(ram_wing, perturbed_theta, perturbed_delta; smooth=false)
-            init!(body_aero; init_aero=false, va=perturbed_va, omega=perturbed_omega)
-            
-            # Get nonlinear solution
-            nonlin_res = VortexStepMethod.solve!(solver, body_aero; log=false)
-            nonlin_res = [solver.sol.force; solver.sol.moment; solver.sol.group_moment_dist]
-            
-            # Compute linearized prediction
-            lin_prediction = lin_res + jac * perturbation
-            
-            # Calculate error ratio
-            prediction_error = norm(lin_prediction - nonlin_res)
-            baseline_difference = norm(lin_res - nonlin_res)
-            error_ratio = prediction_error / baseline_difference
-            
-            @info "$combo_name error ratio: $error_ratio"
-            
-            # Test combinations
-            @test lin_res ≉ nonlin_res atol=1e-2
-            @test lin_prediction ≈ nonlin_res rtol=0.2 atol=0.05
-            @test error_ratio < 2e-3
+        for (combo_name, active_indices, perturbation, idx_mappings) in combination_tests
+            @testset "$combo_name" begin
+                # Start with a fresh model for each combination test
+                VortexStepMethod.group_deform!(ram_wing, theta, delta; smooth=false)
+                init!(body_aero; init_aero=false, va, omega)
+                
+                # Create the appropriate input vector for this combination
+                input_vec = Vector{Float64}(undef, length(active_indices))
+                
+                # Fill with base values
+                if !isnothing(idx_mappings.theta_idxs)
+                    input_vec[idx_mappings.theta_idxs] .= theta
+                end
+                if !isnothing(idx_mappings.va_idxs)
+                    input_vec[idx_mappings.va_idxs] .= va
+                end
+                if !isnothing(idx_mappings.omega_idxs)
+                    input_vec[idx_mappings.omega_idxs] .= omega
+                end
+                if !isnothing(idx_mappings.delta_idxs)
+                    input_vec[idx_mappings.delta_idxs] .= delta
+                end
+                
+                # Get the Jacobian and linearization result for this specific combination
+                jac_combo, lin_res_combo = VortexStepMethod.linearize(
+                    solver, 
+                    body_aero, 
+                    input_vec; 
+                    idx_mappings...
+                )
+                
+                # Get baseline results
+                baseline_res = VortexStepMethod.solve!(solver, body_aero; log=false)
+                baseline_res = [solver.sol.force; solver.sol.moment; solver.sol.group_moment_dist]
+                
+                # Should match the linearization result
+                @test baseline_res ≈ lin_res_combo
+                
+                # Apply perturbation using the appropriate indices
+                perturbed_input = copy(input_vec) + perturbation
+                
+                # Extract components based on the combination being tested
+                perturbed_theta = !isnothing(idx_mappings.theta_idxs) ? 
+                    perturbed_input[idx_mappings.theta_idxs] : theta
+                
+                perturbed_va = !isnothing(idx_mappings.va_idxs) ? 
+                    perturbed_input[idx_mappings.va_idxs] : va
+                
+                perturbed_omega = !isnothing(idx_mappings.omega_idxs) ? 
+                    perturbed_input[idx_mappings.omega_idxs] : omega
+                
+                perturbed_delta = !isnothing(idx_mappings.delta_idxs) ? 
+                    perturbed_input[idx_mappings.delta_idxs] : delta
+                
+                # Apply to nonlinear model
+                VortexStepMethod.group_deform!(ram_wing, perturbed_theta, perturbed_delta; smooth=false)
+                init!(body_aero; init_aero=false, va=perturbed_va, omega=perturbed_omega)
+                
+                # Get nonlinear solution with perturbation
+                nonlin_res = VortexStepMethod.solve!(solver, body_aero; log=false)
+                nonlin_res = [solver.sol.force; solver.sol.moment; solver.sol.group_moment_dist]
+                
+                # Compute linearized prediction using our specialized Jacobian
+                lin_prediction = lin_res_combo + jac_combo * perturbation
+                
+                # Ensure perturbation had an effect
+                @test baseline_res ≉ nonlin_res atol=1e-3
+                
+                # Calculate error ratio
+                prediction_error = norm(lin_prediction - nonlin_res)
+                baseline_difference = norm(baseline_res - nonlin_res)
+                error_ratio = prediction_error / baseline_difference
+                
+                @info "$combo_name error metrics" prediction_error baseline_difference error_ratio
+                
+                # Validate the prediction
+                @test lin_prediction ≈ nonlin_res rtol=0.1 atol=1e-3
+                @test error_ratio < 0.05
+            end
         end
     end
 end