diff --git a/src/Field/BareField.h b/src/Field/BareField.h index 72af46b6a8a66fd5eb71ed81a320dd9d058542e0..9609d30f3e6a5681cce1664d81a87a73ac3c7a8d 100644 --- a/src/Field/BareField.h +++ b/src/Field/BareField.h @@ -62,6 +62,7 @@ namespace ippl { //! View type storing the data using view_type = typename detail::ViewType<T, Dim>::view_type; using HostMirror = typename view_type::host_mirror_type; + using policy_type = typename detail::RangePolicy<Dim>::policy_type; /*! A default constructor, which should be used only if the user calls the @@ -144,16 +145,32 @@ namespace ippl { detail::HaloCells<T, Dim>& getHalo() { return halo_m; } // Assignment from a constant. + template <unsigned dim = Dim, std::enable_if_t<(dim == 2), bool> = true> BareField<T, Dim>& operator=(T x); + template <unsigned dim = Dim, std::enable_if_t<(dim == 3), bool> = true> + BareField<T, Dim>& operator=(T x); + + /*! + * Assign an arbitrary BareField expression (2D version) + * @tparam E expression type + * @tparam N size of the expression, this is necessary for running on the + * device since otherwise it does not allocate enough memory + * @param expr is the expression + */ + template <typename E, size_t N, + unsigned dim = Dim, std::enable_if_t<(dim == 2), bool> = true> + BareField<T, Dim>& operator=(const detail::Expression<E, N>& expr); + /*! - * Assign an arbitrary BareField expression + * Assign an arbitrary BareField expression (3D version) * @tparam E expression type * @tparam N size of the expression, this is necessary for running on the * device since otherwise it does not allocate enough memory * @param expr is the expression */ - template <typename E, size_t N> + template <typename E, size_t N, + unsigned dim = Dim, std::enable_if_t<(dim == 3), bool> = true> BareField<T, Dim>& operator=(const detail::Expression<E, N>& expr); /*! @@ -182,17 +199,44 @@ namespace ippl { return Kokkos::create_mirror(dview_m); } + template <unsigned dim = Dim, std::enable_if_t<(dim == 2), bool> = true> + policy_type getRangePolicy(const int nghost = 0) const { + PAssert_LE(nghost, nghost_m); + const size_t shift = nghost_m - nghost; + return policy_type({shift, shift}, + {dview_m.extent(0) - shift, + dview_m.extent(1) - shift}); + } + + template <unsigned dim = Dim, std::enable_if_t<(dim == 3), bool> = true> + policy_type getRangePolicy(const int nghost = 0) const { + PAssert_LE(nghost, nghost_m); + const size_t shift = nghost_m - nghost; + return policy_type({shift, shift, shift}, + {dview_m.extent(0) - shift, + dview_m.extent(1) - shift, + dview_m.extent(2) - shift}); + } + /*! * Print the BareField. * @param out stream */ void write(std::ostream& out = std::cout) const; - T sum(int nghost = 0); - T max(int nghost = 0); - T min(int nghost = 0); - T prod(int nghost = 0); + #define DefineOperation(fun) \ + template <unsigned dim = Dim, \ + std::enable_if_t<(dim == 2), bool> = true> \ + T fun(int nghost = 0); \ + \ + template <unsigned dim = Dim, \ + std::enable_if_t<(dim == 3), bool> = true> \ + T fun(int nghost = 0); + DefineOperation(sum) + DefineOperation(max) + DefineOperation(min) + DefineOperation(prod) private: //! Number of ghost layers on each field boundary @@ -217,5 +261,6 @@ namespace ippl { } #include "Field/BareField.hpp" +#include "Field/BareFieldOperations.hpp" #endif diff --git a/src/Field/BareField.hpp b/src/Field/BareField.hpp index c6d3f11365706767dd4fb616baaad61b6ac56bb7..de696a04f2eb6dd56666782b687a44fff2088627 100644 --- a/src/Field/BareField.hpp +++ b/src/Field/BareField.hpp @@ -58,13 +58,11 @@ namespace ippl { template <typename T, unsigned Dim> void BareField<T, Dim>::setup() { - static_assert(Dim == 3, "Only 3-dimensional fields supported at the momment!"); + static_assert(Dim == 2 || Dim == 3, "Only 2D and 3D fields supported at the momment!"); owned_m = layout_m->getLocalNDIndex(); - if constexpr(Dim == 1) { - this->resize(owned_m[0].length() + 2 * nghost_m); - } else if constexpr(Dim == 2) { + if constexpr(Dim == 2) { this->resize(owned_m[0].length() + 2 * nghost_m, owned_m[1].length() + 2 * nghost_m); } else if constexpr(Dim == 3) { @@ -98,36 +96,67 @@ namespace ippl { } + template <typename T, unsigned Dim> + template <unsigned dim, std::enable_if_t<(dim == 2), bool>> BareField<T, Dim>& BareField<T, Dim>::operator=(T x) { - using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>; + policy_type policy = getRangePolicy(nghost_m); + Kokkos::parallel_for("BareField::operator=(T)", - mdrange_type({0, 0, 0}, - {dview_m.extent(0), - dview_m.extent(1), - dview_m.extent(2) - }), + policy, + KOKKOS_CLASS_LAMBDA(const size_t i, + const size_t j) + { + dview_m(i, j) = x; + }); + return *this; + } + + + template <typename T, unsigned Dim> + template <unsigned dim, std::enable_if_t<(dim == 3), bool>> + BareField<T, Dim>& BareField<T, Dim>::operator=(T x) { + policy_type policy = getRangePolicy(nghost_m); + + Kokkos::parallel_for("BareField::operator=(T)", + policy, + KOKKOS_CLASS_LAMBDA(const size_t i, + const size_t j, + const size_t k) + { + dview_m(i, j, k) = x; + }); + return *this; + } + + + template <typename T, unsigned Dim> + template <typename E, size_t N, unsigned dim, std::enable_if_t<(dim == 2), bool>> + BareField<T, Dim>& BareField<T, Dim>::operator=(const detail::Expression<E, N>& expr) { + using capture_type = detail::CapturedExpression<E, N>; + capture_type expr_ = reinterpret_cast<const capture_type&>(expr); + + policy_type policy = getRangePolicy(nghost_m); + Kokkos::parallel_for("BareField::operator=(const Expression&)", + policy, KOKKOS_CLASS_LAMBDA(const size_t i, - const size_t j, - const size_t k) + const size_t j) { - dview_m(i, j, k) = x; + dview_m(i, j) = expr_(i, j); }); return *this; } template <typename T, unsigned Dim> - template <typename E, size_t N> + template <typename E, size_t N, unsigned dim, std::enable_if_t<(dim == 3), bool>> BareField<T, Dim>& BareField<T, Dim>::operator=(const detail::Expression<E, N>& expr) { using capture_type = detail::CapturedExpression<E, N>; capture_type expr_ = reinterpret_cast<const capture_type&>(expr); - using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>; + + policy_type policy = getRangePolicy(nghost_m); Kokkos::parallel_for("BareField::operator=(const Expression&)", - mdrange_type({nghost_m, nghost_m, nghost_m}, - {dview_m.extent(0) - nghost_m, - dview_m.extent(1) - nghost_m, - dview_m.extent(2) - nghost_m}), + policy, KOKKOS_CLASS_LAMBDA(const size_t i, const size_t j, const size_t k) @@ -145,32 +174,51 @@ namespace ippl { } - #define DefineReduction(fun, name, op, MPI_Op) \ - template <typename T, unsigned Dim> \ - T BareField<T, Dim>::name(int nghost) { \ - PAssert_LE(nghost, nghost_m); \ - T temp = 0.0; \ - const size_t shift = nghost_m - nghost; \ - Kokkos::parallel_reduce("fun", \ - Kokkos::MDRangePolicy<Kokkos::Rank<3>>({shift, shift, shift}, \ - {dview_m.extent(0) - shift, \ - dview_m.extent(1) - shift, \ - dview_m.extent(2) - shift}), \ - KOKKOS_CLASS_LAMBDA(const size_t i, const size_t j, \ - const size_t k, T& valL) { \ - T myVal = dview_m(i, j, k); \ - op; \ - }, Kokkos::fun<T>(temp)); \ - T globaltemp = 0.0; \ - MPI_Datatype type = get_mpi_datatype<T>(temp); \ - MPI_Allreduce(&temp, &globaltemp, 1, type, MPI_Op, Ippl::getComm()); \ - return globaltemp; \ + #define DefineReduction(fun, name, op, MPI_Op) \ + template <typename T, unsigned Dim> \ + template <unsigned dim, std::enable_if_t<(dim == 2), bool>> \ + T BareField<T, Dim>::name(int nghost) { \ + T temp = 0.0; \ + policy_type policy = getRangePolicy(nghost); \ + Kokkos::parallel_reduce("fun", \ + policy, \ + KOKKOS_CLASS_LAMBDA(const size_t i, \ + const size_t j, \ + T& valL) \ + { \ + T myVal = dview_m(i, j); \ + op; \ + }, Kokkos::fun<T>(temp)); \ + T globaltemp = 0.0; \ + MPI_Datatype type = get_mpi_datatype<T>(temp); \ + MPI_Allreduce(&temp, &globaltemp, 1, type, MPI_Op, Ippl::getComm()); \ + return globaltemp; \ + } \ + \ + \ + template <typename T, unsigned Dim> \ + template <unsigned dim, std::enable_if_t<(dim == 3), bool>> \ + T BareField<T, Dim>::name(int nghost) { \ + T temp = 0.0; \ + policy_type policy = getRangePolicy(nghost); \ + Kokkos::parallel_reduce("fun", \ + policy, \ + KOKKOS_CLASS_LAMBDA(const size_t i, \ + const size_t j, \ + const size_t k, \ + T& valL) \ + { \ + T myVal = dview_m(i, j, k); \ + op; \ + }, Kokkos::fun<T>(temp)); \ + T globaltemp = 0.0; \ + MPI_Datatype type = get_mpi_datatype<T>(temp); \ + MPI_Allreduce(&temp, &globaltemp, 1, type, MPI_Op, Ippl::getComm()); \ + return globaltemp; \ } DefineReduction(Sum, sum, valL += myVal, MPI_SUM) DefineReduction(Max, max, if(myVal > valL) valL = myVal, MPI_MAX) DefineReduction(Min, min, if(myVal < valL) valL = myVal, MPI_MIN) DefineReduction(Prod, prod, valL *= myVal, MPI_PROD) - - } diff --git a/src/Field/BareFieldOperations.hpp b/src/Field/BareFieldOperations.hpp new file mode 100644 index 0000000000000000000000000000000000000000..d75a745462b0738070793f4b39f327709927bd3f --- /dev/null +++ b/src/Field/BareFieldOperations.hpp @@ -0,0 +1,138 @@ +// +// File BareFieldOperations +// Norms, and a scalar product for fields +// +// Copyright (c) 2021 Paul Scherrer Institut, Villigen PSI, Switzerland +// All rights reserved +// +// This file is part of IPPL. +// +// IPPL is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// You should have received a copy of the GNU General Public License +// along with IPPL. If not, see <https://www.gnu.org/licenses/>. +// + +namespace ippl { + #define DefineBinaryBareFieldOperation2D(fun, op, kokkos_op, mpi_op)\ + template <typename T, unsigned Dim, \ + std::enable_if_t<(Dim == 2), bool> = true> \ + T fun(const BareField<T, Dim>& bf1, \ + const BareField<T, Dim>& bf2, \ + const int nghost = 0) \ + { \ + T lres = 0; \ + Kokkos::parallel_reduce("BareField::fun", \ + bf1.getRangePolicy(nghost), \ + KOKKOS_LAMBDA(const size_t i, \ + const size_t j, \ + T& val) \ + { \ + op; \ + }, kokkos_op<T>(lres)); \ + T gres = 0; \ + MPI_Datatype type = get_mpi_datatype<T>(lres); \ + MPI_Allreduce(&lres, &gres, 1, type, mpi_op, Ippl::getComm()); \ + return gres; \ + } + + + #define DefineBinaryBareFieldOperation3D(fun, op, kokkos_op, mpi_op)\ + template <typename T, unsigned Dim, \ + std::enable_if_t<(Dim == 3), bool> = true> \ + T fun(const BareField<T, Dim>& bf1, \ + const BareField<T, Dim>& bf2, \ + const int nghost = 0) \ + { \ + T lres = 0; \ + Kokkos::parallel_reduce("BareField::fun", \ + bf1.getRangePolicy(nghost), \ + KOKKOS_LAMBDA(const size_t i, \ + const size_t j, \ + const size_t k, \ + T& val) \ + { \ + op; \ + }, kokkos_op<T>(lres)); \ + T gres = 0; \ + MPI_Datatype type = get_mpi_datatype<T>(lres); \ + MPI_Allreduce(&lres, &gres, 1, type, mpi_op, Ippl::getComm()); \ + return gres; \ + } + + + #define DefineUnaryBareFieldOperation2D(fun, op, kokkos_op, mpi_op) \ + template <typename T, unsigned Dim, \ + std::enable_if_t<(Dim == 2), bool> = true> \ + T fun(const BareField<T, Dim>& bf, \ + const int nghost = 0) \ + { \ + T lres = 0; \ + Kokkos::parallel_reduce("BareField::fun", \ + bf.getRangePolicy(nghost), \ + KOKKOS_LAMBDA(const size_t i, \ + const size_t j, \ + T& val) \ + { \ + op; \ + }, kokkos_op<T>(lres)); \ + T gres = 0; \ + MPI_Datatype type = get_mpi_datatype<T>(lres); \ + MPI_Allreduce(&lres, &gres, 1, type, mpi_op, Ippl::getComm()); \ + return gres; \ + } + + + #define DefineUnaryBareFieldOperation3D(fun, op, kokkos_op, mpi_op) \ + template <typename T, unsigned Dim, \ + std::enable_if_t<(Dim == 3), bool> = true> \ + T fun(const BareField<T, Dim>& bf, \ + const int nghost = 0) \ + { \ + T lres = 0; \ + auto policy = bf.getRangePolicy(nghost); \ + Kokkos::parallel_reduce("BareField::fun", \ + policy, \ + KOKKOS_LAMBDA(const size_t i, \ + const size_t j, \ + const size_t k, \ + T& val) \ + { \ + op; \ + }, kokkos_op<T>(lres)); \ + T gres = 0; \ + MPI_Datatype type = get_mpi_datatype<T>(lres); \ + MPI_Allreduce(&lres, &gres, 1, type, mpi_op, Ippl::getComm()); \ + return gres; \ + } + + + DefineBinaryBareFieldOperation2D(innerProduct, + val += bf1(i, j) * bf2(i, j);, + Kokkos::Sum, MPI_SUM) + DefineBinaryBareFieldOperation3D(innerProduct, + val += bf1(i, j, k) * bf2(i, j, k);, + Kokkos::Sum, MPI_SUM) + + DefineUnaryBareFieldOperation2D(norm1, + val += std::abs(bf(i, j)), + Kokkos::Sum, MPI_SUM) + DefineUnaryBareFieldOperation3D(norm1, + val += std::abs(bf(i, j, k)), + Kokkos::Sum, MPI_SUM) + + DefineUnaryBareFieldOperation2D(normInf, + val = std::abs(bf(i, j)), + Kokkos::Max, MPI_MAX) + DefineUnaryBareFieldOperation3D(normInf, + val = std::abs(bf(i, j, k)), + Kokkos::Max, MPI_MAX) + + template <typename T, unsigned Dim> + T norm2(const BareField<T, Dim>& bf, const int nghost = 0) { + return std::sqrt(innerProduct(bf, bf, nghost)); + } +} diff --git a/src/Field/CMakeLists.txt b/src/Field/CMakeLists.txt index b8f0e5fdf4a02797c533708825e416ba9caf8d3b..7dbe286a002addd50331231de2e07bdec9f36d3c 100644 --- a/src/Field/CMakeLists.txt +++ b/src/Field/CMakeLists.txt @@ -4,6 +4,7 @@ set (_SRCS set (_HDRS BareField.hpp BareField.h + BareFieldOperations.hpp BConds.h BConds.hpp BcTypes.h diff --git a/src/Field/Field.hpp b/src/Field/Field.hpp index e8fcc660583484e509478799cda25a404737b1e8..7c11c1adca5d8e0e8f43e2a807d175ab054268fd 100644 --- a/src/Field/Field.hpp +++ b/src/Field/Field.hpp @@ -30,7 +30,7 @@ namespace ippl { Field<T,Dim,M,C>::Field() : BareField<T, Dim>() , mesh_m(nullptr) - , bc_m() + , bc_m() { } ////////////////////////////////////////////////////////////////////////// @@ -39,7 +39,7 @@ namespace ippl { Field<T,Dim,M,C>::Field(Mesh_t& m, Layout_t& l, int nghost) : BareField<T,Dim>(l, nghost) , mesh_m(&m) - { + { for (unsigned int face=0; face < 2 * Dim; ++face) { bc_m[face] = std::make_shared<NoBcFace<T, Dim>>(face); } diff --git a/src/Field/FieldOperations.hpp b/src/Field/FieldOperations.hpp index 03f20a7748c0069194f6d342b29301e646ba7293..7819eb539eac512d4a46a0e4aad66bd8d747bc42 100644 --- a/src/Field/FieldOperations.hpp +++ b/src/Field/FieldOperations.hpp @@ -1,6 +1,6 @@ // // File FieldOperations -// Differential operators, norms, and a scalar product for fields +// Differential operators for fields // // Copyright (c) 2021 Paul Scherrer Institut, Villigen PSI, Switzerland // All rights reserved @@ -17,88 +17,6 @@ // namespace ippl { - /*! - * Computes the inner product of two fields - * @param f1 first field - * @param f2 second field - * @return Result of f1^T f2 - */ - template <typename T, unsigned Dim> - T innerProduct(const Field<T, Dim>& f1, const Field<T, Dim>& f2) { - T sum = 0; - const int shift = f1.getNghost(); - auto view1 = f1.getView(); - auto view2 = f2.getView(); - Kokkos::parallel_reduce("Field::innerProduct(Field&, Field&)", - Kokkos::MDRangePolicy<Kokkos::Rank<3>>({shift, shift, shift}, { - view1.extent(0) - shift, - view1.extent(1) - shift, - view1.extent(2) - shift}), - KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, T& val) { - val += view1(i, j, k) * view2(i, j, k); - }, - Kokkos::Sum<T>(sum) - ); - T globalSum = 0; - MPI_Datatype type = get_mpi_datatype<T>(sum); - MPI_Allreduce(&sum, &globalSum, 1, type, MPI_SUM, Ippl::getComm()); - return globalSum; - } - - /*! - * Computes the Lp-norm of a field - * @param field field - * @param p desired norm (default 2) - * @return The desired norm of the field - */ - template<typename T, unsigned Dim, class M, class C> - T norm(const Field<T, Dim, M, C>& field, int p = 2) { - T local = 0; - const int shift = field.getNghost(); - auto view = field.getView(); - switch (p) { - case 0: - { - Kokkos::parallel_reduce("Field::norm(0)", - Kokkos::MDRangePolicy<Kokkos::Rank<3>>({shift, shift, shift}, { - view.extent(0) - shift, - view.extent(1) - shift, - view.extent(2) - shift}), - KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, T& val) { - T myVal = std::abs(view(i, j, k)); - if (myVal > val) - val = myVal; - }, - Kokkos::Max<T>(local) - ); - T globalMax = 0; - MPI_Datatype type = get_mpi_datatype<T>(local); - MPI_Allreduce(&local, &globalMax, 1, type, MPI_MAX, Ippl::getComm()); - return globalMax; - } - case 2: - return std::sqrt(innerProduct(field, field)); - default: - { - Kokkos::parallel_reduce("Field::norm(int) general", - Kokkos::MDRangePolicy<Kokkos::Rank<3>>({shift, shift, shift}, { - view.extent(0) - shift, - view.extent(1) - shift, - view.extent(2) - shift}), - KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k, double& val) { - val += std::pow(std::abs(view(i, j, k)), p); - }, - Kokkos::Sum<T>(local) - ); - T globalSum = 0; - MPI_Datatype type = get_mpi_datatype<T>(local); - MPI_Allreduce(&local, &globalSum, 1, type, MPI_SUM, Ippl::getComm()); - return std::pow(globalSum, 1.0 / p); - } - } - } - - /*! * User interface of gradient in three dimensions. * @param u field diff --git a/src/Types/ViewTypes.h b/src/Types/ViewTypes.h index fd782a3fab94e2f238ca330c31fd716a127f7e2f..1e62b47234c6f3ba2c7d5df2ad3b2f0559d3a658 100644 --- a/src/Types/ViewTypes.h +++ b/src/Types/ViewTypes.h @@ -62,6 +62,23 @@ namespace ippl { }; + /*! + * Muldimensional range policies. + */ + template <unsigned Dim> + struct RangePolicy { + typedef Kokkos::MDRangePolicy<Kokkos::Rank<Dim>> policy_type; + }; + + /*! + * Specialized range policy for one dimension. + */ + template <> + struct RangePolicy<1> { + typedef Kokkos::RangePolicy<> policy_type; + }; + + /*! * Empty function for general write. * @tparam T view data type diff --git a/test/field/CMakeLists.txt b/test/field/CMakeLists.txt index 77431b50708cfc42fddc7eeb21d6ce96301c93f1..51e1b7a5edd11e8d119f4cd9f386c654339673b5 100644 --- a/test/field/CMakeLists.txt +++ b/test/field/CMakeLists.txt @@ -63,14 +63,6 @@ target_link_libraries ( ${Boost_LIBRARIES} ) -add_executable (TestFieldNorm TestFieldNorm.cpp) -target_link_libraries ( - TestFieldNorm - ${IPPL_LIBS} - ${MPI_CXX_LIBRARIES} - ${Boost_LIBRARIES} -) - # vi: set et ts=4 sw=4 sts=4: # Local Variables: diff --git a/test/field/TestFieldNorm.cpp b/test/field/TestFieldNorm.cpp deleted file mode 100644 index 6021ecded936133721d63fc6fb3c85970dde1657..0000000000000000000000000000000000000000 --- a/test/field/TestFieldNorm.cpp +++ /dev/null @@ -1,88 +0,0 @@ -// Usage: ./TestFieldNorm (log_2 field size) --info 10 -#include "Ippl.h" -#include "Utility/IpplTimings.h" - -#include <iostream> -#include <typeinfo> - -#include <cstdlib> - -void checkError(double computed, double correct, int N, int p, - Inform& mout, Inform& merr, double tolerance = 1e-16) { - double absError = fabs(computed - correct); - double relError = absError / correct; - mout << "(" << N << ", L" << p << "): " << absError << "," << relError << endl; - if (relError > tolerance) { - merr << "L" << p << " norm for N = " << N << " does not match.\n\tGot " - << computed << ", expected " << correct << ". Relative error: " << relError << endl; - } -} - -int main(int argc, char *argv[]) { - - Ippl ippl(argc,argv); - - int pt = 4; - - if (argc >= 2) { - pt = 1 << (int)strtol(argv[1], NULL, 10); - } - - constexpr unsigned int dim = 3; - - - ippl::Index I(pt); - ippl::NDIndex<dim> owned(I, I, I); - - ippl::e_dim_tag allParallel[dim]; // Specifies SERIAL, PARALLEL dims - for (unsigned int d=0; d<dim; d++) - allParallel[d] = ippl::PARALLEL; - - // all parallel layout, standard domain, normal axis order - ippl::FieldLayout<dim> layout(owned,allParallel); - - double dx = 1.0 / double(pt); - ippl::Vector<double, 3> hx = {dx, dx, dx}; - ippl::Vector<double, 3> origin = {0, 0, 0}; - ippl::UniformCartesian<double, 3> mesh(owned, hx, origin); - - - typedef ippl::Field<double, dim> field_type; - - field_type field(mesh, layout); - - double pi = acos(-1.0); - - field = pi/4; - - double l2 = pow(pt, 1.5) * pi / 4; - double l1 = pow(pt, 3) * pi / 4; - double linf = pi / 4; - - IpplTimings::TimerRef l2Timer = IpplTimings::getTimer("L2"), - l1Timer = IpplTimings::getTimer("L1"), - l0Timer = IpplTimings::getTimer("Max"); - - IpplTimings::startTimer(l2Timer); - double compute_l2 = ippl::norm(field); - IpplTimings::stopTimer(l2Timer); - - IpplTimings::startTimer(l1Timer); - double compute_l1 = ippl::norm(field, 1); - IpplTimings::stopTimer(l1Timer); - - IpplTimings::startTimer(l0Timer); - double compute_linf = ippl::norm(field, 0); - IpplTimings::stopTimer(l0Timer); - - IpplTimings::print("timings.dat"); - - Inform m1("DATA"); - Inform m2("Deviation", std::cerr); - - checkError(compute_l2, l2, pt, 2, m1, m2); - checkError(compute_l1, l1, pt, 1, m1, m2); - checkError(compute_linf, linf, pt, 0, m1, m2); - - return 0; -} diff --git a/unit_tests/Field/CMakeLists.txt b/unit_tests/Field/CMakeLists.txt index ca1673e59ad2cb3fd280aa6b460518d5c00131b2..695388be153f344f3a1b359e5737d0354334b7e4 100644 --- a/unit_tests/Field/CMakeLists.txt +++ b/unit_tests/Field/CMakeLists.txt @@ -17,6 +17,7 @@ add_executable (Field Field.cpp) target_link_libraries ( Field ippl + pthread ${MPI_CXX_LIBRARIES} ${GTEST_BOTH_LIBRARIES} ${Boost_LIBRARIES} @@ -26,6 +27,7 @@ add_executable (FieldBC FieldBC.cpp) target_link_libraries ( FieldBC ippl + pthread ${MPI_CXX_LIBRARIES} ${GTEST_BOTH_LIBRARIES} ${Boost_LIBRARIES} diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp index 35dbcefd2453b7fd70bb49c4904d398b38b981ef..be490629abc73d0a077a9b1022a1005efb31f849 100644 --- a/unit_tests/Field/Field.cpp +++ b/unit_tests/Field/Field.cpp @@ -27,7 +27,7 @@ public: typedef ippl::Field<double, dim> field_type; FieldTest() - : nPoints(100) + : nPoints(8) { setup(); } @@ -36,9 +36,9 @@ public: ippl::Index I(nPoints); ippl::NDIndex<dim> owned(I, I, I); - ippl::e_dim_tag allParallel[dim]; // Specifies SERIAL, PARALLEL dims + ippl::e_dim_tag allParallel[dim]; for (unsigned int d = 0; d < dim; d++) - allParallel[d] = ippl::SERIAL; + allParallel[d] = ippl::PARALLEL; ippl::FieldLayout<dim> layout(owned, allParallel); @@ -56,7 +56,7 @@ public: -TEST_F(FieldTest, FieldSum) { +TEST_F(FieldTest, Sum) { double val = 1.0; *field = val; @@ -67,10 +67,60 @@ TEST_F(FieldTest, FieldSum) { } +TEST_F(FieldTest, Norm1) { + double val = -1.5; + + *field = val; + + double norm1 = ippl::norm1(*field); + + ASSERT_DOUBLE_EQ(-val * std::pow(nPoints, dim), norm1); +} + + +TEST_F(FieldTest, Norm2) { + double val = 1.5; + + *field = val; + + double norm2 = ippl::norm2(*field); + + ASSERT_DOUBLE_EQ(std::sqrt(val * val * std::pow(nPoints, dim)), norm2); +} + +TEST_F(FieldTest, NormInf) { + + const ippl::NDIndex<dim> lDom = field->getLayout().getLocalNDIndex(); + + auto view = field->getView(); + auto policy = field->getRangePolicy(0); + + Kokkos::parallel_for("Assign field", + policy, + KOKKOS_LAMBDA(const int i, + const int j, + const int k) + { + const size_t ig = i + lDom[0].first(); + const size_t jg = j + lDom[1].first(); + const size_t kg = k + lDom[2].first(); + + view(i, j, k) = -1.0 + (ig + jg + kg); + }); + + + double normInf = ippl::normInf(*field); + + double val = -1.0 + 3 * nPoints; + + ASSERT_DOUBLE_EQ(val, normInf); +} + + int main(int argc, char *argv[]) { Ippl ippl(argc,argv); ::testing::InitGoogleTest(&argc, argv); return RUN_ALL_TESTS(); -} \ No newline at end of file +}