[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.

Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/i_coupled_processes.md

-Defines the coupled processes to a process.
\ No newline at end of file

Documentation/ProjectFile/prj/time_loop/processes/process/coupled_processes/t_coupled_process.md

-Specifies a single coupled process to a process by name.
\ No newline at end of file

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*/,

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,

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

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

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
 

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

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

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"

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,

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,

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>::

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;

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

ProcessLib/LiquidFlow/LiquidFlowProcess.h

@@ -88,7 +88,7 @@ public:
         return _material_properties.get();
     }
 
-    void setCoupledSolutionsForStaggeredSchemeToLocalAssemblers() override;
+    void setStaggeredCouplingTermToLocalAssemblers() override;
 
 private:
     void initializeConcreteProcess(

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);
     }
 }
 

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*/)
     {
     }

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

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.