diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 39373dfa025bd5d47702cebabbfcf3eddf8b4da7..f8d562512ea75bd441cd4accb8db74976643defa 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -110,7 +110,7 @@ void ComponentTransportProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void ComponentTransportProcess::setCoupledSolutionsOfPreviousTimeStep()
@@ -141,7 +141,7 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 Eigen::Vector3d ComponentTransportProcess::getFlux(
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index a669c9c68d877ac4dfe8f19a947721ceb829f9e9..4701e5de261637b1e4b1efe94bc08b33d7b21f36 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -115,7 +115,7 @@ void HTProcess::assembleConcreteProcess(const double t, double const dt,
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void HTProcess::assembleWithJacobianConcreteProcess(
@@ -144,7 +144,7 @@ void HTProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 void HTProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 1c6cbfa30daae469a1b1d497843e7c191ba2a13c..7124426f8572959b4c29f2c64f115b487242108d 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -90,7 +90,7 @@ void HeatConductionProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 
     auto const residuum = computeResiduum(*x[0], *xdot[0], M, K, b);
     transformVariableFromGlobalVector(residuum, 0 /*variable id*/,
@@ -114,7 +114,7 @@ void HeatConductionProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     transformVariableFromGlobalVector(b, 0 /*variable id*/,
                                       *_local_to_global_index_map, *_heat_flux,
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
index 7d57f19c365f231bbf1361b8ebb5404fe0496b57..a7f40e2c8bc8f1a307ef690af68234c128b69098 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
@@ -172,7 +172,7 @@ void HeatTransportBHEProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void HeatTransportBHEProcess::assembleWithJacobianConcreteProcess(
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
index 625bb2f4aff25fc6cf5b474f5a87507e23fd15ea..b994ada35a7098faa7fc7d3b937c341f39026e78 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
@@ -398,7 +398,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, xdot, process_id, M, K, b, _coupled_solutions);
+        dof_table, t, dt, x, xdot, process_id, M, K, b);
 }
 
 template <int DisplacementDim>
@@ -443,7 +443,7 @@ void HydroMechanicsProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         if (_use_monolithic_scheme)
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 49c2e6855f110dfec09a07e0f2041eacbaa69571..b28d2aeff669450cb386b5698c2b214927eb4c0e 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -579,7 +579,7 @@ 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, xdot, process_id, M, K, b, _coupled_solutions);
+        dof_table, t, dt, x, xdot, process_id, M, K, b);
 }
 
 template <int GlobalDim>
