diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 2c0a5849f4fbb2f0698987602b9aa895dc852221..cf1d3202ceba566957b2f74dd1635d15ba75c251 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -22,15 +22,15 @@
 namespace NumLib
 {
 void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
-    GlobalVector const& x) const
+    GlobalVector const& x, int const process_id) const
 {
-    _equation_system->assemble(x);
+    _equation_system->assemble(x, process_id);
 }
 
 NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
     GlobalVector& x,
-    std::function<void(int, GlobalVector const&)> const&
-        postIterationCallback)
+    std::function<void(int, GlobalVector const&)> const& postIterationCallback,
+    int const process_id)
 {
     namespace LinAlg = MathLib::LinAlg;
     auto& sys = *_equation_system;
@@ -58,7 +58,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
         time_iteration.start();
 
         timer_dirichlet.start();
-        sys.computeKnownSolutions(x_new);
+        sys.computeKnownSolutions(x_new, process_id);
         sys.applyKnownSolutions(x_new);
         time_dirichlet += timer_dirichlet.elapsed();
 
@@ -66,7 +66,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
 
         BaseLib::RunTime time_assembly;
         time_assembly.start();
-        sys.assemble(x_new);
+        sys.assemble(x_new, process_id);
         sys.getA(A);
         sys.getRhs(rhs);
         INFO("[time] Assembly took %g s.", time_assembly.elapsed());
@@ -177,9 +177,9 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
 }
 
 void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
