From 2ffa66cff21bf404fa37f3db0c12ce379c35dbf7 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 12 Oct 2017 17:38:56 +0200
Subject: [PATCH] [Coupling] Removed the coupled process members

 from the two classes of coupled term.

[Coupling] Changed the members of the coupled term,

 which only contains the solutions of the coupled equations.
---
 .../coupled_processes/i_coupled_processes.md  |   1 -
 .../coupled_processes/t_coupled_process.md    |   1 -
 ProcessLib/AbstractJacobianAssembler.h        |   1 -
 .../ComponentTransportProcess.cpp             |   5 +
 .../CoupledSolutionsForStaggeredScheme.cpp    |  56 +++++----
 .../CoupledSolutionsForStaggeredScheme.h      |  69 +++++------
 ProcessLib/HT/HTProcess.cpp                   |  31 +++++
 ProcessLib/HT/HTProcess.h                     |  15 +++
 .../HeatConduction/HeatConductionFEM-impl.h   |   3 +-
 ProcessLib/HeatConduction/HeatConductionFEM.h |  86 --------------
 .../HeatConduction/HeatConductionProcess.cpp  |  17 ---
 .../HeatConduction/HeatConductionProcess.h    |  10 --
 .../LiquidFlowLocalAssembler-impl.h           |  17 ++-
 .../LiquidFlow/LiquidFlowLocalAssembler.h     |  23 ++--
 ProcessLib/LiquidFlow/LiquidFlowProcess.cpp   |  10 +-
 ProcessLib/LiquidFlow/LiquidFlowProcess.h     |   2 +-
 ProcessLib/LocalAssemblerInterface.cpp        |  36 +++---
 ProcessLib/LocalAssemblerInterface.h          |  22 ++--
 ProcessLib/Process.cpp                        |  16 ++-
 ProcessLib/Process.h                          |  10 +-
 ProcessLib/UncoupledProcessesTimeLoop.cpp     | 107 ++++--------------
 ProcessLib/UncoupledProcessesTimeLoop.h       |   8 +-
 ProcessLib/VectorMatrixAssembler.cpp          | 105 +++++------------
 ProcessLib/VectorMatrixAssembler.h            |   6 +-
 24 files changed, 241 insertions(+), 416 deletions(-)
 delete mode 100644 Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/i_coupled_processes.md
 delete mode 100644 Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/t_coupled_process.md

