Add Taylor–Green NS benchmark studying order reduction induced by time-dependent bcs

[Ouardghi]

This PR adds a 2D Taylor–Green incompressible Navier–Stokes example implemented
with FEniCSx, together with a small runnable example and tests.

The example uses a monolithic velocity–pressure formulation and includes:

• time-dependent and time-independent Dirichlet boundary conditions,
• periodic boundary condition,
• a manufactured analytical solution for verification.

The main purpose of this benchmark is to study the effect of time-dependent
boundary conditions on numerical accuracy and order reduction in time
integration methods.

Added

• Monolithic NSE formulation and Taylor–Green benchmark setup with periodic and non-periodic bcs variants
• Small run example
• Convergence and correctness tests

[Copilot]

Pull request overview

This PR introduces a new 2D Taylor–Green incompressible Navier–Stokes benchmark for the StroemungsRaum_FEniCSx project, intended to study order reduction effects from time-dependent boundary conditions, with a runnable driver script and pytest-based verification.

Changes:

• Added a monolithic velocity–pressure FEniCSx Navier–Stokes problem class with both non-periodic and periodic BC variants.
• Added a runnable script to set up, run, and postprocess the benchmark.
• Added tests covering system solve, RHS evaluation, and an end-to-end run for both BC variants; plus project README and a conda environment spec.

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 4 comments.

Show a summary per file

File	Description
pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py	Implements the new monolithic NSE Taylor–Green benchmark (non-periodic + periodic via MPC).
pySDC/projects/StroemungsRaum_FEniCSx/run_Navier_Stokes_equations_TaylGreen_monolithic_FEniCSx.py	Provides a small runnable setup/simulation/postprocessing entry point for the benchmark.
pySDC/projects/StroemungsRaum_FEniCSx/tests/test_Navier_Stokes_TaylGreen_monolithic_FEniCSx.py	Adds unit/integration tests validating solve, RHS evaluation, and convergence for both BC variants.
pySDC/projects/StroemungsRaum_FEniCSx/README.rst	Adds project-level documentation for the StroemungsRaum_FEniCSx subtree.
pySDC/projects/StroemungsRaum_FEniCSx/environment.yml	Adds a conda environment file intended to provision the FEniCSx stack for this project.

Comments suppressed due to low confidence (7)

pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py:227

• self.factor is initialized to 0.0 but is used in the nonlinear form as 1 / self.factor (see F = (1 / self.factor) * ...). This makes the form ill-defined (division by zero / Inf) before solve_system sets a nonzero value and can lead to runtime failures or NaNs during form compilation/assembly. Consider formulating the residual without dividing by factor (e.g., multiply the whole equation by factor) and/or initialize self.factor to a safe nonzero value while explicitly handling the factor == 0 case.

        # Initialize the nonlinear form for the Navier–Stokes equations
        self.factor = dfx.fem.Constant(domain, PETSc.ScalarType(0.0))

        self.w = dfx.fem.Function(self.W)

        u, p = ufl.split(self.w)
        rhs_u, rhs_p = ufl.split(self.rhs)

        # define the nonlinear form to be compatible with the multi-point constraints
        F = (1 / self.factor) * ufl.dot(u, self.v) * ufl.dx
        F += ufl.dot(ufl.dot(u, ufl.nabla_grad(u)), self.v) * ufl.dx
        F += ufl.inner(self.nu * ufl.nabla_grad(u), ufl.nabla_grad(self.v)) * ufl.dx
        F -= ufl.dot(p, ufl.div(self.v)) * ufl.dx
        F -= ufl.dot(self.g, self.v) * ufl.dx
        F -= ufl.dot(ufl.div(u), self.q) * ufl.dx
        F -= (1 / self.factor) * ufl.dot(rhs_u, self.v) * ufl.dx
        F -= (1 / self.factor) * ufl.dot(rhs_p, self.q) * ufl.dx

pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py:301

