Skip to content

Commit 0c3140e

Browse files
lmoresilmoresi
committed
Phase G: extend to ETD-1 / ETD-2 stage 1 — both blow up
Three stage-1 variants (BDF-1, ETD-1, ETD-2) on the two-Stokes architecture (stage 2 viscous Stokes with frozen η, body force = -∇·σ_VE, J2 radial return on σ_total = σ_VE + 2η·ε̇(v_pl)): | stage 1 | result | runaway step | peak σ_eq | |---------|-----------------|--------------|-----------| | BDF-1 | clean 120 steps | — | 1.68 | | ETD-1 | runaway | step 39 | 131 | | ETD-2 | runaway | step 27 | 102 | Diagnosis from per-step σ_VE vs σ_total ratio at step 5: * BDF-1: 1.09 (stage 2 adds ~9% — small) * ETD-1: 1.29 (stage 2 adds ~30% — large) The η_frozen mismatch: for BDF-1 the effective coupling to ε̇ is η_VE = 0.0476, which is what we use in stage 2 → minimal extra contribution. For ETD, σ_VE already absorbs the strain-rate term via 2·η_raw(1-α)·ε̇ ≈ 0.098·ε̇; stage 2 adding another 2·η_frozen·ε̇ overcounts the response, breaks stress balance, and ETD's analytical α faithfully propagates the imbalance forward into runaway. Lesson: η_frozen has to match the stage-1 integrator's effective coupling to ε̇. For ETD that means η_frozen = 0 (or vanishingly small) — but then stage 2 collapses to "no correction" and we're back to Phase F single-shot. The architecture really wants outer Picard iteration with the body force updated to use the *corrected* (admissible) stress, not σ_VE, so that the corrected state satisfies momentum balance. Single-shot breaks for ETD because radial return changes σ in yielded zones, which momentum-imbalances the result; ETD's analytical history term then carries that imbalance forward exponentially. BDF-1 dissipates the imbalance through its rational α (lesson #13 from PR #161). Per-step trace files captured the divergence cleanly. Next: try outer Picard with σ-admissible-driven body force (advice §11 from the web doc), and/or set η_frozen to match the integrator (η_raw·(1-α) for ETD-1; ETD-2 has no clean equivalent). Underworld development team with AI support from Claude Code (https://claude.com/claude-code)
1 parent 8e42fe8 commit 0c3140e

4 files changed

Lines changed: 137 additions & 16 deletions

File tree

‎docs/developer/design/_phase_g_bdf1.trace.txt‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
# Phase G two-Stokes operator-split: bdf1
2+
# η_eff_frozen=4.761905e-02
3+
# columns: step, t, V_top, snes_VE, snes_pl, sigma_eq_max_VE, sigma_eq_max_total, sigma_eq_max_admissible, u_VE_y_max, u_pl_y_max, yielded_fraction
4+
1 0.0500 +0.4985 1 1 4.128507e-02 8.253547e-02 8.253547e-02 1.648185e-05 4.234912e-05 0.050505
5+
2 0.1000 +0.4938 1 1 1.316340e-01 1.791930e-01 1.791930e-01 3.644807e-03 9.122831e-04 0.074380
6+
3 0.1500 +0.4862 1 1 2.419653e-01 3.118790e-01 3.118790e-01 1.755278e-02 5.023160e-03 0.081726
7+
4 0.2000 +0.4755 1 1 3.602410e-01 4.626600e-01 4.626600e-01 2.919148e-02 1.247657e-02 0.089073
8+
5 0.2500 +0.4619 1 1 4.876179e-01 6.270607e-01 6.270607e-01 3.621142e-02 2.122805e-02 0.089991
9+
6 0.3000 +0.4455 1 1 6.693912e-01 8.405240e-01 8.405240e-01 3.954460e-02 2.975425e-02 0.097337
10+
7 0.3500 +0.4263 1 1 9.004092e-01 1.200959e+00 1.200959e+00 4.157879e-02 3.639124e-02 0.099174
11+
8 0.4000 +0.4045 1 1 1.168381e+00 1.641332e+00 1.641332e+00 4.446546e-02 5.340553e-02 0.100092
12+
9 0.4500 +0.3802 1 1 1.440688e+00 2.099474e+00 2.099474e+00 4.779403e-02 8.345398e-02 0.100092
13+
10 0.5000 +0.3536 1 1 1.671884e+00 2.489029e+00 2.489029e+00 5.038318e-02 1.114998e-01 0.100092
14+
11 0.5500 +0.3247 1 1 2.155106e+00 3.461931e+00 3.461931e+00 8.313991e-02 1.415293e-01 0.100092
15+
12 0.6000 +0.2939 1 1 2.597636e+00 4.298415e+00 4.298415e+00 1.104659e-01 2.062055e-01 0.101928
16+
13 0.6500 +0.2612 1 1 2.797180e+00 4.657680e+00 4.657680e+00 1.564161e-01 2.422354e-01 0.101928
17+
14 0.7000 +0.2270 1 1 2.724759e+00 4.487114e+00 4.487114e+00 2.239714e-01 2.415248e-01 0.100092
18+
15 0.7500 +0.1913 1 1 2.670327e+00 4.262410e+00 4.262410e+00 2.538862e-01 2.115074e-01 0.096419
19+
16 0.8000 +0.1545 1 1 2.485657e+00 3.875198e+00 3.875198e+00 2.473791e-01 1.507465e-01 0.093664
20+
17 0.8500 +0.1167 1 1 2.141240e+00 3.033080e+00 3.033080e+00 1.940982e-01 9.404939e-02 0.089991
21+
18 0.9000 +0.0782 1 1 2.005078e+00 2.761229e+00 2.761229e+00 1.441770e-01 5.113701e-02 0.082645
22+
19 0.9500 +0.0392 1 1 1.679371e+00 2.144277e+00 2.144277e+00 6.941391e-02 3.884767e-02 0.078972
23+
20 1.0000 -0.0000 1 1 1.417554e+00 1.681398e+00 1.681398e+00 4.768886e-02 2.365449e-02 0.065197
24+
21 1.0500 -0.0392 1 1 1.252073e+00 1.435440e+00 1.435440e+00 3.231136e-02 1.154432e-02 0.054178
25+
22 1.1000 -0.0782 1 1 1.126075e+00 1.267879e+00 1.267879e+00 1.779313e-02 1.194893e-02 0.037649
26+
23 1.1500 -0.1167 1 1 1.011877e+00 1.115827e+00 1.115827e+00 1.432856e-02 1.011321e-02 0.006428
27+
24 1.2000 -0.1545 1 1 9.149813e-01 9.944294e-01 9.944294e-01 1.073262e-02 9.832712e-03 0.004591
28+
25 1.2500 -0.1913 1 1 8.318880e-01 8.981140e-01 8.981140e-01 8.755614e-03 9.963349e-03 0.012856
29+
26 1.3000 -0.2270 1 1 7.541963e-01 8.118726e-01 8.118726e-01 8.728405e-03 1.007611e-02 0.038567
30+
27 1.3500 -0.2612 1 1 6.779184e-01 7.304670e-01 7.304670e-01 8.731282e-03 1.014226e-02 0.056015
31+
28 1.4000 -0.2939 1 1 6.006306e-01 6.511986e-01 6.511986e-01 1.197136e-02 1.013980e-02 0.067034
32+
29 1.4500 -0.3247 1 1 5.227237e-01 5.750199e-01 5.750199e-01 1.681685e-02 1.086650e-02 0.076217
33+
30 1.5000 -0.3536 1 1 4.467655e-01 5.059827e-01 5.059827e-01 2.227344e-02 1.279448e-02 0.078972
34+
31 1.5500 -0.3802 1 1 3.889689e-01 5.514624e-01 5.514624e-01 2.729231e-02 1.542516e-02 0.081726
35+
32 1.6000 -0.4045 1 1 5.193100e-01 7.808144e-01 7.808144e-01 3.140208e-02 1.911398e-02 0.087236
36+
33 1.6500 -0.4263 1 1 7.519007e-01 1.169903e+00 1.169903e+00 3.457438e-02 3.482811e-02 0.089073
37+
34 1.7000 -0.4455 1 1 1.327625e+00 2.410562e+00 2.410562e+00 3.965486e-02 7.263343e-02 0.090909
38+
35 1.7500 -0.4619 1 1 3.052793e+00 5.596395e+00 5.596395e+00 6.263999e-02 1.672252e-01 0.093664
39+
36 1.8000 -0.4755 1 1 7.007915e+00 1.291836e+01 1.291836e+01 1.529835e-01 3.853982e-01 0.100092
40+
37 1.8500 -0.4862 1 1 1.618656e+01 2.992956e+01 2.992956e+01 3.857033e-01 1.021598e+00 0.106520
41+
38 1.9000 -0.4938 1 1 3.633115e+01 6.738515e+01 6.738515e+01 8.217616e-01 2.556586e+00 0.112948
42+
39 1.9500 -0.4985 1 1 7.149545e+01 1.329483e+02 1.329483e+02 3.332106e+00 5.614007e+00 0.110193

‎docs/developer/design/_phase_g_etd1.trace.txt‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
# Phase G two-Stokes operator-split: etd1
2+
# η_eff_frozen=4.761905e-02
3+
# columns: step, t, V_top, snes_VE, snes_pl, sigma_eq_max_VE, sigma_eq_max_total, sigma_eq_max_admissible, u_VE_y_max, u_pl_y_max, yielded_fraction
4+
1 0.0500 +0.4985 1 1 4.228343e-02 8.353688e-02 8.353688e-02 1.648185e-05 4.298041e-05 0.050505
5+
2 0.1000 +0.4938 1 1 1.338467e-01 1.816076e-01 1.816076e-01 3.708296e-03 9.478455e-04 0.074380
6+
3 0.1500 +0.4862 1 1 2.452311e-01 3.157243e-01 3.157243e-01 1.747473e-02 5.139399e-03 0.082645
7+
4 0.2000 +0.4755 1 1 3.648172e-01 4.679188e-01 4.679188e-01 2.889278e-02 1.267561e-02 0.089073
8+
5 0.2500 +0.4619 1 1 4.931296e-01 6.339584e-01 6.339584e-01 3.579862e-02 2.149520e-02 0.089991
9+
6 0.3000 +0.4455 1 1 6.767704e-01 8.497245e-01 8.497245e-01 3.897289e-02 3.006446e-02 0.097337
10+
7 0.3500 +0.4263 1 1 9.090430e-01 1.212227e+00 1.212227e+00 4.101584e-02 3.668803e-02 0.099174
11+
8 0.4000 +0.4045 1 1 1.177608e+00 1.653585e+00 1.653585e+00 4.391562e-02 5.371288e-02 0.100092
12+
9 0.4500 +0.3802 1 1 1.449137e+00 2.110172e+00 2.110172e+00 4.712469e-02 8.371426e-02 0.100092
13+
10 0.5000 +0.3536 1 1 1.678351e+00 2.495996e+00 2.495996e+00 4.919130e-02 1.114959e-01 0.100092
14+
11 0.5500 +0.3247 1 1 2.154955e+00 3.456698e+00 3.456698e+00 8.103724e-02 1.412088e-01 0.100092
15+
12 0.6000 +0.2939 1 1 2.587359e+00 4.273581e+00 4.273581e+00 1.076910e-01 2.047298e-01 0.101928
16+
13 0.6500 +0.2612 1 1 2.775767e+00 4.611472e+00 4.611472e+00 1.521676e-01 2.393180e-01 0.101928
17+
14 0.7000 +0.2270 1 1 2.703669e+00 4.442758e+00 4.442758e+00 2.165808e-01 2.377984e-01 0.100092
18+
15 0.7500 +0.1913 1 1 2.634097e+00 4.189665e+00 4.189665e+00 2.437180e-01 2.071882e-01 0.096419
19+
16 0.8000 +0.1545 1 1 2.450862e+00 3.809074e+00 3.809074e+00 2.355259e-01 1.470401e-01 0.093664
20+
17 0.8500 +0.1167 1 1 2.076273e+00 2.907532e+00 2.907532e+00 1.855852e-01 9.025236e-02 0.089991
21+
18 0.9000 +0.0782 1 1 1.938055e+00 2.634946e+00 2.634946e+00 1.349406e-01 4.948152e-02 0.082645
22+
19 0.9500 +0.0392 1 1 1.637303e+00 2.072370e+00 2.072370e+00 6.486592e-02 3.714019e-02 0.077135
23+
20 1.0000 -0.0000 1 1 1.393603e+00 1.648510e+00 1.648510e+00 4.490199e-02 2.217466e-02 0.065197
24+
21 1.0500 -0.0392 1 1 1.227265e+00 1.401392e+00 1.401392e+00 2.989512e-02 1.111236e-02 0.054178
25+
22 1.1000 -0.0782 1 1 1.102303e+00 1.235347e+00 1.235347e+00 1.687911e-02 1.131177e-02 0.037649
26+
23 1.1500 -0.1167 1 1 9.897544e-01 1.086570e+00 1.086570e+00 1.328152e-02 9.557927e-03 0.006428
27+
24 1.2000 -0.1545 1 1 8.944513e-01 9.684075e-01 9.684075e-01 9.791643e-03 9.723986e-03 0.004591
28+
25 1.2500 -0.1913 1 1 8.114818e-01 8.727156e-01 8.727156e-01 8.462881e-03 9.821042e-03 0.011938
29+
26 1.3000 -0.2270 1 1 7.333416e-01 7.861862e-01 7.861862e-01 8.404212e-03 9.907685e-03 0.040404
30+
27 1.3500 -0.2612 1 1 6.561896e-01 7.038870e-01 7.038870e-01 8.380816e-03 9.955806e-03 0.058770
31+
28 1.4000 -0.2939 1 1 5.776439e-01 6.231668e-01 6.231668e-01 1.182627e-02 9.941190e-03 0.067034
32+
29 1.4500 -0.3247 1 1 4.980084e-01 5.448549e-01 5.448549e-01 1.662201e-02 1.068939e-02 0.076217
33+
30 1.5000 -0.3536 1 1 4.198138e-01 4.728574e-01 4.728574e-01 2.200783e-02 1.253828e-02 0.078972
34+
31 1.5500 -0.3802 1 1 3.703329e-01 5.377762e-01 5.377762e-01 2.687670e-02 1.504018e-02 0.084481
35+
32 1.6000 -0.4045 1 1 5.071143e-01 7.587619e-01 7.587619e-01 3.089766e-02 1.851831e-02 0.087236
36+
33 1.6500 -0.4263 1 1 7.310406e-01 1.138017e+00 1.138017e+00 3.400894e-02 3.383885e-02 0.089073
37+
34 1.7000 -0.4455 1 1 1.289662e+00 2.341276e+00 2.341276e+00 3.869493e-02 7.049121e-02 0.090909
38+
35 1.7500 -0.4619 1 1 2.966282e+00 5.434488e+00 5.434488e+00 5.936084e-02 1.622793e-01 0.093664
39+
36 1.8000 -0.4755 1 1 6.785056e+00 1.250391e+01 1.250391e+01 1.443445e-01 3.728532e-01 0.100092
40+
37 1.8500 -0.4862 1 1 1.557334e+01 2.878723e+01 2.878723e+01 3.612655e-01 9.843966e-01 0.106520
41+
38 1.9000 -0.4938 1 1 3.538998e+01 6.561595e+01 6.561595e+01 7.871726e-01 2.473244e+00 0.112029
42+
39 1.9500 -0.4985 1 1 7.063478e+01 1.314397e+02 1.314397e+02 3.060833e+00 5.416528e+00 0.110193

‎docs/developer/design/_phase_g_etd2.trace.txt‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
# Phase G two-Stokes operator-split: etd2
2+
# η_eff_frozen=4.761905e-02
3+
# columns: step, t, V_top, snes_VE, snes_pl, sigma_eq_max_VE, sigma_eq_max_total, sigma_eq_max_admissible, u_VE_y_max, u_pl_y_max, yielded_fraction
4+
1 0.0500 +0.4985 1 1 2.131789e-02 6.250717e-02 6.250717e-02 1.648185e-05 2.972363e-05 0.033976
5+
2 0.1000 +0.4938 1 1 1.058113e-01 1.492579e-01 1.492579e-01 1.754050e-03 2.214438e-04 0.068871
6+
3 0.1500 +0.4862 1 1 2.145201e-01 2.771855e-01 2.771855e-01 2.618482e-02 3.615327e-03 0.081726
7+
4 0.2000 +0.4755 1 1 3.319387e-01 4.281065e-01 4.281065e-01 2.733942e-02 1.083505e-02 0.089073
8+
5 0.2500 +0.4619 1 1 4.595646e-01 5.889546e-01 5.889546e-01 4.185421e-02 1.933509e-02 0.089073
9+
6 0.3000 +0.4455 1 1 6.259885e-01 7.744456e-01 7.744456e-01 3.569678e-02 2.806569e-02 0.096419
10+
7 0.3500 +0.4263 1 1 8.413309e-01 1.102725e+00 1.102725e+00 4.749592e-02 3.462466e-02 0.098255
11+
8 0.4000 +0.4045 1 1 1.104505e+00 1.542842e+00 1.542842e+00 4.621091e-02 4.725708e-02 0.099174
12+
9 0.4500 +0.3802 1 1 1.347161e+00 1.952024e+00 1.952024e+00 5.391612e-02 7.356886e-02 0.100092
13+
10 0.5000 +0.3536 1 1 1.589223e+00 2.343676e+00 2.343676e+00 5.263648e-02 1.025499e-01 0.100092
14+
11 0.5500 +0.3247 1 1 1.965574e+00 3.087602e+00 3.087602e+00 1.090111e-01 1.181559e-01 0.100092
15+
12 0.6000 +0.2939 1 1 2.400881e+00 3.904230e+00 3.904230e+00 8.099861e-02 1.771298e-01 0.100092
16+
13 0.6500 +0.2612 1 1 2.579540e+00 4.216018e+00 4.216018e+00 1.778706e-01 2.116509e-01 0.100092
17+
14 0.7000 +0.2270 1 1 2.347043e+00 3.721574e+00 3.721574e+00 2.365827e-01 1.923570e-01 0.098255
18+
15 0.7500 +0.1913 1 1 2.455868e+00 3.835060e+00 3.835060e+00 2.097092e-01 1.831820e-01 0.098255
19+
16 0.8000 +0.1545 1 1 1.880064e+00 2.601273e+00 2.601273e+00 3.527839e-01 1.033590e-01 0.094582
20+
17 0.8500 +0.1167 1 1 1.788661e+00 2.547532e+00 2.547532e+00 2.528213e-01 8.337054e-02 0.092746
21+
18 0.9000 +0.0782 1 1 1.881565e+00 2.551438e+00 2.551438e+00 3.488893e-01 6.302538e-02 0.081726
22+
19 0.9500 +0.0392 1 1 1.952635e+00 2.986596e+00 2.986596e+00 3.248196e-01 7.449594e-02 0.078053
23+
20 1.0000 -0.0000 1 1 2.492530e+00 4.041323e+00 4.041323e+00 3.488298e-01 1.222615e-01 0.068871
24+
21 1.0500 -0.0392 1 1 2.797503e+00 4.610039e+00 4.610039e+00 3.750499e-01 1.405274e-01 0.060606
25+
22 1.1000 -0.0782 1 1 3.236328e+00 5.470478e+00 5.470478e+00 4.481368e-01 1.802310e-01 0.048669
26+
23 1.1500 -0.1167 1 1 4.582489e+00 8.020845e+00 8.020845e+00 5.779805e-01 3.302763e-01 0.027548
27+
24 1.2000 -0.1545 1 1 6.411467e+00 1.185821e+01 1.185821e+01 1.134602e+00 3.673761e-01 0.021120
28+
25 1.2500 -0.1913 1 1 1.213353e+01 2.213372e+01 2.213372e+01 1.734974e+00 9.365296e-01 0.041322
29+
26 1.3000 -0.2270 1 1 2.391635e+01 4.446023e+01 4.446023e+01 1.255055e+00 2.079103e+00 0.052342
30+
27 1.3500 -0.2612 1 1 5.392796e+01 1.016061e+02 1.016061e+02 3.386313e+00 4.504132e+00 0.078053

‎docs/developer/design/_phase_g_two_stokes.py‎

Lines changed: 23 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -100,19 +100,14 @@ def _build_setup(label):
100100
)
101101

