New coupling assembly computations

Documentation/ProjectFile/prj/processes/process/HT/solid_thermal_expansion/i_solid_thermal_expansion.md

@@ -0,0 +1,2 @@
+Parameters of solid thermal expansion in flow process.
+

Documentation/ProjectFile/prj/processes/process/HT/solid_thermal_expansion/t_thermal_expansion.md

@@ -0,0 +1 @@
+Line thermal expansion of solid.

Documentation/ProjectFile/prj/processes/process/HT/t_coupling_scheme.md

@@ -0,0 +1,2 @@
+An optional input to select the coupling scheme. The variable is 'staggered'.
+If this input is omitted, the coupling scheme is the monolithic scheme by default.

Documentation/ProjectFile/prj/processes/process/HT/thermal_dispersivity/i_thermal_dispersivity.md

@@ -0,0 +1 @@
+Optional input for thermal dispersion properties.

Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_coupling_scheme.md

@@ -0,0 +1,2 @@
+An optional input to select the coupling scheme. So far, only the monolithic
+scheme is available for HYDRO_MECHANICS, and this input can be omitted.

Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS_WITH_LIE/t_coupling_scheme.md

@@ -0,0 +1,2 @@
+An optional input to select the coupling scheme. So far, only the monolithic
+scheme is available for HYDRO_MECHANICS_WITH_LIE, and this input can be omitted.

Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_coupling_scheme.md

@@ -0,0 +1,2 @@
+An optional input to select the coupling scheme. So far, only the monolithic
+scheme is available for PHASE_FIELD, and this input can be omitted.

Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_coupling_scheme.md

@@ -0,0 +1,3 @@
+An optional input to select the coupling scheme. So far, only the monolithic
+scheme is available for RichardsComponentTransport, and this input can be
+omitted.

Documentation/ProjectFile/prj/processes/process/THERMO_MECHANICS/t_coupling_scheme.md

@@ -0,0 +1,2 @@
+An optional input to select the coupling scheme. So far, only the monolithic
+scheme is available for THERMO_MECHANICS, and this input can be omitted.

ProcessLib/AbstractJacobianAssembler.h

@@ -34,7 +34,6 @@ class AbstractJacobianAssembler
     //! \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

@@ -22,13 +22,16 @@ ComponentTransportProcess::ComponentTransportProcess(
     std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     unsigned const integration_order,
-    std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+        process_variables,
     ComponentTransportProcessData&& process_data,
     SecondaryVariableCollection&& secondary_variables,
-    NumLib::NamedFunctionCaller&& named_function_caller)
+    NumLib::NamedFunctionCaller&& named_function_caller,
+    bool const use_monolithic_scheme)
     : Process(mesh, std::move(jacobian_assembler), parameters,
               integration_order, std::move(process_variables),
-              std::move(secondary_variables), std::move(named_function_caller)),
+              std::move(secondary_variables), std::move(named_function_caller),
+              use_monolithic_scheme),
       _process_data(std::move(process_data))
 {
 }