diff --git a/Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/i_coupled_processes.md b/Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/i_coupled_processes.md
deleted file mode 100644
index 6550f6798a1..00000000000
--- a/Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/i_coupled_processes.md
+++ /dev/null
@@ -1 +0,0 @@
-Defines the coupled processes to a process.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/t_coupled_process.md b/Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/t_coupled_process.md
deleted file mode 100644
index 98b97648e20..00000000000
--- a/Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/t_coupled_process.md
+++ /dev/null
@@ -1 +0,0 @@
-Specifies a single coupled process to a process by name.
\ No newline at end of file
diff --git a/ProcessLib/AbstractJacobianAssembler.h b/ProcessLib/AbstractJacobianAssembler.h
index 50dde014d8a..8645a26a952 100644
--- a/ProcessLib/AbstractJacobianAssembler.h
+++ b/ProcessLib/AbstractJacobianAssembler.h
@@ -34,7 +34,6 @@ public:
     //! \f$b\f$ with coupling.
     virtual void assembleWithJacobianAndCoupling(
         LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
-        std::vector<double> const& /*local_x*/,
         std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
         const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
         std::vector<double>& /*local_K_data*/,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index bf446dacc27..68bd1a4a59d 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -59,6 +59,8 @@ void ComponentTransportProcess::assembleConcreteProcess(
     GlobalVector& b)
 {
     DBUG("Assemble ComponentTransportProcess.");
+    if (!_is_monolithic_scheme)
+        setCoupledSolutionsOfPreviousTimeStep();
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
@@ -73,6 +75,9 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian ComponentTransportProcess.");
 
+    if (!_is_monolithic_scheme)
+        setCoupledSolutionsOfPreviousTimeStep();
+
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
index 1f7eecb535d..590af56f026 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.cpp
@@ -18,48 +18,46 @@
 namespace ProcessLib
 {
 CoupledSolutionsForStaggeredScheme::CoupledSolutionsForStaggeredScheme(
-    std::unordered_map<std::type_index, Process const&> const&
-        coupled_processes_,
-    std::unordered_map<std::type_index, GlobalVector const&> const& coupled_xs_,
-    const double dt_)
-    : coupled_processes(coupled_processes_), coupled_xs(coupled_xs_), dt(dt_)
+    std::vector<std::reference_wrapper<GlobalVector const>> const& coupled_xs_,
+    const double dt_, const int process_id_)
+    : coupled_xs(coupled_xs_), dt(dt_), process_id(process_id_)
 {
-    for (auto const& coupled_x_pair : coupled_xs)
+    for (auto const& coupled_x : coupled_xs)
     {
-        auto const& coupled_x = coupled_x_pair.second;
-        MathLib::LinAlg::setLocalAccessibleVector(coupled_x);
+        MathLib::LinAlg::setLocalAccessibleVector(coupled_x.get());
     }
+}
+
+std::vector<std::vector<double>> getPreviousLocalSolutions(
+    const CoupledSolutionsForStaggeredScheme& cpl_xs,
+    const std::vector<GlobalIndexType>& indices)
+{
+    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);
 
-    for (auto const& coupled_process_pair : coupled_processes)
+    for (std::size_t i = 0; i < number_of_coupled_solutions; i++)
     {
-        auto const& coupled_pcs = coupled_process_pair.second;
-        auto const prevous_time_x = coupled_pcs.getPreviousTimeStepSolution();
-        if (prevous_time_x)
-        {
-            MathLib::LinAlg::setLocalAccessibleVector(*prevous_time_x);
-        }
+        auto const& x_t0 = *cpl_xs.coupled_xs_t0[i];
+        local_xs_t0.emplace_back(x_t0.get(indices));
     }
+    return local_xs_t0;
 }
 
-std::unordered_map<std::type_index, const std::vector<double>>
-getCurrentLocalSolutionsOfCoupledProcesses(
-    const std::unordered_map<std::type_index, GlobalVector const&>&
-        global_coupled_xs,
+std::vector<std::vector<double>> getCurrentLocalSolutions(
+    const CoupledSolutionsForStaggeredScheme& cpl_xs,
     const std::vector<GlobalIndexType>& indices)
 {
-    std::unordered_map<std::type_index, const std::vector<double>>
-        local_coupled_xs;
+    const auto number_of_coupled_solutions = cpl_xs.coupled_xs.size();
+    std::vector<std::vector<double>> local_xs_t1;
+    local_xs_t1.reserve(number_of_coupled_solutions);
 
-    // Get local nodal solutions of the coupled equations.
-    for (auto const& global_coupled_x_pair : global_coupled_xs)
+    for (std::size_t i = 0; i < number_of_coupled_solutions; i++)
     {
-        auto const& coupled_x = global_coupled_x_pair.second;
-        auto const local_coupled_x = coupled_x.get(indices);
-        BaseLib::insertIfTypeIndexKeyUniqueElseError(
-            local_coupled_xs, global_coupled_x_pair.first, local_coupled_x,
-            "local_coupled_x");
+        auto const& x_t1 = cpl_xs.coupled_xs[i].get();
+        local_xs_t1.emplace_back(x_t1.get(indices));
     }
-    return local_coupled_xs;
+    return local_xs_t1;
 }
 
 }  // end of ProcessLib
diff --git a/ProcessLib/CoupledSolutionsForStaggeredScheme.h b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
index c2880de217d..9a940d1ebc7 100644
--- a/ProcessLib/CoupledSolutionsForStaggeredScheme.h
+++ b/ProcessLib/CoupledSolutionsForStaggeredScheme.h
@@ -12,18 +12,17 @@
 
 #pragma once
 
-#include <unordered_map>
-#include <typeindex>
+#include <functional>
+#include <utility>
+#include <vector>
 
 #include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
 
 namespace ProcessLib
 {
-class Process;
-
 /**
- *  A struct to keep the references of the coupled processes and the references
- *  of the current solutions of the equations of the coupled processes.
+ *  A struct to keep the references of the current solutions of the equations of
+ *  the coupled processes.
  *
  *  During staggered coupling iteration, an instance of this struct is created
  *  and passed through interfaces to global and local assemblers for each
@@ -32,68 +31,54 @@ class Process;
 struct CoupledSolutionsForStaggeredScheme
 {
     CoupledSolutionsForStaggeredScheme(
-        std::unordered_map<std::type_index, Process const&> const&
-            coupled_processes_,
-        std::unordered_map<std::type_index, GlobalVector const&> const&
+        std::vector<std::reference_wrapper<GlobalVector const>> const&
             coupled_xs_,
-        const double dt_);
-
-    /// References to the coupled processes are distinguished by the keys of
-    /// process types.
-    std::unordered_map<std::type_index, Process const&> const&
-        coupled_processes;
+        const double dt_, const int process_id_);
 
     /// References to the current solutions of the coupled processes.
-    /// The coupled solutions are distinguished by the keys of process types.
-    std::unordered_map<std::type_index, GlobalVector const&> const& coupled_xs;
+    std::vector<std::reference_wrapper<GlobalVector const>> const& coupled_xs;
+
+    /// Pointers to the vector of the solutions of the previous time step.
+    std::vector<GlobalVector*> coupled_xs_t0;
 
     const double dt;  ///< Time step size.
+    const int process_id;
 };
 
 /**
- *  A struct to keep the references of the coupled processes, and local element
- *  solutions  of the current and previous time step solutions of the equations
- *  of the coupled processes.
+ *  A struct to keep the references to the local element solutions  of the
+ *  current and previous time step solutions of the equations of the coupled
+ *  processes.
  *
  *  During the global assembly loop, an instance of this struct is created for
  *  each element and it is then passed to local assemblers.
  */
 struct LocalCoupledSolutions
 {
-    LocalCoupledSolutions(
-        const double dt_,
-        std::unordered_map<std::type_index, Process const&> const&
-            coupled_processes_,
-        std::unordered_map<std::type_index, const std::vector<double>>&&
-            local_coupled_xs0_,
-        std::unordered_map<std::type_index, const std::vector<double>>&&
-            local_coupled_xs_)
+    LocalCoupledSolutions(const double dt_, const int process_id_,
+                          std::vector<std::vector<double>>&& local_coupled_xs0_,
+                          std::vector<std::vector<double>>&& local_coupled_xs_)
         : dt(dt_),
-          coupled_processes(coupled_processes_),
+          process_id(process_id_),
           local_coupled_xs0(std::move(local_coupled_xs0_)),
           local_coupled_xs(std::move(local_coupled_xs_))
     {
     }
 
     const double dt;  ///< Time step size.
-
-    /// References to the coupled processes are distinguished by the keys of
-    /// process types.
-    std::unordered_map<std::type_index, Process const&> const&
-        coupled_processes;
+    const int process_id;
 
     /// Local solutions of the previous time step.
-    std::unordered_map<std::type_index, const std::vector<double>> const
-        local_coupled_xs0;
+    std::vector<std::vector<double>> const local_coupled_xs0;
     /// Local solutions of the current time step.
-    std::unordered_map<std::type_index, const std::vector<double>> const
-        local_coupled_xs;
+    std::vector<std::vector<double>> const local_coupled_xs;
 };
 
-std::unordered_map<std::type_index, const std::vector<double>>
-getCurrentLocalSolutionsOfCoupledProcesses(
-    const std::unordered_map<std::type_index, GlobalVector const&>&
-        global_coupled_xs,
+std::vector<std::vector<double>> getPreviousLocalSolutions(
+    const CoupledSolutionsForStaggeredScheme& cpl_xs,
     const std::vector<GlobalIndexType>& indices);
 
+std::vector<std::vector<double>> getCurrentLocalSolutions(
+    const CoupledSolutionsForStaggeredScheme& cpl_xs,
+    const std::vector<GlobalIndexType>& indices);
 }  // end of ProcessLib
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 3b95d1e54b5..5c634c84e82 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -73,6 +73,9 @@ void HTProcess::assembleConcreteProcess(const double t,
 {
     DBUG("Assemble HTProcess.");
 
+    if (!_is_monolithic_scheme)
+        setCoupledSolutionsOfPreviousTimeStep();
+
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
@@ -86,6 +89,9 @@ void HTProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian HTProcess.");
 
+    if (!_is_monolithic_scheme)
+        setCoupledSolutionsOfPreviousTimeStep();
+
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
@@ -93,6 +99,31 @@ void HTProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
+void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
+                                            const double /*t*/,
+                                            const double /*delta_t*/,
+                                            const int process_id)
+{
+    assert(process_id < 2);
+
+    if (_is_monolithic_scheme)
+        return;
+
+    if (!_xs_previous_timestep[process_id])
+    {
+        _xs_previous_timestep[process_id] =
+            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+    }
+    else
+    {
+        auto& x0 = *_xs_previous_timestep[process_id];
+        MathLib::LinAlg::copy(x, x0);
+    }
+    
+    auto& x0 = *_xs_previous_timestep[process_id];
+    MathLib::LinAlg::setLocalAccessibleVector(x0);
+}
+
 }  // namespace HT
 }  // namespace ProcessLib
 
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 2cf76cc7964..c1d2bb5da0d 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -58,6 +58,13 @@ public:
     bool isLinear() const override { return false; }
     //! @}
 
+    // Get the solution of the previous time step.
+    GlobalVector* getPreviousTimeStepSolution(
+        const int process_id) const override
+    {
+        return _xs_previous_timestep[process_id].get();
+    }
+
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -73,9 +80,17 @@ private:
         const double dxdot_dx, const double dx_dx, 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;
+
     const std::unique_ptr<HTMaterialProperties> _material_properties;
 
     std::vector<std::unique_ptr<HTLocalAssemblerInterface>> _local_assemblers;
+    
+    /// Solutions of the previous time step
+    std::array<std::unique_ptr<GlobalVector>, 2> _xs_previous_timestep;
+
 };
 
 }  // namespace HT
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM-impl.h b/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
index 3367e05f4b9..8c21f2d0364 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
@@ -112,6 +112,7 @@ void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
                             std::vector<double>& /*local_b_data*/,
                             LocalCoupledSolutions const& coupled_term)
 {
+    /*
     for (auto const& coupled_process_pair : coupled_term.coupled_processes)
     {
         if (coupled_process_pair.first ==
@@ -159,7 +160,7 @@ void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
                 "This coupled process is not presented for "
                 "HeatConduction process");
         }
-    }
+    }*/
 }
 
 }  // namespace HeatConduction
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index 3d61130db91..3f8d58c5f3c 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -21,10 +21,6 @@
 #include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Utils/InitShapeMatrices.h"
 
