From bc2e50609490d4b160636494a62fd0c3cbc1b7b7 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri, 10 Feb 2017 13:51:57 +0100
Subject: [PATCH] [Coupling] Implemented the computation of the secondary
 variables with

 coupled processes.
---
 .../GroundwaterFlowProcess.cpp                |  7 +++--
 .../GroundwaterFlow/GroundwaterFlowProcess.h  |  4 ++-
 .../HeatConduction/HeatConductionProcess.cpp  |  7 +++--
 .../HeatConduction/HeatConductionProcess.h    |  5 ++--
 .../HydroMechanics/HydroMechanicsProcess.cpp  |  5 ++--
 .../HydroMechanics/HydroMechanicsProcess.h    |  4 ++-
 .../LiquidFlowLocalAssembler-impl.h           | 29 +++++++++++++++++++
 .../LiquidFlow/LiquidFlowLocalAssembler.h     |  9 +++++-
 ProcessLib/LiquidFlow/LiquidFlowProcess.cpp   |  8 +++--
 ProcessLib/LiquidFlow/LiquidFlowProcess.h     |  6 ++--
 ProcessLib/LocalAssemblerInterface.cpp        | 17 +++++++++--
 ProcessLib/LocalAssemblerInterface.h          | 13 ++++++++-
 ProcessLib/Process.cpp                        |  6 ++--
 ProcessLib/Process.h                          |  7 +++--
 ProcessLib/UncoupledProcessesTimeLoop.cpp     | 17 +++++++----
 15 files changed, 114 insertions(+), 30 deletions(-)

diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index 68a5c72a375..d50fb4e08c6 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -104,13 +104,14 @@ void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
 }
 
 
-void GroundwaterFlowProcess::computeSecondaryVariableConcrete(const double t,
-                                                         GlobalVector const& x)
+void GroundwaterFlowProcess::computeSecondaryVariableConcrete(
+     const double t, GlobalVector const& x,
+     StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("Compute the velocity for GroundwaterFlowProcess.");
     GlobalExecutor::executeMemberOnDereferenced(
             &GroundwaterFlowLocalAssemblerInterface::computeSecondaryVariable,
-            _local_assemblers, *_local_to_global_index_map, t, x);
+            _local_assemblers, *_local_to_global_index_map, t, x, coupled_term);
 }
 
 }   // namespace GroundwaterFlow
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index 694a84aa0fd..eadd1739792 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -88,7 +88,9 @@ public:
     }
 
     void computeSecondaryVariableConcrete(double const t,
-                                          GlobalVector const& x) override;
+                                          GlobalVector const& x,
+                                          StaggeredCouplingTerm const&
+                                          coupled_term) override;
 
 private:
     void initializeConcreteProcess(
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index c32b9ec8c1c..d341c30a4f9 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -115,13 +115,14 @@ void HeatConductionProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, coupling_term);
 }
 
-void HeatConductionProcess::computeSecondaryVariableConcrete(const double t,
-                                                         GlobalVector const& x)
+void HeatConductionProcess::computeSecondaryVariableConcrete(
+    const double t, GlobalVector const& x,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("Compute heat flux for HeatConductionProcess.");
     GlobalExecutor::executeMemberOnDereferenced(
             &HeatConductionLocalAssemblerInterface::computeSecondaryVariable,
-            _local_assemblers, *_local_to_global_index_map, t, x);
+            _local_assemblers, *_local_to_global_index_map, t, x, coupled_term);
 }
 
 }  // namespace HeatConduction
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h
index 06f54405591..ca7c77a8330 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.h
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.h
@@ -38,8 +38,9 @@ public:
 
     bool isLinear() const override { return true; }
 