@@ -59,6 +62,8 @@ void ComponentTransportProcess::assembleConcreteProcess(
     GlobalVector& b)
 {
     DBUG("Assemble ComponentTransportProcess.");
+    if (!_use_monolithic_scheme)
+        setCoupledSolutionsOfPreviousTimeStep();
 
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
@@ -73,6 +78,9 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
 {
     DBUG("AssembleWithJacobian ComponentTransportProcess.");
 
+    if (!_use_monolithic_scheme)
+        setCoupledSolutionsOfPreviousTimeStep();
+
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberDereferenced(
         _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,

ProcessLib/ComponentTransport/ComponentTransportProcess.h

@@ -83,11 +83,12 @@ class ComponentTransportProcess final : public Process
             jacobian_assembler,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters,
         unsigned const integration_order,
-        std::vector<std::reference_wrapper<ProcessVariable>>&&
+        std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
             process_variables,
         ComponentTransportProcessData&& process_data,
         SecondaryVariableCollection&& secondary_variables,
-        NumLib::NamedFunctionCaller&& named_function_caller);
+        NumLib::NamedFunctionCaller&& named_function_caller,
+        bool const use_monolithic_scheme);
 
     //! \name ODESystem interface
     //! @{

ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp

@@ -35,19 +35,43 @@ std::unique_ptr<Process> createComponentTransportProcess(
     config.checkConfigParameter("type", "ComponentTransport");
 
     DBUG("Create ComponentTransportProcess.");
+    auto const staggered_scheme =
+        //! \ogs_file_param{prj__processes__process__coupling_scheme}
+        config.getConfigParameterOptional<std::string>("coupling_scheme");
+    const bool use_monolithic_scheme =
+        (staggered_scheme && (*staggered_scheme == "staggered")) ? false : true;
 
     // Process variable.
 
     //! \ogs_file_param{prj__processes__process__ComponentTransport__process_variables}
     auto const pv_config = config.getConfigSubtree("process_variables");
 
-    auto process_variables = findProcessVariables(
-        variables, pv_config,
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+        process_variables;
+    if (use_monolithic_scheme)  // monolithic scheme.
+    {
+        auto per_process_variables = findProcessVariables(
+            variables, pv_config,
+            {
+            //! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__concentration}
+             "concentration",
+             //! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__pressure}
+             "pressure"});
+        process_variables.push_back(std::move(per_process_variables));
+    }
+    else  // staggered scheme.
+    {
+        std::array<std::string, 2> variable_names = {
+            {"concentration",
+             "pressure"}};  // double-braces required in C++11 (not in C++14)
+
+        for (int i = 0; i < 2; i++)
         {
-        //! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__concentration}
-        "concentration",
-        //! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__pressure}
-        "pressure"});
+            auto per_process_variables =
+                findProcessVariables(variables, pv_config, {variable_names[i]});
+            process_variables.push_back(std::move(per_process_variables));
+        }
+    }
 
     MaterialLib::PorousMedium::PorousMediaProperties porous_media_properties{
         MaterialLib::PorousMedium::createPorousMediaProperties(
@@ -144,7 +168,8 @@ std::unique_ptr<Process> createComponentTransportProcess(
     return std::make_unique<ComponentTransportProcess>(
         mesh, std::move(jacobian_assembler), parameters, integration_order,
         std::move(process_variables), std::move(process_data),
-        std::move(secondary_variables), std::move(named_function_caller));
+        std::move(secondary_variables), std::move(named_function_caller),
+        use_monolithic_scheme);
 }
 
 }  // namespace ComponentTransport

ProcessLib/CoupledSolutionsForStaggeredScheme.cpp

@@ -18,48 +18,44 @@
 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 (auto const& x_t0 : cpl_xs.coupled_xs_t0)
     {
-        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);
-        }
+        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 (auto const& x_t1 : cpl_xs.coupled_xs)
     {
-        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");
+        local_xs_t1.emplace_back(x_t1.get().get(indices));
     }
-    return local_coupled_xs;
+    return local_xs_t1;
 }
 
-}  // end of ProcessLib
+}  // namespace 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/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp

@@ -42,10 +42,13 @@ std::unique_ptr<Process> createGroundwaterFlowProcess(
     //! \ogs_file_param{prj__processes__process__GROUNDWATER_FLOW__process_variables}
     auto const pv_config = config.getConfigSubtree("process_variables");
 
-    auto process_variables = findProcessVariables(
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+        process_variables;
+    auto per_process_variables = findProcessVariables(
         variables, pv_config,
         {//! \ogs_file_param_special{prj__processes__process__GROUNDWATER_FLOW__process_variables__process_variable}
          "process_variable"});
+    process_variables.push_back(std::move(per_process_variables));
 
     // Hydraulic conductivity parameter.
     auto& hydraulic_conductivity = findParameter<double>(

ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp

@@ -22,7 +22,8 @@ GroundwaterFlowProcess::GroundwaterFlowProcess(
     std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     unsigned const integration_order,
-    std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+        process_variables,
     GroundwaterFlowProcessData&& process_data,
     SecondaryVariableCollection&& secondary_variables,
     NumLib::NamedFunctionCaller&& named_function_caller,

ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h

@@ -33,7 +33,7 @@ class GroundwaterFlowProcess final : public Process
             jacobian_assembler,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters,
         unsigned const integration_order,
-        std::vector<std::reference_wrapper<ProcessVariable>>&&
+        std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
             process_variables,
         GroundwaterFlowProcessData&& process_data,
         SecondaryVariableCollection&& secondary_variables,