-// For coupling
-#include "ProcessLib/LiquidFlow/LiquidFlowProcess.h"
-#include "ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h"
-
 namespace ProcessLib
 {
 namespace HeatConduction
@@ -131,12 +127,6 @@ public:
         }
     }
 
-    void assembleWithCoupledTerm(
-        double const t, std::vector<double> const& local_x,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
-        std::vector<double>& /*local_b_data*/,
-        LocalCoupledSolutions const& coupled_term) override;
-
     void computeSecondaryVariableConcrete(
         const double t, std::vector<double> const& local_x) override
     {
@@ -216,83 +206,7 @@ private:
         _shape_matrices;
 
     std::vector<std::vector<double>> _heat_fluxes;
-
-    /**
-     * @brief Assemble local matrices and vectors of the equations of
-     *        heat transport process in porous media with full saturated liquid
-     *        flow.
-     *
-     * @param t                          Time.
-     * @param material_id                Material ID of the element.
-     * @param pos                        Position (t, x) of the element.
-     * @param gravitational_axis_id      The ID of the axis, which indicates the
-     *                                   direction of gravity portion of the
-     *                                   Darcy's velocity.
-     * @param gravitational_acceleration Gravitational acceleration, 9.81 in the
-     *                                   SI unit standard.
-     * @param permeability               Intrinsic permeability of liquid
-     *                                   in porous media.
-     * @param liquid_flow_properties     Liquid flow properties.
-     * @param local_x                    Local vector of unknowns, e.g.
-     *                                   nodal temperatures of the element.
-     * @param local_p                    Local vector of nodal pore pressures.
-     * @param local_M_data               Local mass matrix.
-     * @param local_K_data               Local Laplace matrix.
-     */
-    template <typename LiquidFlowVelocityCalculator>
-    void assembleHeatTransportLiquidFlow(
-        double const t, int const material_id, SpatialPosition& pos,
-        int const gravitational_axis_id,
-        double const gravitational_acceleration,
-        Eigen::MatrixXd const& permeability,
-        ProcessLib::LiquidFlow::LiquidFlowMaterialProperties const&
-            liquid_flow_properties,
-        std::vector<double> const& local_x, std::vector<double> const& local_p,
-        std::vector<double>& local_M_data, std::vector<double>& local_K_data);
-
-    /// Calculator of liquid flow velocity for isotropic permeability
-    /// tensor
-    struct IsotropicLiquidFlowVelocityCalculator
-    {
-        static GlobalDimVectorType calculateVelocity(
-            Eigen::Map<const NodalVectorType> const& local_p,
-            ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
-            double const mu, double const rho_g,
-            int const gravitational_axis_id)
-        {
-            const double K = permeability(0, 0) / mu;
-            // Compute the velocity
-            GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p;
-            // gravity term
-            if (gravitational_axis_id >= 0)
-                darcy_velocity[gravitational_axis_id] -= K * rho_g;
-            return darcy_velocity;
-        }
-    };
-
-    /// Calculator of liquid flow velocity for anisotropic permeability
-    /// tensor
-    struct AnisotropicLiquidFlowVelocityCalculator
-    {
-        static GlobalDimVectorType calculateVelocity(
-            Eigen::Map<const NodalVectorType> const& local_p,
-            ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
-            double const mu, double const rho_g,
-            int const gravitational_axis_id)
-        {
-            GlobalDimVectorType darcy_velocity =
-                -permeability * sm.dNdx * local_p / mu;
-            if (gravitational_axis_id >= 0)
-            {
-                darcy_velocity.noalias() -=
-                    rho_g * permeability.col(gravitational_axis_id) / mu;
-            }
-            return darcy_velocity;
-        }
-    };
 };
 
 }  // namespace HeatConduction
 }  // namespace ProcessLib
