Aitken acceleration not working / working well for >=3 subdomains

[ikalash]

@ikalash discovered in her studies that there is a lack of convergence for Aitken acceleration for >= 3 subdomains.

For 3 subdomains, sometimes things work (laser weld) but sometimes the method is slower than classical relaxation (cuboid) - see attached image. For >3 subdomains, there is no convergence, whereas there is some convergence for classical relaxation with certain relaxation parameters (also in the attached image). FWIW, this result does seem to be consistent with the literature, but it may be worth verifying that there is no bug in the code.

Attached are tarballs for the cuboids problem to reproduce the issue.

3sd.tar.gz
4sd.tar.gz
5sd.tar.gz

[ikalash]

@lxmota I looked more carefully at my runs, and in the cuboid >2 subdomain runs, I was using the Aug. 4 version of Norma from right before your fix for the multi-domain case. I was using this sha: 8155ccd (this is before your fix, right?). For laser weld, I used the Aug. 6 version, so the results there are correct.

I think before you dig into this more, I need to rerun this with the recent bug fixes, and the earlier multi-subdomain bug fix. Please stay tuned.

[lxmota]

@ikalash I reproduced the cuboid cases and there is a bug, and it is worse than slow convergence: on these decompositions aitken secant can report convergence after a single Schwarz sweep and return a wrong answer with failed = false.

What happens

Running your 3sd case on current main, one step:

fixed theta=0.1     sweeps=76
aitken secant       sweeps=1

Comparing the two displacement fields:

subdomain 1: max|u_secant - u_fixed| = 2.045e-03   max|u_fixed| = 2.045e-03   relative = 1.000
subdomain 2: max|u_secant - u_fixed| = 4.095e-03   max|u_fixed| = 4.095e-03   relative = 1.000
subdomain 3: max|u_secant - u_fixed| = 2.053e-03   max|u_fixed| = 6.156e-03   relative = 0.334

The error is as large as the solution itself in two of the three subdomains. The log shows the mechanism:

[SCHWARZ] Aitken secant θ[iter=0] = 1.0000e-01   (once per interface)
[SCHWARZ] Criterion [0] |ΔU| = 7.59e-01 : [INITIAL]
[SCHWARZ] Aitken secant θ[iter=1] = 1.0000e-01
[SCHWARZ] Aitken secant θ[iter=1] = -0.0000e+00
[SCHWARZ] Criterion [1] |ΔU| = 3.25e-13 : |ΔU|/|U| = 3.02e-15 : [CONVERGED]

Why

On the DN and overlap path the first sweep of a stop finds an empty relaxation slot and falls back to λ_prev = interp_disp. That makes the residual identically zero, stores the incoming datum as the previous iterate, and writes back θ·interp + (1-θ)·interp, which equals the incoming datum for any θ. On the next sweep the iterate has therefore not moved, so d = g^(n) - g^(n-1) is exactly the zero vector and the secant factor -d·δ/‖δ‖² is exactly zero. The -0.0000e+00 in the log is that exact zero, not a small number.

A zero factor freezes the interface iterate. In a chain of three or more subdomains every subdomain then re-solves against unchanged data and reproduces its previous solution, so the whole sweep is a fixed point, |ΔU| collapses to 1e-13, and the displacement-based criterion reads it as convergence. The criterion cannot distinguish "nothing changed because we converged" from "nothing changed because the update was multiplied by zero".

Scope, established by experiment rather than argument:

• It is not the interface predictor. With interface predictor: false the 3sd secant run still stops after one sweep.
• It is not new at three subdomains. The two-subdomain static DN test produces the same exact-zero factor; it escapes only because there the freeze does not propagate to the whole system, so |ΔU| stays at 2.67e+00 and the run continues. More subdomains make the freeze total.
• Robin-Robin and impedance are immune. That path seeds its stored iterate with zeros rather than with the incoming datum, so there is no degenerate pair. Measured on the cantilever Robin-Robin benchmark, no exact-zero factors occur and Aitken is worth about a factor of ten: 599 sweeps for fixed theta = 0.5 against 56 for recursive and 58 for secant over 10 steps.
• Aitken recursive is immune, because below N0 it clears the stored residual and so never forms a factor from the degenerate pair.

This explains the N0 = 2 result in #215

Setting N0 = 2 makes the gate iter >= N0 skip the degenerate pair, so the first computed factor uses two genuine iterates. Measured on your 3sd case before the fix:

secant N0=1    1 sweep    (false convergence, wrong answer)
secant N0=2  106 sweeps   (converges properly)

So N0 = 2 was working around this bug rather than tuning a method parameter. I would hold #215 until you have re-measured on the fixed code; the reason to prefer 2 may simply be gone.

It also corrects something I posted in #215. The Dirichlet-Neumann cantilever number I gave there, where the secant took 20 sweeps against 40 for the recursive form, was itself produced by this bug: the zero factor discarded the first interface iterate, which happened to help that particular problem. With the fix, that case is 41 against 40 and N0 = 1 and N0 = 2 give identical counts. My argument that N0 = 2 "gives back a large advantage" on that case does not survive.

Fixes

One commit on main, d682a515, covering two things.

The degenerate pair now falls back to the input theta instead of producing a zero factor, so the first genuine secant factor forms one sweep later from two real iterates. Nothing else about the method changes.

