diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index 11e916d309c2f363840a9dd4f64b4701961250cc..460a0bd06146143fdf8f421d2a61940f0397de9e 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -54,7 +54,9 @@ public:
         // Therefore there is nothing to do here.
     }
 
-    virtual void preTimestep(const double /*t*/, GlobalVector const& /*x*/)
+    virtual void preTimestep(const double /*t*/,
+                             std::vector<GlobalVector*> const& /*x*/,
+                             int const /*process_id*/)
     {
         // A hook added for solution dependent dirichlet
     }
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
index b1c57e3b14a85f7b60acae05b03f471680387b37..f9844d7558c653373e623c521b0aebd9cad5f41d 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
@@ -53,11 +53,12 @@ public:
         _boundary_conditions.push_back(std::move(bc));
     }
 
-    void preTimestep(const double t, GlobalVector const& x)
+    void preTimestep(const double t, std::vector<GlobalVector*> const& x,
+                     int const process_id)
     {
         for (auto const& bc_ptr : _boundary_conditions)
         {
-            bc_ptr->preTimestep(t, x);
+            bc_ptr->preTimestep(t, x, process_id);
         }
     }
 
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
index f3dc9a7e9cf5828d8ccdb287d6e49abc8f4fb099..587170d86728deb18c199d781fed61bb6650fe23 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
@@ -106,17 +106,18 @@ ConstraintDirichletBoundaryCondition::ConstraintDirichletBoundaryCondition(
         _integration_order, bulk_mesh, _bulk_ids);
 }
 