-
-#include "HeatConductionFEM-impl.h"
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 986aeb31e46..4644395a05f 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -33,23 +33,6 @@ HeatConductionProcess::HeatConductionProcess(
 {
 }
 
-void HeatConductionProcess::preTimestepConcreteProcess(GlobalVector const& x,
-                                            const double /*t*/,
-                                            const double /*delta_t*/,
-                                            const int /*process_id*/)
-{
-    if (!_x_previous_timestep)
-    {
-        _x_previous_timestep =
-            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
-    }
-    else
-    {
-        auto& x0 = *_x_previous_timestep;
-        MathLib::LinAlg::copy(x, x0);
-    }
-}
-
 void HeatConductionProcess::initializeConcreteProcess(
     NumLib::LocalToGlobalIndexMap const& dof_table,
     MeshLib::Mesh const& mesh,
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index 3965ef37f5e..e81e93a9ab5 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -41,16 +41,6 @@ public:
     void computeSecondaryVariableConcrete(
         double const t, GlobalVector const& x) override;
 
-    void preTimestepConcreteProcess(GlobalVector const& x, const double t,
-                                    const double delta_t,
-                                    const int process_id) override;
-
-    // Get the solution of the previous time step.
-    GlobalVector* getPreviousTimeStepSolution() const override
-    {
-        return _x_previous_timestep.get();
-    }
-
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 6e6c923c5dd..2d6262194aa 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -53,6 +53,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
             pos, permeability);
 }
 
+/*
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
@@ -60,7 +61,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
                             std::vector<double>& local_M_data,
                             std::vector<double>& local_K_data,
                             std::vector<double>& local_b_data,
-                            LocalCoupledSolutions const& coupled_term)
+                            LocalCouplingTerm const& coupled_term)
 {
     SpatialPosition pos;
     pos.setElementID(_element.getID());
@@ -104,6 +105,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         }
     }
 }
+*/
 
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
@@ -166,6 +168,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     }
 }
 
+/*
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 template <typename LaplacianGravityVelocityCalculator>
@@ -240,6 +243,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
             eff_thermal_expansion * (T - T0) * integration_factor * sm.N / dt;
     }
 }
+*/
 
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
@@ -268,22 +272,23 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     //  the assert must be changed to perm.rows() == _element->getDimension()
     assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
 
-    if (!_coupled_solutions)
+    //if (!_coupling_term)
     {
         computeDarcyVelocity(permeability, local_x, veloctiy_cache_vectors);
     }
+    /*
     else
     {
         auto const local_coupled_xs =
             getCurrentLocalSolutionsOfCoupledProcesses(
-                _coupled_solutions->coupled_xs, indices);
+                _coupling_term->coupled_xs, indices);
         if (local_coupled_xs.empty())
             computeDarcyVelocity(permeability, local_x, veloctiy_cache_vectors);
         else
             computeDarcyVelocityWithCoupling(permeability, local_x,
                                              local_coupled_xs,
                                              veloctiy_cache_vectors);
-    }
+    }*/
 
     return veloctiy_cache;
 }
@@ -341,6 +346,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     }
 }
 
+/*
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
@@ -362,6 +368,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
                                    darcy_velocity_at_ips);
 }
 
+
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 template <typename LaplacianGravityVelocityCalculator>
@@ -399,7 +406,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
             _gravitational_axis_id, darcy_velocity_at_ips);
     }
 }
-
+*/
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index a5b29cf9281..b7ea92d5cfd 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -42,11 +42,7 @@ class LiquidFlowLocalAssemblerInterface
       public NumLib::ExtrapolatableElement
 {
 public:
-    LiquidFlowLocalAssemblerInterface(
-        CoupledSolutionsForStaggeredScheme* const coupled_solutions)
-        : _coupled_solutions(coupled_solutions)
-    {
-    }
+    LiquidFlowLocalAssemblerInterface() = default;
 
     void setCoupledSolutionsForStaggeredScheme(
         std::size_t const /*mesh_item_id*/,
@@ -92,9 +88,8 @@ public:
         int const gravitational_axis_id,
         double const gravitational_acceleration,
         double const reference_temperature,
-        LiquidFlowMaterialProperties const& material_propertries,
-        CoupledSolutionsForStaggeredScheme* coupled_solutions)
-        : LiquidFlowLocalAssemblerInterface(coupled_solutions),
+        LiquidFlowMaterialProperties const& material_propertries)
+        : LiquidFlowLocalAssemblerInterface(),
           _element(element),
           _integration_method(integration_order),
           _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
@@ -112,11 +107,12 @@ public:
                   std::vector<double>& local_K_data,
                   std::vector<double>& local_b_data) override;
 
+    /*
     void assembleWithCoupledTerm(
         double const t, std::vector<double> const& local_x,
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data,
-        LocalCoupledSolutions const& coupled_term) override;
+        LocalCouplingTerm const& coupled_term) override;*/
 
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
@@ -191,24 +187,26 @@ private:
                                  SpatialPosition const& pos,
                                  Eigen::MatrixXd const& permeability);
 
+    /*
     template <typename LaplacianGravityVelocityCalculator>
     void assembleWithCoupledWithHeatTransport(
         const int material_id, double const t, double const dt,
         std::vector<double> const& local_x, std::vector<double> const& local_T0,
         std::vector<double> const& local_T1, std::vector<double>& local_M_data,
         std::vector<double>& local_K_data, std::vector<double>& local_b_data,
-        SpatialPosition const& pos, Eigen::MatrixXd const& permeability);
+        SpatialPosition const& pos, Eigen::MatrixXd const& permeability);*/
 
     void computeDarcyVelocity(
         Eigen::MatrixXd const& permeability,
         std::vector<double> const& local_x,
         MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
 
+    /*
     void computeDarcyVelocityWithCoupling(
         Eigen::MatrixXd const& permeability, std::vector<double> const& local_x,
         std::unordered_map<std::type_index, const std::vector<double>> const&
             coupled_local_solutions,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
+        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;*/
 
     template <typename LaplacianGravityVelocityCalculator>
     void computeDarcyVelocityLocal(
@@ -216,12 +214,13 @@ private:
         Eigen::MatrixXd const& permeability,
         MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
 
+    /*
     template <typename LaplacianGravityVelocityCalculator>
     void computeDarcyVelocityCoupledWithHeatTransportLocal(
         std::vector<double> const& local_x,
         std::vector<double> const& local_T,
         Eigen::MatrixXd const& permeability,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
+        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;*/
 
     const int _gravitational_axis_id;
     const double _gravitational_acceleration;
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 2d2d0740801..87f3b55ec60 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -66,7 +66,7 @@ void LiquidFlowProcess::initializeConcreteProcess(
         pv.getShapeFunctionOrder(), _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _gravitational_axis_id,
         _gravitational_acceleration, _reference_temperature,
-        *_material_properties, _coupled_solutions);
+        *_material_properties);
 
     _secondary_variables.addSecondaryVariable(
         "darcy_velocity",
@@ -112,13 +112,13 @@ void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
         _coupled_solutions);
 }
 
