diff --git a/pySDC/playgrounds/FEniCS/allen_cahn_1d/__init__.py b/pySDC/playgrounds/FEniCS/allen_cahn_1d/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/pySDC/playgrounds/FEniCS/allen_cahn_1d/problem_classes.py b/pySDC/playgrounds/FEniCS/allen_cahn_1d/problem_classes.py new file mode 100644 index 0000000000..d1e607ef60 --- /dev/null +++ b/pySDC/playgrounds/FEniCS/allen_cahn_1d/problem_classes.py @@ -0,0 +1,486 @@ +""" +Problem classes for the 1D Allen-Cahn FEniCS playground. + +Demonstrates SDC order reduction caused by inhomogeneous time-dependent +Dirichlet boundary conditions, and how boundary lifting can restore the +full convergence order. + +Two classes are provided: + +* :class:`fenics_allencahn_imex_timebc` — the Allen-Cahn equation with + naive time-dependent BC imposition (causes order reduction). +* :class:`fenics_allencahn_imex_timebc_lift` — the same equation + reformulated via boundary lifting so that the transformed variable + satisfies **homogeneous** BCs (restores full SDC order). +""" + +import logging + +import dolfin as df +import numpy as np + +from pySDC.core.problem import Problem +from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh + + +class fenics_allencahn_imex_timebc(Problem): + r""" + One-dimensional Allen-Cahn equation on :math:`[0, 1]` with inhomogeneous + time-dependent Dirichlet boundary conditions, solved with an IMEX split + in space (diffusion implicit, nonlinear reaction explicit) and FEniCS as + the spatial discretisation. + + **Equation** + + .. math:: + \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + - \frac{2}{\varepsilon^2} u (1 - u)(1 - 2u) - 6 d_w u (1 - u), + + with exact (traveling-wave) solution + + .. math:: + u(x, t) = \tfrac{1}{2}\!\left(1 + + \tanh\!\left(\frac{x - 0.5 - v\,t}{\sqrt{2}\,\varepsilon}\right)\right), + \qquad v = 3\sqrt{2}\,\varepsilon\,d_w. + + The spatial domain is mapped from :math:`[-0.5, 0.5]` to :math:`[0, 1]` + by the substitution :math:`x \to x - 0.5`. + + **Boundary conditions** are inhomogeneous and time-dependent (exact + solution evaluated at :math:`x = 0` and :math:`x = 1`). + + **Order reduction mechanism** + + In :meth:`solve_system`, the time-dependent Dirichlet BC is enforced via + + .. code-block:: python + + self.bc.apply(T, b.values.vector()) + + Dolfin's ``apply(A, b)`` call modifies the system matrix *and* overwrites + the boundary DOFs of ``b`` with the Dirichlet value :math:`u_D(t)`. This + replaces whatever the SDC iteration had accumulated at the boundary nodes + with the exact data at the *new* time, disrupting the implicit sweeper's + fixed-point equation and leading to **SDC order reduction**. + + **IMEX split** + + * Implicit part: :math:`f_{\mathrm{impl}} = K u` + (assembled stiffness matrix applied to :math:`u`). + * Explicit part: :math:`f_{\mathrm{expl}} = M(-N(u))` + where :math:`N(u) = \frac{2}{\varepsilon^2} u(1-u)(1-2u) + 6 d_w u(1-u)` + is the positive reaction function, so that the PDE reads + :math:`u_t = u_{xx} - N(u)`. + + Parameters + ---------- + c_nvars : int, optional + Number of mesh intervals (spatial resolution). Default ``128``. + t0 : float, optional + Starting time. Default ``0.0``. + family : str, optional + FEM element family. Default ``'CG'``. + order : int, optional + FEM element polynomial order. Default ``4``. + refinements : int, optional + Number of mesh refinements. Default ``1``. + eps : float, optional + Interface parameter :math:`\varepsilon`. Default ``0.3``. + dw : float, optional + Driving force :math:`d_w`. Default ``-0.04``. + + Attributes + ---------- + V : dolfin.FunctionSpace + M : dolfin.Matrix + Mass matrix :math:`\int_\Omega u\,v\,dx`. + K : dolfin.Matrix + Stiffness matrix :math:`-\int_\Omega \nabla u \cdot \nabla v\,dx`. + bc : dolfin.DirichletBC + Time-dependent inhomogeneous Dirichlet BC. + bc_hom : dolfin.DirichletBC + Homogeneous Dirichlet BC (used for residual fixing and mass inversion). + v_speed : float + Front propagation speed :math:`v = 3\sqrt{2}\,\varepsilon\,d_w`. + fix_bc_for_residual : bool + Always ``True``; instructs the sweeper to call :meth:`fix_residual`. + """ + + dtype_u = fenics_mesh + dtype_f = rhs_fenics_mesh + + def __init__( + self, + c_nvars=128, + t0=0.0, + family='CG', + order=4, + refinements=1, + eps=0.3, + dw=-0.04, + ): + """Initialisation routine.""" + + def Boundary(x, on_boundary): + return on_boundary + + logging.getLogger('FFC').setLevel(logging.WARNING) + logging.getLogger('UFL').setLevel(logging.WARNING) + + df.parameters["form_compiler"]["optimize"] = True + df.parameters["form_compiler"]["cpp_optimize"] = True + df.parameters['allow_extrapolation'] = True + + mesh = df.UnitIntervalMesh(c_nvars) + for _ in range(refinements): + mesh = df.refine(mesh) + + self.V = df.FunctionSpace(mesh, family, order) + tmp = df.Function(self.V) + print('DoFs on this level:', len(tmp.vector()[:])) + + super().__init__(self.V) + self._makeAttributeAndRegister( + 'c_nvars', + 't0', + 'family', + 'order', + 'refinements', + 'eps', + 'dw', + localVars=locals(), + readOnly=True, + ) + + # Front speed: v = 3*sqrt(2)*eps*dw + self.v_speed = 3.0 * np.sqrt(2.0) * eps * dw + + # Assemble mass and stiffness matrices + u_trial = df.TrialFunction(self.V) + v_test = df.TestFunction(self.V) + self.M = df.assemble(u_trial * v_test * df.dx) + self.K = df.assemble(-1.0 * df.inner(df.nabla_grad(u_trial), df.nabla_grad(v_test)) * df.dx) + + # Time-dependent Dirichlet BC (exact traveling-wave evaluated at boundary) + self.u_D = df.Expression( + '0.5*(1 + tanh((x[0] - 0.5 - vsp*t) / (s2*eps)))', + vsp=self.v_speed, + t=t0, + s2=np.sqrt(2.0), + eps=eps, + degree=order, + ) + self.bc = df.DirichletBC(self.V, self.u_D, Boundary) + self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary) + + self.fix_bc_for_residual = True + + def solve_system(self, rhs, factor, u0, t): + r""" + Solve the linear system :math:`(M - \text{factor}\,K)\,u = \text{rhs}`. + + The time-dependent BC is imposed via + + .. code-block:: python + + self.bc.apply(T, b.values.vector()) + + Dolfin's ``apply(A, b)`` call simultaneously modifies the system matrix + *and* **overwrites** the boundary DOFs of ``b`` with the Dirichlet value + :math:`u_D(t)`. This overwrites whatever values the SDC iteration had + accumulated at the boundary nodes, breaking the fixed-point property of + the implicit sweep and causing **SDC order reduction**. + + Parameters + ---------- + rhs : fenics_mesh + factor : float + u0 : fenics_mesh + t : float + + Returns + ------- + fenics_mesh + """ + u = self.dtype_u(u0) + T = self.M - factor * self.K + b = self.dtype_u(rhs) + + # Update expression time. + # bc.apply(T, b) modifies the system matrix AND overwrites boundary DOFs + # of b with u_D(t) — this is the source of SDC order reduction. + self.u_D.t = t + self.bc.apply(T, b.values.vector()) + + df.solve(T, u.values.vector(), b.values.vector()) + + return u + + def eval_f(self, u, t): + r""" + Evaluate both parts of the right-hand side. + + * ``f.impl = K u`` — stiffness matrix applied to :math:`u`. + * ``f.expl = M(-N(u))`` — mass-matrix-weighted negated reaction, where + :math:`N(u) = \tfrac{2}{\varepsilon^2}u(1-u)(1-2u) + 6d_w u(1-u)` + so that the Allen-Cahn PDE reads :math:`u_t = u_{xx} - N(u)`. + + Parameters + ---------- + u : fenics_mesh + t : float + + Returns + ------- + rhs_fenics_mesh + """ + f = self.dtype_f(self.V) + + # Implicit part: stiffness applied to u + self.K.mult(u.values.vector(), f.impl.values.vector()) + + # Explicit part: -N(u), mass-weighted + # N(u) = (2/eps^2)*u*(1-u)*(1-2u) + 6*dw*u*(1-u) [positive reaction] + u_vec = u.values.vector()[:] + eps2 = self.eps**2 + pos_N = (2.0 / eps2) * u_vec * (1.0 - u_vec) * (1.0 - 2.0 * u_vec) + 6.0 * self.dw * u_vec * (1.0 - u_vec) + expl_fn = self.dtype_u(self.V) + expl_fn.values.vector()[:] = -pos_N + f.expl = self.apply_mass_matrix(expl_fn) + + return f + + def apply_mass_matrix(self, u): + r"""Apply the mass matrix :math:`M` to ``u``.""" + me = self.dtype_u(self.V) + self.M.mult(u.values.vector(), me.values.vector()) + return me + + def fix_residual(self, res): + """Zero out boundary DOFs of the residual (required by the sweeper).""" + self.bc_hom.apply(res.values.vector()) + return None + + def u_exact(self, t): + r""" + Return the exact traveling-wave solution at time :math:`t`. + + Parameters + ---------- + t : float + + Returns + ------- + fenics_mesh + """ + u0_expr = df.Expression( + '0.5*(1 + tanh((x[0] - 0.5 - vsp*t) / (s2*eps)))', + vsp=self.v_speed, + t=t, + s2=np.sqrt(2.0), + eps=self.eps, + degree=self.order, + ) + return self.dtype_u(df.interpolate(u0_expr, self.V), val=self.V) + + +class fenics_allencahn_imex_timebc_lift(fenics_allencahn_imex_timebc): + r""" + Same Allen-Cahn equation using **boundary lifting** to restore the full + SDC convergence order. + + **Idea** — decompose :math:`u = v + E` where the *lift* + + .. math:: + E(x, t) = u_L(t)\,(1 - x) + u_R(t)\,x + + interpolates linearly between the two time-dependent boundary values + + .. math:: + u_L(t) = \tfrac{1}{2}\!\left(1 + + \tanh\!\left(\frac{-0.5 - v\,t}{\sqrt{2}\,\varepsilon}\right)\right), + \qquad + u_R(t) = \tfrac{1}{2}\!\left(1 + + \tanh\!\left(\frac{0.5 - v\,t}{\sqrt{2}\,\varepsilon}\right)\right). + + The transformed variable :math:`v = u - E` satisfies **homogeneous** + Dirichlet BCs :math:`v(0, t) = v(1, t) = 0` at all times. Because + :meth:`solve_system` only applies ``bc_hom`` (no overwrite of non-zero + values), the SDC fixed-point property is preserved and the full + collocation order is recovered. + + **Transformed equation** + + Since the lift is linear in :math:`x` (:math:`E_{xx} = 0`), the equation + for :math:`v` reads + + .. math:: + v_t = v_{xx} - N(v + E) - E_t, + + where :math:`E_t(x, t) = u_L'(t)(1 - x) + u_R'(t)\,x`. + + **Modified IMEX split** + + * Implicit part: :math:`f_{\mathrm{impl}} = K v` (unchanged). + * Explicit part: :math:`f_{\mathrm{expl}} = M(-N(v + E) - E_t)`. + + Parameters + ---------- + Same as :class:`fenics_allencahn_imex_timebc`. + """ + + def __init__( + self, + c_nvars=128, + t0=0.0, + family='CG', + order=4, + refinements=1, + eps=0.3, + dw=-0.04, + ): + """Initialisation routine.""" + super().__init__(c_nvars, t0, family, order, refinements, eps, dw) + + # ------------------------------------------------------------------ + # helpers for the lift and its time derivative + # ------------------------------------------------------------------ + + def _get_lift_values(self, t): + r""" + Nodal values of :math:`E(x, t)` at the DOF coordinates. + + Returns + ------- + numpy.ndarray shape ``(n_dofs,)`` + """ + vsp, eps, s2 = self.v_speed, self.eps, np.sqrt(2.0) + uL = 0.5 * (1.0 + np.tanh((-0.5 - vsp * t) / (s2 * eps))) + uR = 0.5 * (1.0 + np.tanh((0.5 - vsp * t) / (s2 * eps))) + x = self.V.tabulate_dof_coordinates()[:, 0] + return uL * (1.0 - x) + uR * x + + def _get_lift_dt_values(self, t): + r""" + Nodal values of :math:`\partial_t E(x, t)` at the DOF coordinates. + + Returns + ------- + numpy.ndarray shape ``(n_dofs,)`` + """ + vsp, eps, s2 = self.v_speed, self.eps, np.sqrt(2.0) + sL = (-0.5 - vsp * t) / (s2 * eps) + sR = (0.5 - vsp * t) / (s2 * eps) + duL_dt = -vsp / (2.0 * s2 * eps * np.cosh(sL) ** 2) + duR_dt = -vsp / (2.0 * s2 * eps * np.cosh(sR) ** 2) + x = self.V.tabulate_dof_coordinates()[:, 0] + return duL_dt * (1.0 - x) + duR_dt * x + + # ------------------------------------------------------------------ + # overridden problem interface + # ------------------------------------------------------------------ + + def solve_system(self, rhs, factor, u0, t): + r""" + Solve :math:`(M - \text{factor}\,K)\,v = \text{rhs}` with + **homogeneous** Dirichlet BCs. + + Only ``bc_hom.apply(T, b.values.vector())`` is called — no boundary + DOFs are overwritten with non-zero values, so the SDC fixed-point + property for :math:`v` is unaffected. + + Parameters + ---------- + rhs : fenics_mesh + factor : float + u0 : fenics_mesh + t : float + + Returns + ------- + fenics_mesh + """ + u = self.dtype_u(u0) + T = self.M - factor * self.K + b = self.dtype_u(rhs) + + self.bc_hom.apply(T, b.values.vector()) + + df.solve(T, u.values.vector(), b.values.vector()) + + return u + + def eval_f(self, v, t): + r""" + Evaluate both parts of the right-hand side for the lifted variable. + + * ``f.impl = K v`` — stiffness applied to :math:`v`. + * ``f.expl = M(-N(v + E) - E_t)`` — mass-weighted negated reaction + :math:`N(v+E)` plus negated lift derivative :math:`E_t`, where + :math:`N(u) = \tfrac{2}{\varepsilon^2}u(1-u)(1-2u) + 6d_w u(1-u)`. + + Parameters + ---------- + v : fenics_mesh + t : float + + Returns + ------- + rhs_fenics_mesh + """ + f = self.dtype_f(self.V) + + # Implicit part: stiffness of v + self.K.mult(v.values.vector(), f.impl.values.vector()) + + # Lift and its time derivative at current time + E_vals = self._get_lift_values(t) + Et_vals = self._get_lift_dt_values(t) + + # u = v + E at DOF locations + u_vals = v.values.vector()[:] + E_vals + + # N(v+E) = (2/eps^2)*u*(1-u)*(1-2u) + 6*dw*u*(1-u) [positive reaction] + eps2 = self.eps**2 + pos_N = (2.0 / eps2) * u_vals * (1.0 - u_vals) * (1.0 - 2.0 * u_vals) + 6.0 * self.dw * u_vals * ( + 1.0 - u_vals + ) + + # Explicit forcing for v: -N(v+E) - E_t (matches lifted equation v_t = v_xx - N(v+E) - E_t) + expl_nodal = -pos_N - Et_vals + + expl_fn = self.dtype_u(self.V) + expl_fn.values.vector()[:] = expl_nodal + f.expl = self.apply_mass_matrix(expl_fn) + + return f + + def u_exact(self, t): + r""" + Exact solution of the transformed variable :math:`v = u - E` at + time :math:`t`. + + Parameters + ---------- + t : float + + Returns + ------- + fenics_mesh + """ + # Interpolate the exact traveling wave + u_expr = df.Expression( + '0.5*(1 + tanh((x[0] - 0.5 - vsp*t) / (s2*eps)))', + vsp=self.v_speed, + t=t, + s2=np.sqrt(2.0), + eps=self.eps, + degree=self.order, + ) + u_fn = df.interpolate(u_expr, self.V) + + # Subtract the lift: v = u - E + v_fn = df.Function(self.V) + v_fn.vector()[:] = u_fn.vector()[:] - self._get_lift_values(t) + + return self.dtype_u(v_fn, val=self.V) diff --git a/pySDC/playgrounds/FEniCS/allen_cahn_1d/run_convergence.py b/pySDC/playgrounds/FEniCS/allen_cahn_1d/run_convergence.py new file mode 100644 index 0000000000..c065eddecb --- /dev/null +++ b/pySDC/playgrounds/FEniCS/allen_cahn_1d/run_convergence.py @@ -0,0 +1,171 @@ +""" +Convergence study for SDC applied to the 1D Allen-Cahn equation with +inhomogeneous time-dependent Dirichlet boundary conditions, using FEM (FEniCS). + +Two test cases are compared: + +- **Naive** (``fenics_allencahn_imex_timebc``): + :math:`u(x,t) = \\tfrac{1}{2}(1 + \\tanh((x - 0.5 - vt)/(\\sqrt{2}\\,\\varepsilon)))` + with time-dependent Dirichlet BCs enforced via ``bc.apply(b.values.vector())`` + inside ``solve_system``. This is the standard FEniCS mechanism but it + **causes order reduction**: the observed convergence order is lower than the + theoretical SDC order :math:`2M - 1`. + +- **Lifted** (``fenics_allencahn_imex_timebc_lift``): + The solution is decomposed as :math:`u = v + E` where + :math:`E(x,t) = u_L(t)(1-x) + u_R(t)\\,x` is a linear lift matching the + time-dependent BCs. The transformed variable :math:`v` satisfies + **homogeneous** BCs, and ``solve_system`` only applies ``bc_hom``. This + **restores the full SDC order**. + +Usage +----- +Run directly:: + + 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_mass import imex_1st_order_mass +from pySDC.playgrounds.FEniCS.allen_cahn_1d.problem_classes import ( + fenics_allencahn_imex_timebc, + fenics_allencahn_imex_timebc_lift, +) + + +def build_description(problem_class, num_nodes, dt, t0=0.0, c_nvars=64, eps=0.3, dw=-0.04): + """ + Build the pySDC description dictionary for a single run. + + Parameters + ---------- + problem_class : type + Either :class:`fenics_allencahn_imex_timebc` or + :class:`fenics_allencahn_imex_timebc_lift`. + num_nodes : int + Number of SDC collocation nodes (RADAU-RIGHT). + dt : float + Time-step size. + t0 : float, optional + Start time. Default ``0.0``. + c_nvars : int, optional + Spatial degrees of freedom. Default ``64``. + eps : float, optional + Interface parameter :math:`\\varepsilon`. Default ``0.3``. + dw : float, optional + Driving force. Default ``-0.04``. + + Returns + ------- + description : dict + controller_params : dict + """ + level_params = {'restol': 1e-12, 'dt': dt} + step_params = {'maxiter': 50} + sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': num_nodes} + problem_params = { + 'c_nvars': c_nvars, + 't0': t0, + 'family': 'CG', + 'order': 4, + 'refinements': 1, + 'eps': eps, + 'dw': dw, + } + controller_params = {'logger_level': 30} + + description = { + 'problem_class': problem_class, + 'problem_params': problem_params, + 'sweeper_class': imex_1st_order_mass, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + } + return description, controller_params + + +def run_sdc(problem_class, dt, num_nodes=3, t0=0.0, Tend=1.0, c_nvars=64, eps=0.3, dw=-0.04): + """ + Run a single SDC solve and return the relative :math:`L^2` error at ``Tend``. + + Parameters + ---------- + problem_class : type + dt : float + num_nodes : int, optional + t0 : float, optional + Tend : float, optional + c_nvars : int, optional + eps : float, optional + dw : float, optional + + Returns + ------- + float + Relative error :math:`\\|u_h - u_{\\text{exact}}\\| / \\|u_{\\text{exact}}\\|` + at ``Tend``. + """ + description, controller_params = build_description( + problem_class, num_nodes, dt, t0=t0, c_nvars=c_nvars, eps=eps, dw=dw + ) + + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + uend, _ = controller.run(u0=uinit, t0=t0, Tend=Tend) + uex = P.u_exact(Tend) + return float(abs(uex - uend) / abs(uex)) + + +def compute_order(problem_class, dts, **kwargs): + """ + Compute errors and least-squares convergence order estimate. + + Parameters + ---------- + problem_class : type + dts : list of float + **kwargs + Forwarded to :func:`run_sdc`. + + Returns + ------- + errors : list of float + order : float + """ + errors = [run_sdc(problem_class, dt, **kwargs) for dt in dts] + order = float(np.polyfit(np.log(dts), np.log(errors), 1)[0]) + return errors, order + + +def main(): + """Run the convergence study and print results.""" + num_nodes = 3 + # Large-dt regime where temporal errors dominate over the FEM spatial error + dts = [0.5 / 2**k for k in range(3)] + Tend = 1.0 + + print("=" * 70) + print("SDC convergence study: order reduction with time-dependent BCs") + print(" Allen-Cahn equation in FEniCS (IMEX, diffusion implicit,") + print(f" reaction explicit), RADAU-RIGHT M={num_nodes} nodes, Tend={Tend}") + print("=" * 70) + + for problem_class, label in [ + (fenics_allencahn_imex_timebc, "Naive (time-dependent BCs, order reduction)"), + (fenics_allencahn_imex_timebc_lift, "Lifted (boundary lifting, full order restored)"), + ]: + errors, order = compute_order(problem_class, dts, num_nodes=num_nodes, Tend=Tend) + print(f"\n{label}") + print(f" Estimated order : {order:.2f} (expected ~{2*num_nodes-1} without reduction)") + for dt, err in zip(dts, errors): + print(f" dt = {dt:.5f} rel. error = {err:.3e}") + + print("\n" + "=" * 70) + + +if __name__ == '__main__': + main() diff --git a/pySDC/playgrounds/FEniCS/allen_cahn_1d/test_order_reduction.py b/pySDC/playgrounds/FEniCS/allen_cahn_1d/test_order_reduction.py new file mode 100644 index 0000000000..9c1eb8f452 --- /dev/null +++ b/pySDC/playgrounds/FEniCS/allen_cahn_1d/test_order_reduction.py @@ -0,0 +1,118 @@ +""" +Tests verifying SDC order-reduction behaviour with inhomogeneous time-dependent +Dirichlet boundary conditions for the 1D Allen-Cahn equation, using FEniCS. + +Two scenarios are tested: + +1. **Naive** (``fenics_allencahn_imex_timebc``, time-dependent BCs): + The standard FEniCS BC imposition via ``bc.apply(b.values.vector())`` in + ``solve_system`` causes order reduction — the observed convergence order + is lower than the theoretical SDC order :math:`2M - 1`. + +2. **Lifted** (``fenics_allencahn_imex_timebc_lift``): + Boundary lifting decomposes :math:`u = v + E`, where :math:`E` is a + linear lift matching the time-dependent BCs. The transformed variable + :math:`v` satisfies homogeneous BCs and SDC applied to :math:`v` + **restores the full convergence order**. + +All tests are marked ``@pytest.mark.fenics`` because they require FEniCS/dolfin. + +Parameters +---------- +_DTS : list of float + Three step sizes in the large-dt asymptotic regime ``[0.5, 0.25, 0.125]``. + At these step sizes temporal errors dominate over the high-order FEM + spatial discretisation error (CG-4 on a refined mesh). +_TEND : float + End time ``1.0``. +_C_NVARS : int + Coarse spatial resolution ``32`` — keeps tests fast while preserving + the temporal convergence signal. +""" + +import numpy as np +import pytest + +# dt values in the large-dt asymptotic regime where temporal errors dominate +_DTS = [0.5 / 2**k for k in range(3)] +_TEND = 1.0 +_C_NVARS = 32 # coarse mesh keeps tests fast; CG-4 spatial error is negligible + + +@pytest.mark.fenics +def test_allencahn_order_reduction(): + """ + Naive Allen-Cahn (time-dependent BCs) must exhibit order reduction. + + The observed convergence order with ``fenics_allencahn_imex_timebc`` must + be strictly less than the theoretical SDC order :math:`2M - 1 - 0.3`, + confirming that the standard FEniCS BC imposition causes a measurable + loss of temporal accuracy. + + For RADAU-RIGHT with M=3 nodes the theoretical order is 5; the observed + order should be noticeably lower. + """ + from pySDC.playgrounds.FEniCS.allen_cahn_1d.run_convergence import compute_order + from pySDC.playgrounds.FEniCS.allen_cahn_1d.problem_classes import fenics_allencahn_imex_timebc + + num_nodes = 3 + _, order = compute_order(fenics_allencahn_imex_timebc, _DTS, num_nodes=num_nodes, Tend=_TEND, c_nvars=_C_NVARS) + + theoretical_order = 2 * num_nodes - 1 # = 5 for M=3 + assert order < theoretical_order - 0.3, ( + f"Naive timebc case (M={num_nodes}): expected order < {theoretical_order - 0.3:.1f} " + f"(order reduction), got {order:.2f}" + ) + + +@pytest.mark.fenics +@pytest.mark.parametrize("num_nodes", [2, 3]) +def test_lifting_restores_full_order(num_nodes): + """ + Allen-Cahn with boundary lifting must recover the full SDC order. + + The lifting approach (``fenics_allencahn_imex_timebc_lift``) transforms + the problem so that the solver always sees homogeneous BCs. This + eliminates the source of order reduction and the observed convergence + order must be at least :math:`2M - 2`. + """ + from pySDC.playgrounds.FEniCS.allen_cahn_1d.run_convergence import compute_order + from pySDC.playgrounds.FEniCS.allen_cahn_1d.problem_classes import fenics_allencahn_imex_timebc_lift + + _, order = compute_order( + fenics_allencahn_imex_timebc_lift, _DTS, num_nodes=num_nodes, Tend=_TEND, c_nvars=_C_NVARS + ) + + expected = 2 * num_nodes - 1 + assert order >= expected - 1.0, ( + f"Lifting case (M={num_nodes}): expected order >= {expected - 1.0:.1f} " + f"(full order restored), got {order:.2f}" + ) + + +@pytest.mark.fenics +def test_lifting_has_higher_order_than_naive(): + """ + The boundary lifting case must converge faster than the naive case. + + This verifies that the lifting correction actually fixes the order + reduction. With RADAU-RIGHT M=3 nodes, the lifted order must exceed + the naive order by at least 0.3. + """ + from pySDC.playgrounds.FEniCS.allen_cahn_1d.run_convergence import compute_order + from pySDC.playgrounds.FEniCS.allen_cahn_1d.problem_classes import ( + fenics_allencahn_imex_timebc, + fenics_allencahn_imex_timebc_lift, + ) + + num_nodes = 3 + _, order_naive = compute_order( + fenics_allencahn_imex_timebc, _DTS, num_nodes=num_nodes, Tend=_TEND, c_nvars=_C_NVARS + ) + _, order_lift = compute_order( + fenics_allencahn_imex_timebc_lift, _DTS, num_nodes=num_nodes, Tend=_TEND, c_nvars=_C_NVARS + ) + + assert order_lift > order_naive + 0.3, ( + f"Lifting order ({order_lift:.2f}) should be > naive timebc order ({order_naive:.2f}) + 0.3" + )