-    void computeSecondaryVariableConcrete(double const t,
-                                          GlobalVector const& x) override;
+    void computeSecondaryVariableConcrete(
+        double const t, GlobalVector const& x,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     void preTimestepConcreteProcess(GlobalVector const& x, const double t,
                                     const double delta_t) override;
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 209c6169aa8..9b6be35e0d9 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -349,12 +349,13 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
 
 template <unsigned GlobalDim>
 void HydroMechanicsProcess<GlobalDim>::computeSecondaryVariableConcrete(
-    const double t, GlobalVector const& x)
+    const double t, GlobalVector const& x,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("Compute the secondary variables for HydroMechanicsProcess.");
     GlobalExecutor::executeMemberOnDereferenced(
         &HydroMechanicsLocalAssemblerInterface::computeSecondaryVariable,
-        _local_assemblers, *_local_to_global_index_map, t, x);
+        _local_assemblers, *_local_to_global_index_map, t, x, coupled_term);
 }
 
 // ------------------------------------------------------------------------------------
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index 92a2f3043b8..a3ecb7442f5 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -53,7 +53,9 @@ public:
     //! @}
 
     void computeSecondaryVariableConcrete(double const t,
-                                          GlobalVector const& x) override;
+                                          GlobalVector const& x,
+                                          StaggeredCouplingTerm const&
+                                          coupled_term) override;
 
 private:
     using LocalAssemblerInterface = HydroMechanicsLocalAssemblerInterface;
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index e460bd7075c..06c45aae362 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -306,6 +306,35 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     }
 }
 
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+    computeSecondaryVariableWithCoupledProcessConcrete(
+        double const t, std::vector<double> const& local_x,
+        std::unordered_map<std::type_index, const std::vector<double>> const&
+            coupled_local_solutions)
+{
+    SpatialPosition pos;
+    pos.setElementID(_element.getID());
+    const int material_id = _material_properties.getMaterialID(pos);
+    const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+        material_id, t, pos, _element.getDimension());
+
+    const auto local_T = coupled_local_solutions.at(std::type_index(
+        typeid(ProcessLib::HeatConduction::HeatConductionProcess)));
+
+    // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
+    //  the assert must be changed to perm.rows() == _element->getDimension()
+    assert(perm.rows() == GlobalDim || perm.rows() == 1);
+
+    if (perm.size() == 1)  // isotropic or 1D problem.
+        computeSecondaryVariableCoupledWithHeatTransportLocal<
+            IsotropicCalculator>(t, local_x, local_T, pos, perm);
+    else
+        computeSecondaryVariableCoupledWithHeatTransportLocal<
+            AnisotropicCalculator>(t, local_x, local_T, pos, perm);
+}
+
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 template <typename LaplacianGravityVelocityCalculator>
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 6b2bf9e5ce3..ad3aded431f 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -12,8 +12,10 @@
 
 #pragma once
 
+#include <map>
+#include <unordered_map>
 #include <vector>
-#include <memory>
+#include <typeindex>
 
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
@@ -95,6 +97,11 @@ public:
     void computeSecondaryVariableConcrete(
         double const /*t*/, std::vector<double> const& local_x) override;
 
+    void computeSecondaryVariableWithCoupledProcessConcrete(
+        double const t, std::vector<double> const& local_x,
+        std::unordered_map<std::type_index, const std::vector<double>> const&
+            coupled_local_solutions) override;
+
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 4ab25b277b4..74c9bde7697 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -115,13 +115,15 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, coupling_term);
 }
 
-void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
-                                                         GlobalVector const& x)
+void LiquidFlowProcess::computeSecondaryVariableConcrete(
+    const double t,
+    GlobalVector const& x,
+    StaggeredCouplingTerm const& coupled_term)
 {
     DBUG("Compute the velocity for LiquidFlowProcess.");
     GlobalExecutor::executeMemberOnDereferenced(
         &LiquidFlowLocalAssemblerInterface::computeSecondaryVariable,
-        _local_assemblers, *_local_to_global_index_map, t, x);
+        _local_assemblers, *_local_to_global_index_map, t, x, coupled_term);
 }
 
 }  // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index fdb2092e790..dd5836da3ac 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -72,8 +72,10 @@ public:
         double const gravitational_acceleration,
         BaseLib::ConfigTree const& config);
 