-void LiquidFlowProcess::setCoupledSolutionsForStaggeredSchemeToLocalAssemblers()
+void LiquidFlowProcess::setStaggeredCouplingTermToLocalAssemblers()
 {
     DBUG("Compute the velocity for LiquidFlowProcess.");
+    /*
     GlobalExecutor::executeMemberOnDereferenced(
-        &LiquidFlowLocalAssemblerInterface::
-            setCoupledSolutionsForStaggeredScheme,
-        _local_assemblers, _coupled_solutions);
+        &LiquidFlowLocalAssemblerInterface::setStaggeredCouplingTerm,
+        _local_assemblers, _coupled_solutions);*/
 }
 
 }  // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index 1e6fca223ef..0ab5f4ac69b 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -88,7 +88,7 @@ public:
         return _material_properties.get();
     }
 
-    void setCoupledSolutionsForStaggeredSchemeToLocalAssemblers() override;
+    void setStaggeredCouplingTermToLocalAssemblers() override;
 
 private:
     void initializeConcreteProcess(
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 910a73b311a..1fb8dc5081f 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -15,8 +15,18 @@
 
 namespace ProcessLib
 {
+void LocalAssemblerInterface::assemble(double const /*t*/,
+                                       std::vector<double> const& /*local_x*/,
+                                       std::vector<double>& /*local_M_data*/,
+                                       std::vector<double>& /*local_K_data*/,
+                                       std::vector<double>& /*local_b_data*/)
+{
+    OGS_FATAL(
+        "The assemble() function is not implemented in the local assembler.");
+}
+
 void LocalAssemblerInterface::assembleWithCoupledTerm(
-    double const /*t*/, std::vector<double> const& /*local_x*/,
+    double const /*t*/,
     std::vector<double>& /*local_M_data*/,
     std::vector<double>& /*local_K_data*/,
     std::vector<double>& /*local_b_data*/,
@@ -41,13 +51,13 @@ void LocalAssemblerInterface::assembleWithJacobian(
 }
 
 void LocalAssemblerInterface::assembleWithJacobianAndCoupling(
-    double const /*t*/, std::vector<double> const& /*local_x*/,
-    std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
-    const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
+    double const /*t*/, std::vector<double> const& /*local_xdot*/,
+    const double /*dxdot_dx*/, const double /*dx_dx*/,
+    std::vector<double>& /*local_M_data*/,
     std::vector<double>& /*local_K_data*/,
     std::vector<double>& /*local_b_data*/,
     std::vector<double>& /*local_Jac_data*/,
-    LocalCoupledSolutions const& /*coupled_solutions*/)
+    LocalCoupledSolutions const& /*local_coupled_solutions*/)
 {
     OGS_FATAL(
         "The assembleWithJacobianAndCoupling() function is not implemented in"
@@ -57,26 +67,20 @@ void LocalAssemblerInterface::assembleWithJacobianAndCoupling(
 void LocalAssemblerInterface::computeSecondaryVariable(
     std::size_t const mesh_item_id,
     NumLib::LocalToGlobalIndexMap const& dof_table, double const t,
-    GlobalVector const& x,
-    CoupledSolutionsForStaggeredScheme const* coupled_term)
+    GlobalVector const& x, CoupledSolutionsForStaggeredScheme const* coupled_xs)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
-    auto const local_x = x.get(indices);
 
-    if (!coupled_term)
+    if (!coupled_xs)
     {
+        auto const local_x = x.get(indices);
         computeSecondaryVariableConcrete(t, local_x);
     }
     else
     {
         auto const local_coupled_xs =
-            getCurrentLocalSolutionsOfCoupledProcesses(coupled_term->coupled_xs,
-                                                       indices);
-        if (!local_coupled_xs.empty())
-            computeSecondaryVariableWithCoupledProcessConcrete(
-                t, local_x, local_coupled_xs);
-        else
-            computeSecondaryVariableConcrete(t, local_x);
+            getCurrentLocalSolutions(*coupled_xs, indices);
+        computeSecondaryVariableWithCoupledProcessConcrete(t, local_coupled_xs);
     }
 }
 
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index e3fc5cb7983..6769abcaa40 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -9,9 +9,6 @@
 
 #pragma once
 
-#include <unordered_map>
-#include <typeindex>
-
 #include "NumLib/NumericsConfig.h"
 #include "MathLib/Point3d.h"
 
@@ -41,11 +38,10 @@ public:
     virtual void assemble(double const t, std::vector<double> const& local_x,
                           std::vector<double>& local_M_data,
                           std::vector<double>& local_K_data,
-                          std::vector<double>& local_b_data) = 0;
+                          std::vector<double>& local_b_data);
 
     virtual void assembleWithCoupledTerm(
         double const t,
-        std::vector<double> const& local_x,
         std::vector<double>& local_M_data,
         std::vector<double>& local_K_data,
         std::vector<double>& local_b_data,
@@ -61,18 +57,17 @@ public:
                                       std::vector<double>& local_Jac_data);
 
     virtual void assembleWithJacobianAndCoupling(
-        double const t, std::vector<double> const& local_x,
-        std::vector<double> const& local_xdot, const double dxdot_dx,
-        const double dx_dx, std::vector<double>& local_M_data,
-        std::vector<double>& local_K_data, std::vector<double>& local_b_data,
-        std::vector<double>& local_Jac_data,
-        LocalCoupledSolutions const& coupled_solutions);
+        double const t, std::vector<double> const& local_xdot,
+        const double dxdot_dx, const double dx_dx,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
+        LocalCoupledSolutions const& local_coupled_solutions);
 
     virtual void computeSecondaryVariable(
         std::size_t const mesh_item_id,
         NumLib::LocalToGlobalIndexMap const& dof_table, const double t,
         GlobalVector const& x,
-        CoupledSolutionsForStaggeredScheme const* coupled_term);
+        CoupledSolutionsForStaggeredScheme const* coupled_xs);
 
     virtual void preTimestep(std::size_t const mesh_item_id,
                              NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -105,8 +100,7 @@ private:
     }
 
     virtual void computeSecondaryVariableWithCoupledProcessConcrete(
-        double const /*t*/, std::vector<double> const& /*local_x*/,
-        std::unordered_map<std::type_index, const std::vector<double>> const&
+        double const /*t*/, std::vector<std::vector<double>> const&
         /*coupled_local_solutions*/)
     {
     }
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 43003cedc30..874ac8d7ad0 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -7,6 +7,8 @@
  *
  */
 
+#include <bits/stl_vector.h>
+
 #include "Process.h"
 
 #include "BaseLib/Functional.h"