-void ConstraintDirichletBoundaryCondition::preTimestep(double t,
-                                                       GlobalVector const& x)
+void ConstraintDirichletBoundaryCondition::preTimestep(
+    double t, std::vector<GlobalVector*> const& x, int const process_id)
 {
     DBUG(
         "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
         "constraints");
+    auto const& solution = *x[process_id];
     for (auto const* boundary_element : _bc_mesh.getElements())
     {
         _flux_values[boundary_element->getID()] =
             _local_assemblers[boundary_element->getID()]->integrate(
-                x, t,
+                solution, t,
                 [this](std::size_t const element_id,
                        MathLib::Point3d const& pnt, double const t,
                        GlobalVector const& x) {
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
index ee9af4b741d6777c128efca9a1a2cc70447c5e08..186264607d394acd34e8e4be91e6aa84d789a89c 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
@@ -61,7 +61,8 @@ public:
                                       GlobalVector const&)>
             getFlux);
 
-    void preTimestep(double t, GlobalVector const& x) override;
+    void preTimestep(double t, std::vector<GlobalVector*> const& x,
+                     int const process_id) override;
 
     void getEssentialBCValues(
         const double t, const GlobalVector& x,
diff --git a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
index d217b698e55f8cc265853bd7b9a53208415f3e1d..2cd43338fb607bf80cf0a4c788eaa704aa742098 100644
--- a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
@@ -38,7 +38,8 @@ void PhaseFieldIrreversibleDamageOracleBoundaryCondition::getEssentialBCValues(
 
 // update new values and corresponding indices.
 void PhaseFieldIrreversibleDamageOracleBoundaryCondition::preTimestep(
-    const double /*t*/, const GlobalVector& x)
+    const double /*t*/, std::vector<GlobalVector*> const& x,
+    int const process_id)
 {
     // phase-field variable is considered irreversible if it loses more than 95%
     // of the stiffness, which is a widely used threshold.
@@ -56,7 +57,7 @@ void PhaseFieldIrreversibleDamageOracleBoundaryCondition::preTimestep(
         const auto g_idx =
             _dof_table.getGlobalIndex(l, _variable_id, _component_id);
 
-        if (x[node_id] <= irreversibleDamage)
+        if ((*x[process_id])[node_id] <= irreversibleDamage)
         {
             _bc_values.ids.emplace_back(g_idx);
             _bc_values.values.emplace_back(0.0);
diff --git a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
index b436480096d3b527fa3e9683cf919fe0ab1684bd..cd661a1db364045535d2b1a112fcb2b7724deacf 100644
--- a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
@@ -47,7 +47,8 @@ public:
         const double t, const GlobalVector& x,
         NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
 
-    void preTimestep(const double t, const GlobalVector& x) override;
+    void preTimestep(const double t, std::vector<GlobalVector*> const& x,
+                     int const process_id) override;
 
 private:
     NumLib::LocalToGlobalIndexMap const& _dof_table;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 6ed8579335b7c9b7960fb8eac0f8113897f28be5..85ab201ca55e2f393fdf6b239789d01d4cfbccda 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -186,8 +186,8 @@ void ComponentTransportProcess::
 }
 
 void ComponentTransportProcess::preTimestepConcreteProcess(
-    GlobalVector const& x, const double /*t*/, const double /*delta_t*/,
-    int const process_id)
+    std::vector<GlobalVector*> const& x, const double /*t*/,
+    const double /*delta_t*/, int const process_id)
 {
     if (_use_monolithic_scheme)
     {
@@ -197,17 +197,18 @@ void ComponentTransportProcess::preTimestepConcreteProcess(
     if (!_xs_previous_timestep[process_id])
     {
         _xs_previous_timestep[process_id] =
-            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+                *x[process_id]);
     }
     else
     {
         auto& x0 = *_xs_previous_timestep[process_id];
-        MathLib::LinAlg::copy(x, x0);
+        MathLib::LinAlg::copy(*x[process_id], x0);
     }
 }
 
 void ComponentTransportProcess::postTimestepConcreteProcess(
-    GlobalVector const& x,
+    std::vector<GlobalVector*> const& x,
     const double t,
     const double /*delta_t*/,
     int const process_id)
@@ -233,8 +234,9 @@ void ComponentTransportProcess::postTimestepConcreteProcess(
 
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
-    _surfaceflux->integrate(x, t, *this, process_id, _integration_order,
-                            _mesh, pv.getActiveElementIDs());
+    _surfaceflux->integrate(*x[process_id], t, *this, process_id,
+                            _integration_order, _mesh,
+                            pv.getActiveElementIDs());
     _surfaceflux->save(t);
 }
 
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 675d4d265984c7d842c34eeb2ead188aded0b05c..745b59a4cd06e93398b7cea915411692f423eb9c 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -126,11 +126,12 @@ public:
     void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
         int const process_id) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, const double /*t*/,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    const double /*t*/,
                                     const double /*delta_t*/,
                                     int const process_id) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
                                      const double t,
                                      const double delta_t,
                                      int const process_id) override;
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index b77a0b6b2f3b4537cbea8b03626ba4f46e2be26c..a5dc28189c59f2c890c005c5524aa84a82c03c4b 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -61,7 +61,7 @@ public:
         return _local_assemblers[element_id]->getFlux(p, t, local_x);
     }
 
-    void postTimestepConcreteProcess(GlobalVector const& x,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
                                      const double t,
                                      const double /*delta_t*/,
                                      int const process_id) override
@@ -79,8 +79,9 @@ public:
 
         ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
-        _surfaceflux->integrate(x, t, *this, process_id, _integration_order,
-                                _mesh, pv.getActiveElementIDs());
+        _surfaceflux->integrate(*x[process_id], t, *this, process_id,
+                                _integration_order, _mesh,
+                                pv.getActiveElementIDs());
         _surfaceflux->save(t);
     }
 
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index ef9fb38cd41233a9798dab6fd12a9204092dae1a..084d516f73ffb38df9dc42e1e411977c22d77d54 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -154,7 +154,7 @@ void HTProcess::assembleWithJacobianConcreteProcess(
         dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
-void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
+void HTProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
                                            const double /*t*/,
                                            const double /*delta_t*/,
                                            const int process_id)
@@ -169,12 +169,13 @@ void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
     if (!_xs_previous_timestep[process_id])
     {
         _xs_previous_timestep[process_id] =
-            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+                *x[process_id]);
     }
     else
     {
         auto& x0 = *_xs_previous_timestep[process_id];
-        MathLib::LinAlg::copy(x, x0);
+        MathLib::LinAlg::copy(*x[process_id], x0);
     }
 
     auto& x0 = *_xs_previous_timestep[process_id];
@@ -254,7 +255,7 @@ Eigen::Vector3d HTProcess::getFlux(std::size_t element_id,
 }
 
 // this is almost a copy of the implementation in the GroundwaterFlow
-void HTProcess::postTimestepConcreteProcess(GlobalVector const& x,
+void HTProcess::postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
                                             const double t,
                                             const double /*delta_t*/,
                                             int const process_id)
