From 43e9bbfd05c65df2308c8db47937414d8921165f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 17 Feb 2022 14:23:10 +0100
Subject: [PATCH] Add min, max and sum for ItemArray
- These function returns the min, max or the sum of all values (for
each array of each item)
- These values are computed in parallel: ghost item values are not
taken into account
---
src/mesh/ItemArrayUtils.hpp | 296 ++++++++++++++++++++++++++++++++++
tests/test_ItemArrayUtils.cpp | 270 +++++++++++++++++++++++++++++++
2 files changed, 566 insertions(+)
diff --git a/src/mesh/ItemArrayUtils.hpp b/src/mesh/ItemArrayUtils.hpp
index 6025bd0f3..ca339de9a 100644
--- a/src/mesh/ItemArrayUtils.hpp
+++ b/src/mesh/ItemArrayUtils.hpp
@@ -10,6 +10,302 @@
#include <iostream>
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+min(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
+{
+ using ItemArrayType = ItemArray<DataType, item_type, ConnectivityPtr>;
+ using ItemIsOwnedType = ItemValue<const bool, item_type>;
+ using data_type = std::remove_const_t<typename ItemArrayType::data_type>;
+ using index_type = typename ItemArrayType::index_type;
+
+ static_assert(std::is_arithmetic_v<data_type>, "min cannot be called on non-arithmetic data");
+ static_assert(not std::is_same_v<data_type, bool>, "min cannot be called on boolean data");
+
+ class ItemArrayMin
+ {
+ private:
+ const ItemArrayType& m_item_value;
+ const ItemIsOwnedType m_is_owned;
+
+ public:
+ PUGS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_item_value.numberOfItems(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ PUGS_INLINE
+ void
+ operator()(const index_type& i_item, data_type& data) const
+ {
+ if (m_is_owned[i_item]) {
+ Array array = m_item_value[i_item];
+ for (size_t i = 0; i < m_item_value.sizeOfArrays(); ++i) {
+ if (array[i] < data) {
+ data = array[i];
+ }
+ }
+ }
+ }
+
+ PUGS_INLINE
+ void
+ join(volatile data_type& dst, const volatile data_type& src) const
+ {
+ if (src < dst) {
+ // cannot be reached if initial value is the min
+ dst = src; // LCOV_EXCL_LINE
+ }
+ }
+
+ PUGS_INLINE
+ void
+ init(data_type& value) const
+ {
+ value = std::numeric_limits<data_type>::max();
+ }
+
+ PUGS_INLINE
+ ItemArrayMin(const ItemArrayType& item_value)
+ : m_item_value(item_value), m_is_owned([&](const IConnectivity& connectivity) {
+ Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+ "unexpected connectivity dimension");
+
+ switch (connectivity.dimension()) {
+ case 1: {
+ const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+ return connectivity_1d.isOwned<item_type>();
+ break;
+ }
+ case 2: {
+ const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+ return connectivity_2d.isOwned<item_type>();
+ break;
+ }
+ case 3: {
+ const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+ return connectivity_3d.isOwned<item_type>();
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("unexpected dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }(*item_value.connectivity_ptr()))
+ {
+ ;
+ }
+
+ PUGS_INLINE
+ ~ItemArrayMin() = default;
+ };
+
+ const DataType local_min = ItemArrayMin{item_value};
+ return parallel::allReduceMin(local_min);
+}
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+max(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
+{
+ using ItemArrayType = ItemArray<DataType, item_type, ConnectivityPtr>;
+ using ItemIsOwnedType = ItemValue<const bool, item_type>;
+ using data_type = std::remove_const_t<typename ItemArrayType::data_type>;
+ using index_type = typename ItemArrayType::index_type;
+
+ static_assert(std::is_arithmetic_v<data_type>, "max cannot be called on non-arithmetic data");
+ static_assert(not std::is_same_v<data_type, bool>, "max cannot be called on boolean data");
+
+ class ItemArrayMax
+ {
+ private:
+ const ItemArrayType& m_item_value;
+ const ItemIsOwnedType m_is_owned;
+
+ public:
+ PUGS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_item_value.numberOfItems(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ PUGS_INLINE
+ void
+ operator()(const index_type& i_item, data_type& data) const
+ {
+ if (m_is_owned[i_item]) {
+ Array array = m_item_value[i_item];
+ for (size_t i = 0; i < m_item_value.sizeOfArrays(); ++i) {
+ if (array[i] > data) {
+ data = array[i];
+ }
+ }
+ }
+ }
+
+ PUGS_INLINE
+ void
+ join(volatile data_type& dst, const volatile data_type& src) const
+ {
+ if (src > dst) {
+ // cannot be reached if initial value is the max
+ dst = src; // LCOV_EXCL_LINE
+ }
+ }
+
+ PUGS_INLINE
+ void
+ init(data_type& value) const
+ {
+ value = std::numeric_limits<data_type>::min();
+ }
+
+ PUGS_INLINE
+ ItemArrayMax(const ItemArrayType& item_value)
+ : m_item_value(item_value), m_is_owned([&](const IConnectivity& connectivity) {
+ Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+ "unexpected connectivity dimension");
+
+ switch (connectivity.dimension()) {
+ case 1: {
+ const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+ return connectivity_1d.isOwned<item_type>();
+ break;
+ }
+ case 2: {
+ const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+ return connectivity_2d.isOwned<item_type>();
+ break;
+ }
+ case 3: {
+ const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+ return connectivity_3d.isOwned<item_type>();
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("unexpected dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }(*item_value.connectivity_ptr()))
+ {
+ ;
+ }
+
+ PUGS_INLINE
+ ~ItemArrayMax() = default;
+ };
+
+ const DataType local_max = ItemArrayMax{item_value};
+ return parallel::allReduceMax(local_max);
+}
+
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+sum(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
+{
+ using ItemArrayType = ItemArray<DataType, item_type, ConnectivityPtr>;
+ using ItemIsOwnedType = ItemValue<const bool, item_type>;
+ using data_type = std::remove_const_t<typename ItemArrayType::data_type>;
+ using index_type = typename ItemArrayType::index_type;
+
+ static_assert(not std::is_same_v<data_type, bool>, "sum cannot be called on boolean data");
+
+ class ItemArraySum
+ {
+ private:
+ const ItemArrayType& m_item_value;
+ const ItemIsOwnedType m_is_owned;
+
+ public:
+ PUGS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_item_value.numberOfItems(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ PUGS_INLINE
+ void
+ operator()(const index_type& i_item, data_type& data) const
+ {
+ if (m_is_owned[i_item]) {
+ Array array = m_item_value[i_item];
+ for (size_t i = 0; i < m_item_value.sizeOfArrays(); ++i) {
+ data += array[i];
+ }
+ }
+ }
+
+ PUGS_INLINE
+ void
+ join(volatile data_type& dst, const volatile data_type& src) const
+ {
+ dst += src;
+ }
+
+ PUGS_INLINE
+ void
+ init(data_type& value) const
+ {
+ if constexpr (std::is_arithmetic_v<data_type>) {
+ value = 0;
+ } else {
+ static_assert(is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>, "invalid data type");
+ value = zero;
+ }
+ }
+
+ PUGS_INLINE
+ ItemArraySum(const ItemArrayType& item_value)
+ : m_item_value(item_value), m_is_owned([&](const IConnectivity& connectivity) {
+ Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+ "unexpected connectivity dimension");
+
+ switch (connectivity.dimension()) {
+ case 1: {
+ const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+ return connectivity_1d.isOwned<item_type>();
+ break;
+ }
+ case 2: {
+ const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+ return connectivity_2d.isOwned<item_type>();
+ break;
+ }
+ case 3: {
+ const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+ return connectivity_3d.isOwned<item_type>();
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("unexpected dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }(*item_value.connectivity_ptr()))
+ {
+ ;
+ }
+
+ PUGS_INLINE
+ ~ItemArraySum() = default;
+ };
+
+ const DataType local_sum = ItemArraySum{item_value};
+ return parallel::allReduceSum(local_sum);
+}
+
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void
synchronize(ItemArray<DataType, item_type, ConnectivityPtr>& item_array)
diff --git a/tests/test_ItemArrayUtils.cpp b/tests/test_ItemArrayUtils.cpp
index 85e2e6e07..f4ab1f47d 100644
--- a/tests/test_ItemArrayUtils.cpp
+++ b/tests/test_ItemArrayUtils.cpp
@@ -839,4 +839,274 @@ TEST_CASE("ItemArrayUtils", "[mesh]")
}
}
}
+
+ SECTION("min")
+ {
+ SECTION("1D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+ CellArray<int> cell_array{connectivity, 3};
+ cell_array.fill(-1);
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ mesh_1d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+ cell_array[cell_id][i] = 10 + parallel::rank() + i;
+ }
+ }
+ });
+
+ REQUIRE(min(cell_array) == 10);
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ CellArray<int> cell_array{connectivity, 3};
+ cell_array.fill(-1);
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ mesh_2d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+ cell_array[cell_id][i] = 10 + parallel::rank() - i;
+ }
+ }
+ });
+
+ REQUIRE(min(cell_array) == 8);
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+ CellArray<int> cell_array{connectivity, 3};
+ cell_array.fill(-1);
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ mesh_3d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+ cell_array[cell_id][i] = 10 + parallel::rank() + 2 - i;
+ }
+ }
+ });
+
+ REQUIRE(min(cell_array) == 10);
+ }
+ }
+ }
+ }
+
+ SECTION("max")
+ {
+ SECTION("1D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+ CellArray<size_t> cell_array{connectivity, 3};
+ cell_array.fill(std::numeric_limits<size_t>::max());
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ mesh_1d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+ cell_array[cell_id][i] = parallel::rank() + 1 + i;
+ }
+ }
+ });
+
+ REQUIRE(max(cell_array) == parallel::size() + 2);
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ CellArray<size_t> cell_array{connectivity, 3};
+ cell_array.fill(std::numeric_limits<size_t>::max());
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ mesh_2d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+ cell_array[cell_id][i] = parallel::rank() + i;
+ }
+ }
+ });
+
+ REQUIRE(max(cell_array) == parallel::size() + 1);
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+ CellArray<size_t> cell_array{connectivity, 2};
+ cell_array.fill(std::numeric_limits<size_t>::max());
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ mesh_3d->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ for (size_t i = 0; i < cell_array.sizeOfArrays(); ++i) {
+ cell_array[cell_id][i] = parallel::rank() + 1 + i;
+ }
+ }
+ });
+
+ REQUIRE(max(cell_array) == parallel::size() + 1);
+ }
+ }
+ }
+ }
+
+ SECTION("sum")
+ {
+ SECTION("1D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+ CellArray<size_t> cell_array{connectivity, 3};
+ cell_array.fill(5);
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+
+ const size_t global_number_of_cells = [&] {
+ size_t number_of_cells = 0;
+ for (CellId cell_id = 0; cell_id < cell_is_owned.numberOfItems(); ++cell_id) {
+ number_of_cells += cell_is_owned[cell_id];
+ }
+ return parallel::allReduceSum(number_of_cells);
+ }();
+
+ REQUIRE(sum(cell_array) == 5 * cell_array.sizeOfArrays() * global_number_of_cells);
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ FaceArray<size_t> face_array{connectivity, 3};
+ face_array.fill(2);
+
+ auto face_is_owned = connectivity.faceIsOwned();
+
+ const size_t global_number_of_faces = [&] {
+ size_t number_of_faces = 0;
+ for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+ number_of_faces += face_is_owned[face_id];
+ }
+ return parallel::allReduceSum(number_of_faces);
+ }();
+
+ REQUIRE(sum(face_array) == 2 * face_array.sizeOfArrays() * global_number_of_faces);
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+ NodeArray<size_t> node_array{connectivity, 3};
+ node_array.fill(3);
+
+ auto node_is_owned = connectivity.nodeIsOwned();
+
+ const size_t global_number_of_nodes = [&] {
+ size_t number_of_nodes = 0;
+ for (NodeId node_id = 0; node_id < node_is_owned.numberOfItems(); ++node_id) {
+ number_of_nodes += node_is_owned[node_id];
+ }
+ return parallel::allReduceSum(number_of_nodes);
+ }();
+
+ REQUIRE(sum(node_array) == 3 * node_array.sizeOfArrays() * global_number_of_nodes);
+ }
+ }
+ }
+ }
}
--
GitLab