From ea24e9c7ba38ce773dce11d96b275a7db14da7c4 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 19 Jan 2017 14:01:43 +0100
Subject: [PATCH] [Coupling] Added functions to fetch the coupled solutions of
 the previous time step

---
 ProcessLib/LocalAssemblerInterface.h      | 11 +--
 ProcessLib/Process.h                      |  6 ++
 ProcessLib/StaggeredCouplingTerm.cpp      |  2 +-
 ProcessLib/StaggeredCouplingTerm.h        | 16 +++-
 ProcessLib/UncoupledProcessesTimeLoop.cpp |  2 +-
 ProcessLib/VectorMatrixAssembler.cpp      | 90 ++++++++++++++++-------
 ProcessLib/VectorMatrixAssembler.h        |  2 +-
 7 files changed, 89 insertions(+), 40 deletions(-)

diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index aeb7329343d..3cb9d7b64da 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -36,11 +36,12 @@ public:
         std::vector<double>& local_M_data, std::vector<double>& local_K_data,
         std::vector<double>& local_b_data) = 0;
 
-    virtual void coupling_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,
-        LocalCouplingTerm const& coupled_term);
+    virtual void coupling_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,
+                                   LocalCouplingTerm const& coupled_term);
 
     virtual void assembleWithJacobian(double const t,
                                       std::vector<double> const& local_x,
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 1f55b9e9b64..73ab1bc0395 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -111,6 +111,12 @@ public:
         return _secondary_variables;
     }
 
+    // Get the solution of the previous time step.
+    virtual GlobalVector* getPreviousTimeStepSolution() const
+    {
+       return nullptr;
+    }
+
     // Used as a call back for CalculateSurfaceFlux process.
     virtual std::vector<double> getFlux(std::size_t /*element_id*/,
                                         MathLib::Point3d const& /*p*/,
diff --git a/ProcessLib/StaggeredCouplingTerm.cpp b/ProcessLib/StaggeredCouplingTerm.cpp
index 5cf8ed0aa77..76c77a949a9 100644
--- a/ProcessLib/StaggeredCouplingTerm.cpp
+++ b/ProcessLib/StaggeredCouplingTerm.cpp
@@ -20,7 +20,7 @@ const StaggeredCouplingTerm createVoidStaggeredCouplingTerm()
     std::map<ProcessType, Process const&> coupled_pcsss;
     std::map<ProcessType, GlobalVector const*> coupled_xs;
     const bool empty = true;
-    return StaggeredCouplingTerm(coupled_pcsss, coupled_xs, empty);
+    return StaggeredCouplingTerm(coupled_pcsss, coupled_xs, 0.0, empty);
 }
 
 } // end of ProcessLib
diff --git a/ProcessLib/StaggeredCouplingTerm.h b/ProcessLib/StaggeredCouplingTerm.h
index c488fa38eb7..24bf9eceb66 100644
--- a/ProcessLib/StaggeredCouplingTerm.h
+++ b/ProcessLib/StaggeredCouplingTerm.h
@@ -26,28 +26,36 @@ struct StaggeredCouplingTerm
     StaggeredCouplingTerm(
         std::map<ProcessType, Process const&> const& coupled_processes_,
         std::map<ProcessType, GlobalVector const*> const& coupled_xs_,
-        const bool empty_ = false)
+        const double dt_, const bool empty_ = false)
     : coupled_processes(coupled_processes_), coupled_xs(coupled_xs_),
-      empty(empty_)
+      dt(dt_), empty(empty_)
     {
     }
 
     std::map<ProcessType, Process const&> const& coupled_processes;
     std::map<ProcessType, GlobalVector const*> const& coupled_xs;