-    GlobalVector const& x) const
+    GlobalVector const& x, int const process_id) const
 {
-    _equation_system->assemble(x);
+    _equation_system->assemble(x, process_id);
     // TODO if the equation system would be reset to nullptr after each
     //      assemble() or solve() call, the user would be forced to set the
     //      equation every time and could not forget it.
@@ -187,8 +187,8 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
 
 NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
     GlobalVector& x,
-    std::function<void(int, GlobalVector const&)> const&
-        postIterationCallback)
+    std::function<void(int, GlobalVector const&)> const& postIterationCallback,
+    int const process_id)
 {
     namespace LinAlg = MathLib::LinAlg;
     auto& sys = *_equation_system;
@@ -220,7 +220,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
         time_iteration.start();
 
         timer_dirichlet.start();
-        sys.computeKnownSolutions(x);
+        sys.computeKnownSolutions(x, process_id);
         sys.applyKnownSolutions(x);
         time_dirichlet += timer_dirichlet.elapsed();
 
@@ -228,7 +228,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
 
         BaseLib::RunTime time_assembly;
         time_assembly.start();
-        sys.assemble(x);
+        sys.assemble(x, process_id);
         sys.getResidual(x, res);
         sys.getJacobian(J);
         INFO("[time] Assembly took %g s.", time_assembly.elapsed());
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index a702cb60947808192d67370a37f9b393ecbb2df1..4ec739bbcc01e09c87a57691d271114e2b83638d 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -44,13 +44,16 @@ public:
      *
      * \param x             the state at which the equation system will be
      *                      assembled.
+     * \param process_id    the process' id which will be assembled.
      */
-    virtual void assemble(GlobalVector const& x) const = 0;
+    virtual void assemble(GlobalVector const& x,
+                          int const process_id) const = 0;
 
     /*! Assemble and solve the equation system.
      *
      * \param x   in: the initial guess, out: the solution.
      * \param postIterationCallback called after each iteration if set.
+     * \param process_id usually used in staggered schemes.
      *
      * \retval true if the equation system could be solved
      * \retval false otherwise
@@ -58,7 +61,8 @@ public:
     virtual NonlinearSolverStatus solve(
         GlobalVector& x,
         std::function<void(int, GlobalVector const&)> const&
-            postIterationCallback) = 0;
+            postIterationCallback,
+        int const process_id) = 0;
 
     virtual ~NonlinearSolverBase() = default;
 };
@@ -106,12 +110,13 @@ public:
         _convergence_criterion = &conv_crit;
     }
 
-    void assemble(GlobalVector const& x) const override;
+    void assemble(GlobalVector const& x, int const process_id) const override;
 
     NonlinearSolverStatus solve(
         GlobalVector& x,
         std::function<void(int, GlobalVector const&)> const&
-            postIterationCallback) override;
+            postIterationCallback,
+        int const process_id) override;
 
 private:
     GlobalLinearSolver& _linear_solver;
@@ -166,12 +171,13 @@ public:
         _convergence_criterion = &conv_crit;
     }
 
-    void assemble(GlobalVector const& x) const override;
+    void assemble(GlobalVector const& x, int const process_id) const override;
 
     NonlinearSolverStatus solve(
         GlobalVector& x,
         std::function<void(int, GlobalVector const&)> const&
-            postIterationCallback) override;
+            postIterationCallback,
+        int const process_id) override;
 
 private:
     GlobalLinearSolver& _linear_solver;
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index c009918599b2e154da017d74ee5cbdfb732902d4..d5367eae3cd81075e24b247a6af0031983de2cc7 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -37,7 +37,7 @@ public:
     //! Assembles the linearized equation system at the point \c x.
     //! The linearized system is \f$A(x) \cdot x = b(x)\f$. Here the matrix
     //! \f$A(x)\f$ and the vector \f$b(x)\f$ are assembled.
-    virtual void assemble(GlobalVector const& x) = 0;
+    virtual void assemble(GlobalVector const& x, int const process_id) = 0;
 
     /*! Writes the residual at point \c x to \c res.
      *
@@ -55,7 +55,8 @@ public:
     virtual void getJacobian(GlobalMatrix& Jac) const = 0;
 
     //! Pre-compute known solutions and possibly store them internally.
-    virtual void computeKnownSolutions(GlobalVector const& x) = 0;
+    virtual void computeKnownSolutions(GlobalVector const& x,
+                                       int const process_id) = 0;
 
     //! Apply known solutions to the solution vector \c x.
     //! \pre computeKnownSolutions() must have been called before.
@@ -82,7 +83,7 @@ public:
     //! Assembles the linearized equation system at the point \c x.
     //! The linearized system is \f$J(x) \cdot \Delta x = (x)\f$. Here the
     //! residual vector \f$r(x)\f$ and its Jacobian \f$J(x)\f$ are assembled.
-    virtual void assemble(GlobalVector const& x) = 0;
+    virtual void assemble(GlobalVector const& x, int const process_id) = 0;
 
     //! Writes the linearized equation system matrix to \c A.
     //! \pre assemble() must have been called before.
@@ -93,7 +94,8 @@ public:
     virtual void getRhs(GlobalVector& rhs) const = 0;
 
     //! Pre-compute known solutions and possibly store them internally.
-    virtual void computeKnownSolutions(GlobalVector const& x) = 0;
+    virtual void computeKnownSolutions(GlobalVector const& x,
+                                       int const process_id) = 0;
 
     //! Apply known solutions to the solution vector \c x.
     //! \pre computeKnownSolutions() must have been called before.
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index 5e20904c26926e23f3ca8d494b7e9a1ea823ae3e..d3c3df2099759569119a94d0bd5dd483f46fd184 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -51,18 +51,18 @@ public:
 
     //! Assemble \c M, \c K and \c b at the provided state (\c t, \c x).
     virtual void assemble(const double t, double const dt,
-                          GlobalVector const& x, GlobalMatrix& M,
-                          GlobalMatrix& K, GlobalVector& b) = 0;
+                          GlobalVector const& x, int const process_id,
+                          GlobalMatrix& M, GlobalMatrix& K,
+                          GlobalVector& b) = 0;
 
     using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
 
     //! Provides known solutions (Dirichlet boundary conditions) vector for
-    //! the ode system at the given time \c t.
+    //! the ode system at the given time \c t and \c process_id.
     virtual std::vector<NumLib::IndexValueVector<Index>> const*
-    getKnownSolutions(double const t, GlobalVector const& x) const
+    getKnownSolutions(double const /*t*/, GlobalVector const& /*x*/,
+                      int const /*process_id*/) const
     {
-        (void)t;
-        (void)x;
         return nullptr;  // by default there are no known solutions
     }
 };
@@ -131,8 +131,9 @@ public:
                                       GlobalVector const& x,
                                       GlobalVector const& xdot,
                                       const double dxdot_dx, const double dx_dx,
