Add linearization feature and nonlinear solve (#140)

* Add smooth deform function

* Reduce allocs

* Add linearize function and docs

* Working linearize example

* Add nonlin solver which added allocations

* Move bench2 test to bench

* Fix allocations

* Add nonlinear solver

* Add nonlin test

* Remove a lot of unnecessary function signature types

* Use a lot more mvec

* Nonlin solve for jacobian

* Move cache inside struct

* Fix example

* Fix docstring

* Fix docstring

* Add types back to struct

* Non-allocating smooth deform

* Improve jacobian calculation

* Less allocs

* Add result tests

* Rename wa to body_aero

* Fix unchanged benchmark

* Add docs to linearize

* Fix create interp bug

* Add linearization tests

* Fix test for windows

* Fix tests for windows

* Remove unused P

* Add back types

* Remove function arg type

* Back to old types

Author: 1-Bart-1

Project.toml

@@ -7,6 +7,8 @@ version = "1.1.2"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DefaultApplication = "3f0dd361-4fe0-5fc6-8523-80b14ec94d85"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -17,6 +19,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -39,8 +42,10 @@ ControlPlots = "0.2.5"
 DataFrames = "1.7"
 DefaultApplication = "1"
 DelimitedFiles = "1"
+DifferentiationInterface = "0.6.48"
 Distributed = "1.10"
 Documenter = "1.8"
+FiniteDiff = "2.27.0"
 Interpolations = "0.15"
 LaTeXStrings = "1"
 LinearAlgebra = "1"
@@ -51,6 +56,7 @@ Parameters = "0.12"
 Pkg = "1"
 PreallocationTools = "0.4.25"
 PrecompileTools = "1.2.1"
+SciMLBase = "2.79.0"
 Serialization = "1"
 SharedArrays = "1"
 StaticArrays = "1"

README.md

@@ -117,7 +117,7 @@ add_section!(wing,
     INVISCID)
 
 # Step 3: Initialize aerodynamics
-wa = BodyAerodynamics([wing])
+body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a

docs/src/examples.md

@@ -54,19 +54,19 @@ add_section!(wing,
 
 #### Step 4: Initialize aerodynamics
 ```
-wa = BodyAerodynamics([wing])
+body_aero = BodyAerodynamics([wing])
 ```
 We need to pass here an array of wing objects, because a body can have multiple wings.
 
 ###### Set inflow conditions
 ```julia
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(wa, vel_app, [0, 0, 0.1])
+set_va!(body_aero, vel_app, [0, 0, 0.1])
 ```
 #### Step 5: Plot the geometry
 ```julia
 plot_geometry(
-      wa,
+      body_aero,
       "Rectangular_wing_geometry";
       data_type=".pdf",
       save_path=".",
@@ -86,8 +86,8 @@ vsm_solver = Solver(aerodynamic_model_type=VSM)
 
 #### Step 7: Solve using both methods
 ```
-results_llt = solve(llt_solver, wa)
-results_vsm = solve(vsm_solver, wa)
+results_llt = solve(llt_solver, body_aero)
+results_vsm = solve(vsm_solver, body_aero)
 ```
 
 ##### Print results comparison
@@ -103,7 +103,7 @@ println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²"
 
 #### Step 8: Plot spanwise distributions
 ```julia
-y_coordinates = [panel.aero_center[2] for panel in wa.panels]
+y_coordinates = [panel.aero_center[2] for panel in body_aero.panels]
 
 plot_distribution(
     [y_coordinates, y_coordinates],
@@ -121,7 +121,7 @@ You should see a plot like this:
 angle_range = range(0, 20, 20)
 plot_polars(
     [llt_solver, vsm_solver],
-    [wa, wa],
+    [body_aero, body_aero],
     ["LLT", "VSM"];
     angle_range,
     angle_type="angle_of_attack",

docs/src/index.md

@@ -115,7 +115,7 @@ add_section!(wing,
     INVISCID)
 
 # Step 3: Initialize aerodynamics
-wa = BodyAerodynamics([wing])
+body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a

docs/src/types.md

@@ -28,7 +28,7 @@ A body is constructed of one or more abstract wings. An abstract wing can be a W
 A Wing/ RamAirWing has one or more sections.
 ```@docs
 Section
-Section(LE_point::Vector{Float64}, TE_point::Vector{Float64}, aero_model=nothing, aero_data=nothing)
+Section(LE_point::PosVector, TE_point::PosVector, aero_model)
 Wing
 Wing(n_panels::Int; spanwise_distribution::PanelDistribution=LINEAR,
      spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]))

examples/bench.jl

@@ -32,29 +32,28 @@ add_section!(wing,
     INVISCID)
 
 # Step 3: Initialize aerodynamics
-wa = BodyAerodynamics([wing])
+body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(wa, vel_app)
+set_va!(body_aero, vel_app)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
-P = length(wa.panels)
-llt_solver = Solver{P}(aerodynamic_model_type=LLT)
-vsm_solver = Solver{P}(aerodynamic_model_type=VSM)
+llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
+vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
-results_vsm = solve(vsm_solver, wa)
-sol = solve!(vsm_solver, wa)
-results_vsm_base = solve_base!(vsm_solver, wa)
+results_vsm = solve(vsm_solver, body_aero)
+sol = solve!(vsm_solver, body_aero)
+results_vsm_base = solve_base!(vsm_solver, body_aero)
 println("Rectangular wing, solve_base!:")
-@time results_vsm_base = solve_base!(vsm_solver, wa)
+@time results_vsm_base = solve_base!(vsm_solver, body_aero)
 # time Python: 32.0  ms Ryzen 7950x
 # time Julia:   0.42 ms Ryzen 7950x
 println("Rectangular wing, solve!:")
-@time sol = solve!(vsm_solver, wa)
+@time sol = solve!(vsm_solver, body_aero)
 println("Rectangular wing, solve:")
-@time solve(vsm_solver, wa)
+@time solve(vsm_solver, body_aero)
 
 # Create wing geometry
 wing = RamAirWing("data/ram_air_kite_body.obj", "data/ram_air_kite_foil.dat")
@@ -64,7 +63,8 @@ body_aero = BodyAerodynamics([wing])
 vsm_solver = Solver(
     body_aero;
     aerodynamic_model_type=VSM,
-    is_with_artificial_damping=false
+    is_with_artificial_damping=false,
+    solver_type=NONLIN,
 )
 
 # Setting velocity conditions
@@ -82,8 +82,8 @@ set_va!(body_aero, vel_app)
 
 # Solving and plotting distributions
 results = solve(vsm_solver, body_aero)
-results_base = solve_base!(vsm_solver, body_aero)
+solve_base!(vsm_solver, body_aero)
 println("RAM-air kite:")
-@time results_base = solve_base!(vsm_solver, body_aero)
+@time solve_base!(vsm_solver, body_aero)
 
 nothing

examples/ram_air_kite.jl

@@ -2,9 +2,10 @@ using ControlPlots
 using VortexStepMethod
 using LinearAlgebra
 
-PLOT = false
+PLOT = true
 USE_TEX = false
-DEFORM = true
+DEFORM = false
+LINEARIZE = false
 
 # Create wing geometry
 wing = RamAirWing("data/ram_air_kite_body.obj", "data/ram_air_kite_foil.dat")
@@ -14,23 +15,18 @@ println("First init")
 
 if DEFORM
     # Linear interpolation of alpha from 10° at one tip to 0° at the other
-    n_panels = wing.n_panels
-    theta_start = deg2rad(10)
-    theta_end = deg2rad(0)
-    delta_start = deg2rad(10)
-    delta_end = deg2rad(0)
-    theta_dist = [theta_start - i * (theta_start - theta_end)/(n_panels-1) for i in 0:(n_panels-1)]
-    delta_dist = [delta_start - i * (delta_start - delta_end)/(n_panels-1) for i in 0:(n_panels-1)]
     println("Deform")
-    @time VortexStepMethod.deform!(wing, theta_dist, delta_dist)
+    @time VortexStepMethod.smooth_deform!(wing, deg2rad.([10,20,10,0]), deg2rad.([-10,0,-10,0]))
     println("Deform init")
     @time VortexStepMethod.init!(body_aero; init_aero=false)
 end
 
 # Create solvers
-vsm_solver = Solver(body_aero;
+solver = Solver(body_aero;
     aerodynamic_model_type=VSM,
-    is_with_artificial_damping=false
+    is_with_artificial_damping=false,
+    rtol=1e-5,
+    solver_type=NONLIN
 )
 
 # Setting velocity conditions
@@ -46,6 +42,26 @@ vel_app = [
 ] * v_a
 set_va!(body_aero, vel_app)
 
+if LINEARIZE
+    println("Linearize")
+    jac, res = VortexStepMethod.linearize(
+        solver, 
+        body_aero, 
+        [zeros(4); vel_app; zeros(3)];
+        theta_idxs=1:4,
+        va_idxs=5:7,
+        omega_idxs=8:10,
+        moment_frac=0.1)
+    @time jac, res = VortexStepMethod.linearize(
+        solver, 
+        body_aero, 
+        [zeros(4); vel_app; zeros(3)]; 
+        theta_idxs=1:4, 
+        va_idxs=5:7, 
+        omega_idxs=8:10,
+        moment_frac=0.1)
+end
+
 # Plotting polar data
 PLOT && plot_polar_data(body_aero)
 
@@ -63,8 +79,9 @@ PLOT && plot_geometry(
 )
 
 # Solving and plotting distributions
-results = solve(vsm_solver, body_aero; log=true)
-@time results = solve(vsm_solver, body_aero; log=true)
+println("Solve")
+results = VortexStepMethod.solve(solver, body_aero; log=true)
+@time results = solve(solver, body_aero; log=true)
 
 body_y_coordinates = [panel.aero_center[2] for panel in body_aero.panels]
 
@@ -80,7 +97,7 @@ PLOT && plot_distribution(
 )
 
 PLOT && plot_polars(
-    [vsm_solver],
+    [solver],
     [body_aero],
     [
         "VSM from Ram Air Kite OBJ and DAT file",

examples/rectangular_wing.jl

@@ -28,22 +28,21 @@ add_section!(wing,
     INVISCID)
 
 # Step 3: Initialize aerodynamics
-wa = BodyAerodynamics([wing])
+body_aero = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(wa, vel_app, [0, 0, 0.1])
+set_va!(body_aero, vel_app, [0, 0, 0.1])
 
 # Step 4: Initialize solvers for both LLT and VSM methods
-P = length(wa.panels)
-llt_solver = Solver{P}(aerodynamic_model_type=LLT)
-vsm_solver = Solver{P}(aerodynamic_model_type=VSM)
+llt_solver = Solver(body_aero; aerodynamic_model_type=LLT)
+vsm_solver = Solver(body_aero; aerodynamic_model_type=VSM)
 
 # Step 5: Solve using both methods
-results_llt = solve(llt_solver, wa)
-@time results_llt = solve(llt_solver, wa)
-results_vsm = solve(vsm_solver, wa)
-@time results_vsm = solve(vsm_solver, wa)
+results_llt = solve(llt_solver, body_aero)
+@time results_llt = solve(llt_solver, body_aero)
+results_vsm = solve(vsm_solver, body_aero)
+@time results_vsm = solve(vsm_solver, body_aero)
 
 # Print results comparison
 println("\nLifting Line Theory Results:")
@@ -56,7 +55,7 @@ println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²"
 
 # Step 6: Plot geometry
 PLOT && plot_geometry(
-      wa,
+      body_aero,
       "Rectangular_wing_geometry";
       data_type=".pdf",
       save_path=".",
@@ -66,7 +65,7 @@ PLOT && plot_geometry(
 )
 
 # Step 7: Plot spanwise distributions
-y_coordinates = [panel.aero_center[2] for panel in wa.panels]
+y_coordinates = [panel.aero_center[2] for panel in body_aero.panels]
 
 PLOT && plot_distribution(
     [y_coordinates, y_coordinates],
@@ -80,7 +79,7 @@ PLOT && plot_distribution(
 angle_range = range(0, 20, 20)
 PLOT && plot_polars(
     [llt_solver, vsm_solver],
-    [wa, wa],
+    [body_aero, body_aero],
     ["LLT", "VSM"];
     angle_range,
     angle_type="angle_of_attack",

examples/stall_model.jl

@@ -42,15 +42,14 @@ CAD_wing = Wing(n_panels; spanwise_distribution)
 for rib in rib_list
     add_section!(CAD_wing, rib[1], rib[2], rib[3], rib[4])
 end
-body_aero_CAD_19ribs = BodyAerodynamics([CAD_wing])
+body_aero = BodyAerodynamics([CAD_wing])
 
 # Create solvers
-P = length(body_aero_CAD_19ribs.panels)
-vsm_solver = Solver{P}(
+vsm_solver = Solver(body_aero; 
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false
 )
-VSM_with_stall_correction = Solver{P}(
+VSM_with_stall_correction = Solver(body_aero; 
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=true
 )
@@ -66,11 +65,11 @@ vel_app = [
     sin(side_slip),
     sin(aoa_rad)
 ] * v_a
-set_va!(body_aero_CAD_19ribs, vel_app)
+set_va!(body_aero, vel_app)
 
 # Plotting geometry
 PLOT && plot_geometry(
-    body_aero_CAD_19ribs,
+    body_aero,
     "";
     data_type=".svg",
     save_path="",
@@ -82,11 +81,11 @@ PLOT && plot_geometry(
 )
 
 # Solving and plotting distributions
-results = solve(vsm_solver, body_aero_CAD_19ribs)
-@time results_with_stall = solve(VSM_with_stall_correction, body_aero_CAD_19ribs)
-@time results_with_stall = solve(VSM_with_stall_correction, body_aero_CAD_19ribs)
+results = solve(vsm_solver, body_aero)
+@time results_with_stall = solve(VSM_with_stall_correction, body_aero)
+@time results_with_stall = solve(VSM_with_stall_correction, body_aero)
 
-CAD_y_coordinates = [panel.aero_center[2] for panel in body_aero_CAD_19ribs.panels]
+CAD_y_coordinates = [panel.aero_center[2] for panel in body_aero.panels]
 
 PLOT && plot_distribution(
     [CAD_y_coordinates, CAD_y_coordinates],
@@ -112,7 +111,7 @@ path_cfd_lebesque = joinpath(
 
 PLOT && plot_polars(
     [vsm_solver, VSM_with_stall_correction],
-    [body_aero_CAD_19ribs, body_aero_CAD_19ribs],
+    [body_aero, body_aero],
     [
         "VSM CAD 19ribs",
         "VSM CAD 19ribs , with stall correction",