From 574995195094857bd160f8ee1bb5a6fb60fe8a4c Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 15 Feb 2019 18:37:21 +0100
Subject: [PATCH] Rework reductions
- define min, max and sum functions for generic array types. These are
local (i.e. non-parallel) reductions
- define min, max and sum functions for item value types. These are
parallel reductions
Parallel versions are not finished since
- ghost items need to be ignored
- mpi sum reduction needs to take into account non scalar
types (i.e. TinyVector's and TinyMatrix's for instance)
---
src/algebra/TinyMatrix.hpp | 9 ++
src/algebra/TinyVector.hpp | 9 ++
src/mesh/GmshReader.cpp | 10 +-
src/mesh/ItemValueUtils.hpp | 191 ++++++++++++++++++++++++
src/output/VTKWriter.hpp | 7 -
src/scheme/AcousticSolver.hpp | 3 +-
src/utils/ArrayUtils.hpp | 263 +++++++++++++++++++---------------
src/utils/Messenger.hpp | 31 ++++
tests/test_ArrayUtils.cpp | 107 ++++++--------
9 files changed, 438 insertions(+), 192 deletions(-)
create mode 100644 src/mesh/ItemValueUtils.hpp
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 453063b76..bafd84bdc 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -173,6 +173,15 @@ public:
return *this;
}
+ PASTIS_INLINE
+ constexpr volatile TinyMatrix& operator+=(const volatile TinyMatrix& A) volatile
+ {
+ for (size_t i=0; i<N*N; ++i) {
+ m_values[i] += A.m_values[i];
+ }
+ return *this;
+ }
+
PASTIS_INLINE
constexpr TinyMatrix& operator-=(const TinyMatrix& A)
{
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index c39bd5ccc..417d26fad 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -150,6 +150,15 @@ class TinyVector
return *this;
}
+ PASTIS_INLINE
+ constexpr volatile TinyVector& operator+=(const volatile TinyVector& v) volatile
+ {
+ for (size_t i=0; i<N; ++i) {
+ m_values[i] += v.m_values[i];
+ }
+ return *this;
+ }
+
PASTIS_INLINE
constexpr TinyVector& operator-=(const TinyVector& v)
{
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 0ff0f36ad..f3aba4e18 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -576,7 +576,7 @@ void GmshReader::_dispatch()
std::vector<Array<int>> recv_cell_node_number_by_proc(parallel::size());
for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
recv_cell_node_number_by_proc[i_rank]
- = Array<int>(Sum(recv_number_of_node_per_cell_by_proc[i_rank]));
+ = Array<int>(sum(recv_number_of_node_per_cell_by_proc[i_rank]));
}
parallel::exchange(cell_node_number_to_send_by_proc, recv_cell_node_number_by_proc);
@@ -765,7 +765,7 @@ void GmshReader::_dispatch()
std::vector<Array<int>> recv_cell_face_number_by_proc(parallel::size());
for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
recv_cell_face_number_by_proc[i_rank]
- = Array<int>(Sum(recv_number_of_face_per_cell_by_proc[i_rank]));
+ = Array<int>(sum(recv_number_of_face_per_cell_by_proc[i_rank]));
}
parallel::exchange(cell_face_number_to_send_by_proc, recv_cell_face_number_by_proc);
@@ -932,7 +932,7 @@ void GmshReader::_dispatch()
std::vector<Array<int>> recv_face_node_number_by_proc(parallel::size());
for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
recv_face_node_number_by_proc[i_rank]
- = Array<int>(Sum(recv_face_number_of_node_by_proc[i_rank]));
+ = Array<int>(sum(recv_face_number_of_node_by_proc[i_rank]));
}
parallel::exchange(face_node_number_to_send_by_proc, recv_face_node_number_by_proc);
@@ -1009,7 +1009,7 @@ void GmshReader::_dispatch()
parallel::broadcast(ref_name_size_list, sender_rank);
// sending references name size
- Array<RefId::TagNameType::value_type> ref_name_cat{Sum{ref_name_size_list}};
+ Array<RefId::TagNameType::value_type> ref_name_cat{sum(ref_name_size_list)};
if (parallel::rank() == sender_rank){
size_t i_char=0;
for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
@@ -1131,7 +1131,7 @@ void GmshReader::_dispatch()
Array<const size_t> number_of_ref_node_list_per_proc
= parallel::allGather(mesh.connectivity().numberOfRefNodeList());
- if (Max(number_of_ref_node_list_per_proc) > 0) {
+ if (max(number_of_ref_node_list_per_proc) > 0) {
perr() << __FILE__ << ':' << __LINE__ << ": "
<< rang::fgB::red << "exchange of node ref list not implemented yet!"
<< rang::fg::reset << '\n';
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
new file mode 100644
index 000000000..138137397
--- /dev/null
+++ b/src/mesh/ItemValueUtils.hpp
@@ -0,0 +1,191 @@
+#ifndef ITEM_VALUE_UTILS_HPP
+#define ITEM_VALUE_UTILS_HPP
+
+#include <Messenger.hpp>
+#include <ItemValue.hpp>
+
+template <typename DataType,
+ ItemType item_type>
+std::remove_const_t<DataType>
+min(const ItemValue<DataType, item_type>& item_value)
+{
+ using ItemValueType = ItemValue<DataType, item_type>;
+ using data_type = std::remove_const_t<typename ItemValueType::data_type>;
+ using index_type = typename ItemValueType::index_type;
+
+ class ItemValueMin
+ {
+ private:
+ const ItemValueType& m_item_value;
+
+ public:
+ PASTIS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_item_value.size(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ PASTIS_INLINE
+ void operator()(const index_type& i, data_type& data) const
+ {
+ if (m_item_value[i] < data) {
+ data = m_item_value[i];
+ }
+ }
+
+ PASTIS_INLINE
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ if (src < dst) {
+ dst = src;
+ }
+ }
+
+ PASTIS_INLINE
+ void init(data_type& value) const
+ {
+ value = std::numeric_limits<data_type>::max();
+ }
+
+ PASTIS_INLINE
+ ItemValueMin(const ItemValueType& item_value)
+ : m_item_value(item_value)
+ {
+ ;
+ }
+
+ PASTIS_INLINE
+ ~ItemValueMin() = default;
+ };
+
+ const DataType local_min = ItemValueMin{item_value};
+ return parallel::allReduceMin(local_min);
+}
+
+template <typename DataType,
+ ItemType item_type>
+std::remove_const_t<DataType>
+max(const ItemValue<DataType, item_type>& item_value)
+{
+ using ItemValueType = ItemValue<DataType, item_type>;
+ using data_type = std::remove_const_t<typename ItemValueType::data_type>;
+ using index_type = typename ItemValueType::index_type;
+
+ class ItemValueMax
+ {
+ private:
+ const ItemValueType& m_item_value;
+
+ public:
+ PASTIS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_item_value.size(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ PASTIS_INLINE
+ void operator()(const index_type& i, data_type& data) const
+ {
+ if (m_item_value[i] > data) {
+ data = m_item_value[i];
+ }
+ }
+
+ PASTIS_INLINE
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ if (src > dst) {
+ dst = src;
+ }
+ }
+
+ PASTIS_INLINE
+ void init(data_type& value) const
+ {
+ value = std::numeric_limits<data_type>::min();
+ }
+
+ PASTIS_INLINE
+ ItemValueMax(const ItemValueType& item_value)
+ : m_item_value(item_value)
+ {
+ ;
+ }
+
+ PASTIS_INLINE
+ ~ItemValueMax() = default;
+ };
+
+ const DataType local_max = ItemValueMax{item_value};
+ return parallel::allReduceMax(local_max);
+}
+
+
+template <typename DataType,
+ ItemType item_type>
+std::remove_const_t<DataType>
+sum(const ItemValue<DataType, item_type>& item_value)
+{
+ using ItemValueType = ItemValue<DataType, item_type>;
+ using data_type = std::remove_const_t<typename ItemValueType::data_type>;
+ using index_type = typename ItemValueType::index_type;
+
+ class ItemValueSum
+ {
+ private:
+ const ItemValueType& m_item_value;
+
+ public:
+ PASTIS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_item_value.size(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ PASTIS_INLINE
+ void operator()(const index_type& i, data_type& data) const
+ {
+ data += m_item_value[i];
+ }
+
+ PASTIS_INLINE
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ dst += src;
+ }
+
+ PASTIS_INLINE
+ void init(data_type& value) const
+ {
+ if constexpr (std::is_arithmetic_v<data_type>) {
+ value = 0;
+ } else {
+ value = zero;
+ }
+ }
+
+ PASTIS_INLINE
+ ItemValueSum(const ItemValueType& item_value)
+ : m_item_value(item_value)
+ {
+ ;
+ }
+
+ PASTIS_INLINE
+ ~ItemValueSum() = default;
+ };
+
+ const DataType local_sum = ItemValueSum{item_value};
+ return parallel::allReduceSum(local_sum);
+}
+
+#endif // ITEM_VALUE_UTILS_HPP
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index f26e0b4fd..dc3d1bb12 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -361,13 +361,6 @@ class VTKWriter
fout << "</VTKFile>\n";
}
m_file_number++;
-
- if (parallel::size() > 1) {
- perr() << __FILE__ << ':' << __LINE__ << ": stopping parallel execution\n";
- parallel::barrier();
- parallel::messenger().destroy();
- std::exit(0);
- }
}
VTKWriter(const std::string& base_filename,
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 3cb320957..fe5e8bf18 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -17,6 +17,7 @@
#include <SubItemValuePerItem.hpp>
#include <Messenger.hpp>
+#include <ItemValueUtils.hpp>
template<typename MeshData>
class AcousticSolver
@@ -265,7 +266,7 @@ class AcousticSolver
m_Vj_over_cj[j] = 2*Vj[j]/(S*cj[j]);
});
- return Min(m_Vj_over_cj);
+ return min(m_Vj_over_cj);
}
void computeNextStep(const double&, const double& dt,
diff --git a/src/utils/ArrayUtils.hpp b/src/utils/ArrayUtils.hpp
index 6392eddd6..51e423094 100644
--- a/src/utils/ArrayUtils.hpp
+++ b/src/utils/ArrayUtils.hpp
@@ -4,156 +4,189 @@
#include <PastisMacros.hpp>
#include <PastisUtils.hpp>
-template <typename ArrayType>
-class Min
+#include <Types.hpp>
+
+template <typename DataType,
+ template <typename> typename ArrayT>
+std::remove_const_t<DataType>
+min(const ArrayT<DataType>& array)
{
- private:
- const ArrayType& m_array;
+ using ArrayType = ArrayT<DataType>;
using data_type = std::remove_const_t<typename ArrayType::data_type>;
using index_type = typename ArrayType::index_type;
- public:
- PASTIS_INLINE
- operator data_type()
+ class ArrayMin
{
- data_type reduced_value;
- parallel_reduce(m_array.size(), *this, reduced_value);
- return reduced_value;
- }
- PASTIS_INLINE
- void operator()(const index_type& i, data_type& data) const
- {
- if (m_array[i] < data) {
- data = m_array[i];
+ private:
+ const ArrayType m_array;
+
+ public:
+ PASTIS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_array.size(), *this, reduced_value);
+ return reduced_value;
}
- }
- PASTIS_INLINE
- void join(volatile data_type& dst,
- const volatile data_type& src) const
- {
- if (src < dst) {
- dst = src;
+ PASTIS_INLINE
+ void operator()(const index_type& i, data_type& data) const
+ {
+ if (m_array[i] < data) {
+ data = m_array[i];
+ }
}
- }
- PASTIS_INLINE
- void init(data_type& value) const
- {
- value = std::numeric_limits<data_type>::max();
- }
+ PASTIS_INLINE
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ if (src < dst) {
+ dst = src;
+ }
+ }
- PASTIS_INLINE
- Min(const ArrayType& array)
- : m_array(array)
- {
- ;
- }
+ PASTIS_INLINE
+ void init(data_type& value) const
+ {
+ value = std::numeric_limits<data_type>::max();
+ }
- PASTIS_INLINE
- ~Min() = default;
-};
+ PASTIS_INLINE
+ ArrayMin(const ArrayType& array)
+ : m_array(array)
+ {
+ ;
+ }
+ PASTIS_INLINE
+ ~ArrayMin() = default;
+ };
-template <typename ArrayType>
-class Max
+
+ return ArrayMin(array);
+}
+
+template <typename DataType,
+ template <typename> typename ArrayT>
+std::remove_const_t<DataType>
+max(const ArrayT<DataType>& array)
{
- private:
- const ArrayType& m_array;
+ using ArrayType = ArrayT<DataType>;
using data_type = std::remove_const_t<typename ArrayType::data_type>;
using index_type = typename ArrayType::index_type;
- public:
- PASTIS_INLINE
- operator data_type()
+ class ArrayMax
{
- data_type reduced_value;
- parallel_reduce(m_array.size(), *this, reduced_value);
- return reduced_value;
- }
+ private:
+ const ArrayType m_array;
+
+ public:
+ PASTIS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_array.size(), *this, reduced_value);
+ return reduced_value;
+ }
- PASTIS_INLINE
- void operator()(const index_type& i, data_type& data) const
- {
- if (m_array[i] > data) {
- data = m_array[i];
+ PASTIS_INLINE
+ void operator()(const index_type& i, data_type& data) const
+ {
+ if (m_array[i] > data) {
+ data = m_array[i];
+ }
}
- }
- PASTIS_INLINE
- void join(volatile data_type& dst,
- const volatile data_type& src) const
- {
- if (src > dst) {
- dst = src;
+ PASTIS_INLINE
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ if (src > dst) {
+ dst = src;
+ }
}
- }
- PASTIS_INLINE
- void init(data_type& value) const
- {
- value = -std::numeric_limits<data_type>::max();
- }
+ PASTIS_INLINE
+ void init(data_type& value) const
+ {
+ value = std::numeric_limits<data_type>::min();
+ }
- PASTIS_INLINE
- Max(const ArrayType& array)
- : m_array(array)
- {
- ;
- }
+ PASTIS_INLINE
+ ArrayMax(const ArrayType& array)
+ : m_array(array)
+ {
+ ;
+ }
- PASTIS_INLINE
- ~Max() = default;
-};
+ PASTIS_INLINE
+ ~ArrayMax() = default;
+ };
+ return ArrayMax(array);
+}
-template <typename DataType>
-class Sum
+template <typename DataType,
+ template <typename> typename ArrayT>
+std::remove_const_t<DataType>
+sum(const ArrayT<DataType>& array)
{
- private:
- using data_type = DataType;
- const Array<data_type>& m_array;
- using index_type = typename Array<DataType>::index_type;
-
- public:
- PASTIS_INLINE
- operator data_type()
- {
- data_type reduced_value;
- parallel_reduce(m_array.size(), *this, reduced_value);
- return reduced_value;
- }
+ using ArrayType = ArrayT<DataType>;
+ using data_type = std::remove_const_t<DataType>;
+ using index_type = typename ArrayType::index_type;
- PASTIS_INLINE
- void operator()(const index_type& i, data_type& data) const
+ class ArraySum
{
- data += m_array[i];
- }
+ private:
+ const ArrayType& m_array;
+
+ public:
+ PASTIS_INLINE
+ operator data_type()
+ {
+ data_type reduced_value;
+ parallel_reduce(m_array.size(), *this, reduced_value);
+ return reduced_value;
+ }
- PASTIS_INLINE
- void join(volatile data_type& dst,
- const volatile data_type& src) const
- {
- dst += src;
- }
+ PASTIS_INLINE
+ void operator()(const index_type& i, data_type& data) const
+ {
+ data += m_array[i];
+ }
- PASTIS_INLINE
- void init(data_type& value) const
- {
- value = 0;
- }
+ PASTIS_INLINE
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ dst += src;
+ }
- PASTIS_INLINE
- Sum(const Array<DataType>& array)
- : m_array(array)
- {
- ;
- }
+ PASTIS_INLINE
+ void init(data_type& value) const
+ {
+ if constexpr (std::is_arithmetic_v<data_type>) {
+ value = 0;
+ } else {
+ value = zero;
+ }
+ }
- PASTIS_INLINE
- ~Sum() = default;
-};
+ PASTIS_INLINE
+ ArraySum(const ArrayType& array)
+ : m_array(array)
+ {
+ ;
+ }
+
+ PASTIS_INLINE
+ ~ArraySum() = default;
+ };
+
+ return ArraySum(array);
+}
template <template <typename ...SourceT> typename SourceArray,
template <typename ...ImageT> typename ImageArray,
diff --git a/src/utils/Messenger.hpp b/src/utils/Messenger.hpp
index 12c237893..328be198b 100644
--- a/src/utils/Messenger.hpp
+++ b/src/utils/Messenger.hpp
@@ -3,6 +3,7 @@
#include <PastisMacros.hpp>
#include <PastisAssert.hpp>
+#include <PastisOStream.hpp>
#include <Array.hpp>
#include <CastArray.hpp>
@@ -328,6 +329,29 @@ class Messenger
#endif // PASTIS_HAS_MPI
}
+ template <typename DataType>
+ DataType allReduceSum(const DataType& data) const
+ {
+#ifdef PASTIS_HAS_MPI
+ static_assert(not std::is_const_v<DataType>);
+ static_assert(std::is_arithmetic_v<DataType>);
+
+ MPI_Datatype mpi_datatype
+ = Messenger::helper::mpiType<DataType>();
+
+ perr() << __FILE__ << ':' << __LINE__ << ": Implementation not finished\n";
+ std::terminate();
+
+
+ DataType data_sum = data;
+ MPI_Allreduce(&data, &data_sum, 1, mpi_datatype, MPI_SUM, MPI_COMM_WORLD);
+
+ return data_sum;
+#else // PASTIS_HAS_MPI
+ return data;
+#endif // PASTIS_HAS_MPI
+ }
+
template <typename DataType>
PASTIS_INLINE
Array<DataType>
@@ -538,6 +562,13 @@ DataType allReduceMin(const DataType& data)
return messenger().allReduceMin(data);
}
+template <typename DataType>
+PASTIS_INLINE
+DataType allReduceSum(const DataType& data)
+{
+ return messenger().allReduceSum(data);
+}
+
template <typename DataType>
PASTIS_INLINE
Array<DataType>
diff --git a/tests/test_ArrayUtils.cpp b/tests/test_ArrayUtils.cpp
index 80ead5291..9d5666785 100644
--- a/tests/test_ArrayUtils.cpp
+++ b/tests/test_ArrayUtils.cpp
@@ -4,15 +4,15 @@
#include <Array.hpp>
#include <ArrayUtils.hpp>
+#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
+
// Instantiate to ensure full coverage is performed
template class Array<int>;
TEST_CASE("ArrayUtils", "[utils]") {
- SECTION("checking for reductions") {
- using ArrayType = Array<int>;
- using data_type = ArrayType::data_type;
-
+ SECTION("checking for Array reductions") {
Array<int> a(10);
a[0] =13;
a[1] = 1;
@@ -26,71 +26,50 @@ TEST_CASE("ArrayUtils", "[utils]") {
a[9] = 9;
SECTION("Min") {
- Min r_min(a);
- data_type data;
-
- r_min.init(data);
- REQUIRE(data == std::numeric_limits<data_type>::max());
-
- r_min(2,data);
- REQUIRE(data == 8);
-
- r_min(9,data);
- REQUIRE(data == 8);
-
- r_min(5,data);
- REQUIRE(data == -1);
-
- {
- data = 0;
- data_type r_value{3};
- r_min.join(data, r_value);
- REQUIRE(data == 0);
- }
-
- {
- data = 0;
- data_type r_value{-5};
- r_min.join(data, r_value);
- REQUIRE(data == -5);
- }
-
- REQUIRE((r_min == -3));
-
- REQUIRE((Min(a) == -3));
+ REQUIRE((min(a) == -3));
}
- SECTION("ReduceMax") {
- Max r_max(a);
- data_type data;
-
- r_max.init(data);
- REQUIRE(data == -std::numeric_limits<data_type>::max());
-
- r_max(3,data);
- REQUIRE(data == -3);
- r_max(5,data);
- REQUIRE(data == -1);
- r_max(8,data);
- REQUIRE(data == 12);
-
- {
- data = 0;
- data_type r_value{3};
- r_max.join(data, r_value);
- REQUIRE(data == 3);
- }
+ SECTION("Max") {
+ REQUIRE((max(a) == 23));
+ }
- {
- data = 0;
- data_type r_value{-5};
- r_max.join(data, r_value);
- REQUIRE(data == 0);
- }
+ SECTION("Sum") {
+ REQUIRE((sum(a) == 75));
+ }
- REQUIRE((r_max == 23));
+ SECTION("TinyVector Sum") {
+ using N2 = TinyVector<2,int>;
+ Array<N2> b(10);
+ b[0] ={13, 2};
+ b[1] = {1, 3};
+ b[2] = {8, -2};
+ b[3] ={-3, 2};
+ b[4] ={23, 4};
+ b[5] ={-1, -3};
+ b[6] ={13, 17};
+ b[7] = {0, 9};
+ b[8] ={12, 13};
+ b[9] = {9, -17};
+
+ REQUIRE((sum(b) == N2{75, 28}));
+ }
- REQUIRE((Max(a) == 23));
+ SECTION("TinyMatrix Sum") {
+ using N22 = TinyMatrix<2,int>;
+ Array<N22> b(10);
+ b[0] ={13, 2, 0, 1};
+ b[1] = {1, 3, 6, 3};
+ b[2] = {8, -2, -1, 21};
+ b[3] ={-3, 2, 5, 12};
+ b[4] ={23, 4, 7, 1};
+ b[5] ={-1, -3, 33, 11};
+ b[6] ={13, 17, 12, 13};
+ b[7] = {0, 9, 1, 14};
+ b[8] ={12, 13, -3,-71};
+ b[9] = {9,-17, 0, 16};
+
+ REQUIRE((sum(b) == N22{75, 28, 60, 21}));
}
+
}
}
--
GitLab