-                                      GlobalMatrix& M, GlobalMatrix& K,
-                                      GlobalVector& b, GlobalMatrix& Jac) = 0;
+                                      int const process_id, GlobalMatrix& M,
+                                      GlobalMatrix& K, GlobalVector& b,
+                                      GlobalMatrix& Jac) = 0;
 };
 
 //! @}
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index f8ccbad684fcb96fa9f88c4491612798285d4d7f..d57c682f41219d6809c587f72a682ec3fb388506 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -71,7 +71,8 @@ TimeDiscretizedODESystem<
 
 void TimeDiscretizedODESystem<
     ODESystemTag::FirstOrderImplicitQuasilinear,
-    NonlinearSolverTag::Newton>::assemble(const GlobalVector& x_new_timestep)
+    NonlinearSolverTag::Newton>::assemble(const GlobalVector& x_new_timestep,
+                                          int const process_id)
 {
     namespace LinAlg = MathLib::LinAlg;
 
@@ -90,8 +91,8 @@ void TimeDiscretizedODESystem<
     _Jac->setZero();
 
     _ode.preAssemble(t, dt, x_curr);
-    _ode.assembleWithJacobian(t, dt, x_curr, xdot, dxdot_dx, dx_dx, *_M, *_K,
-                              *_b, *_Jac);
+    _ode.assembleWithJacobian(t, dt, x_curr, xdot, dxdot_dx, dx_dx, process_id,
+                              *_M, *_K, *_b, *_Jac);
 
     LinAlg::finalizeAssembly(*_M);
     LinAlg::finalizeAssembly(*_K);
@@ -126,9 +127,11 @@ void TimeDiscretizedODESystem<
 
 void TimeDiscretizedODESystem<
     ODESystemTag::FirstOrderImplicitQuasilinear,
-    NonlinearSolverTag::Newton>::computeKnownSolutions(GlobalVector const& x)
+    NonlinearSolverTag::Newton>::computeKnownSolutions(GlobalVector const& x,
+                                                       int const process_id)
 {
-    _known_solutions = _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
+    _known_solutions =
+        _ode.getKnownSolutions(_time_disc.getCurrentTime(), x, process_id);
 }
 
 void TimeDiscretizedODESystem<
@@ -187,7 +190,8 @@ TimeDiscretizedODESystem<
 
 void TimeDiscretizedODESystem<
     ODESystemTag::FirstOrderImplicitQuasilinear,
-    NonlinearSolverTag::Picard>::assemble(const GlobalVector& x_new_timestep)
+    NonlinearSolverTag::Picard>::assemble(const GlobalVector& x_new_timestep,
+                                          int const process_id)
 {
     namespace LinAlg = MathLib::LinAlg;
 
@@ -200,7 +204,7 @@ void TimeDiscretizedODESystem<
     _b->setZero();
 
     _ode.preAssemble(t, dt, x_curr);
-    _ode.assemble(t, dt, x_curr, *_M, *_K, *_b);
+    _ode.assemble(t, dt, x_curr, process_id, *_M, *_K, *_b);
 
     LinAlg::finalizeAssembly(*_M);
     LinAlg::finalizeAssembly(*_K);
@@ -209,9 +213,11 @@ void TimeDiscretizedODESystem<
 
 void TimeDiscretizedODESystem<
     ODESystemTag::FirstOrderImplicitQuasilinear,
-    NonlinearSolverTag::Picard>::computeKnownSolutions(GlobalVector const& x)
+    NonlinearSolverTag::Picard>::computeKnownSolutions(GlobalVector const& x,
+                                                       int const process_id)
 {
-    _known_solutions = _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
+    _known_solutions =
+        _ode.getKnownSolutions(_time_disc.getCurrentTime(), x, process_id);
 }
 
 void TimeDiscretizedODESystem<
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 2a38e5f0840a5326174acc09c7749719d92b99e3..9fa14e095fe7000f4ed7cb3e990b2de79317c573 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -83,14 +83,16 @@ public:
 
     ~TimeDiscretizedODESystem() override;
 
-    void assemble(const GlobalVector& x_new_timestep) override;
+    void assemble(const GlobalVector& x_new_timestep,
+                  int const process_id) override;
 
     void getResidual(GlobalVector const& x_new_timestep,
                      GlobalVector& res) const override;
 
     void getJacobian(GlobalMatrix& Jac) const override;
 
-    void computeKnownSolutions(GlobalVector const& x) override;
+    void computeKnownSolutions(GlobalVector const& x,
+                               int const process_id) override;
 
     void applyKnownSolutions(GlobalVector& x) const override;
 
@@ -179,7 +181,8 @@ public:
 
     ~TimeDiscretizedODESystem() override;
 
-    void assemble(const GlobalVector& x_new_timestep) override;
+    void assemble(const GlobalVector& x_new_timestep,
+                  int const process_id) override;
 
     void getA(GlobalMatrix& A) const override
     {
@@ -191,7 +194,8 @@ public:
         _mat_trans->computeRhs(*_M, *_K, *_b, rhs);
     }
 
-    void computeKnownSolutions(GlobalVector const& x) override;
+    void computeKnownSolutions(GlobalVector const& x,
+                               int const process_id) override;
 
     void applyKnownSolutions(GlobalVector& x) const override;
 
diff --git a/ProcessLib/AbstractJacobianAssembler.h b/ProcessLib/AbstractJacobianAssembler.h
index 89a656e86a45480fdd6a638b2ede3b09a0dc84bc..e11cb1917c75945caa7930cad3ba2dc7c8609ea5 100644
--- a/ProcessLib/AbstractJacobianAssembler.h
+++ b/ProcessLib/AbstractJacobianAssembler.h
@@ -40,7 +40,7 @@ public:
         LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
         double const /*dt*/, std::vector<double> const& /*local_xdot*/,
         const double /*dxdot_dx*/, const double /*dx_dx*/,
-        std::vector<double>& /*local_M_data*/,
+        int const /*process_id*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
         std::vector<double>& /*local_b_data*/,
         std::vector<double>& /*local_Jac_data*/,
diff --git a/ProcessLib/AnalyticalJacobianAssembler.cpp b/ProcessLib/AnalyticalJacobianAssembler.cpp
index 97992726372afe2e8d106e610ce9450a82025a0c..1959c2ef151139493faf15ef3ea8516982fd1c06 100644
--- a/ProcessLib/AnalyticalJacobianAssembler.cpp
+++ b/ProcessLib/AnalyticalJacobianAssembler.cpp
@@ -29,14 +29,14 @@ void AnalyticalJacobianAssembler::assembleWithJacobian(
 void AnalyticalJacobianAssembler::assembleWithJacobianForStaggeredScheme(
     LocalAssemblerInterface& local_assembler, double const t, double const dt,
     std::vector<double> const& local_xdot, const double dxdot_dx,
-    const double dx_dx, std::vector<double>& local_M_data,
+    const double dx_dx, int const process_id, std::vector<double>& local_M_data,
     std::vector<double>& local_K_data, std::vector<double>& local_b_data,
     std::vector<double>& local_Jac_data,
     LocalCoupledSolutions const& local_coupled_solutions)
 {
     local_assembler.assembleWithJacobianForStaggeredScheme(
-        t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
-        local_b_data, local_Jac_data, local_coupled_solutions);
+        t, dt, local_xdot, dxdot_dx, dx_dx, process_id, local_M_data,
+        local_K_data, local_b_data, local_Jac_data, local_coupled_solutions);
 }
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/AnalyticalJacobianAssembler.h b/ProcessLib/AnalyticalJacobianAssembler.h
index 181a595547002b5b9a21305754a564e01b6d6a21..814f1166f0937c7d084005b786f9b5878bc70095 100644
--- a/ProcessLib/AnalyticalJacobianAssembler.h
+++ b/ProcessLib/AnalyticalJacobianAssembler.h
@@ -44,7 +44,7 @@ public:
     void assembleWithJacobianForStaggeredScheme(
         LocalAssemblerInterface& local_assembler, double const t,
         double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
+        const double dxdot_dx, const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions) override;
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
index 9401ee461903a3809db04fb8798501ac1d76a389..767e03cc0075db00d4e54460bf8af57161808cd3 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
@@ -18,9 +18,8 @@
 namespace ProcessLib
 {
 CoupledSolutionsForStaggeredScheme::CoupledSolutionsForStaggeredScheme(
-    std::vector<std::reference_wrapper<GlobalVector const>> const& coupled_xs_,
-    const int process_id_)
-    : coupled_xs(coupled_xs_), process_id(process_id_)
+    std::vector<std::reference_wrapper<GlobalVector const>> const& coupled_xs_)
+    : coupled_xs(coupled_xs_)
 {
     for (auto const& coupled_x : coupled_xs)
     {
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.h b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
index 98528d45e3fb94dad8db8309cba9a91c8f0e0507..63f605d0098008ce03b34a9d8298a30dca447bb7 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.h
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
@@ -31,16 +31,13 @@ struct CoupledSolutionsForStaggeredScheme
 {
     CoupledSolutionsForStaggeredScheme(
         std::vector<std::reference_wrapper<GlobalVector const>> const&
-            coupled_xs_,
-        const int process_id_);
+            coupled_xs_);
 
     /// References to the current solutions of the coupled processes.
     std::vector<std::reference_wrapper<GlobalVector const>> const& coupled_xs;
 
     /// Pointers to the vector of the solutions of the previous time step.
     std::vector<GlobalVector*> coupled_xs_t0;
-
-    const int process_id;
 };
 
 /**
@@ -53,17 +50,13 @@ struct CoupledSolutionsForStaggeredScheme
  */
 struct LocalCoupledSolutions
 {
-    LocalCoupledSolutions(const int process_id_,
-                          std::vector<std::vector<double>>&& local_coupled_xs0_,
+    LocalCoupledSolutions(std::vector<std::vector<double>>&& local_coupled_xs0_,
                           std::vector<std::vector<double>>&& local_coupled_xs_)
-        : process_id(process_id_),
-          local_coupled_xs0(std::move(local_coupled_xs0_)),
+        : local_coupled_xs0(std::move(local_coupled_xs0_)),
           local_coupled_xs(std::move(local_coupled_xs_))
     {
     }
 
-    const int process_id;
-
     /// Local solutions of the previous time step.
     std::vector<std::vector<double>> const local_coupled_xs0;
     /// Local solutions of the current time step.
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 8c568dbca9df5a29aa4fe088dcaa5c2dec31b3a0..be82ee86765cb793d403a9724c7a2d61eaea65b2 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -28,7 +28,7 @@ void LocalAssemblerInterface::assemble(double const /*t*/,
 }
 
 void LocalAssemblerInterface::assembleForStaggeredScheme(
-    double const /*t*/, double const /*dt*/,
+    double const /*t*/, double const /*dt*/, int const /*process_id*/,
     std::vector<double>& /*local_M_data*/,
     std::vector<double>& /*local_K_data*/,
     std::vector<double>& /*local_b_data*/,
@@ -56,7 +56,8 @@ void LocalAssemblerInterface::assembleWithJacobian(
 void LocalAssemblerInterface::assembleWithJacobianForStaggeredScheme(
     double const /*t*/, double const /*dt*/,
     std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
-    const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
+    const double /*dx_dx*/, int const /*process_id*/,
+    std::vector<double>& /*local_M_data*/,
     std::vector<double>& /*local_K_data*/,
     std::vector<double>& /*local_b_data*/,
     std::vector<double>& /*local_Jac_data*/,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 7db6e22ef6187026d23f0cb4ba40cdba5f187d91..7d447480c6ff3534dae6d8424c8da076bdf15119 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -53,8 +53,9 @@ public:
                           std::vector<double>& local_b_data);
 
     virtual void assembleForStaggeredScheme(
-        double const t, double const dt, std::vector<double>& local_M_data,
-        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
+        double const t, double const dt, int const process_id,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+        std::vector<double>& local_b_data,
         LocalCoupledSolutions const& coupled_solutions);
 
     virtual void assembleWithJacobian(double const t, double const dt,
@@ -68,7 +69,7 @@ public:
 
     virtual void assembleWithJacobianForStaggeredScheme(
         double const t, double const dt, std::vector<double> const& local_xdot,
-        const double dxdot_dx, const double dx_dx,
+        const double dxdot_dx, const double dx_dx, int const process_id,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
         LocalCoupledSolutions const& local_coupled_solutions);
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index ac52ba1b1c53dde5ffdad3963276f8e9bd8aa52e..7e8815d8209a7acd1951723d1d11e19387dcbdd7 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -196,41 +196,39 @@ void Process::preAssemble(const double t, double const dt,
 }
 
 void Process::assemble(const double t, double const dt, GlobalVector const& x,
-                       GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+                       int const process_id, GlobalMatrix& M, GlobalMatrix& K,
+                       GlobalVector& b)
 {
     MathLib::LinAlg::setLocalAccessibleVector(x);
 
-    assembleConcreteProcess(t, dt, x, M, K, b);
+    assembleConcreteProcess(t, dt, x, process_id, M, K, b);
 
-    const auto pcs_id =
-        (_coupled_solutions) != nullptr ? _coupled_solutions->process_id : 0;
     // the last argument is for the jacobian, nullptr is for a unused jacobian
-    _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b, nullptr);
+    _boundary_conditions[process_id].applyNaturalBC(t, x, K, b, nullptr);
 
     // the last argument is for the jacobian, nullptr is for a unused jacobian
-    _source_term_collections[pcs_id].integrate(t, x, b, nullptr);
+    _source_term_collections[process_id].integrate(t, x, b, nullptr);
 }
 
 void Process::assembleWithJacobian(const double t, double const dt,
                                    GlobalVector const& x,
                                    GlobalVector const& xdot,
                                    const double dxdot_dx, const double dx_dx,
-                                   GlobalMatrix& M, GlobalMatrix& K,
-                                   GlobalVector& b, GlobalMatrix& Jac)
+                                   int const process_id, GlobalMatrix& M,
+                                   GlobalMatrix& K, GlobalVector& b,
+                                   GlobalMatrix& Jac)
 {
     MathLib::LinAlg::setLocalAccessibleVector(x);
     MathLib::LinAlg::setLocalAccessibleVector(xdot);
 
-    assembleWithJacobianConcreteProcess(t, dt, x, xdot, dxdot_dx, dx_dx, M, K,
-                                        b, Jac);
+    assembleWithJacobianConcreteProcess(t, dt, x, xdot, dxdot_dx, dx_dx,
+                                        process_id, M, K, b, Jac);
 
     // TODO: apply BCs to Jacobian.
-    const auto pcs_id =
-        (_coupled_solutions) != nullptr ? _coupled_solutions->process_id : 0;
-    _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b, &Jac);
+    _boundary_conditions[process_id].applyNaturalBC(t, x, K, b, &Jac);
 
     // the last argument is for the jacobian, nullptr is for a unused jacobian
-    _source_term_collections[pcs_id].integrate(t, x, b, &Jac);
+    _source_term_collections[process_id].integrate(t, x, b, &Jac);
 }
 
 void Process::constructDofTable()
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 4e6b7a3ec3f5836e94735f65f0bd222609f33962..addd0df65ad2107b975ceb12082036a48d7a27bf 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -96,24 +96,28 @@ public:
     void updateDeactivatedSubdomains(double const time, const int process_id);
 
     bool isMonolithicSchemeUsed() const { return _use_monolithic_scheme; }
-    virtual void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() {}
+    virtual void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
+        int const /*process_id*/)
+    {
+    }
     void preAssemble(const double t, double const dt,
                      GlobalVector const& x) final;
     void assemble(const double t, double const dt, GlobalVector const& x,
-                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) final;
+                  int const process_id, GlobalMatrix& M, GlobalMatrix& K,
+                  GlobalVector& b) final;
 
     void assembleWithJacobian(const double t, double const dt,
                               GlobalVector const& x, GlobalVector const& xdot,
                               const double dxdot_dx, const double dx_dx,
-                              GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+                              int const process_id, GlobalMatrix& M,
+                              GlobalMatrix& K, GlobalVector& b,
                               GlobalMatrix& Jac) final;
 
     std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
-    getKnownSolutions(double const t, GlobalVector const& x) const final
+    getKnownSolutions(double const t, GlobalVector const& x,
+                      int const process_id) const final
     {
-        const auto pcs_id =
-            (_coupled_solutions) ? _coupled_solutions->process_id : 0;
-        return _boundary_conditions[pcs_id].getKnownSolutions(t, x);
+        return _boundary_conditions[process_id].getKnownSolutions(t, x);
     }
 
     virtual NumLib::LocalToGlobalIndexMap const& getDOFTable(
@@ -192,13 +196,14 @@ private:
     }
 
     virtual void assembleConcreteProcess(const double t, double const dt,
-                                         GlobalVector const& x, GlobalMatrix& M,
+                                         GlobalVector const& x,
+                                         int const process_id, GlobalMatrix& M,
                                          GlobalMatrix& K, GlobalVector& b) = 0;
 
     virtual void assembleWithJacobianConcreteProcess(
         const double t, double const dt, GlobalVector const& x,
         GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
-        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) = 0;
 
     virtual void preTimestepConcreteProcess(GlobalVector const& /*x*/,
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 065f4e42616baebdd04266338129e0f485f92c58..822adbc1777f4e7b980a2d3d1bcbbcdeac56c69c 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -157,7 +157,7 @@ std::vector<GlobalVector*> setInitialConditions(
 
             auto const nl_tag = process_data->nonlinear_solver_tag;
             setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
-            nonlinear_solver.assemble(x0);
+            nonlinear_solver.assemble(x0, process_id);
             time_disc.pushState(
                 t0, x0,
                 mat_strg);  // TODO: that might do duplicate work
@@ -173,6 +173,7 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
     Output& output_control)
 {
     auto& process = process_data.process;
+    int const process_id = process_data.process_id;
     auto& time_disc = *process_data.time_disc;
     auto& conv_crit = *process_data.conv_crit;
     auto& ode_sys = *process_data.tdisc_ode_sys;
@@ -190,16 +191,16 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
 
     auto const post_iteration_callback = [&](int iteration,
                                              GlobalVector const& x) {
-        output_control.doOutputNonlinearIteration(
-            process, process_data.process_id, timestep, t, x, iteration);
+        output_control.doOutputNonlinearIteration(process, process_id, timestep,
+                                                  t, x, iteration);
     };
 
     auto const nonlinear_solver_status =
-        nonlinear_solver.solve(x, post_iteration_callback);
+        nonlinear_solver.solve(x, post_iteration_callback, process_id);
 
     if (nonlinear_solver_status.error_norms_met)
     {
-        process.postNonLinearSolver(x, t, delta_t, process_data.process_id);
+        process.postNonLinearSolver(x, t, delta_t, process_id);
     }
 
     return nonlinear_solver_status;
@@ -574,7 +575,7 @@ void postTimestepForAllProcesses(
         if (is_staggered_coupling)
         {
             CoupledSolutionsForStaggeredScheme coupled_solutions(
-                solutions_of_coupled_processes, process_id);
+                solutions_of_coupled_processes);
             pcs.setCoupledSolutionsForStaggeredScheme(&coupled_solutions);
         }
         auto& x = *_process_solutions[process_id];
@@ -662,7 +663,7 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
             auto& x = *_process_solutions[process_id];
 
             CoupledSolutionsForStaggeredScheme coupled_solutions(
-                _solutions_of_coupled_processes, process_id);
+                _solutions_of_coupled_processes);
 
             process_data->process.setCoupledSolutionsForStaggeredScheme(
                 &coupled_solutions);
@@ -787,12 +788,13 @@ void TimeLoop::outputSolutions(bool const output_initial_condition,
         if (is_staggered_coupling)
         {
             CoupledSolutionsForStaggeredScheme coupled_solutions(
-                _solutions_of_coupled_processes, process_id);
+                _solutions_of_coupled_processes);
 
             process_data->process.setCoupledSolutionsForStaggeredScheme(
                 &coupled_solutions);
             process_data->process
-                .setCoupledTermForTheStaggeredSchemeToLocalAssemblers();
+                .setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
+                    process_id);
         }
         (output_object.*output_class_member)(pcs, process_id, timestep, t, x);
     }
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index c4b039614e6f483450bfc365ba6c5346f4ab9e08..68268ff6096c3ab5f3eb21e723216c4bf6612716 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -43,8 +43,8 @@ void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
         dof_tables,
-    const double t, double const dt, const GlobalVector& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b,
+    const double t, double const dt, const GlobalVector& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
     CoupledSolutionsForStaggeredScheme const* const cpl_xs)
 {
     std::vector<std::vector<GlobalIndexType>> indices_of_processes;
@@ -55,9 +55,7 @@ void VectorMatrixAssembler::assemble(
             NumLib::getIndices(mesh_item_id, dof_table.get()));
     }
 
-    auto const& indices = (cpl_xs == nullptr)
-                              ? indices_of_processes[0]
-                              : indices_of_processes[cpl_xs->process_id];
+    auto const& indices = indices_of_processes[process_id];
     _local_M_data.clear();
     _local_K_data.clear();
     _local_b_data.clear();
@@ -76,12 +74,11 @@ void VectorMatrixAssembler::assemble(
             getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
 
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-            cpl_xs->process_id, std::move(local_coupled_xs0),
-            std::move(local_coupled_xs));
+            std::move(local_coupled_xs0), std::move(local_coupled_xs));
 
