From a92be2dc46cf846c32d34dbb08d26b254adffade Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 06:53:36 +0000 Subject: [PATCH 1/9] Redesign: use genericImplicitConstrained sweeper with MeshDAE state Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/96b37f2c-762b-4b92-a75a-9eca004dca16 --- .../Stokes_Poiseuille_1D_FD.py | 416 ++++++++++++++++++ .../Stokes_Poiseuille_1D_FD/__init__.py | 0 .../run_convergence.py | 212 +++++++++ 3 files changed, 628 insertions(+) create mode 100644 pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py create mode 100644 pySDC/playgrounds/Stokes_Poiseuille_1D_FD/__init__.py create mode 100644 pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py new file mode 100644 index 0000000000..428c2f87b4 --- /dev/null +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -0,0 +1,416 @@ +r""" +1-D unsteady Stokes / Poiseuille problem – semi-explicit index-1 DAE +===================================================================== + +Problem +------- +The 1-D unsteady Stokes equations on :math:`y \in [0, 1]` with a global +incompressibility constraint read + +.. math:: + + u_t = \nu\,u_{yy} + G(t) + f(y, t), \qquad u(0,t)=u(1,t)=0, + +.. math:: + + \int_0^1 u\,\mathrm{d}y = q(t), + +where :math:`G(t)` is the unknown (spatially uniform) pressure gradient and +:math:`q(t)` is a prescribed flow-rate. After finite-difference discretisation +this becomes the index-1 semi-explicit DAE + +.. math:: + + \mathbf{u}' = \nu A\,\mathbf{u} + G\,\mathbf{1} + \mathbf{f}(t), + +.. math:: + + 0 = B\,\mathbf{u} - q(t), + +where :math:`A` is the FD Laplacian, +:math:`B = h\,\mathbf{1}^T` (rectangle-rule integral, :math:`h = 1/(N+1)`), +:math:`\mathbf{1}` is the vector of ones, and :math:`G(t)` is the Lagrange +multiplier that enforces the flow-rate constraint. + +Manufactured solution +--------------------- +.. math:: + + u_\text{ex}(y, t) = \sin(\pi y)\,\sin(t), \quad + G_\text{ex}(t) = \cos(t). + +The corresponding manufactured forcing that makes the system consistent is + +.. math:: + + f(y, t) = \sin(\pi y)\cos(t) + \nu\pi^2\sin(\pi y)\sin(t) - \cos(t). + +IMEX split and saddle-point solve +---------------------------------- +The IMEX split used with the ``imex_1st_order`` sweeper is: + +* :math:`f_\text{impl} = \nu A\,\mathbf{u}` – autonomous diffusion (stiff). +* :math:`f_\text{expl} = G\,\mathbf{1} + \mathbf{f}(t)` – pressure gradient + and manufactured forcing (non-stiff). + +The pressure :math:`G` is the algebraic variable of the DAE. It is updated +at each SDC node during ``solve_system`` via a **Schur-complement saddle-point +solve**: + +1. Solve :math:`(I - \alpha\nu A)\,\mathbf{w} = \mathbf{r}` (standard diffusion + solve). +2. Solve :math:`(I - \alpha\nu A)\,\mathbf{v} = \mathbf{1}` (unit-load solve). +3. Enforce the constraint: + :math:`G = \bigl(q(t) - B\,\mathbf{w}\bigr) /\! \bigl(\alpha\,B\,\mathbf{v}\bigr)`. +4. Return :math:`\mathbf{u} = \mathbf{w} + \alpha\,G\,\mathbf{v}`. + +Between successive SDC nodes :math:`G` is stored in ``self._G_last`` and used +in the subsequent ``eval_f`` call. + +Spatial discretisation +----------------------- +A **fourth-order** FD Laplacian is used to push the spatial error floor to +:math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}` with ``nvars = 1023``. +Because the exact velocity profile satisfies homogeneous Dirichlet BCs no +boundary-correction vector :math:`b_\text{bc}` is needed. + +Classes +------- +stokes_poiseuille_1d_fd + Single-variant problem class. +""" + +import numpy as np +import scipy.sparse as sp +from scipy.sparse.linalg import spsolve + +from pySDC.core.problem import Problem, WorkCounter +from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh + + +# --------------------------------------------------------------------------- +# Fourth-order FD Laplacian (same stencil as AllenCahn_1D_FD) +# --------------------------------------------------------------------------- + +def _build_laplacian(n, dx): + r""" + Assemble the fourth-order FD Laplacian on *n* interior points with + **zero** Dirichlet boundary conditions. + + Interior rows use the centred stencil + :math:`(-u_{k-2}+16u_{k-1}-30u_k+16u_{k+1}-u_{k+2})/(12\Delta x^2)`; + the outermost rows use a 6-point one-sided stencil to maintain + fourth-order accuracy without ghost points. + + Parameters + ---------- + n : int + Number of interior grid points (must be ≥ 5). + dx : float + Grid spacing. + + Returns + ------- + scipy.sparse.csc_matrix + """ + inv12dx2 = 1.0 / (12.0 * dx**2) + L = sp.lil_matrix((n, n)) + + # Row 0: 6-point one-sided stencil (u_L = 0 → no boundary correction). + L[0, 0] = -15 + L[0, 1] = -4 + L[0, 2] = 14 + L[0, 3] = -6 + L[0, 4] = 1 + + # Row 1: standard centred stencil; u_{-1} = u_L = 0. + L[1, 0] = 16 + L[1, 1] = -30 + L[1, 2] = 16 + L[1, 3] = -1 + + # Interior rows 2 … n-3. + for k in range(2, n - 2): + L[k, k - 2] = -1 + L[k, k - 1] = 16 + L[k, k] = -30 + L[k, k + 1] = 16 + L[k, k + 2] = -1 + + # Row n-2: standard centred stencil; u_n = u_R = 0. + L[n - 2, n - 4] = -1 + L[n - 2, n - 3] = 16 + L[n - 2, n - 2] = -30 + L[n - 2, n - 1] = 16 + + # Row n-1: mirror of row 0 (u_R = 0 → no boundary correction). + L[n - 1, n - 5] = 1 + L[n - 1, n - 4] = -6 + L[n - 1, n - 3] = 14 + L[n - 1, n - 2] = -4 + L[n - 1, n - 1] = -15 + + return (L * inv12dx2).tocsc() + + +# --------------------------------------------------------------------------- +# Problem class +# --------------------------------------------------------------------------- + +class stokes_poiseuille_1d_fd(Problem): + r""" + 1-D unsteady Stokes / Poiseuille problem discretised with fourth-order + finite differences. + + **State variable**: velocity :math:`\mathbf{u} \in \mathbb{R}^N` + (N interior grid points, zero Dirichlet BCs). + + **Algebraic variable**: pressure gradient :math:`G(t) \in \mathbb{R}` + (Lagrange multiplier for the flow-rate constraint). + + **Exact solution** (manufactured): + + .. math:: + + u_\text{ex}(y_i, t) = \sin(\pi y_i)\,\sin(t), \quad + G_\text{ex}(t) = \cos(t). + + **Flow-rate constraint RHS**: + + .. math:: + + q(t) = B\,\mathbf{u}_\text{ex}(t) = C_B\,\sin(t), + \quad C_B = h\sum_i \sin(\pi y_i). + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity :math:`\nu` (default 1.0). + + Attributes + ---------- + dx : float + Grid spacing :math:`1/(N+1)`. + xvalues : numpy.ndarray + Interior grid-point :math:`y`-coordinates. + A : scipy.sparse.csc_matrix + :math:`\nu`-scaled fourth-order FD Laplacian. + ones : numpy.ndarray + Vector of ones (pressure-gradient coupling). + C_B : float + :math:`h\sum_i \sin(\pi y_i)` – coefficient in :math:`q(t)`. + s : float + :math:`h N` – discrete integral of the unit vector (:math:`B\mathbf{1}`). + _G_last : float + Most recently computed algebraic variable :math:`G`. + """ + + dtype_u = mesh + dtype_f = imex_mesh + + def __init__(self, nvars=127, nu=1.0): + if nvars < 5: + raise ValueError(f'nvars must be >= 5 for the 4th-order FD Laplacian; got {nvars}') + super().__init__((nvars, None, np.dtype('float64'))) + self._makeAttributeAndRegister('nvars', 'nu', localVars=locals(), readOnly=True) + + self.dx = 1.0 / (nvars + 1) + self.xvalues = np.linspace(self.dx, 1.0 - self.dx, nvars) + + # nu-scaled Laplacian (4th-order FD, zero BCs → no b_bc correction) + self.A = nu * _build_laplacian(nvars, self.dx) + + # Discrete-integral operator B = h * 1^T as a plain vector + self.ones = np.ones(nvars) + self.s = self.dx * nvars # B * 1 = h*N + self.C_B = self.dx * np.sum(np.sin(np.pi * self.xvalues)) # B * sin(pi*y) + + # Initialise algebraic variable to G_ex(0) = cos(0) = 1 + self._G_last = 1.0 + + self.work_counters['rhs'] = WorkCounter() + + # ----------------------------------------------------------------------- + # Manufactured-solution helpers + # ----------------------------------------------------------------------- + + def _q(self, t): + r""" + Flow-rate constraint RHS: :math:`q(t) = C_B\,\sin(t)`. + + Parameters + ---------- + t : float + + Returns + ------- + float + """ + return self.C_B * np.sin(t) + + def _forcing(self, t): + r""" + Manufactured forcing :math:`f(y_i, t)` that makes the system + consistent with :math:`u_\text{ex}` and :math:`G_\text{ex}`: + + .. math:: + + f(y, t) = \sin(\pi y)\cos(t) + + \nu\pi^2\sin(\pi y)\sin(t) - \cos(t). + + Parameters + ---------- + t : float + + Returns + ------- + numpy.ndarray, shape (nvars,) + """ + y = self.xvalues + return ( + np.sin(np.pi * y) * np.cos(t) + + self.nu * np.pi**2 * np.sin(np.pi * y) * np.sin(t) + - np.cos(t) * self.ones + ) + + # ----------------------------------------------------------------------- + # pySDC interface + # ----------------------------------------------------------------------- + + def eval_f(self, u, t): + r""" + Evaluate the IMEX right-hand side. + + * :math:`f_\text{impl} = \nu A\,\mathbf{u}` (stiff diffusion). + * :math:`f_\text{expl} = G\,\mathbf{1} + \mathbf{f}(t)` where + :math:`G` is taken from ``self._G_last`` (set by ``solve_system`` + immediately before this call, or initialised to :math:`G_\text{ex}(0)`). + + Parameters + ---------- + u : dtype_u + Current velocity. + t : float + Current time. + + Returns + ------- + f : imex_mesh + """ + f = self.dtype_f(self.init) + f.impl[:] = self.A.dot(np.asarray(u)) + f.expl[:] = self._G_last * self.ones + self._forcing(t) + self.work_counters['rhs']() + return f + + def solve_system(self, rhs, factor, u0, t): + r""" + Solve the saddle-point system via Schur-complement elimination. + + Given the SDC right-hand side :math:`\mathbf{r}` and step-size + prefactor :math:`\alpha = \Delta t\,\tilde{q}_{mm}`, find + :math:`(\mathbf{u}, G)` satisfying + + .. math:: + + (I - \alpha\nu A)\,\mathbf{u} = \mathbf{r} + \alpha\,G\,\mathbf{1}, + + .. math:: + + B\,\mathbf{u} = q(t). + + **Algorithm** (Schur complement): + + 1. :math:`\mathbf{w} = (I - \alpha\nu A)^{-1}\,\mathbf{r}`. + 2. :math:`\mathbf{v} = (I - \alpha\nu A)^{-1}\,\mathbf{1}`. + 3. :math:`G = \bigl(q(t) - B\,\mathbf{w}\bigr) / + \bigl(\alpha\,B\,\mathbf{v}\bigr)`. + 4. :math:`\mathbf{u} = \mathbf{w} + \alpha\,G\,\mathbf{v}`. + + The computed :math:`G` is stored in ``self._G_last`` for use in the + immediately following ``eval_f`` call. + + Parameters + ---------- + rhs : dtype_u + SDC right-hand side. + factor : float + Implicit step-size prefactor :math:`\alpha`. + u0 : dtype_u + Initial guess (unused; direct solver). + t : float + Current time. + + Returns + ------- + me : dtype_u + Updated velocity satisfying :math:`B\,\mathbf{u} = q(t)`. + """ + me = self.dtype_u(self.init) + rhs_np = np.asarray(rhs) + + if factor == 0.0: + # Trivial step: no diffusion, just copy rhs. + # Enforce constraint by computing the required G correction. + Br = self.dx * np.sum(rhs_np) + q_t = self._q(t) + # With factor=0 there is no v to scale, so we handle + # the constraint by projecting rhs onto the constraint manifold. + G_corr = (q_t - Br) / self.s + me[:] = rhs_np + G_corr * self.ones + self._G_last = G_corr + return me + + # Assemble and factor the diffusion operator. + M = sp.eye(self.nvars, format='csc') - factor * self.A + + # Step 1: diffusion-only solve. + w = spsolve(M, rhs_np) + + # Step 2: unit-load solve (effect of a unit G contribution). + v = spsolve(M, self.ones) + + # Step 3: enforce flow-rate constraint B*(w + factor*G*v) = q(t). + Bw = self.dx * np.sum(w) # B.dot(w) = h * sum(w) + Bv = self.dx * np.sum(v) # B.dot(v) = h * sum(v) > 0 + G_new = (self._q(t) - Bw) / (factor * Bv) + + # Step 4: assemble velocity. + me[:] = w + factor * G_new * v + + # Persist G for the next eval_f call. + self._G_last = float(G_new) + return me + + def u_exact(self, t): + r""" + Exact velocity :math:`u_\text{ex}(y_i, t) = \sin(\pi y_i)\,\sin(t)`. + + Parameters + ---------- + t : float + + Returns + ------- + me : dtype_u + """ + me = self.dtype_u(self.init, val=0.0) + me[:] = np.sin(np.pi * self.xvalues) * np.sin(t) + return me + + @staticmethod + def G_exact(t): + r""" + Exact pressure gradient :math:`G_\text{ex}(t) = \cos(t)`. + + Parameters + ---------- + t : float + + Returns + ------- + float + """ + return float(np.cos(t)) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/__init__.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py new file mode 100644 index 0000000000..0ab5d7cf54 --- /dev/null +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -0,0 +1,212 @@ +r""" +Temporal order of convergence for the 1-D Stokes/Poiseuille index-1 DAE +======================================================================== + +This script tests the temporal order of convergence of IMEX-SDC applied to +the 1-D unsteady Stokes / Poiseuille problem, formulated as a semi-explicit +index-1 DAE: + +.. math:: + + \mathbf{u}' = \nu A\,\mathbf{u} + G(t)\,\mathbf{1} + \mathbf{f}(t), + +.. math:: + + 0 = B\,\mathbf{u} - q(t), + +where :math:`G(t)` is the algebraic variable (pressure gradient) and the +constraint is enforced via a Schur-complement saddle-point solve inside +``solve_system``. + +Two quantities are tracked at :math:`T_\text{end}`: + +* **Velocity error** :math:`\|\mathbf{u}_h - \mathbf{u}_\text{ex}\|_\infty` + – the differential variable. +* **Pressure error** :math:`|G_h - G_\text{ex}(T_\text{end})|` + – the algebraic variable. + +With ``restol = 1e-13`` (fully converged SDC) and :math:`M = 3` +RADAU-RIGHT nodes the expected result is: + +* **Velocity**: converges at the full collocation order + :math:`2M - 1 = 5`. +* **Pressure gradient**: may show **order reduction** compared to the + velocity. For an index-1 DAE the algebraic variable is determined via a + Schur-complement involving a :math:`1/\Delta t` factor, which can degrade + the convergence order relative to the differential variable. + +**Spatial resolution**: ``nvars = 1023`` interior points with a +fourth-order FD Laplacian (:math:`\Delta x = 1/1024`, spatial error floor +:math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}`). + +Usage:: + + python run_convergence.py +""" + +import numpy as np + +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order +from pySDC.playgrounds.Stokes_Poiseuille_1D_FD.Stokes_Poiseuille_1D_FD import stokes_poiseuille_1d_fd + +# --------------------------------------------------------------------------- +# Fixed parameters +# --------------------------------------------------------------------------- + +# 4th-order FD Laplacian: spatial error ~ O(dx^4). +# With nvars=1023 (dx=1/1024), spatial floor ~ 1e-12. +_NVARS = 1023 +_NU = 1.0 +_TEND = 0.5 +_NUM_NODES = 3 +_RESTOL = 1e-13 + +_SWEEPER_PARAMS = { + 'quad_type': 'RADAU-RIGHT', + 'num_nodes': _NUM_NODES, + 'QI': 'LU', + 'QE': 'EE', + 'initial_guess': 'spread', +} + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _run(dt, restol=_RESTOL, max_iter=50): + """ + Run one simulation and return ``(uend_array, problem_instance)``. + + Parameters + ---------- + dt : float + Time-step size. + restol : float + Residual tolerance for SDC convergence. + max_iter : int + Maximum number of SDC iterations per step. + + Returns + ------- + u : numpy.ndarray + Numerical velocity at :math:`T_\text{end}`. + P : stokes_poiseuille_1d_fd + Problem instance (gives access to ``_G_last`` and exact solutions). + """ + desc = { + 'problem_class': stokes_poiseuille_1d_fd, + 'problem_params': {'nvars': _NVARS, 'nu': _NU}, + 'sweeper_class': imex_1st_order, + 'sweeper_params': _SWEEPER_PARAMS, + 'level_params': {'restol': restol, 'dt': dt}, + 'step_params': {'maxiter': max_iter}, + } + ctrl = controller_nonMPI(num_procs=1, controller_params={'logger_level': 40}, description=desc) + P = ctrl.MS[0].levels[0].prob + uend, _ = ctrl.run(u0=P.u_exact(0.0), t0=0.0, Tend=_TEND) + return np.asarray(uend).copy(), P + + +def _print_table(dts, vel_errs, pres_errs, vel_order, pres_order): + """Print a two-column convergence table.""" + header = ( + f' {"dt":>10} {"vel error":>14} {"vel ord":>8} {"exp":>4}' + f' {"pres error":>14} {"pres ord":>9} {"exp":>4}' + ) + print(header) + for i, dt in enumerate(dts): + ve = vel_errs[i] + pe = pres_errs[i] + if i > 0 and vel_errs[i - 1] > 0.0 and ve > 0.0: + vo = np.log(vel_errs[i - 1] / ve) / np.log(dts[i - 1] / dt) + vo_str = f'{vo:>8.2f}' + else: + vo_str = f'{"---":>8}' + if i > 0 and pres_errs[i - 1] > 0.0 and pe > 0.0: + po = np.log(pres_errs[i - 1] / pe) / np.log(dts[i - 1] / dt) + po_str = f'{po:>9.2f}' + else: + po_str = f'{"---":>9}' + print( + f' {dt:>10.5f} {ve:>14.4e} {vo_str} {vel_order:>4d}' + f' {pe:>14.4e} {po_str} {pres_order:>4d}' + ) + + +# --------------------------------------------------------------------------- +# Main convergence study +# --------------------------------------------------------------------------- + +def main(): + r""" + Temporal convergence study for the Stokes/Poiseuille index-1 DAE. + + Fixed parameters: + + * ``restol = 1e-13``, :math:`\nu = 1.0`, :math:`M = 3` RADAU-RIGHT. + * ``nvars = 1023`` (4th-order FD, spatial floor :math:`\approx 10^{-12}`). + * :math:`T_\text{end} = 0.5`. + + Expected collocation order :math:`2M-1 = 5` for the velocity. + Check whether the pressure gradient :math:`G` shows order reduction. + """ + max_order = 2 * _NUM_NODES - 1 # = 5 + + # dt range: temporal errors dominate for coarser dt; 6 halvings from T/2. + dts = [_TEND / (2**k) for k in range(1, 7)] # 0.25 … 0.0078 + + print(f'\nFully-converged IMEX-SDC (restol={_RESTOL:.0e}, ν={_NU}, M={_NUM_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'Error vs. exact analytical solution at T={_TEND}') + print(f'Pressure G is extracted from the Schur-complement solve\n') + + vel_errs = [] + pres_errs = [] + + for dt in dts: + u, P = _run(dt) + u_ex = np.asarray(P.u_exact(_TEND)).copy() + G_ex = P.G_exact(_TEND) + vel_errs.append(float(np.linalg.norm(u - u_ex, np.inf))) + pres_errs.append(abs(P._G_last - G_ex)) + + # Estimate asymptotic order for pressure from the middle refinement steps + # (avoid the regime where spatial floor dominates or pre-asymptotic). + orders_pres = [] + for i in range(2, len(dts)): + if pres_errs[i] > 0.0 and pres_errs[i - 1] > 0.0: + orders_pres.append( + np.log(pres_errs[i - 1] / pres_errs[i]) / np.log(dts[i - 1] / dts[i]) + ) + pres_est = int(round(np.median(orders_pres))) if orders_pres else 0 + + print('=' * 80) + _print_table(dts, vel_errs, pres_errs, max_order, pres_est) + print('=' * 80) + + print(f'\nSummary') + print(f' Velocity: converging to collocation order {max_order}') + print(f' Pressure gradient: observed order ≈ {pres_est}') + if pres_est < max_order: + print( + f' → Order reduction in the algebraic variable: {max_order} → {pres_est}.' + ) + print( + ' This is expected for index-1 DAEs: the Schur-complement solve' + ' introduces a 1/Δt factor that degrades the convergence of G by one order.' + ) + else: + print( + f' → No order reduction detected: both variables converge at order {max_order}.' + ) + print( + f' (Convergence may plateau at the 4th-order spatial floor ~1e-12' + f'\n once temporal errors fall below O(dx⁴).)' + ) + + +if __name__ == '__main__': + main() From 18f3a8864b3ba56729330e58ed7c5654ef49fb9c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 07:08:01 +0000 Subject: [PATCH 2/9] Implement Stokes/Poiseuille DAE playground with SemiImplicitDAE sweeper Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/96b37f2c-762b-4b92-a75a-9eca004dca16 --- .../Stokes_Poiseuille_1D_FD.py | 304 ++++++++++-------- .../run_convergence.py | 109 ++++--- 2 files changed, 228 insertions(+), 185 deletions(-) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py index 428c2f87b4..0e1769ee9c 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -15,9 +15,9 @@ \int_0^1 u\,\mathrm{d}y = q(t), -where :math:`G(t)` is the unknown (spatially uniform) pressure gradient and -:math:`q(t)` is a prescribed flow-rate. After finite-difference discretisation -this becomes the index-1 semi-explicit DAE +where :math:`G(t)` is the unknown pressure gradient and :math:`q(t)` is a +prescribed flow-rate. After finite-difference discretisation this becomes +the semi-explicit index-1 DAE .. math:: @@ -27,65 +27,83 @@ 0 = B\,\mathbf{u} - q(t), -where :math:`A` is the FD Laplacian, +where :math:`A` is the 4th-order FD Laplacian, :math:`B = h\,\mathbf{1}^T` (rectangle-rule integral, :math:`h = 1/(N+1)`), :math:`\mathbf{1}` is the vector of ones, and :math:`G(t)` is the Lagrange -multiplier that enforces the flow-rate constraint. +multiplier for the flow-rate constraint. Manufactured solution --------------------- .. math:: u_\text{ex}(y, t) = \sin(\pi y)\,\sin(t), \quad - G_\text{ex}(t) = \cos(t). + G_\text{ex}(t) = \cos(t). -The corresponding manufactured forcing that makes the system consistent is +The manufactured forcing consistent with this exact solution is .. math:: f(y, t) = \sin(\pi y)\cos(t) + \nu\pi^2\sin(\pi y)\sin(t) - \cos(t). -IMEX split and saddle-point solve ----------------------------------- -The IMEX split used with the ``imex_1st_order`` sweeper is: +State and sweeper +----------------- +The state variable uses :class:`~pySDC.projects.DAE.misc.meshDAE.MeshDAE` +initialised with ``nvars`` interior points: -* :math:`f_\text{impl} = \nu A\,\mathbf{u}` – autonomous diffusion (stiff). -* :math:`f_\text{expl} = G\,\mathbf{1} + \mathbf{f}(t)` – pressure gradient - and manufactured forcing (non-stiff). +* ``u.diff[:]`` – velocity on :math:`N` interior grid points. +* ``u.alg[0]`` – pressure gradient :math:`G`; ``u.alg[1:]`` is unused. -The pressure :math:`G` is the algebraic variable of the DAE. It is updated -at each SDC node during ``solve_system`` via a **Schur-complement saddle-point -solve**: +The compatible sweeper is +:class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE`, +which keeps the *U-formulation*: ``L.f[m].diff`` stores the velocity +derivative :math:`U_m = u'(t_m)`, and ``L.u[m].diff`` is reconstructed by +time-integration. Only the differential part is integrated; the algebraic +variable :math:`G` is enforced at every SDC node via ``solve_system``. -1. Solve :math:`(I - \alpha\nu A)\,\mathbf{w} = \mathbf{r}` (standard diffusion - solve). -2. Solve :math:`(I - \alpha\nu A)\,\mathbf{v} = \mathbf{1}` (unit-load solve). -3. Enforce the constraint: - :math:`G = \bigl(q(t) - B\,\mathbf{w}\bigr) /\! \bigl(\alpha\,B\,\mathbf{v}\bigr)`. -4. Return :math:`\mathbf{u} = \mathbf{w} + \alpha\,G\,\mathbf{v}`. +Saddle-point solve +------------------ +``solve_system`` bypasses Newton and solves the linear saddle-point system -Between successive SDC nodes :math:`G` is stored in ``self._G_last`` and used -in the subsequent ``eval_f`` call. +.. math:: + + (I - \alpha\nu A)\,\mathbf{U} - G\,\mathbf{1} + = \nu A\,\mathbf{u}_\text{approx} + \mathbf{f}(t), + +.. math:: + + B\bigl(\mathbf{u}_\text{approx} + \alpha\,\mathbf{U}\bigr) = q(t), + +directly by Schur-complement elimination: + +1. :math:`\mathbf{r}_\text{eff} = \nu A\,\mathbf{u}_\text{approx} + \mathbf{f}(t)`. +2. Solve :math:`(I - \alpha\nu A)\,\mathbf{w} = \mathbf{r}_\text{eff}`. +3. Solve :math:`(I - \alpha\nu A)\,\mathbf{v}_0 = \mathbf{1}`. +4. :math:`G = \bigl(q(t) - B\,\mathbf{u}_\text{approx} - \alpha\,B\,\mathbf{w}\bigr) + /\!\bigl(\alpha\,B\,\mathbf{v}_0\bigr)`. +5. Return :math:`\mathbf{U} = \mathbf{w} + G\,\mathbf{v}_0`, :math:`G`. + +At the SDC fixed point this reproduces the DAE collocation solution. Spatial discretisation ----------------------- -A **fourth-order** FD Laplacian is used to push the spatial error floor to +A **fourth-order** FD Laplacian pushes the spatial error floor to :math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}` with ``nvars = 1023``. -Because the exact velocity profile satisfies homogeneous Dirichlet BCs no -boundary-correction vector :math:`b_\text{bc}` is needed. +With homogeneous Dirichlet BCs the boundary-correction vector +:math:`b_\text{bc}` vanishes. Classes ------- stokes_poiseuille_1d_fd - Single-variant problem class. + Problem class (inherits from + :class:`~pySDC.projects.DAE.misc.problemDAE.ProblemDAE`). """ import numpy as np import scipy.sparse as sp from scipy.sparse.linalg import spsolve -from pySDC.core.problem import Problem, WorkCounter -from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh +from pySDC.projects.DAE.misc.problemDAE import ProblemDAE +from pySDC.projects.DAE.misc.meshDAE import MeshDAE # --------------------------------------------------------------------------- @@ -99,8 +117,9 @@ def _build_laplacian(n, dx): Interior rows use the centred stencil :math:`(-u_{k-2}+16u_{k-1}-30u_k+16u_{k+1}-u_{k+2})/(12\Delta x^2)`; - the outermost rows use a 6-point one-sided stencil to maintain - fourth-order accuracy without ghost points. + the outermost rows use a 6-point one-sided stencil. With homogeneous + BCs the boundary-correction vector :math:`b_\text{bc}` is identically + zero. Parameters ---------- @@ -116,14 +135,14 @@ def _build_laplacian(n, dx): inv12dx2 = 1.0 / (12.0 * dx**2) L = sp.lil_matrix((n, n)) - # Row 0: 6-point one-sided stencil (u_L = 0 → no boundary correction). + # Row 0: 6-point one-sided stencil (u_L = 0). L[0, 0] = -15 L[0, 1] = -4 L[0, 2] = 14 L[0, 3] = -6 L[0, 4] = 1 - # Row 1: standard centred stencil; u_{-1} = u_L = 0. + # Row 1: standard centred stencil; u_{-1} = 0. L[1, 0] = 16 L[1, 1] = -30 L[1, 2] = 16 @@ -137,13 +156,13 @@ def _build_laplacian(n, dx): L[k, k + 1] = 16 L[k, k + 2] = -1 - # Row n-2: standard centred stencil; u_n = u_R = 0. + # Row n-2: standard centred stencil; u_n = 0. L[n - 2, n - 4] = -1 L[n - 2, n - 3] = 16 L[n - 2, n - 2] = -30 L[n - 2, n - 1] = 16 - # Row n-1: mirror of row 0 (u_R = 0 → no boundary correction). + # Row n-1: mirror of row 0 (u_R = 0). L[n - 1, n - 5] = 1 L[n - 1, n - 4] = -6 L[n - 1, n - 3] = 14 @@ -157,16 +176,21 @@ def _build_laplacian(n, dx): # Problem class # --------------------------------------------------------------------------- -class stokes_poiseuille_1d_fd(Problem): +class stokes_poiseuille_1d_fd(ProblemDAE): r""" 1-D unsteady Stokes / Poiseuille problem discretised with fourth-order finite differences. - **State variable**: velocity :math:`\mathbf{u} \in \mathbb{R}^N` - (N interior grid points, zero Dirichlet BCs). + The state variable is a :class:`~pySDC.projects.DAE.misc.meshDAE.MeshDAE` + with shape ``(2, nvars)`` (two components of length ``nvars`` each): - **Algebraic variable**: pressure gradient :math:`G(t) \in \mathbb{R}` - (Lagrange multiplier for the flow-rate constraint). + * ``u.diff[:]`` – velocity on :math:`N` interior grid points. + * ``u.alg[0]`` – pressure gradient :math:`G`; ``u.alg[1:]`` is unused. + + This class is compatible with the + :class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE` + sweeper, which expects ``eval_f(u, du, t)`` to return the fully-implicit + DAE residual. **Exact solution** (manufactured): @@ -175,19 +199,15 @@ class stokes_poiseuille_1d_fd(Problem): u_\text{ex}(y_i, t) = \sin(\pi y_i)\,\sin(t), \quad G_\text{ex}(t) = \cos(t). - **Flow-rate constraint RHS**: - - .. math:: - - q(t) = B\,\mathbf{u}_\text{ex}(t) = C_B\,\sin(t), - \quad C_B = h\sum_i \sin(\pi y_i). - Parameters ---------- nvars : int Number of interior grid points (default 127; must be ≥ 5). nu : float Kinematic viscosity :math:`\nu` (default 1.0). + newton_tol : float + Tolerance passed to ``ProblemDAE``; unused since ``solve_system`` + is overridden with a direct solver (default 1e-10). Attributes ---------- @@ -198,22 +218,19 @@ class stokes_poiseuille_1d_fd(Problem): A : scipy.sparse.csc_matrix :math:`\nu`-scaled fourth-order FD Laplacian. ones : numpy.ndarray - Vector of ones (pressure-gradient coupling). + Vector of ones (pressure-gradient coupling, shape ``(nvars,)``). C_B : float - :math:`h\sum_i \sin(\pi y_i)` – coefficient in :math:`q(t)`. - s : float - :math:`h N` – discrete integral of the unit vector (:math:`B\mathbf{1}`). - _G_last : float - Most recently computed algebraic variable :math:`G`. + :math:`h\sum_i \sin(\pi y_i)` – coefficient in + :math:`q(t) = C_B\sin(t)`. """ - dtype_u = mesh - dtype_f = imex_mesh - - def __init__(self, nvars=127, nu=1.0): + def __init__(self, nvars=127, nu=1.0, newton_tol=1e-10): if nvars < 5: - raise ValueError(f'nvars must be >= 5 for the 4th-order FD Laplacian; got {nvars}') - super().__init__((nvars, None, np.dtype('float64'))) + raise ValueError( + f'nvars must be >= 5 for the 4th-order FD Laplacian; got {nvars}' + ) + # ProblemDAE: super().__init__((nvars, None, dtype)) → MeshDAE shape (2, nvars) + super().__init__(nvars=nvars, newton_tol=newton_tol) self._makeAttributeAndRegister('nvars', 'nu', localVars=locals(), readOnly=True) self.dx = 1.0 / (nvars + 1) @@ -222,15 +239,9 @@ def __init__(self, nvars=127, nu=1.0): # nu-scaled Laplacian (4th-order FD, zero BCs → no b_bc correction) self.A = nu * _build_laplacian(nvars, self.dx) - # Discrete-integral operator B = h * 1^T as a plain vector + # Discrete-integral operator B = h * 1^T (used as a plain vector) self.ones = np.ones(nvars) - self.s = self.dx * nvars # B * 1 = h*N - self.C_B = self.dx * np.sum(np.sin(np.pi * self.xvalues)) # B * sin(pi*y) - - # Initialise algebraic variable to G_ex(0) = cos(0) = 1 - self._G_last = 1.0 - - self.work_counters['rhs'] = WorkCounter() + self.C_B = self.dx * float(np.sum(np.sin(np.pi * self.xvalues))) # ----------------------------------------------------------------------- # Manufactured-solution helpers @@ -252,8 +263,8 @@ def _q(self, t): def _forcing(self, t): r""" - Manufactured forcing :math:`f(y_i, t)` that makes the system - consistent with :math:`u_\text{ex}` and :math:`G_\text{ex}`: + Manufactured forcing consistent with :math:`u_\text{ex}` and + :math:`G_\text{ex} = \cos(t)`: .. math:: @@ -276,117 +287,123 @@ def _forcing(self, t): ) # ----------------------------------------------------------------------- - # pySDC interface + # pySDC / SemiImplicitDAE interface # ----------------------------------------------------------------------- - def eval_f(self, u, t): + def eval_f(self, u, du, t): r""" - Evaluate the IMEX right-hand side. + Evaluate the fully-implicit DAE residual + :math:`F(\mathbf{u}, \mathbf{u}', t) = 0`: + + .. math:: + + F_\text{diff} = \mathbf{u}' - \nu A\,\mathbf{u} + - G\,\mathbf{1} - \mathbf{f}(t), - * :math:`f_\text{impl} = \nu A\,\mathbf{u}` (stiff diffusion). - * :math:`f_\text{expl} = G\,\mathbf{1} + \mathbf{f}(t)` where - :math:`G` is taken from ``self._G_last`` (set by ``solve_system`` - immediately before this call, or initialised to :math:`G_\text{ex}(0)`). + .. math:: + + F_\text{alg} = B\,\mathbf{u} - q(t). Parameters ---------- - u : dtype_u - Current velocity. + u : MeshDAE + Current state. ``u.diff[:]`` = velocity, ``u.alg[0]`` = G. + du : MeshDAE + Current derivative. ``du.diff[:]`` = :math:`\mathbf{u}'`. t : float Current time. Returns ------- - f : imex_mesh + f : MeshDAE + Residual: ``f.diff[:]`` = :math:`F_\text{diff}`, + ``f.alg[0]`` = :math:`F_\text{alg}` (scalar; zero at the + exact solution). """ - f = self.dtype_f(self.init) - f.impl[:] = self.A.dot(np.asarray(u)) - f.expl[:] = self._G_last * self.ones + self._forcing(t) + f = self.dtype_f(self.init, val=0.0) + u_vel = np.asarray(u.diff) + du_vel = np.asarray(du.diff) + G = float(u.alg[0]) + + f.diff[:] = du_vel - (self.A.dot(u_vel) + G * self.ones + self._forcing(t)) + f.alg[0] = self.dx * float(np.sum(u_vel)) - self._q(t) + self.work_counters['rhs']() return f - def solve_system(self, rhs, factor, u0, t): + def solve_system(self, impl_sys, u_approx, factor, u0, t): r""" - Solve the saddle-point system via Schur-complement elimination. + Direct Schur-complement solve for one SDC node (bypasses Newton). - Given the SDC right-hand side :math:`\mathbf{r}` and step-size - prefactor :math:`\alpha = \Delta t\,\tilde{q}_{mm}`, find - :math:`(\mathbf{u}, G)` satisfying + The ``SemiImplicitDAE`` sweeper calls this method to find + :math:`(\mathbf{U}_m, G_m)` satisfying .. math:: - (I - \alpha\nu A)\,\mathbf{u} = \mathbf{r} + \alpha\,G\,\mathbf{1}, + \mathbf{U}_m - \bigl(\nu A\,(\mathbf{u}_\text{approx} + + \alpha\,\mathbf{U}_m) + G_m\,\mathbf{1} + \mathbf{f}(t)\bigr) = 0, .. math:: - B\,\mathbf{u} = q(t). + B\,(\mathbf{u}_\text{approx} + \alpha\,\mathbf{U}_m) = q(t), - **Algorithm** (Schur complement): + which reduces to the system shown in the module docstring. - 1. :math:`\mathbf{w} = (I - \alpha\nu A)^{-1}\,\mathbf{r}`. - 2. :math:`\mathbf{v} = (I - \alpha\nu A)^{-1}\,\mathbf{1}`. - 3. :math:`G = \bigl(q(t) - B\,\mathbf{w}\bigr) / - \bigl(\alpha\,B\,\mathbf{v}\bigr)`. - 4. :math:`\mathbf{u} = \mathbf{w} + \alpha\,G\,\mathbf{v}`. + After the solve: - The computed :math:`G` is stored in ``self._G_last`` for use in the - immediately following ``eval_f`` call. + * ``me.diff[:]`` = :math:`\mathbf{U}_m` (velocity derivative at + the node, stored in ``L.f[m].diff`` by the sweeper). + * ``me.alg[0]`` = :math:`G_m` (stored in ``L.u[m].alg`` by the + sweeper). Parameters ---------- - rhs : dtype_u - SDC right-hand side. + impl_sys : callable + The :func:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE.F` + function (unused; linear system solved directly). + u_approx : MeshDAE + Current approximation of the velocity at the node. factor : float Implicit step-size prefactor :math:`\alpha`. - u0 : dtype_u + u0 : MeshDAE Initial guess (unused; direct solver). t : float Current time. Returns ------- - me : dtype_u - Updated velocity satisfying :math:`B\,\mathbf{u} = q(t)`. + me : MeshDAE + Contains :math:`(\mathbf{U}_m, G_m)`. """ - me = self.dtype_u(self.init) - rhs_np = np.asarray(rhs) - - if factor == 0.0: - # Trivial step: no diffusion, just copy rhs. - # Enforce constraint by computing the required G correction. - Br = self.dx * np.sum(rhs_np) - q_t = self._q(t) - # With factor=0 there is no v to scale, so we handle - # the constraint by projecting rhs onto the constraint manifold. - G_corr = (q_t - Br) / self.s - me[:] = rhs_np + G_corr * self.ones - self._G_last = G_corr - return me - - # Assemble and factor the diffusion operator. - M = sp.eye(self.nvars, format='csc') - factor * self.A + me = self.dtype_u(self.init, val=0.0) + u_vel_approx = np.asarray(u_approx.diff).copy() - # Step 1: diffusion-only solve. - w = spsolve(M, rhs_np) + # Effective RHS for the velocity block. + rhs_eff = self.A.dot(u_vel_approx) + self._forcing(t) + + # Assemble operator M = I - alpha * nu * A. + M = sp.eye(self.nvars, format='csc') - factor * self.A - # Step 2: unit-load solve (effect of a unit G contribution). - v = spsolve(M, self.ones) + # Step 1: solve without pressure. + w = spsolve(M, rhs_eff) - # Step 3: enforce flow-rate constraint B*(w + factor*G*v) = q(t). - Bw = self.dx * np.sum(w) # B.dot(w) = h * sum(w) - Bv = self.dx * np.sum(v) # B.dot(v) = h * sum(v) > 0 - G_new = (self._q(t) - Bw) / (factor * Bv) + # Step 2: unit-pressure response. + v0 = spsolve(M, self.ones) - # Step 4: assemble velocity. - me[:] = w + factor * G_new * v + # Step 3: enforce constraint B*(u_approx + alpha*(w + G*v0)) = q(t). + Bw = self.dx * float(np.sum(w)) + Bv0 = self.dx * float(np.sum(v0)) # > 0 for alpha > 0 + Bu_approx = self.dx * float(np.sum(u_vel_approx)) + G_new = (self._q(t) - Bu_approx - factor * Bw) / (factor * Bv0) - # Persist G for the next eval_f call. - self._G_last = float(G_new) + # Step 4: assemble velocity derivative and algebraic variable. + me.diff[:] = w + G_new * v0 + me.alg[0] = G_new return me def u_exact(self, t): r""" - Exact velocity :math:`u_\text{ex}(y_i, t) = \sin(\pi y_i)\,\sin(t)`. + Exact solution at time :math:`t`. Parameters ---------- @@ -394,16 +411,18 @@ def u_exact(self, t): Returns ------- - me : dtype_u + me : MeshDAE + ``me.diff[:]`` = :math:`\sin(\pi y_i)\sin(t)`, + ``me.alg[0]`` = :math:`\cos(t)`. """ me = self.dtype_u(self.init, val=0.0) - me[:] = np.sin(np.pi * self.xvalues) * np.sin(t) + me.diff[:] = np.sin(np.pi * self.xvalues) * np.sin(t) + me.alg[0] = np.cos(t) return me - @staticmethod - def G_exact(t): + def du_exact(self, t): r""" - Exact pressure gradient :math:`G_\text{ex}(t) = \cos(t)`. + Exact time derivative at time :math:`t`. Parameters ---------- @@ -411,6 +430,11 @@ def G_exact(t): Returns ------- - float + me : MeshDAE + ``me.diff[:]`` = :math:`\sin(\pi y_i)\cos(t)`, + ``me.alg[0]`` = :math:`-\sin(t)`. """ - return float(np.cos(t)) + me = self.dtype_u(self.init, val=0.0) + me.diff[:] = np.sin(np.pi * self.xvalues) * np.cos(t) + me.alg[0] = -np.sin(t) + return me diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py index 0ab5d7cf54..69fa228688 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -2,7 +2,7 @@ Temporal order of convergence for the 1-D Stokes/Poiseuille index-1 DAE ======================================================================== -This script tests the temporal order of convergence of IMEX-SDC applied to +This script tests the temporal order of convergence of SDC applied to the 1-D unsteady Stokes / Poiseuille problem, formulated as a semi-explicit index-1 DAE: @@ -14,9 +14,13 @@ 0 = B\,\mathbf{u} - q(t), -where :math:`G(t)` is the algebraic variable (pressure gradient) and the -constraint is enforced via a Schur-complement saddle-point solve inside -``solve_system``. +where :math:`G(t)` is the algebraic variable (pressure gradient) enforced via +a Schur-complement saddle-point solve inside ``solve_system``. + +The sweeper used is +:class:`~pySDC.playgrounds.DAE.genericImplicitDAE.genericImplicitConstrained`, +which integrates only the differential components and enforces the algebraic +constraint at every SDC node. Two quantities are tracked at :math:`T_\text{end}`: @@ -29,11 +33,19 @@ RADAU-RIGHT nodes the expected result is: * **Velocity**: converges at the full collocation order - :math:`2M - 1 = 5`. -* **Pressure gradient**: may show **order reduction** compared to the - velocity. For an index-1 DAE the algebraic variable is determined via a - Schur-complement involving a :math:`1/\Delta t` factor, which can degrade - the convergence order relative to the differential variable. + :math:`2M - 1 = 5`. Due to pre-asymptotic effects from the SDC + iteration, orders below 5 are seen at larger :math:`\Delta t`; the + error approaches the 4th-order spatial floor + :math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}` at fine + :math:`\Delta t`. +* **Pressure gradient**: shows **order reduction** to :math:`M = 3`. + The algebraic variable :math:`G` is the Lagrange multiplier determined + at each collocation node by the Schur-complement constraint solve. + At a single node :math:`t_m` the derivative :math:`U_m = u'(t_m)` has + accuracy :math:`\mathcal{O}(\Delta t^M)` (not the superconvergent order + :math:`2M-1`, which applies only to the endpoint value via Gauss + quadrature), and :math:`G_m` inherits that same order. Hence the + pressure at the endpoint converges at order :math:`M = 3`. **Spatial resolution**: ``nvars = 1023`` interior points with a fourth-order FD Laplacian (:math:`\Delta x = 1/1024`, spatial error floor @@ -47,7 +59,7 @@ import numpy as np from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order +from pySDC.projects.DAE.sweepers.semiImplicitDAE import SemiImplicitDAE from pySDC.playgrounds.Stokes_Poiseuille_1D_FD.Stokes_Poiseuille_1D_FD import stokes_poiseuille_1d_fd # --------------------------------------------------------------------------- @@ -66,7 +78,6 @@ 'quad_type': 'RADAU-RIGHT', 'num_nodes': _NUM_NODES, 'QI': 'LU', - 'QE': 'EE', 'initial_guess': 'spread', } @@ -77,7 +88,7 @@ def _run(dt, restol=_RESTOL, max_iter=50): """ - Run one simulation and return ``(uend_array, problem_instance)``. + Run one simulation and return ``(uend, problem_instance)``. Parameters ---------- @@ -90,15 +101,15 @@ def _run(dt, restol=_RESTOL, max_iter=50): Returns ------- - u : numpy.ndarray - Numerical velocity at :math:`T_\text{end}`. + uend : MeshDAE + Numerical solution at :math:`T_\\text{end}`. P : stokes_poiseuille_1d_fd - Problem instance (gives access to ``_G_last`` and exact solutions). + Problem instance (gives access to exact solutions). """ desc = { 'problem_class': stokes_poiseuille_1d_fd, 'problem_params': {'nvars': _NVARS, 'nu': _NU}, - 'sweeper_class': imex_1st_order, + 'sweeper_class': SemiImplicitDAE, 'sweeper_params': _SWEEPER_PARAMS, 'level_params': {'restol': restol, 'dt': dt}, 'step_params': {'maxiter': max_iter}, @@ -106,14 +117,14 @@ def _run(dt, restol=_RESTOL, max_iter=50): ctrl = controller_nonMPI(num_procs=1, controller_params={'logger_level': 40}, description=desc) P = ctrl.MS[0].levels[0].prob uend, _ = ctrl.run(u0=P.u_exact(0.0), t0=0.0, Tend=_TEND) - return np.asarray(uend).copy(), P + return uend, P -def _print_table(dts, vel_errs, pres_errs, vel_order, pres_order): +def _print_table(dts, vel_errs, pres_errs, vel_order, pres_order_str): """Print a two-column convergence table.""" header = ( f' {"dt":>10} {"vel error":>14} {"vel ord":>8} {"exp":>4}' - f' {"pres error":>14} {"pres ord":>9} {"exp":>4}' + f' {"pres error":>14} {"pres ord":>9}' ) print(header) for i, dt in enumerate(dts): @@ -131,7 +142,7 @@ def _print_table(dts, vel_errs, pres_errs, vel_order, pres_order): po_str = f'{"---":>9}' print( f' {dt:>10.5f} {ve:>14.4e} {vo_str} {vel_order:>4d}' - f' {pe:>14.4e} {po_str} {pres_order:>4d}' + f' {pe:>14.4e} {po_str}' ) @@ -146,65 +157,73 @@ def main(): Fixed parameters: * ``restol = 1e-13``, :math:`\nu = 1.0`, :math:`M = 3` RADAU-RIGHT. - * ``nvars = 1023`` (4th-order FD, spatial floor :math:`\approx 10^{-12}`). + * ``nvars = 1023`` (4th-order FD, spatial floor + :math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}`). * :math:`T_\text{end} = 0.5`. - Expected collocation order :math:`2M-1 = 5` for the velocity. - Check whether the pressure gradient :math:`G` shows order reduction. + Expected collocation order :math:`2M-1 = 5` for the velocity at the + endpoint. The pressure gradient :math:`G` (Lagrange multiplier) is + predicted to converge at only order :math:`M = 3`, because the Schur + complement computes :math:`G` from the derivative :math:`U_m`, which + achieves only :math:`\mathcal{O}(\Delta t^M)` accuracy at each internal + collocation node (no super-convergence for internal nodes). """ max_order = 2 * _NUM_NODES - 1 # = 5 - # dt range: temporal errors dominate for coarser dt; 6 halvings from T/2. + # dt range: 6 halvings from T/2 to T/64. dts = [_TEND / (2**k) for k in range(1, 7)] # 0.25 … 0.0078 - print(f'\nFully-converged IMEX-SDC (restol={_RESTOL:.0e}, ν={_NU}, M={_NUM_NODES})') + print(f'\nFully-converged SDC (restol={_RESTOL:.0e}, ν={_NU}, M={_NUM_NODES})') + print(f'Sweeper: SemiImplicitDAE (DAE-specific, constraint enforced at each node)') 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'Error vs. exact analytical solution at T={_TEND}') - print(f'Pressure G is extracted from the Schur-complement solve\n') + print(f'Error vs. exact analytical solution at T={_TEND}\n') vel_errs = [] pres_errs = [] for dt in dts: - u, P = _run(dt) - u_ex = np.asarray(P.u_exact(_TEND)).copy() - G_ex = P.G_exact(_TEND) - vel_errs.append(float(np.linalg.norm(u - u_ex, np.inf))) - pres_errs.append(abs(P._G_last - G_ex)) - - # Estimate asymptotic order for pressure from the middle refinement steps - # (avoid the regime where spatial floor dominates or pre-asymptotic). + uend, P = _run(dt) + u_ex = P.u_exact(_TEND) + vel_errs.append(float(np.linalg.norm(np.asarray(uend.diff) - np.asarray(u_ex.diff), np.inf))) + pres_errs.append(abs(float(uend.alg[0]) - float(u_ex.alg[0]))) + + # Estimate asymptotic order for pressure (skip the first point). orders_pres = [] for i in range(2, len(dts)): if pres_errs[i] > 0.0 and pres_errs[i - 1] > 0.0: orders_pres.append( np.log(pres_errs[i - 1] / pres_errs[i]) / np.log(dts[i - 1] / dts[i]) ) - pres_est = int(round(np.median(orders_pres))) if orders_pres else 0 + pres_est = round(float(np.median(orders_pres)), 1) if orders_pres else float('nan') + pres_order_str = f'≈ {pres_est:.1f}' print('=' * 80) - _print_table(dts, vel_errs, pres_errs, max_order, pres_est) + _print_table(dts, vel_errs, pres_errs, max_order, pres_order_str) print('=' * 80) print(f'\nSummary') print(f' Velocity: converging to collocation order {max_order}') - print(f' Pressure gradient: observed order ≈ {pres_est}') - if pres_est < max_order: + print(f' Pressure gradient: observed order {pres_order_str}') + if pres_est < max_order - 0.5: print( - f' → Order reduction in the algebraic variable: {max_order} → {pres_est}.' + f' → Order reduction in the algebraic variable: {max_order} → {pres_order_str}.' ) print( - ' This is expected for index-1 DAEs: the Schur-complement solve' - ' introduces a 1/Δt factor that degrades the convergence of G by one order.' + ' The pressure gradient G (Lagrange multiplier) converges at order M = 3,\n' + ' not at the superconvergent velocity order 2M-1 = 5.\n' + ' Reason: the Schur complement computes G_m from the velocity derivative\n' + ' U_m = u\'(t_m), which has only O(dt^M) accuracy at internal collocation\n' + ' nodes (super-convergence order 2M-1 applies only to the endpoint value\n' + ' via the Gauss-quadrature formula, not to the nodal derivatives).' ) else: print( - f' → No order reduction detected: both variables converge at order {max_order}.' + f' → No order reduction detected: both variables converge at order ≈ {max_order}.' ) print( - f' (Convergence may plateau at the 4th-order spatial floor ~1e-12' - f'\n once temporal errors fall below O(dx⁴).)' + f' (Convergence may plateau at the 4th-order spatial floor ~1e-12\n' + f' once temporal errors fall below O(dx⁴).)' ) From d6aa8013d2144322236f27563cb965870b938c6e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 07:32:23 +0000 Subject: [PATCH 3/9] Add constraint-lifting variant stokes_poiseuille_1d_fd_lift; update run_convergence.py to compare both Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/fe78add1-f3a2-429c-b4d5-e9257c7405a5 --- .../Stokes_Poiseuille_1D_FD.py | 507 ++++++++++++------ .../run_convergence.py | 257 +++++---- 2 files changed, 486 insertions(+), 278 deletions(-) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py index 0e1769ee9c..451f7c91b6 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -5,7 +5,7 @@ Problem ------- The 1-D unsteady Stokes equations on :math:`y \in [0, 1]` with a global -incompressibility constraint read +flow-rate constraint read .. math:: @@ -36,10 +36,9 @@ --------------------- .. math:: - u_\text{ex}(y, t) = \sin(\pi y)\,\sin(t), \quad - G_\text{ex}(t) = \cos(t). + u_\text{ex}(y, t) = \sin(\pi y)\,\sin(t), \quad G_\text{ex}(t) = \cos(t). -The manufactured forcing consistent with this exact solution is +Forcing: .. math:: @@ -48,54 +47,59 @@ State and sweeper ----------------- The state variable uses :class:`~pySDC.projects.DAE.misc.meshDAE.MeshDAE` -initialised with ``nvars`` interior points: +with ``nvars`` interior points: * ``u.diff[:]`` – velocity on :math:`N` interior grid points. -* ``u.alg[0]`` – pressure gradient :math:`G`; ``u.alg[1:]`` is unused. +* ``u.alg[0]`` – pressure gradient :math:`G` (Lagrange multiplier). The compatible sweeper is -:class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE`, -which keeps the *U-formulation*: ``L.f[m].diff`` stores the velocity -derivative :math:`U_m = u'(t_m)`, and ``L.u[m].diff`` is reconstructed by -time-integration. Only the differential part is integrated; the algebraic -variable :math:`G` is enforced at every SDC node via ``solve_system``. +:class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE`. -Saddle-point solve ------------------- -``solve_system`` bypasses Newton and solves the linear saddle-point system +Without constraint lifting (class :class:`stokes_poiseuille_1d_fd`), the +constraint :math:`B\mathbf{u} = q(t)` has a time-dependent right-hand side, +which causes order reduction in :math:`G` to order :math:`M` (= 3 for +3 RADAU-RIGHT nodes). + +Constraint lifting (class :class:`stokes_poiseuille_1d_fd_lift`) +----------------------------------------------------------------- +Introduce the lifting function .. math:: - (I - \alpha\nu A)\,\mathbf{U} - G\,\mathbf{1} - = \nu A\,\mathbf{u}_\text{approx} + \mathbf{f}(t), + \mathbf{u}_\ell(t) = \frac{q(t)}{s}\,\mathbf{1}, \qquad + s = B\mathbf{1} = h N, + +which satisfies :math:`B\mathbf{u}_\ell(t) = q(t)` exactly. Let +:math:`\tilde{\mathbf{v}} = \mathbf{u} - \mathbf{u}_\ell(t)`. The lifted +variable satisfies the **homogeneous** (autonomous) constraint .. math:: - B\bigl(\mathbf{u}_\text{approx} + \alpha\,\mathbf{U}\bigr) = q(t), + 0 = B\,\tilde{\mathbf{v}}, -directly by Schur-complement elimination: +and evolves according to -1. :math:`\mathbf{r}_\text{eff} = \nu A\,\mathbf{u}_\text{approx} + \mathbf{f}(t)`. -2. Solve :math:`(I - \alpha\nu A)\,\mathbf{w} = \mathbf{r}_\text{eff}`. -3. Solve :math:`(I - \alpha\nu A)\,\mathbf{v}_0 = \mathbf{1}`. -4. :math:`G = \bigl(q(t) - B\,\mathbf{u}_\text{approx} - \alpha\,B\,\mathbf{w}\bigr) - /\!\bigl(\alpha\,B\,\mathbf{v}_0\bigr)`. -5. Return :math:`\mathbf{U} = \mathbf{w} + G\,\mathbf{v}_0`, :math:`G`. +.. math:: -At the SDC fixed point this reproduces the DAE collocation solution. + \tilde{\mathbf{v}}' = \nu A\,\tilde{\mathbf{v}} + + \bigl[\nu A\,\mathbf{u}_\ell(t) + + \mathbf{f}(t) + - \dot{\mathbf{u}}_\ell(t)\bigr] + + G\,\mathbf{1}. -Spatial discretisation ------------------------ -A **fourth-order** FD Laplacian pushes the spatial error floor to -:math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}` with ``nvars = 1023``. -With homogeneous Dirichlet BCs the boundary-correction vector -:math:`b_\text{bc}` vanishes. +Because the constraint :math:`B\tilde{\mathbf{v}} = 0` no longer contains a +time-dependent right-hand side, the Lagrange multiplier :math:`G` is expected +to converge at the full collocation order :math:`2M-1`, matching the velocity. Classes ------- stokes_poiseuille_1d_fd - Problem class (inherits from - :class:`~pySDC.projects.DAE.misc.problemDAE.ProblemDAE`). + No lifting; constraint :math:`B\mathbf{u} = q(t)` (time-dependent). + Pressure converges at order :math:`M`. + +stokes_poiseuille_1d_fd_lift + Constraint lifting; homogeneous :math:`B\tilde{\mathbf{v}} = 0`. + Expected to restore full order :math:`2M-1` in the pressure. """ import numpy as np @@ -103,7 +107,6 @@ from scipy.sparse.linalg import spsolve from pySDC.projects.DAE.misc.problemDAE import ProblemDAE -from pySDC.projects.DAE.misc.meshDAE import MeshDAE # --------------------------------------------------------------------------- @@ -173,41 +176,24 @@ def _build_laplacian(n, dx): # --------------------------------------------------------------------------- -# Problem class +# Shared base class # --------------------------------------------------------------------------- -class stokes_poiseuille_1d_fd(ProblemDAE): +class _StokesBase(ProblemDAE): r""" - 1-D unsteady Stokes / Poiseuille problem discretised with fourth-order - finite differences. - - The state variable is a :class:`~pySDC.projects.DAE.misc.meshDAE.MeshDAE` - with shape ``(2, nvars)`` (two components of length ``nvars`` each): - - * ``u.diff[:]`` – velocity on :math:`N` interior grid points. - * ``u.alg[0]`` – pressure gradient :math:`G`; ``u.alg[1:]`` is unused. - - This class is compatible with the - :class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE` - sweeper, which expects ``eval_f(u, du, t)`` to return the fully-implicit - DAE residual. + Shared setup for the Stokes/Poiseuille problem classes. - **Exact solution** (manufactured): - - .. math:: - - u_\text{ex}(y_i, t) = \sin(\pi y_i)\,\sin(t), \quad - G_\text{ex}(t) = \cos(t). + Builds the 4th-order FD Laplacian, grid coordinates, and the + manufactured-forcing helpers. Parameters ---------- nvars : int - Number of interior grid points (default 127; must be ≥ 5). + Number of interior grid points (must be ≥ 5). nu : float - Kinematic viscosity :math:`\nu` (default 1.0). + Kinematic viscosity. newton_tol : float - Tolerance passed to ``ProblemDAE``; unused since ``solve_system`` - is overridden with a direct solver (default 1e-10). + Tolerance passed to :class:`~pySDC.projects.DAE.misc.problemDAE.ProblemDAE`. Attributes ---------- @@ -218,46 +204,36 @@ class stokes_poiseuille_1d_fd(ProblemDAE): A : scipy.sparse.csc_matrix :math:`\nu`-scaled fourth-order FD Laplacian. ones : numpy.ndarray - Vector of ones (pressure-gradient coupling, shape ``(nvars,)``). + Vector of ones, shape ``(nvars,)``. C_B : float - :math:`h\sum_i \sin(\pi y_i)` – coefficient in - :math:`q(t) = C_B\sin(t)`. + :math:`h\sum_i \sin(\pi y_i)` — coefficient in + :math:`q(t) = C_B \sin(t)`. + s : float + :math:`B \mathbf{1} = h N` — discrete integral of the all-ones vector. """ - def __init__(self, nvars=127, nu=1.0, newton_tol=1e-10): + def __init__(self, nvars, nu, newton_tol): if nvars < 5: raise ValueError( f'nvars must be >= 5 for the 4th-order FD Laplacian; got {nvars}' ) - # ProblemDAE: super().__init__((nvars, None, dtype)) → MeshDAE shape (2, nvars) super().__init__(nvars=nvars, newton_tol=newton_tol) self._makeAttributeAndRegister('nvars', 'nu', localVars=locals(), readOnly=True) self.dx = 1.0 / (nvars + 1) self.xvalues = np.linspace(self.dx, 1.0 - self.dx, nvars) - # nu-scaled Laplacian (4th-order FD, zero BCs → no b_bc correction) + # nu-scaled Laplacian (4th-order FD, zero BCs) self.A = nu * _build_laplacian(nvars, self.dx) # Discrete-integral operator B = h * 1^T (used as a plain vector) self.ones = np.ones(nvars) self.C_B = self.dx * float(np.sum(np.sin(np.pi * self.xvalues))) - - # ----------------------------------------------------------------------- - # Manufactured-solution helpers - # ----------------------------------------------------------------------- + self.s = self.dx * nvars # B * ones = h * N def _q(self, t): r""" Flow-rate constraint RHS: :math:`q(t) = C_B\,\sin(t)`. - - Parameters - ---------- - t : float - - Returns - ------- - float """ return self.C_B * np.sin(t) @@ -270,14 +246,6 @@ def _forcing(self, t): f(y, t) = \sin(\pi y)\cos(t) + \nu\pi^2\sin(\pi y)\sin(t) - \cos(t). - - Parameters - ---------- - t : float - - Returns - ------- - numpy.ndarray, shape (nvars,) """ y = self.xvalues return ( @@ -286,124 +254,244 @@ def _forcing(self, t): - np.cos(t) * self.ones ) - # ----------------------------------------------------------------------- - # pySDC / SemiImplicitDAE interface - # ----------------------------------------------------------------------- + def _B_dot(self, v): + r"""Rectangle-rule integral: :math:`B \mathbf{v} = h\sum_i v_i`.""" + return self.dx * float(np.sum(v)) - def eval_f(self, u, du, t): + def _schur_solve(self, rhs_eff, v_approx, factor, constraint_rhs): r""" - Evaluate the fully-implicit DAE residual - :math:`F(\mathbf{u}, \mathbf{u}', t) = 0`: + Schur-complement saddle-point solve. + + Finds :math:`(\mathbf{U}, G)` satisfying .. math:: - F_\text{diff} = \mathbf{u}' - \nu A\,\mathbf{u} - - G\,\mathbf{1} - \mathbf{f}(t), + (I - \alpha\nu A)\,\mathbf{U} - G\,\mathbf{1} = \mathbf{r}_\text{eff}, .. math:: - F_\text{alg} = B\,\mathbf{u} - q(t). + B(\mathbf{v}_\text{approx} + \alpha\,\mathbf{U}) = c, + + where :math:`c` is ``constraint_rhs`` (``q(t)`` for the standard + formulation, ``0`` for the lifted formulation). Parameters ---------- - u : MeshDAE - Current state. ``u.diff[:]`` = velocity, ``u.alg[0]`` = G. - du : MeshDAE - Current derivative. ``du.diff[:]`` = :math:`\mathbf{u}'`. - t : float - Current time. + rhs_eff : numpy.ndarray + Effective velocity RHS :math:`\mathbf{r}_\text{eff}`. + v_approx : numpy.ndarray + Current velocity approximation at the node. + factor : float + Implicit prefactor :math:`\alpha = \Delta t\,\tilde{q}_{mm}`. + constraint_rhs : float + Right-hand side of the constraint equation. Returns ------- - f : MeshDAE - Residual: ``f.diff[:]`` = :math:`F_\text{diff}`, - ``f.alg[0]`` = :math:`F_\text{alg}` (scalar; zero at the - exact solution). + U : numpy.ndarray + Velocity derivative at the node. + G_new : float + Pressure gradient (Lagrange multiplier). """ - f = self.dtype_f(self.init, val=0.0) - u_vel = np.asarray(u.diff) - du_vel = np.asarray(du.diff) - G = float(u.alg[0]) + M = sp.eye(self.nvars, format='csc') - factor * self.A + w = spsolve(M, rhs_eff) + v0 = spsolve(M, self.ones) - f.diff[:] = du_vel - (self.A.dot(u_vel) + G * self.ones + self._forcing(t)) - f.alg[0] = self.dx * float(np.sum(u_vel)) - self._q(t) + Bw = self._B_dot(w) + # Bv0 is positive because (I - factor*A) is an M-matrix for typical + # factor values, so its inverse has positive row-sums, giving B*v0 > 0. + Bv0 = self._B_dot(v0) + Bv = self._B_dot(v_approx) + G_new = (constraint_rhs - Bv - factor * Bw) / (factor * Bv0) - self.work_counters['rhs']() - return f + U = w + G_new * v0 + return U, float(G_new) - def solve_system(self, impl_sys, u_approx, factor, u0, t): - r""" - Direct Schur-complement solve for one SDC node (bypasses Newton). - The ``SemiImplicitDAE`` sweeper calls this method to find - :math:`(\mathbf{U}_m, G_m)` satisfying +# --------------------------------------------------------------------------- +# Case 1: No lifting – constraint B·u = q(t) +# --------------------------------------------------------------------------- + +class stokes_poiseuille_1d_fd(_StokesBase): + r""" + 1-D Stokes/Poiseuille DAE **without** constraint lifting. + + The constraint :math:`B\mathbf{u} = q(t)` has a time-dependent RHS, + which causes **order reduction** in the pressure gradient :math:`G` to + order :math:`M` (number of collocation nodes) instead of the full + collocation order :math:`2M-1`. + + **Exact solution** (manufactured): + + .. math:: + + u_\text{ex}(y_i, t) = \sin(\pi y_i)\,\sin(t), \quad + G_\text{ex}(t) = \cos(t). + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity (default 1.0). + newton_tol : float + Unused; passed to base class (default 1e-10). + """ + + def __init__(self, nvars=127, nu=1.0, newton_tol=1e-10): + super().__init__(nvars=nvars, nu=nu, newton_tol=newton_tol) + + def eval_f(self, u, du, t): + r""" + Fully-implicit DAE residual: .. math:: - \mathbf{U}_m - \bigl(\nu A\,(\mathbf{u}_\text{approx} - + \alpha\,\mathbf{U}_m) + G_m\,\mathbf{1} + \mathbf{f}(t)\bigr) = 0, + F_\text{diff} = \mathbf{u}' - \nu A\,\mathbf{u} + - G\,\mathbf{1} - \mathbf{f}(t), .. math:: - B\,(\mathbf{u}_\text{approx} + \alpha\,\mathbf{U}_m) = q(t), + F_\text{alg} = B\,\mathbf{u} - q(t). + """ + f = self.dtype_f(self.init, val=0.0) + u_vel = np.asarray(u.diff) + du_vel = np.asarray(du.diff) + G = float(u.alg[0]) - which reduces to the system shown in the module docstring. + f.diff[:] = du_vel - (self.A.dot(u_vel) + G * self.ones + self._forcing(t)) + f.alg[0] = self._B_dot(u_vel) - self._q(t) - After the solve: + self.work_counters['rhs']() + return f - * ``me.diff[:]`` = :math:`\mathbf{U}_m` (velocity derivative at - the node, stored in ``L.f[m].diff`` by the sweeper). - * ``me.alg[0]`` = :math:`G_m` (stored in ``L.u[m].alg`` by the - sweeper). + def solve_system(self, impl_sys, u_approx, factor, u0, t): + r""" + Schur-complement solve with constraint :math:`B\mathbf{u} = q(t)`. Parameters ---------- impl_sys : callable - The :func:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE.F` - function (unused; linear system solved directly). + Unused; system solved directly. u_approx : MeshDAE - Current approximation of the velocity at the node. + Current velocity approximation at the node. factor : float - Implicit step-size prefactor :math:`\alpha`. + Implicit prefactor :math:`\alpha`. u0 : MeshDAE - Initial guess (unused; direct solver). + Unused (direct solver). t : float Current time. Returns ------- me : MeshDAE - Contains :math:`(\mathbf{U}_m, G_m)`. + ``me.diff[:]`` = velocity derivative :math:`\mathbf{U}_m`, + ``me.alg[0]`` = pressure gradient :math:`G_m`. """ me = self.dtype_u(self.init, val=0.0) - u_vel_approx = np.asarray(u_approx.diff).copy() + v_approx = np.asarray(u_approx.diff).copy() - # Effective RHS for the velocity block. - rhs_eff = self.A.dot(u_vel_approx) + self._forcing(t) + rhs_eff = self.A.dot(v_approx) + self._forcing(t) + U, G_new = self._schur_solve(rhs_eff, v_approx, factor, self._q(t)) - # Assemble operator M = I - alpha * nu * A. - M = sp.eye(self.nvars, format='csc') - factor * self.A + me.diff[:] = U + me.alg[0] = G_new + return me - # Step 1: solve without pressure. - w = spsolve(M, rhs_eff) + def u_exact(self, t): + r""" + Exact solution: ``diff`` = :math:`\sin(\pi y)\sin(t)`, + ``alg[0]`` = :math:`\cos(t)`. + """ + me = self.dtype_u(self.init, val=0.0) + me.diff[:] = np.sin(np.pi * self.xvalues) * np.sin(t) + me.alg[0] = np.cos(t) + return me - # Step 2: unit-pressure response. - v0 = spsolve(M, self.ones) + def du_exact(self, t): + r""" + Exact time derivative: ``diff`` = :math:`\sin(\pi y)\cos(t)`, + ``alg[0]`` = :math:`-\sin(t)`. + """ + me = self.dtype_u(self.init, val=0.0) + me.diff[:] = np.sin(np.pi * self.xvalues) * np.cos(t) + me.alg[0] = -np.sin(t) + return me - # Step 3: enforce constraint B*(u_approx + alpha*(w + G*v0)) = q(t). - Bw = self.dx * float(np.sum(w)) - Bv0 = self.dx * float(np.sum(v0)) # > 0 for alpha > 0 - Bu_approx = self.dx * float(np.sum(u_vel_approx)) - G_new = (self._q(t) - Bu_approx - factor * Bw) / (factor * Bv0) - # Step 4: assemble velocity derivative and algebraic variable. - me.diff[:] = w + G_new * v0 - me.alg[0] = G_new - return me +# --------------------------------------------------------------------------- +# Case 2: Constraint lifting – homogeneous constraint B·ṽ = 0 +# --------------------------------------------------------------------------- - def u_exact(self, t): +class stokes_poiseuille_1d_fd_lift(_StokesBase): + r""" + 1-D Stokes/Poiseuille DAE **with** constraint lifting. + + **Lifting function** + + .. math:: + + \mathbf{u}_\ell(t) = \frac{q(t)}{s}\,\mathbf{1}, \qquad + s = B\mathbf{1} = h N, + + satisfies :math:`B\mathbf{u}_\ell = q(t)` exactly. The lifted + variable :math:`\tilde{\mathbf{v}} = \mathbf{u} - \mathbf{u}_\ell(t)` + satisfies the **homogeneous** constraint + + .. math:: + + 0 = B\,\tilde{\mathbf{v}}, + + and evolves according to + + .. math:: + + \tilde{\mathbf{v}}' = \nu A\,\tilde{\mathbf{v}} + G\,\mathbf{1} + + \bigl[\nu A\,\mathbf{u}_\ell(t) + + \mathbf{f}(t) + - \dot{\mathbf{u}}_\ell(t)\bigr]. + + Because :math:`B\tilde{\mathbf{v}} = 0` is autonomous (no time + dependence), the Lagrange multiplier :math:`G` is expected to converge + at the full collocation order :math:`2M-1`, matching the velocity. + + **State variable** ``u``: + + * ``u.diff[:]`` = :math:`\tilde{\mathbf{v}}` (lifted velocity). + * ``u.alg[0]`` = :math:`G`. + + **Exact lifted solution**: + + .. math:: + + \tilde{\mathbf{v}}_\text{ex}(y_i, t) = + \sin(\pi y_i)\sin(t) - \mathbf{u}_\ell(t), + \quad G_\text{ex}(t) = \cos(t). + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity (default 1.0). + newton_tol : float + Unused; passed to base class (default 1e-10). + """ + + def __init__(self, nvars=127, nu=1.0, newton_tol=1e-10): + super().__init__(nvars=nvars, nu=nu, newton_tol=newton_tol) + # Precompute A*ones (needed in eval_f and solve_system). + self._A_ones = np.asarray(self.A.dot(self.ones)).copy() + + # ------------------------------------------------------------------ + # Lifting helpers + # ------------------------------------------------------------------ + + def lift(self, t): r""" - Exact solution at time :math:`t`. + Lifting function :math:`\mathbf{u}_\ell(t) = (q(t)/s)\,\mathbf{1}`. + + Satisfies :math:`B\mathbf{u}_\ell(t) = q(t)` exactly. Parameters ---------- @@ -411,30 +499,121 @@ def u_exact(self, t): Returns ------- - me : MeshDAE - ``me.diff[:]`` = :math:`\sin(\pi y_i)\sin(t)`, - ``me.alg[0]`` = :math:`\cos(t)`. + numpy.ndarray, shape (nvars,) """ - me = self.dtype_u(self.init, val=0.0) - me.diff[:] = np.sin(np.pi * self.xvalues) * np.sin(t) - me.alg[0] = np.cos(t) - return me + return (self.C_B * np.sin(t) / self.s) * self.ones - def du_exact(self, t): + def _lift_dot(self, t): + r"""Time derivative :math:`\dot{\mathbf{u}}_\ell(t) = (C_B\cos(t)/s)\,\mathbf{1}`.""" + return (self.C_B * np.cos(t) / self.s) * self.ones + + def _A_lift(self, t): + r""" + :math:`\nu A\,\mathbf{u}_\ell(t) = (C_B\sin(t)/s)\,A\mathbf{1}`. + + Non-zero only at boundary rows (boundary-correction effect of the + spatially-constant lift profile). + """ + return (self.C_B * np.sin(t) / self.s) * self._A_ones + + def _net_forcing(self, t): + r""" + Net explicit forcing for the lifted system: + + .. math:: + + \mathbf{F}_\text{net}(t) = \nu A\,\mathbf{u}_\ell(t) + + \mathbf{f}(t) + - \dot{\mathbf{u}}_\ell(t). + """ + return self._A_lift(t) + self._forcing(t) - self._lift_dot(t) + + # ------------------------------------------------------------------ + # pySDC / SemiImplicitDAE interface + # ------------------------------------------------------------------ + + def eval_f(self, v, dv, t): r""" - Exact time derivative at time :math:`t`. + Fully-implicit DAE residual for the lifted variable + :math:`\tilde{\mathbf{v}}`: + + .. math:: + + F_\text{diff} = \tilde{\mathbf{v}}' - \nu A\,\tilde{\mathbf{v}} + - G\,\mathbf{1} - \mathbf{F}_\text{net}(t), + + .. math:: + + F_\text{alg} = B\,\tilde{\mathbf{v}} \quad + (\text{homogeneous, zero at exact solution}). + """ + f = self.dtype_f(self.init, val=0.0) + v_vel = np.asarray(v.diff) + dv_vel = np.asarray(dv.diff) + G = float(v.alg[0]) + + f.diff[:] = dv_vel - ( + self.A.dot(v_vel) + G * self.ones + self._net_forcing(t) + ) + f.alg[0] = self._B_dot(v_vel) # B·ṽ (should vanish at exact solution) + + self.work_counters['rhs']() + return f + + def solve_system(self, impl_sys, v_approx, factor, u0, t): + r""" + Schur-complement solve with **homogeneous** constraint + :math:`B(\tilde{\mathbf{v}}_\text{approx} + \alpha\,\tilde{\mathbf{U}}) = 0`. Parameters ---------- + impl_sys : callable + Unused; system solved directly. + v_approx : MeshDAE + Current lifted velocity approximation at the node. + factor : float + Implicit prefactor :math:`\alpha`. + u0 : MeshDAE + Unused (direct solver). t : float + Current time. Returns ------- me : MeshDAE - ``me.diff[:]`` = :math:`\sin(\pi y_i)\cos(t)`, - ``me.alg[0]`` = :math:`-\sin(t)`. + ``me.diff[:]`` = lifted velocity derivative :math:`\tilde{\mathbf{U}}_m`, + ``me.alg[0]`` = pressure gradient :math:`G_m`. """ me = self.dtype_u(self.init, val=0.0) - me.diff[:] = np.sin(np.pi * self.xvalues) * np.cos(t) + vv = np.asarray(v_approx.diff).copy() + + rhs_eff = self.A.dot(vv) + self._net_forcing(t) + U, G_new = self._schur_solve(rhs_eff, vv, factor, 0.0) + + me.diff[:] = U + me.alg[0] = G_new + return me + + def u_exact(self, t): + r""" + Exact lifted solution at time :math:`t`: + + ``diff`` = :math:`\sin(\pi y)\sin(t) - \mathbf{u}_\ell(t)`, + ``alg[0]`` = :math:`\cos(t)`. + """ + me = self.dtype_u(self.init, val=0.0) + me.diff[:] = np.sin(np.pi * self.xvalues) * np.sin(t) - self.lift(t) + me.alg[0] = np.cos(t) + return me + + def du_exact(self, t): + r""" + Exact time derivative of the lifted solution: + + ``diff`` = :math:`\sin(\pi y)\cos(t) - \dot{\mathbf{u}}_\ell(t)`, + ``alg[0]`` = :math:`-\sin(t)`. + """ + me = self.dtype_u(self.init, val=0.0) + me.diff[:] = np.sin(np.pi * self.xvalues) * np.cos(t) - self._lift_dot(t) me.alg[0] = -np.sin(t) return me diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py index 69fa228688..53031f660e 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -12,40 +12,41 @@ .. math:: - 0 = B\,\mathbf{u} - q(t), - -where :math:`G(t)` is the algebraic variable (pressure gradient) enforced via -a Schur-complement saddle-point solve inside ``solve_system``. - -The sweeper used is -:class:`~pySDC.playgrounds.DAE.genericImplicitDAE.genericImplicitConstrained`, -which integrates only the differential components and enforces the algebraic -constraint at every SDC node. - -Two quantities are tracked at :math:`T_\text{end}`: - -* **Velocity error** :math:`\|\mathbf{u}_h - \mathbf{u}_\text{ex}\|_\infty` - – the differential variable. -* **Pressure error** :math:`|G_h - G_\text{ex}(T_\text{end})|` - – the algebraic variable. - -With ``restol = 1e-13`` (fully converged SDC) and :math:`M = 3` -RADAU-RIGHT nodes the expected result is: - -* **Velocity**: converges at the full collocation order - :math:`2M - 1 = 5`. Due to pre-asymptotic effects from the SDC - iteration, orders below 5 are seen at larger :math:`\Delta t`; the - error approaches the 4th-order spatial floor - :math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}` at fine - :math:`\Delta t`. -* **Pressure gradient**: shows **order reduction** to :math:`M = 3`. - The algebraic variable :math:`G` is the Lagrange multiplier determined - at each collocation node by the Schur-complement constraint solve. - At a single node :math:`t_m` the derivative :math:`U_m = u'(t_m)` has - accuracy :math:`\mathcal{O}(\Delta t^M)` (not the superconvergent order - :math:`2M-1`, which applies only to the endpoint value via Gauss - quadrature), and :math:`G_m` inherits that same order. Hence the - pressure at the endpoint converges at order :math:`M = 3`. + 0 = B\,\mathbf{u} - q(t). + +The sweeper is +:class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE`; +the saddle-point solve at each node is a direct Schur-complement +factorisation that bypasses Newton. + +Two formulations are compared +------------------------------- +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`. + +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). **Spatial resolution**: ``nvars = 1023`` interior points with a fourth-order FD Laplacian (:math:`\Delta x = 1/1024`, spatial error floor @@ -60,14 +61,15 @@ from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.projects.DAE.sweepers.semiImplicitDAE import SemiImplicitDAE -from pySDC.playgrounds.Stokes_Poiseuille_1D_FD.Stokes_Poiseuille_1D_FD import stokes_poiseuille_1d_fd +from pySDC.playgrounds.Stokes_Poiseuille_1D_FD.Stokes_Poiseuille_1D_FD import ( + stokes_poiseuille_1d_fd, + stokes_poiseuille_1d_fd_lift, +) # --------------------------------------------------------------------------- # Fixed parameters # --------------------------------------------------------------------------- -# 4th-order FD Laplacian: spatial error ~ O(dx^4). -# With nvars=1023 (dx=1/1024), spatial floor ~ 1e-12. _NVARS = 1023 _NU = 1.0 _TEND = 0.5 @@ -86,29 +88,17 @@ # Helpers # --------------------------------------------------------------------------- -def _run(dt, restol=_RESTOL, max_iter=50): +def _run(problem_class, dt, restol=_RESTOL, max_iter=50, nvars=_NVARS): """ Run one simulation and return ``(uend, problem_instance)``. - Parameters - ---------- - dt : float - Time-step size. - restol : float - Residual tolerance for SDC convergence. - max_iter : int - Maximum number of SDC iterations per step. - - Returns - ------- - uend : MeshDAE - Numerical solution at :math:`T_\\text{end}`. - P : stokes_poiseuille_1d_fd - Problem instance (gives access to exact solutions). + For the lifted variant ``uend.diff`` contains the *lifted* velocity + :math:`\\tilde{\\mathbf{v}}`; call :meth:`P.lift` to recover the + physical velocity. """ desc = { - 'problem_class': stokes_poiseuille_1d_fd, - 'problem_params': {'nvars': _NVARS, 'nu': _NU}, + 'problem_class': problem_class, + 'problem_params': {'nvars': nvars, 'nu': _NU}, 'sweeper_class': SemiImplicitDAE, 'sweeper_params': _SWEEPER_PARAMS, 'level_params': {'restol': restol, 'dt': dt}, @@ -120,7 +110,29 @@ def _run(dt, restol=_RESTOL, max_iter=50): return uend, P -def _print_table(dts, vel_errs, pres_errs, vel_order, pres_order_str): +def _errors(uend, P): + """ + Compute ``(vel_err, pres_err)`` compared to the exact analytical solution. + + For the lifted variant the physical velocity is reconstructed before + computing the error. + """ + u_ex = P.u_exact(_TEND) + + if isinstance(P, stokes_poiseuille_1d_fd_lift): + # Recover physical velocity from lifted variable + lift at T_end. + u_phys = np.asarray(uend.diff) + P.lift(_TEND) + u_ex_phys = np.asarray(u_ex.diff) + P.lift(_TEND) + else: + u_phys = np.asarray(uend.diff) + u_ex_phys = np.asarray(u_ex.diff) + + vel_err = float(np.linalg.norm(u_phys - u_ex_phys, np.inf)) + pres_err = abs(float(uend.alg[0]) - float(u_ex.alg[0])) + return vel_err, pres_err + + +def _print_table(dts, vel_errs, pres_errs, vel_order): """Print a two-column convergence table.""" header = ( f' {"dt":>10} {"vel error":>14} {"vel ord":>8} {"exp":>4}' @@ -146,13 +158,22 @@ def _print_table(dts, vel_errs, pres_errs, vel_order, pres_order_str): ) +def _asymptotic_order(dts, errs, skip=2): + """Estimate the asymptotic convergence order.""" + orders = [] + for i in range(skip, len(dts)): + if errs[i] > 0.0 and errs[i - 1] > 0.0: + orders.append(np.log(errs[i - 1] / errs[i]) / np.log(dts[i - 1] / dts[i])) + return round(float(np.median(orders)), 1) if orders else float('nan') + + # --------------------------------------------------------------------------- # Main convergence study # --------------------------------------------------------------------------- def main(): r""" - Temporal convergence study for the Stokes/Poiseuille index-1 DAE. + Compare the standard and lifted Stokes/Poiseuille formulations. Fixed parameters: @@ -160,70 +181,78 @@ def main(): * ``nvars = 1023`` (4th-order FD, spatial floor :math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}`). * :math:`T_\text{end} = 0.5`. - - Expected collocation order :math:`2M-1 = 5` for the velocity at the - endpoint. The pressure gradient :math:`G` (Lagrange multiplier) is - predicted to converge at only order :math:`M = 3`, because the Schur - complement computes :math:`G` from the derivative :math:`U_m`, which - achieves only :math:`\mathcal{O}(\Delta t^M)` accuracy at each internal - collocation node (no super-convergence for internal nodes). """ max_order = 2 * _NUM_NODES - 1 # = 5 - - # dt range: 6 halvings from T/2 to T/64. - dts = [_TEND / (2**k) for k in range(1, 7)] # 0.25 … 0.0078 + dts = [_TEND / (2**k) for k in range(1, 7)] # 0.25 … 0.0078 print(f'\nFully-converged SDC (restol={_RESTOL:.0e}, ν={_NU}, M={_NUM_NODES})') - print(f'Sweeper: SemiImplicitDAE (DAE-specific, constraint enforced at each node)') + 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'Error vs. exact analytical solution at T={_TEND}\n') - - vel_errs = [] - pres_errs = [] - - for dt in dts: - uend, P = _run(dt) - u_ex = P.u_exact(_TEND) - vel_errs.append(float(np.linalg.norm(np.asarray(uend.diff) - np.asarray(u_ex.diff), np.inf))) - pres_errs.append(abs(float(uend.alg[0]) - float(u_ex.alg[0]))) - - # Estimate asymptotic order for pressure (skip the first point). - orders_pres = [] - for i in range(2, len(dts)): - if pres_errs[i] > 0.0 and pres_errs[i - 1] > 0.0: - orders_pres.append( - np.log(pres_errs[i - 1] / pres_errs[i]) / np.log(dts[i - 1] / dts[i]) - ) - pres_est = round(float(np.median(orders_pres)), 1) if orders_pres else float('nan') - pres_order_str = f'≈ {pres_est:.1f}' - - print('=' * 80) - _print_table(dts, vel_errs, pres_errs, max_order, pres_order_str) - print('=' * 80) - - print(f'\nSummary') - print(f' Velocity: converging to collocation order {max_order}') - print(f' Pressure gradient: observed order {pres_order_str}') - if pres_est < max_order - 0.5: - print( - f' → Order reduction in the algebraic variable: {max_order} → {pres_order_str}.' - ) - print( - ' The pressure gradient G (Lagrange multiplier) converges at order M = 3,\n' - ' not at the superconvergent velocity order 2M-1 = 5.\n' - ' Reason: the Schur complement computes G_m from the velocity derivative\n' - ' U_m = u\'(t_m), which has only O(dt^M) accuracy at internal collocation\n' - ' nodes (super-convergence order 2M-1 applies only to the endpoint value\n' - ' via the Gauss-quadrature formula, not to the nodal derivatives).' - ) - else: - print( - f' → No order reduction detected: both variables converge at order ≈ {max_order}.' - ) + print(f'Error vs. exact analytical solution at T={_TEND}') + + cases = [ + (stokes_poiseuille_1d_fd, 'Standard (B·u = q(t), time-dependent constraint)'), + (stokes_poiseuille_1d_fd_lift, 'Lifted (B·ṽ = 0, homogeneous constraint) '), + ] + + results = {} + for cls, label in cases: + print() + print('=' * 72) + print(f' {label}') + print('=' * 72) + + vel_errs, pres_errs = [], [] + for dt in dts: + uend, P = _run(cls, dt) + ve, pe = _errors(uend, P) + vel_errs.append(ve) + pres_errs.append(pe) + + _print_table(dts, vel_errs, pres_errs, max_order) + results[cls.__name__] = (vel_errs, pres_errs) + + # ---- Summary ---- + print() + print('=' * 72) + print(' Summary') + print('=' * 72) + for cls, label in cases: + 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})' + else: + status = f'order reduced to {pres_ord:.1f}' + print(f' Pressure order ≈ {pres_ord:.1f} → {status}') + else: + print(f' Pressure order ≈ {pres_ord:.1f} → full collocation 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).' + ) + print( + f' • Lifted: pressure converges at increasing order {pres_ord_lft:.1f}+' + f' (constraint B·ṽ = 0 is autonomous; order reduction removed).' + ) + print( + ' • The velocity accuracy is identical in both cases, confirming that\n' + ' the lifting only affects the algebraic variable.' + ) print( - f' (Convergence may plateau at the 4th-order spatial floor ~1e-12\n' - f' once temporal errors fall below O(dx⁴).)' + ' • Convergence may plateau at the 4th-order spatial floor ~1e-12\n' + ' once temporal errors fall below O(dx⁴) at fine Δt.' ) From 289e58188d3dfa73824e8bc387dc6d3595f34a71 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 08:06:52 +0000 Subject: [PATCH 4/9] Use nu=0.1, Tend=1.0 for clean M+1=4 convergence orders; document U-formulation order limit Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/0a7bcad2-2358-4ceb-b9be-ee9049ffa9cf --- .../run_convergence.py | 156 ++++++++++++------ 1 file changed, 108 insertions(+), 48 deletions(-) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py index 53031f660e..1537ce41c6 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py +++ b/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.' ) From 6d21c46cdef5cfacebdcd812ba246c2e8d12a2e0 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 09:04:40 +0000 Subject: [PATCH 5/9] Add FullyImplicitDAE-compatible Stokes classes; update run_convergence.py with 4-variant comparison Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/caaa10d2-854c-4d03-8f8f-986aa9a48220 --- .../Stokes_Poiseuille_1D_FD.py | 266 +++++++++++++++++- .../run_convergence.py | 149 +++++----- 2 files changed, 339 insertions(+), 76 deletions(-) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py index 451f7c91b6..9aedb5c012 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -52,13 +52,21 @@ * ``u.diff[:]`` – velocity on :math:`N` interior grid points. * ``u.alg[0]`` – pressure gradient :math:`G` (Lagrange multiplier). -The compatible sweeper is -:class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE`. +Two sweepers are supported: + +* :class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE` + (U-formulation): velocity order :math:`M+1`, pressure order :math:`M` + (or increasing toward :math:`M+1` with lifting). + +* :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE` + (collocation-consistent): treats :math:`G` as a differential variable + updated by the same quadrature as the velocity, restoring the full + collocation order :math:`2M-1` for both. Without constraint lifting (class :class:`stokes_poiseuille_1d_fd`), the constraint :math:`B\mathbf{u} = q(t)` has a time-dependent right-hand side, which causes order reduction in :math:`G` to order :math:`M` (= 3 for -3 RADAU-RIGHT nodes). +3 RADAU-RIGHT nodes) when using :class:`SemiImplicitDAE`. Constraint lifting (class :class:`stokes_poiseuille_1d_fd_lift`) ----------------------------------------------------------------- @@ -95,11 +103,13 @@ ------- stokes_poiseuille_1d_fd No lifting; constraint :math:`B\mathbf{u} = q(t)` (time-dependent). - Pressure converges at order :math:`M`. + Compatible with **SemiImplicitDAE** (pressure order :math:`M`) and + **FullyImplicitDAE** (pressure order :math:`2M-1`). stokes_poiseuille_1d_fd_lift Constraint lifting; homogeneous :math:`B\tilde{\mathbf{v}} = 0`. - Expected to restore full order :math:`2M-1` in the pressure. + Compatible with both sweepers; FullyImplicitDAE restores full + :math:`2M-1` order cleanly. """ import numpy as np @@ -307,6 +317,72 @@ def _schur_solve(self, rhs_eff, v_approx, factor, constraint_rhs): U = w + G_new * v0 return U, float(G_new) + def _schur_solve_full_implicit(self, rhs_eff, v_approx, factor, constraint_rhs): + r""" + Schur-complement saddle-point solve for use with + :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE`. + + Finds :math:`(\mathbf{U}, G')` satisfying + + .. math:: + + (I - \alpha\nu A)\,\mathbf{U} - \alpha G'\,\mathbf{1} + = \mathbf{r}_\text{eff}, + + .. math:: + + B(\mathbf{v}_\text{approx} + \alpha\,\mathbf{U}) = c, + + where :math:`G' = \mathrm{d}G/\mathrm{d}t` is the **derivative** of + the pressure gradient and :math:`c` is ``constraint_rhs``. + + Compared to :meth:`_schur_solve`, here ``rhs_eff`` must already + include the :math:`G_0\,\mathbf{1}` term (the current pressure + estimate from ``u_approx.alg[0]``), and the denominator carries an + extra :math:`\alpha` factor: + + .. math:: + + G' = \frac{c - B\mathbf{v} - \alpha B\mathbf{w}} + {\alpha^2 B\mathbf{v}_0}, + + Parameters + ---------- + rhs_eff : numpy.ndarray + Effective velocity RHS :math:`\nu A\mathbf{v} + G_0\mathbf{1} + + \mathbf{f}_\text{net}`. + v_approx : numpy.ndarray + Current velocity approximation at the node. + factor : float + Implicit prefactor :math:`\alpha = \Delta t\,\tilde{q}_{mm}`. + constraint_rhs : float + RHS of the constraint :math:`c` (``q(t)`` standard; ``0`` lifted). + + Returns + ------- + U : numpy.ndarray + Velocity derivative at the node. + G_prime : float + Derivative of the pressure gradient :math:`G'`. + """ + M = sp.eye(self.nvars, format='csc') - factor * self.A + w = spsolve(M, rhs_eff) + v0 = spsolve(M, self.ones) + + Bw = self._B_dot(w) + Bv0 = self._B_dot(v0) + Bv = self._B_dot(v_approx) + # Extra factor of alpha in denominator vs _schur_solve + denom = factor**2 * Bv0 + assert abs(denom) > 0.0, ( + f'_schur_solve_full_implicit: denominator factor²·B·v₀ = {denom:.3e} is zero; ' + f'factor = {factor}, B·v₀ = {Bv0:.3e}' + ) + G_prime = (constraint_rhs - Bv - factor * Bw) / denom + + U = w + factor * G_prime * v0 + return U, float(G_prime) + # --------------------------------------------------------------------------- # Case 1: No lifting – constraint B·u = q(t) @@ -617,3 +693,183 @@ def du_exact(self, t): me.diff[:] = np.sin(np.pi * self.xvalues) * np.cos(t) - self._lift_dot(t) me.alg[0] = -np.sin(t) return me + + +# --------------------------------------------------------------------------- +# Case 3: FullyImplicitDAE – no lifting +# --------------------------------------------------------------------------- + +class stokes_poiseuille_1d_fd_full(stokes_poiseuille_1d_fd): + r""" + 1-D Stokes/Poiseuille DAE **without** constraint lifting, compatible with + :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE`. + + The key difference from :class:`stokes_poiseuille_1d_fd` is in + ``solve_system``: here the unknown is the full derivative + :math:`(\mathbf{U}, G') = (\mathbf{u}', G')`, i.e. ``me.alg[0]`` is the + *derivative* of the pressure gradient. The pressure gradient at the + node is then recovered by quadrature: + + .. math:: + + G_m = G_0 + \Delta t \sum_{j=1}^{M} Q_{mj}\,G'_j, + + which gives the pressure the same collocation structure as + the velocity (both recovered via the quadrature formula). + + .. note:: + + In practice, :class:`FullyImplicitDAE` and + :class:`~.semiImplicitDAE.SemiImplicitDAE` converge to the + **same collocation fixed point** for this index-1 DAE: the + :math:`\mathcal{O}(\Delta t^M)` errors in the stage pressure + values break the :math:`2M-1` superconvergence, and both sweepers + achieve velocity order :math:`M+1` and pressure order :math:`M` + (standard) or increasing toward :math:`M+1` (lifted). The + :class:`FullyImplicitDAE` formulation is included for completeness + and pedagogical comparison. + + The Schur-complement solve for the unknown :math:`(\mathbf{U}, G')`: + + .. math:: + + (I - \alpha\nu A)\,\mathbf{U} - \alpha G'\,\mathbf{1} + = \nu A\mathbf{v} + G_0\mathbf{1} + \mathbf{f}(t), + + .. math:: + + B(\mathbf{v} + \alpha\,\mathbf{U}) = q(t), + + yields + + .. math:: + + G' = \frac{q(t) - B\mathbf{v} - \alpha B\mathbf{w}}{\alpha^2 B\mathbf{v}_0}, + \quad + \mathbf{U} = \mathbf{w} + \alpha G'\,\mathbf{v}_0, + + where :math:`\mathbf{w} = (I-\alpha\nu A)^{-1}(\nu A\mathbf{v} + + G_0\mathbf{1} + \mathbf{f})` and + :math:`\mathbf{v}_0 = (I-\alpha\nu A)^{-1}\mathbf{1}`. + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity (default 1.0). + newton_tol : float + Unused; passed to base class (default 1e-10). + """ + + def solve_system(self, impl_sys, u_approx, factor, u0, t): + r""" + Schur-complement solve for + :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE`. + + Returns the **derivative** :math:`(\mathbf{U}, G')`. + ``me.alg[0]`` = :math:`G'` (pressure-gradient time derivative). + + Parameters + ---------- + impl_sys : callable + Unused; system solved directly. + u_approx : MeshDAE + Approximation :math:`(\mathbf{v}, G_0)` at the current node. + factor : float + Implicit prefactor :math:`\alpha`. + u0 : MeshDAE + Unused (direct solver). + t : float + Current time. + + Returns + ------- + me : MeshDAE + ``me.diff[:]`` = :math:`\mathbf{U}`, + ``me.alg[0]`` = :math:`G'`. + """ + me = self.dtype_u(self.init, val=0.0) + v_approx = np.asarray(u_approx.diff).copy() + G0 = float(u_approx.alg[0]) + + # rhs_eff includes current G0: FullyImplicitDAE treats G as a differential + # variable (G = G0 + factor*G'), so G0 enters the RHS (unlike SemiImplicitDAE + # where u_approx.alg = 0 and G is purely algebraic in the local solve). + rhs_eff = self.A.dot(v_approx) + G0 * self.ones + self._forcing(t) + U, G_prime = self._schur_solve_full_implicit(rhs_eff, v_approx, factor, self._q(t)) + + me.diff[:] = U + me.alg[0] = G_prime + return me + + +# --------------------------------------------------------------------------- +# Case 4: FullyImplicitDAE – with lifting +# --------------------------------------------------------------------------- + +class stokes_poiseuille_1d_fd_lift_full(stokes_poiseuille_1d_fd_lift): + r""" + 1-D Stokes/Poiseuille DAE **with** constraint lifting, compatible with + :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE`. + + Combines the homogeneous constraint :math:`B\tilde{\mathbf{v}} = 0` + from :class:`stokes_poiseuille_1d_fd_lift` with the + :class:`FullyImplicitDAE`-consistent ``solve_system`` from + :class:`stokes_poiseuille_1d_fd_full`. + + .. note:: + + :class:`FullyImplicitDAE` and :class:`SemiImplicitDAE` converge + to the same collocation fixed point for this DAE: velocity order + :math:`M+1` and pressure order increasing toward :math:`M+1` + (same as :class:`stokes_poiseuille_1d_fd_lift`). This class is + provided for pedagogical comparison. + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity (default 1.0). + newton_tol : float + Unused; passed to base class (default 1e-10). + """ + + def solve_system(self, impl_sys, v_approx_mesh, factor, u0, t): + r""" + Schur-complement solve for + :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE` + with the **homogeneous** constraint :math:`B\tilde{\mathbf{v}} = 0`. + + Returns :math:`(\tilde{\mathbf{U}}, G')`. + + Parameters + ---------- + impl_sys : callable + Unused; system solved directly. + v_approx_mesh : MeshDAE + Approximation :math:`(\tilde{\mathbf{v}}, G_0)` at the node. + factor : float + Implicit prefactor :math:`\alpha`. + u0 : MeshDAE + Unused (direct solver). + t : float + Current time. + + Returns + ------- + me : MeshDAE + ``me.diff[:]`` = :math:`\tilde{\mathbf{U}}`, + ``me.alg[0]`` = :math:`G'`. + """ + me = self.dtype_u(self.init, val=0.0) + vv = np.asarray(v_approx_mesh.diff).copy() + G0 = float(v_approx_mesh.alg[0]) + + rhs_eff = self.A.dot(vv) + G0 * self.ones + self._net_forcing(t) + U, G_prime = self._schur_solve_full_implicit(rhs_eff, vv, factor, 0.0) + + me.diff[:] = U + me.alg[0] = G_prime + return me diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py index 1537ce41c6..a551b06100 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -14,54 +14,54 @@ 0 = B\,\mathbf{u} - q(t). -The sweeper is -:class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE`; -the saddle-point solve at each node is a direct Schur-complement -factorisation that bypasses Newton. - -Two formulations are compared -------------------------------- -1. **Standard** (:class:`~.Stokes_Poiseuille_1D_FD.stokes_poiseuille_1d_fd`) - — constraint :math:`B\mathbf{u} = q(t)` has a time-dependent RHS. - 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`) - — 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. +Two sweepers are compared +-------------------------- +1. :class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE` + (U-formulation): stores and integrates velocity derivatives + :math:`U_m = u'(\tau_m)`. +2. :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE` + (collocation-consistent): treats the full state :math:`(u, G)` as a + differential system with derivative :math:`(U_m, G'_m)`. + +Two constraint formulations are compared +----------------------------------------- +1. **Standard** (no lifting): constraint :math:`B\mathbf{u} = q(t)` has a + time-dependent RHS. Causes order reduction in :math:`G` to order :math:`M`. +2. **Lifted**: lifting :math:`\mathbf{u}_\ell(t) = (q(t)/s)\,\mathbf{1}` makes + the constraint **homogeneous**: :math:`B\tilde{\mathbf{v}} = 0`. 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`: +polynomial achieves the superconvergent order :math:`2M-1` at the endpoint. + +For this **DAE**, both SemiImplicitDAE and FullyImplicitDAE converge to the +**same collocation fixed point**. At that fixed point the stage values of the +pressure gradient :math:`G_m` have only :math:`\mathcal{O}(\Delta t^M)` +accuracy (the constraint at each node is exact, but the velocity +:math:`u_m = u_0 + \Delta t\sum_j Q_{mj} U_j` at intermediate nodes has +only stage-order accuracy :math:`\mathcal{O}(\Delta t^{M+1})`). These +:math:`\mathcal{O}(\Delta t^M)` errors in :math:`G_m` feed back into the +velocity derivatives :math:`U_m = A u_m + G_m \mathbf{1} + f(\tau_m)`, +**breaking the superconvergence** and limiting: + +* velocity to :math:`\mathcal{O}(\Delta t^{M+1})` — one higher than the + stage order; +* pressure :math:`G` to :math:`\mathcal{O}(\Delta t^M)` (standard) or + increasing toward :math:`M+1` (lifted). + +This is confirmed across :math:`M = 2, 3, 4` RADAU-RIGHT and is independent +of the sweeper (SemiImplicitDAE = FullyImplicitDAE at the fixed point): * :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`) +Achieving the full collocation order :math:`2M-1` for **both** velocity and +pressure would require solving all :math:`M` RADAU stages simultaneously as +a coupled system (standard RADAU-IIA), which goes beyond the node-by-node +SDC sweep used here. + Observed results (:math:`\nu = 0.1`, ``nvars = 1023``, ``restol = 1e-13``, :math:`M = 3`) --------------------------------------------------------------------------- @@ -78,6 +78,8 @@ * **Lifted**: velocity at :math:`M+1 = 4` (unchanged); pressure order increases monotonically, approaching :math:`M+1 = 4` (homogeneous constraint removes the order reduction). +* **FullyImplicitDAE** (standard or lifted): **identical results** to + SemiImplicitDAE — both converge to the same collocation fixed point. **Spatial resolution**: ``nvars = 1023`` interior points with a fourth-order FD Laplacian (:math:`\Delta x = 1/1024`, spatial error floor @@ -92,9 +94,12 @@ from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.projects.DAE.sweepers.semiImplicitDAE import SemiImplicitDAE +from pySDC.projects.DAE.sweepers.fullyImplicitDAE import FullyImplicitDAE from pySDC.playgrounds.Stokes_Poiseuille_1D_FD.Stokes_Poiseuille_1D_FD import ( stokes_poiseuille_1d_fd, stokes_poiseuille_1d_fd_lift, + stokes_poiseuille_1d_fd_full, + stokes_poiseuille_1d_fd_lift_full, ) # --------------------------------------------------------------------------- @@ -124,7 +129,7 @@ # Helpers # --------------------------------------------------------------------------- -def _run(problem_class, dt, restol=_RESTOL, max_iter=50, nvars=_NVARS): +def _run(problem_class, sweeper_class, dt, restol=_RESTOL, max_iter=50, nvars=_NVARS): """ Run one simulation and return ``(uend, problem_instance)``. @@ -135,7 +140,7 @@ def _run(problem_class, dt, restol=_RESTOL, max_iter=50, nvars=_NVARS): desc = { 'problem_class': problem_class, 'problem_params': {'nvars': nvars, 'nu': _NU}, - 'sweeper_class': SemiImplicitDAE, + 'sweeper_class': sweeper_class, 'sweeper_params': _SWEEPER_PARAMS, 'level_params': {'restol': restol, 'dt': dt}, 'step_params': {'maxiter': max_iter}, @@ -155,7 +160,7 @@ def _errors(uend, P): """ u_ex = P.u_exact(_TEND) - if isinstance(P, stokes_poiseuille_1d_fd_lift): + if hasattr(P, 'lift'): # Recover physical velocity from lifted variable + lift at T_end. u_phys = np.asarray(uend.diff) + P.lift(_TEND) u_ex_phys = np.asarray(u_ex.diff) + P.lift(_TEND) @@ -209,7 +214,7 @@ def _asymptotic_order(dts, errs, skip=2): def main(): r""" - Compare the standard and lifted Stokes/Poiseuille formulations. + Compare four formulations: two sweepers × two constraint treatments. Fixed parameters: @@ -220,22 +225,23 @@ def main(): Expected orders (see module docstring for derivation): - * **Velocity**: :math:`M+1 = 4` for both formulations. + * **Velocity**: :math:`M+1 = 4` for all four 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). + * **SemiImplicitDAE vs FullyImplicitDAE**: identical results — both + converge to the same collocation fixed point for this DAE. """ - vel_order = _NUM_NODES + 1 # M+1 = 4 (U-formulation of SemiImplicitDAE) + vel_order = _NUM_NODES + 1 # M+1 = 4 (U-formulation limit) 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 (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)') + f'(pure-ODE collocation order 2M-1 = {2*_NUM_NODES-1} not achieved; see module docstring)') 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)') @@ -244,12 +250,18 @@ def main(): print(f'Error vs. exact analytical solution at T={_TEND}') cases = [ - (stokes_poiseuille_1d_fd, 'Standard (B·u = q(t), time-dependent constraint)'), - (stokes_poiseuille_1d_fd_lift, 'Lifted (B·ṽ = 0, homogeneous constraint) '), + (stokes_poiseuille_1d_fd, SemiImplicitDAE, + 'SemiImplicitDAE, standard (B·u = q(t))'), + (stokes_poiseuille_1d_fd_lift, SemiImplicitDAE, + 'SemiImplicitDAE, lifted (B·ṽ = 0) '), + (stokes_poiseuille_1d_fd_full, FullyImplicitDAE, + 'FullyImplicitDAE, standard (B·u = q(t)) [same fixed pt as SemiImplicit]'), + (stokes_poiseuille_1d_fd_lift_full, FullyImplicitDAE, + 'FullyImplicitDAE, lifted (B·ṽ = 0) [same fixed pt as SemiImplicit]'), ] results = {} - for cls, label in cases: + for cls, sweeper, label in cases: print() print('=' * 72) print(f' {label}') @@ -257,26 +269,26 @@ def main(): vel_errs, pres_errs = [], [] for dt in dts: - uend, P = _run(cls, dt) + uend, P = _run(cls, sweeper, dt) ve, pe = _errors(uend, P) vel_errs.append(ve) pres_errs.append(pe) _print_table(dts, vel_errs, pres_errs, vel_order) - results[cls.__name__] = (vel_errs, pres_errs) + results[label] = (vel_errs, pres_errs) # ---- Summary ---- print() print('=' * 72) print(' Summary') print('=' * 72) - for cls, label in cases: - vel_errs, pres_errs = results[cls.__name__] + for cls, sweeper, label in cases: + vel_errs, pres_errs = results[label] vel_ord = _asymptotic_order(dts, vel_errs) pres_ord = _asymptotic_order(dts, pres_errs) - is_lift = isinstance(cls(), stokes_poiseuille_1d_fd_lift) + is_lift = hasattr(cls(), 'lift') exp_pres = vel_order if is_lift else pres_order - print(f'\n {cls.__name__}:') + print(f'\n {label}:') print(f' Velocity order ≈ {vel_ord:.1f} (expected M+1 = {vel_order})') if pres_ord < vel_order - 0.4: if is_lift: @@ -289,32 +301,27 @@ def main(): 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' • Velocity order M+1 = {vel_order} confirmed for both formulations.' + f' • All four formulations confirm velocity order M+1 = {vel_order}.' ) print( - f' • Standard: pressure at order {pres_ord_std:.1f} = M ' - f'(order reduction from time-dependent constraint).' + f' • Standard: pressure at order M = {pres_order} (order reduction from time-dependent constraint).' ) print( - f' • Lifted: pressure at increasing order {pres_ord_lft:.1f}+' - f' (heading to M+1 = {vel_order}; autonomous constraint removes reduction).' + f' • Lifted: pressure increasing toward M+1 = {vel_order} (autonomous constraint removes reduction).' ) print( - ' • 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.)' + ' • SemiImplicitDAE and FullyImplicitDAE converge to the SAME collocation\n' + ' fixed point for this DAE: identical velocity and pressure errors.\n' + ' The O(dt^M) stage pressure errors break the 2M-1 superconvergence\n' + ' of both sweepers; achieving 2M-1 would require solving all M RADAU\n' + ' stages simultaneously (full RADAU-IIA, beyond SDC node-by-node sweeps).' ) print( - 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.' + f' • Spatial floor ~1e-12 (nvars={_NVARS}) may cause order plateau at fine Δt.' ) if __name__ == '__main__': main() + From 870f3424f1dceab08e88c01603de641b7f02e6de Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 10:02:07 +0000 Subject: [PATCH 6/9] Add differentiated-constraint classes achieving full 2M-1 RADAU order for velocity and pressure Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/93f40bfb-629c-47d8-8e13-7f5d9f451f22 --- .../Stokes_Poiseuille_1D_FD.py | 341 +++++++++++++++++- .../run_convergence.py | 225 ++++++------ 2 files changed, 428 insertions(+), 138 deletions(-) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py index 9aedb5c012..68a428cc3e 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -56,17 +56,15 @@ * :class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE` (U-formulation): velocity order :math:`M+1`, pressure order :math:`M` - (or increasing toward :math:`M+1` with lifting). + (standard), approaching :math:`M+1` (lifted), or full :math:`2M-1` + (differentiated constraint). -* :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE` - (collocation-consistent): treats :math:`G` as a differential variable - updated by the same quadrature as the velocity, restoring the full - collocation order :math:`2M-1` for both. +Three constraint treatments are provided — see Classes section below. Without constraint lifting (class :class:`stokes_poiseuille_1d_fd`), the constraint :math:`B\mathbf{u} = q(t)` has a time-dependent right-hand side, which causes order reduction in :math:`G` to order :math:`M` (= 3 for -3 RADAU-RIGHT nodes) when using :class:`SemiImplicitDAE`. +3 RADAU-RIGHT nodes). Constraint lifting (class :class:`stokes_poiseuille_1d_fd_lift`) ----------------------------------------------------------------- @@ -95,21 +93,47 @@ - \dot{\mathbf{u}}_\ell(t)\bigr] + G\,\mathbf{1}. -Because the constraint :math:`B\tilde{\mathbf{v}} = 0` no longer contains a -time-dependent right-hand side, the Lagrange multiplier :math:`G` is expected -to converge at the full collocation order :math:`2M-1`, matching the velocity. +Because the constraint :math:`B\tilde{\mathbf{v}} = 0` is autonomous, the +order reduction in :math:`G` is reduced (pressure approaches :math:`M+1`) +but not fully eliminated in the SDC context. + +Differentiated constraint (class :class:`stokes_poiseuille_1d_fd_diffconstr`) +------------------------------------------------------------------------------ +Instead of enforcing the algebraic constraint +:math:`B\mathbf{u}_m = q(\tau_m)` at each SDC stage, enforce its time +derivative :math:`B\mathbf{U}_m = q'(\tau_m)`. The Schur formula becomes + +.. math:: + + G_m = \frac{q'(\tau_m) - B\mathbf{w}}{B\mathbf{v}_0}, + +and the stage pressure error reduces from :math:`\mathcal{O}(\Delta t^M)` +to :math:`\mathcal{O}(\Delta t^{M+1})`, restoring the full RADAU +collocation order :math:`2M-1` for **both** velocity and pressure. Classes ------- stokes_poiseuille_1d_fd - No lifting; constraint :math:`B\mathbf{u} = q(t)` (time-dependent). - Compatible with **SemiImplicitDAE** (pressure order :math:`M`) and - **FullyImplicitDAE** (pressure order :math:`2M-1`). + No lifting; algebraic constraint :math:`B\mathbf{u} = q(t)`. + SemiImplicitDAE: vel :math:`M+1`, pres :math:`M`. stokes_poiseuille_1d_fd_lift Constraint lifting; homogeneous :math:`B\tilde{\mathbf{v}} = 0`. - Compatible with both sweepers; FullyImplicitDAE restores full - :math:`2M-1` order cleanly. + SemiImplicitDAE: vel :math:`M+1`, pres approaching :math:`M+1`. + +stokes_poiseuille_1d_fd_diffconstr + Differentiated constraint :math:`B\mathbf{U}_m = q'(\tau_m)`. + SemiImplicitDAE: **vel and pres both at** :math:`2M-1`. + +stokes_poiseuille_1d_fd_lift_diffconstr + Lifting + differentiated :math:`B\tilde{\mathbf{U}} = 0`. + Equivalent to :class:`stokes_poiseuille_1d_fd_lift` at the fixed point. + +stokes_poiseuille_1d_fd_full + No lifting, FullyImplicitDAE (same fixed point as standard). + +stokes_poiseuille_1d_fd_lift_full + Lifting, FullyImplicitDAE (same fixed point as lifted). """ import numpy as np @@ -247,6 +271,12 @@ def _q(self, t): """ return self.C_B * np.sin(t) + def _q_prime(self, t): + r""" + Time derivative of the flow-rate RHS: :math:`q'(t) = C_B\,\cos(t)`. + """ + return self.C_B * np.cos(t) + def _forcing(self, t): r""" Manufactured forcing consistent with :math:`u_\text{ex}` and @@ -383,6 +413,68 @@ def _schur_solve_full_implicit(self, rhs_eff, v_approx, factor, constraint_rhs): U = w + factor * G_prime * v0 return U, float(G_prime) + def _schur_solve_diffconstr(self, rhs_eff, factor, q_prime_val): + r""" + Schur-complement solve using the **differentiated constraint** + :math:`B\mathbf{U} = q'(t)`. + + Instead of enforcing the original algebraic constraint + :math:`B(\mathbf{v} + \alpha\mathbf{U}) = q(t)`, this method uses + the differentiated form :math:`B\mathbf{U} = q'(t)`. The key + formula is + + .. math:: + + G = \frac{q'(t) - B\mathbf{w}}{B\mathbf{v}_0}, + + where + :math:`\mathbf{w} = (I - \alpha\nu A)^{-1}\mathbf{r}_\text{eff}` and + :math:`\mathbf{v}_0 = (I - \alpha\nu A)^{-1}\mathbf{1}`. + + **Why this gives higher-order pressure**: At the SDC fixed point the + stage velocities satisfy :math:`\mathbf{u}_m - \mathbf{u}(\tau_m) = + \mathcal{O}(\Delta t^{M+1})` (collocation accuracy). The + differentiated-constraint error propagates as + + .. math:: + + e_{G_m} = G_m - G(\tau_m) = -\frac{B A\,e_{\mathbf{u}_m}}{s} + = \mathcal{O}(\Delta t^{M+1}), + + whereas the original algebraic constraint gives only + :math:`\mathcal{O}(\Delta t^M)`. + + Parameters + ---------- + rhs_eff : numpy.ndarray + Effective velocity RHS + :math:`\nu A\mathbf{v}_\text{approx} + \mathbf{f}_\text{net}(t)`. + factor : float + Implicit prefactor :math:`\alpha = \Delta t\,\tilde{q}_{mm}`. + q_prime_val : float + Value of :math:`q'(t)` at the current stage time. + + Returns + ------- + U : numpy.ndarray + Velocity derivative at the node. + G_new : float + Pressure gradient satisfying :math:`B\mathbf{U} = q'(t)`. + """ + M_mat = sp.eye(self.nvars, format='csc') - factor * self.A + w = spsolve(M_mat, rhs_eff) + v0 = spsolve(M_mat, self.ones) + + Bw = self._B_dot(w) + Bv0 = self._B_dot(v0) + assert abs(Bv0) > 0.0, ( + f'_schur_solve_diffconstr: B·v₀ = {Bv0:.3e} is zero; factor = {factor}' + ) + G_new = (q_prime_val - Bw) / Bv0 + + U = w + G_new * v0 + return U, float(G_new) + # --------------------------------------------------------------------------- # Case 1: No lifting – constraint B·u = q(t) @@ -873,3 +965,224 @@ def solve_system(self, impl_sys, v_approx_mesh, factor, u0, t): me.diff[:] = U me.alg[0] = G_prime return me + + +# --------------------------------------------------------------------------- +# Case 5: No lifting – differentiated constraint B·U = q'(t) +# --------------------------------------------------------------------------- + +class stokes_poiseuille_1d_fd_diffconstr(stokes_poiseuille_1d_fd): + r""" + 1-D Stokes/Poiseuille DAE using the **differentiated constraint** + :math:`B\mathbf{u}' = q'(t)` at every SDC stage. + + This is a remedy for the order reduction in the pressure gradient. + Rather than enforcing the original algebraic constraint + :math:`B\mathbf{u}_m = q(\tau_m)` (which limits :math:`G` to order + :math:`M`), each stage solve uses + + .. math:: + + B\,\mathbf{U}_m = q'(\tau_m), + + where :math:`\mathbf{U}_m = \mathbf{u}'(\tau_m)` is the velocity + derivative. The corresponding Schur-complement formula is + + .. math:: + + G_m = \frac{q'(\tau_m) - B\mathbf{w}}{B\mathbf{v}_0}, + + giving pressure error + :math:`e_{G_m} = -BA\,e_{\mathbf{u}_m}/s = \mathcal{O}(\Delta t^{M+1})` + (one order higher than the algebraic constraint). + + ``eval_f`` is also modified to check :math:`B\mathbf{u}' - q'(t) = 0` + (so that the SDC residual converges to machine precision at the fixed + point of the differentiated-constraint iteration). + + .. note:: + + The original constraint :math:`B\mathbf{u} = q(t)` is satisfied + approximately: since :math:`B\mathbf{u}(t_n) = q(t_n)` (consistent + IC) and :math:`B\mathbf{U}_m = q'(\tau_m)` at every stage, + the quadrature formula gives + :math:`B\mathbf{u}_{n+1} - q(t_{n+1}) = \mathcal{O}(\Delta t^{2M})` + at the endpoint (Gauss quadrature error) and + :math:`\mathcal{O}(\Delta t^{M+1})` at interior nodes. + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity (default 1.0). + newton_tol : float + Unused; passed to base class (default 1e-10). + """ + + def eval_f(self, u, du, t): + r""" + Fully-implicit DAE residual using the **differentiated constraint**: + + .. math:: + + F_\text{diff} = \mathbf{u}' - \nu A\,\mathbf{u} + - G\,\mathbf{1} - \mathbf{f}(t), + + .. math:: + + F_\text{alg} = B\,\mathbf{u}' - q'(t). + + The second equation uses the velocity **derivative** ``du.diff`` + instead of the velocity state ``u.diff``, making the SDC residual + exactly zero (machine precision) when :meth:`solve_system` has + enforced :math:`B\mathbf{U}_m = q'(\tau_m)`. + """ + f = self.dtype_f(self.init, val=0.0) + u_vel = np.asarray(u.diff) + du_vel = np.asarray(du.diff) + G = float(u.alg[0]) + + f.diff[:] = du_vel - (self.A.dot(u_vel) + G * self.ones + self._forcing(t)) + f.alg[0] = self._B_dot(du_vel) - self._q_prime(t) + + self.work_counters['rhs']() + return f + + def solve_system(self, impl_sys, u_approx, factor, u0, t): + r""" + Schur-complement solve using the differentiated constraint + :math:`B\mathbf{U} = q'(t)`. + + Parameters + ---------- + impl_sys : callable + Unused; system solved directly. + u_approx : MeshDAE + Current velocity approximation at the node. + factor : float + Implicit prefactor :math:`\alpha`. + u0 : MeshDAE + Unused (direct solver). + t : float + Current time. + + Returns + ------- + me : MeshDAE + ``me.diff[:]`` = velocity derivative :math:`\mathbf{U}_m` + satisfying :math:`B\mathbf{U}_m = q'(t)`, + ``me.alg[0]`` = pressure gradient :math:`G_m`. + """ + me = self.dtype_u(self.init, val=0.0) + v_approx = np.asarray(u_approx.diff).copy() + + rhs_eff = self.A.dot(v_approx) + self._forcing(t) + U, G_new = self._schur_solve_diffconstr(rhs_eff, factor, self._q_prime(t)) + + me.diff[:] = U + me.alg[0] = G_new + return me + + +# --------------------------------------------------------------------------- +# Case 6: Lifting + differentiated constraint B·Ũ = 0 +# --------------------------------------------------------------------------- + +class stokes_poiseuille_1d_fd_lift_diffconstr(stokes_poiseuille_1d_fd_lift): + r""" + 1-D Stokes/Poiseuille DAE with **constraint lifting** and the + **differentiated constraint** :math:`B\tilde{\mathbf{u}}' = 0`. + + Combines: + + * The homogeneous constraint :math:`B\tilde{\mathbf{v}} = 0` from + :class:`stokes_poiseuille_1d_fd_lift` (reduces order reduction). + * The differentiated constraint :math:`B\tilde{\mathbf{U}} = 0` in + the stage solve, analogous to :class:`stokes_poiseuille_1d_fd_diffconstr`. + + .. note:: + + For the lifted problem the original constraint is + :math:`B\tilde{\mathbf{v}} = 0` (homogeneous, constant in time). + Its time derivative is :math:`B\tilde{\mathbf{v}}' = 0`, which is + the same condition. The differentiated-constraint Schur solve + therefore reduces to + :math:`G = -B\mathbf{w}_\text{net} / (B\mathbf{v}_0)`, + which is **equivalent to the original lifted Schur solve at the fixed + point** (when :math:`B\tilde{\mathbf{v}} \approx 0`). Both this class + and :class:`stokes_poiseuille_1d_fd_lift` therefore converge to the + same fixed point and give identical convergence orders. This class is + included for completeness and to confirm the equivalence. + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity (default 1.0). + newton_tol : float + Unused; passed to base class (default 1e-10). + """ + + def eval_f(self, u, du, t): + r""" + Fully-implicit DAE residual using the differentiated homogeneous + constraint :math:`B\tilde{\mathbf{u}}' = 0`: + + .. math:: + + F_\text{diff} = \tilde{\mathbf{u}}' - \nu A\,\tilde{\mathbf{v}} + - G\,\mathbf{1} - \mathbf{f}_\text{net}(t), + + .. math:: + + F_\text{alg} = B\,\tilde{\mathbf{u}}' - 0. + """ + f = self.dtype_f(self.init, val=0.0) + v_tilde = np.asarray(u.diff) + dv_tilde = np.asarray(du.diff) + G = float(u.alg[0]) + + f.diff[:] = dv_tilde - (self.A.dot(v_tilde) + G * self.ones + self._net_forcing(t)) + f.alg[0] = self._B_dot(dv_tilde) # B·ũ' = 0 + + self.work_counters['rhs']() + return f + + def solve_system(self, impl_sys, u_approx, factor, u0, t): + r""" + Schur-complement solve using the differentiated homogeneous + constraint :math:`B\tilde{\mathbf{U}} = 0`. + + Parameters + ---------- + impl_sys : callable + Unused; system solved directly. + u_approx : MeshDAE + Current velocity approximation :math:`(\tilde{\mathbf{v}}, G)`. + factor : float + Implicit prefactor :math:`\alpha`. + u0 : MeshDAE + Unused (direct solver). + t : float + Current time. + + Returns + ------- + me : MeshDAE + ``me.diff[:]`` = lifted velocity derivative + :math:`\tilde{\mathbf{U}}_m` (satisfies + :math:`B\tilde{\mathbf{U}}_m = 0`), + ``me.alg[0]`` = pressure gradient :math:`G_m`. + """ + me = self.dtype_u(self.init, val=0.0) + v_approx = np.asarray(u_approx.diff).copy() + + rhs_eff = self.A.dot(v_approx) + self._net_forcing(t) + # Differentiated homogeneous constraint: B·U_tilde = 0 → q'_eff = 0 + U, G_new = self._schur_solve_diffconstr(rhs_eff, factor, 0.0) + + me.diff[:] = U + me.alg[0] = G_new + return me diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py index a551b06100..268609d002 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -14,76 +14,67 @@ 0 = B\,\mathbf{u} - q(t). -Two sweepers are compared --------------------------- -1. :class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE` - (U-formulation): stores and integrates velocity derivatives - :math:`U_m = u'(\tau_m)`. -2. :class:`~pySDC.projects.DAE.sweepers.fullyImplicitDAE.FullyImplicitDAE` - (collocation-consistent): treats the full state :math:`(u, G)` as a - differential system with derivative :math:`(U_m, G'_m)`. - -Two constraint formulations are compared +Three constraint treatments are compared ----------------------------------------- -1. **Standard** (no lifting): constraint :math:`B\mathbf{u} = q(t)` has a - time-dependent RHS. Causes order reduction in :math:`G` to order :math:`M`. -2. **Lifted**: lifting :math:`\mathbf{u}_\ell(t) = (q(t)/s)\,\mathbf{1}` makes - the constraint **homogeneous**: :math:`B\tilde{\mathbf{v}} = 0`. - -Why M+1 (not 2M-1) for the velocity? --------------------------------------- -For a pure ODE discretised with RADAU-RIGHT :math:`M` nodes, the collocation -polynomial achieves the superconvergent order :math:`2M-1` at the endpoint. - -For this **DAE**, both SemiImplicitDAE and FullyImplicitDAE converge to the -**same collocation fixed point**. At that fixed point the stage values of the -pressure gradient :math:`G_m` have only :math:`\mathcal{O}(\Delta t^M)` -accuracy (the constraint at each node is exact, but the velocity -:math:`u_m = u_0 + \Delta t\sum_j Q_{mj} U_j` at intermediate nodes has -only stage-order accuracy :math:`\mathcal{O}(\Delta t^{M+1})`). These -:math:`\mathcal{O}(\Delta t^M)` errors in :math:`G_m` feed back into the -velocity derivatives :math:`U_m = A u_m + G_m \mathbf{1} + f(\tau_m)`, -**breaking the superconvergence** and limiting: - -* velocity to :math:`\mathcal{O}(\Delta t^{M+1})` — one higher than the - stage order; -* pressure :math:`G` to :math:`\mathcal{O}(\Delta t^M)` (standard) or - increasing toward :math:`M+1` (lifted). - -This is confirmed across :math:`M = 2, 3, 4` RADAU-RIGHT and is independent -of the sweeper (SemiImplicitDAE = FullyImplicitDAE at the fixed point): - -* :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`) - -Achieving the full collocation order :math:`2M-1` for **both** velocity and -pressure would require solving all :math:`M` RADAU stages simultaneously as -a coupled system (standard RADAU-IIA), which goes beyond the node-by-node -SDC sweep used here. - -Observed results (:math:`\nu = 0.1`, ``nvars = 1023``, ``restol = 1e-13``, +1. **Standard** (no lifting, algebraic constraint): constraint + :math:`B\mathbf{u} = q(t)` has a time-dependent RHS. Causes order + reduction in :math:`G` to order :math:`M` (= 3 for 3 RADAU nodes). + +2. **Lifted** (homogeneous constraint): lifting + :math:`\mathbf{u}_\ell(t) = (q(t)/s)\,\mathbf{1}` makes the constraint + **homogeneous**: :math:`B\tilde{\mathbf{v}} = 0`. Pressure order + increases toward :math:`M+1`. + +3. **Differentiated constraint** :math:`B\mathbf{u}' = q'(t)`: replaces the + algebraic constraint with its time derivative. Converts the DAE to an + equivalent index-0 ODE system, allowing RADAU to achieve the full + collocation order :math:`2M-1` for **both** velocity and pressure. + +Key findings (: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). -* **FullyImplicitDAE** (standard or lifted): **identical results** to - SemiImplicitDAE — both converge to the same collocation fixed point. - -**Spatial resolution**: ``nvars = 1023`` interior points with a -fourth-order FD Laplacian (:math:`\Delta x = 1/1024`, spatial error floor -:math:`\mathcal{O}(\Delta x^4) \approx 10^{-12}`). +------------------------------------------------------------------------ + +* **Standard (algebraic)**: velocity :math:`\to M+1 = 4`, pressure + :math:`\to M = 3`. The :math:`\mathcal{O}(\Delta t^M)` stage pressure + errors feed back into the velocity derivatives and break the :math:`2M-1` + superconvergence. + +* **Lifted (homogeneous)**: velocity :math:`\to M+1 = 4` (unchanged); + pressure order increases monotonically toward :math:`M+1 = 4`. + +* **Differentiated constraint**: velocity :math:`\to 2M-1 = 5` ✓, + pressure :math:`\to 2M-1 = 5` ✓. By enforcing :math:`B\mathbf{U}_m = + q'(\tau_m)` at each stage instead of :math:`B\mathbf{u}_m = q(\tau_m)`, + the stage pressure errors are reduced to + :math:`e_{G_m} = -BAe_{\mathbf{u}_m}/s = \mathcal{O}(\Delta t^{M+1})`, + eliminating the feedback that breaks superconvergence. The index-reduced + system is essentially an ODE and RADAU achieves its full order :math:`2M-1`. + +Why the differentiated constraint works +---------------------------------------- +For the original algebraic constraint at each node: + +.. math:: + + B(\mathbf{v}_m + \alpha\,\mathbf{U}_m) = q(\tau_m) + \quad\Rightarrow\quad + G_m = \frac{q(\tau_m) - B\mathbf{v}_m - \alpha B\mathbf{w}} + {\alpha\,B\mathbf{v}_0} + = \mathcal{O}(\Delta t^M). + +For the differentiated constraint: + +.. math:: + + B\,\mathbf{U}_m = q'(\tau_m) + \quad\Rightarrow\quad + G_m = \frac{q'(\tau_m) - B\mathbf{w}}{B\mathbf{v}_0} + = G(\tau_m) - \frac{B A\,e_{\mathbf{u}_m}}{s} + = \mathcal{O}(\Delta t^{M+1}), + +since the stage velocity error :math:`e_{\mathbf{u}_m} = \mathcal{O}(\Delta +t^{M+1})` at collocation nodes. With :math:`G_m` at order :math:`M+1`, the +velocity endpoint recovers the full RADAU superconvergence :math:`2M-1`. Usage:: @@ -94,12 +85,11 @@ from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.projects.DAE.sweepers.semiImplicitDAE import SemiImplicitDAE -from pySDC.projects.DAE.sweepers.fullyImplicitDAE import FullyImplicitDAE from pySDC.playgrounds.Stokes_Poiseuille_1D_FD.Stokes_Poiseuille_1D_FD import ( stokes_poiseuille_1d_fd, stokes_poiseuille_1d_fd_lift, - stokes_poiseuille_1d_fd_full, - stokes_poiseuille_1d_fd_lift_full, + stokes_poiseuille_1d_fd_diffconstr, + stokes_poiseuille_1d_fd_lift_diffconstr, ) # --------------------------------------------------------------------------- @@ -130,13 +120,6 @@ # --------------------------------------------------------------------------- def _run(problem_class, sweeper_class, dt, restol=_RESTOL, max_iter=50, nvars=_NVARS): - """ - Run one simulation and return ``(uend, problem_instance)``. - - For the lifted variant ``uend.diff`` contains the *lifted* velocity - :math:`\\tilde{\\mathbf{v}}`; call :meth:`P.lift` to recover the - physical velocity. - """ desc = { 'problem_class': problem_class, 'problem_params': {'nvars': nvars, 'nu': _NU}, @@ -214,7 +197,7 @@ def _asymptotic_order(dts, errs, skip=2): def main(): r""" - Compare four formulations: two sweepers × two constraint treatments. + Compare three constraint treatments using :class:`SemiImplicitDAE`. Fixed parameters: @@ -225,43 +208,40 @@ def main(): Expected orders (see module docstring for derivation): - * **Velocity**: :math:`M+1 = 4` for all four 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). - * **SemiImplicitDAE vs FullyImplicitDAE**: identical results — both - converge to the same collocation fixed point for this DAE. + * **Standard (algebraic)**: velocity :math:`M+1 = 4`, pressure + :math:`M = 3`. + * **Lifted (homogeneous)**: velocity :math:`M+1 = 4`, pressure + approaches :math:`M+1 = 4`. + * **Differentiated constraint**: both velocity and pressure + approach full collocation order :math:`2M-1 = 5`. + * **Lifted + differentiated**: same as lifted (the differentiated + homogeneous constraint is equivalent to the original). """ - vel_order = _NUM_NODES + 1 # M+1 = 4 (U-formulation limit) - pres_order = _NUM_NODES # M = 3 (algebraic variable at each node) + colloc_order = 2 * _NUM_NODES - 1 # 2M-1 = 5 # 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'Expected velocity order M+1 = {vel_order} ' - f'(pure-ODE collocation order 2M-1 = {2*_NUM_NODES-1} not achieved; see module docstring)') - print(f'Expected pressure order M = {pres_order} (no-lift) ' - f'/ approaches M+1 = {vel_order} (lifted)') + print(f'Full RADAU collocation order: 2M-1 = {colloc_order}') 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 = [ - (stokes_poiseuille_1d_fd, SemiImplicitDAE, - 'SemiImplicitDAE, standard (B·u = q(t))'), - (stokes_poiseuille_1d_fd_lift, SemiImplicitDAE, - 'SemiImplicitDAE, lifted (B·ṽ = 0) '), - (stokes_poiseuille_1d_fd_full, FullyImplicitDAE, - 'FullyImplicitDAE, standard (B·u = q(t)) [same fixed pt as SemiImplicit]'), - (stokes_poiseuille_1d_fd_lift_full, FullyImplicitDAE, - 'FullyImplicitDAE, lifted (B·ṽ = 0) [same fixed pt as SemiImplicit]'), + (stokes_poiseuille_1d_fd, + 'Standard (B·u = q(t), algebraic constraint)'), + (stokes_poiseuille_1d_fd_lift, + 'Lifted (B·ṽ = 0, homogeneous constraint)'), + (stokes_poiseuille_1d_fd_diffconstr, + 'Diff.constr. (B·U = q\'(t), differentiated) ← remedy'), + (stokes_poiseuille_1d_fd_lift_diffconstr, + 'Lifted+diff. (B·Ũ = 0, equiv. to lifted)'), ] results = {} - for cls, sweeper, label in cases: + for cls, label in cases: print() print('=' * 72) print(f' {label}') @@ -269,12 +249,12 @@ def main(): vel_errs, pres_errs = [], [] for dt in dts: - uend, P = _run(cls, sweeper, dt) + uend, P = _run(cls, SemiImplicitDAE, dt) ve, pe = _errors(uend, P) vel_errs.append(ve) pres_errs.append(pe) - _print_table(dts, vel_errs, pres_errs, vel_order) + _print_table(dts, vel_errs, pres_errs, colloc_order) results[label] = (vel_errs, pres_errs) # ---- Summary ---- @@ -282,43 +262,40 @@ def main(): print('=' * 72) print(' Summary') print('=' * 72) - for cls, sweeper, label in cases: + for cls, label in cases: vel_errs, pres_errs = results[label] vel_ord = _asymptotic_order(dts, vel_errs) pres_ord = _asymptotic_order(dts, pres_errs) - is_lift = hasattr(cls(), 'lift') - exp_pres = vel_order if is_lift else pres_order print(f'\n {label}:') - 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} = M (time-dependent constraint)' - print(f' Pressure order ≈ {pres_ord:.1f} → {status}') - else: - print(f' Pressure order ≈ {pres_ord:.1f} → full order M+1 = {vel_order} ✓') + print(f' Velocity order ≈ {vel_ord:.1f}') + print(f' Pressure order ≈ {pres_ord:.1f}') print() print(' Conclusion:') print( - f' • All four formulations confirm velocity order M+1 = {vel_order}.' + f' • Standard (algebraic B·u=q): vel→M+1={_NUM_NODES+1}, pres→M={_NUM_NODES}.' + ) + print( + f' • Lifted (B·ṽ=0): vel→M+1, pres increasing toward M+1 (pre-asymptotic).' ) print( - f' • Standard: pressure at order M = {pres_order} (order reduction from time-dependent constraint).' + f' • Differentiated constraint (B·U=q\'): vel→2M-1={colloc_order}, ' + f'pres→2M-1={colloc_order}.' ) print( - f' • Lifted: pressure increasing toward M+1 = {vel_order} (autonomous constraint removes reduction).' + ' Enforcing B·U_m = q\'(τ_m) at each stage reduces the stage pressure\n' + ' error from O(dt^M) to O(dt^(M+1)), restoring the full RADAU\n' + ' superconvergence order 2M-1 for both velocity and pressure.\n' + ' The differentiated constraint converts the semi-explicit index-1\n' + ' DAE to an equivalent index-0 ODE system at each stage.' ) print( - ' • SemiImplicitDAE and FullyImplicitDAE converge to the SAME collocation\n' - ' fixed point for this DAE: identical velocity and pressure errors.\n' - ' The O(dt^M) stage pressure errors break the 2M-1 superconvergence\n' - ' of both sweepers; achieving 2M-1 would require solving all M RADAU\n' - ' stages simultaneously (full RADAU-IIA, beyond SDC node-by-node sweeps).' + f' • Lifted+differentiated: same as lifted (the differentiated\n' + f' homogeneous constraint B·Ũ=0 is equivalent to the original B·ṽ=0\n' + f' at the SDC fixed point; no additional improvement).' ) print( - f' • Spatial floor ~1e-12 (nvars={_NVARS}) may cause order plateau at fine Δt.' + f' • Spatial floor ~1e-12 (nvars={_NVARS}) limits the finest accessible Δt.' ) From 6e79a2c80a35cdf0737c583c535b31d29fa89b58 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 10:21:27 +0000 Subject: [PATCH 7/9] Fix diffconstr order claim: M+2 (not 2M-1); for M=3 M+2=5=2M-1 is coincidental Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/62cb79d8-8b37-46d3-89b3-20e84e919ae3 --- .../Stokes_Poiseuille_1D_FD.py | 17 +++-- .../run_convergence.py | 66 +++++++++++-------- 2 files changed, 50 insertions(+), 33 deletions(-) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py index 68a428cc3e..85e5adfa3e 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -56,8 +56,9 @@ * :class:`~pySDC.projects.DAE.sweepers.semiImplicitDAE.SemiImplicitDAE` (U-formulation): velocity order :math:`M+1`, pressure order :math:`M` - (standard), approaching :math:`M+1` (lifted), or full :math:`2M-1` - (differentiated constraint). + (standard), approaching :math:`M+1` (lifted), or order :math:`M+2` + (differentiated constraint — coincides with :math:`2M-1` only for + :math:`M = 3`). Three constraint treatments are provided — see Classes section below. @@ -108,8 +109,13 @@ G_m = \frac{q'(\tau_m) - B\mathbf{w}}{B\mathbf{v}_0}, and the stage pressure error reduces from :math:`\mathcal{O}(\Delta t^M)` -to :math:`\mathcal{O}(\Delta t^{M+1})`, restoring the full RADAU -collocation order :math:`2M-1` for **both** velocity and pressure. +to :math:`\mathcal{O}(\Delta t^{M+1})`. The U-formulation quadrature then +gives endpoint error :math:`\Delta t \cdot \mathcal{O}(\Delta t^{M+1}) = +\mathcal{O}(\Delta t^{M+2})`. **For :math:`M = 3`, :math:`M+2 = 5 = 2M-1` +coincidentally equals the full RADAU collocation order.** For :math:`M = 4`, +the order is :math:`M+2 = 6 \neq 2M-1 = 7`, as confirmed numerically. +Achieving :math:`2M-1` for :math:`M \geq 4` requires the y-formulation +(standard RADAU-IIA), not the U-formulation used here. Classes ------- @@ -123,7 +129,8 @@ stokes_poiseuille_1d_fd_diffconstr Differentiated constraint :math:`B\mathbf{U}_m = q'(\tau_m)`. - SemiImplicitDAE: **vel and pres both at** :math:`2M-1`. + SemiImplicitDAE: vel and pres both at :math:`M+2` (= :math:`2M-1` + only for :math:`M = 3`). stokes_poiseuille_1d_fd_lift_diffconstr Lifting + differentiated :math:`B\tilde{\mathbf{U}} = 0`. diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py index 268609d002..bcaebd9741 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -26,9 +26,14 @@ increases toward :math:`M+1`. 3. **Differentiated constraint** :math:`B\mathbf{u}' = q'(t)`: replaces the - algebraic constraint with its time derivative. Converts the DAE to an - equivalent index-0 ODE system, allowing RADAU to achieve the full - collocation order :math:`2M-1` for **both** velocity and pressure. + algebraic constraint with its time derivative, raising the stage pressure + accuracy by one order from :math:`\mathcal{O}(\Delta t^M)` to + :math:`\mathcal{O}(\Delta t^{M+1})`. The endpoint error in the + **U-formulation** then becomes + :math:`\Delta t\,\sum_j Q_{Mj}\,\mathcal{O}(\Delta t^{M+1}) = + \mathcal{O}(\Delta t^{M+2})`. For :math:`M = 3` this gives + :math:`M+2 = 5 = 2M-1` (coincidence!). For :math:`M \geq 4` the order + is :math:`M+2`, **not** :math:`2M-1`. Key findings (:math:`\nu = 0.1`, ``nvars = 1023``, ``restol = 1e-13``, :math:`M = 3`) @@ -36,23 +41,20 @@ * **Standard (algebraic)**: velocity :math:`\to M+1 = 4`, pressure :math:`\to M = 3`. The :math:`\mathcal{O}(\Delta t^M)` stage pressure - errors feed back into the velocity derivatives and break the :math:`2M-1` - superconvergence. + errors feed back into the velocity derivatives and break superconvergence. * **Lifted (homogeneous)**: velocity :math:`\to M+1 = 4` (unchanged); pressure order increases monotonically toward :math:`M+1 = 4`. -* **Differentiated constraint**: velocity :math:`\to 2M-1 = 5` ✓, - pressure :math:`\to 2M-1 = 5` ✓. By enforcing :math:`B\mathbf{U}_m = - q'(\tau_m)` at each stage instead of :math:`B\mathbf{u}_m = q(\tau_m)`, - the stage pressure errors are reduced to - :math:`e_{G_m} = -BAe_{\mathbf{u}_m}/s = \mathcal{O}(\Delta t^{M+1})`, - eliminating the feedback that breaks superconvergence. The index-reduced - system is essentially an ODE and RADAU achieves its full order :math:`2M-1`. +* **Differentiated constraint**: velocity :math:`\to M+2 = 5`, + pressure :math:`\to M+2 = 5`. For :math:`M = 3`, :math:`M+2 = 5 = 2M-1` + coincidentally equals the full RADAU collocation order. For :math:`M = 4`, + the observed order is :math:`M+2 = 6 \neq 2M-1 = 7`, confirming that the + improvement is one order (+1) above the standard, not a jump to :math:`2M-1`. -Why the differentiated constraint works ----------------------------------------- -For the original algebraic constraint at each node: +Why the differentiated constraint improves order by +1 +------------------------------------------------------- +For the original algebraic constraint at each SDC node: .. math:: @@ -74,7 +76,12 @@ since the stage velocity error :math:`e_{\mathbf{u}_m} = \mathcal{O}(\Delta t^{M+1})` at collocation nodes. With :math:`G_m` at order :math:`M+1`, the -velocity endpoint recovers the full RADAU superconvergence :math:`2M-1`. +stage velocity derivative errors :math:`e_{\mathbf{U}_m} = \mathcal{O}(\Delta +t^{M+1})`, and the U-formulation quadrature gives endpoint order +:math:`\Delta t \cdot \mathcal{O}(\Delta t^{M+1}) = \mathcal{O}(\Delta +t^{M+2})`. Achieving the full collocation order :math:`2M-1 \geq M+2` for +:math:`M \geq 4` would require the **y-formulation** (standard RADAU-IIA), not +the U-formulation. Usage:: @@ -212,18 +219,21 @@ def main(): :math:`M = 3`. * **Lifted (homogeneous)**: velocity :math:`M+1 = 4`, pressure approaches :math:`M+1 = 4`. - * **Differentiated constraint**: both velocity and pressure - approach full collocation order :math:`2M-1 = 5`. + * **Differentiated constraint**: both velocity and pressure at + order :math:`M+2 = 5`. For :math:`M = 3`, :math:`M+2 = 5 = 2M-1` + coincidentally equals the full RADAU collocation order; the improvement + is one order (+1) above the standard, not a jump to :math:`2M-1`. * **Lifted + differentiated**: same as lifted (the differentiated homogeneous constraint is equivalent to the original). """ - colloc_order = 2 * _NUM_NODES - 1 # 2M-1 = 5 + mplus2 = _NUM_NODES + 2 # M+2 = 5 for M=3 (coincides with 2M-1=5) # 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'Full RADAU collocation order: 2M-1 = {colloc_order}') + print(f'Standard vel order: M+1={_NUM_NODES+1}; diffconstr vel order: M+2={mplus2}') + print(f' (For M=3: M+2=5=2M-1=5 — coincidence; for M=4: M+2=6≠2M-1=7)') 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') @@ -235,7 +245,7 @@ def main(): (stokes_poiseuille_1d_fd_lift, 'Lifted (B·ṽ = 0, homogeneous constraint)'), (stokes_poiseuille_1d_fd_diffconstr, - 'Diff.constr. (B·U = q\'(t), differentiated) ← remedy'), + f'Diff.constr. (B·U = q\'(t)) order M+2={mplus2} ← remedy'), (stokes_poiseuille_1d_fd_lift_diffconstr, 'Lifted+diff. (B·Ũ = 0, equiv. to lifted)'), ] @@ -254,7 +264,7 @@ def main(): vel_errs.append(ve) pres_errs.append(pe) - _print_table(dts, vel_errs, pres_errs, colloc_order) + _print_table(dts, vel_errs, pres_errs, mplus2) results[label] = (vel_errs, pres_errs) # ---- Summary ---- @@ -279,15 +289,15 @@ def main(): f' • Lifted (B·ṽ=0): vel→M+1, pres increasing toward M+1 (pre-asymptotic).' ) print( - f' • Differentiated constraint (B·U=q\'): vel→2M-1={colloc_order}, ' - f'pres→2M-1={colloc_order}.' + f' • Differentiated constraint (B·U=q\'): vel→M+2={mplus2}, ' + f'pres→M+2={mplus2}.' ) print( ' Enforcing B·U_m = q\'(τ_m) at each stage reduces the stage pressure\n' - ' error from O(dt^M) to O(dt^(M+1)), restoring the full RADAU\n' - ' superconvergence order 2M-1 for both velocity and pressure.\n' - ' The differentiated constraint converts the semi-explicit index-1\n' - ' DAE to an equivalent index-0 ODE system at each stage.' + ' error from O(dt^M) to O(dt^{M+1}), giving endpoint order O(dt^{M+2})\n' + ' in the U-formulation. For M=3: M+2=5=2M-1 (coincidence!). For M=4:\n' + ' M+2=6, which is observed, confirming the improvement is +1 order\n' + ' over the standard, not a jump to the full collocation order 2M-1.' ) print( f' • Lifted+differentiated: same as lifted (the differentiated\n' From 0aad6709a89da7b2cdad317861c5286c562e966b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 11:45:16 +0000 Subject: [PATCH 8/9] Add stokes_poiseuille_1d_fd_coupled: (N+1)x(N+1) block solve with optional projection Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/a4d9fa91-6674-4917-9903-3f9123628b48 --- .../Stokes_Poiseuille_1D_FD.py | 269 ++++++++++++++++++ .../run_convergence.py | 41 ++- 2 files changed, 301 insertions(+), 9 deletions(-) diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py index 85e5adfa3e..ef2ae44421 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -141,6 +141,16 @@ stokes_poiseuille_1d_fd_lift_full Lifting, FullyImplicitDAE (same fixed point as lifted). + +stokes_poiseuille_1d_fd_coupled + Explicit :math:`(N+1)\times(N+1)` block solve for the differentiated + constraint (mathematically equivalent to diffconstr Schur). + Optional ``project=True`` also enforces :math:`B\mathbf{u}_m = q(\tau_m)` + via a post-solve projection, but this **degrades** convergence because + it creates an inconsistency between ``solve_system`` (which after + projection enforces the algebraic constraint) and ``eval_f`` (which + checks the differentiated constraint). Provided for pedagogical + comparison only. """ import numpy as np @@ -420,6 +430,116 @@ def _schur_solve_full_implicit(self, rhs_eff, v_approx, factor, constraint_rhs): U = w + factor * G_prime * v0 return U, float(G_prime) + def _coupled_block_solve(self, rhs_eff, v_approx, factor, q_prime_val, + project=False, q_val=None): + r""" + Fully-coupled :math:`(N+1)\times(N+1)` block solve for the + differentiated constraint :math:`B\mathbf{U} = q'(t)`. + + Assembles and solves the full block system + + .. math:: + + \begin{pmatrix} + I - \alpha\nu A & -\mathbf{1} \\ + h\,\mathbf{1}^T & 0 + \end{pmatrix} + \begin{pmatrix} \mathbf{U} \\ G \end{pmatrix} + = + \begin{pmatrix} \mathbf{r}_\text{eff} \\ q'(t) \end{pmatrix} + + This is **mathematically equivalent** to + :meth:`_schur_solve_diffconstr` (which reduces the same system to a + scalar Schur equation), but solves the full :math:`(N+1)` sparse + system directly. + + If ``project=True``, a post-solve projection step enforces the + **original** algebraic constraint :math:`B(\mathbf{v}+\alpha\mathbf{U}) + = q(t)` in addition to the differentiated one. The minimal-norm + correction + + .. math:: + + \delta = -\frac{B(\mathbf{v}+\alpha\mathbf{U}) - q(t)}{s} + \,\mathbf{1} + + is added to :math:`\mathbf{v}+\alpha\mathbf{U}`, giving + + .. math:: + + \mathbf{U}_\text{proj} = + \mathbf{U} + \frac{\delta}{\alpha}. + + The pressure :math:`G` is then updated by solving the Schur formula + :math:`G = (q'(t) - B\mathbf{U}_\text{proj}) / (B\mathbf{v}_0)` so + that :math:`B\mathbf{U}_\text{proj}` uses the corrected derivative. + + Parameters + ---------- + rhs_eff : numpy.ndarray + Effective velocity RHS + :math:`\nu A\mathbf{v}_\text{approx} + \mathbf{f}(t)`. + v_approx : numpy.ndarray + Current velocity approximation at the node. + factor : float + Implicit prefactor :math:`\alpha = \Delta t\,\tilde{q}_{mm}`. + q_prime_val : float + Value of :math:`q'(t)` at the current stage time. + project : bool, optional + If ``True``, apply the projection step that also enforces + :math:`B(\mathbf{v}+\alpha\mathbf{U}) = q(t)`. Default ``False``. + q_val : float or None + Value of :math:`q(t)` required when ``project=True``. + + Returns + ------- + U : numpy.ndarray + Velocity derivative at the node. + G_new : float + Pressure gradient satisfying :math:`B\mathbf{U} = q'(t)` (before + projection) or consistent with the projected velocity (after). + """ + n = self.nvars + top_left = sp.eye(n, format='csc') - factor * self.A + top_right = sp.csc_matrix(-self.ones.reshape(-1, 1)) + bot_left = sp.csc_matrix(self.dx * self.ones.reshape(1, -1)) + bot_right = sp.csc_matrix(np.zeros((1, 1))) + + K = sp.bmat( + [[top_left, top_right], [bot_left, bot_right]], + format='csc', + ) + rhs = np.concatenate([rhs_eff, [q_prime_val]]) + sol = spsolve(K, rhs) + + U = sol[:n].copy() + G = float(sol[n]) + + if project and q_val is not None: + # Post-solve projection: enforce B*(v + factor*U) = q(t). + # WARNING: this projection changes U_m in a way that is + # INCONSISTENT with eval_f, which checks B*u' - q'(t) = 0. + # After projection, B*U_proj = B*U - violation/factor ≠ q'(t) + # in general (since B*U = q'(t) before projection but violation ≠ 0). + # As a result, the SDC residual no longer converges cleanly and + # the sweep converges to a DIFFERENT (worse) fixed point. + # The projection is provided here for pedagogical comparison only; + # it should NOT be used in practice with the differentiated-constraint + # eval_f. + u_m = v_approx + factor * U + violation = self._B_dot(u_m) - q_val + if abs(violation) > 0.0: + delta = -(violation / self.s) * self.ones + U = U + delta / factor + # Update G from the momentum residual. + # At the unperturbed solution, (I-factor*A)*U - G*ones = rhs_eff, + # so G*ones = (I-factor*A)*U - rhs_eff. After projection, U changes + # and all N components of (I-factor*A)*U_proj - rhs_eff theoretically + # equal the same G value; we take the mean for numerical stability. + G = float(np.mean(top_left.dot(U) - rhs_eff)) + + return U, G + def _schur_solve_diffconstr(self, rhs_eff, factor, q_prime_val): r""" Schur-complement solve using the **differentiated constraint** @@ -1193,3 +1313,152 @@ def solve_system(self, impl_sys, u_approx, factor, u0, t): me.diff[:] = U me.alg[0] = G_new return me + + +# --------------------------------------------------------------------------- +# Case 7: Coupled block solve + projection +# --------------------------------------------------------------------------- + +class stokes_poiseuille_1d_fd_coupled(stokes_poiseuille_1d_fd): + r""" + 1-D Stokes/Poiseuille DAE using an **explicit** :math:`(N+1)\times(N+1)` + block solve for the differentiated constraint, with an optional + post-solve projection that also enforces the original algebraic constraint. + + Two sub-variants are provided via the ``project`` constructor parameter: + + **Block solve only** (``project=False``, default) + Assembles and solves the full :math:`(N+1)\times(N+1)` sparse system + + .. math:: + + \begin{pmatrix} + I - \alpha\nu A & -\mathbf{1} \\ + h\,\mathbf{1}^T & 0 + \end{pmatrix} + \begin{pmatrix} \mathbf{U} \\ G \end{pmatrix} + = + \begin{pmatrix} + \nu A\mathbf{v} + \mathbf{f}(\tau_m) \\ + q'(\tau_m) + \end{pmatrix} + + This is **mathematically equivalent** to + :class:`stokes_poiseuille_1d_fd_diffconstr` (which reduces the same + system to a scalar Schur equation). Convergence orders are the same: + velocity :math:`M+2`, pressure :math:`M+2` (:math:`= 2M-1` for + :math:`M = 3` by coincidence). + + **Block solve + projection** (``project=True``) + After the block solve, a minimal-norm correction enforces the + **original** algebraic constraint :math:`B(\mathbf{v}+\alpha\mathbf{U}) + = q(\tau_m)` as well. + + .. warning:: + + This variant gives **worse** results than the plain block solve. + The root cause is an inconsistency between ``solve_system`` (which + after projection enforces :math:`B(\mathbf{v}+\alpha\mathbf{U}) + = q(t)`) and ``eval_f`` (which checks the **differentiated** + constraint :math:`B\mathbf{u}' - q'(t) = 0`). Because the + projection changes :math:`\mathbf{U}` so that + :math:`B\mathbf{U} \neq q'(t)` any more, the SDC residual + never converges cleanly and the sweep converges to a different, + lower-accuracy fixed point. Numerically, the projection variant + achieves only velocity :math:`M+1 \approx 4`, pressure + :math:`M \approx 3` — the same as the standard algebraic formulation. + + The lesson is that **self-consistency between** ``solve_system`` **and** + ``eval_f`` **is essential**: both must enforce the same constraint + (either the algebraic :math:`B\mathbf{u}=q` or the differentiated + :math:`B\mathbf{u}'=q'`). Mixing the two degrades convergence. + + The ``eval_f`` uses the differentiated constraint + :math:`F_\text{alg} = B\mathbf{u}' - q'(t) = 0`, matching + :class:`stokes_poiseuille_1d_fd_diffconstr`. + + Parameters + ---------- + nvars : int + Number of interior grid points (default 127; must be ≥ 5). + nu : float + Kinematic viscosity (default 1.0). + newton_tol : float + Unused; passed to base class (default 1e-10). + project : bool + If ``True``, apply the post-solve projection step that also enforces + :math:`B(\mathbf{v}+\alpha\mathbf{U}) = q(\tau_m)`. Default ``False``. + See warning above — this is provided for pedagogical comparison only. + """ + + def __init__(self, nvars=127, nu=1.0, newton_tol=1e-10, project=False): + super().__init__(nvars=nvars, nu=nu, newton_tol=newton_tol) + self._makeAttributeAndRegister('project', localVars=locals(), readOnly=True) + + def eval_f(self, u, du, t): + r""" + Fully-implicit DAE residual using the **differentiated constraint**: + + .. math:: + + F_\text{diff} = \mathbf{u}' - \nu A\,\mathbf{u} + - G\,\mathbf{1} - \mathbf{f}(t), + + .. math:: + + F_\text{alg} = B\,\mathbf{u}' - q'(t). + + Identical to :class:`stokes_poiseuille_1d_fd_diffconstr`. + """ + f = self.dtype_f(self.init, val=0.0) + u_vel = np.asarray(u.diff) + du_vel = np.asarray(du.diff) + G = float(u.alg[0]) + + f.diff[:] = du_vel - (self.A.dot(u_vel) + G * self.ones + self._forcing(t)) + f.alg[0] = self._B_dot(du_vel) - self._q_prime(t) + + self.work_counters['rhs']() + return f + + def solve_system(self, impl_sys, u_approx, factor, u0, t): + r""" + Coupled :math:`(N+1)\times(N+1)` block solve with the differentiated + constraint :math:`B\mathbf{U} = q'(t)`. + + Optionally applies a post-solve projection onto + :math:`B(\mathbf{v}+\alpha\mathbf{U}) = q(t)` if ``self.project`` + is ``True``. + + Parameters + ---------- + impl_sys : callable + Unused; system solved directly. + u_approx : MeshDAE + Current velocity approximation at the node. + factor : float + Implicit prefactor :math:`\alpha`. + u0 : MeshDAE + Unused (direct solver). + t : float + Current time. + + Returns + ------- + me : MeshDAE + ``me.diff[:]`` = velocity derivative :math:`\mathbf{U}_m`, + ``me.alg[0]`` = pressure gradient :math:`G_m`. + """ + me = self.dtype_u(self.init, val=0.0) + v_approx = np.asarray(u_approx.diff).copy() + + rhs_eff = self.A.dot(v_approx) + self._forcing(t) + U, G_new = self._coupled_block_solve( + rhs_eff, v_approx, factor, self._q_prime(t), + project=self.project, + q_val=self._q(t) if self.project else None, + ) + + me.diff[:] = U + me.alg[0] = G_new + return me diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py index bcaebd9741..b65f7f6aa3 100644 --- a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -97,6 +97,7 @@ stokes_poiseuille_1d_fd_lift, stokes_poiseuille_1d_fd_diffconstr, stokes_poiseuille_1d_fd_lift_diffconstr, + stokes_poiseuille_1d_fd_coupled, ) # --------------------------------------------------------------------------- @@ -126,10 +127,14 @@ # Helpers # --------------------------------------------------------------------------- -def _run(problem_class, sweeper_class, dt, restol=_RESTOL, max_iter=50, nvars=_NVARS): +def _run(problem_class, sweeper_class, dt, restol=_RESTOL, max_iter=50, nvars=_NVARS, + extra_problem_params=None): + prob_params = {'nvars': nvars, 'nu': _NU} + if extra_problem_params: + prob_params.update(extra_problem_params) desc = { 'problem_class': problem_class, - 'problem_params': {'nvars': nvars, 'nu': _NU}, + 'problem_params': prob_params, 'sweeper_class': sweeper_class, 'sweeper_params': _SWEEPER_PARAMS, 'level_params': {'restol': restol, 'dt': dt}, @@ -241,17 +246,27 @@ def main(): cases = [ (stokes_poiseuille_1d_fd, - 'Standard (B·u = q(t), algebraic constraint)'), + 'Standard (B·u = q(t), algebraic constraint)', + {}), (stokes_poiseuille_1d_fd_lift, - 'Lifted (B·ṽ = 0, homogeneous constraint)'), + 'Lifted (B·ṽ = 0, homogeneous constraint)', + {}), (stokes_poiseuille_1d_fd_diffconstr, - f'Diff.constr. (B·U = q\'(t)) order M+2={mplus2} ← remedy'), + f'Diff.constr. (B·U = q\'(t)) order M+2={mplus2} ← remedy', + {}), (stokes_poiseuille_1d_fd_lift_diffconstr, - 'Lifted+diff. (B·Ũ = 0, equiv. to lifted)'), + 'Lifted+diff. (B·Ũ = 0, equiv. to lifted)', + {}), + (stokes_poiseuille_1d_fd_coupled, + f'Coupled block (N+1)×(N+1) solve, no projection [≡ diffconstr]', + {'project': False}), + (stokes_poiseuille_1d_fd_coupled, + f'Coupled block + projection B·u_m=q(τ_m) [degrades!]', + {'project': True}), ] results = {} - for cls, label in cases: + for cls, label, extra_params in cases: print() print('=' * 72) print(f' {label}') @@ -259,7 +274,7 @@ def main(): vel_errs, pres_errs = [], [] for dt in dts: - uend, P = _run(cls, SemiImplicitDAE, dt) + uend, P = _run(cls, SemiImplicitDAE, dt, extra_problem_params=extra_params) ve, pe = _errors(uend, P) vel_errs.append(ve) pres_errs.append(pe) @@ -272,7 +287,7 @@ def main(): print('=' * 72) print(' Summary') print('=' * 72) - for cls, label in cases: + for cls, label, extra_params in cases: vel_errs, pres_errs = results[label] vel_ord = _asymptotic_order(dts, vel_errs) pres_ord = _asymptotic_order(dts, pres_errs) @@ -304,6 +319,14 @@ def main(): f' homogeneous constraint B·Ũ=0 is equivalent to the original B·ṽ=0\n' f' at the SDC fixed point; no additional improvement).' ) + print( + f' • Coupled block (N+1)×(N+1) solve (no proj): mathematically equivalent\n' + f' to diffconstr (same Schur system, different implementation). Same orders.\n' + f' • Coupled + projection B·u_m=q(τ_m): DEGRADES convergence! The projection\n' + f' changes U so that B·U ≠ q\'(t) any more, creating an inconsistency with\n' + f' eval_f (which checks B·u\'-q\'=0). Lesson: solve_system and eval_f must\n' + f' enforce the SAME constraint — mixing algebraic and differentiated degrades.' + ) print( f' • Spatial floor ~1e-12 (nvars={_NVARS}) limits the finest accessible Δt.' ) From 2016d9028907f76ccf83a415a58a6a1b84219660 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 26 Mar 2026 15:05:18 +0000 Subject: [PATCH 9/9] Add run_monolithic.py: monolithic RADAU-IIA collocation (algebraic and differentiated-constraint variants) Co-authored-by: pancetta <7158893+pancetta@users.noreply.github.com> Agent-Logs-Url: https://github.com/Parallel-in-Time/pySDC/sessions/bc544fa4-6fe9-4053-95f8-9c6e2c1e884f --- .../Stokes_Poiseuille_1D_FD/run_monolithic.py | 414 ++++++++++++++++++ 1 file changed, 414 insertions(+) create mode 100644 pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_monolithic.py diff --git a/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_monolithic.py b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_monolithic.py new file mode 100644 index 0000000000..0a37eb896f --- /dev/null +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_monolithic.py @@ -0,0 +1,414 @@ +r""" +Monolithic RADAU-IIA collocation solver for the 1-D Stokes/Poiseuille DAE +========================================================================== + +This script implements a **monolithic all-at-once** RADAU-IIA collocation +solver for the 1-D Stokes/Poiseuille index-1 DAE: + +.. math:: + + \mathbf{u}' = \nu A\,\mathbf{u} + G(t)\,\mathbf{1} + \mathbf{f}(t), + \quad + 0 = B\,\mathbf{u} - q(t). + +In contrast to the node-by-node SDC sweep in :mod:`run_convergence`, the +monolithic solver assembles all :math:`M` RADAU stages simultaneously into a +single :math:`M(N+1) \times M(N+1)` sparse saddle-point system and solves it +in one shot. The endpoint is taken as the last stage value +:math:`\mathbf{u}_{n+1} = \mathbf{u}_M` (the **y-formulation**, since +RADAU-RIGHT has :math:`c_M = 1`). + +Key finding: monolithic vs.\ SDC standard are **identical** +------------------------------------------------------------ +For the **original algebraic constraint** :math:`B\mathbf{u}_m = q(\tau_m)`, +the monolithic RADAU-IIA collocation system is **mathematically identical** to +the SDC collocation fixed point. Both solvers compute the same stage values +:math:`\mathbf{u}_m`, and both produce the same endpoint errors. The observed +orders are :math:`M+1` for velocity and :math:`M` for pressure — the **same +order reduction** seen in the SDC standard sweep. + +This is not a coincidence. For a semi-explicit index-1 DAE with a time- +dependent algebraic constraint, the theory (Hairer & Wanner) predicts that +RADAU-IIA with the algebraic constraint at each stage achieves only the +**stage order** :math:`M+1` at the endpoint (not the full :math:`2M-1` +order), because the stage accuracy is limited by the constraint. Switching +from iterative SDC sweeping to a single direct solve does not change this. + +To recover full :math:`2M-1` order it is necessary to replace the algebraic +constraint :math:`B\mathbf{u}_m = q(\tau_m)` by its **differentiated form** +:math:`B\mathbf{u}_m' = q'(\tau_m)` (index reduction). This is implemented +in :func:`_monolithic_step_diffconstr` and achieves :math:`2M-1` order for +both velocity and pressure. + +Monolithic collocation system — two constraint variants +-------------------------------------------------------- +For each stage :math:`m = 1,\ldots,M`: + +.. math:: + + \mathbf{u}_m + = \mathbf{u}_0 + \Delta t\sum_{j=1}^M Q_{mj} + \bigl(\nu A\,\mathbf{u}_j + G_j\,\mathbf{1} + \mathbf{f}(\tau_j) + \bigr). + +**Algebraic constraint** (order :math:`M+1` / :math:`M`): + +.. math:: + + B\,\mathbf{u}_m = q(\tau_m). + +Block system :math:`K x = b`: + +.. math:: + + \underbrace{\begin{pmatrix} + I_{MN} - \Delta t\,(Q \otimes \nu A) & -\Delta t\,(Q \otimes \mathbf{1}) \\ + I_M \otimes B & 0_{M \times M} + \end{pmatrix}}_{K_\text{alg}} + \begin{pmatrix} \mathbf{u}_\text{vec} \\ G_\text{vec} \end{pmatrix} + = + \begin{pmatrix} + \mathbf{1}_M \otimes \mathbf{u}_0 + + \Delta t\,(Q \otimes I_N)\,\mathbf{f}_\text{vec} \\ + \mathbf{q}_\text{vec} + \end{pmatrix}. + +**Differentiated constraint** (order :math:`2M-1` / :math:`2M-1`): + +.. math:: + + B\,\mathbf{u}_m' = q'(\tau_m), + \quad + \text{i.e.,}\quad + (B\nu A)\,\mathbf{u}_m + s\,G_m = q'(\tau_m) - B\,\mathbf{f}(\tau_m). + +Same velocity block :math:`(I_{MN} - \Delta t(Q \otimes \nu A),\; +-\Delta t(Q \otimes \mathbf{1}))`, but the constraint rows change to: + +.. math:: + + \underbrace{\begin{pmatrix} + I_{MN} - \Delta t\,(Q \otimes \nu A) & -\Delta t\,(Q \otimes \mathbf{1}) \\ + I_M \otimes (B\nu A) & s\,I_M + \end{pmatrix}}_{K_\text{dc}} + \begin{pmatrix} \mathbf{u}_\text{vec} \\ G_\text{vec} \end{pmatrix} + = + \begin{pmatrix} + \mathbf{1}_M \otimes \mathbf{u}_0 + + \Delta t\,(Q \otimes I_N)\,\mathbf{f}_\text{vec} \\ + [q'(\tau_m) - B\mathbf{f}(\tau_m)]_{m=1}^M + \end{pmatrix}. + +For the differentiated-constraint system the algebraic variable is determined +by an explicit formula once :math:`\mathbf{u}_m` is known, and RADAU-IIA can +achieve its full stage accuracy :math:`2M-1`. + +Summary of orders +----------------- ++-----------------------------------+----------------+----------------+ +| Method | Velocity order | Pressure order | ++===================================+================+================+ +| SDC standard (U-form., alg.) | :math:`M+1` | :math:`M` | ++-----------------------------------+----------------+----------------+ +| SDC diffconstr (U-form., dc) | :math:`M+2` | :math:`M+2` | ++-----------------------------------+----------------+----------------+ +| Monolithic alg. ≡ SDC std fixed pt| :math:`M+1` | :math:`M` | ++-----------------------------------+----------------+----------------+ +| **Monolithic diffconstr (y-form)**| :math:`2M-1` | :math:`2M-1` | ++-----------------------------------+----------------+----------------+ + +For :math:`M = 3`: :math:`M+2 = 5 = 2M-1`, so SDC diffconstr and monolithic +diffconstr give the same asymptotic order (with different constants). +For :math:`M = 4`: monolithic diffconstr gives :math:`7 = 2M-1`, while SDC +diffconstr gives only :math:`6 = M+2`. + +Usage:: + + python run_monolithic.py +""" + +import numpy as np +import scipy.sparse as sp +from scipy.sparse.linalg import spsolve + +from pySDC.core.collocation import CollBase +from pySDC.playgrounds.Stokes_Poiseuille_1D_FD.Stokes_Poiseuille_1D_FD import ( + stokes_poiseuille_1d_fd, +) + +# --------------------------------------------------------------------------- +# Parameters +# --------------------------------------------------------------------------- + +_NVARS = 1023 +_NU = 0.1 +_TEND = 1.0 + + +# --------------------------------------------------------------------------- +# Helpers: build common blocks +# --------------------------------------------------------------------------- + +def _common_velocity_blocks(P, u0, t0, dt, Q, nodes): + """Build the velocity RHS and the velocity-velocity / velocity-pressure blocks.""" + N = P.nvars + M = len(nodes) + A = P.A + + tau = t0 + dt * nodes + f_stages = np.array([P._forcing(tau[m]) for m in range(M)]) + forcing_contrib = dt * (Q @ f_stages) + rhs_vel = np.tile(u0, M).reshape(M, N) + forcing_contrib + rhs_vel_flat = rhs_vel.flatten() + + KUU = ( + sp.eye(M * N, format='csc') + - dt * sp.kron(sp.csc_matrix(Q), A, format='csc') + ) + ones_col = np.ones((N, 1)) + KUG = -dt * sp.kron( + sp.csc_matrix(Q), + sp.csc_matrix(ones_col), + format='csc', + ) + return tau, f_stages, rhs_vel_flat, KUU, KUG + + +# --------------------------------------------------------------------------- +# Monolithic step — algebraic constraint B·u_m = q(τ_m) +# --------------------------------------------------------------------------- + +def _monolithic_step_algebraic(P, u0, t0, dt, Q, nodes): + r""" + Single RADAU-IIA step with the **algebraic constraint** + :math:`B\mathbf{u}_m = q(\tau_m)`. + + Produces the same collocation fixed point as the SDC standard sweep. + Returns ``(u_{n+1}, G_{n+1})`` where :math:`u_{n+1} = \mathbf{u}_M` + (last stage, y-formulation). + """ + N = P.nvars + M = len(nodes) + dx = P.dx + + tau, _f, rhs_vel_flat, KUU, KUG = _common_velocity_blocks( + P, u0, t0, dt, Q, nodes + ) + + rhs_constr = np.array([P._q(tau[m]) for m in range(M)]) + rhs = np.concatenate([rhs_vel_flat, rhs_constr]) + + B_row = dx * np.ones((1, N)) + KCU = sp.kron(sp.eye(M, format='csc'), sp.csc_matrix(B_row), format='csc') + KCG = sp.csc_matrix((M, M)) + + K = sp.bmat([[KUU, KUG], [KCU, KCG]], format='csc') + x = spsolve(K, rhs) + + u_end = x[(M - 1) * N: M * N].copy() + G_end = float(x[M * N + (M - 1)]) + return u_end, G_end + + +# --------------------------------------------------------------------------- +# Monolithic step — differentiated constraint B·u_m' = q'(τ_m) +# --------------------------------------------------------------------------- + +def _monolithic_step_diffconstr(P, u0, t0, dt, Q, nodes): + r""" + Single RADAU-IIA step with the **differentiated constraint** + :math:`B\mathbf{u}'_m = q'(\tau_m)`. + + Substituting :math:`\mathbf{u}'_m = \nu A\mathbf{u}_m + G_m\mathbf{1} + + \mathbf{f}(\tau_m)`, the constraint becomes + + .. math:: + + (B\nu A)\,\mathbf{u}_m + s\,G_m = q'(\tau_m) - B\mathbf{f}(\tau_m), + + where :math:`s = B\mathbf{1}`. The modified constraint block makes + :math:`K_\text{dc}` non-singular and the y-formulation endpoint + :math:`\mathbf{u}_M` achieves the full :math:`2M-1` collocation order. + """ + N = P.nvars + M = len(nodes) + A = P.A + s = P.s # B·1 = dx·N + + tau, f_stages, rhs_vel_flat, KUU, KUG = _common_velocity_blocks( + P, u0, t0, dt, Q, nodes + ) + + # Differentiated constraint RHS: q'(τ_m) − B·f(τ_m) + rhs_constr = np.array([ + P._q_prime(tau[m]) - P.dx * float(np.sum(f_stages[m])) + for m in range(M) + ]) + rhs = np.concatenate([rhs_vel_flat, rhs_constr]) + + # KCU: I_M ⊗ (B·νA) — (B·νA) is a 1×N row = dx·1ᵀ·A + BA_vec = P.dx * (P.ones @ A) # (N,) dense, B·νA row vector + BA = sp.csc_matrix(BA_vec.reshape(1, N)) # (1, N) sparse + KCU = sp.kron( + sp.eye(M, format='csc'), + sp.csc_matrix(BA), + format='csc', + ) + + # KCG: s·I_M (not zero — algebraic constraint involves G_m) + KCG = s * sp.eye(M, format='csc') + + K = sp.bmat([[KUU, KUG], [KCU, KCG]], format='csc') + x = spsolve(K, rhs) + + u_end = x[(M - 1) * N: M * N].copy() + G_end = float(x[M * N + (M - 1)]) + return u_end, G_end + + +# --------------------------------------------------------------------------- +# Time-stepping loop +# --------------------------------------------------------------------------- + +def _run_monolithic(step_fn, M, dt, nvars=_NVARS, nu=_NU, Tend=_TEND): + """Integrate using *step_fn* and return ``(vel_err, pres_err)``.""" + coll = CollBase( + num_nodes=M, tleft=0, tright=1, + node_type='LEGENDRE', quad_type='RADAU-RIGHT', + ) + Q = coll.Qmat[1:, 1:] + nodes = coll.nodes + + P = stokes_poiseuille_1d_fd(nvars=nvars, nu=nu) + t = 0.0 + u = np.sin(np.pi * P.xvalues) * np.sin(t) + + n_steps = int(round(Tend / dt)) + for _ in range(n_steps): + u, G = step_fn(P, u, t, dt, Q, nodes) + t += dt + + u_ex = np.sin(np.pi * P.xvalues) * np.sin(Tend) + G_ex = np.cos(Tend) + return float(np.max(np.abs(u - u_ex))), float(abs(G - G_ex)) + + +# --------------------------------------------------------------------------- +# Convergence table helper +# --------------------------------------------------------------------------- + +def _print_table(dts, vel_errs, pres_errs, exp_ord): + print( + f' {"dt":>10} {"vel error":>14} {"vel ord":>8} {"exp":>4}' + f' {"pres error":>14} {"pres ord":>9}' + ) + for i, dt in enumerate(dts): + ve, pe = vel_errs[i], pres_errs[i] + vo_str = ( + f'{np.log(vel_errs[i-1]/ve)/np.log(dts[i-1]/dt):>8.2f}' + if i > 0 and vel_errs[i-1] > 0 and ve > 0 else f'{"---":>8}' + ) + po_str = ( + f'{np.log(pres_errs[i-1]/pe)/np.log(dts[i-1]/dt):>9.2f}' + if i > 0 and pres_errs[i-1] > 0 and pe > 0 else f'{"---":>9}' + ) + print( + f' {dt:>10.5f} {ve:>14.4e} {vo_str} {exp_ord:>4d}' + f' {pe:>14.4e} {po_str}' + ) + + +# --------------------------------------------------------------------------- +# Main convergence study +# --------------------------------------------------------------------------- + +def main(): + r""" + Convergence study comparing algebraic and differentiated-constraint + monolithic RADAU-IIA solvers. + + Key findings: + + 1. **Algebraic constraint** — monolithic :math:`\equiv` SDC standard: + The monolithic y-formulation and the fully-converged SDC standard sweep + produce **identical** errors (same collocation fixed point). Both + achieve :math:`M+1` for velocity and :math:`M` for pressure. + **The order reduction is inherent in the DAE, not due to SDC sweeping.** + + 2. **Differentiated constraint** — monolithic achieves :math:`2M-1`: + Replacing :math:`B\mathbf{u}_m = q(\tau_m)` with the differentiated + form :math:`B\mathbf{u}'_m = q'(\tau_m)` converts the constraint block + and allows RADAU-IIA to reach its full :math:`2M-1` collocation order + for **both** velocity and pressure. + + For :math:`M = 3`, both SDC diffconstr (:math:`M+2 = 5`) and monolithic + diffconstr (:math:`2M-1 = 5`) give the same asymptotic order (the + coincidence :math:`M+2 = 2M-1` for :math:`M = 3`). For :math:`M = 4`, + monolithic diffconstr achieves :math:`2M-1 = 7`, while SDC diffconstr + gives only :math:`M+2 = 6`. + """ + print('Monolithic RADAU-IIA collocation (y-formulation)') + print(f'ν = {_NU}, nvars = {_NVARS}, T_end = {_TEND}\n') + + dts = [_TEND / (2 ** k) for k in range(1, 8)] + + for M in [3, 4]: + two_m_minus1 = 2 * M - 1 + + # ── Algebraic ──────────────────────────────────────────────────────── + print('=' * 72) + print(f' M = {M} Algebraic constraint B·u_m = q(τ_m)') + print(f' Expected order: M+1 = {M+1} (vel), M = {M} (pres)') + print(f' [Same as SDC standard — order reduction inherent in DAE]') + print('=' * 72) + + vel_errs, pres_errs = [], [] + for dt in dts: + ve, pe = _run_monolithic(_monolithic_step_algebraic, M, dt) + vel_errs.append(ve) + pres_errs.append(pe) + _print_table(dts, vel_errs, pres_errs, M + 1) + print() + + # ── Differentiated ─────────────────────────────────────────────────── + print('=' * 72) + print(f' M = {M} Differentiated constraint B·u_m\' = q\'(τ_m)') + print(f' Expected order: 2M-1 = {two_m_minus1} (vel and pres)') + print(f' [Full RADAU-IIA collocation order via index reduction]') + print('=' * 72) + + vel_errs, pres_errs = [], [] + for dt in dts: + ve, pe = _run_monolithic(_monolithic_step_diffconstr, M, dt) + vel_errs.append(ve) + pres_errs.append(pe) + _print_table(dts, vel_errs, pres_errs, two_m_minus1) + print() + + print('=' * 72) + print(' Conclusion') + print('=' * 72) + print( + '\n Algebraic constraint B·u_m = q(τ_m):' + '\n Monolithic y-formulation produces IDENTICAL results to SDC standard.' + '\n Both achieve vel → M+1, pres → M.' + '\n Order reduction is INHERENT in the DAE structure, not due to SDC.' + '\n' + '\n Differentiated constraint B·u_m\' = q\'(τ_m):' + '\n Monolithic achieves full 2M-1 collocation order for BOTH vel and pres.' + '\n This beats SDC diffconstr (U-form.) for M ≥ 4: 2M-1 vs. M+2.' + '\n For M=3: 2M-1=5=M+2 (coincidence); for M=4: 2M-1=7 > M+2=6.' + '\n' + f'\n Method comparison (M=3, M=4): vel | pres' + f'\n SDC standard (U-form., alg. constr.): M+1 | M' + f'\n SDC diffconstr (U-form., diff. constr.): M+2 | M+2' + f'\n Monolithic alg. (y-form.) ≡ SDC standard: M+1 | M' + f'\n Monolithic diffconstr (y-form.): 2M-1 | 2M-1 ← best' + '\n' + '\n Spatial floor ~1e-12 (nvars=1023) may limit the finest accessible Δt,' + '\n especially for the higher-order methods that converge more rapidly.' + ) + + +if __name__ == '__main__': + main()