@@ -278,7 +279,8 @@ void HTProcess::postTimestepConcreteProcess(GlobalVector const& x,
 
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
 
-    _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
+    _surfaceflux->integrate(*x[process_id], t, *this, process_id,
+                            _integration_order, _mesh,
                             pv.getActiveElementIDs());
     _surfaceflux->save(t);
 }
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 67c0905fe9f6ed10f7b5301d45ceb32a673e7fbf..00e87483be6a303e1051bc8a8ae159dfbb3d51c9 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -81,7 +81,7 @@ public:
     void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
         int const process_id) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
                                      const double t,
                                      const double delta_t,
                                      int const process_id) override;
@@ -103,8 +103,8 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
                                     const int process_id) override;
 
     void setCoupledSolutionsOfPreviousTimeStepPerProcess(const int process_id);
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
index 50580ee1dffa1aec5edae7df9764b7528928f648..21ac46ffd983eb7906ba6ef5e8053eaeaaa85f4d 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
@@ -455,7 +455,7 @@ void HydroMechanicsProcess<DisplacementDim>::
 
 template <int DisplacementDim>
 void HydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PreTimestep HydroMechanicsProcess.");
@@ -466,21 +466,22 @@ void HydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
             getProcessVariables(process_id)[0];
         GlobalExecutor::executeSelectedMemberOnDereferenced(
             &LocalAssemblerInterface::preTimestep, _local_assemblers,
-            pv.getActiveElementIDs(), *_local_to_global_index_map, x, t,
-            dt);
+            pv.getActiveElementIDs(), *_local_to_global_index_map,
+            *x[process_id], t, dt);
     }
 }
 
 template <int DisplacementDim>
 void HydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PostTimestep HydroMechanicsProcess.");
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::postTimestep, _local_assemblers,
-        pv.getActiveElementIDs(), getDOFTable(process_id), x, t, dt);
+        pv.getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
+        dt);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index f47eebb68741c33caecd14a583f70d513ffc2720..7db3cb33b77d6ca02aaa21e4cc409e19fb792357 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -83,12 +83,12 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
                                     const int process_id) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                     const double delta_t,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     const double t, const double delta_t,
                                      int const process_id) override;
 
     void postNonLinearSolverConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 5185fd833994fcbfa764145d7d34e786e308861f..733ba79e49481a3926b7d23bf31c3fc1e9621ddc 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -455,7 +455,7 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
 
 template <int GlobalDim>
 void HydroMechanicsProcess<GlobalDim>::postTimestepConcreteProcess(
-    GlobalVector const& x, const double t, double const dt,
+    std::vector<GlobalVector*> const& x, const double t, double const dt,
     int const process_id)
 {
     DBUG("Compute the secondary variables for HydroMechanicsProcess.");
@@ -467,7 +467,8 @@ void HydroMechanicsProcess<GlobalDim>::postTimestepConcreteProcess(
 
         GlobalExecutor::executeSelectedMemberOnDereferenced(
             &HydroMechanicsLocalAssemblerInterface::postTimestep,
-            _local_assemblers, pv.getActiveElementIDs(), dof_table, x, t, dt);
+            _local_assemblers, pv.getActiveElementIDs(), dof_table,
+            *x[process_id], t, dt);
     }
 
     // Copy displacement jumps in a solution vector to mesh property
@@ -491,7 +492,7 @@ void HydroMechanicsProcess<GlobalDim>::postTimestepConcreteProcess(
         g_variable_id = static_cast<int>(std::distance(pvs.begin(), it));
     }
 
-    MathLib::LinAlg::setLocalAccessibleVector(x);
+    MathLib::LinAlg::setLocalAccessibleVector(*x[process_id]);
 
     const int monolithic_process_id = 0;
     ProcessVariable& pv_g =
@@ -512,7 +513,7 @@ void HydroMechanicsProcess<GlobalDim>::postTimestepConcreteProcess(
             auto const global_index =
                 dof_table.getGlobalIndex(l, g_variable_id, component_id);
             mesh_prop_g[node->getID() * num_comp + component_id] =
-                x[global_index];
+                (*x[process_id])[global_index];
         }
     }
 
@@ -610,8 +611,8 @@ void HydroMechanicsProcess<GlobalDim>::assembleWithJacobianConcreteProcess(
 
 template <int GlobalDim>
 void HydroMechanicsProcess<GlobalDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
-    const int process_id)
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
+    int const process_id)
 {
     DBUG("PreTimestep HydroMechanicsProcess.");
 
@@ -619,8 +620,8 @@ void HydroMechanicsProcess<GlobalDim>::preTimestepConcreteProcess(
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &HydroMechanicsLocalAssemblerInterface::preTimestep, _local_assemblers,
-        pv.getActiveElementIDs(), *_local_to_global_index_map,
-        x, t, dt);
+        pv.getActiveElementIDs(), *_local_to_global_index_map, *x[process_id],
+        t, dt);
 }
 
 // ------------------------------------------------------------------------------------
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index 068b854886293ce03f778019868a943fc7015212..0497533dc00bf8761a627348774799024b2a420e 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -51,8 +51,8 @@ public:
     bool isLinear() const override;
     //! @}
 