-        local_assembler.assembleForStaggeredScheme(t, dt, _local_M_data,
-                                                   _local_K_data, _local_b_data,
-                                                   local_coupled_solutions);
+        local_assembler.assembleForStaggeredScheme(
+            t, dt, process_id, _local_M_data, _local_K_data, _local_b_data,
+            local_coupled_solutions);
     }
 
     auto const num_r_c = indices.size();
@@ -111,8 +108,8 @@ void VectorMatrixAssembler::assembleWithJacobian(
         dof_tables,
     const double t, double const dt, GlobalVector const& x,
     GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
-    GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
-    CoupledSolutionsForStaggeredScheme const* const cpl_xs)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac, CoupledSolutionsForStaggeredScheme const* const cpl_xs)
 {
     std::vector<std::vector<GlobalIndexType>> indices_of_processes;
     indices_of_processes.reserve(dof_tables.size());
@@ -122,9 +119,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
             NumLib::getIndices(mesh_item_id, dof_table.get()));
     }
 
-    auto const& indices = (cpl_xs == nullptr)
-                              ? indices_of_processes[0]
-                              : indices_of_processes[cpl_xs->process_id];
+    auto const& indices = indices_of_processes[process_id];
     auto const local_xdot = xdot.get(indices);
 
     _local_M_data.clear();
