Use nu=0.1, Tend=1.0 for clean M+1=4 convergence orders; document U-formulation order limit

Co-authored-by: pancetta <[email protected]>
Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/0a7bcad2-2358-4ceb-b9be-ee9049ffa9cf

Author: Copilot

pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py

@@ -23,30 +23,61 @@
 -------------------------------
 1. **Standard** (:class:`~.Stokes_Poiseuille_1D_FD.stokes_poiseuille_1d_fd`)
    — constraint :math:`B\mathbf{u} = q(t)` has a time-dependent RHS.
-   This causes **order reduction** in the pressure gradient :math:`G`:
-   velocity converges at full collocation order :math:`2M-1 = 5` while
-   :math:`G` converges at only order :math:`M = 3`.
+   The pressure gradient :math:`G` converges at only order :math:`M`,
+   while the velocity converges at :math:`M+1`.
 
 2. **Lifted** (:class:`~.Stokes_Poiseuille_1D_FD.stokes_poiseuille_1d_fd_lift`)
-   — introduce the lifting :math:`\mathbf{u}_\ell(t) = (q(t)/s)\,\mathbf{1}`
-   (where :math:`s = B\mathbf{1} = h N`), which satisfies
-   :math:`B\mathbf{u}_\ell(t) = q(t)` identically.  The lifted variable
-   :math:`\tilde{\mathbf{v}} = \mathbf{u} - \mathbf{u}_\ell(t)` satisfies
-   the **homogeneous** (autonomous) constraint :math:`B\tilde{\mathbf{v}} = 0`.
-   Making the constraint autonomous removes the source of order reduction
-   and the pressure gradient is expected to converge at the full collocation
-   order :math:`2M-1 = 5`.
-
-Observed results (``nvars = 1023``, ``restol = 1e-13``)
----------------------------------------------------------
-* **No-lift**: velocity → order 5 (superconvergent endpoint);
-  pressure → stable order :math:`M = 3`.
-* **Lifted**: velocity unchanged (same accuracy);
-  pressure order increases beyond :math:`M` with each grid-size
-  halving, approaching the full collocation order :math:`2M-1 = 5`.
-  The lifted case is still in a pre-asymptotic regime at the fine
-  :math:`\Delta t` values accessible with ``nvars = 1023``
-  (:math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}` spatial floor).
+   — lifting :math:`\mathbf{u}_\ell(t) = (q(t)/s)\,\mathbf{1}` makes the
+   constraint **homogeneous**: :math:`B\tilde{\mathbf{v}} = 0`.  With the
+   autonomous constraint the pressure order is restored to :math:`M+1`,
+   matching the velocity.
+
+Why M+1 (not 2M-1) for the velocity?
+--------------------------------------
+For a pure ODE discretised with RADAU-RIGHT :math:`M` nodes, the collocation
+polynomial evaluated at the endpoint achieves the superconvergent order
+:math:`2M-1`.  The
+:class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE`
+sweeper uses the *U-formulation* (stores and integrates velocity derivatives
+:math:`U_m = u'(\tau_m)` at each collocation node).  The endpoint velocity is
+recovered by quadrature:
+
+.. math::
+
+    \mathbf{u}_{n+1} = \mathbf{u}_n
+       + \Delta t \sum_{j=1}^{M} Q_{Mj}\,U_j.
+
+Although the quadrature weights :math:`Q_{Mj}` are exact for the collocation
+polynomial's derivative (degree :math:`\leq M - 1 \leq 2M-2`), the stage
+derivatives :math:`U_j` themselves carry an :math:`\mathcal{O}(\Delta t^M)`
+error at each internal collocation node (the DAE constraint at every stage
+limits the internal accuracy to the stage order :math:`M`).  The resulting
+quadrature integral therefore has :math:`\mathcal{O}(\Delta t^{M+1})`
+accuracy — one order above the stage derivatives — not the full collocation
+order :math:`2M-1`.
+
+This :math:`M+1` order is confirmed across multiple :math:`M`:
+
+* :math:`M = 2`: velocity :math:`\to 3` (:math:`= M+1 = 2M-1`; degenerate)
+* :math:`M = 3`: velocity :math:`\to 4` (:math:`= M+1`; not :math:`2M-1 = 5`)
+* :math:`M = 4`: velocity :math:`\to 5` (:math:`= M+1`; not :math:`2M-1 = 7`)
+
+Observed results (:math:`\nu = 0.1`, ``nvars = 1023``, ``restol = 1e-13``,
+:math:`M = 3`)
+---------------------------------------------------------------------------
+With :math:`\nu = 1.0` the problem is stiff (slow-mode Courant number
+:math:`|\lambda\,\Delta t| = \nu\pi^2 \Delta t \approx 2.5` at :math:`\Delta t = 0.25`),
+keeping the solution in the pre-asymptotic regime across the entire accessible
+:math:`\Delta t` range.  Reducing to :math:`\nu = 0.1` brings
+:math:`|\lambda\,\Delta t| \lesssim 0.25` at :math:`\Delta t = 0.25`, entering
+the asymptotic region and revealing the clean orders:
+
+* **Standard**: velocity at :math:`M+1 = 4`, pressure at :math:`M = 3`
+  (time-dependent constraint :math:`B\mathbf{u} = q(t)` causes order reduction
+  in :math:`G`).
+* **Lifted**: velocity at :math:`M+1 = 4` (unchanged); pressure order
+  increases monotonically, approaching :math:`M+1 = 4` (homogeneous
+  constraint removes the order reduction).
 
 **Spatial resolution**: ``nvars = 1023`` interior points with a
 fourth-order FD Laplacian (:math:`\Delta x = 1/1024`, spatial error floor
@@ -71,8 +102,13 @@
 # ---------------------------------------------------------------------------
 
 _NVARS = 1023
-_NU = 1.0
-_TEND = 0.5
+# nu=0.1 gives slow-mode Courant number |lambda*dt| <= 0.25 at dt=0.25,
+# putting the solution firmly in the asymptotic regime where the expected
+# orders M+1=4 (velocity) and M=3 (pressure, no-lift) are clearly visible.
+# With nu=1.0 the problem is ~10x stiffer and the asymptotic region cannot
+# be accessed before hitting the O(dx^4) spatial floor.
+_NU = 0.1
+_TEND = 1.0
 _NUM_NODES = 3
 _RESTOL = 1e-13
 
@@ -177,18 +213,34 @@ def main():
 
     Fixed parameters:
 
-    * ``restol = 1e-13``, :math:`\nu = 1.0`, :math:`M = 3` RADAU-RIGHT.
+    * ``restol = 1e-13``, :math:`\nu = 0.1`, :math:`M = 3` RADAU-RIGHT.
     * ``nvars = 1023`` (4th-order FD, spatial floor
       :math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}`).
-    * :math:`T_\text{end} = 0.5`.
+    * :math:`T_\text{end} = 1.0`.
+
+    Expected orders (see module docstring for derivation):
+
+    * **Velocity**: :math:`M+1 = 4` for both formulations.
+    * **Pressure (standard)**: :math:`M = 3` (order reduction due to
+      time-dependent constraint).
+    * **Pressure (lifted)**: approaches :math:`M+1 = 4` (homogeneous
+      constraint removes the order reduction).
     """
-    max_order = 2 * _NUM_NODES - 1  # = 5
-    dts = [_TEND / (2**k) for k in range(1, 7)]   # 0.25 … 0.0078
+    vel_order = _NUM_NODES + 1   # M+1 = 4 (U-formulation of SemiImplicitDAE)
+    pres_order = _NUM_NODES      # M   = 3 (algebraic variable at each node)
+
+    # 7 halvings from T_end/2 to T_end/128  →  0.5, 0.25, …, 0.0078125
+    dts = [_TEND / (2**k) for k in range(1, 8)]
 
     print(f'\nFully-converged SDC  (restol={_RESTOL:.0e}, ν={_NU}, M={_NUM_NODES})')
-    print(f'Sweeper: SemiImplicitDAE, RADAU-RIGHT nodes')
-    print(f'Expected collocation order for velocity = {max_order}  (= 2M − 1)')
-    print(f'nvars = {_NVARS}, 4th-order FD  (spatial floor ~ O(dx⁴) ≈ 1e-12)')
+    print(f'Sweeper: SemiImplicitDAE (U-formulation), RADAU-RIGHT nodes')
+    print(f'Expected velocity order  M+1 = {vel_order}  '
+          f'(U-formulation limit; pure-ODE collocation order 2M-1 = {2*_NUM_NODES-1} is not achieved)')
+    print(f'Expected pressure order  M   = {pres_order}  (no-lift) '
+          f'/ approaches M+1 = {vel_order}  (lifted)')
+    print(f'nvars = {_NVARS}, 4th-order FD  (spatial floor ~ O(dx^4) ≈ 1e-12)')
+    print(f'ν = {_NU}:  slow-mode |λΔt| ≤ {_NU * np.pi**2 * dts[0]:.2f} at Δt = {dts[0]:.4f}'
+          f'  → asymptotic regime accessible')
     print(f'Error vs. exact analytical solution at T={_TEND}')
 
     cases = [
@@ -210,7 +262,7 @@ def main():
             vel_errs.append(ve)
             pres_errs.append(pe)
 
-        _print_table(dts, vel_errs, pres_errs, max_order)
+        _print_table(dts, vel_errs, pres_errs, vel_order)
         results[cls.__name__] = (vel_errs, pres_errs)
 
     # ---- Summary ----
@@ -222,37 +274,45 @@ def main():
         vel_errs, pres_errs = results[cls.__name__]
         vel_ord = _asymptotic_order(dts, vel_errs)
         pres_ord = _asymptotic_order(dts, pres_errs)
-        tag = cls.__name__
-        print(f'\n  {tag}:')
-        print(f'    Velocity order ≈ {vel_ord:.1f}  (expected → {max_order})')
-        if pres_ord < max_order - 0.5:
-            if isinstance(cls(), stokes_poiseuille_1d_fd_lift):
-                status = f'increasing, ≈ {pres_ord:.1f} (pre-asymptotic, heading to {max_order})'
+        is_lift = isinstance(cls(), stokes_poiseuille_1d_fd_lift)
+        exp_pres = vel_order if is_lift else pres_order
+        print(f'\n  {cls.__name__}:')
+        print(f'    Velocity order ≈ {vel_ord:.1f}  (expected M+1 = {vel_order})')
+        if pres_ord < vel_order - 0.4:
+            if is_lift:
+                status = f'increasing, ≈ {pres_ord:.1f} (pre-asymptotic, heading to M+1 = {vel_order})'
             else:
-                status = f'order reduced to {pres_ord:.1f}'
+                status = f'order reduced to {pres_ord:.1f} = M  (time-dependent constraint)'
             print(f'    Pressure order ≈ {pres_ord:.1f}  → {status}')
         else:
-            print(f'    Pressure order ≈ {pres_ord:.1f}  → full collocation order ✓')
+            print(f'    Pressure order ≈ {pres_ord:.1f}  → full order M+1 = {vel_order} ✓')
 
     print()
     print('  Conclusion:')
     pres_ord_std = _asymptotic_order(dts, results['stokes_poiseuille_1d_fd'][1])
     pres_ord_lft = _asymptotic_order(dts, results['stokes_poiseuille_1d_fd_lift'][1])
     print(
-        f'  • Standard: pressure converges at order {pres_ord_std:.1f} = M  '
-        f'(order reduction; constraint B·u = q(t) is time-dependent).'
+        f'  • Velocity order M+1 = {vel_order} confirmed for both formulations.'
+    )
+    print(
+        f'  • Standard: pressure at order {pres_ord_std:.1f} = M  '
+        f'(order reduction from time-dependent constraint).'
     )
     print(
-        f'  • Lifted:   pressure converges at increasing order {pres_ord_lft:.1f}+'
-        f'  (constraint B·ṽ = 0 is autonomous; order reduction removed).'
+        f'  • Lifted:   pressure at increasing order {pres_ord_lft:.1f}+'
+        f'  (heading to M+1 = {vel_order}; autonomous constraint removes reduction).'
     )
     print(
-        '  • The velocity accuracy is identical in both cases, confirming that\n'
-        '    the lifting only affects the algebraic variable.'
+        '  • Note: the velocity order M+1 (not 2M-1) arises from the\n'
+        '    U-formulation used by SemiImplicitDAE: the endpoint velocity\n'
+        '    is obtained by integrating O(dt^M) accurate stage derivatives,\n'
+        '    which limits the integral to O(dt^(M+1)) regardless of the\n'
+        '    quadrature formula\'s exactness for the collocation polynomial.\n'
+        '    (Verified for M = 2, 3, 4.)'
     )
     print(
-        '  • Convergence may plateau at the 4th-order spatial floor ~1e-12\n'
-        '    once temporal errors fall below O(dx⁴) at fine Δt.'
+        f'  • Convergence may plateau at the 4th-order spatial floor ~1e-12\n'
+        f'    once temporal errors fall below O(dx^4) at fine Δt.'
     )