diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 453063b76a8977b28a6ff28f7264dedf2b930444..bafd84bdc281f58c9e62e4fd312fe8dce183838b 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 c39bd5ccc75abea27ce652fbb1f1ab72157d5b7b..417d26fadfa8badf208410910661dddb5f163d09 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 0ff0f36ad3594c26567314c305e3f3b84f28ee0b..f3aba4e18f18c308e0f16654221393a0238c3927 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 0000000000000000000000000000000000000000..1381373975405a1bb669b6eb717501d0847b303b
--- /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 f26e0b4fdedecd1c04299dbd16a44f6d87bcccc6..dc3d1bb12018184c6b75a3997a1f889d1f3cc3d3 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 3cb320957e2175da47059eeb215dd8e1a6e34592..fe5e8bf182d0ba02c08acb9f8d406e96244c7fad 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 6392eddd630bf4e22ce613348051aac93ff6c217..51e42309473cd069f0a8946277ae0630665ae686 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 12c2378939c2512d55815fc8dd57e27ead44de83..328be198bdf8e52ebf2f67ba41fbed6c6162cb01 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 80ead52915147e40a0ef5236cfe7972d079168e2..9d56667850c6ca727a154c6ef20fa40b6ba01ac1 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}));
}
+
}
}