@@ -32,7 +34,7 @@ Process::Process(
       _named_function_caller(std::move(named_function_caller)),
       _global_assembler(std::move(jacobian_assembler)),
       _is_monolithic_scheme(true),
-      _coupling_solutions(nullptr),
+      _coupled_solutions(nullptr),
       _integration_order(integration_order),
       _process_variables(std::move(process_variables)),
       _boundary_conditions(parameters)
@@ -329,4 +331,16 @@ NumLib::IterationResult Process::postIteration(const GlobalVector& x)
     return postIterationConcreteProcess(x);
 }
 
+void Process::setCoupledSolutionsOfPreviousTimeStep()
+{
+    const auto number_of_coupled_solutions =
+        _coupled_solutions->coupled_xs.size();
+    _coupled_solutions->coupled_xs_t0.reserve(number_of_coupled_solutions);
+    for (std::size_t i = 0; i < number_of_coupled_solutions; i++)
+    {
+        _coupled_solutions->coupled_xs_t0.emplace_back(
+            getPreviousTimeStepSolution(i));
+    }
+}
+
 }  // namespace ProcessLib
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index e5d9bd3a3b6..5ca6992f8fa 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -119,7 +119,8 @@ public:
     }
 
     // Get the solution of the previous time step.
-    virtual GlobalVector* getPreviousTimeStepSolution() const
+    virtual GlobalVector* getPreviousTimeStepSolution(
+        const int /*process_id*/) const
     {
         return nullptr;
     }
@@ -220,10 +221,13 @@ protected:
     mutable bool _is_monolithic_scheme;
 
     /// Pointer to CoupledSolutionsForStaggeredScheme, which contains the
-    /// references to the coupled processes and the references to the
-    /// solutions of the coupled processes.
+    /// references to the solutions of the coupled processes.
     CoupledSolutionsForStaggeredScheme* _coupled_solutions;
 
+    /// Set the solutions of the previous time step to the coupled term.
+    /// It only performs for the staggered scheme.
+    void setCoupledSolutionsOfPreviousTimeStep();
+
     /// Order of the integration method for element-wise integration.
     /// The Gauss-Legendre integration method and available orders is
     /// implemented in MathLib::GaussLegendre.
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index b9d45047765..92f8a612967 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -18,6 +18,7 @@
 #include "NumLib/TimeStepping/CreateTimeStepper.h"
 
 #include "MathLib/LinAlg/LinAlg.h"
+#include "CoupledSolutionsForStaggeredScheme.h"
 
 std::unique_ptr<ProcessLib::Output> createOutput(
     BaseLib::ConfigTree const& config, std::string const& output_directory)
@@ -109,8 +110,6 @@ struct SingleProcessData
         std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
         std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
         Process& process_,
-        std::unordered_map<std::type_index, Process const&>&&
-            coupled_processes_,
         ProcessOutput&& process_output_);
 
     SingleProcessData(SingleProcessData&& spd);
@@ -136,8 +135,6 @@ struct SingleProcessData
     NumLib::InternalMatrixStorage* mat_strg = nullptr;
 
     Process& process;
-    /// Coupled processes.
-    std::unordered_map<std::type_index, Process const&> const coupled_processes;
     ProcessOutput process_output;
 };
 
@@ -148,7 +145,6 @@ SingleProcessData::SingleProcessData(
     std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
     std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
     Process& process_,
-    std::unordered_map<std::type_index, Process const&>&& coupled_processes_,
     ProcessOutput&& process_output_)
     : timestepper(std::move(timestepper_)),
       nonlinear_solver_tag(NLTag),
@@ -157,7 +153,6 @@ SingleProcessData::SingleProcessData(
       conv_crit(std::move(conv_crit_)),
       time_disc(std::move(time_disc_)),
       process(process_),
-      coupled_processes(coupled_processes_),
       process_output(std::move(process_output_))
 {
 }
@@ -172,7 +167,6 @@ SingleProcessData::SingleProcessData(SingleProcessData&& spd)
       tdisc_ode_sys(std::move(spd.tdisc_ode_sys)),
       mat_strg(spd.mat_strg),
       process(spd.process),
-      coupled_processes(spd.coupled_processes),
       process_output(std::move(spd.process_output))
 {
     spd.mat_strg = nullptr;
@@ -237,7 +231,6 @@ std::unique_ptr<SingleProcessData> makeSingleProcessData(
     Process& process,
     std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
     std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit,
-    std::unordered_map<std::type_index, Process const&>&& coupled_processes,
     ProcessOutput&& process_output)
 {
     using Tag = NumLib::NonlinearSolverTag;
@@ -249,7 +242,7 @@ std::unique_ptr<SingleProcessData> makeSingleProcessData(
         return std::make_unique<SingleProcessData>(
             std::move(timestepper), *nonlinear_solver_picard,
             std::move(conv_crit), std::move(time_disc), process,
-            std::move(coupled_processes), std::move(process_output));
+            std::move(process_output));
     }
     if (auto* nonlinear_solver_newton =
             dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
@@ -258,7 +251,7 @@ std::unique_ptr<SingleProcessData> makeSingleProcessData(
         return std::make_unique<SingleProcessData>(
             std::move(timestepper), *nonlinear_solver_newton,
             std::move(conv_crit), std::move(time_disc), process,
-            std::move(coupled_processes), std::move(process_output));
+            std::move(process_output));
     }
 
     OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
@@ -300,39 +293,12 @@ std::vector<std::unique_ptr<SingleProcessData>> createPerProcessData(
             //! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion}
             pcs_config.getConfigSubtree("convergence_criterion"));
 
-        auto const& coupled_process_tree
-            //! \ogs_file_param{prj__time_loop__processes__process__coupled_processes}
-            = pcs_config.getConfigSubtreeOptional("coupled_processes");
-        std::unordered_map<std::type_index, Process const&> coupled_processes;
-        if (coupled_process_tree)
-        {
-            for (
-                auto const cpl_pcs_name :
-                //! \ogs_file_param{prj__time_loop__processes__process__coupled_processes__coupled_process}
-                coupled_process_tree->getConfigParameterList<std::string>(
-                    "coupled_process"))
-            {
-                auto const& coupled_process = *BaseLib::getOrError(
-                    processes, cpl_pcs_name,
-                    "A process with the given name has not been defined.");
-
-                auto const inserted = coupled_processes.emplace(
-                    std::type_index(typeid(coupled_process)), coupled_process);
-                if (!inserted.second)
-                {  // insertion failed, i.e., key already exists
-                    OGS_FATAL("Coupled process `%s' already exists.",
-                              cpl_pcs_name.data());
-                }
-            }
-        }
-
         //! \ogs_file_param{prj__time_loop__processes__process__output}
         ProcessOutput process_output{pcs_config.getConfigSubtree("output")};
 
         per_process_data.emplace_back(makeSingleProcessData(
             std::move(timestepper), nl_slv, pcs, std::move(time_disc),
-            std::move(conv_crit), std::move(coupled_processes),
-            std::move(process_output)));
+            std::move(conv_crit), std::move(process_output)));
     }
 
     if (per_process_data.size() != processes.size())