-    void postTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                     double const dt,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     double const t, double const dt,
                                      int const process_id) override;
 
 private:
@@ -75,9 +75,9 @@ private:
         GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
-                                    const int /*process_id*/) override;
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
+                                    int const process_id) override;
 
 private:
     HydroMechanicsProcessData<GlobalDim> _process_data;
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index 18a9aff117a17e698d06c9841e63d0fa07fe867f..96e9897755eb14c8bbcab4ecaf35f4107b87cd70 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -579,7 +579,7 @@ void SmallDeformationProcess<DisplacementDim>::
 }
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PreTimestep SmallDeformationProcess.");
@@ -588,7 +588,7 @@ void SmallDeformationProcess<DisplacementDim>::preTimestepConcreteProcess(
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &SmallDeformationLocalAssemblerInterface::preTimestep,
         _local_assemblers, pv.getActiveElementIDs(),
-        *_local_to_global_index_map, x, t, dt);
+        *_local_to_global_index_map, *x[process_id], t, dt);
 }
 // ------------------------------------------------------------------------------------
 // template instantiation
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
index 41355af7d0b728a0f3b31de50c4991a3376bc618..9fc97d66e6b4184cfdc7566629cf00aa518ea464 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.h
@@ -72,9 +72,9 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
-                                    const int /*process_id*/) override;
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
+                                    const int process_id) override;
 
 private:
     SmallDeformationProcessData<DisplacementDim> _process_data;
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.cpp b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
index 37d31666d9582fe96d27602b711773b5beb76a99..0a3df108894419314acd23bf9cbf14b165ecabb2 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
@@ -239,7 +239,7 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
 
 template <int DisplacementDim>
 void PhaseFieldProcess<DisplacementDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PreTimestep PhaseFieldProcess %d.", process_id);
@@ -250,13 +250,14 @@ void PhaseFieldProcess<DisplacementDim>::preTimestepConcreteProcess(
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::preTimestep, _local_assemblers,
-        pv.getActiveElementIDs(), getDOFTable(process_id), x, t, dt);
+        pv.getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
+        dt);
 }
 
 template <int DisplacementDim>
 void PhaseFieldProcess<DisplacementDim>::postTimestepConcreteProcess(
-    GlobalVector const& x, const double t, const double /*delta_t*/,
-    int const process_id)
+    std::vector<GlobalVector*> const& x, const double t,
+    const double /*delta_t*/, int const process_id)
 {
     if (isPhaseFieldProcess(process_id))
     {
@@ -277,7 +278,7 @@ void PhaseFieldProcess<DisplacementDim>::postTimestepConcreteProcess(
 
         GlobalExecutor::executeSelectedMemberOnDereferenced(
             &LocalAssemblerInterface::computeEnergy, _local_assemblers,
-            pv.getActiveElementIDs(), dof_tables, x, t,
+            pv.getActiveElementIDs(), dof_tables, *x[process_id], t,
             _process_data.elastic_energy, _process_data.surface_energy,
             _process_data.pressure_work, _coupled_solutions);
 
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.h b/ProcessLib/PhaseField/PhaseFieldProcess.h
index 066502e3c6ce556a2f6c92d416aa9b01fbab561c..8868d7d34989e30d565c79c381e7766182b15b26 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.h
@@ -98,12 +98,12 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
                                     const int process_id) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                     const double delta_t,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     const double t, const double delta_t,
                                      int const process_id) override;
 
     void postNonLinearSolverConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 7e8815d8209a7acd1951723d1d11e19387dcbdd7..17fd5c5feb19c903d0f0f467e0ae82d9096a26fc 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -378,7 +378,7 @@ void Process::computeSparsityPattern()
         NumLib::computeSparsityPattern(*_local_to_global_index_map, _mesh);
 }
 