@@ -147,12 +142,11 @@ void VectorMatrixAssembler::assembleWithJacobian(
             getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
 
         ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-            cpl_xs->process_id, std::move(local_coupled_xs0),
-            std::move(local_coupled_xs));
+            std::move(local_coupled_xs0), std::move(local_coupled_xs));
 
         _jacobian_assembler->assembleWithJacobianForStaggeredScheme(
-            local_assembler, t, dt, local_xdot, dxdot_dx, dx_dx, _local_M_data,
-            _local_K_data, _local_b_data, _local_Jac_data,
+            local_assembler, t, dt, local_xdot, dxdot_dx, dx_dx, process_id,
+            _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
             local_coupled_solutions);
     }
 
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index 1e79299c34f114cb5d2722ef580d2bc9c074f70b..d0a450d71979098d5cde5ce7b04bbc00591984f2 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -48,7 +48,8 @@ public:
                   std::vector<std::reference_wrapper<
                       NumLib::LocalToGlobalIndexMap>> const& dof_tables,
                   double const t, double const dt, GlobalVector const& x,
-                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+                  int const process_id, GlobalMatrix& M, GlobalMatrix& K,
+                  GlobalVector& b,
                   CoupledSolutionsForStaggeredScheme const* const cpl_xs);
 
     //! Assembles \c M, \c K, \c b, and the Jacobian \c Jac of the residual.
@@ -61,7 +62,8 @@ public:
             dof_tables,
         const double t, double const dt, GlobalVector const& x,
         GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
-        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac,
         CoupledSolutionsForStaggeredScheme const* const cpl_xs);
 
 private: