From 70046c8fc3639e4517c956d33780f369b51672e9 Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Thu, 10 Sep 2026 16:57:10 +0200 Subject: [PATCH 1/7] Add Taylor-Green order-reduction benchmark using FEniCS Reimplements the benchmark from PR #650 with legacy FEniCS in the existing StroemungsRaum project instead of FEniCSx in a new project. Periodicity comes from a constrained_domain on the function space, so dolfinx_mpc is not needed, and the project is already covered by the CI matrix and the `fenics` marker. The manufactured solution and forcing term of #650 were verified symbolically and carried over unchanged; the implementation issues raised in review were not: - No mass matrix inversion. The right-hand side is subtracted from the Newton residual as a vector, so the Mf^-1 / Mf round trip in #650 (and in the merged NavierStokes_2D_monolithic_FEniCS.py) disappears together with its per-node solve and its accuracy floor. - The residual is only fixed on the actual Dirichlet boundary. In #650 the periodic variant zeroed the residual on the periodic boundaries too, which let it stop iterating earlier than the Dirichlet variant and biased exactly the comparison the benchmark exists for. - The form, its Jacobian, the Newton problem and the solver are built once, with the step size as a Constant, rather than rebuilt per solve_system call. - Scaling by `factor` instead of dividing by it, so there is no factor == 0 case and the residual scale, and with it snes/newton tolerances, no longer depends on dt. - Domain extents are not derived from rank-local coordinates. Adds the order study that was missing: with periodic conditions the method attains its design order 2M-1 = 5, with time-dependent Dirichlet conditions it does not. At CI-affordable resolution the order gap is modest (~4.7 vs ~4.3 in the pressure, the Radau IIA reduction from 2M-1 to M+1); the error constant gap is close to an order of magnitude, and that is what the test asserts. Co-Authored-By: Claude Opus 5 --- pySDC/projects/StroemungsRaum/README.rst | 11 +- ...Stokes_2D_TaylorGreen_monolithic_FEniCS.py | 368 ++++++++++++++++++ .../run_Navier_Stokes_TaylorGreen_FEniCS.py | 243 ++++++++++++ .../test_Navier_Stokes_TaylorGreen_FEniCS.py | 181 +++++++++ 4 files changed, 802 insertions(+), 1 deletion(-) create mode 100644 pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py create mode 100644 pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py create mode 100644 pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py diff --git a/pySDC/projects/StroemungsRaum/README.rst b/pySDC/projects/StroemungsRaum/README.rst index c66176ca33..33054ea1ab 100644 --- a/pySDC/projects/StroemungsRaum/README.rst +++ b/pySDC/projects/StroemungsRaum/README.rst @@ -30,7 +30,7 @@ Implemented examples and test cases include: - Heat equation - Convection–diffusion and nonlinear convection–diffusion problems -- Incompressible Navier–Stokes equations, using: +- Incompressible Navier–Stokes equations, using: - Projection methods - Monolithic formulations - DAE- and PDE sweepers @@ -38,6 +38,15 @@ Implemented examples and test cases include: These serve as benchmarks and demonstrators for scalable space–time CFD simulations. +Order reduction from time-dependent boundary conditions +------------------------------------------------------- +``run_Navier_Stokes_TaylorGreen_FEniCS.py`` runs a manufactured Taylor–Green +solution on :math:`[-0.5, 0.5]^2` that is exactly one-periodic in :math:`x` and +constant on the top and bottom boundary. The *same* solution can therefore be +computed with time-dependent Dirichlet conditions in :math:`x` or with periodic +ones, and the difference in the observed temporal order isolates the order +reduction caused by the time-dependent boundary data alone. + Funding ------- Funded by the **German Federal Ministry of Education and Research (BMBF)** under diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py new file mode 100644 index 0000000000..62160ad78f --- /dev/null +++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py @@ -0,0 +1,368 @@ +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 + + +class _PeriodicX(df.SubDomain): + """ + Identifies the right boundary :math:`x = 0.5` with the left boundary :math:`x = -0.5`. + + Passed to ``FunctionSpace`` as ``constrained_domain``, so periodicity is handled by + the dof map itself and the periodic dofs never enter the linear system. + """ + + def inside(self, x, on_boundary): + return bool(df.near(x[0], -0.5) and on_boundary) + + def map(self, x, y): + y[0] = x[0] - 1.0 + y[1] = x[1] + + +class _NewtonStep(df.NonlinearProblem): + r""" + Newton problem for the node-to-node step :math:`M w + \Delta t_{QI} N(w) = rhs`. + + ``rhs`` is kept as an assembled vector and subtracted from the residual here, rather + than being written into the variational form as :math:`\int_\Omega rhs \cdot v\,dx`. + That form would apply the mass matrix to it, which would then have to be undone by a + mass matrix solve beforehand -- an exact round trip that costs a solve per node per + sweep and caps the attainable accuracy at the tolerance of that solve. + + Parameters + ---------- + F : UFL form + Residual form, without the right-hand side term. + J : UFL form + Jacobian of ``F``. + bcs : list of DirichletBC + Boundary conditions, applied in residual form (``bc.apply(b, x)``). + """ + + def __init__(self, F, J, bcs): + super().__init__() + self.F_form = F + self.J_form = J + self.bcs = bcs + self.rhs = None + + def F(self, b, x): + df.assemble(self.F_form, tensor=b) + b.axpy(-1.0, self.rhs) + for bc in self.bcs: + bc.apply(b, x) + + def J(self, A, x): + df.assemble(self.J_form, tensor=A) + for bc in self.bcs: + bc.apply(A) + + +class fenics_NSE_2D_TaylorGreen(Problem): + r""" + Forced two-dimensional incompressible Navier-Stokes equations on :math:`\Omega = [-0.5, 0.5]^2`, + set up to expose the order reduction caused by time-dependent Dirichlet boundary conditions. + + .. math:: + \frac{\partial u}{\partial t} = - u \cdot \nabla u + \nu \Delta u - \nabla p + g, + \qquad \nabla \cdot u = 0 + + The forcing :math:`g` is manufactured from the analytical solution + + .. math:: + u(x, y, t) &= 1 - e^{-8\pi^2\nu t}\sin(2\pi(x - t))\sin(\pi y)\cos(\pi y) \\ + v(x, y, t) &= - e^{-8\pi^2\nu t}\cos(2\pi(x - t))\cos^2(\pi y) \\ + p(x, y, t) &= 1 + \frac{4}{17}e^{-16\pi^2\nu t}\cos(4\pi(x - t))\cos(\pi y) + + which is divergence free and exactly one-periodic in :math:`x`. On :math:`y = \pm 0.5` it + collapses to the constants :math:`u = (1, 0)`, :math:`p = 1`, so the top and bottom boundary + data is time-independent. + + Because the solution is genuinely periodic in :math:`x`, the *same* exact solution satisfies + both variants selected by ``periodic``: + + - ``periodic=False``: time-dependent Dirichlet conditions on :math:`x = \pm 0.5`, + - ``periodic=True``: periodic conditions on :math:`x = \pm 0.5`. + + The only difference between the two runs is therefore the presence of time-dependent boundary + data, which is what isolates the order reduction. + + The problem is discretized in space with Taylor-Hood elements on a mixed velocity-pressure + space and solved monolithically, so the semi-discrete system is the differential-algebraic + system :math:`M \dot{w} = f(w, t)` with the singular mass matrix :math:`M = \mathrm{diag}(M_v, 0)`. + It therefore requires ``generic_implicit_mass`` as sweeper, which applies :math:`M` where needed + instead of inverting it. + + Parameters + ---------- + nelems : int, optional + Number of elements per spatial direction. + t0 : float, optional + Starting time. + order : int, optional + Polynomial degree of the velocity space; the pressure space uses ``order - 1``. + nu : float, optional + Kinematic viscosity :math:`\nu`. + periodic : bool, optional + Use periodic instead of time-dependent Dirichlet conditions on :math:`x = \pm 0.5`. + Sol_tol : float, optional + Absolute tolerance for the Newton solver. + + Attributes + ---------- + V : FunctionSpace + Velocity space. + Q : FunctionSpace + Pressure space. + W : FunctionSpace + Mixed velocity-pressure space. + M : Matrix + The velocity mass matrix :math:`\mathrm{diag}(M_v, 0)` on the mixed space. + g : Expression + Manufactured forcing term. + bc : list of DirichletBC + Dirichlet boundary conditions, time-dependent unless ``periodic``. + bc_hom : list of DirichletBC + Homogeneous conditions on the Dirichlet part of the boundary only, used to fix the residual. + fix_bc_for_residual : bool + Flag indicating that the residual requires special treatment due to boundary conditions. + + References + ---------- + .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, + C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015). + """ + + dtype_u = fenics_mesh + dtype_f = fenics_mesh + + df.set_log_active(False) + + def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, Sol_tol=1e-10): + + # set logger level for FFC and dolfin + logging.getLogger('FFC').setLevel(logging.WARNING) + logging.getLogger('UFL').setLevel(logging.WARNING) + + # set solver and form parameters + df.parameters["form_compiler"]["optimize"] = True + df.parameters["form_compiler"]["cpp_optimize"] = True + + mesh = df.RectangleMesh(df.Point(-0.5, -0.5), df.Point(0.5, 0.5), nelems, nelems) + + # define function spaces (Taylor-Hood); periodicity is baked into the dof map + P2 = df.VectorElement("P", mesh.ufl_cell(), order) + P1 = df.FiniteElement("P", mesh.ufl_cell(), order - 1) + constraint = _PeriodicX() if periodic else None + self.W = df.FunctionSpace(mesh, df.MixedElement([P2, P1]), constrained_domain=constraint) + self.V = df.FunctionSpace(mesh, P2, constrained_domain=constraint) + self.Q = df.FunctionSpace(mesh, P1, constrained_domain=constraint) + + super().__init__(self.W) + self._makeAttributeAndRegister( + 'nelems', 't0', 'order', 'nu', 'periodic', 'Sol_tol', localVars=locals(), readOnly=True + ) + + self.logger.debug('DoFs on this level: %d', self.W.dim()) + + # trial and test functions on the mixed space + self.u, self.p = df.TrialFunctions(self.W) + self.v, self.q = df.TestFunctions(self.W) + + # velocity mass matrix on the mixed space, i.e. diag(M_v, 0) + self.M = df.assemble(df.inner(self.u, self.v) * df.dx) + + # manufactured solution and the forcing term derived from it + self.u_ex = df.Expression( + ( + '1.0 - exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])', + '-exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])', + ), + pi=np.pi, + nu=nu, + t=t0, + degree=order + 2, + ) + self.p_ex = df.Expression( + '1.0 + (4.0/17.0)*exp(-16*pi*pi*nu*t)*cos(4*pi*(x[0] - t))*cos(pi*x[1])', + pi=np.pi, + nu=nu, + t=t0, + degree=order + 2, + ) + self.g = df.Expression( + ( + 'pi/34.0*exp(-16*pi*pi*nu*t)*sin(4*pi*(t - x[0]))*cos(pi*x[1])*(32.0 - 17.0*cos(pi*x[1]))', + '2*pi*pi*nu*exp(-8*pi*pi*nu*t)*cos(2*pi*(t - x[0]))' + ' - pi*exp(-16*pi*pi*nu*t)*sin(pi*x[1])' + '*(2.0*pow(cos(pi*x[1]), 3) + 4.0/17.0*cos(4*pi*(t - x[0])))', + ), + pi=np.pi, + nu=nu, + t=t0, + degree=order + 2, + ) + + # on y = +-0.5 the exact solution is constant in space and time + top_bottom = 'near(x[1], -0.5) || near(x[1], 0.5)' + left_right = 'near(x[0], -0.5) || near(x[0], 0.5)' + self.bc = [ + df.DirichletBC(self.W.sub(0), df.Constant((1.0, 0.0)), top_bottom), + df.DirichletBC(self.W.sub(1), df.Constant(1.0), top_bottom), + ] + if not periodic: + self.bc += [ + df.DirichletBC(self.W.sub(0), self.u_ex, left_right), + df.DirichletBC(self.W.sub(1), self.p_ex, left_right), + ] + + # the residual is meaningless where the solution is prescribed, but only there: with + # periodicity the dofs on x = +-0.5 are unknowns and their residual has to be kept + dirichlet = top_bottom if periodic else 'on_boundary' + self.bc_hom = [ + df.DirichletBC(self.W.sub(0), df.Constant((0.0, 0.0)), dirichlet), + df.DirichletBC(self.W.sub(1), df.Constant(0.0), dirichlet), + ] + self.fix_bc_for_residual = True + + # residual form for a single node-to-node step, assembled once; `factor` and the + # boundary/forcing expressions carry the time dependence + self.factor = df.Constant(0.0) + self.w = df.Function(self.W) + u, p = df.split(self.w) + + F = df.dot(u, self.v) * df.dx + F += self.factor * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx + F += self.factor * self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx + F -= self.factor * df.dot(p, df.div(self.v)) * df.dx + F -= self.factor * df.dot(self.g, self.v) * df.dx + F -= self.factor * df.dot(df.div(u), self.q) * df.dx + + self.step = _NewtonStep(F, df.derivative(F, self.w), self.bc) + self.newton = df.NewtonSolver() + self.newton.parameters['absolute_tolerance'] = Sol_tol + self.newton.parameters['relative_tolerance'] = Sol_tol + self.newton.parameters['maximum_iterations'] = 20 + + def solve_system(self, rhs, factor, u0, t): + r""" + Newton solver for :math:`M w + factor \cdot N(w, t) = rhs`, where :math:`N` collects the + convective, viscous, pressure and divergence terms and ``rhs`` is the mass-weighted + right-hand side assembled by the sweeper. + + Parameters + ---------- + rhs : dtype_f + Right-hand side for the nonlinear system. + factor : float + Abbrev. for the node-to-node stepsize (or any other factor required). + u0 : dtype_u + Initial guess for the iterative solver. + t : float + Current time. + + Returns + ------- + w : dtype_u + Solution. + """ + self.factor.assign(factor) + self.u_ex.t = t + self.p_ex.t = t + self.g.t = t + + self.w.vector()[:] = u0.values.vector()[:] + self.step.rhs = rhs.values.vector() + self.newton.solve(self.step, self.w.vector()) + + me = self.dtype_u(self.W) + me.values.vector()[:] = self.w.vector()[:] + return me + + def eval_f(self, w, t): + r""" + Routine to evaluate the right-hand side of the problem in weak form, i.e. *without* + applying :math:`M^{-1}`. + + Parameters + ---------- + w : dtype_u + Current values of the numerical solution. + t : float + Current time at which the numerical solution is computed. + + Returns + ------- + f : dtype_f + The right-hand side. + """ + u, p = df.split(w.values) + self.g.t = t + + F = -df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx + F -= self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx + F += df.dot(p, df.div(self.v)) * df.dx + F += df.dot(self.g, self.v) * df.dx + F += df.dot(df.div(u), self.q) * df.dx + + f = self.dtype_f(self.W) + df.assemble(F, tensor=f.values.vector()) + return f + + def apply_mass_matrix(self, w): + r""" + Routine to apply the velocity mass matrix. + + Parameters + ---------- + w : dtype_u + Current values of the numerical solution. + + Returns + ------- + me : dtype_u + The product :math:`M \vec{w}`. + """ + me = self.dtype_u(self.W) + self.M.mult(w.values.vector(), me.values.vector()) + return me + + def u_exact(self, t): + r""" + Routine to compute the exact solution at time :math:`t`. + + Parameters + ---------- + t : float + Time of the exact solution. + + Returns + ------- + me : dtype_u + Exact solution. + """ + self.u_ex.t = t + self.p_ex.t = t + + me = self.dtype_u(self.W) + df.assign(me.values.sub(0), df.interpolate(self.u_ex, self.V)) + df.assign(me.values.sub(1), df.interpolate(self.p_ex, self.Q)) + return me + + def fix_residual(self, res): + """ + Applies homogeneous Dirichlet boundary conditions to the residual, on the Dirichlet part + of the boundary only. + + Parameters + ---------- + res : dtype_u + Residual. + """ + for bc in self.bc_hom: + bc.apply(res.values.vector()) + return None diff --git a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py new file mode 100644 index 0000000000..03dbee5d3f --- /dev/null +++ b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py @@ -0,0 +1,243 @@ +import dolfin as df +import numpy as np + +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import ( + fenics_NSE_2D_TaylorGreen, +) +from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import generic_implicit_mass + + +def setup(t0=0.0, dt=0.1, periodic=False, nelems=24, nu=0.1, num_nodes=3, maxiter=40, restol=1e-12): + """ + Helper routine to set up parameters + + Args: + t0: float, + initial time + dt: float, + time step size + periodic: bool, + use periodic instead of time-dependent Dirichlet conditions in x + nelems: int, + number of elements per spatial direction + nu: float, + kinematic viscosity + num_nodes: int, + number of collocation nodes + maxiter: int, + maximum number of SDC iterations + restol: float, + residual tolerance + + Returns: + description: dict, + pySDC description dictionary containing problem and method parameters. + controller_params: dict, + Parameters for the pySDC controller. + """ + # initialize level parameters + level_params = dict() + level_params['restol'] = restol + level_params['dt'] = dt + + # initialize step parameters + step_params = dict() + step_params['maxiter'] = maxiter + + # initialize sweeper parameters + sweeper_params = dict() + sweeper_params['quad_type'] = 'RADAU-RIGHT' + sweeper_params['num_nodes'] = num_nodes + sweeper_params['QI'] = 'LU' + + # initialize problem parameters + problem_params = dict() + problem_params['nelems'] = nelems + problem_params['t0'] = t0 + problem_params['order'] = 2 + problem_params['nu'] = nu + problem_params['periodic'] = periodic + problem_params['Sol_tol'] = 1e-13 + + # initialize controller parameters + controller_params = dict() + controller_params['logger_level'] = 30 + + # Fill description dictionary + description = dict() + description['problem_class'] = fenics_NSE_2D_TaylorGreen + description['sweeper_class'] = generic_implicit_mass + description['problem_params'] = problem_params + description['sweeper_params'] = sweeper_params + description['level_params'] = level_params + description['step_params'] = step_params + + return description, controller_params + + +def run_simulation(description, controller_params, Tend): + """ + Run the time integration for the 2D Taylor-Green Navier-Stokes benchmark. + + Args: + description: dict, + pySDC problem and method description. + controller_params: dict, + Parameters for the pySDC controller. + Tend: float, + Final simulation time. + + Returns: + P: problem instance, + Problem instance holding the function spaces and the exact solution. + stats: dict, + Collected runtime statistics. + uend: dtype_u, + Final solution at time Tend. + """ + t0 = description['problem_params']['t0'] + + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + P = controller.MS[0].levels[0].prob + uend, stats = controller.run(u0=P.u_exact(t0), t0=t0, Tend=Tend) + + return P, stats, uend + + +def relative_errors(u, uref): + """ + Relative L2 errors in velocity and pressure between two solutions on the same space. + + Args: + u: dtype_u, + Numerical solution. + uref: dtype_u, + Reference solution. + + Returns: + tuple of float: relative L2 error in velocity and in pressure. + """ + un, pn = u.values.split(deepcopy=True) + ur, pr = uref.values.split(deepcopy=True) + + return ( + df.errornorm(ur, un, 'L2') / df.norm(ur, 'L2'), + df.errornorm(pr, pn, 'L2') / df.norm(pr, 'L2'), + ) + + +def run_postprocessing(P, uend, Tend): + """ + Compute relative L2 errors between the numerical and the exact solution at the final time. + + Args: + P: problem instance, + Problem instance holding the exact solution. + uend: dtype_u, + Final solution at time Tend. + Tend: float, + Final simulation time. + + Returns: + tuple of float: relative L2 error in velocity and in pressure. + """ + return relative_errors(uend, P.u_exact(Tend)) + + +def order_study(dts, Tend, dt_ref=None, periodic=False, **kwargs): + r""" + Measure the observed temporal order of convergence. + + Errors are *not* taken against the exact solution: the spatial discretization error + dominates it for any affordable mesh, which hides the temporal order completely. Instead + two variants are offered, both of which cancel the spatial error exactly because every run + uses the same mesh: + + - ``dt_ref`` given: compare against a reference run with that much smaller step size, + - ``dt_ref`` omitted: compare consecutive step sizes with each other (Richardson). This + needs no reference run and is therefore a lot cheaper, at the cost of one order estimate. + + Args: + dts: list of float, + Step sizes to run, largest first, each one half of the previous. + Tend: float, + Final simulation time; must be an integer multiple of every step size. + dt_ref: float, + Step size for the reference run, or ``None`` to compare consecutive step sizes. + periodic: bool, + Use periodic instead of time-dependent Dirichlet conditions in x. + kwargs: + Passed on to :func:`setup`. + + Returns: + dts_out: list of float, + Step sizes the errors belong to; one shorter than ``dts`` without a reference. + errors_u: list of float, + Relative L2 velocity error per step size. + errors_p: list of float, + Relative L2 pressure error per step size. + """ + solutions = [] + for dt in dts: + description, controller_params = setup(dt=dt, periodic=periodic, **kwargs) + solutions.append(run_simulation(description, controller_params, Tend)[2]) + + if dt_ref is None: + pairs = list(zip(solutions[:-1], solutions[1:], strict=True)) + dts_out = dts[:-1] + else: + description, controller_params = setup(dt=dt_ref, periodic=periodic, **kwargs) + uref = run_simulation(description, controller_params, Tend)[2] + pairs = [(u, uref) for u in solutions] + dts_out = list(dts) + + errors = [relative_errors(u, ref) for u, ref in pairs] + + return dts_out, [e[0] for e in errors], [e[1] for e in errors] + + +def observed_order(dts, errors): + """ + Observed order of convergence between consecutive step sizes. + + Args: + dts: list of float, + Step sizes. + errors: list of float, + Corresponding errors. + + Returns: + list of float: observed orders, one shorter than the inputs. + """ + return [np.log(errors[i] / errors[i + 1]) / np.log(dts[i] / dts[i + 1]) for i in range(len(dts) - 1)] + + +def main(): + """ + Run the order study for both boundary condition variants and report the observed orders. + """ + Tend = 0.2 + dts = [0.2, 0.1, 0.05, 0.025] + + results = {} + for periodic in (True, False): + dts_out, errors_u, errors_p = order_study(dts, Tend, periodic=periodic) + results[periodic] = (dts_out, errors_u, errors_p) + + label = 'periodic' if periodic else 'time-dependent Dirichlet' + print(f'\n{label} boundary conditions in x:') + print(f'{"dt":>10} {"err(u)":>12} {"order(u)":>9} {"err(p)":>12} {"order(p)":>9}') + orders_u = [None] + observed_order(dts_out, errors_u) + orders_p = [None] + observed_order(dts_out, errors_p) + for dt, eu, ou, ep, op in zip(dts_out, errors_u, orders_u, errors_p, orders_p, strict=True): + su = ' --- ' if ou is None else f'{ou:9.2f}' + sp = ' --- ' if op is None else f'{op:9.2f}' + print(f'{dt:10.5f} {eu:12.4e} {su} {ep:12.4e} {sp}') + + return results + + +if __name__ == "__main__": + main() diff --git a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py new file mode 100644 index 0000000000..6f1f5b6dd1 --- /dev/null +++ b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py @@ -0,0 +1,181 @@ +import numpy as np +import pytest + + +@pytest.mark.fenics +def test_exact_solution_is_periodic(): + """ + The benchmark compares time-dependent Dirichlet against periodic conditions in x using the + *same* manufactured solution. That is only meaningful because the solution is genuinely + one-periodic in x, and constant on the top and bottom boundary. Check both, since every + conclusion drawn from this benchmark rests on it. + """ + import dolfin as df + from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import ( + fenics_NSE_2D_TaylorGreen, + ) + + prob = fenics_NSE_2D_TaylorGreen(nelems=8, t0=0.0, order=2, nu=0.05) + + for t in (0.0, 0.13, 0.4): + prob.u_ex.t = t + prob.p_ex.t = t + for y in np.linspace(-0.5, 0.5, 11): + left = np.hstack([prob.u_ex(-0.5, y), prob.p_ex(-0.5, y)]) + right = np.hstack([prob.u_ex(0.5, y), prob.p_ex(0.5, y)]) + assert np.allclose(left, right, atol=1e-12), f"solution is not periodic in x at t={t}, y={y}" + + for x in np.linspace(-0.5, 0.5, 11): + for y in (-0.5, 0.5): + value = np.hstack([prob.u_ex(x, y), prob.p_ex(x, y)]) + assert np.allclose(value, [1.0, 0.0, 1.0], atol=1e-12), f"top/bottom data varies at t={t}" + + # the periodic function space must not lose anything when representing this solution + prob_periodic = fenics_NSE_2D_TaylorGreen(nelems=8, t0=0.0, order=2, nu=0.05, periodic=True) + assert prob_periodic.W.dim() < prob.W.dim(), "periodic space should have fewer dofs" + + u, p = prob.u_exact(0.3).values.split(deepcopy=True) + up, pp = prob_periodic.u_exact(0.3).values.split(deepcopy=True) + assert abs(df.norm(u, 'L2') - df.norm(up, 'L2')) < 1e-12 + assert abs(df.norm(p, 'L2') - df.norm(pp, 'L2')) < 1e-12 + + +@pytest.mark.fenics +def test_eval_f(): + """ + Anchor ``eval_f`` on the analytical solution: the semi-discrete system is M u' = f(u, t), + so evaluating f at the exact solution must reproduce M du/dt on the interior dofs, and do + so with the second order expected from the interpolation of the data. + + Compared against du/dt rather than against ``solve_system`` on purpose -- a sign error + shared by both would pass a consistency check between them. + """ + import dolfin as df + from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import ( + fenics_NSE_2D_TaylorGreen, + ) + + t, nu = 0.3, 0.05 + errors = [] + for nelems in (16, 32): + prob = fenics_NSE_2D_TaylorGreen(nelems=nelems, t0=0.0, order=2, nu=nu) + + dudt = df.Expression( + ( + '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])' + ' + 2*pi*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])', + '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])' + ' - 2*pi*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])', + ), + pi=np.pi, + nu=nu, + t=t, + degree=prob.order + 2, + ) + + ut = prob.dtype_u(prob.W) + df.assign(ut.values.sub(0), df.interpolate(dudt, prob.V)) + expected = prob.apply_mass_matrix(ut) + f = prob.eval_f(prob.u_exact(t), t) + + # eval_f integrates by parts, so the two can only agree away from the boundary + prob.fix_residual(f) + prob.fix_residual(expected) + + velocity_dofs = np.array(prob.W.sub(0).dofmap().dofs()) + a = f.values.vector()[velocity_dofs] + b = expected.values.vector()[velocity_dofs] + errors.append(np.linalg.norm(a - b) / np.linalg.norm(b)) + + assert errors[0] < 2e-2, f"eval_f does not match M du/dt: relative error {errors[0]:.3e}" + order = np.log2(errors[0] / errors[1]) + assert order > 1.7, f"eval_f converges at order {order:.2f}, expected second order" + + +@pytest.mark.fenics +def test_order_reduction(): + """ + The point of the whole benchmark: the same exact solution, computed with periodic + conditions in x, reaches the design order 2M-1 = 5 of RADAU-RIGHT with M = 3, while with + time-dependent Dirichlet conditions in x it does not. + + What is asserted here is the robust part of that. At a mesh resolution CI can afford, the + difference in observed *order* is modest (roughly 4.7 against 4.3 in the pressure, the + reduction Radau IIA is known for, 2M-1 down to M+1), and too small a margin to assert on. + The difference in the error *constant* is not: the time-dependent boundary data costs + close to an order of magnitude in the pressure at every step size tested. Refining the + mesh deepens both effects, because the reduction is driven by stiffness -- see the + docstring of ``order_study`` and the numbers printed by running the script directly. + """ + from pySDC.projects.StroemungsRaum.run_Navier_Stokes_TaylorGreen_FEniCS import ( + order_study, + observed_order, + ) + + Tend, dts = 0.2, [0.2, 0.1, 0.05] + errors, orders = {}, {} + for periodic in (True, False): + dts_out, errors_u, errors_p = order_study(dts, Tend, periodic=periodic) + errors[periodic] = (errors_u, errors_p) + orders[periodic] = (observed_order(dts_out, errors_u)[0], observed_order(dts_out, errors_p)[0]) + + # with periodic conditions there is no time-dependent boundary data and the method + # attains (close to) its design order + assert orders[True][0] > 4.5, f"periodic velocity order {orders[True][0]:.2f} below design order" + + # time-dependent Dirichlet data costs roughly an order of magnitude in the pressure + for i, dt in enumerate(dts[:-1]): + ratio = errors[False][1][i] / errors[True][1][i] + assert ratio > 3.0, f"pressure error ratio at dt={dt} is only {ratio:.1f}, expected a clear gap" + + # and it does not reach the order the periodic variant does + assert orders[False][1] < orders[True][1], ( + f"pressure order with Dirichlet data ({orders[False][1]:.2f}) is not below " + f"the periodic one ({orders[True][1]:.2f})" + ) + + +@pytest.mark.fenics +@pytest.mark.parametrize("periodic", [False, True]) +def test_run_benchmark(periodic): + """ + End-to-end smoke test: a couple of SDC steps have to land close to the exact solution. + The tolerance is set by the spatial discretization error on this coarse mesh, not by the + time integration, so it says nothing about the temporal order -- see test_order_reduction. + """ + from pySDC.projects.StroemungsRaum.run_Navier_Stokes_TaylorGreen_FEniCS import ( + setup, + run_simulation, + run_postprocessing, + ) + + description, controller_params = setup(dt=0.05, periodic=periodic, nelems=16, nu=0.05, num_nodes=2) + P, _, uend = run_simulation(description, controller_params, Tend=0.1) + error_u, error_p = run_postprocessing(P, uend, Tend=0.1) + + assert error_u < 1e-3, f"relative velocity error {error_u:.3e} exceeds tolerance" + assert error_p < 1e-2, f"relative pressure error {error_p:.3e} exceeds tolerance" + + +@pytest.mark.fenics +@pytest.mark.parametrize("periodic", [False, True]) +def test_solve_system(periodic): + """ + ``solve_system`` solves M w - factor * f(w) = rhs, so feeding it the right-hand side built + from the exact solution must return the exact solution. Newton is started away from the + answer so this actually exercises the solve, the Jacobian and the boundary conditions. + """ + from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import ( + fenics_NSE_2D_TaylorGreen, + ) + + t, factor = 0.3, 0.01 + prob = fenics_NSE_2D_TaylorGreen(nelems=16, t0=0.0, order=2, nu=0.05, periodic=periodic) + + uex = prob.u_exact(t) + rhs = prob.apply_mass_matrix(uex) - factor * prob.eval_f(uex, t) + + w = prob.solve_system(rhs, factor, prob.dtype_u(prob.W), t) + + rel_err = abs(w - uex) / abs(uex) + assert rel_err < 1e-9, f"solve_system did not recover the exact solution: {rel_err:.3e}" From bd05ff4fe0142c948d68575a1554d6e204928640 Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Thu, 10 Sep 2026 18:07:49 +0200 Subject: [PATCH 2/7] Remove the mass matrix inversion from the monolithic NSE problem solve_system applied Mf^-1 to the right-hand side handed over by the sweeper, only for the residual form to apply Mf straight back to it via the dot(rhs, v)*dx terms. The round trip is exact, so it bought nothing while costing an LU solve per node per sweep and capping attainable accuracy at the tolerance of that solve. It was there because dolfin's solve(F == 0, ...) accepts a form and not a vector. Subtract the right-hand side vector in a NonlinearProblem instead. The same pattern is needed by the Taylor-Green benchmark, so it lives in newton_step.py and both problem classes use it. Building the Newton problem this way also means the form, its Jacobian and the solver are set up once, with the step size carried by a Constant, rather than rebuilt on every solve_system call. Previously each distinct value of `factor` was baked into the form as a literal, giving a different form signature per step size (verified) and so a separate compilation. FFC caches those across runs, so this shows up on a cold cache or when new step sizes keep appearing, not in steady-state runtime: on a warm cache the drag/lift benchmark test is unchanged at 2.2s against 2.1s before. Numerically neutral: drag and lift in that benchmark agree with the previous implementation to 13 significant digits (relative difference 9e-15 and 2e-12), which is the round-off of the LU solve that has been removed. Co-Authored-By: Claude Opus 5 --- ...Stokes_2D_TaylorGreen_monolithic_FEniCS.py | 43 +--------- .../NavierStokes_2D_monolithic_FEniCS.py | 81 +++++++------------ .../problem_classes/newton_step.py | 50 ++++++++++++ 3 files changed, 84 insertions(+), 90 deletions(-) create mode 100644 pySDC/projects/StroemungsRaum/problem_classes/newton_step.py diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py index 62160ad78f..e9f568625e 100644 --- a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py @@ -5,6 +5,7 @@ from pySDC.core.problem import Problem from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh +from pySDC.projects.StroemungsRaum.problem_classes.newton_step import NewtonStep class _PeriodicX(df.SubDomain): @@ -23,45 +24,6 @@ def map(self, x, y): y[1] = x[1] -class _NewtonStep(df.NonlinearProblem): - r""" - Newton problem for the node-to-node step :math:`M w + \Delta t_{QI} N(w) = rhs`. - - ``rhs`` is kept as an assembled vector and subtracted from the residual here, rather - than being written into the variational form as :math:`\int_\Omega rhs \cdot v\,dx`. - That form would apply the mass matrix to it, which would then have to be undone by a - mass matrix solve beforehand -- an exact round trip that costs a solve per node per - sweep and caps the attainable accuracy at the tolerance of that solve. - - Parameters - ---------- - F : UFL form - Residual form, without the right-hand side term. - J : UFL form - Jacobian of ``F``. - bcs : list of DirichletBC - Boundary conditions, applied in residual form (``bc.apply(b, x)``). - """ - - def __init__(self, F, J, bcs): - super().__init__() - self.F_form = F - self.J_form = J - self.bcs = bcs - self.rhs = None - - def F(self, b, x): - df.assemble(self.F_form, tensor=b) - b.axpy(-1.0, self.rhs) - for bc in self.bcs: - bc.apply(b, x) - - def J(self, A, x): - df.assemble(self.J_form, tensor=A) - for bc in self.bcs: - bc.apply(A) - - class fenics_NSE_2D_TaylorGreen(Problem): r""" Forced two-dimensional incompressible Navier-Stokes equations on :math:`\Omega = [-0.5, 0.5]^2`, @@ -242,7 +204,7 @@ def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, Sol_tol= F -= self.factor * df.dot(self.g, self.v) * df.dx F -= self.factor * df.dot(df.div(u), self.q) * df.dx - self.step = _NewtonStep(F, df.derivative(F, self.w), self.bc) + self.step = NewtonStep(F, df.derivative(F, self.w)) self.newton = df.NewtonSolver() self.newton.parameters['absolute_tolerance'] = Sol_tol self.newton.parameters['relative_tolerance'] = Sol_tol @@ -277,6 +239,7 @@ def solve_system(self, rhs, factor, u0, t): self.w.vector()[:] = u0.values.vector()[:] self.step.rhs = rhs.values.vector() + self.step.bcs = self.bc self.newton.solve(self.step, self.w.vector()) me = self.dtype_u(self.W) diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py index eb34279281..c81289115b 100755 --- a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py @@ -6,6 +6,7 @@ from pySDC.core.problem import Problem from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh +from pySDC.projects.StroemungsRaum.problem_classes.newton_step import NewtonStep class fenics_NSE_2D_Monolithic(Problem): @@ -58,8 +59,6 @@ class fenics_NSE_2D_Monolithic(Problem): Defines the mixed function space for the coupled velocity-pressure system. M : scalar, vector, matrix or higher rank tensor Denotes the expression :math:`\int_\Omega u_t v\,dx`. - Mf : scalar, vector, matrix or higher rank tensor - Denotes the expression :math:`\int_\Omega u v\,dx + \int_\Omega p q\,dx`. g : Expression The forcing term :math:`f` in the Navier-Stokes momentum equation. bc : DirichletBC @@ -120,10 +119,6 @@ def __init__(self, t0=0.0, order=2, nu=0.001, Sol_tol=1e-10): a_M = df.inner(self.u, self.v) * df.dx self.M = df.assemble(a_M) - # full mass matrix - a_Mf = df.inner(self.u, self.v) * df.dx + df.inner(self.p, self.q) * df.dx - Mf = df.assemble(a_Mf) - # define the time-dependent inflow profile as an Expression Uin = '4.0*1.5*sin(pi*t/8)*x[1]*(0.41 - x[1]) / pow(0.41, 2)' self.u_in = df.Expression((Uin, '0'), pi=np.pi, t=t0, degree=self.order) @@ -160,9 +155,25 @@ def __init__(self, t0=0.0, order=2, nu=0.001, Sol_tol=1e-10): self.xdmffile_p = None self.xdmffile_u = None - # set up linear solver for the inversion of the mass matrix - self.solver = df.LUSolver(Mf) - # self.solver.parameters['reuse_factorization'] = True + # interpolated forcing term referenced by the residual form + self.g_h = df.Function(self.V) + + # residual form for a single node-to-node step, assembled once; the step size is a + # Constant and the boundary and forcing data carry the time dependence + self.factor = df.Constant(0.0) + self.w = df.Function(self.W) + u, p = df.split(self.w) + + F = df.dot(u, self.v) * df.dx + F += self.factor * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx + F += self.factor * self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx + F -= self.factor * df.dot(p, df.div(self.v)) * df.dx + F -= self.factor * df.dot(self.g_h, self.v) * df.dx + F -= self.factor * df.dot(df.div(u), self.q) * df.dx + + self.step = NewtonStep(F, df.derivative(F, self.w)) + self.newton = df.NewtonSolver() + self.newton.parameters['absolute_tolerance'] = Sol_tol def solve_system(self, rhs, factor, u0, t): r""" @@ -186,33 +197,22 @@ def solve_system(self, rhs, factor, u0, t): w : dtype_u Solution. """ - # introduce the coupled solution vector for velocity and pressure - w = self.dtype_u(u0) - u, p = df.split(w.values) - - # get the SDC right-hand side - rhs = self.__invert_mass_matrix(rhs) - rhs_u, rhs_p = df.split(rhs.values) - - # update time in boundary conditions + # update time in boundary conditions and in the forcing term self.u_in.t = t - - # get the forcing term self.g.t = t - g = df.interpolate(self.g, self.V) + self.g_h.interpolate(self.g) + self.factor.assign(factor) - # build the variational form for the coupled system - F = df.dot(u, self.v) * df.dx - F += factor * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx - F += factor * self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx - F -= factor * df.dot(p, df.div(self.v)) * df.dx - F -= factor * df.dot(g, self.v) * df.dx - F -= factor * df.dot(df.div(u), self.q) * df.dx - F -= df.dot(rhs_u, self.v) * df.dx - F -= df.dot(rhs_p, self.q) * df.dx + # the SDC right-hand side enters the residual as a vector, no mass matrix involved + self.w.vector()[:] = u0.values.vector()[:] + self.step.rhs = rhs.values.vector() + self.step.bcs = self.bc # solve the nonlinear system using Newton's method - df.solve(F == 0, w.values, self.bc, solver_parameters={"newton_solver": {"absolute_tolerance": self.Sol_tol}}) + self.newton.solve(self.step, self.w.vector()) + + w = self.dtype_u(self.W) + w.values.vector()[:] = self.w.vector()[:] return w @@ -271,25 +271,6 @@ def apply_mass_matrix(self, w): return me - def __invert_mass_matrix(self, w): - r""" - Helper routine to invert the full mass matrix Mf. - - Parameters - ---------- - w : dtype_u - Current values of the numerical solution. - - Returns - ------- - me : dtype_u - The product :math:`Mf^{-1} \vec{w}`. - """ - - me = self.dtype_u(self.W) - self.solver.solve(me.values.vector(), w.values.vector()) - return me - def u_exact(self, t): r""" Routine to compute the exact solution at time :math:`t`. diff --git a/pySDC/projects/StroemungsRaum/problem_classes/newton_step.py b/pySDC/projects/StroemungsRaum/problem_classes/newton_step.py new file mode 100644 index 0000000000..1722516f68 --- /dev/null +++ b/pySDC/projects/StroemungsRaum/problem_classes/newton_step.py @@ -0,0 +1,50 @@ +import dolfin as df + + +class NewtonStep(df.NonlinearProblem): + r""" + Newton problem for a single SDC node-to-node step :math:`M w + \Delta t_{QI} N(w) = rhs`. + + The right-hand side handed over by the sweeper is an assembled vector, and it is subtracted + from the residual here, at the algebraic level. The alternative -- writing it into the + variational form as :math:`\int_\Omega rhs \cdot v\,dx` so that dolfin's high level + ``solve(F == 0, ...)`` interface can be used -- applies the mass matrix to it, which then + has to be undone by a mass matrix solve beforehand. That round trip is exact in exact + arithmetic, so it buys nothing, while costing a solve per node per sweep and capping the + attainable accuracy at the tolerance of that solve. + + ``rhs`` and ``bcs`` are set per solve rather than at construction, so that the form and its + Jacobian can be compiled once and reused with the step size carried by a ``Constant``. + + Parameters + ---------- + F : UFL form + Residual form of the step, *without* the right-hand side term. + J : UFL form + Jacobian of ``F``. + + Attributes + ---------- + rhs : GenericVector + Right-hand side vector for the current solve, subtracted from the residual. + bcs : list of DirichletBC + Boundary conditions for the current solve, applied in residual form. + """ + + def __init__(self, F, J): + super().__init__() + self.F_form = F + self.J_form = J + self.rhs = None + self.bcs = [] + + def F(self, b, x): + df.assemble(self.F_form, tensor=b) + b.axpy(-1.0, self.rhs) + for bc in self.bcs: + bc.apply(b, x) + + def J(self, A, x): + df.assemble(self.J_form, tensor=A) + for bc in self.bcs: + bc.apply(A) From daa5fda021f1aed0b9c1a7ddae0ab41b6deb2dca Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Thu, 10 Sep 2026 18:53:16 +0200 Subject: [PATCH 3/7] Make the order reduction visible: use four collocation nodes RADAU-RIGHT with M nodes has design order 2M-1 and falls back to the stiff order M+1 when the boundary data is time dependent, so the gap this benchmark can show is M-2. That is *identically zero for M = 2*, which is what the original FEniCSx version ran with: that configuration cannot exhibit the phenomenon it was written to demonstrate, whatever else is fixed. Measured at nelems=24, nu=0.1, Tend=0.2, orders from consecutive step sizes: M design stiff periodic order(p) Dirichlet order(p) 2 3 3 2.85 2.80 3 5 4 4.69 4.28 4 7 5 7.35 4.96 M = 4 is the cheapest setting that separates the two unmistakably: order 7 against 5, with the pressure error 52x larger at the finest step size. The study script and the order test now use it, and the rule is recorded in the problem class docstring and the project README so the M = 2 trap is not stepped in again. The test asserts the gap and the error ratio rather than absolute orders, and still runs in about a minute. Co-Authored-By: Claude Opus 5 --- pySDC/projects/StroemungsRaum/README.rst | 6 ++ ...Stokes_2D_TaylorGreen_monolithic_FEniCS.py | 6 ++ .../run_Navier_Stokes_TaylorGreen_FEniCS.py | 10 +++- .../test_Navier_Stokes_TaylorGreen_FEniCS.py | 59 ++++++++++--------- 4 files changed, 50 insertions(+), 31 deletions(-) diff --git a/pySDC/projects/StroemungsRaum/README.rst b/pySDC/projects/StroemungsRaum/README.rst index 33054ea1ab..cf61d58488 100644 --- a/pySDC/projects/StroemungsRaum/README.rst +++ b/pySDC/projects/StroemungsRaum/README.rst @@ -47,6 +47,12 @@ computed with time-dependent Dirichlet conditions in :math:`x` or with periodic ones, and the difference in the observed temporal order isolates the order reduction caused by the time-dependent boundary data alone. +The number of collocation nodes decides whether the effect is visible: RADAU-RIGHT +with :math:`M` nodes drops from its design order :math:`2M-1` to the stiff order +:math:`M+1`, so the gap is :math:`M-2` and vanishes for :math:`M = 2`. With +:math:`M = 4` the measured pressure orders are 7 with periodic and 5 with +time-dependent Dirichlet conditions. + Funding ------- Funded by the **German Federal Ministry of Education and Research (BMBF)** under diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py index e9f568625e..3c9cd2c697 100644 --- a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py @@ -53,6 +53,12 @@ class fenics_NSE_2D_TaylorGreen(Problem): The only difference between the two runs is therefore the presence of time-dependent boundary data, which is what isolates the order reduction. + Note that the number of collocation nodes decides whether anything can be seen at all. + RADAU-RIGHT with :math:`M` nodes has design order :math:`2M-1` and falls back to the stiff + order :math:`M+1` in the presence of time-dependent boundary data, so the gap on offer is + :math:`M-2`: **zero for M = 2**, where both are 3. Use :math:`M \geq 4`; at :math:`M = 4` + the measured pressure orders are 7 (periodic) against 5 (Dirichlet). + The problem is discretized in space with Taylor-Hood elements on a mixed velocity-pressure space and solved monolithically, so the semi-discrete system is the differential-algebraic system :math:`M \dot{w} = f(w, t)` with the singular mass matrix :math:`M = \mathrm{diag}(M_v, 0)`. diff --git a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py index 03dbee5d3f..adc9d9cbdc 100644 --- a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py @@ -8,7 +8,7 @@ from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import generic_implicit_mass -def setup(t0=0.0, dt=0.1, periodic=False, nelems=24, nu=0.1, num_nodes=3, maxiter=40, restol=1e-12): +def setup(t0=0.0, dt=0.1, periodic=False, nelems=24, nu=0.1, num_nodes=4, maxiter=40, restol=1e-12): """ Helper routine to set up parameters @@ -215,8 +215,14 @@ def observed_order(dts, errors): def main(): - """ + r""" Run the order study for both boundary condition variants and report the observed orders. + + RADAU-RIGHT with M nodes has design order :math:`2M-1` and, on a stiff problem with + time-dependent boundary data, drops to the stiff order :math:`M+1`. The gap the benchmark + can show is therefore :math:`M-2`, and **nothing at all is visible for M = 2**, where the + two coincide at 3. M = 4 is used here because it is the cheapest setting that makes the + reduction unmistakable: order 7 against 5 in the pressure. """ Tend = 0.2 dts = [0.2, 0.1, 0.05, 0.025] diff --git a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py index 6f1f5b6dd1..2a8ac5a2ae 100644 --- a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py @@ -94,45 +94,46 @@ def test_eval_f(): @pytest.mark.fenics def test_order_reduction(): - """ + r""" The point of the whole benchmark: the same exact solution, computed with periodic - conditions in x, reaches the design order 2M-1 = 5 of RADAU-RIGHT with M = 3, while with - time-dependent Dirichlet conditions in x it does not. - - What is asserted here is the robust part of that. At a mesh resolution CI can afford, the - difference in observed *order* is modest (roughly 4.7 against 4.3 in the pressure, the - reduction Radau IIA is known for, 2M-1 down to M+1), and too small a margin to assert on. - The difference in the error *constant* is not: the time-dependent boundary data costs - close to an order of magnitude in the pressure at every step size tested. Refining the - mesh deepens both effects, because the reduction is driven by stiffness -- see the - docstring of ``order_study`` and the numbers printed by running the script directly. + conditions in x, reaches the design order 2M-1 of RADAU-RIGHT, while with time-dependent + Dirichlet conditions in x it drops to the stiff order M+1. + + M = 4 is deliberate. The gap the benchmark can show is (2M-1) - (M+1) = M-2, so it is + *identically zero for M = 2*, where both orders are 3 -- a setup with two nodes cannot + exhibit this phenomenon no matter what else is done. At M = 4 the two orders are 7 and 5 + and the separation is unmistakable. """ from pySDC.projects.StroemungsRaum.run_Navier_Stokes_TaylorGreen_FEniCS import ( order_study, observed_order, ) - Tend, dts = 0.2, [0.2, 0.1, 0.05] + Tend, dts, num_nodes = 0.2, [0.2, 0.1, 0.05, 0.025], 4 errors, orders = {}, {} for periodic in (True, False): - dts_out, errors_u, errors_p = order_study(dts, Tend, periodic=periodic) + dts_out, errors_u, errors_p = order_study(dts, Tend, periodic=periodic, num_nodes=num_nodes) errors[periodic] = (errors_u, errors_p) - orders[periodic] = (observed_order(dts_out, errors_u)[0], observed_order(dts_out, errors_p)[0]) - - # with periodic conditions there is no time-dependent boundary data and the method - # attains (close to) its design order - assert orders[True][0] > 4.5, f"periodic velocity order {orders[True][0]:.2f} below design order" - - # time-dependent Dirichlet data costs roughly an order of magnitude in the pressure - for i, dt in enumerate(dts[:-1]): - ratio = errors[False][1][i] / errors[True][1][i] - assert ratio > 3.0, f"pressure error ratio at dt={dt} is only {ratio:.1f}, expected a clear gap" - - # and it does not reach the order the periodic variant does - assert orders[False][1] < orders[True][1], ( - f"pressure order with Dirichlet data ({orders[False][1]:.2f}) is not below " - f"the periodic one ({orders[True][1]:.2f})" - ) + # the finest pair is the most asymptotic estimate + orders[periodic] = (observed_order(dts_out, errors_u)[-1], observed_order(dts_out, errors_p)[-1]) + + design, stiff = 2 * num_nodes - 1, num_nodes + 1 + + # without time-dependent boundary data the method attains its design order + assert ( + orders[True][0] > design - 1.0 + ), f"periodic velocity order {orders[True][0]:.2f} is not close to the design order {design}" + + # with it, the pressure drops towards the stiff order and stays well clear of the design one + assert ( + orders[False][1] < (design + stiff) / 2 + ), f"pressure order with Dirichlet data is {orders[False][1]:.2f}, expected near {stiff}" + gap = orders[True][1] - orders[False][1] + assert gap > 1.0, f"pressure order gap is only {gap:.2f}, expected close to {design - stiff}" + + # and the accumulated error differs by more than an order of magnitude at the finest step + ratio = errors[False][1][-1] / errors[True][1][-1] + assert ratio > 10.0, f"pressure error ratio at the finest step size is only {ratio:.1f}" @pytest.mark.fenics From 6db3c8309d61f966f4c6c3bbc4ef9bcb5fc2682e Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Thu, 10 Sep 2026 18:58:51 +0200 Subject: [PATCH 4/7] Document that the prescribed boundary pressure acts as a partial lifting Removing the time-dependent pressure condition on x = +-0.5 drops the observed pressure order from M+1 to M and grows the error by ~24x (measured at M=4: 4.96 -> 4.11). Prescribing it therefore supplies constraint information rather than being neutral, which places this benchmark at the "constraint lifting" rung of PR #641's ladder rather than the plain algebraic one, and means the 7-vs-5 gap understates a setup without prescribed boundary pressure. Co-Authored-By: Claude Opus 5 --- .../NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py index 3c9cd2c697..2b9aed8c42 100644 --- a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py @@ -53,6 +53,12 @@ class fenics_NSE_2D_TaylorGreen(Problem): The only difference between the two runs is therefore the presence of time-dependent boundary data, which is what isolates the order reduction. + On :math:`x = \pm 0.5` the *pressure* is prescribed from the exact solution as well. That is + not neutral: it acts as a partial lifting of the algebraic constraint, and it lifts the + observed pressure order from :math:`M` to :math:`M+1`. Dropping it gives order :math:`M` + and a roughly 20 times larger error, so the gap measured here understates what a setup + without prescribed boundary pressure would show. + Note that the number of collocation nodes decides whether anything can be seen at all. RADAU-RIGHT with :math:`M` nodes has design order :math:`2M-1` and falls back to the stiff order :math:`M+1` in the presence of time-dependent boundary data, so the gap on offer is From 2cd08a7b77d03599005304508a32341cd73a260f Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Thu, 10 Sep 2026 19:43:30 +0200 Subject: [PATCH 5/7] Add differentiated boundary conditions to recover the lost order Ships the remedy prototyped against the order reduction the benchmark measures. Instead of evaluating the time-dependent boundary data pointwise at the node, u_B(tau_m) = g(tau_m), it is imposed on the derivative and the stage value is recovered by the collocation quadrature, u_B(tau_m) = g(t0) + dt * sum_j Q[m,j] * gdot(tau_j). The two differ by the quadrature error O(dt^(M+1)), but the second is consistent with the collocation polynomial rather than pointwise exact. This is the boundary-condition analogue of the differentiated-constraint remedy explored for a time-dependent constraint in #641. Measured at M=4, nelems=24, nu=0.1, pressure order and error at dt=0.1: periodic (best possible) 6.32 7.4e-08 pointwise 5.74 9.7e-07 differentiated 6.30 1.3e-07 so it removes most of the penalty, leaving a factor of 1.8 against the periodic case that is a constant rather than a rate. No new sweeper was needed beyond a hook: generic_implicit_mass is already the y-formulation, so generic_implicit_mass_diffbc only hands the problem the collocation data of the step, which a problem class cannot see on its own. The self-consistency trap #641 warns about does not apply here, because eval_f never sees the boundary condition -- fix_residual zeroes those rows. The observed orders are pre-asymptotic: the periodic reference does not reach its design order 7 either, and on finer step sizes every variant collapses against a solver floor near 1e-10. The tests therefore assert the error rather than the order, and the docstrings say so; separating 2M-1 from 2M-2 needs a better conditioned testbed than a 2D nonlinear NSE benchmark. Verified that the differentiated variant solves the same problem: its solution converges to the pointwise one at O(dt^5) = O(dt^(M+1)), the size of the perturbation. Basing the quadrature on the exact g(t0) or chaining it from the incoming numerical value gives results identical to three digits. Co-Authored-By: Claude Opus 5 --- pySDC/projects/StroemungsRaum/README.rst | 7 + ...Stokes_2D_TaylorGreen_monolithic_FEniCS.py | 148 +++++++++++++++++- .../run_Navier_Stokes_TaylorGreen_FEniCS.py | 36 ++++- .../sweepers/generic_implicit_mass.py | 26 +++ .../test_Navier_Stokes_TaylorGreen_FEniCS.py | 55 +++++++ 5 files changed, 259 insertions(+), 13 deletions(-) diff --git a/pySDC/projects/StroemungsRaum/README.rst b/pySDC/projects/StroemungsRaum/README.rst index cf61d58488..340c2c77ed 100644 --- a/pySDC/projects/StroemungsRaum/README.rst +++ b/pySDC/projects/StroemungsRaum/README.rst @@ -53,6 +53,13 @@ with :math:`M` nodes drops from its design order :math:`2M-1` to the stiff order :math:`M = 4` the measured pressure orders are 7 with periodic and 5 with time-dependent Dirichlet conditions. +The third variant, ``differentiated_bc``, imposes the boundary data on its time +derivative and recovers the stage values by collocation quadrature instead of +evaluating the data pointwise at the nodes. This is the boundary-condition +analogue of the differentiated-constraint remedy explored in pull request #641, +and it removes most of the penalty: at :math:`M = 4` the pressure error drops by +roughly an order of magnitude, to within a small factor of the periodic case. + Funding ------- Funded by the **German Federal Ministry of Education and Research (BMBF)** under diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py index 2b9aed8c42..a03c082025 100644 --- a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py @@ -59,6 +59,11 @@ class fenics_NSE_2D_TaylorGreen(Problem): and a roughly 20 times larger error, so the gap measured here understates what a setup without prescribed boundary pressure would show. + Setting ``differentiated_bc`` imposes the time-dependent data in differentiated form and + recovers most of the lost order, following the remedy explored for a time-dependent + *constraint* in pull request #641. It requires the ``generic_implicit_mass_diffbc`` sweeper. + See :meth:`prepare_step` for the construction and its measured effect. + Note that the number of collocation nodes decides whether anything can be seen at all. RADAU-RIGHT with :math:`M` nodes has design order :math:`2M-1` and falls back to the stiff order :math:`M+1` in the presence of time-dependent boundary data, so the gap on offer is @@ -83,6 +88,9 @@ class fenics_NSE_2D_TaylorGreen(Problem): Kinematic viscosity :math:`\nu`. periodic : bool, optional Use periodic instead of time-dependent Dirichlet conditions on :math:`x = \pm 0.5`. + differentiated_bc : bool, optional + Impose the time-dependent boundary data in differentiated form; needs ``periodic=False`` + and the ``generic_implicit_mass_diffbc`` sweeper. Sol_tol : float, optional Absolute tolerance for the Newton solver. @@ -116,7 +124,7 @@ class fenics_NSE_2D_TaylorGreen(Problem): df.set_log_active(False) - def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, Sol_tol=1e-10): + def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, differentiated_bc=False, Sol_tol=1e-10): # set logger level for FFC and dolfin logging.getLogger('FFC').setLevel(logging.WARNING) @@ -138,7 +146,15 @@ def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, Sol_tol= super().__init__(self.W) self._makeAttributeAndRegister( - 'nelems', 't0', 'order', 'nu', 'periodic', 'Sol_tol', localVars=locals(), readOnly=True + 'nelems', + 't0', + 'order', + 'nu', + 'periodic', + 'differentiated_bc', + 'Sol_tol', + localVars=locals(), + readOnly=True, ) self.logger.debug('DoFs on this level: %d', self.W.dim()) @@ -183,17 +199,26 @@ def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, Sol_tol= # on y = +-0.5 the exact solution is constant in space and time top_bottom = 'near(x[1], -0.5) || near(x[1], 0.5)' - left_right = 'near(x[0], -0.5) || near(x[0], 0.5)' - self.bc = [ + self.left_right = 'near(x[0], -0.5) || near(x[0], 0.5)' + self.bc_fixed = [ df.DirichletBC(self.W.sub(0), df.Constant((1.0, 0.0)), top_bottom), df.DirichletBC(self.W.sub(1), df.Constant(1.0), top_bottom), ] + self.bc = list(self.bc_fixed) if not periodic: self.bc += [ - df.DirichletBC(self.W.sub(0), self.u_ex, left_right), - df.DirichletBC(self.W.sub(1), self.p_ex, left_right), + df.DirichletBC(self.W.sub(0), self.u_ex, self.left_right), + df.DirichletBC(self.W.sub(1), self.p_ex, self.left_right), ] + # boundary conditions per collocation node, filled in by prepare_step + self._node_times = None + self._node_bcs = None + if differentiated_bc: + if periodic: + raise ValueError('differentiated_bc has no effect without time-dependent boundary data') + self.u_dot, self.p_dot = self._boundary_derivatives(nu, order, t0) + # the residual is meaningless where the solution is prescribed, but only there: with # periodicity the dofs on x = +-0.5 are unknowns and their residual has to be kept dirichlet = top_bottom if periodic else 'on_boundary' @@ -222,6 +247,109 @@ def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, Sol_tol= self.newton.parameters['relative_tolerance'] = Sol_tol self.newton.parameters['maximum_iterations'] = 20 + @staticmethod + def _boundary_derivatives(nu, order, t0): + r""" + Time derivatives of the boundary data, needed to impose it in differentiated form. + + Returns + ------- + u_dot, p_dot : Expression + :math:`\partial_t u` and :math:`\partial_t p` of the manufactured solution. + """ + kwargs = dict(pi=np.pi, nu=nu, t=t0, degree=order + 2) + u_dot = df.Expression( + ( + '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])' + ' + 2*pi*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])', + '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])' + ' - 2*pi*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])', + ), + **kwargs, + ) + p_dot = df.Expression( + '(4.0/17.0)*cos(pi*x[1])*(' + '-16*pi*pi*nu*exp(-16*pi*pi*nu*t)*cos(4*pi*(x[0] - t))' + ' + 4*pi*exp(-16*pi*pi*nu*t)*sin(4*pi*(x[0] - t)))', + **kwargs, + ) + return u_dot, p_dot + + def prepare_step(self, t0, dt, coll): + r""" + Build the differentiated boundary conditions for every collocation node of a step. + + Rather than evaluating the boundary data pointwise at the node, :math:`u_B(\tau_m) = + g(\tau_m)`, the condition is imposed on the *derivative* and the stage value recovered + by the collocation quadrature, + + .. math:: + u_B(\tau_m) = g(t_0) + \Delta t \sum_j Q_{mj}\, \dot{g}(\tau_j). + + The two differ by the quadrature error :math:`O(\Delta t^{M+1})`, but the second is + consistent with the collocation polynomial instead of pointwise exact, which is what + recovers the order lost to time-dependent boundary data. + + Measured at :math:`M = 4`, ``nelems=24``, ``nu=0.1``, orders and errors in the pressure + taken from consecutive step sizes: + + ========================= ========== ===================== + boundary condition order error at ``dt = 0.1`` + ========================= ========== ===================== + periodic (best possible) 6.32 7.4e-08 + pointwise 5.74 9.7e-07 + differentiated 6.30 1.3e-07 + ========================= ========== ===================== + + The remaining factor of 1.8 against the periodic case is a constant, not a rate. Note + that the observed orders here are pre-asymptotic -- the periodic reference does not + reach its design order 7 either -- so these numbers show that the remedy works, not + that it restores exactly :math:`2M-1`. + + Called once per step by :class:`generic_implicit_mass_diffbc`; ``solve_system`` then + picks the condition belonging to the node it is asked to solve at. + + Parameters + ---------- + t0 : float + Left end of the step. + dt : float + Step size. + coll : pySDC.core.collocation.CollBase + Collocation rule of the sweeper, supplying the nodes and the matrix Q. + """ + M = coll.num_nodes + Q = coll.Qmat[1:, 1:] + self._node_times = t0 + dt * np.asarray(coll.nodes) + + u_rate, p_rate = [], [] + for j in range(M): + self.u_dot.t = self._node_times[j] + self.p_dot.t = self._node_times[j] + u_rate.append(df.interpolate(self.u_dot, self.V)) + p_rate.append(df.interpolate(self.p_dot, self.Q)) + + self.u_ex.t = t0 + self.p_ex.t = t0 + u_base = df.interpolate(self.u_ex, self.V) + p_base = df.interpolate(self.p_ex, self.Q) + + self._node_bcs = [] + for m in range(M): + gu, gp = df.Function(self.V), df.Function(self.Q) + gu.assign(u_base) + gp.assign(p_base) + for j in range(M): + gu.vector().axpy(dt * Q[m, j], u_rate[j].vector()) + gp.vector().axpy(dt * Q[m, j], p_rate[j].vector()) + self._node_bcs.append( + self.bc_fixed + + [ + df.DirichletBC(self.W.sub(0), gu, self.left_right), + df.DirichletBC(self.W.sub(1), gp, self.left_right), + ] + ) + def solve_system(self, rhs, factor, u0, t): r""" Newton solver for :math:`M w + factor \cdot N(w, t) = rhs`, where :math:`N` collects the @@ -249,6 +377,14 @@ def solve_system(self, rhs, factor, u0, t): self.p_ex.t = t self.g.t = t + if self.differentiated_bc: + if self._node_bcs is None: + raise RuntimeError( + 'differentiated_bc requires the generic_implicit_mass_diffbc sweeper, ' + 'which calls prepare_step once per step' + ) + self.bc = self._node_bcs[int(np.argmin(np.abs(self._node_times - t)))] + self.w.vector()[:] = u0.values.vector()[:] self.step.rhs = rhs.values.vector() self.step.bcs = self.bc diff --git a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py index adc9d9cbdc..0c97f03e53 100644 --- a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py @@ -5,10 +5,23 @@ from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import ( fenics_NSE_2D_TaylorGreen, ) -from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import generic_implicit_mass +from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import ( + generic_implicit_mass, + generic_implicit_mass_diffbc, +) -def setup(t0=0.0, dt=0.1, periodic=False, nelems=24, nu=0.1, num_nodes=4, maxiter=40, restol=1e-12): +def setup( + t0=0.0, + dt=0.1, + periodic=False, + differentiated_bc=False, + nelems=24, + nu=0.1, + num_nodes=4, + maxiter=40, + restol=1e-12, +): """ Helper routine to set up parameters @@ -19,6 +32,9 @@ def setup(t0=0.0, dt=0.1, periodic=False, nelems=24, nu=0.1, num_nodes=4, maxite time step size periodic: bool, use periodic instead of time-dependent Dirichlet conditions in x + differentiated_bc: bool, + impose the time-dependent boundary data in differentiated form, which recovers + the order it otherwise costs; requires periodic=False nelems: int, number of elements per spatial direction nu: float, @@ -58,6 +74,7 @@ def setup(t0=0.0, dt=0.1, periodic=False, nelems=24, nu=0.1, num_nodes=4, maxite problem_params['order'] = 2 problem_params['nu'] = nu problem_params['periodic'] = periodic + problem_params['differentiated_bc'] = differentiated_bc problem_params['Sol_tol'] = 1e-13 # initialize controller parameters @@ -67,7 +84,7 @@ def setup(t0=0.0, dt=0.1, periodic=False, nelems=24, nu=0.1, num_nodes=4, maxite # Fill description dictionary description = dict() description['problem_class'] = fenics_NSE_2D_TaylorGreen - description['sweeper_class'] = generic_implicit_mass + description['sweeper_class'] = generic_implicit_mass_diffbc if differentiated_bc else generic_implicit_mass description['problem_params'] = problem_params description['sweeper_params'] = sweeper_params description['level_params'] = level_params @@ -227,12 +244,17 @@ def main(): Tend = 0.2 dts = [0.2, 0.1, 0.05, 0.025] + cases = [ + ('periodic', dict(periodic=True)), + ('time-dependent Dirichlet', dict(periodic=False)), + ('time-dependent Dirichlet, differentiated', dict(periodic=False, differentiated_bc=True)), + ] + results = {} - for periodic in (True, False): - dts_out, errors_u, errors_p = order_study(dts, Tend, periodic=periodic) - results[periodic] = (dts_out, errors_u, errors_p) + for label, kwargs in cases: + dts_out, errors_u, errors_p = order_study(dts, Tend, **kwargs) + results[label] = (dts_out, errors_u, errors_p) - label = 'periodic' if periodic else 'time-dependent Dirichlet' print(f'\n{label} boundary conditions in x:') print(f'{"dt":>10} {"err(u)":>12} {"order(u)":>9} {"err(p)":>12} {"order(p)":>9}') orders_u = [None] + observed_order(dts_out, errors_u) diff --git a/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py b/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py index d98e7d1735..302d401f32 100755 --- a/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py +++ b/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py @@ -144,3 +144,29 @@ def compute_residual(self, stage=None): L.status.updated = False return None + + +class generic_implicit_mass_diffbc(generic_implicit_mass): + """ + Variant of ``generic_implicit_mass`` for problems that impose their boundary conditions in + differentiated form. + + Such a problem needs the collocation data of the current step to build the boundary values + of its stages by quadrature, and a problem class cannot see that data on its own: it only + ever learns the time of the node it is asked to solve at. This sweeper hands it over once + per step, before sweeping. + + The problem class must provide ``prepare_step(t0, dt, coll)``. + """ + + def update_nodes(self): + """ + Supply the collocation data of this step to the problem, then sweep as usual. + + Returns: + None + """ + L = self.level + L.prob.prepare_step(L.time, L.dt, self.coll) + + return super().update_nodes() diff --git a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py index 2a8ac5a2ae..97238cc696 100644 --- a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py @@ -180,3 +180,58 @@ def test_solve_system(periodic): rel_err = abs(w - uex) / abs(uex) assert rel_err < 1e-9, f"solve_system did not recover the exact solution: {rel_err:.3e}" + + +@pytest.mark.fenics +def test_differentiated_boundary_condition(): + r""" + Imposing the time-dependent boundary data in differentiated form recovers most of the order + it otherwise costs. + + Asserted on the error rather than on the observed order: the order estimates are not clean + enough on a mesh CI can afford. The *periodic* reference itself only reaches about 6.3 here + instead of its design order 7, and on finer step sizes it collapses to 5 against a solver + floor near 1e-10, so separating 2M-1 from 2M-2 is beyond what this benchmark resolves. The + error, on the other hand, is unambiguous: with the differentiated condition it lands close + to the periodic case, while the pointwise one is an order of magnitude away. + """ + from pySDC.projects.StroemungsRaum.run_Navier_Stokes_TaylorGreen_FEniCS import order_study + + Tend, dts = 0.4, [0.2, 0.1] + errors = {} + for label, kwargs in [ + ('periodic', dict(periodic=True)), + ('pointwise', dict(periodic=False)), + ('differentiated', dict(periodic=False, differentiated_bc=True)), + ]: + _, _, errors_p = order_study(dts, Tend, num_nodes=4, restol=1e-13, **kwargs) + errors[label] = errors_p[0] + + # the differentiated condition must be a clear improvement on the pointwise one ... + gain = errors['pointwise'] / errors['differentiated'] + assert gain > 2.5, f"differentiated boundary condition only improves the error by {gain:.1f}x" + + # ... and land near the periodic case, which is the best this discretization can do + remaining = errors['differentiated'] / errors['periodic'] + assert remaining < 4.0, f"differentiated error is still {remaining:.1f}x the periodic one" + + # sanity: the pointwise variant is the one that is far off + assert errors['pointwise'] / errors['periodic'] > 5.0, "pointwise variant unexpectedly accurate" + + +@pytest.mark.fenics +def test_differentiated_boundary_condition_needs_its_sweeper(): + """ + ``differentiated_bc`` silently doing nothing would be worse than failing, since the run + would look fine and just be less accurate. Check both guards. + """ + from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import ( + fenics_NSE_2D_TaylorGreen, + ) + + prob = fenics_NSE_2D_TaylorGreen(nelems=8, nu=0.05, differentiated_bc=True) + with pytest.raises(RuntimeError, match='prepare_step'): + prob.solve_system(prob.u_exact(0.0), 0.01, prob.dtype_u(prob.W), 0.0) + + with pytest.raises(ValueError, match='time-dependent'): + fenics_NSE_2D_TaylorGreen(nelems=8, nu=0.05, periodic=True, differentiated_bc=True) From 698b6d8b0562b06ab5075a28459cf46497cf90b1 Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Thu, 10 Sep 2026 19:50:23 +0200 Subject: [PATCH 6/7] Rename the funding ministry from BMBF to BMFTR The ministry is now the Federal Ministry of Research, Technology and Space. The main README already carried the BMFTR logo while its text still said BMBF, so both it and the StroemungsRaum project README are updated, along with the link, which now points at bmftr.bund.de (verified to resolve). Grant numbers are unchanged. Co-Authored-By: Claude Opus 5 --- README.md | 4 ++-- pySDC/projects/StroemungsRaum/README.rst | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 76d4d78c75..124b8a23f6 100644 --- a/README.md +++ b/README.md @@ -119,8 +119,8 @@ and grant agreement No 101118139. The JU receives support from the European Union's Horizon 2020 research and innovation programme and Belgium, France, Germany, and Switzerland. This project also received funding from the [German Federal Ministry of -Education and Research](https://www.bmbf.de/bmbf/en/home/home_node.html) -(BMBF) grants 16HPC047, 16ME0708 and 16ME0679K. Supported by the European Union - NextGenerationEU. +Research, Technology and Space](https://www.bmftr.bund.de/EN/Home/home_node.html) +(BMFTR) grants 16HPC047, 16ME0708 and 16ME0679K. Supported by the European Union - NextGenerationEU. The project also received help from the [Joint Lab "Helmholtz Information - Research Software Engineering" (HiRSE)](https://www.helmholtz-hirse.de/).

