Add nonlinear solver

Author: 1-Bart-1

Project.toml

@@ -19,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"
@@ -55,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"

examples/ram_air_kite.jl

@@ -25,7 +25,8 @@ end
 vsm_solver = Solver(body_aero;
     aerodynamic_model_type=VSM,
     is_with_artificial_damping=false,
-    allowed_error=1e-8
+    allowed_error=1e-8,
+    solver_type=NONLIN
 )
 
 # Setting velocity conditions

src/VortexStepMethod.jl

@@ -19,6 +19,7 @@ using PreallocationTools
 using PrecompileTools
 using Pkg
 using DifferentiationInterface
+import SciMLBase: succesful_retcode
 
 # Export public interface
 export Wing, Section, RamAirWing
@@ -33,6 +34,7 @@ export AeroModel, LEI_AIRFOIL_BREUKELS, POLAR_VECTORS, POLAR_MATRICES, INVISCID
 export PanelDistribution, LINEAR, COSINE, COSINE_VAN_GARREL, SPLIT_PROVIDED, UNCHANGED
 export InitialGammaDistribution, ELLIPTIC, ZEROS
 export SolverStatus, FEASIBLE, INFEASIBLE, FAILURE
+export SolverType, LOOP, NONLIN
 
 export plot_geometry, plot_distribution, plot_circulation_distribution, plot_geometry, plot_polars, save_plot, show_plot, plot_polar_data
 

src/solver.jl

@@ -402,7 +402,7 @@ function solve_base!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribu
     # Run main iteration loop
     gamma_loop!(solver, body_aero, panels, relaxation_factor; log)
     # Try again with reduced relaxation factor if not converged
-    if ! solver.lr.converged && relaxation_factor > 1e-3
+    if solver.solver_type == LOOP && !solver.lr.converged && relaxation_factor > 1e-3
         log && @warn "Running again with half the relaxation_factor = $(relaxation_factor/2)"
         solver.lr.gamma_new .= gamma_initial
         gamma_loop!(solver, body_aero, panels, relaxation_factor/2; log)
@@ -453,9 +453,7 @@ function gamma_loop!(
     velocity_view_y = @view induced_velocity_all[:, 2]
     velocity_view_z = @view induced_velocity_all[:, 3]
 
-    function calc_gamma_new!(gamma_new, gamma)
-        gamma .= gamma_new
-        
+    function calc_gamma_new!(gamma_new, gamma)     
         # Calculate induced velocities
         mul!(velocity_view_x, AIC_x, gamma)
         mul!(velocity_view_y, AIC_y, gamma)
@@ -493,15 +491,16 @@ function gamma_loop!(
             d_gamma .= solver.lr.gamma_new .- gamma
             nothing
         end
-        d_gamma = abs_gamma_new
-        prob = NonlinearProblem(f_nonlin!, d_gamma, solver.lr.gamma_new, p)
-        sol = solve(prob, NewtonRaphson())
+        prob = NonlinearProblem(f_nonlin!, solver.lr.gamma_new, nothing)
+        sol = NonlinearSolve.solve(prob, NewtonRaphson(autodiff=AutoFiniteDiff()))
         gamma .= sol.u
         solver.lr.gamma_new .= sol.u
+        solver.lr.converged = succesful_retcode(sol)
     end
 
     if solver.solver_type == LOOP
         function f_loop!(gamma_new, gamma, damp)
+            gamma .= gamma_new
             calc_gamma_new!(gamma_new, gamma)
             # Update gamma with relaxation and damping
             @. gamma_new = (1 - relaxation_factor) * gamma + 

test/bench.jl

@@ -140,7 +140,7 @@ using LinearAlgebra
                     0.5;
                     log = false
                 ) samples=1 evals=1
-                @test result.allocs ≤ 20
+                @test result.allocs ≤ 10
                 @info "Model: $model \t Aero_model: $aero_model \t Allocations: $(result.allocs) Memory: $(result.memory)"
             end
         end

test/runtests.jl

@@ -21,7 +21,6 @@ println("Running tests...")
 @testset verbose = true "Testing VortexStepMethod..." begin
     if build_is_production_build
         include("bench.jl")
-        include("bench2.jl")
     end
     include("test_bound_filament.jl")
     include("test_panel.jl")