diff --git a/src/core/data/particles/particle.hpp b/src/core/data/particles/particle.hpp index af1c66f07..badab19ab 100644 --- a/src/core/data/particles/particle.hpp +++ b/src/core/data/particles/particle.hpp @@ -61,26 +61,17 @@ struct Particle std::array delta = ConstArray(); std::array v = ConstArray(); - double Ex = 0, Ey = 0, Ez = 0; - double Bx = 0, By = 0, Bz = 0; - NO_DISCARD bool operator==(Particle const& that) const { return (this->weight == that.weight) && // (this->charge == that.charge) && // (this->iCell == that.iCell) && // (this->delta == that.delta) && // - (this->v == that.v) && // - (this->Ex == that.Ex) && // - (this->Ey == that.Ey) && // - (this->Ez == that.Ez) && // - (this->Bx == that.Bx) && // - (this->By == that.By) && // - (this->Bz == that.Bz); + (this->v == that.v); } template - friend std::ostream& operator<<(std::ostream& out, const Particle& particle); + friend std::ostream& operator<<(std::ostream& out, Particle const& particle); }; template @@ -102,8 +93,6 @@ std::ostream& operator<<(std::ostream& out, Particle const& particle) out << v << ","; } out << "), charge : " << particle.charge << ", weight : " << particle.weight; - out << ", Exyz : " << particle.Ex << "," << particle.Ey << "," << particle.Ez; - out << ", Bxyz : " << particle.Bx << "," << particle.By << "," << particle.Bz; out << '\n'; return out; } @@ -132,9 +121,9 @@ inline constexpr auto is_phare_particle_type template typename ParticleA, template typename ParticleB> -NO_DISCARD typename std::enable_if_t> - and is_phare_particle_type>, - bool> +NO_DISCARD typename std::enable_if_t< + is_phare_particle_type> and is_phare_particle_type>, + bool> operator==(ParticleA const& particleA, ParticleB const& particleB) { return particleA.weight == particleB.weight and // diff --git a/src/core/data/particles/particle_array.hpp b/src/core/data/particles/particle_array.hpp index 0c4e42384..75ba5f396 100644 --- a/src/core/data/particles/particle_array.hpp +++ b/src/core/data/particles/particle_array.hpp @@ -231,6 +231,8 @@ class ParticleArray NO_DISCARD auto& vector() { return particles_; } NO_DISCARD auto& vector() const { return particles_; } + auto& box() const { return box_; } + private: Vector particles_; box_t box_; diff --git a/src/core/numerics/interpolator/interpolator.hpp b/src/core/numerics/interpolator/interpolator.hpp index c3a014ec2..8e9e15f17 100644 --- a/src/core/numerics/interpolator/interpolator.hpp +++ b/src/core/numerics/interpolator/interpolator.hpp @@ -9,12 +9,16 @@ #include "core/data/grid/gridlayout.hpp" #include "core/data/vecfield/vecfield_component.hpp" #include "core/utilities/point/point.hpp" -#include "core/def.hpp" +#include "core/def.hpp" #include "core/logger.hpp" + + namespace PHARE::core { + + //! return the size of the index and weights arrays for //! interpolation at a given order // Number of points where the interpolator of order interpOrder @@ -457,7 +461,8 @@ class Interpolator : private Weighter public: auto static constexpr interp_order = interpOrder; auto static constexpr dimension = dim; - /**\brief interpolate electromagnetic fields on all particles in the range + + /**\brief interpolate electromagnetic fields on a particle and return particles EB * * For each particle : * - The function first calculates the startIndex and weights for interpolation at @@ -465,18 +470,13 @@ class Interpolator : private Weighter * - then it uses Interpol<> to calculate the interpolation of E and B components * onto the particle. */ - template - inline void operator()(ParticleRange&& particleRange, Electromag const& Em, - GridLayout const& layout) + template + inline auto operator()(Particle_t& currPart, Electromag const& Em, GridLayout const& layout) { PHARE_LOG_SCOPE("Interpolator::operator()"); - auto begin = particleRange.begin(); - auto end = particleRange.end(); - - using Scalar = HybridQuantity::Scalar; - auto const& [Ex, Ey, Ez] = Em.E(); - auto const& [Bx, By, Bz] = Em.B(); + using E_B_tuple = std::tuple, std::array>; + using Scalar = HybridQuantity::Scalar; // for each particle, first calculate the startIndex and weights for dual and // primal quantities. then, knowing the centering (primal or dual) of each @@ -485,33 +485,36 @@ class Interpolator : private Weighter // calculated twice, and not for each E,B component. PHARE_LOG_START("MeshToParticle::operator()"); - for (auto currPart = begin; currPart != end; ++currPart) - { - auto& iCell = currPart->iCell; - auto& delta = currPart->delta; - indexAndWeights_(layout, iCell, delta); - indexAndWeights_(layout, iCell, delta); - - auto indexWeights = std::forward_as_tuple(dual_startIndex_, dual_weights_, - primal_startIndex_, primal_weights_); - - currPart->Ex - = meshToParticle_.template operator()(Ex, indexWeights); - currPart->Ey - = meshToParticle_.template operator()(Ey, indexWeights); - currPart->Ez - = meshToParticle_.template operator()(Ez, indexWeights); - currPart->Bx - = meshToParticle_.template operator()(Bx, indexWeights); - currPart->By - = meshToParticle_.template operator()(By, indexWeights); - currPart->Bz - = meshToParticle_.template operator()(Bz, indexWeights); - } + + auto& iCell = currPart.iCell; + auto& delta = currPart.delta; + indexAndWeights_(layout, iCell, delta); + indexAndWeights_(layout, iCell, delta); + + auto indexWeights = std::forward_as_tuple(dual_startIndex_, dual_weights_, + primal_startIndex_, primal_weights_); + + E_B_tuple particle_EB; + auto& [pE, pB] = particle_EB; + auto& [pEx, pEy, pEz] = pE; + auto& [pBx, pBy, pBz] = pB; + + auto const& [Ex, Ey, Ez] = Em.E(); + auto const& [Bx, By, Bz] = Em.B(); + + pEx = meshToParticle_.template operator()(Ex, indexWeights); + pEy = meshToParticle_.template operator()(Ey, indexWeights); + pEz = meshToParticle_.template operator()(Ez, indexWeights); + pBx = meshToParticle_.template operator()(Bx, indexWeights); + pBy = meshToParticle_.template operator()(By, indexWeights); + pBz = meshToParticle_.template operator()(Bz, indexWeights); + PHARE_LOG_STOP("MeshToParticle::operator()"); + return particle_EB; } + /**\brief interpolate electromagnetic fields on all particles in the range * * For each particle : @@ -521,7 +524,7 @@ class Interpolator : private Weighter * onto the particle. */ template - inline void operator()(ParticleRange&& particleRange, Field& density, VecField& flux, + inline void operator()(ParticleRange& particleRange, Field& density, VecField& flux, GridLayout const& layout, double coef = 1.) { auto begin = particleRange.begin(); @@ -548,6 +551,12 @@ class Interpolator : private Weighter } PHARE_LOG_STOP("ParticleToMesh::operator()"); } + template + inline void operator()(ParticleRange&& range, Field& density, VecField& flux, + GridLayout const& layout, double coef = 1.) + { + (*this)(range, density, flux, layout, coef); + } /** diff --git a/src/core/numerics/pusher/boris.hpp b/src/core/numerics/pusher/boris.hpp index 7faf26194..593516f39 100644 --- a/src/core/numerics/pusher/boris.hpp +++ b/src/core/numerics/pusher/boris.hpp @@ -15,6 +15,8 @@ namespace PHARE::core { + + template class BorisPusher @@ -85,29 +87,32 @@ class BorisPusher // rangeIn : t=n, rangeOut : t=n+1/2 // Do not partition on this step - this is to keep all domain and ghost // particles consistent. see: https://github.com/PHAREHUB/PHARE/issues/571 - pushStep_(rangeIn, rangeOut, PushStep::PrePush); + prePushStep_(rangeIn, rangeOut); rangeOut = firstSelector(rangeOut); - // get electromagnetic fields interpolated on the particles of rangeOut - // stop at newEnd. - interpolator(rangeOut, emFields, layout); - // get the particle velocity from t=n to t=n+1 - accelerate_(rangeOut, rangeOut, mass); + double const dto2m = 0.5 * dt_ / mass; + for (auto idx = rangeOut.ibegin(); idx < rangeOut.iend(); ++idx) + { + auto& currPart = rangeOut.array()[idx]; - // now advance the particles from t=n+1/2 to t=n+1 using v_{n+1} just calculated - // and get a pointer to the first leaving particle - pushStep_(rangeOut, rangeOut, PushStep::PostPush); + // get electromagnetic fields interpolated on the particles of rangeOut stop at newEnd. + // get the particle velocity from t=n to t=n+1 + accelerate_(currPart, interpolator(currPart, emFields, layout), dto2m); + + // now advance the particles from t=n+1/2 to t=n+1 using v_{n+1} just calculated + // and get a pointer to the first leaving particle + postPushStep_(rangeOut, idx); + } return secondSelector(rangeOut); } - /** see Pusher::move() documentation*/ - virtual void setMeshAndTimeStep(std::array ms, double ts) override + void setMeshAndTimeStep(std::array ms, double ts) override { std::transform(std::begin(ms), std::end(ms), std::begin(halfDtOverDl_), [ts](double& x) { return 0.5 * ts / x; }); @@ -117,8 +122,6 @@ class BorisPusher private: - enum class PushStep { PrePush, PostPush }; - /** move the particle partIn of half a time step and store it in partOut */ template @@ -148,7 +151,7 @@ class BorisPusher * @return the function returns and iterator on the first leaving particle, as * detected by the ParticleSelector */ - void pushStep_(ParticleRange const& rangeIn, ParticleRange& rangeOut, PushStep step) + void prePushStep_(ParticleRange const& rangeIn, ParticleRange& rangeOut) { auto& inParticles = rangeIn.array(); auto& outParticles = rangeOut.array(); @@ -164,88 +167,89 @@ class BorisPusher // not strictly speaking this function's business // to take advantage that we're already looping // over rangeIn particles. - if (step == PushStep::PrePush) - { - outParticles[outIdx].charge = inParticles[inIdx].charge; - outParticles[outIdx].weight = inParticles[inIdx].weight; - outParticles[outIdx].v = inParticles[inIdx].v; - } + + outParticles[outIdx].charge = inParticles[inIdx].charge; + outParticles[outIdx].weight = inParticles[inIdx].weight; + outParticles[outIdx].v = inParticles[inIdx].v; + auto newCell = advancePosition_(inParticles[inIdx], outParticles[outIdx]); if (newCell != inParticles[inIdx].iCell) outParticles.change_icell(newCell, outIdx); } } + void postPushStep_(ParticleRange& range, std::size_t idx) + { + auto& particles = range.array(); + auto newCell = advancePosition_(particles[idx], particles[idx]); + if (newCell != particles[idx].iCell) + particles.change_icell(newCell, idx); + } /** Accelerate the particles in rangeIn and put the new velocity in rangeOut */ - void accelerate_(ParticleRange rangeIn, ParticleRange rangeOut, double mass) + template + void accelerate_(Particle_t& part, ParticleEB const& particleEB, double const& dto2m) { - double dto2m = 0.5 * dt_ / mass; + auto& [pE, pB] = particleEB; + auto& [pEx, pEy, pEz] = pE; + auto& [pBx, pBy, pBz] = pB; - auto& inParticles = rangeIn.array(); - auto& outParticles = rangeOut.array(); - for (auto inIdx = rangeIn.ibegin(), outIdx = rangeOut.ibegin(); inIdx < rangeIn.iend(); - ++inIdx, ++outIdx) - { - auto& inPart = inParticles[inIdx]; - auto& outPart = outParticles[inIdx]; - double coef1 = inPart.charge * dto2m; + double const coef1 = part.charge * dto2m; - // We now apply the 3 steps of the BORIS PUSHER + // We now apply the 3 steps of the BORIS PUSHER - // 1st half push of the electric field - double velx1 = inPart.v[0] + coef1 * inPart.Ex; - double vely1 = inPart.v[1] + coef1 * inPart.Ey; - double velz1 = inPart.v[2] + coef1 * inPart.Ez; + // 1st half push of the electric field + double velx1 = part.v[0] + coef1 * pEx; + double vely1 = part.v[1] + coef1 * pEy; + double velz1 = part.v[2] + coef1 * pEz; - // preparing variables for magnetic rotation - double const rx = coef1 * inPart.Bx; - double const ry = coef1 * inPart.By; - double const rz = coef1 * inPart.Bz; + // preparing variables for magnetic rotation + double const rx = coef1 * pBx; + double const ry = coef1 * pBy; + double const rz = coef1 * pBz; - double const rx2 = rx * rx; - double const ry2 = ry * ry; - double const rz2 = rz * rz; - double const rxry = rx * ry; - double const rxrz = rx * rz; - double const ryrz = ry * rz; + double const rx2 = rx * rx; + double const ry2 = ry * ry; + double const rz2 = rz * rz; + double const rxry = rx * ry; + double const rxrz = rx * rz; + double const ryrz = ry * rz; - double const invDet = 1. / (1. + rx2 + ry2 + rz2); + double const invDet = 1. / (1. + rx2 + ry2 + rz2); - // preparing rotation matrix due to the magnetic field - // m = invDet*(I + r*r - r x I) - I where x denotes the cross product - double const mxx = 1. + rx2 - ry2 - rz2; - double const mxy = 2. * (rxry + rz); - double const mxz = 2. * (rxrz - ry); + // preparing rotation matrix due to the magnetic field + // m = invDet*(I + r*r - r x I) - I where x denotes the cross product + double const mxx = 1. + rx2 - ry2 - rz2; + double const mxy = 2. * (rxry + rz); + double const mxz = 2. * (rxrz - ry); - double const myx = 2. * (rxry - rz); - double const myy = 1. + ry2 - rx2 - rz2; - double const myz = 2. * (ryrz + rx); + double const myx = 2. * (rxry - rz); + double const myy = 1. + ry2 - rx2 - rz2; + double const myz = 2. * (ryrz + rx); - double const mzx = 2. * (rxrz + ry); - double const mzy = 2. * (ryrz - rx); - double const mzz = 1. + rz2 - rx2 - ry2; + double const mzx = 2. * (rxrz + ry); + double const mzy = 2. * (ryrz - rx); + double const mzz = 1. + rz2 - rx2 - ry2; - // magnetic rotation - double const velx2 = (mxx * velx1 + mxy * vely1 + mxz * velz1) * invDet; - double const vely2 = (myx * velx1 + myy * vely1 + myz * velz1) * invDet; - double const velz2 = (mzx * velx1 + mzy * vely1 + mzz * velz1) * invDet; + // magnetic rotation + double const velx2 = (mxx * velx1 + mxy * vely1 + mxz * velz1) * invDet; + double const vely2 = (myx * velx1 + myy * vely1 + myz * velz1) * invDet; + double const velz2 = (mzx * velx1 + mzy * vely1 + mzz * velz1) * invDet; - // 2nd half push of the electric field - velx1 = velx2 + coef1 * inPart.Ex; - vely1 = vely2 + coef1 * inPart.Ey; - velz1 = velz2 + coef1 * inPart.Ez; + // 2nd half push of the electric field + velx1 = velx2 + coef1 * pEx; + vely1 = vely2 + coef1 * pEy; + velz1 = velz2 + coef1 * pEz; - // Update particle velocity - outPart.v[0] = velx1; - outPart.v[1] = vely1; - outPart.v[2] = velz1; - } + // Update particle velocity + part.v[0] = velx1; + part.v[1] = vely1; + part.v[2] = velz1; } diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index c6579299f..fd3e1fe09 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -25,6 +25,9 @@ class box_iterator; template struct Box { + static const size_t dimension = dim; + + Point lower; Point upper; @@ -257,8 +260,8 @@ NO_DISCARD bool isIn(Point const& point, } -template -Box grow(Box const& box, Type const& size) +template +Box grow(Box const& box, OType const& size) { auto copy{box}; copy.grow(size); diff --git a/src/core/utilities/point/point.hpp b/src/core/utilities/point/point.hpp index 4ad05e902..754b0a12e 100644 --- a/src/core/utilities/point/point.hpp +++ b/src/core/utilities/point/point.hpp @@ -85,7 +85,7 @@ namespace core NO_DISCARD bool operator!=(Point const& other) const { return !(*this == other); } - template + template NO_DISCARD auto toArray() const { std::array destArray; diff --git a/tests/amr/data/particles/copy/test_particledata_copyNd.cpp b/tests/amr/data/particles/copy/test_particledata_copyNd.cpp index 27132c515..f7f5e823a 100644 --- a/tests/amr/data/particles/copy/test_particledata_copyNd.cpp +++ b/tests/amr/data/particles/copy/test_particledata_copyNd.cpp @@ -121,12 +121,6 @@ TYPED_TEST(AParticlesDataND, PreservesAllParticleAttributesAfterCopy) Pointwise(DoubleEq(), this->particle.delta)); EXPECT_THAT(this->destData.domainParticles[0].weight, DoubleEq(this->particle.weight)); EXPECT_THAT(this->destData.domainParticles[0].charge, DoubleEq(this->particle.charge)); - EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Ex, this->particle.Ex); - EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Ey, this->particle.Ey); - EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Ez, this->particle.Ez); - EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Bx, this->particle.Bx); - EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].By, this->particle.By); - EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Bz, this->particle.Bz); // particle is in the domain of the source patchdata // and in last ghost of the destination patchdata @@ -142,12 +136,6 @@ TYPED_TEST(AParticlesDataND, PreservesAllParticleAttributesAfterCopy) Pointwise(DoubleEq(), this->particle.delta)); EXPECT_THAT(this->destData.patchGhostParticles[0].weight, DoubleEq(this->particle.weight)); EXPECT_THAT(this->destData.patchGhostParticles[0].charge, DoubleEq(this->particle.charge)); - EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Ex, this->particle.Ex); - EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Ey, this->particle.Ey); - EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Ez, this->particle.Ez); - EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Bx, this->particle.Bx); - EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].By, this->particle.By); - EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Bz, this->particle.Bz); } diff --git a/tests/amr/data/particles/copy_overlap/test_particledata_copy_periodicNd.cpp b/tests/amr/data/particles/copy_overlap/test_particledata_copy_periodicNd.cpp index 09d19ba1d..31d53cf05 100644 --- a/tests/amr/data/particles/copy_overlap/test_particledata_copy_periodicNd.cpp +++ b/tests/amr/data/particles/copy_overlap/test_particledata_copy_periodicNd.cpp @@ -139,12 +139,6 @@ TYPED_TEST(twoParticlesDataNDTouchingPeriodicBorders, preserveParticleAttributes EXPECT_THAT(this->destPdat.patchGhostParticles[0].delta, Eq(this->particle.delta)); EXPECT_THAT(this->destPdat.patchGhostParticles[0].weight, Eq(this->particle.weight)); EXPECT_THAT(this->destPdat.patchGhostParticles[0].charge, Eq(this->particle.charge)); - EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Ex, this->particle.Ex); - EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Ey, this->particle.Ey); - EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Ez, this->particle.Ez); - EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Bx, this->particle.Bx); - EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].By, this->particle.By); - EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Bz, this->particle.Bz); } diff --git a/tests/amr/data/particles/stream_pack/test_main.cpp b/tests/amr/data/particles/stream_pack/test_main.cpp index a2bfe65fd..7905eb579 100644 --- a/tests/amr/data/particles/stream_pack/test_main.cpp +++ b/tests/amr/data/particles/stream_pack/test_main.cpp @@ -228,12 +228,6 @@ TYPED_TEST(StreamPackTest, EXPECT_THAT(destData.domainParticles[0].delta, Eq(particle.delta)); EXPECT_THAT(destData.domainParticles[0].weight, Eq(particle.weight)); EXPECT_THAT(destData.domainParticles[0].charge, Eq(particle.charge)); - EXPECT_DOUBLE_EQ(destData.domainParticles[0].Ex, particle.Ex); - EXPECT_DOUBLE_EQ(destData.domainParticles[0].Ey, particle.Ey); - EXPECT_DOUBLE_EQ(destData.domainParticles[0].Ez, particle.Ez); - EXPECT_DOUBLE_EQ(destData.domainParticles[0].Bx, particle.Bx); - EXPECT_DOUBLE_EQ(destData.domainParticles[0].By, particle.By); - EXPECT_DOUBLE_EQ(destData.domainParticles[0].Bz, particle.Bz); } @@ -271,12 +265,6 @@ TYPED_TEST(StreamPackTest, EXPECT_THAT(destData.patchGhostParticles[0].delta, Eq(particle.delta)); EXPECT_THAT(destData.patchGhostParticles[0].weight, Eq(particle.weight)); EXPECT_THAT(destData.patchGhostParticles[0].charge, Eq(particle.charge)); - EXPECT_DOUBLE_EQ(destData.patchGhostParticles[0].Ex, particle.Ex); - EXPECT_DOUBLE_EQ(destData.patchGhostParticles[0].Ey, particle.Ey); - EXPECT_DOUBLE_EQ(destData.patchGhostParticles[0].Ez, particle.Ez); - EXPECT_DOUBLE_EQ(destData.patchGhostParticles[0].Bx, particle.Bx); - EXPECT_DOUBLE_EQ(destData.patchGhostParticles[0].By, particle.By); - EXPECT_DOUBLE_EQ(destData.patchGhostParticles[0].Bz, particle.Bz); } diff --git a/tests/core/data/gridlayout/gridlayout_test.hpp b/tests/core/data/gridlayout/gridlayout_test.hpp index 7f1897fdc..9ecd8c280 100644 --- a/tests/core/data/gridlayout/gridlayout_test.hpp +++ b/tests/core/data/gridlayout/gridlayout_test.hpp @@ -5,6 +5,7 @@ #include "core/data/grid/gridlayoutimplyee.hpp" #include "core/utilities/types.hpp" +#include "test_gridlayout.hpp" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -23,22 +24,5 @@ class GridLayoutTest : public ::testing::TestWithParam> Param param; }; -template -class TestGridLayout : public GridLayout -{ // to expose a default constructor -public: - auto static constexpr dim = GridLayout::dimension; - - TestGridLayout() = default; - - TestGridLayout(std::uint32_t cells) - : GridLayout{PHARE::core::ConstArray(1.0 / cells), - PHARE::core::ConstArray(cells), - PHARE::core::Point{PHARE::core::ConstArray(0)}} - { - } - - auto static make(std::uint32_t cells) { return TestGridLayout{cells}; } -}; #endif diff --git a/tests/core/data/gridlayout/test_gridlayout.hpp b/tests/core/data/gridlayout/test_gridlayout.hpp new file mode 100644 index 000000000..d31f630a6 --- /dev/null +++ b/tests/core/data/gridlayout/test_gridlayout.hpp @@ -0,0 +1,27 @@ +#ifndef TESTS_CORE_DATA_GRIDLAYOUT_TEST_GRIDLAYOUT_HPP +#define TESTS_CORE_DATA_GRIDLAYOUT_TEST_GRIDLAYOUT_HPP + +#include "core/data/grid/gridlayout.hpp" +#include "core/data/grid/gridlayoutimplyee.hpp" +#include "core/utilities/types.hpp" + + +template +class TestGridLayout : public GridLayout +{ // to expose a default constructor +public: + auto static constexpr dim = GridLayout::dimension; + + TestGridLayout() = default; + + TestGridLayout(std::uint32_t cells) + : GridLayout{PHARE::core::ConstArray(1.0 / cells), + PHARE::core::ConstArray(cells), + PHARE::core::Point{PHARE::core::ConstArray(0)}} + { + } + + auto static make(std::uint32_t cells) { return TestGridLayout{cells}; } +}; + +#endif /*TESTS_CORE_DATA_GRIDLAYOUT_TEST_GRIDLAYOUT_HPP*/ diff --git a/tests/core/data/particles/test_main.cpp b/tests/core/data/particles/test_main.cpp index a4771ec47..f3ad991c1 100644 --- a/tests/core/data/particles/test_main.cpp +++ b/tests/core/data/particles/test_main.cpp @@ -38,17 +38,6 @@ TEST_F(AParticle, ParticleChargeIsInitiliazedOK) EXPECT_DOUBLE_EQ(1., part.charge); } -TEST_F(AParticle, ParticleFieldsAreInitializedToZero) -{ - EXPECT_DOUBLE_EQ(0.0, part.Ex); - EXPECT_DOUBLE_EQ(0.0, part.Ey); - EXPECT_DOUBLE_EQ(0.0, part.Ez); - - EXPECT_DOUBLE_EQ(0.0, part.Bx); - EXPECT_DOUBLE_EQ(0.0, part.By); - EXPECT_DOUBLE_EQ(0.0, part.Bz); -} - TEST_F(AParticle, ParticleVelocityIsInitializedOk) { diff --git a/tests/core/numerics/interpolator/test_main.cpp b/tests/core/numerics/interpolator/test_main.cpp index 9d3be8768..dc0f01635 100644 --- a/tests/core/numerics/interpolator/test_main.cpp +++ b/tests/core/numerics/interpolator/test_main.cpp @@ -298,32 +298,18 @@ TYPED_TEST(A1DInterpolator, canComputeAllEMfieldsAtParticle) this->em.B.setBuffer("EM_B_y", &this->by1d_); this->em.B.setBuffer("EM_B_z", &this->bz1d_); - this->interp(makeIndexRange(this->particles), this->em, this->layout); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ex - this->ex0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ey - this->ey0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ez - this->ez0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Bx - this->bx0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.By - this->by0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Bz - this->bz0) < 1e-8; })); - + for (auto const& part : this->particles) + { + auto const [E, B] = this->interp(part, this->em, this->layout); + auto const& [Ex, Ey, Ez] = E; + auto const& [Bx, By, Bz] = B; + EXPECT_NEAR(Ex, this->ex0, 1e-8); + EXPECT_NEAR(Ey, this->ey0, 1e-8); + EXPECT_NEAR(Ez, this->ez0, 1e-8); + EXPECT_NEAR(Bx, this->bx0, 1e-8); + EXPECT_NEAR(By, this->by0, 1e-8); + EXPECT_NEAR(Bz, this->bz0, 1e-8); + } this->em.E.setBuffer("EM_E_x", nullptr); this->em.E.setBuffer("EM_E_y", nullptr); @@ -421,32 +407,18 @@ TYPED_TEST(A2DInterpolator, canComputeAllEMfieldsAtParticle) this->em.B.setBuffer("EM_B_y", &this->by_); this->em.B.setBuffer("EM_B_z", &this->bz_); - this->interp(makeIndexRange(this->particles), this->em, this->layout); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ex - this->ex0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ey - this->ey0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ez - this->ez0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Bx - this->bx0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.By - this->by0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Bz - this->bz0) < 1e-8; })); - + for (auto const& part : this->particles) + { + auto const [E, B] = this->interp(part, this->em, this->layout); + auto const& [Ex, Ey, Ez] = E; + auto const& [Bx, By, Bz] = B; + EXPECT_NEAR(Ex, this->ex0, 1e-8); + EXPECT_NEAR(Ey, this->ey0, 1e-8); + EXPECT_NEAR(Ez, this->ez0, 1e-8); + EXPECT_NEAR(Bx, this->bx0, 1e-8); + EXPECT_NEAR(By, this->by0, 1e-8); + EXPECT_NEAR(Bz, this->bz0, 1e-8); + } this->em.E.setBuffer("EM_E_x", nullptr); this->em.E.setBuffer("EM_E_y", nullptr); @@ -549,31 +521,19 @@ TYPED_TEST(A3DInterpolator, canComputeAllEMfieldsAtParticle) this->em.B.setBuffer("EM_B_y", &this->by_); this->em.B.setBuffer("EM_B_z", &this->bz_); - this->interp(makeIndexRange(this->particles), this->em, this->layout); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ex - this->ex0) < 1e-8; })); - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ey - this->ey0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Ez - this->ez0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Bx - this->bx0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.By - this->by0) < 1e-8; })); - - EXPECT_TRUE( - std::all_of(std::begin(this->particles), std::end(this->particles), - [this](auto const& part) { return std::abs(part.Bz - this->bz0) < 1e-8; })); + for (auto const& part : this->particles) + { + auto const [E, B] = this->interp(part, this->em, this->layout); + auto const& [Ex, Ey, Ez] = E; + auto const& [Bx, By, Bz] = B; + EXPECT_NEAR(Ex, this->ex0, 1e-8); + EXPECT_NEAR(Ey, this->ey0, 1e-8); + EXPECT_NEAR(Ez, this->ez0, 1e-8); + EXPECT_NEAR(Bx, this->bx0, 1e-8); + EXPECT_NEAR(By, this->by0, 1e-8); + EXPECT_NEAR(Bz, this->bz0, 1e-8); + } this->em.E.setBuffer("EM_E_x", nullptr); diff --git a/tests/core/numerics/pusher/test_pusher.cpp b/tests/core/numerics/pusher/test_pusher.cpp index 729848977..8bad3a9d1 100644 --- a/tests/core/numerics/pusher/test_pusher.cpp +++ b/tests/core/numerics/pusher/test_pusher.cpp @@ -20,6 +20,7 @@ using namespace PHARE::core; + struct Trajectory { std::vector x, y, z; @@ -70,19 +71,25 @@ Trajectory readExpectedTrajectory() // to test the pusher. class Interpolator { + using E_B_tuple = std::tuple, std::array>; + public: - template - void operator()(ParticleRange particles, Electromag const&, GridLayout&) + template + auto operator()(Particle_t& particle, Electromag const&, GridLayout&) { - for (auto currPart = std::begin(particles); currPart != std::end(particles); ++currPart) - { - currPart->Ex = 0.01; - currPart->Ey = -0.05; - currPart->Ez = 0.05; - currPart->Bx = 1.; - currPart->By = 1.; - currPart->Bz = 1.; - } + E_B_tuple eb_interop; + auto& [pE, pB] = eb_interop; + auto& [pEx, pEy, pEz] = pE; + auto& [pBx, pBy, pBz] = pB; + + pEx = 0.01; + pEy = -0.05; + pEz = 0.05; + pBx = 1.; + pBy = 1.; + pBz = 1.; + + return eb_interop; } }; @@ -132,6 +139,7 @@ class APusher : public ::testing::Test , tstart{0} , tend{10} , nt{static_cast((tend - tstart) / dt + 1)} + { particlesIn.emplace_back( Particle{1., 1., ConstArray(5), ConstArray(0.), {0., 10., 0.}}); @@ -180,9 +188,7 @@ TEST_F(APusher3D, trajectoryIsOk) actual[1][i] = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * dxyz[1]; actual[2][i] = (particlesOut[0].iCell[2] + particlesOut[0].delta[2]) * dxyz[2]; - pusher->move( - rangeIn, rangeOut, em, mass, interpolator, layout, - [](decltype(rangeIn)& rge) { return rge; }, selector); + pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); } @@ -206,9 +212,7 @@ TEST_F(APusher2D, trajectoryIsOk) actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0]; actual[1][i] = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * dxyz[1]; - pusher->move( - rangeIn, rangeOut, em, mass, interpolator, layout, [](auto& rge) { return rge; }, - selector); + pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); } diff --git a/tools/bench/core/bench.hpp b/tools/bench/core/bench.hpp index 932b1d650..df84a33bb 100644 --- a/tools/bench/core/bench.hpp +++ b/tools/bench/core/bench.hpp @@ -2,7 +2,7 @@ #define PHARE_BENCH_CORE_BENCH_H #include "phare_core.hpp" -#include "benchmark/benchmark.hpp" +#include "benchmark/benchmark.h" namespace PHARE::core::bench @@ -166,4 +166,69 @@ class Electromag : public PHARE::core::Electromag +class LocalisedCellFlattener +{ +public: + static constexpr std::size_t dim = Box_t::dimension; + + LocalisedCellFlattener(Box_t const& b) + : box{b} + , shape{box.shape().toArray()} + { + } + + RValue operator()(std::array icell) const + { + for (std::size_t i = 0; i < dim; ++i) + icell[i] -= box.lower[i]; + + if constexpr (dim == 2) + return icell[1] + icell[0] * shape[1]; + if constexpr (dim == 3) + return icell[2] + icell[1] * shape[2] + icell[0] * shape[1] * shape[2]; + return icell[0]; + } + template + RValue operator()(Particle const& particle) const + { + return (*this)(particle.iCell); + } + + + Box_t const box; + +private: + std::array const shape = box.shape().toArray(); +}; +} // namespace PHARE::core + + +namespace std +{ + +template +auto& sort(PHARE::core::ParticleArray& particles) +{ + using box_t = typename PHARE::core::ParticleArray::box_t; + PHARE::core::LocalisedCellFlattener cell_flattener{grow(particles.box(), 1)}; + std::sort(particles.vector().begin(), particles.vector().end(), + [&](auto const& a, auto const& b) { + return cell_flattener(a.iCell) < cell_flattener(b.iCell); + }); + return particles; +} + +template +auto sort(PHARE::core::ParticleArray&& particles) +{ + return sort(particles); +} + + +} // namespace std + + #endif /*PHARE_BENCH_CORE_BENCH_H*/ diff --git a/tools/bench/core/numerics/pusher/pusher.cpp b/tools/bench/core/numerics/pusher/pusher.cpp index 2bd49b791..55cef42df 100644 --- a/tools/bench/core/numerics/pusher/pusher.cpp +++ b/tools/bench/core/numerics/pusher/pusher.cpp @@ -1,92 +1,60 @@ +#include "tools/bench/core/bench.hpp" +#include "tests/core/data/gridlayout/test_gridlayout.hpp" -#include "benchmark/benchmark.hpp" - -#include "phare_core.hpp" #include "core/numerics/pusher/boris.hpp" #include "core/numerics/ion_updater/ion_updater.hpp" -template -using Field = PHARE::core::Field, - typename PHARE::core::HybridQuantity::Scalar>; -template -using VecField - = PHARE::core::VecField, typename PHARE::core::HybridQuantity>; - -template -PHARE::core::Particle particle() -{ - return {// - /*.weight = */ 0, - /*.charge = */ 1, - /*.iCell = */ PHARE::core::ConstArray(35), - /*.delta = */ PHARE::core::ConstArray(.01), - /*.v = */ {{0, 10., 0}}}; -} - -template -Field field(std::string key, Quantity type, GridLayout const& layout) -{ - Field feeld{key, type, layout.allocSize(type)}; - std::fill(feeld.begin(), feeld.end(), 1); - return feeld; -} +using namespace PHARE; template void push(benchmark::State& state) { - constexpr std::uint32_t cells = 65; - constexpr std::uint32_t parts = 1e7; + constexpr std::uint32_t cells = 65; + constexpr std::uint32_t n_parts = 1e7; using PHARE_Types = PHARE::core::PHARE_Types; + using GridLayout_t = TestGridLayout; using Interpolator = PHARE::core::Interpolator; using BoundaryCondition = PHARE::core::BoundaryCondition; + using Electromag_t = core::bench::Electromag; using Ions_t = typename PHARE_Types::Ions_t; - using Electromag_t = typename PHARE_Types::Electromag_t; - using GridLayout_t = typename PHARE_Types::GridLayout_t; using ParticleArray = typename Ions_t::particle_array_type; - using PartIterator = typename ParticleArray::iterator; + using Particle_t = typename ParticleArray::value_type; + using ParticleRange = core::IndexRange; - using BorisPusher_t = PHARE::core::BorisPusher; - Interpolator interpolator; - ParticleArray domainParticles{parts, particle()}; - ParticleArray tmpDomain{domainParticles.size(), particle()}; - auto rangeIn = PHARE::core::makeRange(domainParticles); - auto rangeOut = PHARE::core::makeRange(tmpDomain); + GridLayout_t layout{cells}; + Electromag_t em{layout}; - auto meshSize = PHARE::core::ConstArray(1.0 / cells); - auto nCells = PHARE::core::ConstArray(cells); - auto origin = PHARE::core::Point{PHARE::core::ConstArray(0)}; - GridLayout_t layout{meshSize, nCells, origin}; + ParticleArray domainParticles{layout.AMRBox()}; + domainParticles.vector() = std::vector(n_parts, core::bench::particle()); + core::bench::disperse(domainParticles, 0, cells - 1, 13337); + // std::sort(domainParticles); - Field bx = field("Bx", PHARE::core::HybridQuantity::Scalar::Bx, layout); - Field by = field("By", PHARE::core::HybridQuantity::Scalar::By, layout); - Field bz = field("Bz", PHARE::core::HybridQuantity::Scalar::Bz, layout); + ParticleArray tmpDomain{layout.AMRBox()}; + tmpDomain.vector() = domainParticles.vector(); - Field ex = field("Ex", PHARE::core::HybridQuantity::Scalar::Ex, layout); - Field ey = field("Ey", PHARE::core::HybridQuantity::Scalar::Ey, layout); - Field ez = field("Ez", PHARE::core::HybridQuantity::Scalar::Ez, layout); - - PHARE::core::Electromag> emFields{std::string{"EM"}}; - emFields.B.setBuffer("EM_B_x", &bx); - emFields.B.setBuffer("EM_B_y", &by); - emFields.B.setBuffer("EM_B_z", &bz); - emFields.E.setBuffer("EM_E_x", &ex); - emFields.E.setBuffer("EM_E_y", &ey); - emFields.E.setBuffer("EM_E_z", &ez); + auto rangeIn = PHARE::core::makeIndexRange(domainParticles); + auto rangeOut = PHARE::core::makeIndexRange(tmpDomain); BorisPusher_t pusher; pusher.setMeshAndTimeStep(layout.meshSize(), .001); - while (state.KeepRunning()) + Interpolator interpolator; + auto const no_op = [](auto& particleRange) { return particleRange; }; + while (state.KeepRunningBatch(1)) { pusher.move( - /*ParticleRange const&*/ rangeIn, /*ParticleRange&*/ rangeOut, - /*Electromag const&*/ emFields, /*double mass*/ 1, /*Interpolator&*/ interpolator, - /*ParticleSelector const&*/ [](auto const& /*part*/) { return true; }, - /*GridLayout const&*/ layout); + /*ParticleRange const&*/ rangeIn, // + /*ParticleRange& */ rangeOut, // + /*Electromag const& */ em, // + /*double mass */ 1, + /*Interpolator&*/ interpolator, + /*GridLayout const&*/ layout, // + no_op, no_op); } } BENCHMARK_TEMPLATE(push, /*dim=*/1, /*interp=*/1)->Unit(benchmark::kMicrosecond); diff --git a/tools/bench/hi5/write_particles.cpp b/tools/bench/hi5/write_particles.cpp index 545438e73..0d4a6e798 100644 --- a/tools/bench/hi5/write_particles.cpp +++ b/tools/bench/hi5/write_particles.cpp @@ -1,5 +1,5 @@ -#include "benchmark/benchmark.hpp" +#include "benchmark/benchmark.h" #define PHARE_DIAG_DOUBLES 0 #include "diagnostic/detail/h5writer.hpp" #include "diagnostic/detail/h5_utils.hpp"