diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
index 1c5c1fbc2f2239eaf01f4a3e5bd5aa31748f07dc..a594ca6f3906b53e62af33c479fc0ad1504d36e4 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
@@ -34,6 +34,9 @@ std::vector<std::vector<double>> getPreviousLocalSolutions(
         std::reference_wrapper<const std::vector<GlobalIndexType>>>&
         indices)
 {
+    if (cpl_xs.coupled_xs_t0.empty())
+        return {};
+
     const auto number_of_coupled_solutions = cpl_xs.coupled_xs.size();
     std::vector<std::vector<double>> local_xs_t0;
     local_xs_t0.reserve(number_of_coupled_solutions);
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 598968e1c2d6aca9ea730d0bc928c570628e3f34..5567b0ffd4af4d3b7e72c2e0bbd70ac10831d262 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -197,17 +197,13 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         auto C = _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u);
 
-        local_Jac
-            .template block<displacement_size, displacement_size>(
-                displacement_index, displacement_index)
-            .noalias() += B.transpose() * C * B * w;
+        local_Jac.noalias() += B.transpose() * C * B * w;
 
         double p_at_xi = 0.;
         NumLib::shapeFunctionInterpolate(local_p, N_p, p_at_xi);
 
         double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
-        local_rhs.template segment<displacement_size>(displacement_index)
-            .noalias() -=
+        local_rhs.noalias() -=
             (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
              N_u_op.transpose() * rho * b) *
             w;
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index 56c17d0ebd1f2c030161f2516f24b533e2cfb4b1..d196a80401431dd224298b7607ab662bc92bed53 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -658,6 +658,7 @@ bool UncoupledProcessesTimeLoop::loop()
             auto& pcs = spd->process;
             _output->addProcess(pcs, pcs_idx);
 