-void Process::preTimestep(GlobalVector const& x, const double t,
+void Process::preTimestep(std::vector<GlobalVector*> const& x, const double t,
                           const double delta_t, const int process_id)
 {
     for (auto& cached_var : _cached_secondary_variables)
@@ -386,16 +386,18 @@ void Process::preTimestep(GlobalVector const& x, const double t,
         cached_var->setTime(t);
     }
 
-    MathLib::LinAlg::setLocalAccessibleVector(x);
+    for (auto* const solution : x)
+        MathLib::LinAlg::setLocalAccessibleVector(*solution);
     preTimestepConcreteProcess(x, t, delta_t, process_id);
 
-    _boundary_conditions[process_id].preTimestep(t, x);
+    _boundary_conditions[process_id].preTimestep(t, x, process_id);
 }
 
-void Process::postTimestep(GlobalVector const& x, const double t,
+void Process::postTimestep(std::vector<GlobalVector*> const& x, const double t,
                            const double delta_t, int const process_id)
 {
-    MathLib::LinAlg::setLocalAccessibleVector(x);
+    for (auto* const solution : x)
+        MathLib::LinAlg::setLocalAccessibleVector(*solution);
     postTimestepConcreteProcess(x, t, delta_t, process_id);
 }
 
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index addd0df65ad2107b975ceb12082036a48d7a27bf..2a08e5b3b4dbf2066872f512dc0fac61249bad04 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -59,11 +59,11 @@ public:
             const bool use_monolithic_scheme = true);
 
     /// Preprocessing before starting assembly for new timestep.
-    void preTimestep(GlobalVector const& x, const double t,
+    void preTimestep(std::vector<GlobalVector*> const& x, const double t,
                      const double delta_t, const int process_id);
 
     /// Postprocessing after a complete timestep.
-    void postTimestep(GlobalVector const& x, const double t,
+    void postTimestep(std::vector<GlobalVector*> const& x, const double t,
                       const double delta_t, int const process_id);
 
     /// Calculates secondary variables, e.g. stress and strain for deformation
@@ -206,17 +206,19 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) = 0;
 
-    virtual void preTimestepConcreteProcess(GlobalVector const& /*x*/,
-                                            const double /*t*/,
-                                            const double /*delta_t*/,
-                                            const int /*process_id*/)
+    virtual void preTimestepConcreteProcess(
+        std::vector<GlobalVector*> const& /*x*/,
+        const double /*t*/,
+        const double /*delta_t*/,
+        const int /*process_id*/)
     {
     }
 
-    virtual void postTimestepConcreteProcess(GlobalVector const& /*x*/,
-                                             const double /*t*/,
-                                             const double /*delta_t*/,
-                                             int const /*process_id*/)
+    virtual void postTimestepConcreteProcess(
+        std::vector<GlobalVector*> const& /*x*/,
+        const double /*t*/,
+        const double /*delta_t*/,
+        int const /*process_id*/)
     {
     }
 
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index acadc4d3673e68955cd2fe13ed49528721f2bf60..1a15b2618a8a758f492d8fbfbdf9f79e072cd72a 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -324,7 +324,7 @@ void RichardsMechanicsProcess<DisplacementDim>::
 
 template <int DisplacementDim>
 void RichardsMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PreTimestep RichardsMechanicsProcess.");
@@ -335,8 +335,8 @@ void RichardsMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
             getProcessVariables(process_id)[0];
         GlobalExecutor::executeSelectedMemberOnDereferenced(
             &LocalAssemblerInterface::preTimestep, _local_assemblers,
-            pv.getActiveElementIDs(), *_local_to_global_index_map, x,
-            t, dt);
+            pv.getActiveElementIDs(), *_local_to_global_index_map,
+            *x[process_id], t, dt);
     }
 }
 
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
index eaa4e943008508e758e55c978265ecd9f860f668..6943ac5170f923e59ff22fa34cf41ec2e34bb21f 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
@@ -85,8 +85,8 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
                                     const int process_id) override;
 
     void postNonLinearSolverConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
