[PL] Extend submesh initialization for staggered

... scheme. Passing a vector of vectors (outer for processes,
inner over the process variables) of residuum names for each
process variable.

ProcessLib/AssemblyMixin.cpp

@@ -10,6 +10,10 @@
 
 #include "AssemblyMixin.h"
 
+#include <range/v3/view/enumerate.hpp>
+#include <range/v3/view/for_each.hpp>
+#include <range/v3/view/transform.hpp>
+#include <range/v3/view/zip.hpp>
 #include <unordered_set>
 
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
@@ -18,30 +22,79 @@
 
 namespace
 {
-std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>
-createResiduumVectors(
-    MeshLib::Mesh& mesh,
-    std::vector<std::string> const& residuum_names,
-    std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
-        pvs)
-{
-    auto const num_residua = residuum_names.size();
 
-    std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>
-        residuum_vectors;
-    residuum_vectors.reserve(num_residua);
+void checkResiduumNamesVsProcessVariables(
+    std::vector<std::vector<std::string>> const& per_process_residuum_names,
+    std::vector<
+        std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>> const&
+        per_process_pvs)
+{
+    if (per_process_pvs.size() != per_process_residuum_names.size())
+    {
+        OGS_FATAL(
+            "The number of passed residuum names ({}) does not match the "
+            "number of processes ({}).",
+            per_process_residuum_names.size(), per_process_pvs.size());
+    }
 
-    for (std::size_t i = 0; i < num_residua; ++i)
+    auto check_sizes = [](std::size_t const process_id, auto const& rns_pvs)
     {
-        auto const& residuum_name = residuum_names[i];
-        auto const& pv = pvs[i].get();
+        auto const& [rns, pvs] = rns_pvs;
 
-        residuum_vectors.emplace_back(*MeshLib::getOrCreateMeshProperty<double>(
-            mesh, residuum_name, MeshLib::MeshItemType::Node,
-            pv.getNumberOfGlobalComponents()));
+        if (rns.size() != pvs.size())
+        {
+            OGS_FATAL(
+                "The number of passed residuum names ({}) does not match the "
+                "number of process variables ({}) for process {}.",
+                rns.size(), pvs.size(), process_id);
+        }
+    };
+    for (auto const& [process_id, rns_pvs] :
+         ranges::views::zip(per_process_residuum_names, per_process_pvs) |
+             ranges::views::enumerate)
+    {
+        check_sizes(process_id, rns_pvs);
     }
+}
 
-    return residuum_vectors;
+std::vector<
+    std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>>
+createResiduumVectors(
+    MeshLib::Mesh& mesh,
+    std::vector<std::vector<std::string>> const& per_process_residuum_names,
+    std::vector<
+        std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>> const&
+        per_process_pvs)
+{
+    checkResiduumNamesVsProcessVariables(per_process_residuum_names,
+                                         per_process_pvs);
+
+    auto create_mesh_property_for_residuum =
+        [&mesh](std::pair<std::string const&,
+                          ProcessLib::ProcessVariable const&> const&
+                    residuum_name_process_variable)
+        -> std::reference_wrapper<MeshLib::PropertyVector<double>>
+    {
+        auto const& [rn, pv] = residuum_name_process_variable;
+        return *MeshLib::getOrCreateMeshProperty<double>(
+            mesh, rn, MeshLib::MeshItemType::Node,
+            pv.getNumberOfGlobalComponents());
+    };
+
+    auto create_mesh_properties =
+        [&create_mesh_property_for_residuum](auto const& rns_pvs)
+    {
+        auto const& [rns, pvs] = rns_pvs;
+        return ranges::views::zip(rns, pvs) |
+               ranges::views::transform(create_mesh_property_for_residuum) |
+               ranges::to<std::vector<
+                   std::reference_wrapper<MeshLib::PropertyVector<double>>>>;
+    };
+
+    return ranges::views::zip(per_process_residuum_names, per_process_pvs) |
+           ranges::views::transform(create_mesh_properties) |
+           ranges::to<std::vector<std::vector<
+               std::reference_wrapper<MeshLib::PropertyVector<double>>>>>;
 }
 }  // namespace
 
@@ -49,7 +102,8 @@ namespace ProcessLib
 {
 SubmeshAssemblyData::SubmeshAssemblyData(
     MeshLib::Mesh const& mesh,
-    std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>&&
+    std::vector<
+        std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>>&&
         residuum_vectors)
     : bulk_element_ids{*MeshLib::bulkElementIDs(mesh)},
       bulk_node_ids{*MeshLib::bulkNodeIDs(mesh)},
@@ -58,24 +112,14 @@ SubmeshAssemblyData::SubmeshAssemblyData(
 }
 
 void AssemblyMixinBase::initializeAssemblyOnSubmeshes(
-    const int process_id,
     MeshLib::Mesh& bulk_mesh,
     std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes,
-    std::vector<std::string> const& residuum_names,
-    std::vector<std::reference_wrapper<ProcessVariable>> const& pvs)
+    std::vector<std::vector<std::string>> const& residuum_names,
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const&
+        pvs)
 {
     DBUG("AssemblyMixinBase initializeSubmeshOutput().");
 
-    auto const num_residua = residuum_names.size();
-
-    if (pvs.size() != num_residua)
-    {
-        OGS_FATAL(
-            "The number of passed residuum names does not match the number "
-            "of process variables for process id {}: {} != {}",
-            process_id, num_residua, pvs.size());
-    }
-
     submesh_assembly_data_.reserve(submeshes.size());
     for (auto& mesh_ref : submeshes)
     {
@@ -184,16 +228,18 @@ void AssemblyMixinBase::copyResiduumVectorsToBulkMesh(
 }
 
 void AssemblyMixinBase::copyResiduumVectorsToSubmesh(
+    int const process_id,
     GlobalVector const& rhs,
     NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
     SubmeshAssemblyData const& sad)
 {
-    for (std::size_t variable_id = 0; variable_id < sad.residuum_vectors.size();
+    auto const& residuum_vectors = sad.residuum_vectors[process_id];
+    for (std::size_t variable_id = 0; variable_id < residuum_vectors.size();
          ++variable_id)
     {
         transformVariableFromGlobalVector(
             rhs, variable_id, local_to_global_index_map,
-            sad.residuum_vectors[variable_id].get(), sad.bulk_node_ids,
+            residuum_vectors[variable_id].get(), sad.bulk_node_ids,
             std::negate<double>());
     }
 }

ProcessLib/AssemblyMixin.h

@@ -24,13 +24,15 @@ struct SubmeshAssemblyData
 {
     explicit SubmeshAssemblyData(
         MeshLib::Mesh const& mesh,
-        std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>&&
+        std::vector<std::vector<
+            std::reference_wrapper<MeshLib::PropertyVector<double>>>>&&
             residuum_vectors);
 
     MeshLib::PropertyVector<std::size_t> const& bulk_element_ids;
     MeshLib::PropertyVector<std::size_t> const& bulk_node_ids;
     std::vector<std::size_t> active_element_ids;
-    std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>
+    std::vector<
+        std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>>
         residuum_vectors;
 };
 
@@ -50,11 +52,11 @@ protected:
         : pvma_{jacobian_assembler} {};
 
     void initializeAssemblyOnSubmeshes(
-        const int process_id,
         MeshLib::Mesh& bulk_mesh,
         std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes,
-        std::vector<std::string> const& residuum_names,
-        std::vector<std::reference_wrapper<ProcessVariable>> const& pvs);
+        std::vector<std::vector<std::string>> const& residuum_names,
+        std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const&
+            pvs);
 
     void updateActiveElements(ProcessLib::Process const& process);
 
@@ -65,6 +67,7 @@ protected:
             residuum_vectors);
 
     static void copyResiduumVectorsToSubmesh(
+        int const process_id,
         GlobalVector const& rhs,
         NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
         SubmeshAssemblyData const& sad);
@@ -74,7 +77,8 @@ private:
 
 protected:
     std::vector<SubmeshAssemblyData> submesh_assembly_data_;
-    std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>
+    std::vector<
+        std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>>
         residuum_vectors_bulk_;
 
     /// ID of the b vector on submeshes, cf. NumLib::VectorProvider.
@@ -122,13 +126,12 @@ public:
      * and staggered coupling.
      */
     void initializeAssemblyOnSubmeshes(
-        const int process_id,
         std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes,
-        std::vector<std::string> const& residuum_names)
+        std::vector<std::vector<std::string>> const& residuum_names)
     {
         AssemblyMixinBase::initializeAssemblyOnSubmeshes(
-            process_id, derived().getMesh(), submeshes, residuum_names,
-            derived().getProcessVariables(process_id));
+            derived().getMesh(), submeshes, residuum_names,
+            derived().getProcessVariables());
     }
 
     void updateActiveElements()