-    void computeSecondaryVariableConcrete(double const t,
-                                          GlobalVector const& x) override;
+    void computeSecondaryVariableConcrete(
+        double const t,
+        GlobalVector const& x,
+        StaggeredCouplingTerm const& coupled_term) override;
 
     bool isLinear() const override { return true; }
     int getGravitationalAxisID() const { return _gravitational_axis_id; }
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 8eb8198e77b..336ee0f67e1 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -56,11 +56,24 @@ void LocalAssemblerInterface::assembleWithJacobianAndCouping(
 void LocalAssemblerInterface::computeSecondaryVariable(
                               std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
-                              double const t, GlobalVector const& x)
+                              double const t, GlobalVector const& x,
+                              StaggeredCouplingTerm const& coupled_term)
 {
     auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
     auto const local_x = x.get(indices);
-    computeSecondaryVariableConcrete(t, local_x);
+
+    if (coupled_term.empty)
+    {
+        computeSecondaryVariableConcrete(t, local_x);
+    }
+    else
+    {
+        auto const local_coupled_xs
+            = getCurrentLocalSolutionsOfCoupledProcesses(
+                    coupled_term.coupled_xs, indices);
+        computeSecondaryVariableWithCoupledProcessConcrete(t, local_x,
+                                                           local_coupled_xs);
+    }
 }
 
 void LocalAssemblerInterface::preTimestep(
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index a70871b4b33..426e5d87e6d 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -9,6 +9,10 @@
 
 #pragma once
 
+
+#include <unordered_map>
+#include <typeindex>
+
 #include "NumLib/NumericsConfig.h"
 #include "MathLib/Point3d.h"
 #include "StaggeredCouplingTerm.h"
@@ -64,7 +68,8 @@ public:
 
     virtual void computeSecondaryVariable(std::size_t const mesh_item_id,
                               NumLib::LocalToGlobalIndexMap const& dof_table,
-                              const double t, GlobalVector const& x);
+                              const double t, GlobalVector const& x,
+                              StaggeredCouplingTerm const& coupled_term);
 
     virtual void preTimestep(std::size_t const mesh_item_id,
                              NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -94,6 +99,12 @@ private:
 
     virtual void computeSecondaryVariableConcrete
                 (double const /*t*/, std::vector<double> const& /*local_x*/) {}
+
+    virtual void computeSecondaryVariableWithCoupledProcessConcrete
+            (double const /*t*/, std::vector<double> const& /*local_x*/,
+             std::unordered_map<std::type_index,
+             const std::vector<double>> const&
+             /*coupled_local_solutions*/) {}
 };
 
 } // namespace ProcessLib
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index addff2becea..0fb1ca21936 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -246,9 +246,11 @@ void Process::postTimestep(GlobalVector const& x)
     postTimestepConcreteProcess(x);
 }
 
-void Process::computeSecondaryVariable( const double t, GlobalVector const& x)
+void Process::computeSecondaryVariable(const double t, GlobalVector const& x,
+                                       StaggeredCouplingTerm const&
+                                       coupled_term)
 {
-    computeSecondaryVariableConcrete(t, x);
+    computeSecondaryVariableConcrete(t, x, coupled_term);
 }
 
 void Process::preIteration(const unsigned iter, const GlobalVector &x)
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 7bec88b5785..674c898b132 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -60,7 +60,8 @@ public:
                       GlobalVector const& x) override final;
 
     /// compute secondary variables for the coupled equations or for output.
-    void computeSecondaryVariable(const double t, GlobalVector const& x);
+    void computeSecondaryVariable(const double t, GlobalVector const& x,
+                                  StaggeredCouplingTerm const& coupled_term);
 
     NumLib::IterationResult postIteration(GlobalVector const& x) override final;
 
@@ -165,7 +166,9 @@ private:
                                              GlobalVector const& /*x*/){}
 
     virtual void computeSecondaryVariableConcrete(const double /*t*/,
-                                                  GlobalVector const& /*x*/) {}
+                                                  GlobalVector const& /*x*/,
+                                                  StaggeredCouplingTerm
+                                                  const& /*coupled_term*/) {}
 
     virtual NumLib::IterationResult postIterationConcreteProcess(
         GlobalVector const& /*x*/)
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index c07fc0eedb9..c7666e8ae0c 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -659,11 +659,14 @@ bool UncoupledProcessesTimeLoop::solveUncoupledEquationSystems(
         auto& x = *_process_solutions[pcs_idx];
         pcs.preTimestep(x, t, dt);
 
-        const auto nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
-            x, timestep_id, t, dt, *spd,
-            ProcessLib::createVoidStaggeredCouplingTerm(), *_output);
+        const auto void_staggered_coupling_term =
+            ProcessLib::createVoidStaggeredCouplingTerm();
+
+        const auto nonlinear_solver_succeeded =
+            solveOneTimeStepOneProcess(x, timestep_id, t, dt, *spd,
+                                       void_staggered_coupling_term, *_output);
         pcs.postTimestep(x);
-        pcs.computeSecondaryVariable(t, x);
+        pcs.computeSecondaryVariable(t, x, void_staggered_coupling_term);
 
         INFO("[time] Solving process #%u took %g s in time step #%u ", pcs_idx,
              time_timestep_process.elapsed(), timestep_id);
@@ -791,7 +794,11 @@ bool UncoupledProcessesTimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
         auto& pcs = spd->process;
         auto& x = *_process_solutions[pcs_idx];
         pcs.postTimestep(x);
-        pcs.computeSecondaryVariable(t, x);
+
+        StaggeredCouplingTerm coupled_term(
+            spd->coupled_processes, _solutions_of_coupled_processes[pcs_idx],
+            0.0);
+        pcs.computeSecondaryVariable(t, x, coupled_term);
 
         _output->doOutput(pcs, spd->process_output, timestep_id, t, x);
         ++pcs_idx;
-- 
GitLab