@@ -599,7 +599,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, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         transformVariableFromGlobalVector(b, variable_id,
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index e3de8775e3087483a6b3e71c5d1b14864999d6f2..26deeeb36693e9bf686b377661a3fb3496343da5 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -562,7 +562,7 @@ void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::
@@ -582,7 +582,7 @@ void SmallDeformationProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 template <int DisplacementDim>
 void SmallDeformationProcess<DisplacementDim>::preTimestepConcreteProcess(
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 25c28c5996c85bb341fd8feaf407f6a68e2a5869..ae54d245ff1d829935e1652332735943f28e08a1 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -83,7 +83,7 @@ void LiquidFlowProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 
     auto const residuum = computeResiduum(*x[0], *xdot[0], M, K, b);
     transformVariableFromGlobalVector(residuum, 0 /*variable id*/,
@@ -107,7 +107,7 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 void LiquidFlowProcess::computeSecondaryVariableConcrete(
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.cpp b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
index b8ac1967353d8e1ac6f36fb4ce3c4c367ff7263d..24914bb733826a6b6016c087a05f4560801757ab 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
@@ -191,7 +191,7 @@ void PhaseFieldProcess<DisplacementDim>::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 template <int DisplacementDim>
@@ -226,7 +226,7 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     if (process_id == 0)
     {
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
index 5bdc63315b0d97e36dcfa2318e9d37c75e6db743..e424e68e34b42eefd3461224579248ab32f90b73 100644
--- a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
@@ -78,7 +78,7 @@ void RichardsComponentTransportProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess(
@@ -97,7 +97,7 @@ void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 }  // namespace RichardsComponentTransport
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index b8ebdab6e0e4e7babc924e672f63a082c7a54058..d6743efd4370002101518908897ce40db6baad17 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -75,7 +75,7 @@ void RichardsFlowProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
@@ -94,7 +94,7 @@ void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 }  // namespace RichardsFlow
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index a8c5b7cdb92b8cbd60abc144ea4e328b6c91cb77..c7e64883de94864c307c8b0b2aed628e1eb17363 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -463,7 +463,7 @@ void RichardsMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 template <int DisplacementDim>
@@ -508,7 +508,7 @@ void RichardsMechanicsProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         if (_use_monolithic_scheme)
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
index d7ca5dce3b0aaa0dec84b3722759e43f9a432ec4..b2b04be822b964b98d532c8bf03507230029b523 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.cpp
@@ -212,7 +212,7 @@ void SmallDeformationProcess<DisplacementDim>::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 template <int DisplacementDim>
@@ -233,7 +233,7 @@ void SmallDeformationProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     transformVariableFromGlobalVector(b, 0, *_local_to_global_index_map,
                                       *_nodal_forces, std::negate<double>());
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
index 90b13c5812434ecf8ce11b90496ada44c44d5a2f..5b3db21887d70511abb53b9eb6a7ddfb4542d40a 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalProcess.cpp
@@ -247,7 +247,7 @@ void SmallDeformationNonlocalProcess<DisplacementDim>::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 template <int DisplacementDim>
@@ -287,7 +287,7 @@ void SmallDeformationNonlocalProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     b.copyValues(*_nodal_forces);
     std::transform(_nodal_forces->begin(), _nodal_forces->end(),
diff --git a/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusion.cpp b/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusion.cpp
index d8668153291f31ee6c8cbe9c84a0fe6979f6a5d0..c1ba6edd4b267ec43ff3f8eb058ad6b7472b94e8 100644
--- a/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusion.cpp
+++ b/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusion.cpp
@@ -71,7 +71,7 @@ void SteadyStateDiffusion::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void SteadyStateDiffusion::assembleWithJacobianConcreteProcess(
@@ -89,7 +89,7 @@ void SteadyStateDiffusion::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 }  // namespace SteadyStateDiffusion
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 67bbb6a9658ac10aadcaa175f3ae28ec760ebc6f..835d07eb3cfbac2c5c3516d8663d3e2f863fe82e 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -223,7 +223,7 @@ void TESProcess::assembleConcreteProcess(const double t, double const dt,
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void TESProcess::assembleWithJacobianConcreteProcess(
@@ -240,7 +240,7 @@ void TESProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 void TESProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
index 2b36ca9ed85e83adee9226727b6defd2016bfc2a..576355aa0bd871faffd8b9a909aacf1487f8f81e 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
@@ -88,7 +88,7 @@ void ThermalTwoPhaseFlowWithPPProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void ThermalTwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
@@ -107,7 +107,7 @@ void ThermalTwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 void ThermalTwoPhaseFlowWithPPProcess::preTimestepConcreteProcess(
     std::vector<GlobalVector*> const& x, double const t, double const delta_t,
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
index 769e0cd06aa5d11c823c442eb1f8eb234928c37e..33d4a079f9c92573a85f5de92bb9e94db1ff908f 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
@@ -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, xdot, process_id, M, K, b, _coupled_solutions);
+        dof_table, t, dt, x, xdot, process_id, M, K, b);
 }
 
 template <int DisplacementDim>
@@ -304,7 +304,7 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, dof_tables, t, dt, x, xdot, dxdot_dx, dx_dx,
-        process_id, M, K, b, Jac, _coupled_solutions);
+        process_id, M, K, b, Jac);
 
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
         if (_use_monolithic_scheme)
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
index 390f01e96720e7795ca4907e4744798f6acf448d..99ab6af2b926cc4f974c2a8acbc40540c43bdecb 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldProcess.cpp
@@ -204,7 +204,7 @@ void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 template <int DisplacementDim>
@@ -254,7 +254,7 @@ void ThermoMechanicalPhaseFieldProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 template <int DisplacementDim>
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
index 91c2df8034f48829faa03a28f6da2bcd0ff24b09..b37b060e9395715349ea06e6ec032e3db0509c40 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.cpp
@@ -300,7 +300,7 @@ void ThermoMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 template <int DisplacementDim>
@@ -362,7 +362,7 @@ void ThermoMechanicsProcess<DisplacementDim>::
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 
     // TODO (naumov): Refactor the copy rhs part. This is copy from HM.
     auto copyRhs = [&](int const variable_id, auto& output_vector) {
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 330cd43408e66ff23c4b0d4a42eb4ac0590e8529..486c490e069d1310d4a39247e62e6e2f4ee2264b 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -770,6 +770,10 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
             BaseLib::RunTime time_timestep_process;
             time_timestep_process.start();
 
+            // The following setting of coupled_solutions can be removed only if
+            // the CoupledSolutionsForStaggeredScheme and related functions are
+            // removed totally from the computation of the secondary variable
+            // and from post-time functions.
             CoupledSolutionsForStaggeredScheme coupled_solutions(
                 _process_solutions);
 
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 878bd2074064c28381a73c14bca9c1349f49a7bf..b0cbbe940ec91acd26a669db95097fa8900b00d2 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -78,7 +78,7 @@ void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
@@ -97,7 +97,7 @@ void TwoPhaseFlowWithPPProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 
 }  // namespace TwoPhaseFlowWithPP
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
index 275df38732421bc5ff8e45643fc86a4e7d089f05..a40a2acb8d9cc5ed6ab410eb25dd705d943f6ef7 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -79,7 +79,7 @@ void TwoPhaseFlowWithPrhoProcess::assembleConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
         pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
-        b, _coupled_solutions);
+        b);
 }
 
 void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