@@ -504,56 +470,27 @@ UncoupledProcessesTimeLoop::UncoupledProcessesTimeLoop(
 
 bool UncoupledProcessesTimeLoop::setCoupledSolutions()
 {
-    // Do nothing if process are not coupled
-    if ((!_global_coupling_conv_crit) || _global_coupling_max_iterations == 1)
+    // All _per_process_data share one process
+    const bool use_monolithic_scheme =
+        _per_process_data[0]->process.useMonolithicScheme();
+    if (use_monolithic_scheme)
         return false;
 
-    unsigned pcs_idx = 0;
     _solutions_of_coupled_processes.reserve(_per_process_data.size());
-    for (auto& spd : _per_process_data)
+    for (unsigned pcs_idx = 0; pcs_idx < _per_process_data.size(); pcs_idx++)
     {
-        auto const& coupled_processes = spd->coupled_processes;
-        std::unordered_map<std::type_index, GlobalVector const&> coupled_xs;
-        for (auto const& coupled_process_pair : coupled_processes)
-        {
-            ProcessLib::Process const& coupled_process =
-                coupled_process_pair.second;
-            auto const found_item = std::find_if(
-                _per_process_data.begin(),
-                _per_process_data.end(),
-                [&coupled_process](
-                    std::unique_ptr<SingleProcessData> const& item) {
-                    auto const& item_process = item->process;
-                    return std::type_index(typeid(coupled_process)) ==
-                           std::type_index(typeid(item_process));
-                });
-
-            if (found_item != _per_process_data.end())
-            {
-                // Id of the coupled process:
-                const std::size_t c_id =
-                    std::distance(_per_process_data.begin(), found_item);
-
-                BaseLib::insertIfTypeIndexKeyUniqueElseError(
-                    coupled_xs, coupled_process_pair.first,
-                    *_process_solutions[c_id], "global_coupled_x");
-            }
-        }
-        _solutions_of_coupled_processes.emplace_back(coupled_xs);
-
         auto const& x = *_process_solutions[pcs_idx];
+        _solutions_of_coupled_processes.emplace_back(x);
 
         // Create a vector to store the solution of the last coupling iteration
-        auto& x_coupling0 = NumLib::GlobalVectorProvider::provider.getVector(x);
-        MathLib::LinAlg::copy(x, x_coupling0);
+        auto& x0 = NumLib::GlobalVectorProvider::provider.getVector(x);
+        MathLib::LinAlg::copy(x, x0);
 
         // append a solution vector of suitable size
-        _solutions_of_last_cpl_iteration.emplace_back(&x_coupling0);
-
-        ++pcs_idx;
-    }  // end of for (auto& spd : _per_process_data)
+        _solutions_of_last_cpl_iteration.emplace_back(&x0);
+    }
 
-    return true;
+    return true;  // use staggered scheme.
 }
 
 double UncoupledProcessesTimeLoop::computeTimeStepping(
@@ -933,8 +870,7 @@ bool UncoupledProcessesTimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
             }
 
             CoupledSolutionsForStaggeredScheme coupled_solutions(
-                spd->coupled_processes,
-                _solutions_of_coupled_processes[pcs_idx], dt);
+                _solutions_of_coupled_processes, dt, pcs_idx);
 
             spd->process.setCoupledSolutionsForStaggeredScheme(
                 &coupled_solutions);
@@ -1051,14 +987,11 @@ void UncoupledProcessesTimeLoop::outputSolutions(
                             spd->timestepper->getTimeStep().dt(), pcs_idx);
         if (is_staggered_coupling)
         {
-            CoupledSolutionsForStaggeredScheme coupled_solutions(
-                spd->coupled_processes,
-                _solutions_of_coupled_processes[pcs_idx], 0.0);
+            CoupledSolutionsForStaggeredScheme coupled_xs(
+                _solutions_of_coupled_processes, 0.0, pcs_idx);
 
-            spd->process.setCoupledSolutionsForStaggeredScheme(
-                &coupled_solutions);
-            spd->process
-                .setCoupledSolutionsForStaggeredSchemeToLocalAssemblers();
+            spd->process.setCoupledSolutionsForStaggeredScheme(&coupled_xs);
+            spd->process.setStaggeredCouplingTermToLocalAssemblers();
             (output_object.*output_class_member)(pcs, spd->process_output,
                                                  timestep, t, x);
         }
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.h b/ProcessLib/UncoupledProcessesTimeLoop.h
index 6f1d7e03d7d..0125ba1d658 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.h
+++ b/ProcessLib/UncoupledProcessesTimeLoop.h
@@ -10,9 +10,7 @@
 #pragma once
 
 #include <memory>
-#include <unordered_map>
-#include <typeinfo>
-#include <typeindex>
+#include <functional>
 
 #include <logog/include/logog.hpp>
 
@@ -77,11 +75,11 @@ private:
     std::unique_ptr<NumLib::ConvergenceCriterion> _global_coupling_conv_crit;
 
     /**
-     *  Vector of solutions of coupled processes of processes.
+     *  Vector of solutions of the coupled processes.
      *  Each vector element stores the references of the solution vectors
      *  (stored in _process_solutions) of the coupled processes of a process.
      */
-    std::vector<std::unordered_map<std::type_index, GlobalVector const&>>
+    std::vector<std::reference_wrapper<GlobalVector const>>
         _solutions_of_coupled_processes;
 
     /// Solutions of the previous coupling iteration for the convergence
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 690c9bcc23c..74646274f0f 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -15,40 +15,11 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "LocalAssemblerInterface.h"
 
