Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects
Commit b890e8a3 authored by vinciguerra_a's avatar vinciguerra_a
Browse files

Revert 2D/3D support

Revert "Merge branch '87-add-mdrangepolicy-generator-for-arbitrary-dimensions-2' into 'master'"

This reverts commit a752facc, reversing
changes made to 07468628.

The following commits are preserved:
- replace TestFieldNorm with unit test
- more rigorous normInf test
parent 37dbe4f8
No related branches found
No related tags found
1 merge request!91Revert 2D/3D support
......@@ -62,7 +62,6 @@ 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
......@@ -145,32 +144,16 @@ 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 (3D version)
* Assign an arbitrary BareField expression
* @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 == 3), bool> = true>
template <typename E, size_t N>
BareField<T, Dim>& operator=(const detail::Expression<E, N>& expr);
/*!
......@@ -199,44 +182,17 @@ 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;
#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);
T sum(int nghost = 0);
T max(int nghost = 0);
T min(int nghost = 0);
T prod(int nghost = 0);
DefineOperation(sum)
DefineOperation(max)
DefineOperation(min)
DefineOperation(prod)
private:
//! Number of ghost layers on each field boundary
......@@ -261,6 +217,5 @@ namespace ippl {
}
#include "Field/BareField.hpp"
#include "Field/BareFieldOperations.hpp"
#endif
......@@ -58,11 +58,13 @@ namespace ippl {
template <typename T, unsigned Dim>
void BareField<T, Dim>::setup() {
static_assert(Dim == 2 || Dim == 3, "Only 2D and 3D fields supported at the momment!");
static_assert(Dim == 3, "Only 3-dimensional fields supported at the momment!");
owned_m = layout_m->getLocalNDIndex();
if constexpr(Dim == 2) {
if constexpr(Dim == 1) {
this->resize(owned_m[0].length() + 2 * nghost_m);
} else 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) {
......@@ -96,67 +98,36 @@ 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) {
policy_type policy = getRangePolicy(nghost_m);
using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
Kokkos::parallel_for("BareField::operator=(T)",
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,
mdrange_type({0, 0, 0},
{dview_m.extent(0),
dview_m.extent(1),
dview_m.extent(2)
}),
KOKKOS_CLASS_LAMBDA(const size_t i,
const size_t j)
const size_t j,
const size_t k)
{
dview_m(i, j) = expr_(i, j);
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 == 3), bool>>
template <typename E, size_t N>
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);
using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
Kokkos::parallel_for("BareField::operator=(const Expression&)",
policy,
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}),
KOKKOS_CLASS_LAMBDA(const size_t i,
const size_t j,
const size_t k)
......@@ -174,51 +145,32 @@ namespace ippl {
}
#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; \
#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; \
}
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)
}
//
// 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));
}
}
......@@ -4,7 +4,6 @@ set (_SRCS
set (_HDRS
BareField.hpp
BareField.h
BareFieldOperations.hpp
BConds.h
BConds.hpp
BcTypes.h
......
......@@ -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);
}
......
//
// File FieldOperations
// Differential operators for fields
// Differential operators, norms, and a scalar product for fields
//
// Copyright (c) 2021 Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
......@@ -17,6 +17,88 @@
//
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
......
......@@ -62,23 +62,6 @@ 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
......
......@@ -36,9 +36,9 @@ public:
ippl::Index I(nPoints);
ippl::NDIndex<dim> owned(I, I, I);
ippl::e_dim_tag allParallel[dim];
ippl::e_dim_tag allParallel[dim]; // Specifies SERIAL, PARALLEL dims
for (unsigned int d = 0; d < dim; d++)
allParallel[d] = ippl::PARALLEL;
allParallel[d] = ippl::SERIAL;
ippl::FieldLayout<dim> layout(owned, allParallel);
......@@ -72,7 +72,7 @@ TEST_F(FieldTest, Norm1) {
*field = val;
double norm1 = ippl::norm1(*field);
double norm1 = ippl::norm(*field, 1);
ASSERT_DOUBLE_EQ(-val * std::pow(nPoints, dim), norm1);
}
......@@ -83,7 +83,7 @@ TEST_F(FieldTest, Norm2) {
*field = val;
double norm2 = ippl::norm2(*field);
double norm2 = ippl::norm(*field);
ASSERT_DOUBLE_EQ(std::sqrt(val * val * std::pow(nPoints, dim)), norm2);
}
......@@ -93,7 +93,11 @@ TEST_F(FieldTest, NormInf) {
const ippl::NDIndex<dim> lDom = field->getLayout().getLocalNDIndex();
auto view = field->getView();
auto policy = field->getRangePolicy(0);
const int shift = field->getNghost();
auto policy = Kokkos::MDRangePolicy<Kokkos::Rank<3>>(
{shift, shift, shift},
{view.extent(0) - shift, view.extent(1) - shift, view.extent(2) - shift}
);
Kokkos::parallel_for("Assign field",
policy,
......@@ -109,7 +113,7 @@ TEST_F(FieldTest, NormInf) {
});
double normInf = ippl::normInf(*field);
double normInf = ippl::norm(*field, 0);
double val = -1.0 + 3 * nPoints;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment