DeltaSDC: restore the nonlinear FEniCS route now that Gray-Scott works

problems_fenics.py was deleted in the previous commit because it could not
be executed: fenics_grayscott did not run at all against the pinned DOLFIN.
With that fixed, the module can be validated, so it comes back.

fenics_grayscott_delta poses the correction problem variationally, in the
same shape as the stock solve_system, and expands the Gray-Scott reaction
term analytically:

    (w1+d1)(w2+d2)^2 - w1 w2^2 = w1(2 w2 d2 + d2^2) + d1(w2 + d2)^2

so every term carries an explicit factor delta. Assembling both forms and
subtracting instead would cancel two O(|F|) vectors and reinstate an
absolute error of order eps*|F|, which is the failure mode the correction
form exists to avoid.

Verified against the stock path: 3.08e-11 at full precision and 3.09e-11
with emulated float32, both consistent with Gray-Scott's own Newton
tolerance of 1e-9 rather than with any error in the formulation.

The fenics-marked test now covers both routes on this backend: the linear
implicit part through fenics_heat and the IMEX delta form, which needs no
problem-class change, and the nonlinear implicit part through Gray-Scott.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: Robert Speck

pySDC/projects/DeltaSDC/README.rst

@@ -48,6 +48,7 @@ Contents
 ``sweepers_MPI.py``          ``delta_implicit_MPI``, one collocation node per rank
 ``problems.py``              ``allencahn_delta``, an example nonlinear correction solve
 ``problems_petsc.py``        ``petsc_fisher_delta`` (reduced precision emulated)
+``problems_fenics.py``       ``fenics_grayscott_delta`` (reduced precision emulated)
 ``run_demo.py``              runnable demonstration
 ``run_petsc.py``             PETSc entry point, called by the ``petsc``-marked test
 ``run_fenics.py``            FEniCS entry point, called by the ``fenics``-marked test
@@ -67,7 +68,7 @@ MSSDC (single level, multi-step)  covered
 Node-parallel sweeper        covered, spawns ``mpirun`` with one rank per node
 IMEX                         covered
 PETSc                        covered, reduced precision **emulated**
-FEniCS                       covered, on ``fenics_heat`` via the IMEX delta form
+FEniCS                       covered, linear via ``fenics_heat``, nonlinear via Gray-Scott
 ===========================  ====================================================================
 
 PETSc and DOLFIN fix their scalar type at build time, so genuine single precision needs a

pySDC/projects/DeltaSDC/problems_fenics.py