+#include "CoupledSolutionsForStaggeredScheme.h"
 #include "Process.h"
 
 namespace ProcessLib
 {
-static std::unordered_map<std::type_index, const std::vector<double>>
-getPreviousLocalSolutionsOfCoupledProcesses(
-    const CoupledSolutionsForStaggeredScheme& coupled_solutions,
-    const std::vector<GlobalIndexType>& indices)
-{
-    std::unordered_map<std::type_index, const std::vector<double>>
-        local_coupled_xs0;
-
-    for (auto const& coupled_process_pair : coupled_solutions.coupled_processes)
-    {
-        auto const& coupled_pcs = coupled_process_pair.second;
-        auto const prevous_time_x = coupled_pcs.getPreviousTimeStepSolution();
-        if (prevous_time_x)
-        {
-            auto const local_coupled_x0 = prevous_time_x->get(indices);
-            BaseLib::insertIfTypeIndexKeyUniqueElseError(
-                local_coupled_xs0, coupled_process_pair.first, local_coupled_x0,
-                "local_coupled_x0");
-        }
-        else
-        {
-            const std::vector<double> local_coupled_x0;
-            BaseLib::insertIfTypeIndexKeyUniqueElseError(
-                local_coupled_xs0, coupled_process_pair.first, local_coupled_x0,
-                "local_coupled_x0");
-        }
-    }
-    return local_coupled_xs0;
-}
-
 VectorMatrixAssembler::VectorMatrixAssembler(
     std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler)
     : _jacobian_assembler(std::move(jacobian_assembler))
@@ -70,42 +41,32 @@ void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
     const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
     const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-    const CoupledSolutionsForStaggeredScheme* coupled_solutions)
+    const CoupledSolutionsForStaggeredScheme* cpl_xs)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
-    auto const local_x = x.get(indices);
 
     _local_M_data.clear();
     _local_K_data.clear();
     _local_b_data.clear();
 
-    if (!coupled_solutions)
+    if (!cpl_xs)
     {
+        auto const local_x = x.get(indices);
         local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
                                  _local_b_data);
     }
     else
     {
-        auto local_coupled_xs0 = getPreviousLocalSolutionsOfCoupledProcesses(
-            *coupled_solutions, indices);
-        auto local_coupled_xs = getCurrentLocalSolutionsOfCoupledProcesses(
-            coupled_solutions->coupled_xs, indices);
-
-        if (local_coupled_xs0.empty() || local_coupled_xs.empty())
-        {
-            local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
-                                     _local_b_data);
-        }
-        else
-        {
-            ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-                coupled_solutions->dt, coupled_solutions->coupled_processes,
-                std::move(local_coupled_xs0), std::move(local_coupled_xs));
-
-            local_assembler.assembleWithCoupledTerm(
-                t, local_x, _local_M_data, _local_K_data, _local_b_data,
-                local_coupled_solutions);
-        }
+        auto local_coupled_xs0 = getPreviousLocalSolutions(*cpl_xs, indices);
+        auto local_coupled_xs = getCurrentLocalSolutions(*cpl_xs, indices);
+
+        ProcessLib::LocalCoupledSolutions local_coupled_solutions(
+            cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
+            std::move(local_coupled_xs));
+
+        local_assembler.assembleWithCoupledTerm(t, _local_M_data, _local_K_data,
+                                                _local_b_data,
+                                                local_coupled_solutions);
     }
 
     auto const num_r_c = indices.size();
@@ -134,10 +95,9 @@ void VectorMatrixAssembler::assembleWithJacobian(
     NumLib::LocalToGlobalIndexMap const& dof_table, const double t,
     GlobalVector const& x, GlobalVector const& xdot, const double dxdot_dx,
     const double dx_dx, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
-    GlobalMatrix& Jac, const CoupledSolutionsForStaggeredScheme* coupled_solutions)
+    GlobalMatrix& Jac, const CoupledSolutionsForStaggeredScheme* cpl_xs)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
-    auto const local_x = x.get(indices);
     auto const local_xdot = xdot.get(indices);
 
     _local_M_data.clear();
@@ -145,35 +105,26 @@ void VectorMatrixAssembler::assembleWithJacobian(
     _local_b_data.clear();
     _local_Jac_data.clear();
 
-    if (!coupled_solutions)
+    if (!cpl_xs)
     {
+        auto const local_x = x.get(indices);
         _jacobian_assembler->assembleWithJacobian(
             local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
             _local_M_data, _local_K_data, _local_b_data, _local_Jac_data);
     }
     else
     {
-        auto local_coupled_xs0 = getPreviousLocalSolutionsOfCoupledProcesses(
-            *coupled_solutions, indices);
-        auto local_coupled_xs = getCurrentLocalSolutionsOfCoupledProcesses(
-            coupled_solutions->coupled_xs, indices);
-        if (local_coupled_xs0.empty() || local_coupled_xs.empty())
-        {
-            _jacobian_assembler->assembleWithJacobian(
-                local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
-                _local_M_data, _local_K_data, _local_b_data, _local_Jac_data);
-        }
-        else
-        {
-            ProcessLib::LocalCoupledSolutions local_coupled_solutions(
-                coupled_solutions->dt, coupled_solutions->coupled_processes,
-                std::move(local_coupled_xs0), std::move(local_coupled_xs));
-
-            _jacobian_assembler->assembleWithJacobianAndCoupling(
-                local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
-                _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
-                local_coupled_solutions);
-        }
+        auto local_coupled_xs0 = getPreviousLocalSolutions(*cpl_xs, indices);
+        auto local_coupled_xs = getCurrentLocalSolutions(*cpl_xs, indices);
+
+        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();
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index a1203361b4d..38b7b855a7f 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -21,6 +21,8 @@ class LocalToGlobalIndexMap;
 
 namespace ProcessLib
 {
+struct CoupledSolutionsForStaggeredScheme;
+
 class LocalAssemblerInterface;
 
 //! Utility class used to assemble global matrices and vectors.
@@ -45,7 +47,7 @@ public:
                   NumLib::LocalToGlobalIndexMap const& dof_table,
                   double const t, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b,
-                  const CoupledSolutionsForStaggeredScheme* coupled_solutions);
+                  const CoupledSolutionsForStaggeredScheme* cpl_xs);
 
     //! Assembles \c M, \c K, \c b, and the Jacobian \c Jac of the residual.
     //! \note The Jacobian must be assembled.
@@ -57,7 +59,7 @@ public:
                               const double dx_dx, GlobalMatrix& M,
                               GlobalMatrix& K, GlobalVector& b,
                               GlobalMatrix& Jac,
-                              const CoupledSolutionsForStaggeredScheme* coupled_solutions);
+                              const CoupledSolutionsForStaggeredScheme* cpl_xs);
 
 private:
     // temporary data only stored here in order to avoid frequent memory
-- 
GitLab