+            spd->equation_id = pcs_idx;
             setTimeDiscretizedODESystem(*spd);
 
             if (auto* conv_crit =
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 9853293febf9643ba95572590c3205c5647463fa..e167bdee25918b8f1faa53d8c85eb082698873f6 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -115,7 +115,8 @@ void VectorMatrixAssembler::assembleWithJacobian(
     // coupling.
     auto const indices =
         ((base_dof_table == nullptr) ||
-         (cpl_xs->process_id == cpl_xs->coupled_xs.size() - 1 &&
+         (cpl_xs->process_id ==
+              static_cast<int>(cpl_xs->coupled_xs.size()) - 1 &&
           cpl_xs != nullptr))
             ? NumLib::getIndices(mesh_item_id, dof_table)
             : NumLib::getIndices(mesh_item_id, *base_dof_table);
@@ -135,58 +136,32 @@ void VectorMatrixAssembler::assembleWithJacobian(
     }
     else
     {
-        std::vector<std::reference_wrapper<const std::vector<GlobalIndexType>>>
-            indices_of_all_coupled_processes;
-        indices_of_all_coupled_processes.reserve(cpl_xs->coupled_xs.size());
         if (base_dof_table == nullptr)
         {
-            for (std::size_t i = 0; i < cpl_xs->coupled_xs.size(); i++)
-            {
-                indices_of_all_coupled_processes.emplace_back(
-                    std::ref(indices));
-            }
+            localAssembleWithJacobianAndCoupling(t, indices, indices,
+                                                 local_xdot, local_assembler,
+                                                 dxdot_dx, dx_dx, cpl_xs);
         }
         else
         {
-            if (cpl_xs->process_id == cpl_xs->coupled_xs.size() - 1)
+            if (cpl_xs->process_id ==
+                static_cast<int>(cpl_xs->coupled_xs.size()) - 1)
             {
                 const auto base_indices =
                     NumLib::getIndices(mesh_item_id, *base_dof_table);
-                for (std::size_t i = 0; i < cpl_xs->coupled_xs.size() - 1; i++)
-                {
-                    indices_of_all_coupled_processes.emplace_back(
-                        std::ref(base_indices));
-                }
-                indices_of_all_coupled_processes.emplace_back(
-                    std::ref(indices));
+                localAssembleWithJacobianAndCoupling(
+                    t, base_indices, indices, local_xdot, local_assembler,
+                    dxdot_dx, dx_dx, cpl_xs);
             }
             else
             {
                 const auto full_indices =
                     NumLib::getIndices(mesh_item_id, dof_table);
-                for (std::size_t i = 0; i < cpl_xs->coupled_xs.size() - 1; i++)
-                {
-                    indices_of_all_coupled_processes.emplace_back(
-                        std::ref(indices));
-                }
-                indices_of_all_coupled_processes.emplace_back(
-                    std::ref(full_indices));
+                localAssembleWithJacobianAndCoupling(
+                    t, indices, full_indices, local_xdot, local_assembler,
+                    dxdot_dx, dx_dx, cpl_xs);
             }
         }
-
-        auto local_coupled_xs0 = getPreviousLocalSolutions(
-            *cpl_xs, indices_of_all_coupled_processes);
-        auto local_coupled_xs =
-            getCurrentLocalSolutions(*cpl_xs, indices_of_all_coupled_processes);
-
-        ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-            cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
-            std::move(local_coupled_xs));
-
-        _jacobian_assembler->assembleWithJacobianAndCoupling(
-            local_assembler, t, local_xdot, dxdot_dx, dx_dx, _local_M_data,
-            _local_K_data, _local_b_data, _local_Jac_data,
-            local_coupled_solutions);
     }
 
     auto const num_r_c = indices.size();
@@ -222,4 +197,34 @@ void VectorMatrixAssembler::assembleWithJacobian(
     }
 }
 
+void VectorMatrixAssembler::localAssembleWithJacobianAndCoupling(
+    const double t, std::vector<GlobalIndexType> const& base_indices,
+    std::vector<GlobalIndexType> const& full_indices,
+    std::vector<double> const& local_xdot,
+    LocalAssemblerInterface& local_assembler, const double dxdot_dx,
+    const double dx_dx, CoupledSolutionsForStaggeredScheme const* cpl_xs)
+{
+    std::vector<std::reference_wrapper<const std::vector<GlobalIndexType>>>
+        indices_of_all_coupled_processes;
+    indices_of_all_coupled_processes.reserve(cpl_xs->coupled_xs.size());
+    for (std::size_t i = 0; i < cpl_xs->coupled_xs.size() - 1; i++)
+    {
+        indices_of_all_coupled_processes.emplace_back(std::ref(base_indices));
+    }
+    indices_of_all_coupled_processes.emplace_back(std::ref(full_indices));
+
+    auto local_coupled_xs0 =
+        getPreviousLocalSolutions(*cpl_xs, indices_of_all_coupled_processes);
+    auto local_coupled_xs =
+        getCurrentLocalSolutions(*cpl_xs, indices_of_all_coupled_processes);
+
+    ProcessLib::LocalCoupledSolutions local_coupled_solutions(
+        cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
+        std::move(local_coupled_xs));
+
+    _jacobian_assembler->assembleWithJacobianAndCoupling(
+        local_assembler, t, local_xdot, dxdot_dx, dx_dx, _local_M_data,
+        _local_K_data, _local_b_data, _local_Jac_data, local_coupled_solutions);
+}
+
 }  // namespace ProcessLib
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index 733ebde8dcb4ff9d7e839cf944cb40402f5745f6..0f41b856681b80622bef07f2c098e4e5d811931b 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -70,6 +70,13 @@ private:
 
     //! Used to assemble the Jacobian.
     std::unique_ptr<AbstractJacobianAssembler> _jacobian_assembler;
+
+    void localAssembleWithJacobianAndCoupling(
+        const double t, std::vector<GlobalIndexType> const& base_indices,
+        std::vector<GlobalIndexType> const& full_indices,
+        std::vector<double> const& local_xdot,
+        LocalAssemblerInterface& local_assembler, const double dxdot_dx,
+        const double dx_dx, CoupledSolutionsForStaggeredScheme const* cpl_xs);
 };
 
 }  // namespace ProcessLib