diff --git a/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py
new file mode 100755
index 0000000000..86fb74b34a
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/problem_classes/NavierStokes_2D_monolithic_FEniCS.py
@@ -0,0 +1,327 @@
+import logging
+import os
+
+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 fenics_NSE_2D_Monolithic(Problem):
+    r"""
+    Example implementing the forced two-dimensional incompressible Navier-Stokes equations with
+    time-dependent Dirichlet boundary conditions
+
+    .. math::
+        \frac{\partial u}{\partial t} = - u \cdot \nabla u  + \nu \nabla u - \nabla p + f
+                      0 = \nabla \cdot u
+
+    for :math:`x \in \Omega`, where the forcing term :math:`f` is defined by
+
+    .. math::
+        f(x, t) = (0, 0).
+
+    This implementation follows a monolithic approach, where velocity and pressure are
+    solved simultaneously in a coupled system using a mixed finite element formulation.
+
+    Boundary conditions are applied on subsets of the boundary:
+    - no-slip on channel walls and cylinder surface,
+    - a time-dependent inflow profile at the inlet,
+    - pressure condition at the outflow.
+
+    In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem
+    is reformulated to the *weak formulation*
+
+    .. math:
+        \int_\Omega u_t v\,dx = - \int_\Omega u \cdot \nabla u v\,dx -  \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega p \nabla \cdot v\,dx + \int_\Omega f v\,dx
+        \int_\Omega \nabla \cdot u q\,dx = 0
+
+    Parameters
+    ----------
+    t0 : float, optional
+        Starting time.
+    order : int, optional
+        Defines the order of the elements in the function space.
+    nu : float, optional
+        Diffusion coefficient :math:`\nu`.
+
+    Attributes
+    ----------
+    V : FunctionSpace
+        Defines the function space of the trial and test functions for velocity.
+    Q : FunctionSpace
+        Defines the function space of the trial and test functions for pressure.
+    W : FunctionSpace
+        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 heat equation.
+    bc : DirichletBC
+        Denotes the time-dependent Dirichlet boundary conditions.
+    bc_hom : DirichletBC
+        Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual
+    fix_bc_for_residual: boolean
+        flag to indicate 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).
+    .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
+        Wells and others. Springer (2012).
+    """
+
+    dtype_u = fenics_mesh
+    dtype_f = fenics_mesh
+
+    df.set_log_active(False)
+
+    def __init__(self, t0=0.0, order=2, nu=0.001):
+
+        # 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
+        df.parameters['allow_extrapolation'] = True
+
+        # load mesh
+        path = f"{os.path.dirname(__file__)}/../meshs/"
+        mesh = df.Mesh(path + '/cylinder.xml')
+
+        # define function spaces for future reference (Taylor-Hood)
+        P2 = df.VectorElement("P", mesh.ufl_cell(), order)
+        P1 = df.FiniteElement("P", mesh.ufl_cell(), order - 1)
+        TH = df.MixedElement([P2, P1])
+        self.W = df.FunctionSpace(mesh, TH)
+        self.V = df.FunctionSpace(mesh, P2)
+        self.Q = df.FunctionSpace(mesh, P1)
+
+        # print the number of DoFs for debugging purposes
+        tmp = df.Function(self.W)
+        self.logger.debug('DoFs on this level:', len(tmp.vector()[:]))
+
+        super().__init__(self.W)
+        self._makeAttributeAndRegister('t0', 'order', 'nu', localVars=locals(), readOnly=True)
+
+        # Trial and test function for the Mixed FE space
+        self.u, self.p = df.TrialFunctions(self.W)
+        self.v, self.q = df.TestFunctions(self.W)
+
+        # velocity mass matrix
+        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
+        self.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)
+
+        # define boundaries
+        inflow = 'near(x[0], 0)'
+        outflow = 'near(x[0], 2.2)'
+        walls = 'near(x[1], 0) || near(x[1], 0.41)'
+        cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'
+
+        # define boundary conditions
+        bc_in = df.DirichletBC(self.W.sub(0), self.u_in, inflow)
+        bc_out = df.DirichletBC(self.W.sub(1), 0, outflow)
+        bc_walls = df.DirichletBC(self.W.sub(0), (0, 0), walls)
+        bc_cylinder = df.DirichletBC(self.W.sub(0), (0, 0), cylinder)
+        self.bc = [bc_cylinder, bc_walls, bc_out, bc_in]
+
+        # homogeneous boundary conditions for fixing the residual
+        bc_hom_u = df.DirichletBC(self.W.sub(0), df.Constant((0, 0)), 'on_boundary')
+        bc_hom_p = df.DirichletBC(self.W.sub(1), df.Constant(0), 'on_boundary')
+        self.bc_hom = [bc_hom_u, bc_hom_p]
+        self.fix_bc_for_residual = True
+
+        # define measure for drag and lift computation
+        Cylinder = df.CompiledSubDomain(cylinder)
+        CylinderBoundary = df.MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
+        Cylinder.mark(CylinderBoundary, 1)
+        self.dsc = df.Measure("ds", domain=mesh, subdomain_data=CylinderBoundary, subdomain_id=1)
+
+        # set forcing term as expression
+        self.g = df.Expression(('0', '0'), a=np.pi, b=self.nu, t=self.t0, degree=self.order)
+
+        # initialize XDMF files for velocity and pressure if needed
+        self.xdmffile_p = None
+        self.xdmffile_u = None
+
+    def solve_system(self, rhs, factor, u0, t):
+        r"""
+        Dolfin's nonlinear solver for :
+
+        (u,v) + factor * (u \cdot \nabla u, v) + factor * \nu (\nabla u, \nabla v) - factor * (p, \nabla \cdot v) - factor * (g, v) - factor * (div(u), q) = (rhs_u, v) + (rhs_p, q)
+
+        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 (not used here so far).
+        t : float
+            Current time.
+
+        Returns
+        -------
+        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
+        self.u_in.t = t
+
+        # get the forcing term
+        self.g.t = t
+        g = df.interpolate(self.g, self.V)
+
+        # 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
+
+        # solve the nonlinear system using Newton's method
+        df.solve(F == 0, w.values, self.bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1e-15}})
+
+        return w
+
+    def eval_f(self, w, t):
+        """
+        Routine to evaluate both parts of the right-hand side of the problem.
+
+        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.
+        """
+
+        f = self.dtype_f(self.W)
+
+        u, p = df.split(w.values)
+
+        # get the forcing term
+        self.g.t = t
+        g = self.dtype_f(df.interpolate(self.g, self.V), val=self.V)
+
+        F = -1.0 * 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(g.values, self.v) * df.dx
+        F += df.dot(df.div(u), self.q) * df.dx
+
+        f.values.vector()[:] = df.assemble(F)
+
+        return f
+
+    def apply_mass_matrix(self, w):
+        r"""
+        Routine to apply 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 __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)
+        df.solve(self.Mf, me.values.vector(), w.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.
+        """
+        # define the exact solution
+        w = df.Function(self.W)
+
+        # assign the exact solution as an Expression
+        df.assign(w.sub(0), df.interpolate(df.Expression(('0.0', '0.0'), degree=self.order), self.V))
+        df.assign(w.sub(1), df.interpolate(df.Expression('0.0', degree=self.order - 1), self.Q))
+
+        # update time in Boundary conditions
+        self.u_in.t = t
+        [bc.apply(w.vector()) for bc in self.bc]
+
+        me = self.dtype_u(w, val=self.W)
+
+        return me
+
+    def fix_residual(self, res):
+        """
+        Applies homogeneous Dirichlet boundary conditions to the residual
+
+        Parameters
+        ----------
+        res : dtype_u
+              Residual
+        """
+        [bc.apply(res.values.vector()) for bc in self.bc_hom]
+
+        return None
diff --git a/pySDC/projects/StroemungsRaum/run_Navier_Stokes_equations_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_equations_monolithic_FEniCS.py
new file mode 100755
index 0000000000..f6a3351e07
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/run_Navier_Stokes_equations_monolithic_FEniCS.py
@@ -0,0 +1,108 @@
+from pathlib import Path
+import json
+
+from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_monolithic_FEniCS import fenics_NSE_2D_Monolithic
+from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import generic_implicit_mass
+
+
+def setup(t0=0):
+    """
+    Helper routine to set up parameters
+
+    Args:
+        t0: float,
+            initial time
+
+    Returns:
+        description: dict,
+            pySDC description dictionary containing problem and method parameters.
+        controller_params: dict,
+            Parameters for the pySDC controller.
+    """
+    # time step size
+    dt = 1 / 100
+
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1e-10
+    level_params['dt'] = dt
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 10
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['quad_type'] = 'RADAU-RIGHT'
+    sweeper_params['num_nodes'] = 2
+    sweeper_params['QI'] = ['LU']
+
+    # initialize problem parameters
+    problem_params = dict()
+    problem_params['nu'] = 0.001
+    problem_params['t0'] = t0
+    problem_params['order'] = 2
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['logger_level'] = 20
+
+    # Fill description dictionary
+    description = dict()
+    description['problem_class'] = fenics_NSE_2D_Monolithic
+    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 Navier-Stokes equations.
+
+     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 containing the final solution and other problem-related information.
+        stats: dict,
+           collected runtime statistics,
+        uend: FEniCS function,
+           Final solution at time Tend.
+    """
+    # get initial time from description
+    t0 = description['problem_params']['t0']
+
+    # quickly generate block of steps
+    controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
+
+    # get initial values on finest level
+    P = controller.MS[0].levels[0].prob
+    uinit = P.u_exact(t0)
+
+    # call main function to get things done...
+    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+
+    return P, stats, uend
+
+
+if __name__ == "__main__":
+
+    t0 = 3.125e-04  # the FEATFLOW data initial time
+    Tend = 8.0
+
+    # run the setup to get description and controller parameters
+    description, controller_params = setup(t0=t0)
+
+    # run the simulation and get the problem, stats and final solution
+    prob, stats, uend = run_simulation(description, controller_params, Tend)
diff --git a/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py b/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py
new file mode 100755
index 0000000000..5a522ccaca
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/sweepers/generic_implicit_mass.py
@@ -0,0 +1,124 @@
+from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+
+
+class generic_implicit_mass(generic_implicit):
+
+    def update_nodes(self):
+        """
+        This sweeper extends the generic_implicit sweeper from the implementations
+        package for a monolithic discretization of the incompressible Navier-Stokes
+        equations with a mass matrix. It updates the solution and right-hand side
+        values at all collocation nodes during a single sweep.
+
+        Returns:
+            None
+        """
+
+        L = self.level
+        P = L.prob
+
+        # only if the level has been touched before
+        assert L.status.unlocked
+
+        # get number of collocation nodes for easier access
+        M = self.coll.num_nodes
+
+        # update the MIN-SR-FLEX preconditioner
+        if self.params.QI == 'MIN-SR-FLEX':
+            self.QI = self.get_Qdelta_implicit(qd_type="MIN-SR-FLEX", k=L.status.sweep)
+
+        # gather all terms which are known already (e.g. from the previous iteration)
+        # this corresponds to u0 + QF(u^k) - QdF(u^k) + tau
+
+        # get QF(u^k)
+        integral = self.integrate()
+
+        # This is somewhat ugly, but we have to apply the mass matrix on u0 only on the finest level
+        if L.level_index == 0:
+            u0 = P.apply_mass_matrix(L.u[0])
+        else:
+            u0 = L.u[0]
+
+        for m in range(M):
+            # get -QdF(u^k)_m
+            for j in range(1, M + 1):
+                integral[m] -= L.dt * self.QI[m + 1, j] * L.f[j]
+
+            # add initial value
+            integral[m] += u0
+            # add tau if associated
+            if L.tau[m] is not None:
+                integral[m] += L.tau[m]
+
+        # do the sweep
+        for m in range(0, M):
+            # build rhs, consisting of the known values from above and new values from previous nodes (at k+1)
+            rhs = P.dtype_u(integral[m])
+            for j in range(1, m + 1):
+                rhs += L.dt * self.QI[m + 1, j] * L.f[j]
+
+            # implicit solve with prefactor stemming from the diagonal of Qd
+            alpha = L.dt * self.QI[m + 1, m + 1]
+            if alpha == 0:
+                L.u[m + 1] = rhs
+            else:
+                L.u[m + 1] = P.solve_system(rhs, alpha, L.u[m + 1], L.time + L.dt * self.coll.nodes[m])
+            # update function values
+            L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m])
+
+        # indicate presence of new values at this level
+        L.status.updated = True
+
+        return None
+
+    def compute_residual(self, stage=None):
+        """
+        Computation of the residual using the collocation matrix Q
+
+        Args:
+            stage (str): The current stage of the step the level belongs to
+        """
+
+        # get current level and problem description
+        L = self.level
+        P = L.prob
+
+        # Check if we want to skip the residual computation to gain performance
+        # Keep in mind that skipping any residual computation is likely to give incorrect outputs of the residual!
+        if stage in self.params.skip_residual_computation:
+            L.status.residual = 0.0 if L.status.residual is None else L.status.residual
+            return None
+
+        # check if there are new values (e.g. from a sweep)
+        # assert L.status.updated
+
+        # compute the residual for each node
+
+        # build QF(u)
+        res_norm = []
+        res = self.integrate()
+        for m in range(self.coll.num_nodes):
+
+            # This is somewhat ugly, but we have to apply the mass matrix on u0 only on the finest level
+
+            if L.level_index == 0:
+                res[m] += P.apply_mass_matrix(L.u[0] - L.u[m + 1])
+            else:
+                res[m] += L.u[0] - P.apply_mass_matrix(L.u[m + 1])
+            # add tau if associated
+            if L.tau[m] is not None:
+                res[m] += L.tau[m]
+
+            # Due to different boundary conditions we might have to fix the residual
+            if L.prob.fix_bc_for_residual:
+                L.prob.fix_residual(res[m])
+            # use abs function from data type here
+            res_norm.append(abs(res[m]))
+
+        # find maximal residual over the nodes
+        L.status.residual = max(res_norm)
+
+        # indicate that the residual has seen the new values
+        L.status.updated = False
+
+        return None
diff --git a/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_monolithic_FEniCS.py b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_monolithic_FEniCS.py
new file mode 100644
index 0000000000..526da5b400
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/tests/test_Navier_Stokes_monolithic_FEniCS.py
@@ -0,0 +1,176 @@
+import pytest
+import dolfin as df
+
+from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh
+from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_monolithic_FEniCS import fenics_NSE_2D_Monolithic
+
+
+@pytest.mark.fenics
+def test_solve_system():
+
+    # Physical and numerical parameters for the test case
+    nu = 0.001
+    t = 0.0
+    factor = 0.01
+
+    # Create the 2D monolithic Navier-Stokes problem
+    prob = fenics_NSE_2D_Monolithic(t0=0.0, order=2, nu=nu)
+
+    # Exact Taylor-Green vortex velocity used as manufactured solution
+    u_exact = df.Expression(
+        ("-cos(x[0])*sin(x[1])*exp(-2.0*nu*t)", " sin(x[0])*cos(x[1])*exp(-2.0*nu*t)"), degree=prob.order, nu=nu, t=t
+    )
+
+    # Exact pressure corresponding to the Taylor-Green vortex
+    p_exact = df.Expression("-0.25*(cos(2.0*x[0]) + cos(2.0*x[1]))*exp(-4.0*nu*t)", degree=prob.order - 1, nu=nu, t=t)
+
+    # Interpolate the exact velocity and pressure into the mixed space
+    wex = df.Function(prob.W)
+    df.assign(wex.sub(0), df.interpolate(u_exact, prob.V))
+    df.assign(wex.sub(1), df.interpolate(p_exact, prob.Q))
+    wex = prob.dtype_u(wex)
+
+    # Manufactured forcing term chosen so that the exact solution satisfies
+    # the discrete Navier-Stokes system used in the problem class
+    rhs_u = df.Expression(
+        (
+            "(-1.0-2.0*nu*factor)*cos(x[0])*sin(x[1])*exp(-2.0*nu*t)",
+            "(1.0-2.0*nu*factor)*sin(x[0])*cos(x[1])*exp(-2.0*nu*t)",
+        ),
+        degree=prob.order,
+        nu=nu,
+        t=t,
+        factor=factor,
+    )
+
+    # Zero pressure contribution in the right-hand side
+    rhs_p = df.Expression("0.0", degree=prob.order - 1, nu=nu, t=t)
+
+    # Interpolate the manufactured forcing term into the mixed space
+    g = df.Function(prob.W)
+    df.assign(g.sub(0), df.interpolate(rhs_u, prob.V))
+    df.assign(g.sub(1), df.interpolate(rhs_p, prob.Q))
+
+    # Dirichlet boundary conditions taken from the exact solution
+    bc_u = df.DirichletBC(prob.W.sub(0), u_exact, 'on_boundary')
+    bc_p = df.DirichletBC(prob.W.sub(1), p_exact, 'on_boundary')
+    prob.bc = [bc_u, bc_p]
+
+    # Apply the mass matrix to obtain the final right-hand side vector
+    rhs = prob.apply_mass_matrix(prob.dtype_u(g))
+
+    # Solve the linear system
+    w = prob.solve_system(rhs=rhs, factor=factor, u0=fenics_mesh(prob.W), t=t)
+
+    # Compute the relative error with respect to the exact solution
+    rel_err = abs(wex - w) / abs(wex)
+    assert rel_err < 5e-5, f"Relative error {rel_err} exceeds tolerance"
+
+
+@pytest.mark.fenics
+def test_eval_f():
+
+    # Physical and numerical parameters for the test case
+    nu = 0.1
+    t = 0.0
+
+    # Create the 2D monolithic Navier-Stokes problem
+    prob = fenics_NSE_2D_Monolithic(t0=0.0, order=2, nu=nu)
+
+    # Exact Taylor-Green vortex velocity used as manufactured solution
+    u_exact = df.Expression(
+        ("-cos(x[0])*sin(x[1])*exp(-2.0*nu*t)", " sin(x[0])*cos(x[1])*exp(-2.0*nu*t)"), degree=prob.order, nu=nu, t=t
+    )
+
+    # Exact pressure corresponding to the Taylor-Green vortex
+    p_exact = df.Expression("-0.25*(cos(2.0*x[0]) + cos(2.0*x[1]))*exp(-4.0*nu*t)", degree=prob.order - 1, nu=nu, t=t)
+
+    # Interpolate the exact velocity and pressure into the mixed space
+    wex = df.Function(prob.W)
+    df.assign(wex.sub(0), df.interpolate(u_exact, prob.V))
+    df.assign(wex.sub(1), df.interpolate(p_exact, prob.Q))
+    wex = prob.dtype_u(wex)
+
+    # Evaluate the right-hand side for the given state and time
+    f = prob.eval_f(w=wex, t=t)
+
+    # Analytical value of the right-hand side for the Taylor-Green vortex
+    fu_exact = df.Expression(
+        ("2.0*nu*cos(x[0])*sin(x[1])*exp(-2.0*nu*t)", "-2.0*nu*sin(x[0])*cos(x[1])*exp(-2.0*nu*t)"),
+        degree=prob.order,
+        nu=nu,
+        t=t,
+    )
+
+    # The pressure component of the right-hand side is zero
+    fp_exact = df.Expression("0.0", degree=prob.order - 1, nu=nu, t=t)
+
+    # Interpolate the analytical right-hand side into the mixed space
+    fex = df.Function(prob.W)
+    df.assign(fex.sub(0), df.interpolate(fu_exact, prob.V))
+    df.assign(fex.sub(1), df.interpolate(fp_exact, prob.Q))
+
+    # Apply the mass matrix to obtain the expected right-hand side vector
+    fex = prob.apply_mass_matrix(prob.dtype_u(fex))
+
+    # Apply boundary conditions to both vectors before comparison
+    bc_u = df.DirichletBC(prob.W.sub(0), fu_exact, 'on_boundary')
+    bc_p = df.DirichletBC(prob.W.sub(1), fp_exact, 'on_boundary')
+    bc = [bc_u, bc_p]
+
+    [bc.apply(f.values.vector()) for bc in bc]
+    [bc.apply(fex.values.vector()) for bc in bc]
+
+    # Compute the relative error between computed and exact right-hand sides
+    rel_err = abs(f - fex) / abs(fex)
+    assert rel_err < 5e-5, f"Relative error {rel_err} exceeds tolerance"
+
+
+@pytest.mark.fenics
+def test_problem_class():
+
+    from pySDC.projects.StroemungsRaum.run_Navier_Stokes_equations_monolithic_FEniCS import (
+        setup,
+        run_simulation,
+    )
+
+    t0 = 3.125e-04
+    Tend = 0.01
+    rho = 1
+
+    # Initialize the simulation setup and controller parameters
+    description, controller_params = setup(t0=t0)
+
+    # Run the simulation up to Tend and retrieve the final solution
+    prob, stats, uend = run_simulation(description, controller_params, Tend)
+
+    # Extract velocity and pressure from the mixed solution
+    u, p = df.split(uend.values)
+
+    # Outward normal vector on the obstacle boundary
+    n = -df.FacetNormal(prob.V.mesh())
+
+    # Tangential component of the velocity on the obstacle boundary
+    u_t = df.inner(df.as_vector((n[1], -n[0])), u)
+
+    # Define the drag coefficient functional
+    drag = df.Form(2 / 0.1 * (prob.nu / rho * df.inner(df.grad(u_t), n) * n[1] - p * n[0]) * prob.dsc)
+
+    # Define the lift coefficient functional
+    lift = df.Form(-2 / 0.1 * (prob.nu / rho * df.inner(df.grad(u_t), n) * n[0] + p * n[1]) * prob.dsc)
+
+    # Evaluate drag and lift coefficients
+    dc = df.assemble(drag)
+    lc = df.assemble(lift)
+
+    # Reference drag and lift coefficients at t = 0.0103125
+    # taken from the FEATFLOW benchmark solution
+    rdc = 1.5689678714e-01
+    rlc = -2.2719989899e-04
+
+    # Compute errors in drag and lift coefficients
+    errors = [abs(dc - rdc), abs(lc - rlc)]
+
+    # Assert that the computed drag and lift coefficients are within the specified tolerances
+    assert errors[0] < 5e-2, f"Error in drag coefficient {errors[0]} exceeds tolerance"
+    assert errors[1] < 5e-2, f"Error in lift coefficient {errors[1]} exceeds tolerance"