@@ -0,0 +1,163 @@
+r"""
+FEniCS problem with a node-local correction solve.
+
+Kept in its own module because importing it requires ``dolfin``, which :mod:`problems` does not.
+
+The correction problem is posed variationally, in the same shape as the stock ``solve_system``:
+find :math:`\delta` such that
+
+.. math::
+    \langle \delta, q\rangle - factor\,\big[F(w+\delta; q) - F(w; q)\big] = \langle r, q\rangle .
+
+The increment :math:`F(w+\delta) - F(w)` is **expanded analytically** rather than written as a
+difference of two assembled forms. Assembling both and subtracting would cancel two
+:math:`\mathcal{O}(|F|)` vectors and reinstate an absolute error of order
+:math:`\varepsilon |F|`, which is exactly the failure mode the correction form exists to avoid.
+
+For Gray-Scott, with :math:`w = (w_1, w_2)` the base and :math:`\delta = (\delta_1, \delta_2)`:
+
+* diffusion is linear, so its increment is the same form evaluated at :math:`\delta`;
+* :math:`A(1 - u_1) \to -A\delta_1` and :math:`B u_2 \to B\delta_2`;
+* the reaction term expands as
+
+  .. math::
+      (w_1+\delta_1)(w_2+\delta_2)^2 - w_1 w_2^2
+        = w_1\left(2 w_2 \delta_2 + \delta_2^2\right) + \delta_1 (w_2 + \delta_2)^2 ,
+
+  in which every term carries an explicit factor :math:`\delta`.
+
+Reduced precision is **emulated**: DOLFIN inherits PETSc's build-time scalar type, so values are
+rounded through the requested working precision and written back, capping the *information* while
+the arithmetic stays at the backend type. That is optimistic about iteration counts and attainable
+accuracy compared with a real single-precision build.
+"""
+
+import dolfin as df
+import numpy as np
+
+from pySDC.implementations.problem_classes.GrayScott_1D_FEniCS_implicit import fenics_grayscott
+
+
+def quantize_function(function, work_precision):
+    """
+    Round a DOLFIN function's coefficients through the working precision, in place.
+
+    Parameters
+    ----------
+    function : dolfin.Function
+        Function whose vector is quantized.
+    work_precision : numpy.dtype or None
+        Working precision. ``None`` leaves the function untouched.
+
+    Returns
+    -------
+    dolfin.Function
+        The same function, for convenience.
+    """
+    if work_precision is None:
+        return function
+    values = function.vector().get_local()
+    function.vector().set_local(values.astype(np.dtype(work_precision)).astype(values.dtype))
+    function.vector().apply('insert')
+    return function
+
+
+class fenics_grayscott_delta(fenics_grayscott):
+    """
+    Gray-Scott exposing ``solve_system_delta`` alongside the stock ``solve_system``.
+
+    Parameters
+    ----------
+    solve_precision : dtype-like or None, optional
+        Working precision to emulate for the node-local solve. ``None`` keeps backend precision.
+    **kwargs
+        Forwarded to :class:`fenics_grayscott`.
+    """
+
+    def __init__(self, solve_precision=None, **kwargs):
+        """Initialization routine"""
+        super().__init__(**kwargs)
+        self.solve_precision = None if solve_precision is None else np.dtype(solve_precision)
+
+        # base state of the correction, assigned per solve
+        self.base = df.Function(self.V)
+        self.delta = df.Function(self.V)
+
+    def _increment_forms(self, test_functions):
+        r"""
+        Build the analytically expanded weak form of :math:`F(w+\delta) - F(w)`.
+
+        Parameters
+        ----------
+        test_functions : tuple
+            The two test functions of the mixed space.
+
+        Returns
+        -------
+        ufl.Form
+            The increment form, every term carrying an explicit factor of the correction.
+        """
+        q1, q2 = test_functions
+        b1, b2 = df.split(self.base)
+        d1, d2 = df.split(self.delta)
+
+        # (b1+d1)(b2+d2)^2 - b1 b2^2, expanded so nothing cancels
+        reaction = b1 * (2 * b2 * d2 + d2**2) + d1 * (b2 + d2) ** 2
+
+        increment1 = (
+            -self.Du * df.inner(df.nabla_grad(d1), df.nabla_grad(q1)) - reaction * q1 - self.A * d1 * q1
+        ) * df.dx
+        increment2 = (
+            -self.Dv * df.inner(df.nabla_grad(d2), df.nabla_grad(q2)) + reaction * q2 - self.B * d2 * q2
+        ) * df.dx
+        return increment1 + increment2
+
+    def solve_system_delta(self, r, factor, base, f_base, t):
+        r"""
+        Solve :math:`\delta - factor\,[f(base+\delta) - f(base)] = r` for the correction.
+
+        Parameters
+        ----------
+        r : dtype_u
+            Right-hand side of the correction equation.
+        factor : float
+            Implicit prefactor assembled by the sweeper.
+        base : dtype_u
+            Base state :math:`w`.
+        f_base : dtype_f
+            ``f`` evaluated at ``base``; accepted so no extra evaluation is needed.
+        t : float
+            Physical time, accepted for interface compatibility.
+
+        Returns
+        -------
+        dtype_u
+            The correction.
+        """
+        self.base.assign(base.values)
+        self.delta.assign(df.Function(self.V))  # start from the zero correction
+
+        q1, q2 = df.TestFunctions(self.V)
+        d1, d2 = df.split(self.delta)
+        r1, r2 = df.split(r.values)
+
+        residual = (d1 * q1 + d2 * q2) * df.dx - factor * self._increment_forms((q1, q2))
+        residual -= (r1 * q1 + r2 * q2) * df.dx
+
+        trial = df.TrialFunction(self.V)
+        jacobian = df.derivative(residual, self.delta, trial)
+
+        problem = df.NonlinearVariationalProblem(residual, self.delta, [], jacobian)
+        solver = df.NonlinearVariationalSolver(problem)
+        prm = solver.parameters['newton_solver']
+        prm['absolute_tolerance'] = 1e-09
+        prm['relative_tolerance'] = 1e-08
+        prm['maximum_iterations'] = 100
+        prm['relaxation_parameter'] = 1.0
+        solver.solve()
+
+        quantize_function(self.delta, self.solve_precision)
+
+        me = self.dtype_u(self.V)
+        me.values.assign(self.delta)
+        return me

pySDC/projects/DeltaSDC/run_fenics.py

@@ -4,8 +4,12 @@
 Runnable entry point in the style of ``pySDC/tutorial/step_7/A_pySDC_with_FEniCS.py``: the logic
 lives here and the FEniCS-marked test simply calls :func:`main`.
 
-Uses the FEniCS heat equation with the IMEX delta-form sweeper, which reproduces the stock
-``imex_1st_order`` path exactly.
+Covers both routes on this backend:
+
+* the **linear implicit** part, via the FEniCS heat equation and the IMEX delta form, which
+  reproduces the stock ``imex_1st_order`` path exactly and needs no problem-class change;
+* the **nonlinear implicit** part, via Gray-Scott and ``fenics_grayscott_delta``, which supplies a
+  variational ``solve_system_delta`` with an analytically expanded increment.
 
 Note what this backend showed about ``linear_implicit=True``: that shortcut reuses the stock
 ``solve_system`` to solve the correction equation directly, which additionally requires the solve to
