LOOP converges on the relaxed step, so rtol is 33x looser than it reads

CHANGELOG.md

@@ -24,6 +24,12 @@
   `test/solver/solver_test_wing.yaml` wing at 26.6° it stopped 3.6% below `LOOP`'s
   peak circulation and reported `FAILURE`, and now lands on the same distribution
   with a fixed-point residual at machine precision.
+- The `LOOP` solver tests convergence on the fixed-point residual instead of on
+  the under-relaxed step, so `rtol` is the tolerance it reads rather than
+  `rtol / relaxation_factor`. Every `LOOP` solve is now tighter by that factor —
+  33x at the defaults — which moves coefficients in the last few digits and costs
+  around 1.4x the iterations, and a solve that misses the tolerances is no longer
+  retried at a tolerance its first attempt would have passed.
 - `plot_slices_3d` in audit mode no longer crashes when a generated deflection
   `.dat` holds no finite coordinates (an all-`NaN` contour from a deflection the
   2D solver converged at no angle): it skips that overlay and warns, naming the

src/settings.jl

@@ -154,7 +154,7 @@ Solver configuration, used within [`VSMSettings`](@ref).
 - `solver_type`: `"LOOP"` or `"NONLIN"` (default `"LOOP"`)
 - `density`: Air density (kg/m^3) (default `1.225`)
 - `max_iterations`: Maximum solver iterations (default `1500`)
-- `rtol`: Relative tolerance (default `1e-5`)
+- `rtol`: Relative tolerance on the fixed-point residual (default `1e-5`)
 - `tol_reference_error`: Reference error tolerance
     (default `0.001`)
 - `relaxation_factor`: Convergence relaxation factor
@@ -187,7 +187,7 @@ Solver configuration, used within [`VSMSettings`](@ref).
     solver_type::String = "LOOP"    # type of solver
     density::Float64 = 1.225                # air density  [kg/m³] 
     max_iterations::Int64 = 1500
-    rtol::Float64 = 1e-5                    # relative error   [-]
+    rtol::Float64 = 1e-5                    # relative residual tolerance [-]
     tol_reference_error::Float64 = 0.001
     relaxation_factor::Float64 = 0.03       # relaxation factor for convergence
     artificial_damping::Bool = false        # whether to apply artificial damping

src/solver.jl

@@ -115,7 +115,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
-- `rtol`::Float64 = 1e-5: relative error
+- `rtol`::Float64 = 1e-5: Relative tolerance on the fixed-point residual
 - `tol_reference_error`::Float64 = 0.001
 - `relaxation_factor`::Float64 = 0.03: Relaxation factor for convergence 
 
@@ -861,9 +861,10 @@ end
 
 Main iteration loop for calculating circulation distribution.
 
-The NONLIN solver is a Newton iteration on the fixed-point residual
-`F(gamma) - gamma` with a finite-difference Jacobian, backtracking along each step
-until it reduces the residual.
+Both solvers converge on the fixed-point residual `F(gamma) - gamma`, measured
+relative to the largest circulation: the LOOP solver takes under-relaxed steps
+towards `F(gamma)`, the NONLIN solver a Newton step with a finite-difference
+Jacobian, backtracked until it reduces the residual.
 
 When `solver.is_with_artificial_viscosity` is set, the LOOP solver replaces the
 explicit target `F(gamma)` with the implicit Li/Gaunaa solution
@@ -1085,9 +1086,10 @@ function gamma_loop!(
             abs_gamma_new .= abs.(solver.lr.gamma_new)
             reference_error = maximum(abs_gamma_new)
             reference_error = max(reference_error, solver.tol_reference_error)
+            # The relaxed step is `relaxation_factor` times the fixed-point residual.
             abs_gamma_new .= abs.(solver.lr.gamma_new .- gamma)
-            error = maximum(abs_gamma_new)
-            normalized_error = error / reference_error
+            residual = maximum(abs_gamma_new) / relaxation_factor
+            normalized_error = residual / reference_error
 
             @debug "Iteration: $i, normalized_error: $normalized_error"
 

test/solver/test_solver.jl

@@ -101,6 +101,45 @@ end
     end
 end
 
+"""
+    unrelaxed_step(body_aero, gamma)
+
+One unrelaxed fixed-point step `F(gamma)` of the LOOP iteration, so that
+`F(gamma) - gamma` is the residual of `gamma`.
+"""
+function unrelaxed_step(body_aero, gamma)
+    # An infinite rtol accepts the single step, so solve_base! skips its retry.
+    probe = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
+        relaxation_factor=1.0, max_iterations=1, rtol=Inf)
+    VortexStepMethod.solve_base!(probe, body_aero, gamma)
+    return copy(probe.lr.gamma_new)
+end
+
+@testset "LOOP converges on the residual, not on the relaxed step" begin
+    settings_file = create_temp_wing_settings(
+        "solver", "solver_test_wing.yaml";
+        alpha=5.0, beta=0.0, wind_speed=10.0,
+    )
+    try
+        settings = VSMSettings(settings_file)
+        wing = Wing(settings)
+        refine!(wing)
+        body_aero = BodyAerodynamics([wing])
+        solver = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
+            type_initial_gamma_distribution=ELLIPTIC)
+
+        for va in ([10.0, 0.0, 0.0], [10.0, 0.0, 5.0])   # 0 deg, and 26.6 deg past stall
+            set_va!(body_aero, va)
+            gamma = copy(solve!(solver, body_aero).gamma_distribution)
+            @test solver.lr.converged
+            residual = maximum(abs, unrelaxed_step(body_aero, gamma) .- gamma)
+            @test residual < solver.rtol * maximum(abs, gamma)
+        end
+    finally
+        rm(settings_file; force=true)
+    end
+end
+
 calc_forces_allocs(solver, body_aero) =
     (calc_forces!(solver, body_aero); @allocated calc_forces!(solver, body_aero))