diff --git a/pySDC/projects/StroemungsRaum/README.rst b/pySDC/projects/StroemungsRaum/README.rst
new file mode 100644
index 0000000000..c66176ca33
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/README.rst
@@ -0,0 +1,45 @@
+StroemungsRaum
+==============
+
+**StroemungsRaum** is a research software project developed within the
+BMBF-funded project
+
+*“StrömungsRaum – Novel Exascale Architectures with Heterogeneous Hardware
+Components for Computational Fluid Dynamics Simulations”*
+(October 2022 – September 2025).
+
+The project addresses the development of scalable numerical methods and
+high-performance algorithms for Computational Fluid Dynamics (CFD) targeting
+future **exascale computing architectures** with heterogeneous hardware.
+
+Scope of This Repository
+------------------------
+This repository contains the **Forschungszentrum Jülich (FZJ)** contribution to
+the StrömungsRaum project, focusing on:
+
+- Parallel-in-time methods
+- Combined space–time parallelization for fluid simulations
+- Algorithmic scalability for time-dependent PDEs
+
+The goal is to expose concurrency beyond spatial parallelism and enable
+efficient execution on large-scale HPC systems.
+
+Model Problems and Methods
+--------------------------
+Implemented examples and test cases include:
+
+- Heat equation
+- Convection–diffusion and nonlinear convection–diffusion problems
+- Incompressible Navier–Stokes equations, using:  
+   - Projection methods
+   - Monolithic formulations
+   - DAE- and PDE sweepers
+
+These serve as benchmarks and demonstrators for scalable space–time CFD
+simulations.
+
+Funding
+-------
+Funded by the **German Federal Ministry of Education and Research (BMBF)** under
+grant number **16ME0708**.
+
diff --git a/pySDC/projects/StroemungsRaum/environment.yml b/pySDC/projects/StroemungsRaum/environment.yml
new file mode 100644
index 0000000000..aa1baeb636
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/environment.yml
@@ -0,0 +1,16 @@
+---
+
+name: pySDC
+channels:
+  - conda-forge
+dependencies:
+  - fenics>=2019.1.0
+  - numpy
+  - scipy>=0.17.1
+  - dill>=0.2.6
+  - matplotlib>=3.0
+  - pytest
+  - pytest-cov
+  - pip
+  - pip:
+    - qmat>=0.1.8
diff --git a/pySDC/projects/StroemungsRaum/problem_classes/HeatEquation_2D_FEniCS.py b/pySDC/projects/StroemungsRaum/problem_classes/HeatEquation_2D_FEniCS.py
new file mode 100644
index 0000000000..0a57d1743e
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/problem_classes/HeatEquation_2D_FEniCS.py
@@ -0,0 +1,240 @@
+import logging
+
+import dolfin as df
+import numpy as np
+
+from pySDC.core.problem import Problem
+from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh
+
+
+class fenics_heat2D_mass(Problem):
+    r"""
+    Example implementing the forced two-dimensional heat equation with Dirichlet boundary conditions
+
+    .. math::
+        \frac{\partial u}{\partial t} = \nu \Delta u + f
+
+    for :math:`x \in \Omega:=[0,1]x[0,1]`, where the forcing term :math:`f` is defined by
+
+    .. math::
+        f(x, y, t) = -\sin(\pi x)\sin(\pi y) (\sin(t) - 2\nu \pi^2 \cos(t)).
+
+    For initial conditions with constant c and
+
+    .. math::
+        u(x, y, 0) = \sin(\pi x)\sin(\pi y) + c
+
+    the exact solution of the problem is given by
+
+    .. math::
+        u(x, y, t) = \sin(\pi x)\sin(\pi y)\cos(t) + c.
+
+    In this class the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem is reformulated to
+    the *weak formulation*
+
+    .. math:
+        \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.
+
+    The part containing the forcing term is treated explicitly, where it is interpolated in the function
+    space. The other part will be treated in an implicit way.
+
+    Parameters
+    ----------
+    c_nvars : int, optional
+        Spatial resolution.
+    t0 : float, optional
+        Starting time.
+    family : str, optional
+        Indicates the family of elements used to create the function space
+        for the trail and test functions. The default is ``'CG'``, which are the class
+        of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
+    order : int, optional
+        Defines the order of the elements in the function space.
+    nu : float, optional
+        Diffusion coefficient :math:`\nu`.
+    c: float, optional
+        Constant for the Dirichlet boundary condition :math: `c`
+
+    Attributes
+    ----------
+    V : FunctionSpace
+        Defines the function space of the trial and test functions.
+    M : scalar, vector, matrix or higher rank tensor
+        Denotes the expression :math:`\int_\Omega u_t v\,dx`.
+    K : scalar, vector, matrix or higher rank tensor
+        Denotes the expression :math:`- \nu \int_\Omega \nabla u \nabla v\,dx`.
+    g : Expression
+        The forcing term :math:`f` in the heat equation.
+    bc : DirichletBC
+        Denotes the Dirichlet 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 = rhs_fenics_mesh
+
+    def __init__(self, c_nvars=64, t0=0.0, family='CG', order=2, nu=0.1, c=0.0):
+
+        # define the Dirichlet boundary
+        def Boundary(x, on_boundary):
+            return on_boundary
+
+        # set solver and form parameters
+        df.parameters["form_compiler"]["optimize"] = True
+        df.parameters["form_compiler"]["cpp_optimize"] = True
+        df.parameters['allow_extrapolation'] = True
+
+        # set mesh and refinement (for multilevel)
+        mesh = df.UnitSquareMesh(c_nvars, c_nvars)
+
+        # define function space for future reference
+        self.V = df.FunctionSpace(mesh, family, order)
+        tmp = df.Function(self.V)
+        self.logger.debug('DoFs on this level:', len(tmp.vector()[:]))
+
+        # invoke super init, passing number of dofs, dtype_u and dtype_f
+        super().__init__(self.V)
+        self._makeAttributeAndRegister('c_nvars', 't0', 'family', 'order', 'nu', 'c', localVars=locals(), readOnly=True)
+
+        # Stiffness term (Laplace)
+        u = df.TrialFunction(self.V)
+        v = df.TestFunction(self.V)
+        a_K = -1.0 * df.inner(df.nabla_grad(u), self.nu * df.nabla_grad(v)) * df.dx
+
+        # Mass term
+        a_M = u * v * df.dx
+
+        self.M = df.assemble(a_M)
+        self.K = df.assemble(a_K)
+
+        # set boundary values
+        self.u_D = df.Expression('sin(a*x[0]) * sin(a*x[1]) * cos(t) + c', c=self.c, a=np.pi, t=t0, degree=self.order)
+        self.bc = df.DirichletBC(self.V, self.u_D, Boundary)
+        self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)
+
+        self.fix_bc_for_residual = True
+
+        # set forcing term as expression
+        self.g = df.Expression(
+            '-sin(a*x[0]) * sin(a*x[1]) * (sin(t) - 2*b*a*a*cos(t))',
+            a=np.pi,
+            b=self.nu,
+            t=self.t0,
+            degree=self.order,
+        )
+
+    def solve_system(self, rhs, factor, u0, t):
+        r"""
+        Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`.
+
+        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
+        -------
+        u : dtype_u
+            Solution.
+        """
+
+        u = self.dtype_u(u0)
+        T = self.M - factor * self.K
+        b = self.dtype_u(rhs)
+
+        self.u_D.t = t
+
+        self.bc.apply(T, b.values.vector())
+        self.bc.apply(b.values.vector())
+
+        df.solve(T, u.values.vector(), b.values.vector())
+
+        return u
+
+    def eval_f(self, u, t):
+        """
+        Routine to evaluate both parts of the right-hand side.
+
+        Parameters
+        ----------
+        u : 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 divided into two parts.
+        """
+
+        f = self.dtype_f(self.V)
+
+        self.K.mult(u.values.vector(), f.impl.values.vector())
+
+        self.g.t = t
+        f.expl = self.apply_mass_matrix(self.dtype_u(df.interpolate(self.g, self.V)))
+        return f
+
+    def apply_mass_matrix(self, u):
+        r"""
+        Routine to apply mass matrix.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Current values of the numerical solution.
+
+        Returns
+        -------
+        me : dtype_u
+            The product :math:`M \vec{u}`.
+        """
+
+        me = self.dtype_u(self.V)
+        self.M.mult(u.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.
+        """
+        u0 = df.Expression('sin(a*x[0]) * sin(a*x[1]) * cos(t) + c', c=self.c, a=np.pi, t=t, degree=self.order)
+        me = self.dtype_u(df.interpolate(u0, self.V), val=self.V)
+
+        return me
+
+    def fix_residual(self, res):
+        """
+        Applies homogeneous Dirichlet boundary conditions to the residual
+
+        Parameters
+        ----------
+        res : dtype_u
+              Residual
+        """
+        self.bc_hom.apply(res.values.vector())
+        return None
diff --git a/pySDC/projects/StroemungsRaum/run_heat_equation_FEniCS.py b/pySDC/projects/StroemungsRaum/run_heat_equation_FEniCS.py
new file mode 100644
index 0000000000..96144797d8
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/run_heat_equation_FEniCS.py
@@ -0,0 +1,175 @@
+import json
+from pathlib import Path
+import dolfin as df
+
+from pySDC.projects.StroemungsRaum.problem_classes.HeatEquation_2D_FEniCS import fenics_heat2D_mass
+from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass
+
+from pySDC.helpers.stats_helper import get_sorted
+from pySDC.implementations.hooks.log_solution import LogSolution
+
+
+def setup(t0=None):
+    """
+    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 = 0.2
+
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1e-12
+    level_params['dt'] = dt
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 20
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['quad_type'] = 'RADAU-RIGHT'
+    sweeper_params['num_nodes'] = 2
+
+    problem_params = dict()
+    problem_params['nu'] = 0.1
+    problem_params['t0'] = t0
+    problem_params['c_nvars'] = 64
+    problem_params['family'] = 'CG'
+    problem_params['order'] = 2
+    problem_params['c'] = 0.0
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['logger_level'] = 20
+    controller_params['hook_class'] = LogSolution
+
+    description = dict()
+
+    description['problem_class'] = fenics_heat2D_mass
+    description['sweeper_class'] = imex_1st_order_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 heat equation problem.
+
+     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,
+        rel_err: float,
+           relative final-time error.
+    """
+    # 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)
+
+    # get exact solution at final time for error calculation
+    uex = P.u_exact(Tend)
+
+    # call main function to get things done...
+    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+
+    # compute relative error at final time
+    rel_err = abs(uex - uend) / abs(uex)
+
+    return P, stats, rel_err
+
+
+def run_postprocessing(description, problem, stats, Tend):
+    """
+    Postprocess and store simulation results for visualization and analysis.
+
+    Args:
+        description: dict,
+            pySDC description containing problem parameters.
+        problem: Problem instance,
+            Problem instance containing the final solution and other problem-related information.
+        stats: dict,
+            collected runtime statistics,
+        Tend: float,
+            Final simulation time.
+
+    Returns: None
+    """
+    # Get the data directory
+    import os
+
+    path = f"{os.path.dirname(__file__)}/data/heat_equation/"
+
+    # If it does not exist, create the 'data' directory at the specified path, including any necessary parent directories
+    Path(path).mkdir(parents=True, exist_ok=True)
+
+    # Save parameters
+    parameters = description['problem_params']
+    parameters.update(description['level_params'])
+    parameters['Tend'] = Tend
+    json.dump(parameters, open(path + "heat_equation_FEniCS_parameters.json", 'w'))
+
+    # Create XDMF file for visualization output
+    xdmffile_u = df.XDMFFile(path + "heat_equation_FEniCS_Temperature.xdmf")
+
+    # Get the solution at every time step, sorted by time
+    Solutions = get_sorted(stats, type='u', sortby='time')
+
+    for i in range(len(Solutions)):
+        time = Solutions[i][0]
+        #
+        un = Solutions[i][1]
+        ux = problem.u_exact(time)
+        #
+        xdmffile_u.write_checkpoint(un.values, "un", time, df.XDMFFile.Encoding.HDF5, True)
+        xdmffile_u.write_checkpoint(ux.values, "ux", time, df.XDMFFile.Encoding.HDF5, True)
+        #
+    xdmffile_u.close()
+
+    return None
+
+
+if __name__ == "__main__":
+
+    t0 = 0.0
+    Tend = 1.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 relative error
+    problem, stats, rel_err = run_simulation(description, controller_params, Tend)
+    print('The relative error at time ', Tend, 'is ', rel_err)
+
+    # run postprocessing to save parameters and solution for visualization
+    run_postprocessing(description, problem, stats, Tend)
diff --git a/pySDC/projects/StroemungsRaum/tests/test_heat_equation_FEniCS.py b/pySDC/projects/StroemungsRaum/tests/test_heat_equation_FEniCS.py
new file mode 100644
index 0000000000..d6edb38f3e
--- /dev/null
+++ b/pySDC/projects/StroemungsRaum/tests/test_heat_equation_FEniCS.py
@@ -0,0 +1,24 @@
+import pytest
+
+
+@pytest.mark.fenics
+@pytest.mark.mpi4py
+def test_problem_class():
+    """
+    This test checks the functionality of the problem class for the heat equation implemented in FEniCS.
+    It runs a short simulation and checks if the relative error at the final time is below a certain threshold,
+    indicating that the problem class is correctly implemented and can be used for time integration.
+    """
+
+    from pySDC.projects.StroemungsRaum.run_heat_equation_FEniCS import setup, run_simulation
+
+    t0 = 0.0
+    Tend = 1.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 relative error
+    _, _, rel_err = run_simulation(description, controller_params, Tend)
+
+    assert rel_err <= 2e-4, 'ERROR: Relative error is too high, got rel_err = {}'.format(rel_err)