From 0e613992698a9f856382859a5e12db5531726d19 Mon Sep 17 00:00:00 2001 From: dhinf Date: Mon, 10 Oct 2022 13:34:54 +0200 Subject: [PATCH] DASH vs ScaFES implementations --- heat_equation/heatDASH_2d_custom_scafes.cpp | 182 +++++++++++++++++++ heat_equation/heatDASH_3d_custom_scafes.cpp | 188 ++++++++++++++++++++ heat_equation/heatScafes_2d.cpp | 41 +++++ heat_equation/heatScafes_2d.hpp | 179 +++++++++++++++++++ heat_equation/heatScafes_3d.cpp | 41 +++++ 5 files changed, 631 insertions(+) create mode 100644 heat_equation/heatDASH_2d_custom_scafes.cpp create mode 100644 heat_equation/heatDASH_3d_custom_scafes.cpp create mode 100644 heat_equation/heatScafes_2d.cpp create mode 100644 heat_equation/heatScafes_2d.hpp create mode 100644 heat_equation/heatScafes_3d.cpp diff --git a/heat_equation/heatDASH_2d_custom_scafes.cpp b/heat_equation/heatDASH_2d_custom_scafes.cpp new file mode 100644 index 0000000..94bd82c --- /dev/null +++ b/heat_equation/heatDASH_2d_custom_scafes.cpp @@ -0,0 +1,182 @@ +#include +#include +#include "dash/Init.h" +#include "dash/halo/Halo.h" +#include "dash/halo/Types.h" + +#include "minimon.h" + +using namespace std; + +constexpr dash::dim_t DIMENSION = 2; +constexpr size_t INIT_WIDTH = 5; + +using PatternT = dash::Pattern; +using MatrixT = dash::Matrix; +using StencilT = dash::halo::StencilPoint; +using StencilSpecT = dash::halo::StencilSpec; +using GlobBoundSpecT = dash::halo::GlobalBoundarySpec; +using HaloMatrixWrapperT = dash::halo::HaloMatrixWrapper; + +using array_t = dash::Array; + +MiniMon minimon{}; + +constexpr double dx = 1.0; +constexpr double dy = 1.0; +constexpr double dt = 0.05; +constexpr double k = 1.0; + +double calcEnergy(const MatrixT& m, array_t &a) { + *a.lbegin() = std::accumulate( m.lbegin(), m.lend(), 0.0); + a.barrier(); + double energy = 0.0; + if(dash::myid() == 0) + energy = std::accumulate(a.begin(), a.end(), 0.0); + a.barrier(); + + return energy; +} + +int main(int argc, char *argv[]) +{ + minimon.enter(); + + if (argc < 3) { + cerr << "Not enough arguments ./ matrix_ext iterations" << endl; + return 1; + } + + auto matrix_ext = std::atoi(argv[1]); + auto iterations = std::atoi(argv[2]); + + minimon.enter(); + dash::init(&argc, &argv); + minimon.leave("initialize"); + + auto myid = dash::myid(); + auto ranks = dash::size(); + + dash::DistributionSpec dist(dash::BLOCKED, dash::BLOCKED); + dash::TeamSpec tspec; + tspec.balance_extents(); + PatternT pattern(dash::SizeSpec(matrix_ext, matrix_ext),dist, tspec, dash::Team::All()); + auto src_matrix = MatrixT(pattern); + auto dst_matrix = MatrixT(pattern); + + StencilSpecT stencil_spec( StencilT(-1, 0), StencilT(1, 0), + StencilT( 0,-1), StencilT(0, 1)); + + // Periodic/cyclic global boundary values for both dimensions + GlobBoundSpecT bound_spec(dash::halo::BoundaryProp::CUSTOM,dash::halo::BoundaryProp::CUSTOM); + + // HaloWrapper for source and destination partitions + HaloMatrixWrapperT src_halo(src_matrix, bound_spec, stencil_spec); + HaloMatrixWrapperT dst_halo(dst_matrix, bound_spec, stencil_spec); + auto src_halo_ptr = &src_halo; + auto dst_halo_ptr = &dst_halo; + + // Stencil specific operator for both partitions + auto src_stencil_op = src_halo_ptr->stencil_operator(stencil_spec); + auto dst_stencil_op = dst_halo_ptr->stencil_operator(stencil_spec); + auto src_op_ptr = &src_stencil_op; + auto dst_op_ptr = &dst_stencil_op; + + for(auto i = 0; i < src_matrix.local.extent(0); ++i) { + for(auto j = 0; j < src_matrix.local.extent(1); ++j) { + if(i < INIT_WIDTH && j < INIT_WIDTH) { + src_matrix.local[i][j] = 1; + dst_matrix.local[i][j] = 1; + } else { + src_matrix.local[i][j] = 0; + dst_matrix.local[i][j] = 0; + } + } + } + + src_halo.set_custom_halos([](const auto& coords) { return 0;}); + dst_halo.set_custom_halos([](const auto& coords) { return 0;}); + + + // initial total energy + auto energy = array_t(ranks); + double initEnergy = calcEnergy(src_halo_ptr->matrix(), energy); + + src_matrix.barrier(); + minimon.enter(); + + for (auto i = 0; i < iterations; ++i) { + minimon.enter(); + auto dst_matrix_lbegin = dst_halo_ptr->matrix().lbegin(); + + minimon.enter(); + src_halo_ptr->update_async(); // start asynchronous Halo data exchange + minimon.leave("async"); + + minimon.enter(); + // Calculation of all inner partition elements + auto iend = src_op_ptr->inner.end(); + for(auto it = src_op_ptr->inner.begin(); it != iend; ++it) { + auto center = *it; + double dtheta = (it.value_at(0) + it.value_at(1) - 2 * center) / (dx * dx) + + (it.value_at(2) + it.value_at(3) - 2 * center) / (dy * dy); + dst_matrix_lbegin[it.lpos()] = center + k * dtheta * dt; + } + minimon.leave("inner"); + + minimon.enter(); + src_halo_ptr->wait(); // Wait until all Halo data exchanges are finished + minimon.leave("wait"); + + minimon.enter(); + // Calculation of all boundary partition elements + auto bend = src_op_ptr->boundary.end(); + for(auto it = src_op_ptr->boundary.begin(); it != bend; ++it) { + auto center = *it; + double dtheta = (it.value_at(0) + it.value_at(1) - 2 * center) / (dx * dx) + + (it.value_at(2) + it.value_at(3) - 2 * center) / (dy * dy); + dst_matrix_lbegin[it.lpos()] = center + k * dtheta * dt; + } + minimon.leave("boundary"); + + // swap source and destination partitions and operators + std::swap(src_halo_ptr, dst_halo_ptr); + std::swap(src_op_ptr, dst_op_ptr); + minimon.leave("calc iter"); + } + + dash::barrier(); + minimon.leave("calc total"); + + // final total energy + double endEnergy = calcEnergy(src_halo_ptr->matrix(), energy); + + // Output + if (myid == 0) { + cout << fixed; + cout.precision(5); + cout << "InitEnergy=" << initEnergy << endl; + cout << "EndEnergy=" << endEnergy << endl; + cout << "DiffEnergy=" << endEnergy - initEnergy << endl; + cout << "Matrixspec: " << matrix_ext << " x " << matrix_ext << endl; + cout << "Iterations: " << iterations << endl; + cout.flush(); + } + + dash::Team::All().barrier(); + minimon.leave("total"); + + if(myid == 0) + std::cout << "# unit_id;function_name;num_calls;avg_runtime;min_runtime;max_runtime" + << std::endl; + + for(auto i = 0; i < dash::size(); ++i) { + if(i == myid) + minimon.print(myid); + dash::Team::All().barrier(); + } + + dash::finalize(); + + return 0; +} diff --git a/heat_equation/heatDASH_3d_custom_scafes.cpp b/heat_equation/heatDASH_3d_custom_scafes.cpp new file mode 100644 index 0000000..2f48b88 --- /dev/null +++ b/heat_equation/heatDASH_3d_custom_scafes.cpp @@ -0,0 +1,188 @@ +#include +#include +#include "dash/Init.h" +#include "dash/halo/Halo.h" +#include "dash/halo/Types.h" + +#include "minimon.h" + +using namespace std; + +constexpr dash::dim_t DIMENSION = 3; +constexpr size_t INIT_WIDTH = 5; + +using PatternT = dash::Pattern; +using MatrixT = dash::Matrix; +using StencilT = dash::halo::StencilPoint; +using StencilSpecT = dash::halo::StencilSpec; +using GlobBoundSpecT = dash::halo::GlobalBoundarySpec; +using HaloMatrixWrapperT = dash::halo::HaloMatrixWrapper; + +using array_t = dash::Array; + +MiniMon minimon{}; + +constexpr double dx = 1.0; +constexpr double dy = 1.0; +constexpr double dz = 1.0; +constexpr double dt = 0.05; +constexpr double k = 1.0; + +double calcEnergy(const MatrixT& m, array_t &a) { + *a.lbegin() = std::accumulate( m.lbegin(), m.lend(), 0.0); + a.barrier(); + double energy = 0.0; + if(dash::myid() == 0) + energy = std::accumulate(a.begin(), a.end(), 0.0); + a.barrier(); + + return energy; +} + +int main(int argc, char *argv[]) +{ + minimon.enter(); + + if (argc < 3) { + cerr << "Not enough arguments ./ matrix_ext iterations" << endl; + return 1; + } + + auto matrix_ext = std::atoi(argv[1]); + auto iterations = std::atoi(argv[2]); + + minimon.enter(); + dash::init(&argc, &argv); + minimon.leave("initialize"); + + auto myid = dash::myid(); + auto ranks = dash::size(); + + dash::DistributionSpec dist(dash::BLOCKED, dash::BLOCKED, dash::BLOCKED); + dash::TeamSpec tspec; + tspec.balance_extents(); + PatternT pattern(dash::SizeSpec(matrix_ext, matrix_ext, matrix_ext),dist, tspec, dash::Team::All()); + auto src_matrix = MatrixT(pattern); + auto dst_matrix = MatrixT(pattern); + + StencilSpecT stencil_spec( StencilT(-1, 0, 0), StencilT(1, 0, 0), + StencilT( 0,-1, 0), StencilT(0, 1, 0), + StencilT( 0, 0,-1), StencilT(0, 0, 1)); + + // Periodic/cyclic global boundary values for both dimensions + GlobBoundSpecT bound_spec(dash::halo::BoundaryProp::CUSTOM,dash::halo::BoundaryProp::CUSTOM,dash::halo::BoundaryProp::CUSTOM); + + // HaloWrapper for source and destination partitions + HaloMatrixWrapperT src_halo(src_matrix, bound_spec, stencil_spec); + HaloMatrixWrapperT dst_halo(dst_matrix, bound_spec, stencil_spec); + auto src_halo_ptr = &src_halo; + auto dst_halo_ptr = &dst_halo; + + // Stencil specific operator for both partitions + auto src_stencil_op = src_halo_ptr->stencil_operator(stencil_spec); + auto dst_stencil_op = dst_halo_ptr->stencil_operator(stencil_spec); + auto src_op_ptr = &src_stencil_op; + auto dst_op_ptr = &dst_stencil_op; + + for(auto i = 0; i < src_matrix.local.extent(0); ++i) { + for(auto j = 0; j < src_matrix.local.extent(1); ++j) { + for(auto k = 0; k < src_matrix.local.extent(2); ++k) { + if(i < INIT_WIDTH && j < INIT_WIDTH && k < INIT_WIDTH) { + src_matrix.local[i][j][k] = 1; + dst_matrix.local[i][j][k] = 1; + } else { + src_matrix.local[i][j][k] = 0; + dst_matrix.local[i][j][k] = 0; + } + } + } + } + + src_halo.set_custom_halos([](const auto& coords) { return 0;}); + dst_halo.set_custom_halos([](const auto& coords) { return 0;}); + + + // initial total energy + auto energy = array_t(ranks); + double initEnergy = calcEnergy(src_halo_ptr->matrix(), energy); + + src_matrix.barrier(); + minimon.enter(); + + for (auto i = 0; i < iterations; ++i) { + minimon.enter(); + auto dst_matrix_lbegin = dst_halo_ptr->matrix().lbegin(); + + minimon.enter(); + src_halo_ptr->update_async(); // start asynchronous Halo data exchange + minimon.leave("async"); + + minimon.enter(); + // Calculation of all inner partition elements + auto iend = src_op_ptr->inner.end(); + for(auto it = src_op_ptr->inner.begin(); it != iend; ++it) { + auto center = *it; + double dtheta = (it.value_at(0) + it.value_at(1) - 2 * center) / (dx * dx) + + (it.value_at(2) + it.value_at(3) - 2 * center) / (dy * dy) + + (it.value_at(4) + it.value_at(5) - 2 * center) / (dz * dz); + dst_matrix_lbegin[it.lpos()] = center + k * dtheta * dt; + } + minimon.leave("inner"); + + minimon.enter(); + src_halo_ptr->wait(); // Wait until all Halo data exchanges are finished + minimon.leave("wait"); + + minimon.enter(); + // Calculation of all boundary partition elements + auto bend = src_op_ptr->boundary.end(); + for(auto it = src_op_ptr->boundary.begin(); it != bend; ++it) { + auto center = *it; + double dtheta = (it.value_at(0) + it.value_at(1) - 2 * center) / (dx * dx) + + (it.value_at(2) + it.value_at(3) - 2 * center) / (dy * dy) + + (it.value_at(4) + it.value_at(5) - 2 * center) / (dz * dz); + dst_matrix_lbegin[it.lpos()] = center + k * dtheta * dt; + } + minimon.leave("boundary"); + + // swap source and destination partitions and operators + std::swap(src_halo_ptr, dst_halo_ptr); + std::swap(src_op_ptr, dst_op_ptr); + minimon.leave("calc iter"); + } + + dash::barrier(); + minimon.leave("calc total"); + + // final total energy + double endEnergy = calcEnergy(src_halo_ptr->matrix(), energy); + + // Output + if (myid == 0) { + cout << fixed; + cout.precision(5); + cout << "InitEnergy=" << initEnergy << endl; + cout << "EndEnergy=" << endEnergy << endl; + cout << "DiffEnergy=" << endEnergy - initEnergy << endl; + cout << "Matrixspec: " << matrix_ext << " x " << matrix_ext <<" x " << matrix_ext << endl; + cout << "Iterations: " << iterations << endl; + cout.flush(); + } + + dash::Team::All().barrier(); + minimon.leave("total"); + + if(myid == 0) + std::cout << "# unit_id;function_name;num_calls;avg_runtime;min_runtime;max_runtime" + << std::endl; + + for(auto i = 0; i < dash::size(); ++i) { + if(i == myid) + minimon.print(myid); + dash::Team::All().barrier(); + } + + dash::finalize(); + + return 0; +} diff --git a/heat_equation/heatScafes_2d.cpp b/heat_equation/heatScafes_2d.cpp new file mode 100644 index 0000000..2375fe5 --- /dev/null +++ b/heat_equation/heatScafes_2d.cpp @@ -0,0 +1,41 @@ +#include "ScaFES.hpp" +#include "heatScafes_2d.hpp" + +#include "minimon.h" + +/** Space dimension of problem. */ +const int DIM = 2; + +/** Main program for HeatEqnFDM. */ +int main(int argc, char *argv[]) { + MiniMon minimon{}; + minimon.enter(); + ScaFES::Parameters paramsCl(argc, argv); + ScaFES::GridGlobal gg(paramsCl); + + std::vector nameDatafield(1); + nameDatafield[0] = "F"; + std::vector stencilWidth(1); + stencilWidth[0] = 1; + std::vector isKnownDf(1); + isKnownDf[0] = false; + std::vector nLayers(1); + nLayers[0] = 0; + std::vector defaultValue(1); + defaultValue[0] = 0.0; + std::vector writeToFile(1, ScaFES::WriteHowOften::LIKE_GIVEN_AT_CL); + std::vector computeError(1); + computeError[0] = false; + + std::vector geomparamsInit; + + HeatEqnFDM ppp(paramsCl, gg, false, nameDatafield, stencilWidth, + isKnownDf, nLayers, defaultValue, writeToFile, + computeError, geomparamsInit); + ppp.iterateOverTime(); + minimon.leave("total"); + + minimon.print(0); + + return 0; +} diff --git a/heat_equation/heatScafes_2d.hpp b/heat_equation/heatScafes_2d.hpp new file mode 100644 index 0000000..72f57e5 --- /dev/null +++ b/heat_equation/heatScafes_2d.hpp @@ -0,0 +1,179 @@ +#include + +#include "ScaFES.hpp" + +constexpr size_t INIT_WIDTH = 5; + +template +class HeatEqnFDM : public ScaFES::Problem, CT, DIM> { + public: + /** All fields which are related to the underlying problem + * are added in terms of an entry of the parameters of + * type \c std::vector. + * @param params Set of ScaFES parameters. + * @param gg Global grid. + * @param useLeapfrog Should the leap frog scheme be used? + * @param nameDatafield Name of the fields. + * @param stencilWidth Stencil width of the fields. + * @param isKnownDf Is the data field are known or unknown one? + * @param nLayers Number of layers at the global boundary. + * @param defaultValue Default value of fields. + * @param writeToFile How often should the data field be written to file. + * @param computeError Should the Linf error between the numerical + * and exact solution be computed? + * @param geomparamsInit Initial guess of geometrical parameters. + */ + HeatEqnFDM(ScaFES::Parameters const& params, + ScaFES::GridGlobal const& gg, + bool useLeapfrog, + std::vector const& nameDatafield, + std::vector const& stencilWidth, + std::vector const& isKnownDf, + std::vector const& nLayers = std::vector(), + std::vector const& defaultValue = std::vector(), + std::vector const& writeToFile + = std::vector(), + std::vector const& computeError = std::vector(), + std::vector const& geomparamsInit = std::vector() ) + : ScaFES::Problem, CT, DIM>(params, gg, useLeapfrog, + nameDatafield, stencilWidth, + isKnownDf, nLayers, + defaultValue, writeToFile, + computeError, geomparamsInit) + { + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + } + + /** Evaluates all fields at one given global inner grid node. + * @param vNew Set of all fields. + * @param idxNode Index of given grid node. + */ + void evalInner(std::vector< ScaFES::DataField >& vNew, + ScaFES::Ntuple const& idxNode, + int const& /*timestep*/) { + } + /** Evaluates all fields at one given global border grid node. + * @param vNew Set of all fields. + * @param idxNode Index of given grid node. + * @param timestep Given time step. + */ + void evalBorder(std::vector< ScaFES::DataField >& vNew, + ScaFES::Ntuple const& idxNode, + int const& timestep) { + } + + /** Initializes all unknown fields at one given global inner grid node. + * @param vNew Set of all unknown fields (return value). + * @param idxNode Index of given grid node. + */ + template + void initInner(std::vector< ScaFES::DataField >& vNew, + std::vector const& /*vOld*/, + ScaFES::Ntuple const& idxNode, + int const& timestep) { + const auto& partition = this->globalGrid().partition(rank); + auto first = partition.idxNodeFirstSub(); + if((idxNode[0] - first[0]) < INIT_WIDTH && (idxNode[1] - first[1]) < INIT_WIDTH) { + vNew[0](idxNode) = 1; + } else { + vNew[0](idxNode) = 0; + } + } + /** Initializes all unknown fields at one given global border grid node. + * @param vNew Set of all unknown fields (return value). + * @param vOld Set of all given fields. + * @param idxNode Index of given grid node. + * @param timestep Given tim step. + */ + template + void initBorder(std::vector< ScaFES::DataField >& vNew, + std::vector const& vOld, + ScaFES::Ntuple const& idxNode, + int const& timestep) { + this->template initInner(vNew, vOld, idxNode, timestep); + } + + /** Updates all unknown fields at one given global inner grid node. + * @param vNew Set of all unknown fields at new time step (return value). + * @param vOld Set of all unknown fields at old time step. + * @param idxNode Index of given grid node. + */ + template + void updateInner(std::vector>& vNew, + std::vector> const& vOld, + ScaFES::Ntuple const& idxNode, + int const& timestep) { + auto center = vOld[0](idxNode); + double dtheta = (vOld[0](this->connect(idxNode, 0)) + vOld[0](this->connect(idxNode, 1)) - 2 * center) / (dx * dx) + + (vOld[0](this->connect(idxNode, 2)) + vOld[0](this->connect(idxNode, 3)) - 2 * center) / (dy * dy); + + vNew[0](idxNode) = center + k * dtheta * dt; + } + /** Updates all unknown fields at one given global border grid node. + * @param vNew Set of all unknown fields at new time step (return value). + * @param idxNode Index of given grid node. + */ + template + void updateBorder(std::vector>& vNew, + std::vector>const& vOld, + ScaFES::Ntuple const& idxNode, + int const& timestep) { + auto center = vOld[0](idxNode); + std::array stencil_values{}; + for(size_t i = 0; i < stencil_values.size(); ++i) { + auto idxNode_tmp = this->connect(idxNode, i); + if(idxNode_tmp[0] < 0 || idxNode_tmp < 0) { + continue; + } + stencil_values[i] = vOld[0](idxNode_tmp); + } + double dtheta = (stencil_values[0] + stencil_values[1] - 2 * center) / (dx * dx) + + (stencil_values[2] + stencil_values[3] - 2 * center) / (dy * dy); + + vNew[0](idxNode) = center + k * dtheta * dt; + } + + /** Updates (2nd cycle) all unknown fields at one given global inner grid node. + * \remarks Only important if leap frog scheme is used. + */ + template + void updateInner2(std::vector>&, + std::vector> const&, + ScaFES::Ntuple const&, + int const&) { } + + /** Updates (2nd cycle) all unknown fields at one given global border + * grid node. + * \remarks Only important if leap frog scheme is used. + */ + template + void updateBorder2(std::vector>&, + std::vector>const&, + ScaFES::Ntuple const&, + int const&) { } + + template + double calc_energy(const std::vector>& v) { + const auto& partition = this->globalGrid().partition(rank); + auto first = partition.idxNodeFirstSub(); + auto last = partition.idxNodeLastSub(); + TT energy{}; + for(int i = first[0]; i <= last[0]; ++i) { + for(int j = first[1]; j <= last[1]; ++j) { + energy += v[0](ScaFES::Ntuple(i,j)); + } + } + + double energy_final = 0; + MPI_Reduce(&energy, &energy_final, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + + return energy_final; + } +private: + static constexpr double dx = 1.0; + static constexpr double dy = 1.0; + static constexpr double dt = 0.05; + static constexpr double k = 1.0; + + int rank = 0; +}; diff --git a/heat_equation/heatScafes_3d.cpp b/heat_equation/heatScafes_3d.cpp new file mode 100644 index 0000000..8d06f45 --- /dev/null +++ b/heat_equation/heatScafes_3d.cpp @@ -0,0 +1,41 @@ +#include "ScaFES.hpp" +#include "heatScafes_3d.hpp" + +#include "minimon.h" + +/** Space dimension of problem. */ +const int DIM = 3; + +/** Main program for HeatEqnFDM. */ +int main(int argc, char *argv[]) { + MiniMon minimon{}; + minimon.enter(); + ScaFES::Parameters paramsCl(argc, argv); + ScaFES::GridGlobal gg(paramsCl); + + std::vector nameDatafield(1); + nameDatafield[0] = "F"; + std::vector stencilWidth(1); + stencilWidth[0] = 1; + std::vector isKnownDf(1); + isKnownDf[0] = false; + std::vector nLayers(1); + nLayers[0] = 0; + std::vector defaultValue(1); + defaultValue[0] = 0.0; + std::vector writeToFile(1, ScaFES::WriteHowOften::LIKE_GIVEN_AT_CL); + std::vector computeError(1); + computeError[0] = false; + + std::vector geomparamsInit; + + HeatEqnFDM ppp(paramsCl, gg, false, nameDatafield, stencilWidth, + isKnownDf, nLayers, defaultValue, writeToFile, + computeError, geomparamsInit); + ppp.iterateOverTime(); + minimon.leave("total"); + + minimon.print(0); + + return 0; +}