+    const double dt; ///< Time step size.
     const bool empty;
 };
 
 struct LocalCouplingTerm
 {
-    LocalCouplingTerm(
+    LocalCouplingTerm(const double dt_,
         std::map<ProcessType, Process const&> const& coupled_processes_,
+        std::map<ProcessType, const std::vector<double>>&& local_coupled_xs0_,
         std::map<ProcessType, const std::vector<double>>&& local_coupled_xs_)
-    : coupled_processes(coupled_processes_),
+    : dt(dt_), coupled_processes(coupled_processes_),
+      local_coupled_xs0(std::move(local_coupled_xs0_)),
       local_coupled_xs(std::move(local_coupled_xs_))
     {
     }
 
+    const double dt; ///< Time step size.
+
     std::map<ProcessType, Process const&> const& coupled_processes;
+    /// Local solutions of the previous time step.
+    std::map<ProcessType, const std::vector<double>> const local_coupled_xs0;
+    /// Local solutions of the current time step.
     std::map<ProcessType, const std::vector<double>> const local_coupled_xs;
 };
 
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index c7489c21488..a6aa357a71d 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -626,7 +626,7 @@ bool UncoupledProcessesTimeLoop::loop()
                         _global_coupling_conv_crit->setNoFirstIteration();
                     StaggeredCouplingTerm coupled_term(
                         spd->coupled_processes,
-                        _solutions_of_coupled_processes[pcs_idx]);
+                        _solutions_of_coupled_processes[pcs_idx], delta_t);
 
                     nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
                         x, timestep, t, delta_t, *spd, coupled_term, *_output);
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 73f6681f00b..08fb297c0c8 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -19,6 +19,56 @@
 
 namespace ProcessLib
 {
+static std::map<ProcessLib::ProcessType, const std::vector<double>>
+getPreviousLocalSolutionsOfCoupledProcesses(
+    const StaggeredCouplingTerm& coupled_term,
+    const std::vector<GlobalIndexType>& indices)
+{
+    std::map<ProcessLib::ProcessType, const std::vector<double>>
+        local_coupled_xs0;
+    auto it = coupled_term.coupled_processes.begin();
+    while (it != coupled_term.coupled_processes.end())
+    {
+        auto const& coupled_pcs = it->second;
+        auto const prevous_time_x = coupled_pcs.getPreviousTimeStepSolution();
+        if (prevous_time_x)
+        {
+            auto const local_coupled_x0 = prevous_time_x->get(indices);
+            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs0, it->first,
+                                                   local_coupled_x0,
+                                                   "local_coupled_x0");
+        }
+        else
+        {
+            const std::vector<double> local_coupled_x0;
+            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs0, it->first,
+                                                   local_coupled_x0,
+                                                   "local_coupled_x0");
+        }
+        it++;
+    }
+    return local_coupled_xs0;
+}
+
+static std::map<ProcessLib::ProcessType, const std::vector<double>>
+getCurrentLocalSolutionsOfCoupledProcesses(
+    const std::map<ProcessType, GlobalVector const*>& global_coupled_xs,
+    const std::vector<GlobalIndexType>& indices)
+{
+    std::map<ProcessLib::ProcessType, const std::vector<double>>
+        local_coupled_xs;
+    auto it = global_coupled_xs.begin();
+    while (it != global_coupled_xs.end())
+    {
+        GlobalVector const* coupled_x = it->second;
+        auto const local_coupled_x = coupled_x->get(indices);
+        BaseLib::insertMapIfKeyUniqueElseError(
+            local_coupled_xs, it->first, local_coupled_x, "local_coupled_x");
+        it++;
+    }
+    return local_coupled_xs;
+}
+
 VectorMatrixAssembler::VectorMatrixAssembler(
     std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler)
     : _jacobian_assembler(std::move(jacobian_assembler))
@@ -45,21 +95,13 @@ void VectorMatrixAssembler::assemble(
     }
     else
     {
-        std::map<ProcessLib::ProcessType, const std::vector<double>>
-            local_coupled_xs;
-        auto it = coupled_term.coupled_xs.begin();
-        while (it != coupled_term.coupled_xs.end())
-        {
-            GlobalVector const* coupled_x = it->second;
-            auto const local_coupled_x = coupled_x->get(indices);
-            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs, it->first,
-                                                   local_coupled_x,
-                                                   "local_coupled_x");
-            it++;
-        }
-
+        auto local_coupled_xs0 =
+            getPreviousLocalSolutionsOfCoupledProcesses(coupled_term, indices);
+        auto local_coupled_xs = getCurrentLocalSolutionsOfCoupledProcesses(
+            coupled_term.coupled_xs, indices);
         ProcessLib::LocalCouplingTerm local_coupling_term(
-            coupled_term.coupled_processes, std::move(local_coupled_xs));
+            coupled_term.dt, coupled_term.coupled_processes,
+            std::move(local_coupled_xs0), std::move(local_coupled_xs));
 
         local_assembler.coupling_assemble(t, local_x, _local_M_data,
                                           _local_K_data, _local_b_data,
@@ -111,21 +153,13 @@ void VectorMatrixAssembler::assembleWithJacobian(
     }
     else
     {
-        std::map<ProcessLib::ProcessType, const std::vector<double>>
-            local_coupled_xs;
-        auto it = coupled_term.coupled_xs.begin();
-        while (it != coupled_term.coupled_xs.end())
-        {
-            GlobalVector const* coupled_x = it->second;
-            auto const local_coupled_x = coupled_x->get(indices);
-            BaseLib::insertMapIfKeyUniqueElseError(local_coupled_xs, it->first,
-                                                   local_coupled_x,
-                                                   "local_coupled_x");
-            it++;
-        }
-
+        auto local_coupled_xs0 =
+            getPreviousLocalSolutionsOfCoupledProcesses(coupled_term, indices);
+        auto local_coupled_xs = getCurrentLocalSolutionsOfCoupledProcesses(
+            coupled_term.coupled_xs, indices);
         ProcessLib::LocalCouplingTerm local_coupling_term(
-            coupled_term.coupled_processes, std::move(local_coupled_xs));
+            coupled_term.dt, coupled_term.coupled_processes,
+            std::move(local_coupled_xs0), std::move(local_coupled_xs));
 
         _jacobian_assembler->coupling_assembleWithJacobian(
             local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index 5466d6d71c0..26046c82ab5 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -66,4 +66,4 @@ private:
     std::unique_ptr<AbstractJacobianAssembler> _jacobian_assembler;
 };
 
-}   // namespace ProcessLib
+}  // namespace ProcessLib
-- 
GitLab