Add nonlin test

Author: 1-Bart-1

examples/ram_air_kite.jl

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

src/VortexStepMethod.jl

@@ -19,7 +19,7 @@ using PreallocationTools
 using PrecompileTools
 using Pkg
 using DifferentiationInterface
-import SciMLBase: succesful_retcode
+import SciMLBase: successful_retcode
 
 # Export public interface
 export Wing, Section, RamAirWing

src/solver.jl

@@ -73,7 +73,7 @@ Main solver structure for the Vortex Step Method.See also: [solve](@ref)
 - `aerodynamic_model_type`::Model = VSM: The model type, see: [Model](@ref)
 - density::Float64 = 1.225: Air density [kg/m³] 
 - `max_iterations`::Int64 = 1500
-- `allowed_error`::Float64 = 1e-5: relative error
+- `rtol`::Float64 = 1e-5: relative error
 - `tol_reference_error`::Float64 = 0.001
 - `relaxation_factor`::Float64 = 0.03: Relaxation factor for convergence 
 
@@ -96,7 +96,8 @@ sol::VSMSolution = VSMSolution(): The result of calling [solve!](@ref)
     aerodynamic_model_type::Model = VSM
     density::Float64 = 1.225
     max_iterations::Int64 = 1500
-    allowed_error::Float64 = 1e-5
+    atol::Float64 = 1e-5
+    rtol::Float64 = 1e-5
     tol_reference_error::Float64 = 0.001
     relaxation_factor::Float64 = 0.03
 
@@ -492,10 +493,10 @@ function gamma_loop!(
             nothing
         end
         prob = NonlinearProblem(f_nonlin!, solver.lr.gamma_new, nothing)
-        sol = NonlinearSolve.solve(prob, NewtonRaphson(autodiff=AutoFiniteDiff()))
+        sol = NonlinearSolve.solve(prob, NewtonRaphson(autodiff=AutoFiniteDiff()); abstol=solver.atol, reltol=solver.rtol)
         gamma .= sol.u
         solver.lr.gamma_new .= sol.u
-        solver.lr.converged = succesful_retcode(sol)
+        solver.lr.converged = SciMLBase.successful_retcode(sol)
     end
 
     if solver.solver_type == LOOP
@@ -532,7 +533,7 @@ function gamma_loop!(
 
             @debug "Iteration: $i, normalized_error: $normalized_error, is_damping_applied: $is_damping_applied"
 
-            if normalized_error < solver.allowed_error
+            if normalized_error < solver.rtol
                 solver.lr.converged = true
                 break
             end

test/test_body_aerodynamics.jl

@@ -146,34 +146,49 @@ end
 
     # Run analysis
     P = length(body_aero.panels)
-    solver_object = Solver{P}(
+    loop_solver = Solver{P}(
         aerodynamic_model_type=model, 
-        core_radius_fraction=core_radius_fraction
+        core_radius_fraction=core_radius_fraction,
+        solver_type=LOOP,
+        atol=1e-8,
+        rtol=1e-8
     )
-    results_NEW = solve(solver_object, body_aero; reference_point=[0,1,0])
+    nonlin_solver = Solver{P}(
+        aerodynamic_model_type=model, 
+        core_radius_fraction=core_radius_fraction,
+        solver_type=NONLIN,
+        atol=1e-8,
+        rtol=1e-8
+    )
+    results_NEW = solve(loop_solver, body_aero; reference_point=[0,1,0])
     # println(results_NEW)
 
     @test results_NEW isa Dict
 
-    sol = solve!(solver_object, body_aero; reference_point=[0,1,0])
+    @testset "Loop and nonlin solve!" begin
+        loop_sol = solve!(loop_solver, body_aero; reference_point=[0,1,0])
+        nonlin_sol = solve!(nonlin_solver, body_aero; reference_point=[0,1,0])
 
-    @test sol.force.x ≈ -117.97225244011436 atol=1e-4
-    @test sol.force.y ≈ 0.0 atol=1e-10
-    @test sol.force.z ≈ 1481.996390329679 atol=1e-4 rtol= 1e-4
+        @test all(isapprox.(nonlin_sol.gamma_distribution, loop_sol.gamma_distribution; atol=1e-4))
 
-    @test sol.moment.x ≈ -1481.996390329678 atol=1e-4 rtol= 1e-4
-    @test sol.moment.y ≈ 0.0 atol=1e-10
-    @test sol.moment.z ≈ -117.97225244011435 atol=1e-4
+        @test loop_sol.force.x ≈ -117.96518414816444 atol=1e-4
+        @test loop_sol.force.y ≈ 0.0 atol=1e-10
+        @test loop_sol.force.z ≈ 1481.996390329679 atol=1e-4 rtol= 1e-4
 
-    @test sol.force_coefficients[1] ≈ -0.039050322560956294 atol=1e-4 # CFx
-    @test sol.force_coefficients[2] ≈ 0.0                   atol=1e-4 # CFy
-    @test sol.force_coefficients[3] ≈ 0.49055973654418716   atol=1e-4 # CFz
-    @test sol.force_coefficients[3] / sol.force_coefficients[1] ≈ sol.force[3] / sol.force[1]
-    @test sol.moment_distribution[1] ≈ -0.0006683746751654071 atol=1e-8
-    @test sol.moment_coeff_dist[1] ≈ -2.212405554436003e-7 atol=1e-10
-    @test sol.moment_distribution[1] / sol.moment_distribution[2] ≈ sol.moment_coeff_dist[1] / sol.moment_coeff_dist[2]
+        @test loop_sol.moment.x ≈ -1481.996390329678 atol=1e-4 rtol= 1e-4
+        @test loop_sol.moment.y ≈ 0.0 atol=1e-10
+        @test loop_sol.moment.z ≈ -117.9651841481644 atol=1e-4
 
-    @test sol.solver_status == FEASIBLE
+        @test loop_sol.force_coefficients[1] ≈ -0.039050322560956294 atol=1e-4 # CFx
+        @test loop_sol.force_coefficients[2] ≈ 0.0                   atol=1e-4 # CFy
+        @test loop_sol.force_coefficients[3] ≈ 0.49055973654418716   atol=1e-4 # CFz
+        @test loop_sol.force_coefficients[3] / loop_sol.force_coefficients[1] ≈ loop_sol.force[3] / loop_sol.force[1]
+        @test loop_sol.moment_distribution[1] ≈ -0.0006683569356186426 atol=1e-8
+        @test loop_sol.moment_coeff_dist[1] ≈ -2.212405554436003e-7 atol=1e-10
+        @test loop_sol.moment_distribution[1] / loop_sol.moment_distribution[2] ≈ loop_sol.moment_coeff_dist[1] / loop_sol.moment_coeff_dist[2]
+
+        @test loop_sol.solver_status == FEASIBLE
+    end
 
     # Calculate forces using uncorrected alpha
     alpha = results_NEW["alpha_uncorrected"]