diff --git a/pySDC/projects/StroemungsRaum/README.rst b/pySDC/projects/StroemungsRaum/README.rst index 340c2c77ed..a4afa2bc29 100644 --- a/pySDC/projects/StroemungsRaum/README.rst +++ b/pySDC/projects/StroemungsRaum/README.rst @@ -2,7 +2,7 @@ StroemungsRaum ============== **StroemungsRaum** is a research software project developed within the -BMBF-funded project +BMFTR-funded project *“StrömungsRaum – Novel Exascale Architectures with Heterogeneous Hardware Components for Computational Fluid Dynamics Simulations”* @@ -62,6 +62,6 @@ roughly an order of magnitude, to within a small factor of the periodic case. Funding ------- -Funded by the **German Federal Ministry of Education and Research (BMBF)** under -grant number **16ME0708**. +Funded by the **German Federal Ministry of Research, Technology and Space (BMFTR)** +under grant number **16ME0708**. From 19edb416bce567a18256ffb06ea1aaa27ff808bd Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Fri, 11 Sep 2026 11:55:54 +0200 Subject: [PATCH 7/7] Address the review of the Taylor-Green benchmark Nine findings from a review of this branch. Boundary conditions per node are now looked up by exact time rather than by nearest match, so a `solve_system` call at anything but a node of the prepared step fails instead of quietly using a neighbouring node's data. The arithmetic on both sides is identical, so the exact lookup costs nothing. `prepare_step` moves from `update_nodes` to `predict`, which the controller calls once per step rather than once per sweep -- what both docstrings already claimed. Measured over three steps at twelve iterations each: three calls instead of thirty-six. The Newton solver drops its `relative_tolerance = Sol_tol` and its reduced iteration budget, both of which pushed it towards raising on a marginal stall with `error_on_nonconvergence` left at its default. It now matches the merged monolithic class. The numbers in the `prepare_step` table were regenerated with this change in place and are unmoved: 7.411e-08 at dt = 0.1, order 6.32. `prepare_step` also gains the guard for the one mismatch that was missing, the `generic_implicit_mass_diffbc` sweeper on a problem with `differentiated_bc` off, which used to die on a missing attribute. `test_eval_f` takes `du/dt` from `_boundary_derivatives` instead of re-typing it, which anchors the velocity derivative the differentiated condition is built from, and checks its pressure half against a finite difference of `p_ex`, which nothing covered before. The rest is removal: the `dt_ref` branch of `order_study` that no caller used, the `results` dictionary `main` returned to nobody, and the `degree_rise` enrichment in `relative_errors`, which is pure cost for two functions in the same space -- verified to leave the result identical to four digits. `Sol_tol` becomes a keyword on `setup` defaulting to one decade below `restol`. A tenth finding, that the table in `prepare_step` mislabels its error column, was wrong: the numbers come from a four-point study whose dt = 0.1 entry is 7.411e-08 with order 6.32, so the label is correct and stands. Co-Authored-By: Claude Opus 5 --- ...Stokes_2D_TaylorGreen_monolithic_FEniCS.py | 16 ++++++-- .../run_Navier_Stokes_TaylorGreen_FEniCS.py | 41 +++++++------------ .../sweepers/generic_implicit_mass.py | 10 +++-- .../test_Navier_Stokes_TaylorGreen_FEniCS.py | 27 ++++++------ 4 files changed, 49 insertions(+), 45 deletions(-) diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py index a03c082025..5d3c8eceed 100644 --- a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_TaylorGreen_monolithic_FEniCS.py @@ -244,8 +244,6 @@ def __init__(self, nelems=32, t0=0.0, order=2, nu=0.02, periodic=False, differen self.step = NewtonStep(F, df.derivative(F, self.w)) self.newton = df.NewtonSolver() self.newton.parameters['absolute_tolerance'] = Sol_tol - self.newton.parameters['relative_tolerance'] = Sol_tol - self.newton.parameters['maximum_iterations'] = 20 @staticmethod def _boundary_derivatives(nu, order, t0): @@ -318,6 +316,12 @@ def prepare_step(self, t0, dt, coll): coll : pySDC.core.collocation.CollBase Collocation rule of the sweeper, supplying the nodes and the matrix Q. """ + if not self.differentiated_bc: + raise RuntimeError( + 'prepare_step builds the differentiated boundary conditions, which this problem ' + 'was not set up for; use generic_implicit_mass or pass differentiated_bc=True' + ) + M = coll.num_nodes Q = coll.Qmat[1:, 1:] self._node_times = t0 + dt * np.asarray(coll.nodes) @@ -383,7 +387,13 @@ def solve_system(self, rhs, factor, u0, t): 'differentiated_bc requires the generic_implicit_mass_diffbc sweeper, ' 'which calls prepare_step once per step' ) - self.bc = self._node_bcs[int(np.argmin(np.abs(self._node_times - t)))] + node = np.flatnonzero(self._node_times == t) + if node.size != 1: + raise RuntimeError( + f'no collocation node of the prepared step is at t = {t}; the prepared step ' + f'covers {self._node_times}' + ) + self.bc = self._node_bcs[node[0]] self.w.vector()[:] = u0.values.vector()[:] self.step.rhs = rhs.values.vector() diff --git a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py index 0c97f03e53..0d46edb290 100644 --- a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py @@ -21,6 +21,7 @@ def setup( num_nodes=4, maxiter=40, restol=1e-12, + Sol_tol=None, ): """ Helper routine to set up parameters @@ -45,6 +46,9 @@ def setup( maximum number of SDC iterations restol: float, residual tolerance + Sol_tol: float, + absolute tolerance of the Newton solve at each node; defaults to one decade below + ``restol``, which SDC cannot converge past Returns: description: dict, @@ -75,7 +79,7 @@ def setup( problem_params['nu'] = nu problem_params['periodic'] = periodic problem_params['differentiated_bc'] = differentiated_bc - problem_params['Sol_tol'] = 1e-13 + problem_params['Sol_tol'] = restol / 10 if Sol_tol is None else Sol_tol # initialize controller parameters controller_params = dict() @@ -140,8 +144,8 @@ def relative_errors(u, uref): ur, pr = uref.values.split(deepcopy=True) return ( - df.errornorm(ur, un, 'L2') / df.norm(ur, 'L2'), - df.errornorm(pr, pn, 'L2') / df.norm(pr, 'L2'), + df.errornorm(ur, un, 'L2', degree_rise=0) / df.norm(ur, 'L2'), + df.errornorm(pr, pn, 'L2', degree_rise=0) / df.norm(pr, 'L2'), ) @@ -163,26 +167,21 @@ def run_postprocessing(P, uend, Tend): return relative_errors(uend, P.u_exact(Tend)) -def order_study(dts, Tend, dt_ref=None, periodic=False, **kwargs): +def order_study(dts, Tend, periodic=False, **kwargs): r""" Measure the observed temporal order of convergence. Errors are *not* taken against the exact solution: the spatial discretization error dominates it for any affordable mesh, which hides the temporal order completely. Instead - two variants are offered, both of which cancel the spatial error exactly because every run - uses the same mesh: - - - ``dt_ref`` given: compare against a reference run with that much smaller step size, - - ``dt_ref`` omitted: compare consecutive step sizes with each other (Richardson). This - needs no reference run and is therefore a lot cheaper, at the cost of one order estimate. + consecutive step sizes are compared with each other (Richardson), which cancels the spatial + error exactly because every run uses the same mesh and needs no reference run, at the cost + of one order estimate. Args: dts: list of float, Step sizes to run, largest first, each one half of the previous. Tend: float, Final simulation time; must be an integer multiple of every step size. - dt_ref: float, - Step size for the reference run, or ``None`` to compare consecutive step sizes. periodic: bool, Use periodic instead of time-dependent Dirichlet conditions in x. kwargs: @@ -190,7 +189,7 @@ def order_study(dts, Tend, dt_ref=None, periodic=False, **kwargs): Returns: dts_out: list of float, - Step sizes the errors belong to; one shorter than ``dts`` without a reference. + Step sizes the errors belong to; one shorter than ``dts``. errors_u: list of float, Relative L2 velocity error per step size. errors_p: list of float, @@ -201,18 +200,10 @@ def order_study(dts, Tend, dt_ref=None, periodic=False, **kwargs): description, controller_params = setup(dt=dt, periodic=periodic, **kwargs) solutions.append(run_simulation(description, controller_params, Tend)[2]) - if dt_ref is None: - pairs = list(zip(solutions[:-1], solutions[1:], strict=True)) - dts_out = dts[:-1] - else: - description, controller_params = setup(dt=dt_ref, periodic=periodic, **kwargs) - uref = run_simulation(description, controller_params, Tend)[2] - pairs = [(u, uref) for u in solutions] - dts_out = list(dts) - + pairs = zip(solutions[:-1], solutions[1:], strict=True) errors = [relative_errors(u, ref) for u, ref in pairs] - return dts_out, [e[0] for e in errors], [e[1] for e in errors] + return dts[:-1], [e[0] for e in errors], [e[1] for e in errors] def observed_order(dts, errors): @@ -250,10 +241,8 @@ def main(): ('time-dependent Dirichlet, differentiated', dict(periodic=False, differentiated_bc=True)), ] - results = {} for label, kwargs in cases: dts_out, errors_u, errors_p = order_study(dts, Tend, **kwargs) - results[label] = (dts_out, errors_u, errors_p) print(f'\n{label} boundary conditions in x:') print(f'{"dt":>10} {"err(u)":>12} {"order(u)":>9} {"err(p)":>12} {"order(p)":>9}') @@ -264,8 +253,6 @@ def main(): sp = ' --- ' if op is None else f'{op:9.2f}' print(f'{dt:10.5f} {eu:12.4e} {su} {ep:12.4e} {sp}') - return results - if __name__ == "__main__": main() diff --git a/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py b/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py index 302d401f32..ace17346c2 100755 --- a/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py +++ b/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py @@ -159,9 +159,13 @@ class generic_implicit_mass_diffbc(generic_implicit_mass): The problem class must provide ``prepare_step(t0, dt, coll)``. """ - def update_nodes(self): + def predict(self): """ - Supply the collocation data of this step to the problem, then sweep as usual. + Supply the collocation data of this step to the problem, then predict as usual. + + ``predict`` rather than ``update_nodes`` because the controller calls it exactly once + per step, while ``update_nodes`` runs once per sweep and would rebuild the same + boundary conditions on every iteration. Returns: None @@ -169,4 +173,4 @@ def update_nodes(self): L = self.level L.prob.prepare_step(L.time, L.dt, self.coll) - return super().update_nodes() + return super().predict() diff --git a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py index 97238cc696..39fc5a29ac 100644 --- a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py +++ b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_TaylorGreen_FEniCS.py @@ -49,6 +49,10 @@ def test_eval_f(): Compared against du/dt rather than against ``solve_system`` on purpose -- a sign error shared by both would pass a consistency check between them. + + ``du/dt`` is taken from ``_boundary_derivatives`` rather than re-typed, so this also anchors + the velocity derivative the differentiated boundary condition is built from. Its pressure + half has no such anchor, so it is checked here against a finite difference of ``p_ex``. """ import dolfin as df from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import ( @@ -60,18 +64,7 @@ def test_eval_f(): for nelems in (16, 32): prob = fenics_NSE_2D_TaylorGreen(nelems=nelems, t0=0.0, order=2, nu=nu) - dudt = df.Expression( - ( - '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])' - ' + 2*pi*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*sin(pi*x[1])*cos(pi*x[1])', - '8*pi*pi*nu*exp(-8*pi*pi*nu*t)*cos(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])' - ' - 2*pi*exp(-8*pi*pi*nu*t)*sin(2*pi*(x[0] - t))*cos(pi*x[1])*cos(pi*x[1])', - ), - pi=np.pi, - nu=nu, - t=t, - degree=prob.order + 2, - ) + dudt, dpdt = prob._boundary_derivatives(nu, prob.order, t) ut = prob.dtype_u(prob.W) df.assign(ut.values.sub(0), df.interpolate(dudt, prob.V)) @@ -87,6 +80,16 @@ def test_eval_f(): b = expected.values.vector()[velocity_dofs] errors.append(np.linalg.norm(a - b) / np.linalg.norm(b)) + # p_dot feeds the differentiated pressure boundary condition and nothing else pins it + h = 1e-6 + for x, y in ((0.5, 0.2), (-0.5, -0.35)): + prob.p_ex.t = t + h + fwd = prob.p_ex(x, y) + prob.p_ex.t = t - h + bwd = prob.p_ex(x, y) + fd = (fwd - bwd) / (2 * h) + assert abs(fd - dpdt(x, y)) < 1e-6, f"p_dot is not dp_ex/dt at ({x}, {y}): {fd} vs {dpdt(x, y)}" + assert errors[0] < 2e-2, f"eval_f does not match M du/dt: relative error {errors[0]:.3e}" order = np.log2(errors[0] / errors[1]) assert order > 1.7, f"eval_f converges at order {order:.2f}, expected second order"