Skip to content

Commit fb1eea7

Browse files
1-Bart-11-Bart-1
authored
Merge pull request #285 from OpenSourceAWE/agent/283-nonlin-misses-the-tolerances-at-26-6-whe
2 parents 3f75701 + 47b799a commit fb1eea7

3 files changed

Lines changed: 70 additions & 18 deletions

File tree

‎CHANGELOG.md‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,14 @@
1111
The default is off: a post-stall solve that misses the tolerances still returns
1212
its `solver_status == FAILURE` solution.
1313

14+
### Fixed
15+
16+
- The `NONLIN` solver backtracks along each Newton step instead of always taking
17+
it whole, so it converges past stall where the full step used to cycle: on the
18+
`test/solver/solver_test_wing.yaml` wing at 26.6° it stopped 3.6% below `LOOP`'s
19+
peak circulation and reported `FAILURE`, and now lands on the same distribution
20+
with a fixed-point residual at machine precision.
21+
1422
## VortexStepMethod v5.0.0 2026-09-07
1523

1624
### Added

‎src/solver.jl‎

Lines changed: 32 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -861,6 +861,10 @@ end
861861
862862
Main iteration loop for calculating circulation distribution.
863863
864+
The NONLIN solver is a Newton iteration on the fixed-point residual
865+
`F(gamma) - gamma` with a finite-difference Jacobian, backtracking along each step
866+
until it reduces the residual.
867+
864868
When `solver.is_with_artificial_viscosity` is set, the LOOP solver replaces the
865869
explicit target `F(gamma)` with the implicit Li/Gaunaa solution
866870
`(I - diag(mu) L) gamma = F(gamma)` before relaxation, stabilizing post-stall
@@ -960,35 +964,45 @@ function gamma_loop!(
960964

961965
_, _, info = LinearAlgebra.LAPACK.getrf!(jac, ipiv; check=false)
962966
info == 0 || break
967+
residual_norm = maximum(abs, residual)
963968
LinearAlgebra.LAPACK.getrs!('N', jac, ipiv, residual)
964969

965-
max_step = 0.0
970+
# Past stall the full Newton step overshoots into a cycle, so take the
971+
# largest of 1, 1/2 ... 1/64 of it that brings the residual down.
972+
step_fraction = 1.0
973+
for _ in 1:7
974+
@inbounds for i in 1:n_panels
975+
gamma_perturbed[i] = gamma_iter[i] - step_fraction * residual[i]
976+
end
977+
update_gamma_candidate!(
978+
residual_perturbed, gamma_perturbed, solver, panels, n_panels,
979+
AIC_x, AIC_y, AIC_z,
980+
velocity_view_x, velocity_view_y, velocity_view_z,
981+
va_array, induced_velocity_all, relative_velocity_array,
982+
y_airf_array, relative_velocity_crossz, v_acrossz_array,
983+
z_airf_array, x_airf_array,
984+
v_normal_array, v_tangential_array,
985+
va_magw_array, cl_dist, chord_array,
986+
)
987+
@inbounds for i in 1:n_panels
988+
residual_perturbed[i] -= gamma_perturbed[i]
989+
end
990+
maximum(abs, residual_perturbed) < residual_norm && break
991+
step_fraction /= 2
992+
end
993+
994+
max_step = maximum(abs, residual)
966995
ref = solver.tol_reference_error
967996
@inbounds for i in 1:n_panels
968-
s = abs(residual[i])
969-
s > max_step && (max_step = s)
970-
gamma_iter[i] -= residual[i]
997+
gamma_iter[i] = gamma_perturbed[i]
998+
residual[i] = residual_perturbed[i]
971999
g = abs(gamma_iter[i])
9721000
g > ref && (ref = g)
9731001
end
9741002
if max_step < solver.atol + solver.rtol * ref
9751003
solver.lr.converged = true
9761004
break
9771005
end
978-
979-
update_gamma_candidate!(
980-
residual, gamma_iter, solver, panels, n_panels,
981-
AIC_x, AIC_y, AIC_z,
982-
velocity_view_x, velocity_view_y, velocity_view_z,
983-
va_array, induced_velocity_all, relative_velocity_array,
984-
y_airf_array, relative_velocity_crossz, v_acrossz_array,
985-
z_airf_array, x_airf_array,
986-
v_normal_array, v_tangential_array,
987-
va_magw_array, cl_dist, chord_array,
988-
)
989-
@inbounds for i in 1:n_panels
990-
residual[i] -= gamma_iter[i]
991-
end
9921006
end
9931007

9941008
gamma .= gamma_iter

‎test/solver/test_solver.jl‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,36 @@ end
7171
end
7272
end
7373

74+
@testset "NONLIN converges past stall, where LOOP already did" begin
75+
settings_file = create_temp_wing_settings(
76+
"solver", "solver_test_wing.yaml";
77+
alpha=5.0, beta=0.0, wind_speed=10.0,
78+
)
79+
try
80+
settings = VSMSettings(settings_file)
81+
wing = Wing(settings)
82+
refine!(wing)
83+
body_aero = BodyAerodynamics([wing])
84+
va = [10.0, 0.0, 5.0] # 26.6 deg angle of attack, past stall
85+
nonlin = Solver(body_aero; solver_type=NONLIN, aerodynamic_model_type=VSM,
86+
type_initial_gamma_distribution=ELLIPTIC)
87+
loop = Solver(body_aero; solver_type=LOOP, aerodynamic_model_type=VSM,
88+
type_initial_gamma_distribution=ELLIPTIC)
89+
90+
set_va!(body_aero, va)
91+
sol_nonlin = solve!(nonlin, body_aero)
92+
gamma_nonlin = copy(sol_nonlin.gamma_distribution)
93+
set_va!(body_aero, va)
94+
sol_loop = solve!(loop, body_aero)
95+
96+
@test sol_nonlin.solver_status == FEASIBLE
97+
@test sol_loop.solver_status == FEASIBLE
98+
@test isapprox(gamma_nonlin, sol_loop.gamma_distribution; rtol=1e-3)
99+
finally
100+
rm(settings_file; force=true)
101+
end
102+
end
103+
74104
calc_forces_allocs(solver, body_aero) =
75105
(calc_forces!(solver, body_aero); @allocated calc_forces!(solver, body_aero))
76106

0 commit comments

Comments
 (0)