• dolfinx.fem.petsc.NonlinearProblem does not provide a .solve() method in the standard DOLFINx API (it is normally solved via dolfinx.nls.petsc.NewtonSolver). As written, problem.solve() is likely to raise AttributeError at runtime. Please instantiate a Newton solver (or SNES solver wrapper, depending on DOLFINx version) and call its solve method with the NonlinearProblem and solution function.

        problem = dolfinx.fem.petsc.NonlinearProblem(
            self.F, self.w, bcs=bcup, petsc_options=self.petsc_options, petsc_options_prefix="nonlinearsolver"
        )

        # solve the nonlinear system using Dolfinx's nonlinear solver
        problem.solve()

pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py:124

• This class prints the number of DoFs unconditionally. In this codebase similar diagnostics are logged via the problem logger (see e.g. pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py:108-111). Please switch this print to self.logger.debug(...) (or remove it) to avoid noisy test output and improve controllability.

        # get number of dofs and set up the problem
        self.tmp_u = dfx.fem.Function(self.W)
        nx = len(self.tmp_u.x.array)
        print('DoFs on this level:', nx)

pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py:585

• xmin/xmax/ymin/ymax are computed via local min()/max() on domain.geometry.x. With MPI, those are rank-local extrema and can differ across ranks, breaking periodic boundary detection/constraint creation. Use communicator reductions (e.g. comm.allreduce(local_min, op=MPI.MIN) / MPI.MAX) to compute global domain extents.

        self.xmin = domain.geometry.x[:, 0].min()
        self.xmax = domain.geometry.x[:, 0].max()
        self.ymin = domain.geometry.x[:, 1].min()
        self.ymax = domain.geometry.x[:, 1].max()
        self.Lx = float(self.xmax - self.xmin)
        self.Ly = float(self.ymax - self.ymin)

pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py:701

• periodic_relation_x unconditionally writes out_x[2] = x[2]. For a 2D mesh created by create_rectangle the coordinate dimension is typically 2, so indexing the third coordinate can raise an IndexError when the relation is evaluated. Please make this dimension-agnostic (only set available components, or guard with if out_x.shape[0] > 2).

        out_x = np.zeros_like(x)
        out_x[0] = x[0] - self.Lx
        out_x[1] = x[1]
        out_x[2] = x[2]
        return out_x

pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py:629

• The Jacobian is derived using ufl.TrialFunction(self.W), but self.w is redefined on W_mpc = self.mpc.function_space. Mixing a trial function from the unconstrained space with a state variable in the constrained space is likely to cause form compilation/assembly errors. Use a TrialFunction from the same space as self.w (e.g. ufl.TrialFunction(W_mpc) / self.w.function_space).

        # compute the Jacobian of the nonlinear form for use in the Newton solver
        Jac = ufl.derivative(F, self.w, ufl.TrialFunction(self.W))

        # assemble the nonlinear problem with the multi-point constraints and configure the Newton solver
        self.problem = dolfinx_mpc.NonlinearProblem(
            F, self.w, mpc=self.mpc, bcs=self.bcs, J=Jac, petsc_options=self.petsc_options
        )

pySDC/projects/StroemungsRaum_FEniCSx/problem_classes/NavierStokes_2D_TaylGreen_monolithic_FEniCSx.py:183

• fix_bc_for_residual is enabled and bc_hom is built as homogeneous Dirichlet conditions on all boundary facets (on_boundary). For the periodic variant, applying Dirichlet BCs on the periodic boundaries in fix_residual will incorrectly zero out residual entries and can harm convergence/accuracy. Consider overriding fix_residual/bc_hom in fenicsx_NSE_periodic_mass (or disabling fix_bc_for_residual) so only true Dirichlet boundaries are constrained and periodicity is handled via the MPC.

        # homogeneous boundary conditions for the velocity and pressure for the residual
        bndry_facets = dfx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
        dofs_bndry_u = dfx.fem.locate_dofs_topological((self.W.sub(0), self.V), fdim, bndry_facets)
        dofs_bndry_p = dfx.fem.locate_dofs_topological((self.W.sub(1), self.Q), fdim, bndry_facets)
        bcu_hom = dfx.fem.dirichletbc(w_zero.sub(0).collapse(), dofs_bndry_u, self.W.sub(0))
        bcp_hom = dfx.fem.dirichletbc(w_zero.sub(1).collapse(), dofs_bndry_p, self.W.sub(1))
        self.bc_hom = [bcu_hom, bcp_hom]
        self.fix_bc_for_residual = True

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.