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/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 91cfac19b9da137e83bac9bf71c04122d76ad325..3b1d21a97944a30188892c84c4a8f9f4371e5cf8 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -398,14 +398,12 @@ public:
     }
 
     void assembleForStaggeredScheme(
-        double const t,
-        double const dt,
-        std::vector<double>& local_M_data,
-        std::vector<double>& local_K_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_xs) override
     {
-        if (coupled_xs.process_id == hydraulic_process_id)
+        if (process_id == hydraulic_process_id)
         {
             assembleHydraulicEquation(t, dt, local_M_data, local_K_data,
                                       local_b_data, coupled_xs);
@@ -413,8 +411,9 @@ public:
         else
         {
             // Go for assembling in an order of transport process id.
-            assembleComponentTransportEquation(
-                t, dt, local_M_data, local_K_data, local_b_data, coupled_xs);
+            assembleComponentTransportEquation(t, dt, local_M_data,
+                                               local_K_data, local_b_data,
+                                               coupled_xs, process_id);
         }
     }
 
@@ -534,9 +533,8 @@ public:
         double const t, double const dt, std::vector<double>& local_M_data,
         std::vector<double>& local_K_data,
         std::vector<double>& /*local_b_data*/,
-        LocalCoupledSolutions const& coupled_xs)
+        LocalCoupledSolutions const& coupled_xs, int const transport_process_id)
     {
-        auto const& transport_process_id = coupled_xs.process_id;
         auto local_C = Eigen::Map<const NodalVectorType>(
             coupled_xs.local_coupled_xs[transport_process_id].data(),
             concentration_size);
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 034ff36795279e4f03dad7f9ea34110020114932..b0776379699329ebb4bbe76fb1315f549af88d8f 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -90,12 +90,11 @@ void ComponentTransportProcess::initializeConcreteProcess(
 }
 
 void ComponentTransportProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble ComponentTransportProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -114,7 +113,7 @@ void ComponentTransportProcess::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_tables, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_tables, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
@@ -134,11 +133,11 @@ void ComponentTransportProcess::setCoupledSolutionsOfPreviousTimeStep()
 void ComponentTransportProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian ComponentTransportProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
@@ -146,7 +145,7 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 Eigen::Vector3d ComponentTransportProcess::getFlux(std::size_t const element_id,
@@ -175,12 +174,10 @@ Eigen::Vector3d ComponentTransportProcess::getFlux(std::size_t const element_id,
 }
 
 void ComponentTransportProcess::
-    setCoupledTermForTheStaggeredSchemeToLocalAssemblers()
+    setCoupledTermForTheStaggeredSchemeToLocalAssemblers(int const process_id)
 {
     DBUG("Set the coupled term for the staggered scheme to local assembers.");
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &ComponentTransportLocalAssemblerInterface::
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 65400f9dc5f2a2b4d9c08bfc20809cc45a4402b8..675d4d265984c7d842c34eeb2ead188aded0b05c 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -123,7 +123,8 @@ public:
         return _process_id_to_component_name_map;
     }
 
-    void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() override;
+    void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
+        int const process_id) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, const double /*t*/,
                                     const double /*delta_t*/,
@@ -141,13 +142,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void setCoupledSolutionsOfPreviousTimeStep();
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/CreateProcessData.cpp b/ProcessLib/CreateProcessData.cpp
index fcf42bafd3ca21c8ea4306624a2138c832ae35d0..4f46cae448f2593ebcf75a016966015269dabdbe 100644
--- a/ProcessLib/CreateProcessData.cpp
+++ b/ProcessLib/CreateProcessData.cpp
@@ -20,6 +20,7 @@ namespace ProcessLib
 static std::unique_ptr<ProcessData> makeProcessData(
     std::unique_ptr<NumLib::TimeStepAlgorithm>&& timestepper,
     NumLib::NonlinearSolverBase& nonlinear_solver,
+    int const process_id,
     Process& process,
     std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
     std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit)
@@ -32,7 +33,7 @@ static std::unique_ptr<ProcessData> makeProcessData(
     {
         return std::make_unique<ProcessData>(
             std::move(timestepper), *nonlinear_solver_picard,
-            std::move(conv_crit), std::move(time_disc), process);
+            std::move(conv_crit), std::move(time_disc), process_id, process);
     }
     if (auto* nonlinear_solver_newton =
             dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
@@ -40,7 +41,7 @@ static std::unique_ptr<ProcessData> makeProcessData(
     {
         return std::make_unique<ProcessData>(
             std::move(timestepper), *nonlinear_solver_newton,
-            std::move(conv_crit), std::move(time_disc), process);
+            std::move(conv_crit), std::move(time_disc), process_id, process);
     }
 
     OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
@@ -53,6 +54,7 @@ std::vector<std::unique_ptr<ProcessData>> createPerProcessData(
         nonlinear_solvers)
 {
     std::vector<std::unique_ptr<ProcessData>> per_process_data;
+    int process_id = 0;
 
     //! \ogs_file_param{prj__time_loop__processes__process}
     for (auto pcs_config : config.getConfigSubtreeList("process"))
@@ -99,9 +101,10 @@ std::vector<std::unique_ptr<ProcessData>> createPerProcessData(
                 "in the current project file!");
         }
 
-        per_process_data.emplace_back(makeProcessData(
-            std::move(timestepper), nl_slv, pcs, std::move(time_disc),
-            std::move(conv_crit)));
+        per_process_data.emplace_back(
+            makeProcessData(std::move(timestepper), nl_slv, process_id, pcs,
+                            std::move(time_disc), std::move(conv_crit)));
+        ++process_id;
     }
 
     if (per_process_data.size() != processes.size())
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 939bf17269637bc59a57df0e550ca5b4b6a439ac..907d031c385eef77a8c5279389bd9445feb38fd1 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -58,30 +58,29 @@ void GroundwaterFlowProcess::initializeConcreteProcess(
 }
 
 void GroundwaterFlowProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble GroundwaterFlowProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void GroundwaterFlowProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian GroundwaterFlowProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
@@ -89,7 +88,7 @@ void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 }   // namespace GroundwaterFlow
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index adc1cc69b4bb31e938bc723923e4a49ca0705bc8..b77a0b6b2f3b4537cbea8b03626ba4f46e2be26c 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -91,13 +91,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     GroundwaterFlowProcessData _process_data;
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 25b86f5eab65a2afe3cceaf680532e7a1a952a8d..ef9fb38cd41233a9798dab6fd12a9204092dae1a 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -87,7 +87,8 @@ void HTProcess::initializeConcreteProcess(
 }
 
 void HTProcess::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)
 {
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -99,7 +100,7 @@ void HTProcess::assembleConcreteProcess(const double t, double const dt,
     }
     else
     {
-        if (_coupled_solutions->process_id == _heat_transport_process_id)
+        if (process_id == _heat_transport_process_id)
         {
             DBUG(
                 "Assemble the equations of heat transport process within "
@@ -116,20 +117,19 @@ void HTProcess::assembleConcreteProcess(const double t, double const dt,
         dof_tables.emplace_back(*_local_to_global_index_map);
     }
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_tables, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_tables, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void HTProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian HTProcess.");
 
@@ -147,13 +147,11 @@ void HTProcess::assembleWithJacobianConcreteProcess(
     }
 
     // Call global assembler for each local assembly item.
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
@@ -183,12 +181,11 @@ void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
     MathLib::LinAlg::setLocalAccessibleVector(x0);
 }
 
-void HTProcess::setCoupledTermForTheStaggeredSchemeToLocalAssemblers()
+void HTProcess::setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
+    int const process_id)
 {
     DBUG("Set the coupled term for the staggered scheme to local assembers.");
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &HTLocalAssemblerInterface::setStaggeredCoupledSolutions,
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 6ddf0c1f9070e70bfc59db456dd9eb52d0d35512..67c0905fe9f6ed10f7b5301d45ceb32a673e7fbf 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -78,7 +78,8 @@ public:
                             double const t,
                             GlobalVector const& x) const override;
 
-    void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() override;
+    void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
+        int const process_id) override;
 
     void postTimestepConcreteProcess(GlobalVector const& x,
                                      const double t,
@@ -92,13 +93,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 1b294c099dd76d3a23a7c4a936b2e7a2b35dcbf2..b5f4865a33f54ddaa38c71b1e0813df8e167d1a4 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -25,12 +25,13 @@ template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
     assembleForStaggeredScheme(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_xs)
 {
-    if (coupled_xs.process_id == _heat_transport_process_id)
+    if (process_id == _heat_transport_process_id)
     {
         assembleHeatTransportEquation(t, local_M_data, local_K_data,
                                       local_b_data, coupled_xs);
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index 13b5b68153f5cf25c25b05bba09202382ec55a61..b1aa5d15fcb169a8f8c8676e10e14321bbcd25d1 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -67,8 +67,9 @@ public:
     }
 
     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_xs) override;
 
     std::vector<double> const& getIntPtDarcyVelocity(
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index b7fe8621c49349632570f76531b567042f916827..a1afe805bba16ebe71451e796a1b7d0f20d7139b 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -73,12 +73,11 @@ void HeatConductionProcess::initializeConcreteProcess(
 }
 
 void HeatConductionProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble HeatConductionProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -86,18 +85,18 @@ void HeatConductionProcess::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void HeatConductionProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian HeatConductionProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -106,7 +105,7 @@ void HeatConductionProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 void HeatConductionProcess::computeSecondaryVariableConcrete(
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index 7126dac0f80d7fb079babb9244a1847c93258c36..5a226946a1899e49b795e868e1c0803518b77a2d 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -51,13 +51,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const /*dt*/,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     HeatConductionProcessData _process_data;
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
index df1dcdc5adcc00dc6ce2f41fb4b76b9aa4cae37c..e4374b836b2a6e1e4ceda604eef5dba640748804 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
@@ -159,12 +159,11 @@ void HeatTransportBHEProcess::initializeConcreteProcess(
 }
 
 void HeatTransportBHEProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble HeatTransportBHE process.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -172,15 +171,15 @@ void HeatTransportBHEProcess::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void HeatTransportBHEProcess::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*/, GlobalMatrix& /*Jac*/)
+    const double /*dx_dx*/, int const /*process_id*/, GlobalMatrix& /*M*/,
+    GlobalMatrix& /*K*/, GlobalVector& /*b*/, GlobalMatrix& /*Jac*/)
 {
     OGS_FATAL(
         "HeatTransportBHE: analytical Jacobian assembly is not implemented");
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
index c01b5eecc757dfa273fa24cfe64b5eed14123b6f..ad26dc925acf41fce17eb55bdba98c4eeebf39b4 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
@@ -55,13 +55,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void createBHEBoundaryConditionTopBottom(
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 88c7b894b19fddf7caca490939611a822f90d6ca..24ec4d958bef31054f0bb9f982454d95b257e6be 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -551,13 +551,13 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                   DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
         const double t, double const dt, const std::vector<double>& 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,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
     // For the equations with pressure
-    if (local_coupled_solutions.process_id == 0)
+    if (process_id == 0)
     {
         assembleWithJacobianForPressureEquations(
             t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index d7390cfa2dee03b3c288431c95bf8fb8bbaa78ff..1afed2c4c5ab63f24a43a9d554547b2332516ae1 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -152,7 +152,7 @@ public:
 
     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) override;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
index f73c76de6e1cb30cb84aef35444dddd09413f6a8..462d6af8198f4a4c1921da50cb14e173f96f269e 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
@@ -278,8 +278,8 @@ void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
 
 template <int DisplacementDim>
 void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble the equations for HydroMechanics");
 
@@ -292,7 +292,7 @@ void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        dof_table, t, dt, x, M, K, b, _coupled_solutions);
+        dof_table, t, dt, x, process_id, M, K, b, _coupled_solutions);
 }
 
 template <int DisplacementDim>
@@ -300,7 +300,8 @@ void HydroMechanicsProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_tables;
@@ -315,7 +316,7 @@ void HydroMechanicsProcess<DisplacementDim>::
     else
     {
         // For the staggered scheme
-        if (_coupled_solutions->process_id == 0)
+        if (process_id == 0)
         {
             DBUG(
                 "Assemble the Jacobian equations of liquid fluid process in "
@@ -331,14 +332,12 @@ void HydroMechanicsProcess<DisplacementDim>::
         dof_tables.emplace_back(*_local_to_global_index_map);
     }
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         if (_use_monolithic_scheme)
@@ -349,16 +348,16 @@ void HydroMechanicsProcess<DisplacementDim>::
         }
         else
         {
-            transformVariableFromGlobalVector(
-                b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
-                std::negate<double>());
+            transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
+                                              output_vector,
+                                              std::negate<double>());
         }
     };
-    if (_use_monolithic_scheme || _coupled_solutions->process_id == 0)
+    if (_use_monolithic_scheme || process_id == 0)
     {
         copyRhs(0, *_hydraulic_flow);
     }
-    if (_use_monolithic_scheme || _coupled_solutions->process_id == 1)
+    if (_use_monolithic_scheme || process_id == 1)
     {
         copyRhs(1, *_nodal_forces);
     }
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index d9007cef85e32c8aae7b6bd2e96bfd98bff70fc5..f47eebb68741c33caecd14a583f70d513ffc2720 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -73,13 +73,14 @@ private:
     void initializeBoundaryConditions() override;
 
     void assembleConcreteProcess(const double t, double const /*dt*/,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 8dd224bf345adf1f7a3dc90be2e93894d653e652..5185fd833994fcbfa764145d7d34e786e308861f 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -566,8 +566,8 @@ bool HydroMechanicsProcess<GlobalDim>::isLinear() const
 
 template <int GlobalDim>
 void HydroMechanicsProcess<GlobalDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble HydroMechanicsProcess.");
 
@@ -576,19 +576,18 @@ void HydroMechanicsProcess<GlobalDim>::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        dof_table, t, dt, x, M, K, b, _coupled_solutions);
+        dof_table, t, dt, x, process_id, M, K, b, _coupled_solutions);
 }
 
 template <int GlobalDim>
 void HydroMechanicsProcess<GlobalDim>::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian HydroMechanicsProcess.");
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
@@ -597,7 +596,7 @@ void HydroMechanicsProcess<GlobalDim>::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         transformVariableFromGlobalVector(b, variable_id,
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index f80940e78d4fe4821bac707becca9dd769d17226..068b854886293ce03f778019868a943fc7015212 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -66,13 +66,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
                                     double const dt,
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index bf6fa8e9939ff4228dc07c8cde6467918485a831..18a9aff117a17e698d06c9841e63d0fa07fe867f 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -542,12 +542,11 @@ bool SmallDeformationProcess<DisplacementDim>::isLinear() const
 
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble SmallDeformationProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -555,7 +554,7 @@ void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 template <int DisplacementDim>
@@ -563,12 +562,12 @@ void SmallDeformationProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
     // Call global assembler for each local assembly item.
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -576,7 +575,7 @@ void SmallDeformationProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::preTimestepConcreteProcess(
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
index f594ef1d52a264d841912327a74728a6d4d2b3a9..41355af7d0b728a0f3b31de50c4991a3376bc618 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
@@ -62,13 +62,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 5c62713e170a13d681926d241c8e050adf9268b0..80009b47e2c10df9460f1b611658b3b36557023b 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -78,45 +78,41 @@ void LiquidFlowProcess::initializeConcreteProcess(
             &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
-void LiquidFlowProcess::assembleConcreteProcess(const double t, double const dt,
-                                                GlobalVector const& x,
-                                                GlobalMatrix& M,
-                                                GlobalMatrix& K,
-                                                GlobalVector& b)
+void LiquidFlowProcess::assembleConcreteProcess(
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble LiquidFlowProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
 
-
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void LiquidFlowProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian LiquidFlowProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index 930e75a6a0c0fedb892ef5698601bb2644ffcbe1..64b66ac3e6e13a6361d7ad9b66d7c0e895f3fe63 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -90,13 +90,14 @@ private:
         MeshLib::Mesh const& mesh, unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     const int _gravitational_axis_id;
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/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 58828fc80b1722b262fcf8cf2217b37540313b12..5098caa42eac07306e43f22fdf3fbbe641a197f4 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -23,13 +23,13 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                               DisplacementDim>::
     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)
 {
     // For the equations with phase field.
-    if (local_coupled_solutions.process_id == 1)
+    if (process_id == 1)
     {
         assembleWithJacobianPhaseFiledEquations(
             t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index 3c1509d5b42cdc8063f1a697adb54b47923eea9b..9727fd3346cc9c50574c0cea406e0dcdfa2d0123 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -204,7 +204,7 @@ public:
 
     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) override;
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.cpp b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
index 4ce9fbc645ba66e44616d05bc38cddb5fdf7d887..3fb0e865729f704403c7f0d1b68b7810ed5b83eb 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
@@ -162,8 +162,8 @@ void PhaseFieldProcess<DisplacementDim>::initializeBoundaryConditions()
 
 template <int DisplacementDim>
 void PhaseFieldProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble PhaseFieldProcess.");
 
@@ -171,7 +171,7 @@ void PhaseFieldProcess<DisplacementDim>::assembleConcreteProcess(
         dof_tables;
 
     // For the staggered scheme
-    if (_coupled_solutions->process_id == 1)
+    if (process_id == 1)
     {
         DBUG(
             "Assemble the equations of phase field in "
@@ -186,13 +186,12 @@ void PhaseFieldProcess<DisplacementDim>::assembleConcreteProcess(
     dof_tables.emplace_back(*_local_to_global_index_map_single_component);
     dof_tables.emplace_back(*_local_to_global_index_map);
 
-    const int process_id = _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_tables, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_tables, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
@@ -200,13 +199,14 @@ template <int DisplacementDim>
 void PhaseFieldProcess<DisplacementDim>::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_tables;
 
     // For the staggered scheme
-    if (_coupled_solutions->process_id == 1)
+    if (process_id == 1)
     {
         DBUG(
             "Assemble the Jacobian equations of phase field in "
@@ -221,16 +221,15 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
     dof_tables.emplace_back(*_local_to_global_index_map);
     dof_tables.emplace_back(*_local_to_global_index_map_single_component);
 
-    const int process_id = _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 
-    if (_coupled_solutions->process_id == 0)
+    if (process_id == 0)
     {
         b.copyValues(*_nodal_forces);
         std::transform(_nodal_forces->begin(), _nodal_forces->end(),
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.h b/ProcessLib/PhaseField/PhaseFieldProcess.h
index 2607f2cd22e29c32c25a26f741471b499c8ec72d..066502e3c6ce556a2f6c92d416aa9b01fbab561c 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.h
@@ -88,13 +88,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
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/ProcessData.h b/ProcessLib/ProcessData.h
index 69415ff11aab1310bacf69bf0e8e48da252bbfcd..ad26f9e75471477e5106661a1000034aa3411619 100644
--- a/ProcessLib/ProcessData.h
+++ b/ProcessLib/ProcessData.h
@@ -27,6 +27,7 @@ struct ProcessData
                 NumLib::NonlinearSolver<NLTag>& nonlinear_solver,
                 std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
                 std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
+                int const process_id_,
                 Process& process_)
         : timestepper(std::move(timestepper_)),
           nonlinear_solver_tag(NLTag),
@@ -34,6 +35,7 @@ struct ProcessData
           nonlinear_solver_status{true, 0},
           conv_crit(std::move(conv_crit_)),
           time_disc(std::move(time_disc_)),
+          process_id(process_id_),
           process(process_)
     {
     }
@@ -47,6 +49,7 @@ struct ProcessData
           time_disc(std::move(pd.time_disc)),
           tdisc_ode_sys(std::move(pd.tdisc_ode_sys)),
           mat_strg(pd.mat_strg),
+          process_id(pd.process_id),
           process(pd.process)
     {
         pd.mat_strg = nullptr;
@@ -67,9 +70,7 @@ struct ProcessData
     //! cast of \c tdisc_ode_sys to NumLib::InternalMatrixStorage
     NumLib::InternalMatrixStorage* mat_strg = nullptr;
 
-    /// Process ID. It is alway 0 when the monolithic scheme is used or
-    /// a single process is modelled.
-    int process_id = 0;
+    int const process_id;
 
     Process& process;
 };
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
index f84286059f709f1bd783cbc1757b0edf91d5a649..ec3302b75a17a5a020323dc0b86fcecdff9ef35f 100644
--- a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
@@ -66,42 +66,39 @@ void RichardsComponentTransportProcess::initializeConcreteProcess(
 }
 
 void RichardsComponentTransportProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble RichardsComponentTransportProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void RichardsComponentTransportProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian RichardsComponentTransportProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 }  // namespace RichardsComponentTransport
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
index 50fb81583c8f7efae576068f56da0b993f6ecd2f..49f659f491dddaa6f548326c9df78fef3d034646 100644
--- a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
@@ -130,13 +130,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     RichardsComponentTransportProcessData _process_data;
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index 1300c4ceeddc0a8dacf45d32c36fc8005b914bfd..b8a534f5fdca9e973da66559246168c485833d90 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -66,40 +66,39 @@ void RichardsFlowProcess::initializeConcreteProcess(
 }
 
 void RichardsFlowProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble RichardsFlowProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void RichardsFlowProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian RichardsFlowProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 }  // namespace RichardsFlow
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.h b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
index e0e26b4e0c33f2c363e10c0f96416db7cd48471b..3f2695f0cc0c7443216e7ec80c93935c331ac843 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
@@ -54,13 +54,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     RichardsFlowProcessData _process_data;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 578dc75aaf4998b3bbbe42b0de73f6db6874a7a9..3c7a31f6c252f8a36a0ae52a9252d33616489e23 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -770,13 +770,13 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                      DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
         const double t, double const dt, const std::vector<double>& 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,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
     // For the equations with pressure
-    if (local_coupled_solutions.process_id == 0)
+    if (process_id == 0)
     {
         assembleWithJacobianForPressureEquations(
             t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index aa9038f8c10c97fe3fc1a673ddd604cd4f5f8172..3623ef73ed9312e63c8e7f599498efe3af7d05f7 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -85,7 +85,7 @@ public:
 
     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) override;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index 38c9b3a009376ce9f3a2c9411b77dc80818411b7..acadc4d3673e68955cd2fe13ed49528721f2bf60 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -234,8 +234,8 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
 
 template <int DisplacementDim>
 void RichardsMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble the equations for RichardsMechanics");
 
@@ -245,14 +245,12 @@ void RichardsMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
@@ -261,7 +259,8 @@ void RichardsMechanicsProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_tables;
@@ -276,7 +275,7 @@ void RichardsMechanicsProcess<DisplacementDim>::
     else
     {
         // For the staggered scheme
-        if (_coupled_solutions->process_id == 0)
+        if (process_id == 0)
         {
             DBUG(
                 "Assemble the Jacobian equations of liquid fluid process in "
@@ -292,14 +291,12 @@ void RichardsMechanicsProcess<DisplacementDim>::
         dof_tables.emplace_back(*_local_to_global_index_map);
     }
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         if (_use_monolithic_scheme)
@@ -310,16 +307,16 @@ void RichardsMechanicsProcess<DisplacementDim>::
         }
         else
         {
-            transformVariableFromGlobalVector(
-                b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
-                std::negate<double>());
+            transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
+                                              output_vector,
+                                              std::negate<double>());
         }
     };
-    if (_use_monolithic_scheme || _coupled_solutions->process_id == 0)
+    if (_use_monolithic_scheme || process_id == 0)
     {
         copyRhs(0, *_hydraulic_flow);
     }
-    if (_use_monolithic_scheme || _coupled_solutions->process_id == 1)
+    if (_use_monolithic_scheme || process_id == 1)
     {
         copyRhs(1, *_nodal_forces);
     }
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
index 85dabe58b786b241b67d3738039f70f459fda6ec..eaa4e943008508e758e55c978265ecd9f860f668 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
@@ -75,13 +75,14 @@ private:
     void initializeBoundaryConditions() override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
index 44fb8146999f4397228ff061c003df32b5e675c5..bda7c7f63dbc0d8e1b2877d5b6f853c9b8351ba1 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
@@ -241,12 +241,11 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
 
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble SmallDeformationProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -254,7 +253,7 @@ void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
@@ -263,11 +262,11 @@ void SmallDeformationProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian SmallDeformationProcess.");
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
@@ -276,7 +275,7 @@ void SmallDeformationProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 
     transformVariableFromGlobalVector(b, 0, *_local_to_global_index_map,
                                       *_nodal_forces, std::negate<double>());
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index 89d5c35540fab66113e464b5e19eee82e8947849..1b54dcd3a48bbad76d62a8c2ee63af32d36ecfac 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -51,13 +51,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void postTimestepConcreteProcess(GlobalVector const& x, const double t,
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
index 28e61f0e163846b9c3a45c9e427b244b12872120..8b01241720a3cc51a7891ab6f489079736a16357 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
@@ -231,21 +231,20 @@ void SmallDeformationNonlocalProcess<DisplacementDim>::
 
 template <int DisplacementDim>
 void SmallDeformationNonlocalProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble SmallDeformationNonlocalProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
@@ -272,21 +271,21 @@ void SmallDeformationNonlocalProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian SmallDeformationNonlocalProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
 
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 
     b.copyValues(*_nodal_forces);
     std::transform(_nodal_forces->begin(), _nodal_forces->end(),
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
index fc8c32f444e7563dfa12e503bb6031414652b926..100f3bac6e7fb54a5d1d985854d73a91d079b198 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
@@ -54,8 +54,9 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     void preAssembleConcreteProcess(const double t, double const dt,
                                     GlobalVector const& x) override;
@@ -63,7 +64,7 @@ private:
     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) override;
 
     void postTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 11d7c85b0b48606495f162200288cfe138eab17c..47a9c9841d6b99861cccbe604e50a5b0e8401121 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -229,40 +229,38 @@ void TESProcess::initializeSecondaryVariables()
 }
 
 void TESProcess::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)
 {
     DBUG("Assemble TESProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void TESProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 void TESProcess::preTimestepConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index 7da1d1ef6dd72fa7e3938bfccce2b80cd12c0903..765ade2fe1da8c6ba7f189a56e68af32f5e57efa 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -60,15 +60,16 @@ private:
         MeshLib::Mesh const& mesh, unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     void initializeSecondaryVariables();
 
     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) override;
 
     GlobalVector const& computeVapourPartialPressure(
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
index 0831bb465bcaaf5313ce2c462da981d0dffcaa01..b5ccf3584c46470a16bf0931920352176eb9bbdb 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
@@ -75,40 +75,39 @@ void ThermalTwoPhaseFlowWithPPProcess::initializeConcreteProcess(
 }
 
 void ThermalTwoPhaseFlowWithPPProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble ThermalTwoPhaseFlowWithPPProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void ThermalTwoPhaseFlowWithPPProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian ThermalTwoPhaseFlowWithPPProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 void ThermalTwoPhaseFlowWithPPProcess::preTimestepConcreteProcess(
     GlobalVector const& x, double const t, double const delta_t,
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h
index d641b584335b9afae57d1c57fc25db78e049f4b8..0087c6921a45f8b9a575edca51c17ead7f91b90c 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h
@@ -62,13 +62,14 @@ private:
         MeshLib::Mesh const& mesh, unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, const double t,
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
index e5e16a1f6bf7cc9f5ee3f0ea806bfe176d9d882f..2f577aacdf593e95624086113ebafdcb9de56985 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
@@ -240,8 +240,8 @@ void ThermoHydroMechanicsProcess<
 
 template <int DisplacementDim>
 void ThermoHydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble the equations for ThermoHydroMechanics");
 
@@ -254,7 +254,7 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        dof_table, t, dt, x, M, K, b, _coupled_solutions);
+        dof_table, t, dt, x, process_id, M, K, b, _coupled_solutions);
 }
 
 template <int DisplacementDim>
@@ -262,7 +262,8 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_tables;
@@ -277,13 +278,13 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::
     else
     {
         // For the staggered scheme
-        if (_coupled_solutions->process_id == 0)
+        if (process_id == 0)
         {
             DBUG(
                 "Assemble the Jacobian equations of heat transport process in "
                 "ThermoHydroMechanics for the staggered scheme.");
         }
-        else if (_coupled_solutions->process_id == 1)
+        else if (process_id == 1)
         {
             DBUG(
                 "Assemble the Jacobian equations of liquid fluid process in "
@@ -302,8 +303,8 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::
 
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
-        _local_assemblers, dof_tables, t, dt, x, xdot, dxdot_dx, dx_dx, M, K, b,
-        Jac, _coupled_solutions);
+        _local_assemblers, dof_tables, t, dt, x, xdot, dxdot_dx, dx_dx,
+        process_id, M, K, b, Jac, _coupled_solutions);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         if (_use_monolithic_scheme)
@@ -314,16 +315,16 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::
         }
         else
         {
-            transformVariableFromGlobalVector(
-                b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
-                std::negate<double>());
+            transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
+                                              output_vector,
+                                              std::negate<double>());
         }
     };
-    if (_use_monolithic_scheme || _coupled_solutions->process_id == 1)
+    if (_use_monolithic_scheme || process_id == 1)
     {
         copyRhs(0, *_hydraulic_flow);
     }
-    if (_use_monolithic_scheme || _coupled_solutions->process_id == 2)
+    if (_use_monolithic_scheme || process_id == 2)
     {
         copyRhs(1, *_nodal_forces);
     }
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
index ad6694deff53883401df7f534e3779faea7ad564..da51d994b9b10a80f62cd4c19afb9444310b2a47 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
@@ -76,13 +76,14 @@ private:
     void initializeBoundaryConditions() override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
index ae7fd38fd97a84c4fad703bfc15c45f8ba6f5d6a..5e2cdf6b12a51a70e3de79cfb95abeae997c0530 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
@@ -25,12 +25,12 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
                                               DisplacementDim>::
     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)
 {
-    if (local_coupled_solutions.process_id == _phase_field_process_id)
+    if (process_id == _phase_field_process_id)
     {
         assembleWithJacobianForPhaseFieldEquations(
             t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
@@ -38,7 +38,7 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         return;
     }
 
-    if (local_coupled_solutions.process_id == _heat_conduction_process_id)
+    if (process_id == _heat_conduction_process_id)
     {
         assembleWithJacobianForHeatConductionEquations(
             t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
index 54104600cde410b1037ff7a6dd8f7f381d3141e6..fff7bf8ecafb28171437546d242b1489d4cdd246 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
@@ -217,7 +217,7 @@ public:
 
     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) override;
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
index fc535f1ec0cc07bd94e523dcac3488d83e653619..2e3ad20164dbe7c9fd0a3b2455e2b6f3849db269 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
@@ -192,6 +192,7 @@ template <int DisplacementDim>
 void ThermoMechanicalPhaseFieldProcess<
     DisplacementDim>::assembleConcreteProcess(const double t, double const dt,
                                               GlobalVector const& x,
+                                              int const process_id,
                                               GlobalMatrix& M, GlobalMatrix& K,
                                               GlobalVector& b)
 {
@@ -199,14 +200,12 @@ void ThermoMechanicalPhaseFieldProcess<
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
@@ -215,12 +214,13 @@ void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_tables;
     // For the staggered scheme
-    if (_coupled_solutions->process_id == _mechanics_related_process_id)
+    if (process_id == _mechanics_related_process_id)
     {
         DBUG(
             "Assemble the Jacobian equations of "
@@ -229,7 +229,7 @@ void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
             "the staggered scheme.");
     }
 
-    if (_coupled_solutions->process_id == _phase_field_process_id)
+    if (process_id == _phase_field_process_id)
     {
         DBUG(
             "Assemble the Jacobian equations of"
@@ -251,14 +251,12 @@ void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
         getDOFTableByProcessID(_mechanics_related_process_id));
     dof_tables.emplace_back(getDOFTableByProcessID(_phase_field_process_id));
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h
index ea6f22f1c5d9ce5ad3ae8c6d53fd8648b13fedc9..c18fd2ae246164a589e06757b774ca3b65763d30 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h
@@ -93,13 +93,14 @@ private:
         unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
index f1cdd5b6327fd1443042f72dfb20eb9d53f07a07..ac1af0754c992afd31085af45e7e1a015d074678 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
@@ -306,14 +306,13 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
                                    DisplacementDim>::
     assembleWithJacobianForStaggeredScheme(
         const double t, double const dt, const std::vector<double>& 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,
         const LocalCoupledSolutions& local_coupled_solutions)
 {
     // For the equations with pressure
-    if (local_coupled_solutions.process_id ==
-        _process_data.heat_conduction_process_id)
+    if (process_id == _process_data.heat_conduction_process_id)
     {
         assembleWithJacobianForHeatConductionEquations(
             t, dt, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index f0408e6c98361d11b128ba6b46ec44bfeb056e4c..dd4397d31c9d3d413ec71465953be79dea5fdf78 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -135,7 +135,7 @@ public:
 
     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) override;
@@ -143,9 +143,9 @@ public:
     void assembleWithJacobian(double const t, double const dt,
                               std::vector<double> const& local_x,
                               std::vector<double> const& local_xdot,
-                              const double /*dxdot_dx*/, const double /*dx_dx*/,
-                              std::vector<double>& /*local_M_data*/,
-                              std::vector<double>& /*local_K_data*/,
+                              const double dxdot_dx, const double dx_dx,
+                              std::vector<double>& local_M_data,
+                              std::vector<double>& local_K_data,
                               std::vector<double>& local_rhs_data,
                               std::vector<double>& local_Jac_data) override;
 
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
index cf37e037d657aad853c067bf8b3d154130b30332..6c149fd75f89b182451f7f1775ffd59c91060874 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -277,21 +277,19 @@ void ThermoMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
 
 template <int DisplacementDim>
 void ThermoMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble ThermoMechanicsProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
         dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
@@ -300,7 +298,8 @@ void ThermoMechanicsProcess<DisplacementDim>::
     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, GlobalMatrix& Jac)
+        int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+        GlobalMatrix& Jac)
 {
     DBUG("AssembleJacobian ThermoMechanicsProcess.");
 
@@ -317,8 +316,7 @@ void ThermoMechanicsProcess<DisplacementDim>::
     else
     {
         // For the staggered scheme
-        if (_coupled_solutions->process_id ==
-            _process_data.heat_conduction_process_id)
+        if (process_id == _process_data.heat_conduction_process_id)
         {
             DBUG(
                 "Assemble the Jacobian equations of heat conduction process in "
@@ -349,14 +347,12 @@ void ThermoMechanicsProcess<DisplacementDim>::
         setCoupledSolutionsOfPreviousTimeStep();
     }
 
-    const int process_id =
-        _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 
     // TODO (naumov): Refactor the copy rhs part. This is copy from HM.
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
@@ -368,19 +364,18 @@ void ThermoMechanicsProcess<DisplacementDim>::
         }
         else
         {
-            transformVariableFromGlobalVector(
-                b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
-                std::negate<double>());
+            transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
+                                              output_vector,
+                                              std::negate<double>());
         }
     };
     if (_use_monolithic_scheme ||
-        _coupled_solutions->process_id ==
-            _process_data.heat_conduction_process_id)
+        process_id == _process_data.heat_conduction_process_id)
     {
         copyRhs(0, *_heat_flux);
     }
     if (_use_monolithic_scheme ||
-        _coupled_solutions->process_id == _process_data.mechanics_process_id)
+        process_id == _process_data.mechanics_process_id)
     {
         copyRhs(1, *_nodal_forces);
     }
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
index 344488ddd60244293877fa95e15d0158e7129a91..24e76adbdaa646f2fd4f19827bb57e5b111dde7f 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -68,13 +68,14 @@ private:
     void initializeBoundaryConditions() override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, double const t,
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 33fe66b1204d579623a1bed51ee925e12750f899..822adbc1777f4e7b980a2d3d1bcbbcdeac56c69c 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -130,10 +130,10 @@ std::vector<GlobalVector*> setInitialConditions(
 {
     std::vector<GlobalVector*> process_solutions;
 
-    int process_id = 0;
     for (auto& process_data : per_process_data)
     {
         auto& pcs = process_data->process;
+        auto const process_id = process_data->process_id;
         auto& time_disc = *process_data->time_disc;
 
         auto& ode_sys = *process_data->tdisc_ode_sys;
@@ -157,24 +157,23 @@ 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
         }
-
-        ++process_id;
     }
 
     return process_solutions;
 }
 
 NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
-    int const process_id, GlobalVector& x, std::size_t const timestep,
-    double const t, double const delta_t, ProcessData const& process_data,
+    GlobalVector& x, std::size_t const timestep, double const t,
+    double const delta_t, ProcessData const& process_data,
     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;
@@ -197,7 +196,7 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
     };
 
     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)
     {
@@ -228,10 +227,9 @@ TimeLoop::TimeLoop(
 void TimeLoop::setCoupledSolutions()
 {
     _solutions_of_coupled_processes.reserve(_per_process_data.size());
-    for (unsigned process_id = 0; process_id < _per_process_data.size();
-         process_id++)
+    for (auto& process_data : _per_process_data)
     {
-        auto const& x = *_process_solutions[process_id];
+        auto const& x = *_process_solutions[process_data->process_id];
         _solutions_of_coupled_processes.emplace_back(x);
 
         // Create a vector to store the solution of the last coupling iteration
@@ -383,13 +381,12 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
 /// initialize output, convergence criterion, etc.
 void TimeLoop::initialize()
 {
-    int process_id = 0;
     for (auto& process_data : _per_process_data)
     {
         auto& pcs = process_data->process;
+        int const process_id = process_data->process_id;
         _output->addProcess(pcs, process_id);
 
-        process_data->process_id = process_id;
         setTimeDiscretizedODESystem(*process_data);
 
         if (auto* conv_crit =
@@ -405,8 +402,6 @@ void TimeLoop::initialize()
         // the fixed times.
         auto& timestepper = process_data->timestepper;
         timestepper->addFixedOutputTimes(_output->getFixedOutputTimes());
-
-        ++process_id;
     }
 
     // init solution storage
@@ -465,11 +460,10 @@ bool TimeLoop::loop()
              timesteps, t, dt);
 
         // Check element deactivation:
-        int process_id = 0;
         for (auto& process_data : _per_process_data)
         {
-            process_data->process.updateDeactivatedSubdomains(t, process_id);
-            ++process_id;
+            process_data->process.updateDeactivatedSubdomains(
+                t, process_data->process_id);
         }
 
         if (is_staggered_coupling)
@@ -549,8 +543,8 @@ static NumLib::NonlinearSolverStatus solveMonolithicProcess(
     BaseLib::RunTime time_timestep_process;
     time_timestep_process.start();
 
-    auto const nonlinear_solver_status = solveOneTimeStepOneProcess(
-        process_data.process_id, x, timestep_id, t, dt, process_data, output);
+    auto const nonlinear_solver_status =
+        solveOneTimeStepOneProcess(x, timestep_id, t, dt, process_data, output);
 
     INFO("[time] Solving process #%u took %g s in time step #%u ",
          process_data.process_id, time_timestep_process.elapsed(), timestep_id);
@@ -581,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];
@@ -598,24 +592,24 @@ NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
     NumLib::NonlinearSolverStatus nonlinear_solver_status;
     for (auto& process_data : _per_process_data)
     {
-        nonlinear_solver_status = solveMonolithicProcess(
-            t, dt, timestep_id, *process_data,
-            *_process_solutions[process_data->process_id], *_output);
+        auto const process_id = process_data->process_id;
+        nonlinear_solver_status =
+            solveMonolithicProcess(t, dt, timestep_id, *process_data,
+                                   *_process_solutions[process_id], *_output);
 
         process_data->nonlinear_solver_status = nonlinear_solver_status;
         if (!nonlinear_solver_status.error_norms_met)
         {
             ERR("The nonlinear solver failed in time step #%u at t = %g s for "
                 "process #%u.",
-                timestep_id, t, process_data->process_id);
+                timestep_id, t, process_id);
 
             if (!process_data->timestepper->isSolutionErrorComputationNeeded())
             {
                 // save unsuccessful solution
-                _output->doOutputAlways(
-                    process_data->process, process_data->process_id,
-                    timestep_id, t,
-                    *_process_solutions[process_data->process_id]);
+                _output->doOutputAlways(process_data->process, process_id,
+                                        timestep_id, t,
+                                        *_process_solutions[process_id]);
                 OGS_FATAL(nonlinear_fixed_dt_fails_info.data());
             }
 
@@ -659,23 +653,23 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
     {
         // TODO(wenqing): use process name
         coupling_iteration_converged = true;
-        int process_id = 0;
         int const last_process_id = _per_process_data.size() - 1;
         for (auto& process_data : _per_process_data)
         {
+            auto const process_id = process_data->process_id;
             BaseLib::RunTime time_timestep_process;
             time_timestep_process.start();
 
             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);
 
             nonlinear_solver_status = solveOneTimeStepOneProcess(
-                process_id, x, timestep_id, t, dt, *process_data, *_output);
+                x, timestep_id, t, dt, *process_data, *_output);
             process_data->nonlinear_solver_status = nonlinear_solver_status;
 
             INFO(
@@ -719,8 +713,6 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
                 }
             }
             MathLib::LinAlg::copy(x, x_old);
-
-            ++process_id;
         }  // end of for (auto& process_data : _per_process_data)
 
         if (coupling_iteration_converged && global_coupling_iteration > 0)
@@ -771,15 +763,14 @@ void TimeLoop::outputSolutions(bool const output_initial_condition,
     bool const is_staggered_coupling =
         !isMonolithicProcess(*_per_process_data[0]);
 
-    unsigned process_id = 0;
     for (auto& process_data : _per_process_data)
     {
+        auto const process_id = process_data->process_id;
         auto& pcs = process_data->process;
         // If nonlinear solver diverged, the solution has already been
         // saved.
         if (!process_data->nonlinear_solver_status.error_norms_met)
         {
-            ++process_id;
             continue;
         }
 
@@ -797,16 +788,15 @@ 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);
-
-        ++process_id;
     }
 }
 
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 49512b1366e0c8c9b7a53b5ba914af04daa8b345..bf5bf0f4c20dae8cd40a4c3520c65250832e977f 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -73,40 +73,39 @@ void TwoPhaseFlowWithPPProcess::initializeConcreteProcess(
 }
 
 void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble TwoPhaseFlowWithPPProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void TwoPhaseFlowWithPPProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPPProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
 }  // namespace TwoPhaseFlowWithPP
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
index acf35353d2ae3b47bd339c66fff01fd47eb22074..83a7355615bd7d791c52876e7a4928fe1c6587ed 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
@@ -61,13 +61,14 @@ private:
         MeshLib::Mesh const& mesh, unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     TwoPhaseFlowWithPPProcessData _process_data;
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
index 376b1e45bc9d7554d4c732b6451b3bfcaa513297..2c1ccb4ecfb51b68f2558ae18a7d389c9069a675 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -73,40 +73,39 @@ void TwoPhaseFlowWithPrhoProcess::initializeConcreteProcess(
 }
 
 void TwoPhaseFlowWithPrhoProcess::assembleConcreteProcess(
-    const double t, double const dt, GlobalVector const& x, GlobalMatrix& M,
-    GlobalMatrix& K, GlobalVector& b)
+    const double t, double const dt, GlobalVector const& x,
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     DBUG("Assemble TwoPhaseFlowWithPrhoProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
-        pv.getActiveElementIDs(), dof_table, t, dt, x, M, K, b,
+        pv.getActiveElementIDs(), dof_table, t, dt, x, process_id, M, K, b,
         _coupled_solutions);
 }
 
 void TwoPhaseFlowWithPrhoProcess::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, GlobalMatrix& Jac)
+    int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+    GlobalMatrix& Jac)
 {
     DBUG("AssembleWithJacobian TwoPhaseFlowWithPrhoProcess.");
 
     std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
        dof_table = {std::ref(*_local_to_global_index_map)};
-    const int process_id = 0;
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
     GlobalVector const& x, double const t, double const dt,
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
index a0547fd48f118a9e1032fa280547671125fa0cb1..bab3f17d6f1dd9a9ec4e78e87c86a284be3b5398 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
@@ -58,13 +58,14 @@ private:
         MeshLib::Mesh const& mesh, unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, double const dt,
-                                 GlobalVector const& x, GlobalMatrix& M,
-                                 GlobalMatrix& K, GlobalVector& b) override;
+                                 GlobalVector const& x, int const process_id,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
 
     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) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, const double t,
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:
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index 17c3072b1d1f023d93bfbd4ace7c61eac6c514e2..28da0e408ed07605532453c4b09ab697af17f72b 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -32,8 +32,8 @@ public:
     }
 
     void assemble(const double /*t*/, double const /*dt*/,
-                  GlobalVector const& /*x*/, GlobalMatrix& M, GlobalMatrix& K,
-                  GlobalVector& b) override
+                  GlobalVector const& /*x*/, int const /*process_id*/,
+                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
     {
         MathLib::setMatrix(M, { 1.0, 0.0,  0.0, 1.0 });
         MathLib::setMatrix(K, { 0.0, 1.0, -1.0, 0.0 });
@@ -45,12 +45,13 @@ public:
                               GlobalVector const& x_curr,
                               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) override
     {
         namespace LinAlg = MathLib::LinAlg;
 
-        assemble(t, dt, x_curr, M, K, b);
+        assemble(t, dt, x_curr, process_id, M, K, b);
 
         // compute Jac = M*dxdot_dx + dx_dx*K
         LinAlg::finalizeAssembly(M);
@@ -117,8 +118,8 @@ public:
     }
 
     void assemble(const double /*t*/, double const /*dt*/,
-                  GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
-                  GlobalVector& b) override
+                  GlobalVector const& x, int const /*process_id*/,
+                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
     {
         MathLib::setMatrix(M, {1.0});
         MathLib::LinAlg::setLocalAccessibleVector(x);
@@ -130,10 +131,11 @@ public:
                               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) override
     {
-        assemble(t, dt, x, M, K, b);
+        assemble(t, dt, x, process_id, M, K, b);
 
         namespace LinAlg = MathLib::LinAlg;
 
@@ -210,8 +212,8 @@ public:
     }
 
     void assemble(const double /*t*/, double const /*dt*/,
-                  GlobalVector const& x_curr, GlobalMatrix& M, GlobalMatrix& K,
-                  GlobalVector& b) override
+                  GlobalVector const& x_curr, int const /*process_id*/,
+                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
     {
         MathLib::LinAlg::setLocalAccessibleVector(x_curr);
         auto const u = x_curr[0];
@@ -226,13 +228,13 @@ public:
     void assembleWithJacobian(const double t, double const dt,
                               GlobalVector const& x_curr,
                               GlobalVector const& xdot, const double dxdot_dx,
-                              const double dx_dx, GlobalMatrix& M,
-                              GlobalMatrix& K, GlobalVector& b,
+                              const double dx_dx, int const process_id,
+                              GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
                               GlobalMatrix& Jac) override
     {
         MathLib::LinAlg::setLocalAccessibleVector(x_curr);
         MathLib::LinAlg::setLocalAccessibleVector(xdot);
-        assemble(t, dt, x_curr, M, K, b);
+        assemble(t, dt, x_curr, process_id, M, K, b);
 
         auto const u = x_curr[0];
         auto const v = x_curr[1];
diff --git a/Tests/NumLib/TimeLoopSingleODE.h b/Tests/NumLib/TimeLoopSingleODE.h
index 55bac5c8f669b95a187c2cffa4cc7e0e37e23f60..c86a9aab4e054d28c1a58d765af8ab57eb4fe387 100644
--- a/Tests/NumLib/TimeLoopSingleODE.h
+++ b/Tests/NumLib/TimeLoopSingleODE.h
@@ -96,8 +96,9 @@ NumLib::NonlinearSolverStatus TimeLoopSingleODE<NLTag>::loop(
 
     if (time_disc.needsPreload())
     {
+        int const process_id = 0;
         MathLib::LinAlg::setLocalAccessibleVector(x);
-        _nonlinear_solver->assemble(x);
+        _nonlinear_solver->assemble(x, process_id);
         time_disc.pushState(t0, x0,
                             _ode_sys);  // TODO: that might do duplicate work
     }
@@ -113,7 +114,9 @@ NumLib::NonlinearSolverStatus TimeLoopSingleODE<NLTag>::loop(
         // INFO("time: %e, delta_t: %e", t, delta_t);
         time_disc.nextTimestep(t, delta_t);
 
-        nonlinear_solver_status = _nonlinear_solver->solve(x, nullptr);
+        int const process_id = 0;
+        nonlinear_solver_status =
+            _nonlinear_solver->solve(x, nullptr, process_id);
         if (!nonlinear_solver_status.error_norms_met)
         {
             break;