102102

103-
def _build_stage1_VE(setup, label):
104-
"""Stage 1: VE Stokes solver (BDF-1, no in-residual yield)."""
103+
def _build_stage1_VE(setup, label, integrator="bdf", order=1):
104+
"""Stage 1: VE Stokes solver (configurable integrator)."""
105105
mesh = setup["mesh"]
106106
u_VE = setup["u_VE"]; p_VE = setup["p_VE"]
107107

108108
stokes_VE = uw.systems.Stokes(mesh, velocityField=u_VE, pressureField=p_VE)
109-
# NOTE: branch is off `development` (pre-PR-#161). The `integrator`
110-
# parameter ships in PR #161 — once that lands, swap to
111-
# `integrator='etd', order=1` (ETD-1) for the higher-accuracy
112-
# exponential factor. For now BDF-1 (the default) is the VE
113-
# predictor; matches ETD-1 numerically anyway in this regime.
114109
stokes_VE.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(
115-
stokes_VE.Unknowns, order=1,
110+
stokes_VE.Unknowns, integrator=integrator, order=order,
116111
)
117112
cm = stokes_VE.constitutive_model
118113
cm.Parameters.shear_viscosity_0 = ETA
@@ -198,7 +193,7 @@ def _eta_eff_bdf1(eta_raw, mu_raw, dt):
198193
return eta_raw * mu_raw * dt / (eta_raw + mu_raw * dt)
199194