index bda7c7f63dbc0d8e1b2877d5b6f853c9b8351ba1..555d0a0889b1460abf118d2671d3ac4d472d52c6 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
@@ -283,7 +283,7 @@ void SmallDeformationProcess<DisplacementDim>::
 
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::postTimestepConcreteProcess(
-    GlobalVector const& x, const double t, const double dt,
+    std::vector<GlobalVector*> const& x, const double t, const double dt,
     int const process_id)
 {
     DBUG("PostTimestep SmallDeformationProcess.");
@@ -292,11 +292,13 @@ void SmallDeformationProcess<DisplacementDim>::postTimestepConcreteProcess(
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::postTimestep, _local_assemblers,
-        pv.getActiveElementIDs(), *_local_to_global_index_map, x, t, dt);
+        pv.getActiveElementIDs(), *_local_to_global_index_map, *x[process_id],
+        t, dt);
 
     std::unique_ptr<GlobalVector> material_forces;
     ProcessLib::SmallDeformation::writeMaterialForces(
-        material_forces, _local_assemblers, *_local_to_global_index_map, x);
+        material_forces, _local_assemblers, *_local_to_global_index_map,
+        *x[process_id]);
 
     material_forces->copyValues(*_material_forces);
 }
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index 7846afd5cccaf9fc1c7db3ba8c13a47532d3ba04..ee555c426bc86c354e52161fe9a4e6bf288dade1 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -61,8 +61,8 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                     const double dt,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     const double t, const double dt,
                                      int const process_id) override;
 
 private:
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
index 8b01241720a3cc51a7891ab6f489079736a16357..7a46b97e4059af641372fec4e1038bc3cfcef8f2 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
@@ -293,11 +293,11 @@ void SmallDeformationNonlocalProcess<DisplacementDim>::
 }
 
 template <int DisplacementDim>
-void SmallDeformationNonlocalProcess<
-    DisplacementDim>::postTimestepConcreteProcess(GlobalVector const& x,
-                                                  double const t,
-                                                  double const dt,
-                                                  int const process_id)
+void SmallDeformationNonlocalProcess<DisplacementDim>::
+    postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                double const t,
+                                double const dt,
+                                int const process_id)
 {
     DBUG("PostTimestep SmallDeformationNonlocalProcess.");
 
@@ -305,7 +305,8 @@ void SmallDeformationNonlocalProcess<
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::postTimestep, _local_assemblers,
-        pv.getActiveElementIDs(), *_local_to_global_index_map, x, t, dt);
+        pv.getActiveElementIDs(), *_local_to_global_index_map, *x[process_id],
+        t, dt);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
index 100f3bac6e7fb54a5d1d985854d73a91d079b198..4450b5c1ae0dc0cbcb1f3380fc4743cfa92444c0 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.h
@@ -67,8 +67,8 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                     double const dt,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     double const t, double const dt,
                                      int const process_id) override;
 
     NumLib::IterationResult postIterationConcreteProcess(
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 47a9c9841d6b99861cccbe604e50a5b0e8401121..39faa68e638606bd66711f7caa3d15d47ff96b05 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -263,10 +263,10 @@ void TESProcess::assembleWithJacobianConcreteProcess(
         dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 
-void TESProcess::preTimestepConcreteProcess(GlobalVector const& x,
+void TESProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
                                             const double t,
                                             const double delta_t,
-                                            const int /*process_id*/)
+                                            const int process_id)
 {
     DBUG("new timestep");
 
@@ -275,7 +275,7 @@ void TESProcess::preTimestepConcreteProcess(GlobalVector const& x,
     ++_assembly_params.timestep;  // TODO remove that
 
     _x_previous_timestep =
-        MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+        MathLib::MatrixVectorTraits<GlobalVector>::newInstance(*x[process_id]);
 }
 
 void TESProcess::preIterationConcreteProcess(const unsigned iter,
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index 765ade2fe1da8c6ba7f189a56e68af32f5e57efa..3f5ef525adcedc8b9a377ae2a29fc7a1d5b8b91d 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -44,8 +44,8 @@ public:
         NumLib::NamedFunctionCaller&& named_function_caller,
         BaseLib::ConfigTree const& config);
 
-    void preTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                    const double delta_t,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    const double t, const double delta_t,
                                     const int process_id) override;
     void preIterationConcreteProcess(const unsigned iter,
                                      GlobalVector const& x) override;
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
index b5ccf3584c46470a16bf0931920352176eb9bbdb..e3541447be01c1ff74aa1784fcb51d782436cc89 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
@@ -110,7 +110,7 @@ void ThermalTwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
         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,
+    std::vector<GlobalVector*> const& x, double const t, double const delta_t,
     const int process_id)
 {
     DBUG("PreTimestep ThermalTwoPhaseFlowWithPPProcess.");
@@ -118,8 +118,8 @@ void ThermalTwoPhaseFlowWithPPProcess::preTimestepConcreteProcess(
     ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::preTimestep, _local_assemblers,
-        pv.getActiveElementIDs(), *_local_to_global_index_map, x, t,
-        delta_t);
+        pv.getActiveElementIDs(), *_local_to_global_index_map, *x[process_id],
+        t, delta_t);
 }
 
 }  // namespace ThermalTwoPhaseFlowWithPP
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h
index 0087c6921a45f8b9a575edca51c17ead7f91b90c..71d7bdb7c67e8818f12d02229d15e81dbf0d4151 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.h
@@ -72,8 +72,8 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                    const double delta_t,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    const double t, const double delta_t,
                                     const int process_id) override;
 
     ThermalTwoPhaseFlowWithPPProcessData _process_data;
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
index 2f577aacdf593e95624086113ebafdcb9de56985..705f7280a6f82c28d4b6ad2b59b4987435b1988f 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
@@ -332,7 +332,7 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::
 
 template <int DisplacementDim>
 void ThermoHydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PreTimestep ThermoHydroMechanicsProcess.");
