From 5597890d95a1b1e9569b522843b66d8df7f09b64 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Mon, 12 Aug 2024 15:02:47 +0200
Subject: [PATCH] [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 | 114 ++++++++++++------
ProcessLib/AssemblyMixin.h | 27 +++--
ProcessLib/CreateTimeLoop.cpp | 3 +-
ProcessLib/Process.h | 7 ++
ProcessLib/SubmeshAssemblySupport.h | 6 +-
ProcessLib/TH2M/TH2MProcess.cpp | 10 +-
ProcessLib/TH2M/TH2MProcess.h | 2 +-
.../ThermoRichardsMechanicsProcess.cpp | 9 +-
.../ThermoRichardsMechanicsProcess.h | 2 +-
9 files changed, 119 insertions(+), 61 deletions(-)
diff --git a/ProcessLib/AssemblyMixin.cpp b/ProcessLib/AssemblyMixin.cpp
index 2b955a23644..bca2c4a6859 100644
--- a/ProcessLib/AssemblyMixin.cpp
+++ b/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>());
}
}
diff --git a/ProcessLib/AssemblyMixin.h b/ProcessLib/AssemblyMixin.h
index c68fb644e36..6b4bc09511f 100644
--- a/ProcessLib/AssemblyMixin.h
+++ b/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:
diff --git a/ProcessLib/CreateTimeLoop.cpp b/ProcessLib/CreateTimeLoop.cpp
index 0b22415787a..49b5bc64e85 100644
--- a/ProcessLib/CreateTimeLoop.cpp
+++ b/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);
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index b8c4fd9ed64..954a71821f0 100644
--- a/ProcessLib/Process.h
+++ b/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
{
diff --git a/ProcessLib/SubmeshAssemblySupport.h b/ProcessLib/SubmeshAssemblySupport.h
index 669b3029c34..e188674ddf6 100644
--- a/ProcessLib/SubmeshAssemblySupport.h
+++ b/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(
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index c7f66f32522..3f7e384cade 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/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;
}
diff --git a/ProcessLib/TH2M/TH2MProcess.h b/ProcessLib/TH2M/TH2MProcess.h
index 6b8942efd44..96ff5ddb1a2 100644
--- a/ProcessLib/TH2M/TH2MProcess.h
+++ b/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;
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
index 3813235dd1b..c763ac46da4 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
+++ b/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;
}
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h
index a13c2b0d338..6ae284ce890 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.h
+++ b/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;
--
GitLab