200195

201-
def run_two_stokes(label, n_periods=1.5):
196+
def run_two_stokes(label, n_periods=1.5, integrator="bdf", order=1):
202197
setup = _build_setup(label)
203198
mesh = setup["mesh"]
204199

@@ -208,7 +203,8 @@ def run_two_stokes(label, n_periods=1.5):
208203
n_x_l = -np.sin(theta); n_y_l = np.cos(theta)
209204

210205
# Stage 1 builds the VE Stokes — psi_star lives on this solver's DDt
211-
stokes_VE, V_top_VE = _build_stage1_VE(setup, label)
206+
stokes_VE, V_top_VE = _build_stage1_VE(setup, label,
207+
integrator=integrator, order=order)
212208
DFDt = stokes_VE.Unknowns.DFDt
213209

214210
sigma_coords = DFDt.psi_star[0].coords
@@ -400,12 +396,23 @@ def run_two_stokes(label, n_periods=1.5):
400396

401397
def main():
402398
os.makedirs(OUT_DIR, exist_ok=True)
403-
cache = os.path.join(OUT_DIR, "phase_g_first_cut.npz")
404-
if os.path.exists(cache):
405-
print(f"=== two-Stokes cache hit: {cache} — skipping ===", flush=True)
406-
return
407-
print("=== Phase G two-Stokes first cut: BDF-1 VE + viscous plastic ===", flush=True)
408-
run_two_stokes("first_cut", n_periods=1.5)
399+
cases = [
400+
# (label, integrator, order)
401+
("bdf1", "bdf", 1),
402+
("etd1", "etd", 1),
403+
("etd2", "etd", 2), # headline experiment
404+
]
405+
for label, integrator, order in cases:
406+
cache = os.path.join(OUT_DIR, f"phase_g_{label}.npz")
407+
if os.path.exists(cache):
408+
print(f"\n=== {label}: cache hit, skipping ===", flush=True)
409+
continue
410+
print(
411+
f"\n=== Phase G two-Stokes ({integrator.upper()}-{order} VE "
412+
f"+ viscous plastic) ===",
413+
flush=True,
414+
)
415+
run_two_stokes(label, n_periods=1.5, integrator=integrator, order=order)
409416

410417

411418
if __name__ == "__main__":

0 commit comments

Comments
 (0)