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..ef2ae44421 --- /dev/null +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/Stokes_Poiseuille_1D_FD.py @@ -0,0 +1,1464 @@ +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 +flow-rate 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 pressure gradient and :math:`q(t)` is a +prescribed flow-rate. After finite-difference discretisation this becomes +the semi-explicit index-1 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 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 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). + +Forcing: + +.. math:: + + f(y, t) = \sin(\pi y)\cos(t) + \nu\pi^2\sin(\pi y)\sin(t) - \cos(t). + +State and sweeper +----------------- +The state variable uses :class:`~pySDC.projects.DAE.misc.meshDAE.MeshDAE` +with ``nvars`` interior points: + +* ``u.diff[:]`` – velocity on :math:`N` interior grid points. +* ``u.alg[0]`` – pressure gradient :math:`G` (Lagrange multiplier). + +Two sweepers are supported: + +* :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 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. + +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:: + + \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:: + + 0 = B\,\tilde{\mathbf{v}}, + +and evolves according to + +.. math:: + + \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}. + +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})`. 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 +------- +stokes_poiseuille_1d_fd + 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`. + 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:`M+2` (= :math:`2M-1` + only for :math:`M = 3`). + +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). + +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 +import scipy.sparse as sp +from scipy.sparse.linalg import spsolve + +from pySDC.projects.DAE.misc.problemDAE import ProblemDAE + + +# --------------------------------------------------------------------------- +# 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. With homogeneous + BCs the boundary-correction vector :math:`b_\text{bc}` is identically + zero. + + 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). + 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} = 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 = 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). + 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() + + +# --------------------------------------------------------------------------- +# Shared base class +# --------------------------------------------------------------------------- + +class _StokesBase(ProblemDAE): + r""" + Shared setup for the Stokes/Poiseuille problem classes. + + Builds the 4th-order FD Laplacian, grid coordinates, and the + manufactured-forcing helpers. + + Parameters + ---------- + nvars : int + Number of interior grid points (must be ≥ 5). + nu : float + Kinematic viscosity. + newton_tol : float + Tolerance passed to :class:`~pySDC.projects.DAE.misc.problemDAE.ProblemDAE`. + + 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, shape ``(nvars,)``. + C_B : float + :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, nu, newton_tol): + if nvars < 5: + raise ValueError( + f'nvars must be >= 5 for the 4th-order FD Laplacian; got {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) + 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))) + 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)`. + """ + 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 + :math:`G_\text{ex} = \cos(t)`: + + .. math:: + + f(y, t) = \sin(\pi y)\cos(t) + + \nu\pi^2\sin(\pi y)\sin(t) - \cos(t). + """ + 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 + ) + + 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 _schur_solve(self, rhs_eff, v_approx, factor, constraint_rhs): + r""" + Schur-complement saddle-point solve. + + Finds :math:`(\mathbf{U}, G)` satisfying + + .. math:: + + (I - \alpha\nu A)\,\mathbf{U} - G\,\mathbf{1} = \mathbf{r}_\text{eff}, + + .. math:: + + 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 + ---------- + 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 + ------- + U : numpy.ndarray + Velocity derivative at the node. + G_new : float + Pressure gradient (Lagrange multiplier). + """ + 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 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) + + 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) + + 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** + :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) +# --------------------------------------------------------------------------- + +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:: + + F_\text{diff} = \mathbf{u}' - \nu A\,\mathbf{u} + - G\,\mathbf{1} - \mathbf{f}(t), + + .. math:: + + 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]) + + 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) + + self.work_counters['rhs']() + return f + + 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 + 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._schur_solve(rhs_eff, v_approx, factor, self._q(t)) + + me.diff[:] = U + me.alg[0] = G_new + return me + + 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 + + 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 + + +# --------------------------------------------------------------------------- +# Case 2: Constraint lifting – homogeneous constraint B·ṽ = 0 +# --------------------------------------------------------------------------- + +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""" + Lifting function :math:`\mathbf{u}_\ell(t) = (q(t)/s)\,\mathbf{1}`. + + Satisfies :math:`B\mathbf{u}_\ell(t) = q(t)` exactly. + + Parameters + ---------- + t : float + + Returns + ------- + numpy.ndarray, shape (nvars,) + """ + return (self.C_B * np.sin(t) / self.s) * self.ones + + 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""" + 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[:]`` = 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) + 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 + + +# --------------------------------------------------------------------------- +# 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 + + +# --------------------------------------------------------------------------- +# 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 + + +# --------------------------------------------------------------------------- +# 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/__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..b65f7f6aa3 --- /dev/null +++ b/pySDC/playgrounds/Stokes_Poiseuille_1D_FD/run_convergence.py @@ -0,0 +1,337 @@ +r""" +Temporal order of convergence for the 1-D Stokes/Poiseuille index-1 DAE +======================================================================== + +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: + +.. math:: + + \mathbf{u}' = \nu A\,\mathbf{u} + G(t)\,\mathbf{1} + \mathbf{f}(t), + +.. math:: + + 0 = B\,\mathbf{u} - q(t). + +Three constraint treatments are compared +----------------------------------------- +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, 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`) +------------------------------------------------------------------------ + +* **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 superconvergence. + +* **Lifted (homogeneous)**: velocity :math:`\to M+1 = 4` (unchanged); + pressure order increases monotonically toward :math:`M+1 = 4`. + +* **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 improves order by +1 +------------------------------------------------------- +For the original algebraic constraint at each SDC 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 +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:: + + python run_convergence.py +""" + +import numpy as np + +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, + stokes_poiseuille_1d_fd_lift, + stokes_poiseuille_1d_fd_diffconstr, + stokes_poiseuille_1d_fd_lift_diffconstr, + stokes_poiseuille_1d_fd_coupled, +) + +# --------------------------------------------------------------------------- +# Fixed parameters +# --------------------------------------------------------------------------- + +_NVARS = 1023 +# 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 + +_SWEEPER_PARAMS = { + 'quad_type': 'RADAU-RIGHT', + 'num_nodes': _NUM_NODES, + 'QI': 'LU', + 'initial_guess': 'spread', +} + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +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': prob_params, + 'sweeper_class': sweeper_class, + '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 uend, P + + +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 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) + 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}' + f' {"pres error":>14} {"pres ord":>9}' + ) + 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}' + ) + + +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""" + Compare three constraint treatments using :class:`SemiImplicitDAE`. + + Fixed parameters: + + * ``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} = 1.0`. + + Expected orders (see module docstring for derivation): + + * **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 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). + """ + 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'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') + print(f'Error vs. exact analytical solution at T={_TEND}') + + cases = [ + (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, + 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)', + {}), + (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, extra_params in cases: + print() + print('=' * 72) + print(f' {label}') + print('=' * 72) + + vel_errs, pres_errs = [], [] + for dt in dts: + uend, P = _run(cls, SemiImplicitDAE, dt, extra_problem_params=extra_params) + ve, pe = _errors(uend, P) + vel_errs.append(ve) + pres_errs.append(pe) + + _print_table(dts, vel_errs, pres_errs, mplus2) + results[label] = (vel_errs, pres_errs) + + # ---- Summary ---- + print() + print('=' * 72) + print(' Summary') + print('=' * 72) + 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) + print(f'\n {label}:') + print(f' Velocity order ≈ {vel_ord:.1f}') + print(f' Pressure order ≈ {pres_ord:.1f}') + + print() + print(' Conclusion:') + print( + 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' • 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}), 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' + 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.' + ) + + +if __name__ == '__main__': + main() + 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()