From 490e859ff48492ce459a1aea71234477c38c53fe Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Fri, 11 Sep 2026 21:32:33 +0200 Subject: [PATCH] Converge LOOP on the residual, not on the relaxed step Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_0123e9orAr632QxtAMPULtSi --- CHANGELOG.md | 6 ++++++ src/settings.jl | 4 ++-- src/solver.jl | 14 ++++++++------ test/solver/test_solver.jl | 39 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 55 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2de9cb76..b039d0bb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,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. ## VortexStepMethod v5.0.0 2026-09-07 diff --git a/src/settings.jl b/src/settings.jl index 5f3b8f8a..d1ae6ef7 100644 --- a/src/settings.jl +++ b/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 diff --git a/src/solver.jl b/src/solver.jl index a0e54110..543c7c9d 100644 --- a/src/solver.jl +++ b/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" diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index c59316be..b940c69e 100644 --- a/test/solver/test_solver.jl +++ b/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))