From 8c8db2619c1aeb4357c26414be4d2c0aac6d2e57 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:18:48 +0200 Subject: [PATCH 1/6] LOOP keeps halving the relaxation factor until it converges; stop at non-finite gamma The :curved reference wing's stiffest mode needs relaxation below 0.0095; the single retry at 0.015 diverged to NaN. Co-Authored-By: Claude Opus 5.5 --- src/solver.jl | 18 ++++++++++-------- .../test_body_aerodynamics.jl | 12 ++++++++++++ 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index 871e75bb..46a27168 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -652,8 +652,9 @@ Converge the circulation distribution and leave it in `solver.lr.gamma_new`, wit turning it into forces. Fills the solver's panel arrays, builds the AIC matrices, starts from `gamma_distribution` (or an elliptical/zero guess when it is `nothing` or `solver.use_gamma_prev` is false) and iterates; a `LOOP` solver that fails to converge -retries once with half the relaxation factor. The circulation half of [`solve!`](@ref), -paired with [`calc_forces!`](@ref). Returns `nothing`. +retries with half the relaxation factor, warning each time, while the factor is above +1e-3. The circulation half of [`solve!`](@ref), paired with [`calc_forces!`](@ref). +Returns `nothing`. """ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_distribution=nothing; log=false) where {P, U, T} @@ -721,11 +722,11 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma solver.lr.gamma_new .= gamma_initial # Run main iteration loop gamma_loop!(solver, body_aero, panels, relaxation_factor; log) - # Try again with reduced relaxation factor if not converged - if solver.solver_type == LOOP && !solver.lr.converged && relaxation_factor > 1e-3 - log && @warn "Running again with half the relaxation_factor = $(relaxation_factor/2)" + while solver.solver_type == LOOP && !solver.lr.converged && relaxation_factor > 1e-3 + relaxation_factor /= 2 + @warn "LOOP did not converge, retrying with relaxation_factor = $relaxation_factor" solver.lr.gamma_new .= gamma_initial - gamma_loop!(solver, body_aero, panels, relaxation_factor/2; log) + gamma_loop!(solver, body_aero, panels, relaxation_factor; log) end nothing @@ -891,8 +892,8 @@ Main iteration loop for calculating circulation distribution. 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. +towards `F(gamma)` and stops at the first non-finite residual, 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 @@ -1096,6 +1097,7 @@ function gamma_loop!( normalized_error = residual / reference_error @debug "Iteration: $i, normalized_error: $normalized_error" + isfinite(normalized_error) || break if normalized_error < solver.rtol solver.lr.converged = true diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index 893853c4..c0be5a6b 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -292,6 +292,18 @@ end end end end + + # CL from the NONLIN solve of the same wing + @testset "VSM LOOP converges on the curved wing at $(alpha_deg)°" for (alpha_deg, CL) in + ((5.7106, 0.422308), (10.0, 0.726321)) + body_aero, _, _, _ = create_geometry(wing_type=:curved) + set_va!(body_aero, 20.0 .* [cosd(alpha_deg), 0.0, sind(alpha_deg)]) + wing = body_aero.wings[1] + solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM) + sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE + @test sol.force_coeffs[3] ≈ CL rtol=1e-4 + end end @testset "Calculate results against output results" begin From c58083755927e4957d074317dd00ba1fa3eab8a6 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:34:32 +0200 Subject: [PATCH 2/6] LOOP relaxes each panel by its self-induced stiffness; unrefined coefficients are area-weighted Narrow tip panels set the LOOP stability limit; scaling their relaxation by median/own stiffness lets the default factor converge on the curved wing and on solver_test_wing up to 40 panels. The solver_test_wing polar gains a -20 deg row so its cl crosses zero and refined tips can unload. cl/cd/cm/alpha_unrefined_dist become area-weighted means, so coefficient times section area is the section load. Co-Authored-By: Claude Opus 5.5 --- src/solver.jl | 74 ++++++++++++++++++++---------- test/solver/solver_test_polar.csv | 1 + test/solver/test_unrefined_dist.jl | 44 +++++++++--------- 3 files changed, 75 insertions(+), 44 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index 46a27168..971c029e 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -27,11 +27,11 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - `moment_dist`::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm] - `moment_coeff_dist`::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-] - `moment_unrefined_dist`::MVector{U, Float64}: Averaged moments for unrefined sections [Nm] -- `cl_unrefined_dist`::MVector{U, Float64}: Averaged lift coefficients for unrefined sections [-] -- `cd_unrefined_dist`::MVector{U, Float64}: Averaged drag coefficients for unrefined sections [-] -- `cm_unrefined_dist`::MVector{U, Float64}: Averaged airfoil moment coefficients for unrefined sections [-] +- `cl_unrefined_dist`::MVector{U, Float64}: Area-weighted mean lift coefficient of each unrefined section [-] +- `cd_unrefined_dist`::MVector{U, Float64}: Area-weighted mean drag coefficient of each unrefined section [-] +- `cm_unrefined_dist`::MVector{U, Float64}: Area-weighted mean airfoil moment coefficient of each unrefined section [-] - `moment_coeff_unrefined_dist`::MVector{U, Float64}: Summed `moment_frac`-referenced pitching-moment coefficient per unrefined section [-] -- `alpha_unrefined_dist`::MVector{U, Float64}: Averaged angles of attack for unrefined sections [rad] +- `alpha_unrefined_dist`::MVector{U, Float64}: Area-weighted mean angle of attack of each unrefined section [rad] - `solver_status`::SolverStatus: enum, see [`SolverStatus`](@ref) """ @with_kw mutable struct VSMSolution{P, U, T} @@ -181,7 +181,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) # Intermediate results lr::LoopResult{P, T} = LoopResult{P, T}() br::BaseResult{P, T} = BaseResult{P, T}() - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:8] + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:10] cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache()] cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4] @@ -498,42 +498,44 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; original_section_idx = wing.refined_panel_mapping[local_panel_idx] target_unrefined_idx = section_range[original_section_idx] + area = panel.chord * panel.width + # Accumulate coefficients and moments moment_unrefined_dist[target_unrefined_idx] += moment_dist[panel_idx] - cl_unrefined_dist[target_unrefined_idx] += solver.sol.cl_dist[panel_idx] - cd_unrefined_dist[target_unrefined_idx] += solver.sol.cd_dist[panel_idx] - cm_unrefined_dist[target_unrefined_idx] += solver.sol.cm_dist[panel_idx] + cl_unrefined_dist[target_unrefined_idx] += solver.sol.cl_dist[panel_idx] * area + cd_unrefined_dist[target_unrefined_idx] += solver.sol.cd_dist[panel_idx] * area + cm_unrefined_dist[target_unrefined_idx] += solver.sol.cm_dist[panel_idx] * area moment_coeff_unrefined_dist[target_unrefined_idx] += solver.sol.moment_coeff_dist[panel_idx] - alpha_unrefined_dist[target_unrefined_idx] += solver.sol.alpha_dist[panel_idx] + alpha_unrefined_dist[target_unrefined_idx] += solver.sol.alpha_dist[panel_idx] * area # Accumulate geometry x_airf_unrefined_dist[target_unrefined_idx] .+= panel.x_airf y_airf_unrefined_dist[target_unrefined_idx] .+= panel.y_airf z_airf_unrefined_dist[target_unrefined_idx] .+= panel.z_airf va_vec_unrefined_dist[target_unrefined_idx] .+= panel.va_vec - chord_unrefined_dist[target_unrefined_idx] += panel.chord + chord_unrefined_dist[target_unrefined_idx] += area width_unrefined_dist[target_unrefined_idx] += panel.width unrefined_count_dist[target_unrefined_idx] += 1 end - # Average coefficients and geometry. width and - # moment_coeff_unrefined_dist stay summed (extensive). + # Coefficients and alpha are area-weighted means, the chord is the + # section area over its width, width and moment_coeff stay summed. for target_unrefined_idx in section_range if unrefined_count_dist[target_unrefined_idx] > 0 count = unrefined_count_dist[target_unrefined_idx] + area = chord_unrefined_dist[target_unrefined_idx] moment_unrefined_dist[target_unrefined_idx] /= count - cl_unrefined_dist[target_unrefined_idx] /= count - cd_unrefined_dist[target_unrefined_idx] /= count - cm_unrefined_dist[target_unrefined_idx] /= count - alpha_unrefined_dist[target_unrefined_idx] /= count + cl_unrefined_dist[target_unrefined_idx] /= area + cd_unrefined_dist[target_unrefined_idx] /= area + cm_unrefined_dist[target_unrefined_idx] /= area + alpha_unrefined_dist[target_unrefined_idx] /= area x_airf_unrefined_dist[target_unrefined_idx] ./= count y_airf_unrefined_dist[target_unrefined_idx] ./= count z_airf_unrefined_dist[target_unrefined_idx] ./= count va_vec_unrefined_dist[target_unrefined_idx] ./= count - chord_unrefined_dist[target_unrefined_idx] /= count - # width_unrefined_dist is NOT averaged - it is the - # sum of panel widths in the unrefined section + chord_unrefined_dist[target_unrefined_idx] = + area / width_unrefined_dist[target_unrefined_idx] end end end @@ -883,6 +885,28 @@ function apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, visco return true end +""" + panel_relaxation_factors!(panel_relaxation, stiffness, relaxation_factor, AIC, + z_airf_dist, chord_dist) + +Fill `panel_relaxation` with the LOOP relaxation factor of each panel: `relaxation_factor`, +scaled down by `median(stiffness) / stiffness[i]` on each panel whose self-induced +stiffness `1 + π c_i |z_airf_i ⋅ AIC[i, i, :]|` exceeds the median. `stiffness` is work +space of the same length. +""" +function panel_relaxation_factors!(panel_relaxation, stiffness, relaxation_factor, AIC, + z_airf_dist, chord_dist) + @inbounds for i in eachindex(panel_relaxation) + normal_self_induction = z_airf_dist[i, 1] * AIC[i, i, 1] + + z_airf_dist[i, 2] * AIC[i, i, 2] + z_airf_dist[i, 3] * AIC[i, i, 3] + panel_relaxation[i] = 1 + π * chord_dist[i] * abs(normal_self_induction) + end + stiffness .= panel_relaxation + reference = median!(stiffness) + @. panel_relaxation = relaxation_factor * min(1, reference / panel_relaxation) + return panel_relaxation +end + """ gamma_loop!(solver::Solver, AIC_x::Matrix{Float64}, AIC_y::Matrix{Float64}, AIC_z::Matrix{Float64}, @@ -926,6 +950,8 @@ function gamma_loop!( relative_velocity_crossz = solver.cache[6][va_vec_dist] v_acrossz_dist = solver.cache[7][va_vec_dist] cl_dist = solver.cache[8][solver.lr.gamma_new] + panel_relaxation = solver.cache[9][solver.lr.gamma_new] + stiffness_work = solver.cache[10][solver.lr.gamma_new] AIC_x = @view body_aero.AIC[:, :, 1] AIC_y = @view body_aero.AIC[:, :, 2] @@ -1043,6 +1069,8 @@ function gamma_loop!( planform_area += panels[i].width * chord_dist[i] end end + panel_relaxation_factors!(panel_relaxation, stiffness_work, relaxation_factor, + body_aero.AIC, z_airf_dist, chord_dist) function f_loop!(gamma_new, gamma) gamma .= gamma_new @@ -1078,7 +1106,7 @@ function gamma_loop!( ) end - @. gamma_new = (1 - relaxation_factor) * gamma + relaxation_factor * gamma_new + @. gamma_new = (1 - panel_relaxation) * gamma + panel_relaxation * gamma_new return nothing end iters = 0 @@ -1091,9 +1119,9 @@ 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) - residual = maximum(abs_gamma_new) / relaxation_factor + # The relaxed step is `panel_relaxation` times the fixed-point residual. + abs_gamma_new .= abs.(solver.lr.gamma_new .- gamma) ./ panel_relaxation + residual = maximum(abs_gamma_new) normalized_error = residual / reference_error @debug "Iteration: $i, normalized_error: $normalized_error" diff --git a/test/solver/solver_test_polar.csv b/test/solver/solver_test_polar.csv index 7e043d9c..dd714098 100644 --- a/test/solver/solver_test_polar.csv +++ b/test/solver/solver_test_polar.csv @@ -1,4 +1,5 @@ alpha,cl,cd,cm +-20.0,-0.8,0.04,-0.02 -10.0,0.1,0.02,-0.01 0.0,1.0,0.01,0.0 10.0,1.8,0.025,0.01 diff --git a/test/solver/test_unrefined_dist.jl b/test/solver/test_unrefined_dist.jl index 9bca8a5c..d284152c 100644 --- a/test/solver/test_unrefined_dist.jl +++ b/test/solver/test_unrefined_dist.jl @@ -2,6 +2,15 @@ using VortexStepMethod using LinearAlgebra using Test +""" + area_weighted_mean(values, panels, indices) + +Mean of `values[indices]` weighted by the area of `panels[indices]`. +""" +area_weighted_mean(values, panels, indices) = + sum(values[i] * panels[i].chord * panels[i].width for i in indices) / + sum(panels[i].chord * panels[i].width for i in indices) + @testset "Unrefined Arrays Tests" begin @testset "Unrefined section array aggregation" begin # Create a simple wing with unrefined sections @@ -28,6 +37,7 @@ using Test va_vec = [10.0, 0.0, 0.0] set_va!(body_aero, va_vec) sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE # Test 1: Unrefined arrays exist and have correct size @test length(sol.cl_unrefined_dist) == wing.n_unrefined_sections @@ -38,20 +48,18 @@ using Test @test !all(sol.cl_unrefined_dist .== 0.0) @test !all(sol.cd_unrefined_dist .== 0.0) - # Test 3: Verify unrefined coefficients are averaged from refined panels + # Test 3: Unrefined coefficients are area-weighted means of their refined panels # refined_panel_mapping maps each refined panel to its unrefined section index for unrefined_idx in 1:wing.n_unrefined_sections # Find all refined panels that map to this unrefined section refined_panel_indices = findall(x -> x == unrefined_idx, wing.refined_panel_mapping) if !isempty(refined_panel_indices) - # Calculate expected average from refined panel coefficients - expected_cl = sum(sol.cl_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cd = sum(sol.cd_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cm = sum(sol.cm_dist[refined_panel_indices]) / length(refined_panel_indices) + expected_cl = area_weighted_mean(sol.cl_dist, body_aero.panels, refined_panel_indices) + expected_cd = area_weighted_mean(sol.cd_dist, body_aero.panels, refined_panel_indices) + expected_cm = area_weighted_mean(sol.cm_dist, body_aero.panels, refined_panel_indices) - # Check if unrefined coefficients match expected averages - # Handle NaN values that can occur in INVISCID models + # Handle NaN values that can occur in INVISCID models if isnan(expected_cl) @test isnan(sol.cl_unrefined_dist[unrefined_idx]) else @@ -102,6 +110,7 @@ using Test va_vec = [10.0, 0.0, 0.0] set_va!(body_aero, va_vec) sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE panels = body_aero.panels for unrefined_idx in 1:wing.n_unrefined_sections @@ -133,17 +142,9 @@ using Test panels[i].width for i in pidxs) @test w_u ≈ expected_width rtol=1e-12 - # Force-like quantities should be approximately - # equal (not exact due to averaging) - if !isnan(cl_u) - @test cl_u * c_u * w_u ≈ sum_cl_cw rtol=0.05 - end - if !isnan(cd_u) - @test cd_u * c_u * w_u ≈ sum_cd_cw rtol=0.05 - end - if !isnan(cm_u) - @test cm_u * c_u * w_u ≈ sum_cm_cw rtol=0.05 - end + @test cl_u * c_u * w_u ≈ sum_cl_cw rtol=1e-12 + @test cd_u * c_u * w_u ≈ sum_cd_cw rtol=1e-12 + @test cm_u * c_u * w_u ≈ sum_cm_cw rtol=1e-12 end finally rm(settings_file; force=true) @@ -176,6 +177,7 @@ using Test va_vec = [10.0, 0.0, 0.0] set_va!(body_aero, va_vec) sol = solve!(solver, body_aero) + @test sol.solver_status == FEASIBLE # Verify arrays have correct size @test length(sol.cl_unrefined_dist) == wing.n_unrefined_sections @@ -187,9 +189,9 @@ using Test refined_panel_indices = findall(x -> x == unrefined_idx, wing.refined_panel_mapping) if !isempty(refined_panel_indices) - expected_cl = sum(sol.cl_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cd = sum(sol.cd_dist[refined_panel_indices]) / length(refined_panel_indices) - expected_cm = sum(sol.cm_dist[refined_panel_indices]) / length(refined_panel_indices) + expected_cl = area_weighted_mean(sol.cl_dist, body_aero.panels, refined_panel_indices) + expected_cd = area_weighted_mean(sol.cd_dist, body_aero.panels, refined_panel_indices) + expected_cm = area_weighted_mean(sol.cm_dist, body_aero.panels, refined_panel_indices) # Handle NaN for all coefficients if isnan(expected_cl) From c2d3f30df994cee5564bf66ad2ac488d90b6051a Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:41:17 +0200 Subject: [PATCH 3/6] Relax each panel against the geometric mean stiffness: no sort, no allocation Co-Authored-By: Claude Opus 5.5 --- src/solver.jl | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index 971c029e..cc527dc9 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -181,7 +181,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) # Intermediate results lr::LoopResult{P, T} = LoopResult{P, T}() br::BaseResult{P, T} = BaseResult{P, T}() - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:10] + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:9] cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache()] cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4] @@ -886,23 +886,24 @@ function apply_artificial_viscosity!(gamma, panels, alpha_dist, laplacian, visco end """ - panel_relaxation_factors!(panel_relaxation, stiffness, relaxation_factor, AIC, - z_airf_dist, chord_dist) + panel_relaxation_factors!(panel_relaxation, relaxation_factor, AIC, z_airf_dist, + chord_dist) Fill `panel_relaxation` with the LOOP relaxation factor of each panel: `relaxation_factor`, -scaled down by `median(stiffness) / stiffness[i]` on each panel whose self-induced -stiffness `1 + π c_i |z_airf_i ⋅ AIC[i, i, :]|` exceeds the median. `stiffness` is work -space of the same length. +scaled by `reference / stiffness_i` on each panel whose self-induced stiffness +`stiffness_i = 1 + π c_i |z_airf_i ⋅ AIC[i, i, :]|` exceeds `reference`, the geometric +mean of the stiffness over all panels. """ -function panel_relaxation_factors!(panel_relaxation, stiffness, relaxation_factor, AIC, - z_airf_dist, chord_dist) +function panel_relaxation_factors!(panel_relaxation, relaxation_factor, AIC, z_airf_dist, + chord_dist) + log_reference = zero(eltype(panel_relaxation)) @inbounds for i in eachindex(panel_relaxation) normal_self_induction = z_airf_dist[i, 1] * AIC[i, i, 1] + z_airf_dist[i, 2] * AIC[i, i, 2] + z_airf_dist[i, 3] * AIC[i, i, 3] panel_relaxation[i] = 1 + π * chord_dist[i] * abs(normal_self_induction) + log_reference += log(panel_relaxation[i]) end - stiffness .= panel_relaxation - reference = median!(stiffness) + reference = exp(log_reference / length(panel_relaxation)) @. panel_relaxation = relaxation_factor * min(1, reference / panel_relaxation) return panel_relaxation end @@ -915,9 +916,10 @@ end Main iteration loop for calculating circulation distribution. 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)` and stops at the first non-finite residual, the NONLIN solver a -Newton step with a finite-difference Jacobian, backtracked until it reduces the residual. +relative to the largest circulation: the LOOP solver steps towards `F(gamma)`, each +panel under-relaxed by [`panel_relaxation_factors!`](@ref), and stops at the first +non-finite residual; the NONLIN solver takes 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 @@ -951,7 +953,6 @@ function gamma_loop!( v_acrossz_dist = solver.cache[7][va_vec_dist] cl_dist = solver.cache[8][solver.lr.gamma_new] panel_relaxation = solver.cache[9][solver.lr.gamma_new] - stiffness_work = solver.cache[10][solver.lr.gamma_new] AIC_x = @view body_aero.AIC[:, :, 1] AIC_y = @view body_aero.AIC[:, :, 2] @@ -1069,8 +1070,8 @@ function gamma_loop!( planform_area += panels[i].width * chord_dist[i] end end - panel_relaxation_factors!(panel_relaxation, stiffness_work, relaxation_factor, - body_aero.AIC, z_airf_dist, chord_dist) + panel_relaxation_factors!(panel_relaxation, relaxation_factor, body_aero.AIC, + z_airf_dist, chord_dist) function f_loop!(gamma_new, gamma) gamma .= gamma_new From df5f85a67cb64eb6e99b5244b455e7d9b67432d0 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:42:50 +0200 Subject: [PATCH 4/6] Cap each panel's LOOP relaxation at its diagonal step instead of a stiffness mean Leaves the relaxation of panels that were already stable unchanged: V3 and pyramid iteration counts match main. Co-Authored-By: Claude Opus 5.5 --- src/solver.jl | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index cc527dc9..fb0c8332 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -890,21 +890,17 @@ end chord_dist) Fill `panel_relaxation` with the LOOP relaxation factor of each panel: `relaxation_factor`, -scaled by `reference / stiffness_i` on each panel whose self-induced stiffness -`stiffness_i = 1 + π c_i |z_airf_i ⋅ AIC[i, i, :]|` exceeds `reference`, the geometric -mean of the stiffness over all panels. +capped at `1 / stiffness_i`, the diagonal step of the panel's self-induced stiffness +`stiffness_i = 1 + π c_i |z_airf_i ⋅ AIC[i, i, :]|`. """ function panel_relaxation_factors!(panel_relaxation, relaxation_factor, AIC, z_airf_dist, chord_dist) - log_reference = zero(eltype(panel_relaxation)) @inbounds for i in eachindex(panel_relaxation) normal_self_induction = z_airf_dist[i, 1] * AIC[i, i, 1] + z_airf_dist[i, 2] * AIC[i, i, 2] + z_airf_dist[i, 3] * AIC[i, i, 3] - panel_relaxation[i] = 1 + π * chord_dist[i] * abs(normal_self_induction) - log_reference += log(panel_relaxation[i]) + stiffness = 1 + π * chord_dist[i] * abs(normal_self_induction) + panel_relaxation[i] = min(relaxation_factor, 1 / stiffness) end - reference = exp(log_reference / length(panel_relaxation)) - @. panel_relaxation = relaxation_factor * min(1, reference / panel_relaxation) return panel_relaxation end From f7afff77a05e4ecbd4db517af096919d48a7655d Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:49:03 +0200 Subject: [PATCH 5/6] Sum moment_unrefined_dist over the section's panels; changelog and docs entry Co-Authored-By: Claude Opus 5.5 --- ...he-curved-reference-wing-does-not-conve.md | 7 ++++ docs/src/private_functions.md | 1 + src/solver.jl | 33 ++++++++++++------- .../test_body_aerodynamics.jl | 7 ++-- test/solver/test_unrefined_dist.jl | 24 ++++++++------ 5 files changed, 47 insertions(+), 25 deletions(-) create mode 100644 changelog.d/395-the-curved-reference-wing-does-not-conve.md diff --git a/changelog.d/395-the-curved-reference-wing-does-not-conve.md b/changelog.d/395-the-curved-reference-wing-does-not-conve.md new file mode 100644 index 00000000..ef5d42a6 --- /dev/null +++ b/changelog.d/395-the-curved-reference-wing-does-not-conve.md @@ -0,0 +1,7 @@ +### Fixed + +- The `LOOP` solver caps each panel's relaxation at `1 / (1 + π c |z_airf ⋅ AIC_ii|)`, the diagonal step of its self-induced stiffness, so narrow panels (cosine-spaced tips) no longer make it diverge to NaN; the curved reference wing and finely refined wings now converge at the default `relaxation_factor`, with unchanged iteration counts on wings that already converged. A `LOOP` solve that still does not converge keeps halving `relaxation_factor`, warning each time, while it is above 1e-3, and each attempt stops at the first non-finite residual. + +### Changed + +- `cl_unrefined_dist`, `cd_unrefined_dist`, `cm_unrefined_dist` and `alpha_unrefined_dist` are area-weighted means over the section's panels instead of plain means, and `chord_unrefined_dist` is the section area over its width, so a coefficient times `chord_unrefined_dist * width_unrefined_dist` is the section's load. `moment_unrefined_dist` is the section's moment, summed over its panels, instead of their mean; `linearize` with `aero_coeffs=false` returns it. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d8b0c836..53d0b9ad 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -17,6 +17,7 @@ solver_kwargs check_dimensions calculate_AIC_matrices! gamma_loop! +panel_relaxation_factors! build_spanwise_laplacian! local_lift_slope! apply_artificial_viscosity! diff --git a/src/solver.jl b/src/solver.jl index fb0c8332..8a3e8b29 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -26,12 +26,17 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - `moment_coeffs`::MVec3: Aerodynamic moment coefficients [CMx, CMy, CMz] [-] - `moment_dist`::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm] - `moment_coeff_dist`::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-] -- `moment_unrefined_dist`::MVector{U, Float64}: Averaged moments for unrefined sections [Nm] -- `cl_unrefined_dist`::MVector{U, Float64}: Area-weighted mean lift coefficient of each unrefined section [-] -- `cd_unrefined_dist`::MVector{U, Float64}: Area-weighted mean drag coefficient of each unrefined section [-] -- `cm_unrefined_dist`::MVector{U, Float64}: Area-weighted mean airfoil moment coefficient of each unrefined section [-] +- `moment_unrefined_dist`::MVector{U, Float64}: Pitching moment per unrefined section, + summed over its panels [Nm] +- `cl_unrefined_dist`::MVector{U, Float64}: Area-weighted mean lift coefficient + per unrefined section [-] +- `cd_unrefined_dist`::MVector{U, Float64}: Area-weighted mean drag coefficient + per unrefined section [-] +- `cm_unrefined_dist`::MVector{U, Float64}: Area-weighted mean airfoil moment coefficient + per unrefined section [-] - `moment_coeff_unrefined_dist`::MVector{U, Float64}: Summed `moment_frac`-referenced pitching-moment coefficient per unrefined section [-] -- `alpha_unrefined_dist`::MVector{U, Float64}: Area-weighted mean angle of attack of each unrefined section [rad] +- `alpha_unrefined_dist`::MVector{U, Float64}: Area-weighted mean angle of attack + per unrefined section [rad] - `solver_status`::SolverStatus: enum, see [`SolverStatus`](@ref) """ @with_kw mutable struct VSMSolution{P, U, T} @@ -181,7 +186,8 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) # Intermediate results lr::LoopResult{P, T} = LoopResult{P, T}() br::BaseResult{P, T} = BaseResult{P, T}() - cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:9] + cache::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = + [LazyBufferCache() for _ in 1:9] cache_base::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache()] cache_lin::Vector{PreallocationTools.LazyBufferCache{typeof(identity), typeof(identity)}} = [LazyBufferCache() for _ in 1:4] @@ -502,11 +508,15 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; # Accumulate coefficients and moments moment_unrefined_dist[target_unrefined_idx] += moment_dist[panel_idx] - cl_unrefined_dist[target_unrefined_idx] += solver.sol.cl_dist[panel_idx] * area - cd_unrefined_dist[target_unrefined_idx] += solver.sol.cd_dist[panel_idx] * area - cm_unrefined_dist[target_unrefined_idx] += solver.sol.cm_dist[panel_idx] * area + cl_unrefined_dist[target_unrefined_idx] += + solver.sol.cl_dist[panel_idx] * area + cd_unrefined_dist[target_unrefined_idx] += + solver.sol.cd_dist[panel_idx] * area + cm_unrefined_dist[target_unrefined_idx] += + solver.sol.cm_dist[panel_idx] * area moment_coeff_unrefined_dist[target_unrefined_idx] += solver.sol.moment_coeff_dist[panel_idx] - alpha_unrefined_dist[target_unrefined_idx] += solver.sol.alpha_dist[panel_idx] * area + alpha_unrefined_dist[target_unrefined_idx] += + solver.sol.alpha_dist[panel_idx] * area # Accumulate geometry x_airf_unrefined_dist[target_unrefined_idx] .+= panel.x_airf @@ -520,12 +530,11 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; end # Coefficients and alpha are area-weighted means, the chord is the - # section area over its width, width and moment_coeff stay summed. + # section area over its width; width and the moments stay summed. for target_unrefined_idx in section_range if unrefined_count_dist[target_unrefined_idx] > 0 count = unrefined_count_dist[target_unrefined_idx] area = chord_unrefined_dist[target_unrefined_idx] - moment_unrefined_dist[target_unrefined_idx] /= count cl_unrefined_dist[target_unrefined_idx] /= area cd_unrefined_dist[target_unrefined_idx] /= area cm_unrefined_dist[target_unrefined_idx] /= area diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index c0be5a6b..10d54532 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -294,12 +294,13 @@ end end # CL from the NONLIN solve of the same wing - @testset "VSM LOOP converges on the curved wing at $(alpha_deg)°" for (alpha_deg, CL) in - ((5.7106, 0.422308), (10.0, 0.726321)) + @testset "VSM LOOP converges on the curved wing at $(alpha_deg)°" for + (alpha_deg, CL) in ((5.7106, 0.422308), (10.0, 0.726321)) body_aero, _, _, _ = create_geometry(wing_type=:curved) set_va!(body_aero, 20.0 .* [cosd(alpha_deg), 0.0, sind(alpha_deg)]) wing = body_aero.wings[1] - solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM) + solver = Solver(wing.n_panels, wing.n_unrefined_sections; + aerodynamic_model_type=VSM) sol = solve!(solver, body_aero) @test sol.solver_status == FEASIBLE @test sol.force_coeffs[3] ≈ CL rtol=1e-4 diff --git a/test/solver/test_unrefined_dist.jl b/test/solver/test_unrefined_dist.jl index d284152c..3e9d3031 100644 --- a/test/solver/test_unrefined_dist.jl +++ b/test/solver/test_unrefined_dist.jl @@ -3,11 +3,11 @@ using LinearAlgebra using Test """ - area_weighted_mean(values, panels, indices) + area_mean(values, panels, indices) Mean of `values[indices]` weighted by the area of `panels[indices]`. """ -area_weighted_mean(values, panels, indices) = +area_mean(values, panels, indices) = sum(values[i] * panels[i].chord * panels[i].width for i in indices) / sum(panels[i].chord * panels[i].width for i in indices) @@ -50,16 +50,17 @@ area_weighted_mean(values, panels, indices) = # Test 3: Unrefined coefficients are area-weighted means of their refined panels # refined_panel_mapping maps each refined panel to its unrefined section index + panels = body_aero.panels for unrefined_idx in 1:wing.n_unrefined_sections # Find all refined panels that map to this unrefined section refined_panel_indices = findall(x -> x == unrefined_idx, wing.refined_panel_mapping) if !isempty(refined_panel_indices) - expected_cl = area_weighted_mean(sol.cl_dist, body_aero.panels, refined_panel_indices) - expected_cd = area_weighted_mean(sol.cd_dist, body_aero.panels, refined_panel_indices) - expected_cm = area_weighted_mean(sol.cm_dist, body_aero.panels, refined_panel_indices) + expected_cl = area_mean(sol.cl_dist, panels, refined_panel_indices) + expected_cd = area_mean(sol.cd_dist, panels, refined_panel_indices) + expected_cm = area_mean(sol.cm_dist, panels, refined_panel_indices) - # Handle NaN values that can occur in INVISCID models + # Handle NaN values that can occur in INVISCID models if isnan(expected_cl) @test isnan(sol.cl_unrefined_dist[unrefined_idx]) else @@ -89,7 +90,7 @@ area_weighted_mean(values, panels, indices) = end end - @testset "Unrefined force consistency: coeff*width*chord" begin + @testset "Unrefined section loads equal the sums of their panel loads" begin # Verify that for each unrefined section: # coeff_unrefined * width_unrefined * chord_unrefined ≈ # sum(coeff_panel * width_panel * chord_panel) for all panels in section @@ -145,6 +146,8 @@ area_weighted_mean(values, panels, indices) = @test cl_u * c_u * w_u ≈ sum_cl_cw rtol=1e-12 @test cd_u * c_u * w_u ≈ sum_cd_cw rtol=1e-12 @test cm_u * c_u * w_u ≈ sum_cm_cw rtol=1e-12 + @test sol.moment_unrefined_dist[unrefined_idx] ≈ + sum(sol.moment_dist[i] for i in pidxs) rtol=1e-12 end finally rm(settings_file; force=true) @@ -185,13 +188,14 @@ area_weighted_mean(values, panels, indices) = @test length(sol.cm_unrefined_dist) == wing.n_unrefined_sections # Verify unrefined coefficients are computed correctly using mapping + panels = body_aero.panels for unrefined_idx in 1:wing.n_unrefined_sections refined_panel_indices = findall(x -> x == unrefined_idx, wing.refined_panel_mapping) if !isempty(refined_panel_indices) - expected_cl = area_weighted_mean(sol.cl_dist, body_aero.panels, refined_panel_indices) - expected_cd = area_weighted_mean(sol.cd_dist, body_aero.panels, refined_panel_indices) - expected_cm = area_weighted_mean(sol.cm_dist, body_aero.panels, refined_panel_indices) + expected_cl = area_mean(sol.cl_dist, panels, refined_panel_indices) + expected_cd = area_mean(sol.cd_dist, panels, refined_panel_indices) + expected_cm = area_mean(sol.cm_dist, panels, refined_panel_indices) # Handle NaN for all coefficients if isnan(expected_cl) From c62bd58cc31cfc4906d62d8d119f56c5aa0adfe2 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 30 Sep 2026 14:59:13 +0200 Subject: [PATCH 6/6] Warn once per solve after halving; drop unreachable NaN branches; unit-test the relaxation cap Review of #402: one warning naming the final relaxation_factor instead of one per halving, the gamma_loop! docstring signature matches the function, and the unrefined-dist testsets no longer branch on NaN now that FEASIBLE is asserted. Co-Authored-By: Claude Opus 5.5 --- ...he-curved-reference-wing-does-not-conve.md | 2 +- src/solver.jl | 19 ++++---- test/solver/test_solver.jl | 12 ++++++ test/solver/test_unrefined_dist.jl | 43 +++---------------- 4 files changed, 30 insertions(+), 46 deletions(-) diff --git a/changelog.d/395-the-curved-reference-wing-does-not-conve.md b/changelog.d/395-the-curved-reference-wing-does-not-conve.md index ef5d42a6..e7e1fa43 100644 --- a/changelog.d/395-the-curved-reference-wing-does-not-conve.md +++ b/changelog.d/395-the-curved-reference-wing-does-not-conve.md @@ -1,6 +1,6 @@ ### Fixed -- The `LOOP` solver caps each panel's relaxation at `1 / (1 + π c |z_airf ⋅ AIC_ii|)`, the diagonal step of its self-induced stiffness, so narrow panels (cosine-spaced tips) no longer make it diverge to NaN; the curved reference wing and finely refined wings now converge at the default `relaxation_factor`, with unchanged iteration counts on wings that already converged. A `LOOP` solve that still does not converge keeps halving `relaxation_factor`, warning each time, while it is above 1e-3, and each attempt stops at the first non-finite residual. +- The `LOOP` solver caps each panel's relaxation at `1 / (1 + π c |z_airf ⋅ AIC_ii|)`, the diagonal step of its self-induced stiffness, so narrow panels (cosine-spaced tips) no longer make it diverge to NaN; the curved reference wing and finely refined wings now converge at the default `relaxation_factor`, with unchanged iteration counts on wings that already converged. A `LOOP` solve that still does not converge keeps halving `relaxation_factor` while it is above 1e-3 and warns once with the factor it stopped at, and each attempt stops at the first non-finite residual. ### Changed diff --git a/src/solver.jl b/src/solver.jl index 8a3e8b29..c507ef10 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -529,8 +529,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; unrefined_count_dist[target_unrefined_idx] += 1 end - # Coefficients and alpha are area-weighted means, the chord is the - # section area over its width; width and the moments stay summed. for target_unrefined_idx in section_range if unrefined_count_dist[target_unrefined_idx] > 0 count = unrefined_count_dist[target_unrefined_idx] @@ -663,9 +661,9 @@ Converge the circulation distribution and leave it in `solver.lr.gamma_new`, wit turning it into forces. Fills the solver's panel arrays, builds the AIC matrices, starts from `gamma_distribution` (or an elliptical/zero guess when it is `nothing` or `solver.use_gamma_prev` is false) and iterates; a `LOOP` solver that fails to converge -retries with half the relaxation factor, warning each time, while the factor is above -1e-3. The circulation half of [`solve!`](@ref), paired with [`calc_forces!`](@ref). -Returns `nothing`. +retries with half the relaxation factor while the factor is above 1e-3, and warns once +with the factor it stopped at. The circulation half of [`solve!`](@ref), paired with +[`calc_forces!`](@ref). Returns `nothing`. """ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma_distribution=nothing; log=false) where {P, U, T} @@ -733,12 +731,16 @@ function solve_base!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics, gamma solver.lr.gamma_new .= gamma_initial # Run main iteration loop gamma_loop!(solver, body_aero, panels, relaxation_factor; log) + initial_relaxation_factor = relaxation_factor while solver.solver_type == LOOP && !solver.lr.converged && relaxation_factor > 1e-3 relaxation_factor /= 2 - @warn "LOOP did not converge, retrying with relaxation_factor = $relaxation_factor" solver.lr.gamma_new .= gamma_initial gamma_loop!(solver, body_aero, panels, relaxation_factor; log) end + if relaxation_factor < initial_relaxation_factor + @warn "LOOP halved relaxation_factor from $initial_relaxation_factor to " * + "$relaxation_factor, converged: $(solver.lr.converged)" + end nothing end @@ -914,9 +916,8 @@ function panel_relaxation_factors!(panel_relaxation, relaxation_factor, AIC, z_a end """ - gamma_loop!(solver::Solver, AIC_x::Matrix{Float64}, - AIC_y::Matrix{Float64}, AIC_z::Matrix{Float64}, - panels::AbstractVector{<:Panel}, relaxation_factor::Float64; log=true) + gamma_loop!(solver::Solver, body_aero::BodyAerodynamics, + panels::AbstractVector{<:Panel}, relaxation_factor; log=true) Main iteration loop for calculating circulation distribution. diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 25ebf514..506cd66d 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -175,6 +175,18 @@ end end end +@testset "LOOP caps only the panels stiffer than 1 / relaxation_factor" begin + AIC = zeros(2, 2, 3) + AIC[1, 1, 3] = -1 / π # stiffness 1 + π * 1.0 * (1/π) = 2 + AIC[2, 2, 3] = 99 / (2π) # stiffness 1 + π * 2.0 * (99/2π) = 100 + z_airf_dist = [0.0 0.0 1.0; 0.0 0.0 1.0] + chord_dist = [1.0, 2.0] + panel_relaxation = VortexStepMethod.panel_relaxation_factors!(zeros(2), 0.03, AIC, + z_airf_dist, chord_dist) + @test panel_relaxation[1] === 0.03 + @test panel_relaxation[2] ≈ 1 / 100 +end + calc_forces_allocs(solver, body_aero) = (calc_forces!(solver, body_aero); @allocated calc_forces!(solver, body_aero)) diff --git a/test/solver/test_unrefined_dist.jl b/test/solver/test_unrefined_dist.jl index 3e9d3031..d601cead 100644 --- a/test/solver/test_unrefined_dist.jl +++ b/test/solver/test_unrefined_dist.jl @@ -60,30 +60,14 @@ area_mean(values, panels, indices) = expected_cd = area_mean(sol.cd_dist, panels, refined_panel_indices) expected_cm = area_mean(sol.cm_dist, panels, refined_panel_indices) - # Handle NaN values that can occur in INVISCID models - if isnan(expected_cl) - @test isnan(sol.cl_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) - end - if isnan(expected_cd) - @test isnan(sol.cd_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) - end - if isnan(expected_cm) - @test isnan(sol.cm_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) - end + @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) + @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) + @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) end end # Test 4: Verify physical consistency (lift coefficients should be positive at positive AoA) - # Skip test if values are NaN - if !any(isnan.(sol.cl_unrefined_dist)) - @test all(sol.cl_unrefined_dist .> 0.0) - end + @test all(sol.cl_unrefined_dist .> 0.0) finally rm(settings_file; force=true) @@ -197,22 +181,9 @@ area_mean(values, panels, indices) = expected_cd = area_mean(sol.cd_dist, panels, refined_panel_indices) expected_cm = area_mean(sol.cm_dist, panels, refined_panel_indices) - # Handle NaN for all coefficients - if isnan(expected_cl) - @test isnan(sol.cl_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) - end - if isnan(expected_cd) - @test isnan(sol.cd_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) - end - if isnan(expected_cm) - @test isnan(sol.cm_unrefined_dist[unrefined_idx]) - else - @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) - end + @test isapprox(sol.cl_unrefined_dist[unrefined_idx], expected_cl, rtol=1e-10) + @test isapprox(sol.cd_unrefined_dist[unrefined_idx], expected_cd, rtol=1e-10) + @test isapprox(sol.cm_unrefined_dist[unrefined_idx], expected_cm, rtol=1e-10) end end