@@ -24,9 +28,12 @@
 import numpy as np
 
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+from pySDC.implementations.problem_classes.GrayScott_1D_FEniCS_implicit import fenics_grayscott
 from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import fenics_heat
+from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
 from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
-from pySDC.projects.DeltaSDC.sweepers import delta_imex_1st_order
+from pySDC.projects.DeltaSDC.problems_fenics import fenics_grayscott_delta
+from pySDC.projects.DeltaSDC.sweepers import delta_imex_1st_order, delta_implicit
 
 T0 = 0.0
 
@@ -91,9 +98,96 @@ def run(sweeper_class, sweeper_params_, maxiter=6, dt=0.2, nsteps=2):
     return uend
 
 
+GRAYSCOTT_PARAMS = {
+    'c_nvars': 64,
+    't0': 0.0,
+    'family': 'CG',
+    'order': 2,
+    'refinements': 1,
+    'Du': 1.0,
+    'Dv': 0.01,
+    'A': 0.09,
+    'B': 0.086,
+}
+
+
+def run_grayscott(problem_class, problem_params, sweeper_class, maxiter=4, dt=1.0, nsteps=2):
+    """
+    Run a short Gray-Scott simulation.
+
+    Parameters
+    ----------
+    problem_class : type
+        Problem class to integrate.
+    problem_params : dict
+        Parameters for the problem class.
+    sweeper_class : type
+        Sweeper class to use.
+    maxiter : int, optional
+        Number of SDC iterations.
+    dt : float, optional
+        Step size.
+    nsteps : int, optional
+        Number of steps.
+
+    Returns
+    -------
+    dtype_u
+        The end value.
+    """
+    description = {
+        'problem_class': problem_class,
+        'problem_params': problem_params,
+        'sweeper_class': sweeper_class,
+        'sweeper_params': {
+            'quad_type': 'RADAU-RIGHT',
+            'node_type': 'LEGENDRE',
+            'num_nodes': 3,
+            'QI': 'LU',
+            'initial_guess': 'spread',
+        },
+        'level_params': {'restol': -1, 'dt': dt},
+        'step_params': {'maxiter': maxiter},
+    }
+    controller = controller_nonMPI(num_procs=1, controller_params={'logger_level': 40}, description=description)
+    prob = controller.MS[0].levels[0].prob
+    uend, _ = controller.run(u0=prob.u_exact(0.0), t0=0.0, Tend=nsteps * dt)
+    return uend
+
+
+def check_nonlinear():
+    """
+    Check the variational correction solve on Gray-Scott, at full and emulated reduced precision.
+
+    Returns
+    -------
+    dict
+        The three end values.
+
+    Raises
+    ------
+    AssertionError
+        If either delta-form variant deviates from the stock path.
+    """
+    results = {}
+    results['stock'] = run_grayscott(fenics_grayscott, GRAYSCOTT_PARAMS, generic_implicit)
+    results['delta'] = run_grayscott(fenics_grayscott_delta, GRAYSCOTT_PARAMS, delta_implicit)
+    results['delta_fp32'] = run_grayscott(
+        fenics_grayscott_delta, dict(GRAYSCOTT_PARAMS, solve_precision=np.float32), delta_implicit
+    )
+
+    # Gray-Scott's own Newton runs to absolute_tolerance 1e-9, which sets the agreement level here
+    for label in ['delta', 'delta_fp32']:
+        deviation = abs(results[label] - results['stock'])
+        print(f'{label:>22} vs stock: {deviation:.3e} (tol 1e-08)')
+        assert deviation < 1e-8, f'{label} deviates from the stock path by {deviation:.3e}'
+
+    return results
+
+
 def main():
     """
-    Check the delta form against the stock IMEX sweeper on a FEniCS backend.
+    Check the delta form against the stock sweepers on a FEniCS backend.
 
     Returns
     -------
@@ -122,6 +216,8 @@ def main():
     else:  # pragma: no cover - guards against a silent regression
         raise AssertionError('correction_precision should not be silently accepted for fenics_mesh')
 
+    check_nonlinear()
+
     return results
 
 

pySDC/tests/test_problems/test_GrayScott_1D_FEniCS.py

@@ -74,9 +74,7 @@ def test_runs_in_sdc():
         'level_params': {'restol': -1, 'dt': 1.0},
         'step_params': {'maxiter': 4},
     }
-    controller = controller_nonMPI(
-        num_procs=1, controller_params={'logger_level': 30}, description=description
-    )
+    controller = controller_nonMPI(num_procs=1, controller_params={'logger_level': 30}, description=description)
     prob = controller.MS[0].levels[0].prob
     u0 = prob.u_exact(0.0)
     uend, _ = controller.run(u0=u0, t0=0.0, Tend=2.0)