@@ -341,19 +341,19 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
     {
         GlobalExecutor::executeMemberOnDereferenced(
             &LocalAssemblerInterface::preTimestep, _local_assemblers,
-            *_local_to_global_index_map, x, t, dt);
+            *_local_to_global_index_map, *x[process_id], t, dt);
     }
 }
 
 template <int DisplacementDim>
 void ThermoHydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PostTimestep ThermoHydroMechanicsProcess.");
     GlobalExecutor::executeMemberOnDereferenced(
         &LocalAssemblerInterface::postTimestep, _local_assemblers,
-        getDOFTable(process_id), x, t, dt);
+        getDOFTable(process_id), *x[process_id], t, dt);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
index da51d994b9b10a80f62cd4c19afb9444310b2a47..f96bff67c79876b5c4ec621071f75f11dac05391 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.h
@@ -86,12 +86,12 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
                                     const int process_id) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                     const double delta_t,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     const double t, const double delta_t,
                                      int const process_id) override;
 
     void postNonLinearSolverConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
index 2e3ad20164dbe7c9fd0a3b2455e2b6f3849db269..09f6e54861a84a642bb14b9d05ed405ce38a4d26 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
@@ -260,11 +260,11 @@ void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
 }
 
 template <int DisplacementDim>
-void ThermoMechanicalPhaseFieldProcess<
-    DisplacementDim>::preTimestepConcreteProcess(GlobalVector const& x,
-                                                 double const t,
-                                                 double const dt,
-                                                 const int process_id)
+void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
+    preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                               double const t,
+                               double const dt,
+                               const int process_id)
 {
     DBUG("PreTimestep ThermoMechanicalPhaseFieldProcess.");
 
@@ -277,16 +277,16 @@ void ThermoMechanicalPhaseFieldProcess<
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &ThermoMechanicalPhaseFieldLocalAssemblerInterface::preTimestep,
-        _local_assemblers, pv.getActiveElementIDs(), getDOFTable(process_id), x,
-        t, dt);
+        _local_assemblers, pv.getActiveElementIDs(), getDOFTable(process_id),
+        *x[process_id], t, dt);
 }
 
 template <int DisplacementDim>
-void ThermoMechanicalPhaseFieldProcess<
-    DisplacementDim>::postTimestepConcreteProcess(GlobalVector const& x,
-                                                  double const t,
-                                                  double const dt,
-                                                  int const process_id)
+void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
+    postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                double const t,
+                                double const dt,
+                                int const process_id)
 {
     DBUG("PostTimestep ThermoMechanicalPhaseFieldProcess.");
 
@@ -294,8 +294,8 @@ void ThermoMechanicalPhaseFieldProcess<
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &ThermoMechanicalPhaseFieldLocalAssemblerInterface::postTimestep,
-        _local_assemblers, pv.getActiveElementIDs(), getDOFTable(process_id), x,
-        t, dt);
+        _local_assemblers, pv.getActiveElementIDs(), getDOFTable(process_id),
+        *x[process_id], t, dt);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h
index c18fd2ae246164a589e06757b774ca3b65763d30..4ef734a38a8776da1d8c45aad47dc6e04d3a68a3 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.h
@@ -103,12 +103,12 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
                                     const int process_id) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                     double const dt,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     double const t, double const dt,
                                      int const process_id) override;
 
     void postNonLinearSolverConcreteProcess(GlobalVector const& x,
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
index 6c149fd75f89b182451f7f1775ffd59c91060874..6b813c5b6f2d0571fb7a143339c7a7266d07d13b 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -383,7 +383,7 @@ void ThermoMechanicsProcess<DisplacementDim>::
 
 template <int DisplacementDim>
 void ThermoMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PreTimestep ThermoMechanicsProcess.");
@@ -397,19 +397,20 @@ void ThermoMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
         GlobalExecutor::executeSelectedMemberOnDereferenced(
             &ThermoMechanicsLocalAssemblerInterface::preTimestep,
             _local_assemblers, pv.getActiveElementIDs(),
-            *_local_to_global_index_map, x, t, dt);
+            *_local_to_global_index_map, *x[process_id], t, dt);
         return;
     }
 
     // For the staggered scheme.
     if (!_previous_T)
     {
-        _previous_T = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+        _previous_T = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+            *x[process_id]);
     }
     else
     {
         auto& x0 = *_previous_T;
-        MathLib::LinAlg::copy(x, x0);
+        MathLib::LinAlg::copy(*x[process_id], x0);
     }
 
     auto& x0 = *_previous_T;