@@ -98,7 +98,7 @@ void TwoPhaseFlowWithPrhoProcess::assembleWithJacobianConcreteProcess(
     GlobalExecutor::executeSelectedMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
         _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
-        dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
+        dxdot_dx, dx_dx, process_id, M, K, b, Jac);
 }
 void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
     std::vector<GlobalVector*> const& x, double const t, double const dt,
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 256cfd7ad6b40bbae640348d547979dea1cd259a..6f5d608486963313ff6fcd0922e5a5cffe96657a 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -45,8 +45,7 @@ void VectorMatrixAssembler::assemble(
         dof_tables,
     const double t, double const dt, std::vector<GlobalVector*> const& x,
     std::vector<GlobalVector*> const& xdot, int const process_id,
-    GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-    CoupledSolutionsForStaggeredScheme const* const /*cpl_xs*/)
+    GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
 {
     std::vector<std::vector<GlobalIndexType>> indices_of_processes;
     indices_of_processes.reserve(dof_tables.size());
@@ -79,18 +78,6 @@ void VectorMatrixAssembler::assemble(
             getCoupledLocalSolutions(xdot, indices_of_processes);
         auto const local_xdot = MathLib::toVector(local_coupled_xdots);
 
-        /*
-        auto local_coupled_xs0 = getCoupledLocalSolutions(cpl_xs->coupled_xs_t0,
-                                                          indices_of_processes);
-
-        auto local_coupled_xs =
-            getCoupledLocalSolutions(x, indices_of_processes);
-
-        auto const local_x = MathLib::toVector(local_coupled_xs);
-
-        ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-            std::move(local_coupled_xs0));
-        */
 
         local_assembler.assembleForStaggeredScheme(
             t, dt, local_x, local_xdot, process_id, _local_M_data,
@@ -125,8 +112,7 @@ void VectorMatrixAssembler::assembleWithJacobian(
     const double t, double const dt, std::vector<GlobalVector*> const& x,
     std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
     const double dx_dx, int const process_id, GlobalMatrix& M, GlobalMatrix& K,
-    GlobalVector& b, GlobalMatrix& Jac,
-    CoupledSolutionsForStaggeredScheme const* const /*cpl_xs*/)
+    GlobalVector& b, GlobalMatrix& Jac)
 {
     std::vector<std::vector<GlobalIndexType>> indices_of_processes;
     indices_of_processes.reserve(dof_tables.size());
@@ -162,18 +148,6 @@ void VectorMatrixAssembler::assembleWithJacobian(
             getCoupledLocalSolutions(xdot, indices_of_processes);
         auto const local_xdot = MathLib::toVector(local_coupled_xdots);
 
-        /*
-        auto local_coupled_xs0 = getCoupledLocalSolutions(cpl_xs->coupled_xs_t0,
-                                                          indices_of_processes);
-
-        auto local_coupled_xs =
-            getCoupledLocalSolutions(x, indices_of_processes);
-
-        auto const local_x = MathLib::toVector(local_coupled_xs);
-
-        ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-            std::move(local_coupled_xs0));
-        */
 
         _jacobian_assembler->assembleWithJacobianForStaggeredScheme(
             local_assembler, t, dt, local_x, local_xdot, dxdot_dx, dx_dx,
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index 1046d9f336329ade6e80db9dcc217fbbdcbf3dde..f9a8ea9f99f8a58e59eaf65e07f1da1898176996 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -50,8 +50,7 @@ public:
                   double const t, double const dt,
                   std::vector<GlobalVector*> const& x,
                   std::vector<GlobalVector*> const& xdot, int const process_id,
-                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-                  CoupledSolutionsForStaggeredScheme const* const cpl_xs);
+                  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b);
 
     //! Assembles \c M, \c K, \c b, and the Jacobian \c Jac of the residual.
     //! \note The Jacobian must be assembled.
@@ -64,8 +63,7 @@ public:
         const double t, double const dt, std::vector<GlobalVector*> const& x,
         std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
         const double dx_dx, int const process_id, GlobalMatrix& M,
-        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac,
-        CoupledSolutionsForStaggeredScheme const* const cpl_xs);
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac);
 
 private:
     // temporary data only stored here in order to avoid frequent memory