Separately, a relaxation factor at or below 1e-12 is now recorded on the controller as having frozen an interface iterate, and the convergence criterion refuses to convert such a sweep into convergence, logging why. With the first fix in place this should not fire, but it turns any future route to a frozen update into a visible non-convergence rather than a silent wrong answer. It is the same blind spot already noted for the paired-impedance neutral mode: the criterion measures the change in the relaxed iterate, not the interface residual, and those differ by a factor of theta.

What is not a bug

Aitken recursive really is slower than fixed relaxation on these cuboid cases: 110 sweeps against 76 at three subdomains, and the 256 cap against 118 at four. That part matches what you reported and I found nothing wrong in the machinery behind it. The relaxation state is keyed per interface, only the Dirichlet side of each interface relaxes, and the state is cleared wholesale each stop, so a subdomain with two interfaces gets two independent factors as intended.

One more thing worth knowing before you rerun: the 8155ccd0 you used is immediately before 21afe107, "schwarz: key the relaxation state per interface, not per partner", from six hours later the same day. Before that commit the state was keyed by partner subdomain, which does not identify an interface, so a subdomain partnering more than one interface shared a single blended iterate. Your three-or-more-subdomain runs predate exactly the fix that matters most for them, as well as the two from yesterday.

[lxmota]

Follow-up with the post-fix numbers, and the 5sd case.

The fix resolves the false convergence

Three subdomains, one step, before and after d682a515:

                 sweeps   subdomain 1   subdomain 2   subdomain 3
before  secant        1       100%          100%          33.4%
after   secant      104     2.7e-05       1.9e-05       2.0e-06
        fixed        76      (reference)

The remaining difference after the fix is the tolerance band, not a discrepancy: the driver sets absolute tolerance: 1.0e-06, and the two runs stop at slightly different points inside it. The 104 sweeps also match the 106 that N0 = 2 needed before the fix, which is the direct confirmation that N0 = 2 was skipping the degenerate pair and that the fix now does that automatically at N0 = 1.

Worth saying plainly: on this problem neither Aitken form beats fixed relaxation, 104 against 76. The earlier appearance that the secant was dramatically better was the bug returning a wrong answer quickly.

Five subdomains

5sd.tar.gz ships no 5sd meshes. It contains the four 4sd meshes plus the 5sd input files and journal, so the case cannot be run as attached. The journal is complete, so I regenerated the meshes with Cubit 17.08: 800, 6400, 800, 6400 and 800 elements, alternating coarse and fine per its fmod(_i,2) branch. These are my meshes rather than yours, so treat the numbers as indicative until you reproduce them.

One step, on current main:

fixed theta=0.1          176 sweeps, converged
aitken recursive N0=1    256 sweeps, hit the cap
aitken recursive N0=2    256 sweeps, hit the cap
aitken secant N0=1       256 sweeps, hit the cap

This reproduces your original report at five subdomains, and it is not explained by the bug. Fixed relaxation converges; every Aitken variant fails to. The difference after the fix is that the secant now fails visibly instead of stopping after one sweep with a wrong answer.

Summarising across the chain lengths, sweeps for one step, the three and four subdomain rows from your meshes:

             fixed   recursive N0=1   recursive N0=2   secant N0=1
3 subdomains    76         110              100            104
4 subdomains   118         256 (cap)        220            n/a
5 subdomains   176         256 (cap)        256 (cap)      256 (cap)

So Aitken degrades monotonically with chain length here and is already losing to fixed relaxation at three subdomains. I found no defect behind that: the relaxation state is keyed per interface, only the Dirichlet side of each interface relaxes, and the state is cleared wholesale each stop, so a subdomain with two interfaces gets two independent factors as intended. It looks like the method property you expected from the literature rather than a coding error, now that the one real bug is out of the way.

Separate issue worth opening

When a step exhausts maximum iterations without converging, the run says nothing. The tail of the 5sd secant log is:

[SCHWARZ] Criterion [256] |ΔU| = 5.32e-02 : |ΔU|/|U| = 3.78e-04 : [CONVERGING]
[SCHWARZ] Performed 256 Schwarz Iterations
[DONE]    Simulation Complete

The relative error is 3.78e-04 against a relative tolerance of 1.0e-10, and the simulation reports success with failed = false. A non-converged Schwarz step should at least warn, and arguably should be able to abort on request, the way the stalled interface jump already does for the impedance condition. That is independent of everything above and I would rather track it on its own.

[ikalash]

Sorry about the issue with the 5sd.tar.gz. Most of the 5sd files are there, just not the .g files. The .jou is there though so these can easily be generated by running CUBIT.

[lxmota]

Yes, I used the Cubit journal files to generate the meshes. This is just FYI.

[ikalash]

@lxmota yes I figured. Sorry about the error, I was doing this quickly at the end of the day.

[lxmota]

@ikalash We are reviewing the open issues. Is this one still relevant? If it has been resolved or superseded, please close it; otherwise a brief status update would help us prioritize it.

[ikalash]

Yes, I think this is still an issue. Below are updated results for cuboids. You can see that Aitken improved (right plot) but is still pretty bad in terms of iterations and not that robust either for >2 subdomains.