@@ -418,7 +419,7 @@ void ThermoMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
 
 template <int DisplacementDim>
 void ThermoMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     int const process_id)
 {
     if (process_id != _process_data.mechanics_process_id)
@@ -433,7 +434,7 @@ void ThermoMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &ThermoMechanicsLocalAssemblerInterface::postTimestep,
         _local_assemblers, pv.getActiveElementIDs(),
-        *_local_to_global_index_map, x, t, dt);
+        *_local_to_global_index_map, *x[process_id], t, dt);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
index 24e76adbdaa646f2fd4f19827bb57e5b111dde7f..1e4026d04cf7fa6f1f068a0e13475caac06defc8 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -78,12 +78,12 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
-                                    double const dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    double const t, double const dt,
                                     const int process_id) override;
 
-    void postTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                     const double delta_t,
+    void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                     const double t, const double delta_t,
                                      int const process_id) override;
 
     NumLib::LocalToGlobalIndexMap const& getDOFTable(
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 2864cfddadc0cc82f31b314c31066249ded76124..96e35a35b795db881f8bef16e7d8fb78aca353e9 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -528,9 +528,8 @@ void preTimestepForAllProcesses(
     for (auto& process_data : per_process_data)
     {
         auto const process_id = process_data->process_id;
-        auto& x = *_process_solutions[process_id];
         auto& pcs = process_data->process;
-        pcs.preTimestep(x, t, dt, process_id);
+        pcs.preTimestep(_process_solutions, t, dt, process_id);
     }
 }
 
@@ -575,7 +574,7 @@ void postTimestepForAllProcesses(
             pcs.setCoupledSolutionsForStaggeredScheme(&coupled_solutions);
         }
         auto& x = *process_solutions[process_id];
-        pcs.postTimestep(x, t, dt, process_id);
+        pcs.postTimestep(process_solutions, t, dt, process_id);
         pcs.computeSecondaryVariable(t, x, process_id);
     }
 }
@@ -772,7 +771,7 @@ void TimeLoop::outputSolutions(bool const output_initial_condition,
 
         if (output_initial_condition)
         {
-            pcs.preTimestep(x, _start_time,
+            pcs.preTimestep(_process_solutions, _start_time,
                             process_data->timestepper->getTimeStep().dt(),
                             process_id);
             // Update secondary variables, which might be uninitialized, before
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
index 2c1ccb4ecfb51b68f2558ae18a7d389c9069a675..c0540947688be9099ef19171ce44b01a143cd3fb 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -108,7 +108,7 @@ void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
         dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
 }
 void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
-    GlobalVector const& x, double const t, double const dt,
+    std::vector<GlobalVector*> const& x, double const t, double const dt,
     const int process_id)
 {
     DBUG("PreTimestep TwoPhaseFlowWithPrhoProcess.");
@@ -117,8 +117,8 @@ void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
 
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerInterface::preTimestep, _local_assemblers,
-        pv.getActiveElementIDs(), *_local_to_global_index_map, x, t,
-        dt);
+        pv.getActiveElementIDs(), *_local_to_global_index_map, *x[process_id],
+        t, dt);
 }
 
 }  // namespace TwoPhaseFlowWithPrho
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
index bab3f17d6f1dd9a9ec4e78e87c86a284be3b5398..4c67a9b1676967a7ffd908fd34126fa3d3004e1e 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
@@ -68,8 +68,8 @@ private:
         int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
         GlobalMatrix& Jac) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                    const double dt,
+    void preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+                                    const double t, const double dt,
                                     const int process_id) override;
 
     TwoPhaseFlowWithPrhoProcessData _process_data;