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/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..86ec08222 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 @@ -91,10 +93,14 @@ class BorisPusher // get electromagnetic fields interpolated on the particles of rangeOut // stop at newEnd. - interpolator(rangeOut, emFields, layout); + dto2m = 0.5 * dt_ / mass; + for (auto& currPart : rangeOut) + { + auto const& particle_EB = interpolator(currPart, emFields, layout); - // get the particle velocity from t=n to t=n+1 - accelerate_(rangeOut, rangeOut, mass); + // get the particle velocity from t=n to t=n+1 + accelerate_(currPart, particle_EB); + } // 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 @@ -148,7 +154,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 pushStep_(ParticleRange const& rangeIn, ParticleRange& rangeOut, PushStep const step) { auto& inParticles = rangeIn.array(); auto& outParticles = rangeOut.array(); @@ -180,72 +186,67 @@ class BorisPusher /** 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 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; } @@ -253,6 +254,7 @@ class BorisPusher std::array halfDtOverDl_; double dt_; + double dto2m; }; } // namespace PHARE::core diff --git a/tests/amr/data/particles/copy/test_particledata_copyNd.cpp b/tests/amr/data/particles/copy/test_particledata_copyNd.cpp index f64bc980d..d26df32c3 100644 --- a/tests/amr/data/particles/copy/test_particledata_copyNd.cpp +++ b/tests/amr/data/particles/copy/test_particledata_copyNd.cpp @@ -117,12 +117,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 @@ -138,12 +132,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 871de56fe..9a2cc79d3 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 @@ -135,12 +135,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 2531652d6..7e00e1e71 100644 --- a/tests/amr/data/particles/stream_pack/test_main.cpp +++ b/tests/amr/data/particles/stream_pack/test_main.cpp @@ -226,12 +226,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); } @@ -269,12 +263,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 b4d2d904f..9ecd8c280 100644 --- a/tests/core/data/gridlayout/gridlayout_test.hpp +++ b/tests/core/data/gridlayout/gridlayout_test.hpp @@ -24,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..85478e0a2 100644 --- a/tests/core/numerics/interpolator/test_main.cpp +++ b/tests/core/numerics/interpolator/test_main.cpp @@ -298,32 +298,20 @@ 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; + EXPECT_TRUE(std::abs(Ex - this->ex0) < 1e-8); + EXPECT_TRUE(std::abs(Ey - this->ey0) < 1e-8); + EXPECT_TRUE(std::abs(Ez - this->ez0) < 1e-8); + + auto const& [Bx, By, Bz] = B; + EXPECT_TRUE(std::abs(Bx - this->bx0) < 1e-8); + EXPECT_TRUE(std::abs(By - this->by0) < 1e-8); + EXPECT_TRUE(std::abs(Bz - this->bz0) < 1e-8); + } this->em.E.setBuffer("EM_E_x", nullptr); this->em.E.setBuffer("EM_E_y", nullptr); @@ -421,32 +409,19 @@ 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; + EXPECT_TRUE(std::abs(Ex - this->ex0) < 1e-8); + EXPECT_TRUE(std::abs(Ey - this->ey0) < 1e-8); + EXPECT_TRUE(std::abs(Ez - this->ez0) < 1e-8); + + auto const& [Bx, By, Bz] = B; + EXPECT_TRUE(std::abs(Bx - this->bx0) < 1e-8); + EXPECT_TRUE(std::abs(By - this->by0) < 1e-8); + EXPECT_TRUE(std::abs(Bz - this->bz0) < 1e-8); + } this->em.E.setBuffer("EM_E_x", nullptr); this->em.E.setBuffer("EM_E_y", nullptr); @@ -549,31 +524,20 @@ 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; + EXPECT_TRUE(std::abs(Ex - this->ex0) < 1e-8); + EXPECT_TRUE(std::abs(Ey - this->ey0) < 1e-8); + EXPECT_TRUE(std::abs(Ez - this->ez0) < 1e-8); + + auto const& [Bx, By, Bz] = B; + EXPECT_TRUE(std::abs(Bx - this->bx0) < 1e-8); + EXPECT_TRUE(std::abs(By - this->by0) < 1e-8); + EXPECT_TRUE(std::abs(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()); }