@@ -182,7 +185,7 @@ public:
                 MathLib::LinAlg::axpy(b, 1.0, b_submesh);
 
                 AssemblyMixinBase::copyResiduumVectorsToSubmesh(
-                    b_submesh, *(dof_tables.front()), sad);
+                    process_id, b_submesh, *(dof_tables[process_id]), sad);
             }
 
             NumLib::GlobalVectorProvider::provider.releaseVector(b_submesh);
@@ -195,7 +198,7 @@ public:
         }
 
         AssemblyMixinBase::copyResiduumVectorsToBulkMesh(
-            b, *(dof_tables.front()), residuum_vectors_bulk_);
+            b, *(dof_tables[process_id]), residuum_vectors_bulk_[process_id]);
     }
 
 private:

ProcessLib/CreateTimeLoop.cpp

@@ -12,6 +12,7 @@
 
 #include <algorithm>
 #include <range/v3/algorithm/any_of.hpp>
+#include <range/v3/view/join.hpp>
 #include <string>
 
 #include "BaseLib/ConfigTree.h"
@@ -61,7 +62,7 @@ std::unique_ptr<TimeLoop> createTimeLoop(
             auto const& residuum_vector_names =
                 process->initializeAssemblyOnSubmeshes(smroc.meshes);
 
-            for (auto const& name : residuum_vector_names)
+            for (auto const& name : residuum_vector_names | ranges::views::join)
             {
                 smroc.output.doNotProjectFromBulkMeshToSubmeshes(
                     name, MeshLib::MeshItemType::Node);

ProcessLib/Process.h

@@ -151,6 +151,13 @@ public:
     }
 
     MeshLib::Mesh& getMesh() const { return _mesh; }
+
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const&
+    getProcessVariables() const
+    {
+        return _process_variables;
+    }
+
     std::vector<std::reference_wrapper<ProcessVariable>> const&
     getProcessVariables(const int process_id) const
     {

ProcessLib/SubmeshAssemblySupport.h

@@ -32,8 +32,10 @@ public:
     /// \attention \c meshes must be a must be a non-overlapping cover of the
     /// entire simulation domain (bulk mesh)!
     ///
-    /// \return The names of the residuum vectors that will be assembled.
-    virtual std::vector<std::string> initializeAssemblyOnSubmeshes(
+    /// \return The names of the residuum vectors that will be assembled for
+    /// each process: outer vector of size 1 for monolithic schemes and greater
+    /// for staggered schemes.
+    virtual std::vector<std::vector<std::string>> initializeAssemblyOnSubmeshes(
         std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes)
     {
         DBUG(

ProcessLib/TH2M/TH2MProcess.cpp

@@ -316,17 +316,17 @@ void TH2MProcess<DisplacementDim>::computeSecondaryVariableConcrete(
 }
 
 template <int DisplacementDim>
-std::vector<std::string>
+std::vector<std::vector<std::string>>
 TH2MProcess<DisplacementDim>::initializeAssemblyOnSubmeshes(
     std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes)
 {
     INFO("TH2M process initializeSubmeshOutput().");
-    const int process_id = 0;
-    std::vector<std::string> residuum_names{
-        "GasMassFlowRate", "LiquidMassFlowRate", "HeatFlowRate", "NodalForces"};
+    std::vector<std::vector<std::string>> residuum_names{
+        {"GasMassFlowRate", "LiquidMassFlowRate", "HeatFlowRate",
+         "NodalForces"}};
 
     AssemblyMixin<TH2MProcess<DisplacementDim>>::initializeAssemblyOnSubmeshes(
-        process_id, meshes, residuum_names);
+        meshes, residuum_names);
 
     return residuum_names;
 }

ProcessLib/TH2M/TH2MProcess.h

@@ -105,7 +105,7 @@ private:
     NumLib::LocalToGlobalIndexMap const& getDOFTable(
         const int process_id) const override;
 
-    std::vector<std::string> initializeAssemblyOnSubmeshes(
+    std::vector<std::vector<std::string>> initializeAssemblyOnSubmeshes(
         std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes)
         override;
 

ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp

@@ -239,19 +239,18 @@ void ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>::
 }
 
 template <int DisplacementDim, typename ConstitutiveTraits>
-std::vector<std::string>
+std::vector<std::vector<std::string>>
 ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>::
     initializeAssemblyOnSubmeshes(
         std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes)
 {
     INFO("TRM process initializeSubmeshOutput().");
-    const int process_id = 0;
-    std::vector<std::string> residuum_names{"HeatFlowRate", "MassFlowRate",
-                                            "NodalForces"};
+    std::vector<std::vector<std::string>> residuum_names{
+        {"HeatFlowRate", "MassFlowRate", "NodalForces"}};
 
     AssemblyMixin<
         ThermoRichardsMechanicsProcess<DisplacementDim, ConstitutiveTraits>>::
-        initializeAssemblyOnSubmeshes(process_id, meshes, residuum_names);
+        initializeAssemblyOnSubmeshes(meshes, residuum_names);
 
     return residuum_names;
 }

ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h

@@ -199,7 +199,7 @@ private:
     NumLib::LocalToGlobalIndexMap const& getDOFTable(
         const int process_id) const override;
 
-    std::vector<std::string> initializeAssemblyOnSubmeshes(
+    std::vector<std::vector<std::string>> initializeAssemblyOnSubmeshes(
         std::vector<std::reference_wrapper<MeshLib::Mesh>> const& meshes)
         override;