From 1a763a04879ab013815d20b3152939377f60b033 Mon Sep 17 00:00:00 2001 From: PhilipDeegan Date: Fri, 27 Feb 2026 17:32:18 +0100 Subject: [PATCH 1/3] Updater does not copy vector --- src/core/CMakeLists.txt | 2 - src/core/data/particles/particle_array.hpp | 53 +- src/core/errors.hpp | 1 + .../numerics/interpolator/interpolator.hpp | 60 +-- src/core/numerics/ion_updater/ion_updater.hpp | 181 ++----- src/core/numerics/pusher/boris.hpp | 452 +++++++++--------- src/core/numerics/pusher/pusher.hpp | 43 -- src/core/numerics/pusher/pusher_factory.hpp | 34 -- src/core/utilities/box/box.hpp | 59 +-- src/core/utilities/cellmap.hpp | 32 ++ src/core/utilities/indexer.hpp | 14 +- .../core/data/gridlayout/test_gridlayout.hpp | 2 +- .../test_ion_population_fixtures.hpp | 24 +- .../numerics/ion_updater/test_updater.cpp | 7 - tests/core/numerics/pusher/test_pusher.cpp | 224 +++++---- .../ion_updater/bench_ion_updater.cpp | 2 +- .../core/numerics/pusher/push_raw_use.cpp | 37 +- tools/bench/core/numerics/pusher/pusher.cpp | 38 +- 18 files changed, 581 insertions(+), 684 deletions(-) delete mode 100644 src/core/numerics/pusher/pusher.hpp delete mode 100644 src/core/numerics/pusher/pusher_factory.hpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 10273d690..8a756b839 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -28,8 +28,6 @@ set( SOURCES_INC numerics/boundary_condition/boundary_condition.hpp numerics/interpolator/interpolator.hpp numerics/pusher/boris.hpp - numerics/pusher/pusher.hpp - numerics/pusher/pusher_factory.hpp numerics/ampere/ampere.hpp numerics/faraday/faraday.hpp numerics/ohm/ohm.hpp diff --git a/src/core/data/particles/particle_array.hpp b/src/core/data/particles/particle_array.hpp index 372b90b80..bc0184c50 100644 --- a/src/core/data/particles/particle_array.hpp +++ b/src/core/data/particles/particle_array.hpp @@ -2,18 +2,19 @@ #define PHARE_CORE_DATA_PARTICLES_PARTICLE_ARRAY_HPP +#include "core/def.hpp" +#include "core/utilities/span.hpp" +#include "core/utilities/box/box.hpp" +#include "core/utilities/cellmap.hpp" +#include "core/utilities/range/range.hpp" + +#include "particle.hpp" + + +#include #include #include -#include -#include "core/utilities/indexer.hpp" -#include "particle.hpp" -#include "core/utilities/point/point.hpp" -#include "core/utilities/cellmap.hpp" -#include "core/logger.hpp" -#include "core/utilities/box/box.hpp" -#include "core/utilities/range/range.hpp" -#include "core/def.hpp" namespace PHARE::core { @@ -146,6 +147,27 @@ class ParticleArray void empty_map() { cellMap_.empty(); } + template + bool swap_last_reduce_by_one(Cell const& oldCell, std::size_t const idx) + { + // swap last to index + // swap idx to last + // --size + // return true if you need to repeat the current index == expected + + bool const idx_is_last = idx == size() - 1; + if (!idx_is_last) + { + cellMap_.swap(oldCell, particles_[size() - 1].iCell, idx, size() - 1); + particles_[idx] = particles_[size() - 1]; + } + + cellMap_.erase(*this, oldCell, size() - 1); // doesn't erase from particles vector + resize(size() - 1); + return !idx_is_last; + } + + NO_DISCARD auto nbr_particles_in(box_t const& box) const { return cellMap_.size(box); } using cell_t = std::array; @@ -180,9 +202,18 @@ class ParticleArray template - void change_icell(Cell const& newCell, std::size_t particleIndex) + void change_icell(Particle_t& /*particle*/, Cell const& oldCell, + std::size_t const particleIndex) + { + if (!box_.isEmpty()) + cellMap_.update(particles_, particleIndex, oldCell); + } + + + template + void change_icell(Cell const& newCell, std::size_t const particleIndex) { - auto oldCell = particles_[particleIndex].iCell; + auto const oldCell = particles_[particleIndex].iCell; particles_[particleIndex].iCell = newCell; if (!box_.isEmpty()) { diff --git a/src/core/errors.hpp b/src/core/errors.hpp index 9a1c6c06c..bfb7cfbdf 100644 --- a/src/core/errors.hpp +++ b/src/core/errors.hpp @@ -20,6 +20,7 @@ class DictionaryException : public std::exception public: DictionaryException() = default; DictionaryException(auto const& k, auto const& v) { (*this)(k, v); } + DictionaryException(auto const& cause) { (*this)("cause", cause); } using Dict_t = cppdict::Dict; diff --git a/src/core/numerics/interpolator/interpolator.hpp b/src/core/numerics/interpolator/interpolator.hpp index e8858d6a4..38ad8c416 100644 --- a/src/core/numerics/interpolator/interpolator.hpp +++ b/src/core/numerics/interpolator/interpolator.hpp @@ -468,41 +468,47 @@ class Interpolator : private Weighter * - then it uses Interpol<> to calculate the interpolation of E and B components * onto the particle. */ + template + void particleToMesh(Particle_t const& particle, Field& particleDensity, Field& chargeDensity, + VecField& flux, GridLayout const& layout, double coef = 1.) + { + auto const& startIndex_ = primal_startIndex_; + auto const& weights_ = primal_weights_; + auto const& [xFlux, yFlux, zFlux] = flux(); + + indexAndWeights_(layout, particle.iCell, + particle.delta); + + particleToMesh_( + particleDensity, particle, [](auto const& /*part*/) { return 1.; }, startIndex_, + weights_, coef); + particleToMesh_( + chargeDensity, particle, [](auto const& part) { return part.charge; }, startIndex_, + weights_, coef); + particleToMesh_( + xFlux, particle, [](auto const& part) { return part.v[0]; }, startIndex_, weights_, + coef); + particleToMesh_( + yFlux, particle, [](auto const& part) { return part.v[1]; }, startIndex_, weights_, + coef); + particleToMesh_( + zFlux, particle, [](auto const& part) { return part.v[2]; }, startIndex_, weights_, + coef); + } + + template inline void operator()(ParticleRange& particleRange, Field& particleDensity, Field& chargeDensity, VecField& flux, GridLayout const& layout, double coef = 1.) { - auto begin = particleRange.begin(); - auto end = particleRange.end(); - auto& startIndex_ = primal_startIndex_; - auto& weights_ = primal_weights_; - auto const& [xFlux, yFlux, zFlux] = flux(); - + auto const end = particleRange.end(); PHARE_LOG_START(3, "ParticleToMesh::operator()"); - for (auto currPart = begin; currPart != end; ++currPart) - { - indexAndWeights_(layout, currPart->iCell, - currPart->delta); - - particleToMesh_( - particleDensity, *currPart, [](auto const& part) { return 1.; }, startIndex_, - weights_, coef); - particleToMesh_( - chargeDensity, *currPart, [](auto const& part) { return part.charge; }, startIndex_, - weights_, coef); - particleToMesh_( - xFlux, *currPart, [](auto const& part) { return part.v[0]; }, startIndex_, weights_, - coef); - particleToMesh_( - yFlux, *currPart, [](auto const& part) { return part.v[1]; }, startIndex_, weights_, - coef); - particleToMesh_( - zFlux, *currPart, [](auto const& part) { return part.v[2]; }, startIndex_, weights_, - coef); - } + for (auto currPart = particleRange.begin(); currPart != end; ++currPart) + particleToMesh(*currPart, particleDensity, chargeDensity, flux, layout, coef); + PHARE_LOG_STOP(3, "ParticleToMesh::operator()"); } template diff --git a/src/core/numerics/ion_updater/ion_updater.hpp b/src/core/numerics/ion_updater/ion_updater.hpp index 57a925916..0d7441d41 100644 --- a/src/core/numerics/ion_updater/ion_updater.hpp +++ b/src/core/numerics/ion_updater/ion_updater.hpp @@ -4,10 +4,8 @@ #include "core/logger.hpp" #include "core/utilities/box/box.hpp" -#include "core/utilities/range/range.hpp" -#include "core/numerics/pusher/pusher.hpp" +#include "core/numerics/pusher/boris.hpp" #include "core/numerics/moments/moments.hpp" -#include "core/numerics/pusher/pusher_factory.hpp" #include "core/numerics/interpolator/interpolator.hpp" #include "core/numerics/boundary_condition/boundary_condition.hpp" @@ -21,39 +19,32 @@ namespace PHARE::core { enum class UpdaterMode { domain_only = 1, all = 2 }; + + template class IonUpdater { - using This = IonUpdater; - public: static constexpr auto dimension = GridLayout::dimension; static constexpr auto interp_order = GridLayout::interp_order; - using Box = PHARE::core::Box; - using Interpolator = PHARE::core::Interpolator; - using VecField = typename Ions::vecfield_type; - using ParticleArray = typename Ions::particle_array_type; - using Particle_t = typename ParticleArray::Particle_t; - using PartIterator = typename ParticleArray::iterator; - using ParticleRange = IndexRange; + using Box = PHARE::core::Box; + using Interpolator = PHARE::core::Interpolator; + using ParticleArray = Ions::particle_array_type; + using Particle_t = ParticleArray::Particle_t; + using PartIterator = ParticleArray::iterator; + using BoundaryCondition = PHARE::core::BoundaryCondition; - using Pusher = PHARE::core::Pusher; + using Pusher = BorisPusher; private: - constexpr static auto makePusher - = PHARE::core::PusherFactory::makePusher; - - std::unique_ptr pusher_; + Pusher pusher_; Interpolator interpolator_; public: - IonUpdater(PHARE::initializer::PHAREDict const& dict) - : pusher_{makePusher(dict["pusher"]["name"].template to())} - { - } + IonUpdater() = default; + IonUpdater(auto const& /*dict*/) {} template void updatePopulations(Ions& ions, Electromag const& em, Boxing_t const& boxing, double dt, @@ -63,28 +54,19 @@ class IonUpdater void updateIons(Ions& ions); - void reset() - { - // clear memory - tmp_particles_ = std::move(ParticleArray{Box{}}); - } + void reset() { /* noop */ } private: template - void updateAndDepositDomain_(Ions& ions, Electromag const& em, Boxing_t const& boxing); + void updateCopy(Ions& ions, Electromag const& em, Boxing_t const& boxing); template - void updateAndDepositAll_(Ions& ions, Electromag const& em, Boxing_t const& boxing); - - - // dealloced on regridding/load balancing coarsest - ParticleArray tmp_particles_{Box{}}; //{std::make_unique(Box{})}; + void updateInplace(Ions& ions, Electromag const& em, Boxing_t const& boxing); }; - template template void IonUpdater::updatePopulations(Ions& ions, Electromag const& em, @@ -94,15 +76,15 @@ void IonUpdater::updatePopulations(Ions& ions, Ele PHARE_LOG_SCOPE(3, "IonUpdater::updatePopulations"); resetMoments(ions); - pusher_->setMeshAndTimeStep(boxing.layout.meshSize(), dt); + pusher_.setMeshAndTimeStep(boxing.layout.meshSize(), dt); if (mode == UpdaterMode::domain_only) { - updateAndDepositDomain_(ions, em, boxing); + updateCopy(ions, em, boxing); } else { - updateAndDepositAll_(ions, em, boxing); + updateInplace(ions, em, boxing); } } @@ -121,8 +103,11 @@ struct UpdaterSelectionBoxing { auto constexpr static partGhostWidth = GridLayout::nbrParticleGhosts(); using GridLayout_t = GridLayout; + using Particle_t = IonUpdater_t::Particle_t; using Box_t = IonUpdater_t::Box; - using Selector_t = IonUpdater_t::Pusher::ParticleSelector; + using ParticleArray_t = IonUpdater_t::ParticleArray; + using ParticleRange = IndexRange; + using Selector_t = std::function; GridLayout_t const layout; std::vector const nonLevelGhostBox; @@ -131,6 +116,10 @@ struct UpdaterSelectionBoxing Selector_t const noop = [](auto& particleRange) { return particleRange; }; + bool isInGhostBox(Particle_t& p) const { return isIn(p, ghostBox); }; + bool isInDomainBox(Particle_t& p) const { return isIn(p, domainBox); }; + bool isInNonLevelGhostBox(auto const& icell) const { return isIn(icell, nonLevelGhostBox); }; + // lambda copy captures to detach from above references in case of class copy construct Selector_t const inDomainBox = [domainBox = domainBox](auto& particleRange) { return particleRange.array().partition( @@ -155,121 +144,39 @@ struct UpdaterSelectionBoxing return isIn(cell, ghostBox) and !isIn(cell, domainBox); }); }; + + Selector_t const outsideGhostBox = [ghostBox = ghostBox](auto& particleRange) { + return particleRange.array().partition( + particleRange, [&](auto const& cell) { return !isIn(cell, ghostBox); }); + }; }; -/** - * @brief IonUpdater::updateAndDepositDomain_ - evolves moments from time n to n+1 without updating particles, which stay at time n - */ + template template -void IonUpdater::updateAndDepositDomain_(Ions& ions, - Electromag const& em, - Boxing_t const& boxing) +void IonUpdater::updateCopy(Ions& ions, Electromag const& em, + Boxing_t const& boxing) { - PHARE_LOG_SCOPE(3, "IonUpdater::updateAndDepositDomain_"); + bool constexpr copy_particle = true; - auto const& layout = boxing.layout; + PHARE_LOG_SCOPE(3, "IonUpdater::updateCopy"); for (auto& pop : ions) - { - auto& domain = (tmp_particles_ = pop.domainParticles()); // make local copy - - // first push all domain particles twice - // accumulate those inNonLevelGhostBox - auto outRange = makeIndexRange(domain); - auto allowed = outRange = pusher_->move(outRange, outRange, em, pop.mass(), interpolator_, - layout, boxing.noop, boxing.inNonLevelGhostBox); - - interpolator_(allowed, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); - - - // push those in the ghostArea (i.e. stop pushing if they're not out of it) - // deposit moments on those which leave to go inDomainBox - - auto pushAndAccumulateGhosts = [&](auto const& inputArray) { - tmp_particles_ = inputArray; // work on local copy - - auto outRange = makeIndexRange(tmp_particles_); - - auto enteredInDomain = pusher_->move(outRange, outRange, em, pop.mass(), interpolator_, - layout, boxing.inGhostBox, boxing.inDomainBox); - - interpolator_(enteredInDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(), - layout); - }; - - // After this function is done domain particles overlaping ghost layers of neighbor patches - // are sent to these neighbor's patchghost particle array. - // After being pushed, some patch ghost particles may enter the domain. These need to be - // copied into the domain array so they are transfered to the neighbor patch - // ghost array and contribute to moments there too. - // On the contrary level ghost particles entering the domain here do not need to be copied - // since they contribute to nodes that are not shared with neighbor patches an since - // level border nodes will receive contributions from levelghost old and new particles - - if (pop.levelGhostParticles().size()) - pushAndAccumulateGhosts(pop.levelGhostParticles()); - } + pusher_.template move(pop, em, interpolator_, boxing); } -/** - * @brief IonUpdater::updateAndDepositDomain_ - evolves moments and particles from time n to n+1 - */ template template -void IonUpdater::updateAndDepositAll_(Ions& ions, - Electromag const& em, - Boxing_t const& boxing) +void IonUpdater::updateInplace(Ions& ions, Electromag const& em, + Boxing_t const& boxing) { - PHARE_LOG_SCOPE(3, "IonUpdater::updateAndDepositAll_"); + bool constexpr copy_particle = false; - auto const& layout = boxing.layout; + PHARE_LOG_SCOPE(3, "IonUpdater::updateInplace"); - // push domain particles, erase from array those leaving domain - // push level ghost particles that are in ghost area (==ghost box without domain) - // copy ghost particles out of ghost area that are in domain, in particle array - // finally all particles in non level ghost box are to be interpolated on mesh. for (auto& pop : ions) - { - auto& domainParticles = pop.domainParticles(); - auto domainPartRange = makeIndexRange(domainParticles); - - auto inDomain = pusher_->move(domainPartRange, domainPartRange, em, pop.mass(), - interpolator_, layout, boxing.noop, boxing.inDomainBox); - - auto now_ghosts = makeRange(domainParticles, inDomain.iend(), domainParticles.size()); - auto const not_level_ghosts = boxing.inNonLevelGhostBox(now_ghosts); - - // copy out new patch ghosts - auto& patchGhost = pop.patchGhostParticles(); - patchGhost.reserve(patchGhost.size() + not_level_ghosts.size()); - std::copy(not_level_ghosts.begin(), not_level_ghosts.end(), std::back_inserter(patchGhost)); - - domainParticles.erase(now_ghosts); // drop all ghosts - - if (pop.levelGhostParticles().size()) - { - auto particleRange = makeIndexRange(pop.levelGhostParticles()); - auto inGhostLayerRange - = pusher_->move(particleRange, particleRange, em, pop.mass(), interpolator_, layout, - boxing.inGhostBox, boxing.inGhostLayer); - - auto& particleArray = particleRange.array(); - particleArray.export_particles( - domainParticles, [&](auto const& cell) { return isIn(cell, boxing.domainBox); }); - - particleArray.erase( - makeRange(particleArray, inGhostLayerRange.iend(), particleArray.size())); - } - - interpolator_( // - domainParticles, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); - interpolator_( // - patchGhost, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); - } + pusher_.template move(pop, em, interpolator_, boxing); } diff --git a/src/core/numerics/pusher/boris.hpp b/src/core/numerics/pusher/boris.hpp index 79b2cc296..6b94ab619 100644 --- a/src/core/numerics/pusher/boris.hpp +++ b/src/core/numerics/pusher/boris.hpp @@ -4,143 +4,133 @@ #include "core/errors.hpp" #include "core/logger.hpp" -#include "core/numerics/pusher/pusher.hpp" +#include "core/utilities/box/box.hpp" #include #include #include -#include #include #include #include +#include -namespace PHARE::core + +namespace PHARE::core::boris { -template -class BorisPusher - : public Pusher +struct MoveTwoCellException : std::exception { - struct MoveTwoCellException : std::exception + MoveTwoCellException(double const d, double const v) + : delta{d} + , vel{v} { - MoveTwoCellException(double const d, double const v) - : delta{d} - , vel{v} - { - } + } - double delta, vel; - }; + double delta, vel; +}; -public: - using Super - = Pusher; +template +auto static advance(Particle_t& p, std::array const& halfDtOverDl) +{ + std::array newCell; + for (std::size_t iDim = 0; iDim < dim; ++iDim) + { + double const delta = p.delta[iDim] + static_cast(halfDtOverDl[iDim] * p.v[iDim]); -private: - using ParticleSelector = typename Super::ParticleSelector; + if (std::abs(delta) > 2) + throw MoveTwoCellException{delta, p.v[iDim]}; -public: - // This move function should be considered when being used so that all particles are pushed - // twice - see: https://github.com/PHAREHUB/PHARE/issues/571 - /** see Pusher::move() documentation*/ -#if 0 - ParticleRange move(ParticleRange const& rangeIn, ParticleRange& rangeOut, - Electromag const& emFields, double mass, Interpolator& interpolator, - ParticleSelector const& particleIsNotLeaving, BoundaryCondition& bc, - GridLayout const& layout) override - { - // push the particles of half a step - // rangeIn : t=n, rangeOut : t=n+1/Z - // get a pointer on the first particle of rangeOut that leaves the patch - auto firstLeaving - = pushStep_(rangeIn, rangeOut, particleIsNotLeaving, PushStep::PrePush); + auto const iCell = static_cast(std::floor(delta)); - // apply boundary condition on the particles in [firstLeaving, rangeOut.end[ - // that actually leave through a physical boundary condition - // get a pointer on the new end of rangeOut. Particles passed newEnd - // are those that have left the patch through a non-physical boundary - // they should be discarded now - auto newEnd = bc.applyOutgoingParticleBC(firstLeaving, rangeOut.end()); + p.delta[iDim] = delta - iCell; + newCell[iDim] = iCell + p.iCell[iDim]; + } + return newCell; +} - rangeOut = makeRange(rangeOut.begin(), std::move(newEnd)); - // get electromagnetic fields interpolated on the particles of rangeOut - // stop at newEnd. - interpolator(rangeOut.begin(), rangeOut.end(), emFields, layout); - // get the particle velocity from t=n to t=n+1 - accelerate_(rangeOut, rangeOut, mass); +/** Accelerate the particles in rangeIn and put the new velocity in rangeOut + */ +template +void accelerate(Particle& p, ParticleEB const& eb, Float const& dto2m) +{ + static constexpr Float one = 1; + static constexpr Float two = 2; - // 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 - firstLeaving = pushStep_(rangeOut, rangeOut, particleIsNotLeaving, PushStep::PostPush); + auto& [pE, pB] = eb; + auto& [pEx, pEy, pEz] = pE; + auto& [pBx, pBy, pBz] = pB; - // apply BC on the leaving particles that leave through physical BC - // and get pointer on new End, discarding particles leaving elsewhere - newEnd = bc.applyOutgoingParticleBC(firstLeaving, rangeOut.end()); + Float const coef1 = p.charge * dto2m; - rangeOut = makeRange(rangeOut.begin(), std::move(newEnd)); + // We now apply the 3 steps of the BORIS PUSHER + // 1st half push of the electric field + Float const velx1 = p.v[0] + coef1 * pEx; + Float const vely1 = p.v[1] + coef1 * pEy; + Float const velz1 = p.v[2] + coef1 * pEz; - return rangeOut.end(); - } -#endif + // preparing variables for magnetic rotation + Float const rx = coef1 * pBx; + Float const ry = coef1 * pBy; + Float const rz = coef1 * pBz; + Float const rx2 = rx * rx; + Float const ry2 = ry * ry; + Float const rz2 = rz * rz; + Float const rxry = rx * ry; + Float const rxrz = rx * rz; + Float const ryrz = ry * rz; - ParticleRange move(ParticleRange const& rangeIn, ParticleRange& rangeOut, - Electromag const& emFields, double mass, Interpolator& interpolator, - GridLayout const& layout, ParticleSelector firstSelector, - ParticleSelector secondSelector) override - { - PHARE_LOG_SCOPE(3, "Boris::move_no_bc"); + Float const invDet = one / (one + rx2 + ry2 + rz2); - // push the particles of half a step - // 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 - prePushStep_(rangeIn, rangeOut); + // preparing rotation matrix due to the magnetic field + // m = invDet*(I + r*r - r x I) - I where x denotes the cross product + Float const mxx = one + rx2 - ry2 - rz2; + Float const mxy = two * (rxry + rz); + Float const mxz = two * (rxrz - ry); - rangeOut = firstSelector(rangeOut); + Float const myx = two * (rxry - rz); + Float const myy = one + ry2 - rx2 - rz2; + Float const myz = two * (ryrz + rx); - double const dto2m = 0.5 * dt_ / mass; - for (auto idx = rangeOut.ibegin(); idx < rangeOut.iend(); ++idx) - { - auto& currPart = rangeOut.array()[idx]; + Float const mzx = two * (rxrz + ry); + Float const mzy = two * (ryrz - rx); + Float const mzz = one + rz2 - rx2 - ry2; - // get electromagnetic fields interpolated on the particles of rangeOut stop at newEnd. - // get the particle velocity from t=n to t=n+1 - auto const& local_em = interpolator(currPart, emFields, layout); - accelerate_(currPart, local_em, dto2m); + // magnetic rotation + Float const velx2 = (mxx * velx1 + mxy * vely1 + mxz * velz1) * invDet; + Float const vely2 = (myx * velx1 + myy * vely1 + myz * velz1) * invDet; + Float const velz2 = (mzx * velx1 + mzy * vely1 + mzz * velz1) * invDet; - // 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 - try - { - postPushStep_(rangeOut, idx); - } - catch (DictionaryException const& bex) - { - auto ex = bex; - auto const& [e, b] = local_em; - for (std::uint16_t i = 0; i < 3; ++i) - ex("E_" + std::to_string(i), std::to_string(e[i])); - for (std::uint16_t i = 0; i < 3; ++i) - ex("B_" + std::to_string(i), std::to_string(b[i])); - ex("level", std::to_string(layout.levelNumber())); - throw ex; - } - } + // 2nd half push of the electric field / Update particle velocity + p.v[0] = velx2 + coef1 * pEx; + p.v[1] = vely2 + coef1 * pEy; + p.v[2] = velz2 + coef1 * pEz; +} + + +} // namespace PHARE::core::boris + + +namespace PHARE::core +{ - return secondSelector(rangeOut); - } +template +class BorisPusher +{ + using Particle_t = ParticleArray_t::value_type; + using OnChangeIcell = std::function const&, std::size_t const)>; - /** see Pusher::move() documentation*/ - void setMeshAndTimeStep(std::array ms, double const ts) override +public: + void setMeshAndTimeStep(std::array ms, double const ts) { std::transform(std::begin(ms), std::end(ms), std::begin(halfDtOverDl_), [ts](double& x) { return 0.5 * ts / x; }); @@ -148,164 +138,194 @@ class BorisPusher } + template + void move(Population& pop, Electromag const& em, Interpolator& interpolator, + Boxing_t const& boxing) + { + move_domain(pop, em, interpolator, boxing); + move_level_ghost(pop, em, interpolator, boxing); + } private: - /** move the particle partIn of half a time step and store it in partOut - */ - template - auto advancePosition_(Particle const& partIn, Particle& partOut) + void move_particle(auto& particle, auto const& em, auto& interpolator, auto const& layout, + auto const& dto2m) { - std::array newCell; - for (std::size_t iDim = 0; iDim < dim; ++iDim) + try { - double const delta - = partIn.delta[iDim] + static_cast(halfDtOverDl_[iDim] * partIn.v[iDim]); - - if (std::abs(delta) > 2) - throw MoveTwoCellException{delta, partIn.v[iDim]}; - - auto const iCell = static_cast(std::floor(delta)); - partOut.delta[iDim] = delta - iCell; - newCell[iDim] = iCell + partIn.iCell[iDim]; + particle.iCell = boris::advance(particle, halfDtOverDl_); } - return newCell; - } - - - /** advance the particles in rangeIn of half a time step and store them - * in rangeOut. - * @return the function returns and iterator on the first leaving particle, as - * detected by the ParticleSelector - */ - void prePushStep_(ParticleRange const& rangeIn, ParticleRange& rangeOut) - { - auto& inParticles = rangeIn.array(); - auto& outParticles = rangeOut.array(); - for (auto inIdx = rangeIn.ibegin(), outIdx = rangeOut.ibegin(); inIdx < rangeIn.iend(); - ++inIdx, ++outIdx) + catch (boris::MoveTwoCellException const& e) { - // in the first push, this is the first time - // we push to rangeOut, which contains crap - // the push will only touch the particle position - // but the next step being the acceleration of - // rangeOut, we need to copy rangeIn weights, charge - // and velocity. This is done here although - // not strictly speaking this function's business - // to take advantage that we're already looping - // over rangeIn particles. - - outParticles[outIdx].charge = inParticles[inIdx].charge; - outParticles[outIdx].weight = inParticles[inIdx].weight; - outParticles[outIdx].v = inParticles[inIdx].v; - - try - { - auto newCell = advancePosition_(inParticles[inIdx], outParticles[outIdx]); - if (newCell != inParticles[inIdx].iCell) - outParticles.change_icell(newCell, outIdx); - } - catch (MoveTwoCellException const& e) - { - std::stringstream ss; - ss << "PrePush Particle moved 2 cells with delta/vel: "; - ss << e.delta << "/" << e.vel << std::endl; - DictionaryException ex{"cause", ss.str()}; - throw ex; - } + std::stringstream ss; + ss << "PrePush Particle moved 2 cells with delta/vel: "; + ss << e.delta << "/" << e.vel << std::endl; + throw DictionaryException{ss.str()}; } - } - void postPushStep_(ParticleRange& range, std::size_t idx) - { + auto const& local_em = interpolator(particle, em, layout); try { - auto& particles = range.array(); - auto newCell = advancePosition_(particles[idx], particles[idx]); - if (newCell != particles[idx].iCell) - particles.change_icell(newCell, idx); + boris::accelerate(particle, local_em, dto2m); + particle.iCell = boris::advance(particle, halfDtOverDl_); } - catch (MoveTwoCellException const& e) + catch (boris::MoveTwoCellException const& e) { std::stringstream ss; ss << "PostPush Particle moved 2 cells with delta/vel: "; ss << e.delta << "/" << e.vel << std::endl; - throw DictionaryException{}("cause", ss.str()); + DictionaryException ex{ss.str()}; + + auto const& [E, B] = local_em; + for (std::uint16_t i = 0; i < 3; ++i) + ex("E_" + std::to_string(i), std::to_string(E[i])); + for (std::uint16_t i = 0; i < 3; ++i) + ex("B_" + std::to_string(i), std::to_string(B[i])); + ex("level", std::to_string(layout.levelNumber())); + throw ex; } - } - + }; - /** Accelerate the particles in rangeIn and put the new velocity in rangeOut - */ - template - void accelerate_(Particle_t& part, ParticleEB const& particleEB, double const& dto2m) + template + void move_level_ghost(Population& pop, Electromag const& em, Interpolator& interpolator, + Boxing_t const& boxing) { - auto& [pE, pB] = particleEB; - auto& [pEx, pEy, pEz] = pE; - auto& [pBx, pBy, pBz] = pB; + PHARE_LOG_SCOPE(3, "Boris::move_level_ghost"); + + using ParticleLoop_t = std::conditional_t; + + auto const& layout = boxing.layout; + double const dto2m = 0.5 * dt_ / pop.mass(); + + auto& domain = pop.domainParticles(); + auto& level_ghost = pop.levelGhostParticles(); + + for (std::size_t i = 0; i < level_ghost.size(); ++i) // size might change on iteration! + { + ParticleLoop_t particle = level_ghost[i]; + auto const oldCell = particle.iCell; + move_particle(particle, em, interpolator, layout, dto2m); - double const coef1 = part.charge * dto2m; + if (oldCell == particle.iCell) + continue; - // We now apply the 3 steps of the BORIS PUSHER + bool const isInDomainBox = boxing.isInDomainBox(particle); + bool const should_interpolate = isInDomainBox; - // 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; + if (should_interpolate) + interpolator.particleToMesh( // + particle, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + if constexpr (!copy_particle) + { + if (isInDomainBox) + { + domain.push_back(particle); + level_ghost.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + continue; + } + + bool const isInNonLevelGhostBox + = isInDomainBox || boxing.isInNonLevelGhostBox(particle.iCell); + bool const isInLevelGhostBox = !isInNonLevelGhostBox; + + if (isInLevelGhostBox) + level_ghost.change_icell(particle, oldCell, i); + else + { + level_ghost.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + } + } + } + } - // preparing variables for magnetic rotation - double const rx = coef1 * pBx; - double const ry = coef1 * pBy; - double const rz = coef1 * pBz; + template + void move_domain(Population& pop, Electromag const& em, Interpolator& interpolator, + Boxing_t const& boxing) + { + PHARE_LOG_SCOPE(3, "Boris::move_domain"); - 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; + using ParticleLoop_t = std::conditional_t; - 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); + auto const& layout = boxing.layout; + double const dto2m = 0.5 * dt_ / pop.mass(); - double const myx = 2. * (rxry - rz); - double const myy = 1. + ry2 - rx2 - rz2; - double const myz = 2. * (ryrz + rx); + auto& domain = pop.domainParticles(); + auto& patch_ghost = pop.patchGhostParticles(); - double const mzx = 2. * (rxrz + ry); - double const mzy = 2. * (ryrz - rx); - double const mzz = 1. + rz2 - rx2 - ry2; + for (std::size_t i = 0; i < domain.size(); ++i) // size might change on iteration! + { + ParticleLoop_t particle = domain[i]; + auto const oldCell = particle.iCell; - // 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; + move_particle(particle, em, interpolator, layout, dto2m); + if (oldCell == particle.iCell) + { + interpolator.particleToMesh( // + particle, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); - // 2nd half push of the electric field - velx1 = velx2 + coef1 * pEx; - vely1 = vely2 + coef1 * pEy; - velz1 = velz2 + coef1 * pEz; + continue; + } - // Update particle velocity - part.v[0] = velx1; - part.v[1] = vely1; - part.v[2] = velz1; - } + bool const isInDomainBox = boxing.isInDomainBox(particle); + bool const isInNonLevelGhostBox + = isInDomainBox || boxing.isInNonLevelGhostBox(particle.iCell); + bool const should_interpolate = isInNonLevelGhostBox; + if (!should_interpolate) + { + if constexpr (!copy_particle) + { + domain.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + } + continue; + } + + interpolator.particleToMesh( // + particle, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + + if constexpr (!copy_particle) + domain.change_icell(particle, oldCell, i); + } + // move to patch_ghost + if constexpr (!copy_particle) + { + auto range = makeIndexRange(domain); + range = domain.partition( + range, [&](auto const& cell) { return core::isIn(cell, boxing.domainBox); }); + + auto const not_in_domain = makeRange(domain, range.iend(), domain.size()); + auto const is_patch_ghost = domain.partition( + not_in_domain, [&](auto const& cell) { return boxing.isInNonLevelGhostBox(cell); }); + + patch_ghost.reserve(patch_ghost.size() + is_patch_ghost.size()); + std::copy(is_patch_ghost.begin(), is_patch_ghost.end(), + std::back_inserter(patch_ghost)); + domain.erase(not_in_domain); + } + PHARE_DEBUG_DO({ + for (auto const& p : domain) + if (!isIn(p, boxing.domainBox)) + throw std::runtime_error("invalid domain"); + for (auto const& p : patch_ghost) + if (isIn(p, boxing.domainBox)) + throw std::runtime_error("invalid patch ghost"); + }) + } std::array halfDtOverDl_; double dt_; }; + + } // namespace PHARE::core diff --git a/src/core/numerics/pusher/pusher.hpp b/src/core/numerics/pusher/pusher.hpp deleted file mode 100644 index 1e094ffd0..000000000 --- a/src/core/numerics/pusher/pusher.hpp +++ /dev/null @@ -1,43 +0,0 @@ -#ifndef PHARE_CORE_NUMERICS_PUSHER_PUSHER_HPP -#define PHARE_CORE_NUMERICS_PUSHER_PUSHER_HPP - -#include -#include -#include -#include - -#include "core/utilities/range/range.hpp" -#include "core/data/particles/particle.hpp" - -namespace PHARE -{ -namespace core -{ - template - class Pusher - { - protected: - static auto constexpr dimension = GridLayout::dimension; - - public: - using ParticleSelector = std::function; - - // TODO : to really be independant on boris which has 2 push steps - // we should have an arbitrary number of selectors, 1 per push step - virtual ParticleRange move(ParticleRange const& rangeIn, ParticleRange& rangeOut, - Electromag const& emFields, double mass, - Interpolator& interpolator, GridLayout const& layout, - ParticleSelector firstSelector, ParticleSelector secondSelector) - = 0; - - - virtual void setMeshAndTimeStep(std::array ms, double ts) = 0; - - virtual ~Pusher() {} - }; - -} // namespace core -} // namespace PHARE - -#endif diff --git a/src/core/numerics/pusher/pusher_factory.hpp b/src/core/numerics/pusher/pusher_factory.hpp deleted file mode 100644 index 28a85ecb2..000000000 --- a/src/core/numerics/pusher/pusher_factory.hpp +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef PHARE_CORE_NUMERIC_PUSHER_PUSHER_FACTORY_HPP -#define PHARE_CORE_NUMERIC_PUSHER_PUSHER_FACTORY_HPP - -#include -#include -#include - -#include "boris.hpp" -#include "pusher.hpp" - -namespace PHARE -{ -namespace core -{ - class PusherFactory - { - public: - template - static auto makePusher(std::string pusherName) - { - if (pusherName == "modified_boris") - { - return std::make_unique>(); - } - - throw std::runtime_error("Error : Invalid Pusher name"); - } - }; -} // namespace core -} // namespace PHARE - -#endif // PUSHER_FACTORY_HPP diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index 57bdd9ecf..d597eda57 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -259,6 +259,22 @@ struct boxes_iterator }; +template typename Point, typename Type, std::size_t SIZE> +NO_DISCARD bool isIn(Point const& point, Box const& box) +{ + for (std::size_t iDim = 0; iDim < SIZE; ++iDim) + if (point[iDim] < box.lower[iDim] || point[iDim] > box.upper[iDim]) + return false; + return true; +} + + +template +NO_DISCARD auto isIn(Particle const& particle, Box const& box) + -> decltype(isIn(particle.iCell, box), bool()) +{ + return isIn(particle.iCell, box); +} /** this overload of isIn takes a Point and a Container of boxes * and returns true if the Point is at least in one of the boxes. @@ -275,48 +291,9 @@ bool isIn(Point const& point, BoxContainer const& boxes) typename BoxContainer::value_type::value_type>::value, "Box and Point should have the same data type"); - - auto isIn1D = [](typename Point::value_type pos, typename Point::value_type lower, - typename Point::value_type upper) { return pos >= lower && pos <= upper; }; - for (auto const& box : boxes) - { - bool pointInBox = true; - - for (auto iDim = 0u; iDim < Point::dimension; ++iDim) - { - pointInBox = pointInBox && isIn1D(point[iDim], box.lower[iDim], box.upper[iDim]); - } - if (pointInBox) - return pointInBox; - } - - return false; -} - -template -NO_DISCARD auto isIn(Particle const& particle, Box const& box) - -> decltype(isIn(particle.iCell, box), bool()) -{ - return isIn(particle.iCell, box); -} - -/** This overload of isIn does the same as the one above but takes only - * one box. - */ -template typename Point, typename Type, std::size_t SIZE> -NO_DISCARD bool isIn(Point const& point, Box const& box) -{ - auto isIn1D = [](auto const pos, auto const lower, auto const upper) { - return pos >= lower && pos <= upper; - }; - - bool pointInBox = true; - - for (auto iDim = 0u; iDim < SIZE; ++iDim) - pointInBox = pointInBox && isIn1D(point[iDim], box.lower[iDim], box.upper[iDim]); - if (pointInBox) - return pointInBox; + if (isIn(point, box)) + return true; return false; } diff --git a/src/core/utilities/cellmap.hpp b/src/core/utilities/cellmap.hpp index 86ab52a48..394b8df6a 100644 --- a/src/core/utilities/cellmap.hpp +++ b/src/core/utilities/cellmap.hpp @@ -180,6 +180,10 @@ class CellMap CellExtractor extract = default_extractor); + template, void>> + void erase(Array const& items, cell_t const oldCell, std::size_t const itemIndex); + + // sort all cell indexes void sort(); @@ -204,6 +208,8 @@ class CellMap return cellIndexes_(local_(cell)); } + void swap(cell_t const a, cell_t const b, std::size_t const i0, std::size_t const i1); + private: template auto local_(Cell const& cell) const @@ -220,6 +226,24 @@ class CellMap }; +template +void CellMap::swap(cell_t const a, cell_t const b, std::size_t const i0, + std::size_t const i1) +{ + auto& indexesA = cellIndexes_(local_(a)); + auto& indexesB = cellIndexes_(local_(b)); + + assert(indexesA.is_indexed(i0)); + assert(indexesB.is_indexed(i1)); + + // find before modify! + auto const posA = indexesA.position(i0); + auto const posB = indexesB.position(i1); + + indexesA.update(posA, i1); + indexesB.update(posB, i0); +} + template template @@ -273,6 +297,14 @@ inline void CellMap::erase(Array const& items, std::size_t it } +template +template +inline void CellMap::erase(Array const& items, cell_t const oldCell, + std::size_t const itemIndex) +{ + auto& blist = cellIndexes_(local_(oldCell)); + blist.remove(itemIndex); +} template diff --git a/src/core/utilities/indexer.hpp b/src/core/utilities/indexer.hpp index f701748bb..94905266b 100644 --- a/src/core/utilities/indexer.hpp +++ b/src/core/utilities/indexer.hpp @@ -49,7 +49,7 @@ class Indexer // void trim(std::size_t max_empty); // to use if an item in an indexed array is moved at another index - void updateIndex(std::size_t oldIndex, std::size_t newIndex) + void updateIndex(std::size_t const oldIndex, std::size_t const newIndex) { // auto it = std::find(std::begin(indexes_), std::end(indexes_), oldIndex); @@ -59,6 +59,18 @@ class Indexer } } + void update(std::size_t const offset, std::size_t const index) { indexes_[offset] = index; } + + auto& findIndex(std::size_t const index) const + { + // assumes exists + return *std::find(std::begin(indexes_), std::end(indexes_), index); + } + auto position(std::size_t const index) const + { + return std::distance(indexes_.data(), &findIndex(index)); + } + // empty the bucketlist, but leaves the capacity untouched void empty() { indexes_.resize(0); }; NO_DISCARD bool is_empty() const { return indexes_.size() == 0; } diff --git a/tests/core/data/gridlayout/test_gridlayout.hpp b/tests/core/data/gridlayout/test_gridlayout.hpp index d31f630a6..4137b9249 100644 --- a/tests/core/data/gridlayout/test_gridlayout.hpp +++ b/tests/core/data/gridlayout/test_gridlayout.hpp @@ -14,7 +14,7 @@ class TestGridLayout : public GridLayout TestGridLayout() = default; - TestGridLayout(std::uint32_t cells) + TestGridLayout(std::uint32_t const cells) : GridLayout{PHARE::core::ConstArray(1.0 / cells), PHARE::core::ConstArray(cells), PHARE::core::Point{PHARE::core::ConstArray(0)}} diff --git a/tests/core/data/ion_population/test_ion_population_fixtures.hpp b/tests/core/data/ion_population/test_ion_population_fixtures.hpp index f8317268c..54b0083f3 100644 --- a/tests/core/data/ion_population/test_ion_population_fixtures.hpp +++ b/tests/core/data/ion_population/test_ion_population_fixtures.hpp @@ -40,7 +40,7 @@ struct UsableIonsDefaultTypes template -class UsableIonsPopulation : public _defaults::IonPopulation_t +class UsableIonsPopulation_ : public _defaults::IonPopulation_t { using GridLayout_t = _defaults::GridLayout_t; using VecField_t = _defaults::VecField_t; @@ -49,7 +49,7 @@ class UsableIonsPopulation : public _defaults::IonPopulation_t using Super = IonPopulation; public: - UsableIonsPopulation(initializer::PHAREDict const& dict, GridLayout_t const& layout) + UsableIonsPopulation_(initializer::PHAREDict const& dict, GridLayout_t const& layout) : Super{dict} , particleDensity{this->name() + "_particleDensity", layout, HybridQuantity::Scalar::rho} , chargeDensity{this->name() + "_chargeDensity", layout, HybridQuantity::Scalar::rho} @@ -76,10 +76,12 @@ class UsableIonsPopulation : public _defaults::IonPopulation_t _defaults::UsableTensorField_t M; UsableParticlesPopulation particles; }; +template +using UsableIonsPopulation = UsableIonsPopulation_>; template -class UsableIons +class UsableIons_ : public Ions { using GridLayout_t = _defaults::GridLayout_t; @@ -105,7 +107,7 @@ class UsableIons } public: - UsableIons(GridLayout_t const& layout, std::vector const& pop_names) + UsableIons_(GridLayout_t const& layout, std::vector const& pop_names) : Super{super(pop_names)} , massDensity{"massDensity", layout, HybridQuantity::Scalar::rho} , chargeDensity{"chargeDensity", layout, HybridQuantity::Scalar::rho} @@ -125,8 +127,8 @@ class UsableIons Super::getRunTimeResourcesViewList().emplace_back(*pop); } - UsableIons(GridLayout_t const& layout, std::string const& pop_name) - : UsableIons{layout, std::vector{pop_name}} + UsableIons_(GridLayout_t const& layout, std::string const& pop_name = "protons") + : UsableIons_{layout, std::vector{pop_name}} { } @@ -138,16 +140,10 @@ class UsableIons _defaults::Grid_t massDensity, chargeDensity; _defaults::UsableVecField_t Vi; _defaults::UsableTensorField_t M; - std::vector> populations; + std::vector> populations; }; - -template -using UsableIonsPopulation_t - = UsableIonsPopulation>; - - template -using UsableIons_t = UsableIons>; +using UsableIons = UsableIons_>; } // namespace PHARE::core diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index e4ffcbc32..4066b332b 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -664,13 +664,6 @@ TYPED_TEST_SUITE(IonUpdaterTest, DimInterps, ); -TYPED_TEST(IonUpdaterTest, ionUpdaterTakesPusherParamsFromPHAREDictAtConstruction) -{ - typename IonUpdaterTest::IonUpdater ionUpdater{ - init_dict["simulation"]["algo"]["ion_updater"]}; -} - - // the following 3 tests are testing the fixture is well configured. diff --git a/tests/core/numerics/pusher/test_pusher.cpp b/tests/core/numerics/pusher/test_pusher.cpp index ceef4c9a5..485126572 100644 --- a/tests/core/numerics/pusher/test_pusher.cpp +++ b/tests/core/numerics/pusher/test_pusher.cpp @@ -2,6 +2,10 @@ #include "gtest/gtest.h" #include +#include +#include +#include +#include #include #include #include @@ -13,10 +17,12 @@ #include "core/data/particles/particle_array.hpp" #include "core/numerics/boundary_condition/boundary_condition.hpp" #include "core/numerics/pusher/boris.hpp" -#include "core/numerics/pusher/pusher_factory.hpp" #include "core/utilities/range/range.hpp" #include "core/utilities/box/box.hpp" +#include "tests/core/data/gridlayout/test_gridlayout.hpp" +#include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" + using namespace PHARE::core; @@ -69,7 +75,7 @@ Trajectory readExpectedTrajectory() // in the python script that generates a trajectory // this way, we don't need a proper Electromag, VecField, Field etc. objects // to test the pusher. -class Interpolator +class TestInterpolator { using E_B_tuple = std::tuple, std::array>; @@ -91,15 +97,39 @@ class Interpolator return eb_interop; } + + void particleToMesh(auto&&... args) {} }; // mock of electromag just so that the Pusher gives something to // the Interpolator -class Electromag +class TestElectromag { }; +template +struct TestIons // Some tests have large domains but no need for fields +{ + using particle_array_type = ParticleArray_; + + struct TestIonPop + { + auto& domainParticles() { return domain; } + auto& levelGhostParticles() { return levelGhost; } + auto& patchGhostParticles() { return patchGhost; } + double mass() const { return 1; } + + GridLayout const& layout; + ParticleArray_ domain{layout.AMRBox()}; + ParticleArray_ patchGhost{layout.AMRBox()}; + ParticleArray_ levelGhost{layout.AMRBox()}; + }; + + GridLayout layout; + TestIonPop pop{layout}; +}; + // with this mock, all particles are found inside class DummySelector @@ -113,27 +143,29 @@ class DummySelector }; -template -struct DummyLayout -{ - static constexpr std::size_t dimension = dimension_; - std::array nbrCells_; - auto nbrCells() const { return nbrCells_; } - auto AMRBox() const { return PHARE::core::emptyBox(); } - auto levelNumber() const { return 0; } -}; template class APusher : public ::testing::Test { + static constexpr auto dimension = dim; + static constexpr auto interp_order = 1; + constexpr static PHARE::SimOpts opts{dimension, interp_order}; + + using PHARE_TYPES = PHARE::core::PHARE_Types; + using Interpolator = TestInterpolator; + using Electromag = TestElectromag; + using GridLayout_t = TestGridLayout; + using Ions_t = TestIons>; + using IonUpdater = typename PHARE::core::IonUpdater; + using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + public: - using Pusher_ = BorisPusher>, Electromag, Interpolator, - BoundaryCondition, DummyLayout>; + using Pusher_ = BorisPusher, Electromag, Interpolator, + BoundaryCondition, GridLayout_t>; APusher() : expectedTrajectory{readExpectedTrajectory()} - , particlesIn{layout.AMRBox()} - , particlesOut(layout.AMRBox()) + , layout{30} , pusher{std::make_unique()} , mass{1} , dt{0.0001} @@ -142,22 +174,22 @@ class APusher : public ::testing::Test , nt{static_cast((tend - tstart) / dt + 1)} { - particlesIn.emplace_back( - Particle{1., 1., ConstArray(5), ConstArray(0.), {0., 10., 0.}}); - particlesOut.emplace_back( + particles.emplace_back( Particle{1., 1., ConstArray(5), ConstArray(0.), {0., 10., 0.}}); dxyz.fill(0.05); for (std::size_t i = 0; i < dim; i++) actual[i].resize(nt, 0.05); pusher->setMeshAndTimeStep(dxyz, dt); + + std::transform(std::begin(dxyz), std::end(dxyz), std::begin(halfDtOverDl), + [dt = this->dt](double& x) { return 0.5 * dt / x; }); } protected: using Particle = typename ParticleArray::Particle_t; Trajectory expectedTrajectory; - DummyLayout layout; - ParticleArray particlesIn; - ParticleArray particlesOut; + GridLayout_t layout; + std::unique_ptr pusher; double mass; double dt; @@ -165,11 +197,26 @@ class APusher : public ::testing::Test std::size_t nt; Electromag em; Interpolator interpolator; - DummySelector selector; // BoundaryCondition bc; std::array, dim> actual; std::array dxyz; + Ions_t ions{layout}; + + ParticleArray& particles = ions.pop.domainParticles(); + + std::array halfDtOverDl; + double const dto2m = 0.5 * dt / mass; + + void move() + { + for (auto& particle : particles) + { + particle.iCell = boris::advance(particle, halfDtOverDl); + boris::accelerate(particle, interpolator(particle, em, layout), dto2m); + particle.iCell = boris::advance(particle, halfDtOverDl); + } + } }; @@ -179,19 +226,13 @@ using APusher3D = APusher<3>; TEST_F(APusher3D, trajectoryIsOk) { - auto rangeIn = makeIndexRange(particlesIn); - auto rangeOut = makeIndexRange(particlesOut); - std::copy(rangeIn.begin(), rangeIn.end(), rangeOut.begin()); - for (decltype(nt) i = 0; i < nt; ++i) { - 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]; - actual[2][i] = (particlesOut[0].iCell[2] + particlesOut[0].delta[2]) * dxyz[2]; - - pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); + actual[0][i] = (particles[0].iCell[0] + particles[0].delta[0]) * dxyz[0]; + actual[1][i] = (particles[0].iCell[1] + particles[0].delta[1]) * dxyz[1]; + actual[2][i] = (particles[0].iCell[2] + particles[0].delta[2]) * dxyz[2]; - std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); + move(); } EXPECT_THAT(actual[0], ::testing::Pointwise(::testing::DoubleNear(1e-5), expectedTrajectory.x)); @@ -204,18 +245,12 @@ TEST_F(APusher3D, trajectoryIsOk) TEST_F(APusher2D, trajectoryIsOk) { - auto rangeIn = makeIndexRange(particlesIn); - auto rangeOut = makeIndexRange(particlesOut); - std::copy(rangeIn.begin(), rangeIn.end(), rangeOut.begin()); - for (decltype(nt) i = 0; i < nt; ++i) { - 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]; + actual[0][i] = (particles[0].iCell[0] + particles[0].delta[0]) * dxyz[0]; + actual[1][i] = (particles[0].iCell[1] + particles[0].delta[1]) * dxyz[1]; - pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); - - std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); + move(); } EXPECT_THAT(actual[0], ::testing::Pointwise(::testing::DoubleNear(1e-5), expectedTrajectory.x)); @@ -226,17 +261,11 @@ TEST_F(APusher2D, trajectoryIsOk) TEST_F(APusher1D, trajectoryIsOk) { - auto rangeIn = makeIndexRange(particlesIn); - auto rangeOut = makeIndexRange(particlesOut); - std::copy(rangeIn.begin(), rangeIn.end(), rangeOut.begin()); - for (decltype(nt) i = 0; i < nt; ++i) { - actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0]; + actual[0][i] = (particles[0].iCell[0] + particles[0].delta[0]) * dxyz[0]; - pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); - - std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); + move(); } EXPECT_THAT(actual[0], ::testing::Pointwise(::testing::DoubleNear(1e-5), expectedTrajectory.x)); @@ -250,21 +279,32 @@ TEST_F(APusher1D, trajectoryIsOk) // and those that stay. class APusherWithLeavingParticles : public ::testing::Test { + using Interpolator = TestInterpolator; + using Electromag = TestElectromag; + + static constexpr auto dimension = 1; + static constexpr auto interp_order = 1; + constexpr static PHARE::SimOpts opts{dimension, interp_order}; + + using PHARE_TYPES = PHARE::core::PHARE_Types; + using ParticleArray_t = ParticleArray<1>; + using GridLayout_t = TestGridLayout; + using Ions_t = PHARE_TYPES::Ions_t; + using IonUpdater = typename PHARE::core::IonUpdater; + using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + public: APusherWithLeavingParticles() - : pusher{std::make_unique< - BorisPusher<1, IndexRange>, Electromag, Interpolator, - BoundaryCondition<1, 1>, DummyLayout<1>>>()} + : pusher{std::make_unique, TestElectromag, TestInterpolator, + BoundaryCondition<1, 1>, GridLayout_t>>()} , mass{1} , dt{0.001} , tstart{0} , tend{10} , nt{static_cast((tend - tstart) / dt + 1)} - , domain{Point{0.}, Point{1.}} - , cells{Point{0}, Point{9}} - , particlesIn{grow(cells, 1)} - , particlesOut1{grow(cells, 1), 1000} - , particlesOut2{grow(cells, 1), 1000} + , layout{10} + // , particlesOut1{grow(domain, 1), 1000} + // , particlesOut2{grow(domain, 1), 1000} { std::random_device rd; std::mt19937 gen(rd()); @@ -281,8 +321,8 @@ class APusherWithLeavingParticles : public ::testing::Test protected: - std::unique_ptr>, Electromag, Interpolator, - BoundaryCondition<1, 1>, DummyLayout<1>>> + std::unique_ptr, Electromag, Interpolator, + BoundaryCondition<1, 1>, GridLayout_t>> pusher; double mass; double dt; @@ -292,12 +332,18 @@ class APusherWithLeavingParticles : public ::testing::Test Electromag em; Interpolator interpolator; double dx = 0.1; - Box domain; - Box cells; - BoundaryCondition<1, 1> bc; - ParticleArray<1> particlesIn; - ParticleArray<1> particlesOut1; - ParticleArray<1> particlesOut2; + GridLayout_t layout; + Box domain = layout.AMRBox(); + // BoundaryCondition<1, 1> bc; + // ParticleArray<1> particlesOut1; + // ParticleArray<1> particlesOut2; + + // TestIons> ions{layout}; + UsableIons ions{layout}; + + ParticleArray& particlesIn = ions[0].domainParticles(); + + Boxing_t const boxing{layout, {layout.AMRBox()}}; }; @@ -305,34 +351,16 @@ class APusherWithLeavingParticles : public ::testing::Test TEST_F(APusherWithLeavingParticles, splitLeavingFromNonLeavingParticles) { - auto rangeIn = makeIndexRange(particlesIn); - auto inDomain = rangeIn; - - auto selector = [this](auto& particleRange) // - { - auto& box = this->cells; - return particleRange.array().partition( - [&](auto const& cell) { return PHARE::core::isIn(Point{cell}, box); }); - }; - for (decltype(nt) i = 0; i < nt; ++i) - { - auto layout = DummyLayout<1>{}; - inDomain - = pusher->move(rangeIn, rangeIn, em, mass, interpolator, layout, selector, selector); - - if (inDomain.end() != std::end(particlesIn)) - { - std::cout << inDomain.iend() << " and " << particlesIn.size() << "\n"; - break; - } - } - EXPECT_TRUE(std::none_of(inDomain.end(), std::end(particlesIn), [this](auto const& particle) { - return PHARE::core::isIn(Point{particle.iCell}, cells); - })); - EXPECT_TRUE(std::all_of(std::begin(inDomain), std::end(inDomain), [this](auto const& particle) { - return PHARE::core::isIn(Point{particle.iCell}, cells); - })); + pusher->move(ions[0], em, interpolator, boxing); + + auto& patchGhost = ions[0].patchGhostParticles(); + EXPECT_TRUE( + std::none_of(std::begin(patchGhost), std::end(patchGhost), + [this](auto const& particle) { return PHARE::core::isIn(particle, domain); })); + EXPECT_TRUE( + std::all_of(std::begin(particlesIn), std::end(particlesIn), + [this](auto const& particle) { return PHARE::core::isIn(particle, domain); })); } @@ -378,16 +406,6 @@ TEST_F(APusherWithLeavingParticles, pusherWithOrWithoutBCReturnsSameNbrOfStaying -TEST(APusherFactory, canReturnABorisPusher) -{ - auto pusher - = PusherFactory::makePusher<1, IndexRange>, Electromag, Interpolator, - BoundaryCondition<1, 1>, DummyLayout<1>>("modified_boris"); - - EXPECT_NE(nullptr, pusher); -} - - int main(int argc, char** argv) { diff --git a/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp b/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp index 0b92fd889..b837af3fd 100644 --- a/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp +++ b/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp @@ -17,7 +17,7 @@ void updater_routine(benchmark::State& state) using Electromag_t = core::UsableElectromag; using ParticleArray = PHARE_Types::ParticleArray_t; using Particle_t = ParticleArray::value_type; - using Ions = PHARE::core::UsableIons_t; + using Ions = PHARE::core::UsableIons; using IonUpdater = core::IonUpdater; using Boxing_t = PHARE::core::UpdaterSelectionBoxing; diff --git a/tools/bench/core/numerics/pusher/push_raw_use.cpp b/tools/bench/core/numerics/pusher/push_raw_use.cpp index 3a3e4e738..ae812f457 100644 --- a/tools/bench/core/numerics/pusher/push_raw_use.cpp +++ b/tools/bench/core/numerics/pusher/push_raw_use.cpp @@ -1,12 +1,14 @@ #include "push_bench.hpp" +#include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" // Does not include google benchmark as we want to see only PHARE operations/instructions/etc template void push() { - auto static constexpr opts = PHARE::SimOpts{dim, interp}; - constexpr std::uint32_t cells = 65; + auto static constexpr opts = PHARE::SimOpts{dim, interp}; + bool static constexpr copy_particle = true; + constexpr std::uint32_t cells = 65; // constexpr std::uint32_t n_parts = 1e7; using PHARE_Types = PHARE::core::PHARE_Types; @@ -16,41 +18,26 @@ void push() using Electromag_t = PHARE::core::UsableElectromag; using Ions_t = PHARE_Types::Ions_t; using ParticleArray = Ions_t::particle_array_type; - using ParticleRange = PHARE::core::IndexRange; - using BorisPusher_t = PHARE::core::BorisPusher; + using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + + using BorisPusher_t = PHARE::core::BorisPusher; GridLayout_t layout{cells}; Electromag_t em{layout}; - - ParticleArray domainParticles{layout.AMRBox()}; + Boxing_t boxing{layout, {layout.AMRBox()}}; + PHARE::core::UsableIons ions{layout}; std::stringstream ss; ss << "unsorted_particles_" << dim << ".raw"; - PHARE::core::bench::read_raw_from_file(domainParticles, ss.str()); - - // std::sort(domainParticles); - - ParticleArray tmpDomain{layout.AMRBox()}; - tmpDomain.vector() = domainParticles.vector(); - - auto rangeIn = PHARE::core::makeIndexRange(domainParticles); - auto rangeOut = PHARE::core::makeIndexRange(tmpDomain); + PHARE::core::bench::read_raw_from_file(ions[0].domainParticles(), ss.str()); BorisPusher_t pusher; pusher.setMeshAndTimeStep(layout.meshSize(), .001); Interpolator interpolator; - auto const no_op = [](auto& particleRange) { return particleRange; }; - - pusher.move( - /*ParticleRange const&*/ rangeIn, // - /*ParticleRange& */ rangeOut, // - /*Electromag const& */ em, // - /*double mass */ 1, - /*Interpolator&*/ interpolator, - /*GridLayout const&*/ layout, // - no_op, no_op); + pusher.template move(ions[0], em, interpolator, boxing); } int main(int /*argc*/, char** /*argv*/) diff --git a/tools/bench/core/numerics/pusher/pusher.cpp b/tools/bench/core/numerics/pusher/pusher.cpp index 64880205b..243c9321f 100644 --- a/tools/bench/core/numerics/pusher/pusher.cpp +++ b/tools/bench/core/numerics/pusher/pusher.cpp @@ -1,11 +1,13 @@ #include "push_bench.hpp" +#include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" template void push(benchmark::State& state) { - auto static constexpr opts = PHARE::SimOpts{dim, interp}; - constexpr std::uint32_t cells = 65; - constexpr std::uint32_t n_parts = 1e7; + auto static constexpr opts = PHARE::SimOpts{dim, interp}; + bool static constexpr copy_particle = true; + constexpr std::uint32_t cells = 65; + constexpr std::uint32_t n_parts = 1e7; using PHARE_Types = PHARE::core::PHARE_Types; using GridLayout_t = TestGridLayout; @@ -15,44 +17,38 @@ void push(benchmark::State& state) using Ions_t = PHARE_Types::Ions_t; using ParticleArray = Ions_t::particle_array_type; using Particle_t = ParticleArray::value_type; - using ParticleRange = PHARE::core::IndexRange; - using BorisPusher_t = PHARE::core::BorisPusher; + using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + + using BorisPusher_t = PHARE::core::BorisPusher; GridLayout_t layout{cells}; + Boxing_t boxing{layout, {layout.AMRBox()}}; + Electromag_t em{layout}; + PHARE::core::UsableIons ions{layout}; - ParticleArray domainParticles{layout.AMRBox()}; + auto& domainParticles = ions[0].domainParticles(); domainParticles.vector() = std::vector(n_parts, PHARE::core::bench::particle()); PHARE::core::bench::disperse(domainParticles, 0, cells - 1, 13337); - // std::sort(domainParticles); - ParticleArray tmpDomain{layout.AMRBox()}; - tmpDomain.vector() = domainParticles.vector(); - auto rangeIn = PHARE::core::makeIndexRange(domainParticles); - auto rangeOut = PHARE::core::makeIndexRange(tmpDomain); BorisPusher_t pusher; pusher.setMeshAndTimeStep(layout.meshSize(), .001); Interpolator interpolator; auto const no_op = [](auto& particleRange) { return particleRange; }; + while (state.KeepRunningBatch(1)) - { - pusher.move( - /*ParticleRange const&*/ rangeIn, // - /*ParticleRange& */ rangeOut, // - /*Electromag const& */ em, // - /*double mass */ 1, - /*Interpolator&*/ interpolator, - /*GridLayout const&*/ layout, // - no_op, no_op); - } + pusher.template move(ions[0], em, interpolator, boxing); } + + BENCHMARK_TEMPLATE(push, /*dim=*/1, /*interp=*/1)->Unit(benchmark::kMicrosecond); BENCHMARK_TEMPLATE(push, /*dim=*/1, /*interp=*/2)->Unit(benchmark::kMicrosecond); BENCHMARK_TEMPLATE(push, /*dim=*/1, /*interp=*/3)->Unit(benchmark::kMicrosecond); From 3a6856b947f84522bc3f8c15809d997392f1f264 Mon Sep 17 00:00:00 2001 From: deegan Date: Tue, 24 Mar 2026 18:02:20 +0100 Subject: [PATCH 2/3] MORE! --- src/amr/resources_manager/amr_utils.hpp | 33 +- src/amr/solvers/solver_ppc.hpp | 10 +- src/amr/utilities/box/amr_box.hpp | 7 +- src/core/data/grid/grid.hpp | 1 + src/core/data/grid/gridlayout.hpp | 6 +- src/core/data/ions/ions.hpp | 10 +- src/core/data/ndarray/ndarray_vector.hpp | 6 + src/core/data/particles/particle_array.hpp | 2 +- .../numerics/interpolator/interpolator.hpp | 1 - src/core/numerics/ion_updater/ion_updater.hpp | 189 +++------- .../ion_updater/ion_updater_impl0.hpp | 327 ++++++++++++++++++ .../ion_updater/ion_updater_impl1.hpp | 312 +++++++++++++++++ .../ion_updater/ion_updater_impl2.hpp | 297 ++++++++++++++++ .../numerics/ion_updater/ion_updater_def.hpp | 78 +++++ src/core/utilities/box/box.hpp | 83 ++++- src/core/utilities/equality.hpp | 26 ++ src/core/utilities/types.hpp | 16 + tests/core/data/field/test_field_fixtures.hpp | 103 ++++++ .../core/data/gridlayout/test_gridlayout.hpp | 39 ++- .../test_ion_population_fixtures.hpp | 200 ++++++++--- .../particles/test_particles_fixtures.hpp | 150 +++++++- .../tensorfield/test_tensorfield_fixtures.hpp | 61 +++- .../core/numerics/ion_updater/CMakeLists.txt | 19 +- .../numerics/ion_updater/test_updater.cpp | 216 ++++++------ .../numerics/ion_updater/test_updaters.cpp | 287 +++++++++++++++ tests/core/numerics/pusher/test_pusher.cpp | 275 ++++++--------- .../ion_updater/bench_ion_updater.cpp | 34 +- .../core/numerics/pusher/push_raw_use.cpp | 100 ++++-- tools/bench/core/numerics/pusher/pusher.cpp | 79 ++++- tools/bench/core/utilities/CMakeLists.txt | 8 + tools/bench/core/utilities/bench_box.hpp | 39 +++ tools/bench/core/utilities/box_iterator.cpp | 45 +++ tools/bench/core/utilities/box_span.cpp | 65 ++++ .../bench/core/utilities/box_span_N_loop.cpp | 85 +++++ 34 files changed, 2646 insertions(+), 563 deletions(-) create mode 100644 src/core/numerics/ion_updater/ion_updater/ion_updater_impl0.hpp create mode 100644 src/core/numerics/ion_updater/ion_updater/ion_updater_impl1.hpp create mode 100644 src/core/numerics/ion_updater/ion_updater/ion_updater_impl2.hpp create mode 100644 src/core/numerics/ion_updater/ion_updater_def.hpp create mode 100644 src/core/utilities/equality.hpp create mode 100644 tests/core/numerics/ion_updater/test_updaters.cpp create mode 100644 tools/bench/core/utilities/CMakeLists.txt create mode 100644 tools/bench/core/utilities/bench_box.hpp create mode 100644 tools/bench/core/utilities/box_iterator.cpp create mode 100644 tools/bench/core/utilities/box_span.cpp create mode 100644 tools/bench/core/utilities/box_span_N_loop.cpp diff --git a/src/amr/resources_manager/amr_utils.hpp b/src/amr/resources_manager/amr_utils.hpp index feaebf92d..f463ca1e5 100644 --- a/src/amr/resources_manager/amr_utils.hpp +++ b/src/amr/resources_manager/amr_utils.hpp @@ -227,29 +227,39 @@ namespace amr noDomainOverlapsOn(hierarchy, iLevel); } - // potentially to replace with SAMRAI coarse to fine boundary stuff template // fow now it gives us a box for only patch ghost layer NO_DISCARD auto makeNonLevelGhostBoxFor(SAMRAI::hier::Patch const& patch, SAMRAI::hier::PatchHierarchy const& hierarchy) { auto constexpr dimension = GridLayoutT::dimension; - auto const lvlNbr = patch.getPatchLevelNumber(); SAMRAI::hier::Box const domain = patch.getBox(); auto const domBox = phare_box_from(domain); auto const particleGhostBox = grow(domBox, GridLayoutT::nbrParticleGhosts()); - - auto const neighbors = getSameLevelNeighbors(patch, hierarchy); + auto const neighbors = getSameLevelNeighbors(patch, hierarchy); std::vector> patchGhostLayerBoxes; patchGhostLayerBoxes.reserve(neighbors.size() + 1); patchGhostLayerBoxes.emplace_back(domBox); for (auto const& neighbox : neighbors) patchGhostLayerBoxes.emplace_back( *(particleGhostBox * phare_box_from(neighbox))); - return patchGhostLayerBoxes; } + template + NO_DISCARD auto patchGhostBoxOverlaps(SAMRAI::hier::Patch const& patch, + SAMRAI::hier::PatchHierarchy const& hierarchy) + { + auto constexpr dimension = GridLayoutT::dimension; + auto const domBox = phare_box_from(patch.getBox()); + auto const particleGhostBox = grow(domBox, GridLayoutT::nbrParticleGhosts()); + return core::generate( + [&](auto& neighbox) { + return *(particleGhostBox * phare_box_from(neighbox)); + }, + getSameLevelNeighbors(patch, hierarchy)); + } + inline auto to_string(auto const& id) { std::stringstream patchID; @@ -340,7 +350,20 @@ namespace amr orMissing(ilvl); } + template + auto removeIntersections(core::Box const& box, core::Box const& remove) + { + return generate([&](auto const& rem) { return phare_box_from(rem); }, + samrai_box_from(box).removeIntersections(samrai_box_from(remove))); + } + template + auto removeIntersections(core::Box const& box, + core::Box const& remove) + { + return generate([&](auto const& rem) { return as_unsigned_phare_box(rem); }, + samrai_box_from(box).removeIntersections(samrai_box_from(remove))); + } } // namespace amr } // namespace PHARE diff --git a/src/amr/solvers/solver_ppc.hpp b/src/amr/solvers/solver_ppc.hpp index 2779052fe..03b2979bd 100644 --- a/src/amr/solvers/solver_ppc.hpp +++ b/src/amr/solvers/solver_ppc.hpp @@ -50,7 +50,7 @@ class SolverPPC : public ISolver using Faraday_t = ModelViews_t::Faraday_t; using Ampere_t = ModelViews_t::Ampere_t; using Ohm_t = ModelViews_t::Ohm_t; - using IonUpdater_t = PHARE::core::IonUpdater; + using IonUpdater_t = PHARE::core::IonUpdater; Electromag electromagPred_{"EMPred"}; Electromag electromagAvg_{"EMAvg"}; @@ -184,7 +184,7 @@ class SolverPPC : public ISolver if (auto [it, suc] = levelBoxing.try_emplace( amr::to_string(patch->getGlobalId()), Boxing_t{amr::layoutFromPatch(*patch), - amr::makeNonLevelGhostBoxFor(*patch, hierarchy)}); + amr::patchGhostBoxOverlaps(*patch, hierarchy)}); !suc) throw std::runtime_error("boxing map insertion failure"); } @@ -198,7 +198,7 @@ class SolverPPC : public ISolver } - using Boxing_t = core::UpdaterSelectionBoxing; + using Boxing_t = core::UpdaterSelectionBoxing; std::unordered_map> boxing; @@ -559,8 +559,8 @@ void SolverPPC::moveIons_(level_t& level, ModelViews_t& auto dt = newTime - currentTime; for (auto& state : views) ionUpdater_.updatePopulations( - state.ions, state.electromagAvg, - levelBoxing.at(amr::to_string(state.patch->getGlobalId())), dt, mode); + mode, state.ions, state.electromagAvg, + levelBoxing.at(amr::to_string(state.patch->getGlobalId())), dt); } catch (core::DictionaryException const& ex) { diff --git a/src/amr/utilities/box/amr_box.hpp b/src/amr/utilities/box/amr_box.hpp index 6badca945..bffd4d65c 100644 --- a/src/amr/utilities/box/amr_box.hpp +++ b/src/amr/utilities/box/amr_box.hpp @@ -1,14 +1,11 @@ #ifndef PHARE_AMR_UTILITIES_BOX_BOX_HPP #define PHARE_AMR_UTILITIES_BOX_BOX_HPP - +#include "core/def.hpp" #include "core/def/phare_mpi.hpp" // IWYU pragma: keep - - -#include "SAMRAI/hier/Box.h" #include "core/utilities/box/box.hpp" -#include "core/def.hpp" +#include "SAMRAI/hier/Box.h" namespace PHARE::amr { diff --git a/src/core/data/grid/grid.hpp b/src/core/data/grid/grid.hpp index e12da2634..80aceaca7 100644 --- a/src/core/data/grid/grid.hpp +++ b/src/core/data/grid/grid.hpp @@ -107,6 +107,7 @@ class Grid : public NdArrayImpl // returns view when getting address of this object, could be misleading, but convenient NO_DISCARD auto operator&() { return &field_; } NO_DISCARD auto operator&() const { return &field_; } + NO_DISCARD auto operator*() const { return field_; } private: std::string name_{"No Name"}; diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index 655d1f38b..f6dbc540c 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -99,7 +99,7 @@ namespace core static constexpr std::size_t interp_order = GridLayoutImpl::interp_order; using This = GridLayout; using implT = GridLayoutImpl; - + using AMRBox_t = Box; /** * @brief Constructor of a GridLayout @@ -109,8 +109,8 @@ namespace core */ GridLayout(std::array const& meshSize, std::array const& nbrCells, - Point const& origin, - Box AMRBox = Box{}, int level_number = 0) + Point const& origin, AMRBox_t const AMRBox = {}, + int level_number = 0) : meshSize_{meshSize} , origin_{origin} , nbrPhysicalCells_{nbrCells} diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 49c8f0ba0..c2c069143 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -27,11 +27,11 @@ namespace core { public: using value_type = IonPopulation; - using field_type = typename IonPopulation::field_type; - using vecfield_type = typename IonPopulation::vecfield_type; - using Float = typename field_type::type; - using tensorfield_type = typename IonPopulation::tensorfield_type; - using particle_array_type = typename IonPopulation::particle_array_type; + using field_type = IonPopulation::field_type; + using vecfield_type = IonPopulation::vecfield_type; + using Float = field_type::type; + using tensorfield_type = IonPopulation::tensorfield_type; + using particle_array_type = IonPopulation::particle_array_type; using gridlayout_type = GridLayout; static constexpr auto dimension = GridLayout::dimension; diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index 5c8b370c3..a9295c2b7 100644 --- a/src/core/data/ndarray/ndarray_vector.hpp +++ b/src/core/data/ndarray/ndarray_vector.hpp @@ -314,6 +314,12 @@ class NdArrayVector return MaskedView{*this, std::forward(mask)}; } + void zero() { fill(0); } + auto& fill(DataType const& v) + { + std::fill(begin(), end(), v); + return *this; + } NO_DISCARD auto& vector() { return data_; } NO_DISCARD auto& vector() const { return data_; } diff --git a/src/core/data/particles/particle_array.hpp b/src/core/data/particles/particle_array.hpp index bc0184c50..9d57e688a 100644 --- a/src/core/data/particles/particle_array.hpp +++ b/src/core/data/particles/particle_array.hpp @@ -42,7 +42,7 @@ class ParticleArray public: - ParticleArray(box_t box) + ParticleArray(box_t box = {}) : box_{box} , cellMap_{box_} { diff --git a/src/core/numerics/interpolator/interpolator.hpp b/src/core/numerics/interpolator/interpolator.hpp index 38ad8c416..2f34e0e95 100644 --- a/src/core/numerics/interpolator/interpolator.hpp +++ b/src/core/numerics/interpolator/interpolator.hpp @@ -459,7 +459,6 @@ class Interpolator : private Weighter } - /**\brief interpolate electromagnetic fields on all particles in the range * * For each particle : diff --git a/src/core/numerics/ion_updater/ion_updater.hpp b/src/core/numerics/ion_updater/ion_updater.hpp index 0d7441d41..de4fe71b6 100644 --- a/src/core/numerics/ion_updater/ion_updater.hpp +++ b/src/core/numerics/ion_updater/ion_updater.hpp @@ -3,181 +3,94 @@ #include "core/logger.hpp" -#include "core/utilities/box/box.hpp" -#include "core/numerics/pusher/boris.hpp" -#include "core/numerics/moments/moments.hpp" -#include "core/numerics/interpolator/interpolator.hpp" -#include "core/numerics/boundary_condition/boundary_condition.hpp" -#include "initializer/data_provider.hpp" +#include "core/numerics/ion_updater/ion_updater/ion_updater_impl0.hpp" +#include "core/numerics/ion_updater/ion_updater/ion_updater_impl1.hpp" +#include "core/numerics/ion_updater/ion_updater/ion_updater_impl2.hpp" +#include -#include +#ifndef PHARE_UPDATER_IMPL +#define PHARE_UPDATER_IMPL 0 +#endif -namespace PHARE::core +namespace PHARE::core::detail { -enum class UpdaterMode { domain_only = 1, all = 2 }; - - -template -class IonUpdater +template +auto reset_updater_impl(Impl_t& impl) -> decltype(impl.reset(), void()) { -public: - static constexpr auto dimension = GridLayout::dimension; - static constexpr auto interp_order = GridLayout::interp_order; - - using Box = PHARE::core::Box; - using Interpolator = PHARE::core::Interpolator; - using ParticleArray = Ions::particle_array_type; - using Particle_t = ParticleArray::Particle_t; - using PartIterator = ParticleArray::iterator; + impl.reset(); +} - using BoundaryCondition = PHARE::core::BoundaryCondition; - using Pusher = BorisPusher; +auto reset_updater_impl(auto&&...) {} -private: - Pusher pusher_; - Interpolator interpolator_; +} // namespace PHARE::core::detail -public: - IonUpdater() = default; - IonUpdater(auto const& /*dict*/) {} +namespace PHARE::core +{ - template - void updatePopulations(Ions& ions, Electromag const& em, Boxing_t const& boxing, double dt, - UpdaterMode = UpdaterMode::all); +template +class IonUpdaterProxy +{ +public: + using Interpolator_t = Impl_t::Interpolator_t; - void updateIons(Ions& ions); + IonUpdaterProxy(auto&&... args) + : impl{args...} + { + } + void updatePopulations(auto&&... args); - void reset() { /* noop */ } + void updateIons(auto& ions); + void reset() { detail::reset_updater_impl(impl); } private: - template - void updateCopy(Ions& ions, Electromag const& em, Boxing_t const& boxing); - - template - void updateInplace(Ions& ions, Electromag const& em, Boxing_t const& boxing); + Impl_t impl; }; - -template -template -void IonUpdater::updatePopulations(Ions& ions, Electromag const& em, - Boxing_t const& boxing, double dt, - UpdaterMode mode) +template +void IonUpdaterProxy::updatePopulations(auto&&... args) { PHARE_LOG_SCOPE(3, "IonUpdater::updatePopulations"); - resetMoments(ions); - pusher_.setMeshAndTimeStep(boxing.layout.meshSize(), dt); - - if (mode == UpdaterMode::domain_only) - { - updateCopy(ions, em, boxing); - } - else - { - updateInplace(ions, em, boxing); - } + impl.updatePopulations(args...); } - -template -void IonUpdater::updateIons(Ions& ions) +template +void IonUpdaterProxy::updateIons(auto& ions) { ions.computeChargeDensity(); ions.computeBulkVelocity(); } -// this is to detach how we partition particles from the updater directly -template -struct UpdaterSelectionBoxing +template +struct IonUpdaterImplResolver { - auto constexpr static partGhostWidth = GridLayout::nbrParticleGhosts(); - using GridLayout_t = GridLayout; - using Particle_t = IonUpdater_t::Particle_t; - using Box_t = IonUpdater_t::Box; - using ParticleArray_t = IonUpdater_t::ParticleArray; - using ParticleRange = IndexRange; - using Selector_t = std::function; - - GridLayout_t const layout; - std::vector const nonLevelGhostBox; - Box_t const domainBox = layout.AMRBox(); - Box_t const ghostBox = grow(domainBox, partGhostWidth); - - Selector_t const noop = [](auto& particleRange) { return particleRange; }; - - bool isInGhostBox(Particle_t& p) const { return isIn(p, ghostBox); }; - bool isInDomainBox(Particle_t& p) const { return isIn(p, domainBox); }; - bool isInNonLevelGhostBox(auto const& icell) const { return isIn(icell, nonLevelGhostBox); }; - - // lambda copy captures to detach from above references in case of class copy construct - Selector_t const inDomainBox = [domainBox = domainBox](auto& particleRange) { - return particleRange.array().partition( - particleRange, [&](auto const& cell) { return core::isIn(cell, domainBox); }); - }; - - Selector_t const inGhostBox = [ghostBox = ghostBox](auto& particleRange) { - return particleRange.array().partition( - particleRange, [&](auto const& cell) { return isIn(cell, ghostBox); }); - }; - - Selector_t const inNonLevelGhostBox - = [nonLevelGhostBox = nonLevelGhostBox](auto& particleRange) { - return particleRange.array().partition(particleRange, [&](auto const& cell) { - return isIn(Point{cell}, nonLevelGhostBox); - }); - }; - - Selector_t const inGhostLayer - = [ghostBox = ghostBox, domainBox = domainBox](auto& particleRange) { - return particleRange.array().partition(particleRange, [&](auto const& cell) { - return isIn(cell, ghostBox) and !isIn(cell, domainBox); - }); - }; - - Selector_t const outsideGhostBox = [ghostBox = ghostBox](auto& particleRange) { - return particleRange.array().partition( - particleRange, [&](auto const& cell) { return !isIn(cell, ghostBox); }); - }; -}; - - -template -template -void IonUpdater::updateCopy(Ions& ions, Electromag const& em, - Boxing_t const& boxing) -{ - bool constexpr copy_particle = true; - - PHARE_LOG_SCOPE(3, "IonUpdater::updateCopy"); - - for (auto& pop : ions) - pusher_.template move(pop, em, interpolator_, boxing); -} - + auto static constexpr updater_impl() + { + if constexpr (PHARE_UPDATER_IMPL == 0) + return static_cast*>(0); + if constexpr (PHARE_UPDATER_IMPL == 1) + return static_cast*>(0); + else + throw std::runtime_error("NO IMPL FOR VALUE"); + } -template -template -void IonUpdater::updateInplace(Ions& ions, Electromag const& em, - Boxing_t const& boxing) -{ - bool constexpr copy_particle = false; + using value_type = std::decay_t; +}; - PHARE_LOG_SCOPE(3, "IonUpdater::updateInplace"); +template +using IonUpdaterImpl_t = IonUpdaterImplResolver::value_type; - for (auto& pop : ions) - pusher_.template move(pop, em, interpolator_, boxing); -} +template +using IonUpdater = IonUpdaterProxy>; diff --git a/src/core/numerics/ion_updater/ion_updater/ion_updater_impl0.hpp b/src/core/numerics/ion_updater/ion_updater/ion_updater_impl0.hpp new file mode 100644 index 000000000..161241f06 --- /dev/null +++ b/src/core/numerics/ion_updater/ion_updater/ion_updater_impl0.hpp @@ -0,0 +1,327 @@ +#ifndef PHARE_CORE_NUMERICS_ION_UPDATER0_HPP +#define PHARE_CORE_NUMERICS_ION_UPDATER0_HPP + +#include "core/errors.hpp" +#include "core/logger.hpp" +#include "core/utilities/box/box.hpp" +#include "core/utilities/range/range.hpp" +#include "core/numerics/pusher/boris.hpp" +#include "core/numerics/moments/moments.hpp" +#include "core/numerics/interpolator/interpolator.hpp" +#include "core/numerics/ion_updater/ion_updater_def.hpp" + + +#include +#include +#include +#include + + +namespace PHARE::core +{ + +template +class IonUpdater0 +{ + using GridLayout_t = Ions::gridlayout_type; + +public: + static constexpr auto dimension = GridLayout_t::dimension; + static constexpr auto interp_order = GridLayout_t::interp_order; + using Interpolator_t = Interpolator; + using ParticleArray = Ions::particle_array_type; + +public: + IonUpdater0(auto&&...) {} + + void updatePopulations( // + UpdaterMode const, Ions& ions, auto const& em, auto const& boxing, double const dt); + + auto static make_pusher(GridLayout_t const& layout, auto const dt, auto const mass) + { + return Pusher{dt, mass, layout}; + } + + void reset() + { + tmp_particles_ = std::move(ParticleArray{}); // clear memory + } + + class Pusher; + +private: + void updateAndDepositDomain_(auto& ions, auto const& em, auto const& boxing, double const dt); + + + void updateAndDepositAll_(auto& ions, auto const& em, auto const& boxing, double const dt); + + // dealloced on regridding/load balancing coarsest + Interpolator_t interpolator; + ParticleArray tmp_particles_{}; +}; + +template +class IonUpdater0::Pusher +{ +public: + auto move(auto const& rangeIn, auto& rangeOut, auto const& emFields, auto& interpolator, + auto const& firstSelector, auto const& secondSelector) const; + + + double const dt; + double const mass; + GridLayout_t const layout; + + double const dto2m = 0.5 * dt / mass; + std::array const halfDtOverDl + = for_N_make_array([&](auto i) { return 0.5 * dt / layout.meshSize()[i]; }); + +private: + void prePushStep_(auto const& rangeIn, auto& rangeOut) const; + void postPushStep_(auto& range, std::size_t const idx) const; +}; + + +template +void IonUpdater0::updatePopulations(UpdaterMode mode, Ions& ions, auto const& em, + auto const& boxing, double const dt) +{ + PHARE_LOG_SCOPE(3, "IonUpdater0::updatePopulations"); + + resetMoments(ions); + + if (mode == UpdaterMode::domain_only) + updateAndDepositDomain_(ions, em, boxing, dt); + + else + updateAndDepositAll_(ions, em, boxing, dt); +} + + + +/** + * @brief IonUpdater0::updateAndDepositDomain_ + evolves moments from time n to n+1 without updating particles, which stay at time n + */ +template +void IonUpdater0::updateAndDepositDomain_( // + auto& ions, auto const& em, auto const& boxing, double const dt) +{ + PHARE_LOG_SCOPE(3, "IonUpdater0::updateAndDepositDomain_"); + + auto const& layout = boxing.layout; + + for (auto& pop : ions) + { + auto const pusher = make_pusher(layout, dt, pop.mass()); + + auto& domain = (tmp_particles_ = pop.domainParticles()); // make local copy + + // first push all domain particles twice + // accumulate those inNonLevelGhostBox + auto outRange = makeIndexRange(domain); + auto allowed = outRange = pusher.move( // + outRange, outRange, em, interpolator, boxing.noop, boxing.inNonLevelGhostBox); + + interpolator(allowed, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + + + // push those in the ghostArea (i.e. stop pushing if they're not out of it) + // deposit moments on those which leave to go inDomainBox + + auto pushAndAccumulateGhosts = [&](auto const& inputArray) { + tmp_particles_ = inputArray; // work on local copy + + auto outRange = makeIndexRange(tmp_particles_); + + auto enteredInDomain = pusher.move( // + outRange, outRange, em, interpolator, boxing.inGhostBox, boxing.inDomainBox); + + interpolator(enteredInDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(), + layout); + }; + + // After this function is done domain particles overlaping ghost layers of neighbor patches + // are sent to these neighbor's patchghost particle array. + // After being pushed, some patch ghost particles may enter the domain. These need to be + // copied into the domain array so they are transfered to the neighbor patch + // ghost array and contribute to moments there too. + // On the contrary level ghost particles entering the domain here do not need to be copied + // since they contribute to nodes that are not shared with neighbor patches an since + // level border nodes will receive contributions from levelghost old and new particles + + if (pop.levelGhostParticles().size()) + pushAndAccumulateGhosts(pop.levelGhostParticles()); + } +} + + +/** + * @brief IonUpdater0::updateAndDepositDomain_ + evolves moments and particles from time n to n+1 + */ +template +void IonUpdater0::updateAndDepositAll_( // + auto& ions, auto const& em, auto const& boxing, double const dt) +{ + PHARE_LOG_SCOPE(3, "IonUpdater0::updateAndDepositAll_"); + + auto const& layout = boxing.layout; + + // push domain particles, erase from array those leaving domain + // push level ghost particles that are in ghost area (==ghost box without domain) + // copy ghost particles out of ghost area that are in domain, in particle array + // finally all particles in non level ghost box are to be interpolated on mesh. + for (auto& pop : ions) + { + auto const pusher = make_pusher(layout, dt, pop.mass()); + auto& domainParticles = pop.domainParticles(); + auto domainPartRange = makeIndexRange(domainParticles); + + auto inDomain = pusher.move( // + domainPartRange, domainPartRange, em, interpolator, boxing.noop, boxing.inDomainBox); + + auto now_ghosts = makeRange(domainParticles, inDomain.iend(), domainParticles.size()); + auto const not_level_ghosts = boxing.inNonLevelGhostBox(now_ghosts); + + // copy out new patch ghosts + auto& patchGhost = pop.patchGhostParticles(); + patchGhost.reserve(patchGhost.size() + not_level_ghosts.size()); + std::copy(not_level_ghosts.begin(), not_level_ghosts.end(), std::back_inserter(patchGhost)); + + domainParticles.erase(now_ghosts); // drop all ghosts + + if (pop.levelGhostParticles().size()) + { + auto particleRange = makeIndexRange(pop.levelGhostParticles()); + auto inGhostLayerRange = pusher.move(particleRange, particleRange, em, interpolator, + boxing.inGhostBox, boxing.inGhostLayer); + + auto& particleArray = particleRange.array(); + particleArray.export_particles( + domainParticles, [&](auto const& cell) { return isIn(cell, boxing.domainBox); }); + + particleArray.erase( + makeRange(particleArray, inGhostLayerRange.iend(), particleArray.size())); + } + + interpolator( // + domainParticles, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + interpolator( // + patchGhost, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + } +} + + + +template +auto IonUpdater0::Pusher::move(auto const& rangeIn, auto& rangeOut, auto const& emFields, + auto& interpolator, auto const& firstSelector, + auto const& secondSelector) const +{ + PHARE_LOG_SCOPE(3, "Boris::move_no_bc"); + + + // push the particles of half a step + // 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 + prePushStep_(rangeIn, rangeOut); + rangeOut = firstSelector(rangeOut); + + for (auto idx = rangeOut.ibegin(); idx < rangeOut.iend(); ++idx) + { + auto& currPart = rangeOut.array()[idx]; + + // get electromagnetic fields interpolated on the particles of rangeOut stop at newEnd. + // get the particle velocity from t=n to t=n+1 + auto const& local_em = interpolator(currPart, emFields, layout); + boris::accelerate(currPart, local_em, 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 + try + { + postPushStep_(rangeOut, idx); + } + catch (DictionaryException const& bex) + { + auto ex = bex; + auto const& [e, b] = local_em; + for (std::uint16_t i = 0; i < 3; ++i) + ex("E_" + std::to_string(i), std::to_string(e[i])); + for (std::uint16_t i = 0; i < 3; ++i) + ex("B_" + std::to_string(i), std::to_string(b[i])); + ex("level", std::to_string(layout.levelNumber())); + throw ex; + } + } + + return secondSelector(rangeOut); +} + + + +template +void IonUpdater0::Pusher::prePushStep_(auto const& rangeIn, auto& rangeOut) const +{ + auto& inParticles = rangeIn.array(); + auto& outParticles = rangeOut.array(); + for (auto inIdx = rangeIn.ibegin(), outIdx = rangeOut.ibegin(); inIdx < rangeIn.iend(); + ++inIdx, ++outIdx) + { + // in the first push, this is the first time + // we push to rangeOut, which contains crap + // the push will only touch the particle position + // but the next step being the acceleration of + // rangeOut, we need to copy rangeIn weights, charge + // and velocity. This is done here although + // not strictly speaking this function's business + // to take advantage that we're already looping + // over rangeIn particles. + + outParticles[outIdx].charge = inParticles[inIdx].charge; + outParticles[outIdx].weight = inParticles[inIdx].weight; + outParticles[outIdx].v = inParticles[inIdx].v; + + try + { + auto const newCell = boris::advance(outParticles[outIdx], halfDtOverDl); + if (newCell != inParticles[inIdx].iCell) + outParticles.change_icell(newCell, outIdx); + } + catch (boris::MoveTwoCellException const& e) + { + std::stringstream ss; + ss << "PrePush Particle moved 2 cells with delta/vel: "; + ss << e.delta << "/" << e.vel << std::endl; + DictionaryException ex{"cause", ss.str()}; + throw ex; + } + } +} + +template +void IonUpdater0::Pusher::postPushStep_(auto& range, std::size_t const idx) const +{ + try + { + auto& particles = range.array(); + auto const newCell = boris::advance(particles[idx], halfDtOverDl); + if (newCell != particles[idx].iCell) + particles.change_icell(newCell, idx); + } + catch (boris::MoveTwoCellException const& e) + { + std::stringstream ss; + ss << "PostPush Particle moved 2 cells with delta/vel: "; + ss << e.delta << "/" << e.vel << std::endl; + throw DictionaryException{}("cause", ss.str()); + } +} + + +} // namespace PHARE::core + + +#endif // ION_UPDATER_HPP diff --git a/src/core/numerics/ion_updater/ion_updater/ion_updater_impl1.hpp b/src/core/numerics/ion_updater/ion_updater/ion_updater_impl1.hpp new file mode 100644 index 000000000..c318709a2 --- /dev/null +++ b/src/core/numerics/ion_updater/ion_updater/ion_updater_impl1.hpp @@ -0,0 +1,312 @@ +#ifndef PHARE_CORE_NUMERICS_ION_UPDATER1_HPP +#define PHARE_CORE_NUMERICS_ION_UPDATER1_HPP + +#include "core/errors.hpp" +#include "core/logger.hpp" +#include "core/utilities/types.hpp" +#include "core/utilities/box/box.hpp" +#include "core/numerics/pusher/boris.hpp" +#include "core/utilities/range/range.hpp" +#include "core/numerics/moments/moments.hpp" +#include "core/numerics/interpolator/interpolator.hpp" +#include "core/numerics/ion_updater/ion_updater_def.hpp" +#include "core/numerics/boundary_condition/boundary_condition.hpp" + + +namespace PHARE::core +{ + +template +class IonUpdater1 +{ + using GridLayout_t = Ions::gridlayout_type; + +public: + static constexpr auto dimension = GridLayout_t::dimension; + static constexpr auto interp_order = GridLayout_t::interp_order; + + using Box = PHARE::core::Box; + using Interpolator_t = Interpolator; + using BoundaryCondition_t = BoundaryCondition; + +public: + IonUpdater1(auto&&... /*dict ?*/) {} + + void updatePopulations(UpdaterMode const mode, auto& ions, auto&&... args); + + auto static make_pusher(GridLayout_t const& layout, auto const dt, auto const mass) + { + return Pusher{dt, mass, layout}; + } + + class Pusher; + +private: + void updateAndDepositDomain_(auto& ions, auto const& em, auto const& boxing, double const dt); + + void updateAndDepositAll_(auto& ions, auto const& em, auto const& boxing, double const dt); + + Interpolator_t interpolator; // default +}; + + +template +class IonUpdater1::Pusher +{ +public: + void move_particle(auto& particle, auto const& em, auto const& layout, + auto& interpolator) const; + + template + void move_interpolate_and_sort(auto& pop, auto&&... args) const + { + move_domain(pop, args...); + move_level_ghost(pop, args...); // might modify domain + } + + double const dt; + double const mass; + GridLayout_t const layout; + + double const dto2m = 0.5 * dt / mass; + std::array const halfDtOverDl + = for_N_make_array([&](auto i) { return 0.5 * dt / layout.meshSize()[i]; }); + +private: + template + void move_domain(auto& pop, auto const& em, auto const& boxing, auto& interpolator) const; + + template + void move_level_ghost(auto& pop, auto const& em, auto const& boxing, auto& interpolator) const; +}; + + +template +void IonUpdater1::updatePopulations(UpdaterMode const mode, auto& ions, auto&&... args) +{ + PHARE_LOG_SCOPE(3, "IonUpdater1::updatePopulations"); + + resetMoments(ions); + + if (mode == UpdaterMode::domain_only) + updateAndDepositDomain_(ions, args...); + else + updateAndDepositAll_(ions, args...); +} + + + +/** + * @brief IonUpdater1::updateAndDepositDomain_ + evolves moments from time n to n+1 without updating particles, which stay at time n + */ +template +void IonUpdater1::updateAndDepositDomain_(auto& ions, auto const& em, auto const& boxing, + double const dt) +{ + bool constexpr copy_particle = true; + + PHARE_LOG_SCOPE(3, "IonUpdater1::updateAndDepositDomain_"); + + for (auto& pop : ions) + Pusher{dt, pop.mass(), boxing.layout}.template // + move_interpolate_and_sort(pop, em, boxing, interpolator); +} + + + +template +void IonUpdater1::updateAndDepositAll_( // + auto& ions, auto const& em, auto const& boxing, double const dt) +{ + bool constexpr copy_particle = false; + + PHARE_LOG_SCOPE(3, "IonUpdater1::updateAndDepositAll_"); + + for (auto& pop : ions) + Pusher{dt, pop.mass(), boxing.layout}.template // + move_interpolate_and_sort(pop, em, boxing, interpolator); +} + + +template +void IonUpdater1::Pusher::move_particle( // + auto& particle, auto const& em, auto const& layout, auto& interpolator) const +{ + try + { + particle.iCell = boris::advance(particle, halfDtOverDl); + } + catch (boris::MoveTwoCellException const& e) + { + std::stringstream ss; + ss << "PrePush Particle moved 2 cells with delta/vel: "; + ss << e.delta << "/" << e.vel << std::endl; + throw DictionaryException{ss.str()}; + } + + auto const& local_em = interpolator(particle, em, layout); + try + { + boris::accelerate(particle, local_em, dto2m); + particle.iCell = boris::advance(particle, halfDtOverDl); + } + catch (boris::MoveTwoCellException const& e) + { + std::stringstream ss; + ss << "PostPush Particle moved 2 cells with delta/vel: "; + ss << e.delta << "/" << e.vel << std::endl; + DictionaryException ex{ss.str()}; + + auto const& [E, B] = local_em; + for (std::uint16_t i = 0; i < 3; ++i) + ex("E_" + std::to_string(i), std::to_string(E[i])); + for (std::uint16_t i = 0; i < 3; ++i) + ex("B_" + std::to_string(i), std::to_string(B[i])); + ex("level", std::to_string(layout.levelNumber())); + throw ex; + } +}; + + + + +template +template +void IonUpdater1::Pusher::move_domain( // + auto& pop, auto const& em, auto const& boxing, auto& interpolator) const +{ + using Particle_t = std::decay_t; + using ParticleLoop_t = std::conditional_t; + + + auto const& layout = boxing.layout; + + auto& domain = pop.domainParticles(); + auto& patch_ghost = pop.patchGhostParticles(); + + for (std::size_t i = 0; i < domain.size(); ++i) // size might change on iteration! + { + ParticleLoop_t particle = domain[i]; + auto const oldCell = particle.iCell; + + move_particle(particle, em, layout, interpolator); + + if (oldCell == particle.iCell) + { + interpolator.particleToMesh( // + particle, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + + continue; + } + + bool const isInDomainBox = boxing.isInDomainBox(particle); + bool const isInNonLevelGhostBox + = isInDomainBox || boxing.isInNonLevelGhostBox(particle.iCell); + bool const should_interpolate = isInNonLevelGhostBox; + + if (!should_interpolate) + { + if constexpr (!copy_particle) + { + domain.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + } + continue; + } + + interpolator.particleToMesh( // + particle, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + + if constexpr (!copy_particle) + domain.change_icell(particle, oldCell, i); + } + + // move to patch_ghost + if constexpr (!copy_particle) + { + auto range = makeIndexRange(domain); + range = domain.partition( + range, [&](auto const& cell) { return core::isIn(cell, boxing.domainBox); }); + + auto const not_in_domain = makeRange(domain, range.iend(), domain.size()); + auto const is_patch_ghost = domain.partition( + not_in_domain, [&](auto const& cell) { return boxing.isInNonLevelGhostBox(cell); }); + + patch_ghost.reserve(patch_ghost.size() + is_patch_ghost.size()); + std::copy(is_patch_ghost.begin(), is_patch_ghost.end(), std::back_inserter(patch_ghost)); + domain.erase(not_in_domain); + } + + // PHARE_DEBUG_DO({ + // for (auto const& p : domain) + // if (!isIn(p, boxing.domainBox)) + // throw std::runtime_error("invalid domain"); + // for (auto const& p : patch_ghost) + // if (isIn(p, boxing.domainBox)) + // throw std::runtime_error("invalid patch ghost"); + // }) +} + + + +template +template +void IonUpdater1::Pusher::move_level_ghost( // + auto& pop, auto const& em, auto const& boxing, auto& interpolator) const +{ + using Particle_t = std::decay_t; + using ParticleLoop_t = std::conditional_t; + + auto const& layout = boxing.layout; + + auto& domain = pop.domainParticles(); + auto& level_ghost = pop.levelGhostParticles(); + + for (std::size_t i = 0; i < level_ghost.size(); ++i) // size might change on iteration! + { + ParticleLoop_t particle = level_ghost[i]; + auto const oldCell = particle.iCell; + + move_particle(particle, em, layout, interpolator); + + if (oldCell == particle.iCell) + continue; + + bool const isInDomainBox = boxing.isInDomainBox(particle); + bool const should_interpolate = isInDomainBox; + + if (should_interpolate) + interpolator.particleToMesh( // + particle, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + + if constexpr (!copy_particle) + { + if (isInDomainBox) + { + domain.push_back(particle); + level_ghost.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + continue; + } + + bool const isInNonLevelGhostBox + = isInDomainBox || boxing.isInNonLevelGhostBox(particle.iCell); + bool const isInLevelGhostBox = !isInNonLevelGhostBox; + + if (isInLevelGhostBox) + level_ghost.change_icell(particle, oldCell, i); + else + { + level_ghost.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + } + } + } +} + + +} // namespace PHARE::core + + +#endif // ION_UPDATER_HPP diff --git a/src/core/numerics/ion_updater/ion_updater/ion_updater_impl2.hpp b/src/core/numerics/ion_updater/ion_updater/ion_updater_impl2.hpp new file mode 100644 index 000000000..1e901bdb5 --- /dev/null +++ b/src/core/numerics/ion_updater/ion_updater/ion_updater_impl2.hpp @@ -0,0 +1,297 @@ +#ifndef PHARE_CORE_NUMERICS_ION_UPDATER2_HPP +#define PHARE_CORE_NUMERICS_ION_UPDATER2_HPP + +#include "core/errors.hpp" +#include "core/logger.hpp" +#include "core/utilities/types.hpp" +#include "core/utilities/box/box.hpp" +#include "core/numerics/pusher/boris.hpp" +#include "core/utilities/range/range.hpp" +#include "core/numerics/moments/moments.hpp" +#include "core/numerics/interpolator/interpolator.hpp" +#include "core/numerics/ion_updater/ion_updater_def.hpp" +#include "core/numerics/boundary_condition/boundary_condition.hpp" + + +namespace PHARE::core +{ + +template +class IonUpdater2 +{ + using GridLayout_t = Ions::gridlayout_type; + +public: + static constexpr auto dimension = GridLayout_t::dimension; + static constexpr auto interp_order = GridLayout_t::interp_order; + + using Box = PHARE::core::Box; + using Interpolator_t = Interpolator; + using BoundaryCondition_t = BoundaryCondition; + +public: + IonUpdater2(auto&&... /*dict ?*/) {} + + void updatePopulations(UpdaterMode const mode, auto& ions, auto&&... args); + + auto static make_pusher(GridLayout_t const& layout, auto const dt, auto const mass) + { + return Pusher{dt, mass, layout}; + } + + class Pusher; + +private: + void updateAndDepositDomain_(auto& ions, auto const& em, auto const& boxing, double const dt); + + void updateAndDepositAll_(auto& ions, auto const& em, auto const& boxing, double const dt); + + Interpolator_t interpolator; // default +}; + + +template +class IonUpdater2::Pusher +{ +public: + void move_particle(auto& particle, auto const& em, auto const& layout, + auto& interpolator) const; + + template + void move_interpolate_and_sort(auto& pop, auto&&... args) const + { + move_domain(pop, args...); + move_level_ghost(pop, args...); // might modify domain + } + + double const dt; + double const mass; + GridLayout_t const layout; + + double const dto2m = 0.5 * dt / mass; + std::array const halfDtOverDl + = for_N_make_array([&](auto i) { return 0.5 * dt / layout.meshSize()[i]; }); + +private: + template + void move_domain(auto& pop, auto const& em, auto const& boxing, auto& interpolator) const; + + template + void move_level_ghost(auto& pop, auto const& em, auto const& boxing, auto& interpolator) const; +}; + + +template +void IonUpdater2::updatePopulations(UpdaterMode const mode, auto& ions, auto&&... args) +{ + PHARE_LOG_SCOPE(3, "IonUpdater2::updatePopulations"); + + resetMoments(ions); + + if (mode == UpdaterMode::domain_only) + updateAndDepositDomain_(ions, args...); + else + updateAndDepositAll_(ions, args...); +} + + + +/** + * @brief IonUpdater2::updateAndDepositDomain_ + evolves moments from time n to n+1 without updating particles, which stay at time n + */ +template +void IonUpdater2::updateAndDepositDomain_(auto& ions, auto const& em, auto const& boxing, + double const dt) +{ + bool constexpr copy_particle = true; + + PHARE_LOG_SCOPE(3, "IonUpdater2::updateAndDepositDomain_"); + + for (auto& pop : ions) + Pusher{dt, pop.mass(), boxing.layout}.template // + move_interpolate_and_sort(pop, em, boxing, interpolator); +} + + + +template +void IonUpdater2::updateAndDepositAll_( // + auto& ions, auto const& em, auto const& boxing, double const dt) +{ + bool constexpr copy_particle = false; + + PHARE_LOG_SCOPE(3, "IonUpdater2::updateAndDepositAll_"); + + for (auto& pop : ions) + Pusher{dt, pop.mass(), boxing.layout}.template // + move_interpolate_and_sort(pop, em, boxing, interpolator); +} + + +template +void IonUpdater2::Pusher::move_particle( // + auto& particle, auto const& em, auto const& layout, auto& interpolator) const +{ + try + { + particle.iCell = boris::advance(particle, halfDtOverDl); + } + catch (boris::MoveTwoCellException const& e) + { + std::stringstream ss; + ss << "PrePush Particle moved 2 cells with delta/vel: "; + ss << e.delta << "/" << e.vel << std::endl; + throw DictionaryException{ss.str()}; + } + + auto const& local_em = interpolator(particle, em, layout); + try + { + boris::accelerate(particle, local_em, dto2m); + particle.iCell = boris::advance(particle, halfDtOverDl); + } + catch (boris::MoveTwoCellException const& e) + { + std::stringstream ss; + ss << "PostPush Particle moved 2 cells with delta/vel: "; + ss << e.delta << "/" << e.vel << std::endl; + DictionaryException ex{ss.str()}; + + auto const& [E, B] = local_em; + for (std::uint16_t i = 0; i < 3; ++i) + ex("E_" + std::to_string(i), std::to_string(E[i])); + for (std::uint16_t i = 0; i < 3; ++i) + ex("B_" + std::to_string(i), std::to_string(B[i])); + ex("level", std::to_string(layout.levelNumber())); + throw ex; + } +}; + + + + +template +template +void IonUpdater2::Pusher::move_domain( // + auto& pop, auto const& em, auto const& boxing, auto& interpolator) const +{ + using Particle_t = std::decay_t; + using ParticleLoop_t = std::conditional_t; + + + auto const& layout = boxing.layout; + auto& rhoP = pop.particleDensity(); + auto& rhoC = pop.chargeDensity(); + auto& F = pop.flux(); + auto& domain = pop.domainParticles(); + auto& patch_ghost = pop.patchGhostParticles(); + patch_ghost.reserve(static_cast(domain.size() / 100)); // do better + + for (std::size_t i = 0; i < domain.size(); ++i) // size might change on iteration! + { + ParticleLoop_t particle = domain[i]; + auto const oldCell = particle.iCell; + + move_particle(particle, em, layout, interpolator); + + if (oldCell == particle.iCell) + { + interpolator.particleToMesh(particle, rhoP, rhoC, F, layout); + continue; + } + + bool const isInDomainBox = boxing.isInDomainBox(particle); + bool const isInPatchGhostBox = boxing.isInPatchGhostBox(particle); + bool const should_interpolate = isInDomainBox or isInPatchGhostBox; + + if (!should_interpolate) + { + if constexpr (!copy_particle) + { + domain.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + } + continue; + } + + interpolator.particleToMesh(particle, rhoP, rhoC, F, layout); + + if constexpr (!copy_particle) + { + if (isInDomainBox) + { + domain.change_icell(particle, oldCell, i); + } + else + { + patch_ghost.push_back(particle); + domain.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + } + } + } +} + + + +template +template +void IonUpdater2::Pusher::move_level_ghost( // + auto& pop, auto const& em, auto const& boxing, auto& interpolator) const +{ + using Particle_t = std::decay_t; + using ParticleLoop_t = std::conditional_t; + + auto const& layout = boxing.layout; + + auto& domain = pop.domainParticles(); + auto& level_ghost = pop.levelGhostParticles(); + + for (std::size_t i = 0; i < level_ghost.size(); ++i) // size might change on iteration! + { + ParticleLoop_t particle = level_ghost[i]; + auto const oldCell = particle.iCell; + + move_particle(particle, em, layout, interpolator); + + if (oldCell == particle.iCell) + continue; + + bool const isInDomainBox = boxing.isInDomainBox(particle); + bool const should_interpolate = isInDomainBox; + + if (should_interpolate) + interpolator.particleToMesh( // + particle, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + + if constexpr (!copy_particle) + { + if (isInDomainBox) + { + domain.push_back(particle); + level_ghost.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + continue; + } + + bool const isInNonLevelGhostBox + = isInDomainBox || boxing.isInNonLevelGhostBox(particle.iCell); + bool const isInLevelGhostBox = !isInNonLevelGhostBox; + + if (isInLevelGhostBox) + level_ghost.change_icell(particle, oldCell, i); + else + { + level_ghost.swap_last_reduce_by_one(oldCell, i); + --i; // redo current index as last is now i + } + } + } +} + + +} // namespace PHARE::core + + +#endif // ION_UPDATER_HPP diff --git a/src/core/numerics/ion_updater/ion_updater_def.hpp b/src/core/numerics/ion_updater/ion_updater_def.hpp new file mode 100644 index 000000000..84298daf5 --- /dev/null +++ b/src/core/numerics/ion_updater/ion_updater_def.hpp @@ -0,0 +1,78 @@ +#ifndef PHARE_CORE_NUMERICS_ION_UPDATER_DEF_HPP +#define PHARE_CORE_NUMERICS_ION_UPDATER_DEF_HPP + + +#include "core/utilities/box/box.hpp" +#include "core/utilities/point/point.hpp" +#include "core/utilities/range/range.hpp" + +#include + +namespace PHARE::core +{ + + +enum class UpdaterMode { domain_only = 1, all = 2 }; + +// this is to detach how we partition particles from the updater directly +template +struct UpdaterSelectionBoxing +{ + auto constexpr static partGhostWidth = GridLayout::nbrParticleGhosts(); + using GridLayout_t = GridLayout; + using Box_t = GridLayout_t::AMRBox_t; + using ParticleRange = IndexRange; + using Selector_t = std::function; + + + GridLayout_t const layout; + std::vector const patchGhostBox; + Box_t const domainBox = layout.AMRBox(); + Box_t const ghostBox = grow(domainBox, partGhostWidth); + + bool isInGhostBox(auto& particle) const { return isIn(particle, ghostBox); }; + bool isInDomainBox(auto& particle) const { return isIn(particle, domainBox); }; + bool isInPatchGhostBox(auto& particle) const { return isIn(particle, patchGhostBox); }; + bool isInLevelGhostBox(auto& particle) const { return !isInNonLevelGhostBox(particle); }; + bool isInNonLevelGhostBox(auto& particle) const + { + return isInDomainBox(particle) or isIn(particle, patchGhostBox); + }; + + Selector_t const noop = [](auto& particleRange) { return particleRange; }; + + // lambda copy captures to detach from above references in case of class copy construct + Selector_t const inDomainBox = [domainBox = domainBox](auto& particleRange) { + return particleRange.array().partition( + particleRange, [&](auto const& particle) { return core::isIn(particle, domainBox); }); + }; + + Selector_t const inGhostBox = [ghostBox = ghostBox](auto& particleRange) { + return particleRange.array().partition( + particleRange, [&](auto const& particle) { return isIn(particle, ghostBox); }); + }; + + Selector_t const inNonLevelGhostBox + = [domainBox = domainBox, patchGhostBox = patchGhostBox](auto& particleRange) { + return particleRange.array().partition(particleRange, [&](auto const& particle) { + return isIn(particle, domainBox) or isIn(particle, patchGhostBox); + }); + }; + + Selector_t const inGhostLayer + = [ghostBox = ghostBox, domainBox = domainBox](auto& particleRange) { + return particleRange.array().partition(particleRange, [&](auto const& particle) { + return isIn(particle, ghostBox) and !isIn(particle, domainBox); + }); + }; + + Selector_t const outsideGhostBox = [ghostBox = ghostBox](auto& particleRange) { + return particleRange.array().partition( + particleRange, [&](auto const& particle) { return !isIn(particle, ghostBox); }); + }; +}; + +} // namespace PHARE::core + + +#endif // PHARE_CORE_NUMERICS_ION_UPDATER_DEF_HPP diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index d597eda57..3215a9d50 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -12,6 +12,7 @@ #include #include #include +#include namespace PHARE::core { @@ -269,8 +270,8 @@ NO_DISCARD bool isIn(Point const& point, Box const& box) } -template -NO_DISCARD auto isIn(Particle const& particle, Box const& box) +template +NO_DISCARD auto isIn(Particle const& particle, Box const& box) -> decltype(isIn(particle.iCell, box), bool()) { return isIn(particle.iCell, box); @@ -280,17 +281,10 @@ NO_DISCARD auto isIn(Particle const& particle, Box co * and returns true if the Point is at least in one of the boxes. * Returns occurs at the first box the point is in. */ -template = dummy::value> -bool isIn(Point const& point, BoxContainer const& boxes) +template +bool isIn(auto const& point, Boxes const& boxes) + requires(is_iterable_v) { - if (boxes.size() == 0) - return false; - - - static_assert(std::is_same::value, - "Box and Point should have the same data type"); - for (auto const& box : boxes) if (isIn(point, box)) return true; @@ -338,6 +332,71 @@ auto& operator<<(std::ostream& os, Box const& box) } + +template +std::vector> remove(Box const box, Box const& to_remove) +{ + using box_t = Box; + using _m = std::unordered_map; + + auto overlap = box * to_remove; + + if (not overlap) + return std::vector{box}; + + + auto copy = [](auto cpy, auto const& replace) { + for (auto const& [i, v] : replace) + cpy[i] = v; + return cpy; + }; + + auto intersection = *overlap; + + + std::unordered_map boxes; + + if (intersection.lower[0] > box.lower[0]) + boxes["left"] = Box(box.lower, copy(box.upper, _m{{0, intersection.lower[0] - 1}})); + if (intersection.upper[0] < box.upper[0]) + boxes["right"] = box_t{copy(box.lower, _m{{0, intersection.upper[0] + 1}}), box.upper}; + + [[maybe_unused]] Type minx = 0, maxx = 0; + if constexpr (dim > 1) + { + minx = boxes.count("left") > 0 ? intersection.lower[0] : box.lower[0]; + maxx = boxes.count("right") > 0 ? intersection.upper[0] : box.upper[0]; + + if (intersection.lower[1] > box.lower[1]) + boxes["down"] = box_t{copy(box.lower, _m{{0, minx}}), + copy(box.upper, _m{{0, maxx}, {1, intersection.lower[1] - 1}})}; + + if (intersection.upper[1] < box.upper[1]) + boxes["up"] = Box(copy(box.lower, _m{{0, minx}, {1, intersection.upper[1] + 1}}), + copy(box.upper, _m{{0, maxx}})); + } + + if constexpr (dim > 2) + { + Type miny = boxes.count("down") > 0 ? intersection.lower[1] : box.lower[1]; + Type maxy = boxes.count("up") > 0 ? intersection.upper[1] : box.upper[1]; + + if (intersection.lower[2] > box.lower[2]) + boxes["back"] = Box(copy(box.lower, _m{{0, minx}, {1, miny}}), + copy(intersection.lower - 1, _m{{0, maxx}, {1, maxy}})); + if (intersection.upper[2] < box.upper[2]) + boxes["front"] = Box(copy(intersection.upper + 1, _m{{0, minx}, {1, miny}}), + copy(box.upper, _m{{0, maxx}, {1, maxy}})); + } + + std::vector remaining; + for (auto const& [key, val] : boxes) + remaining.emplace_back(val); + return remaining; +} + + + } // namespace PHARE::core #endif diff --git a/src/core/utilities/equality.hpp b/src/core/utilities/equality.hpp new file mode 100644 index 000000000..56fa5dbac --- /dev/null +++ b/src/core/utilities/equality.hpp @@ -0,0 +1,26 @@ +#ifndef PHARE_CORE_UTILITIES_EQUALITY_HPP +#define PHARE_CORE_UTILITIES_EQUALITY_HPP + + +#include +#include + +namespace PHARE::core +{ + +struct EqualityReport +{ + auto operator()() const { return eq; } + operator bool() const { return eq; } + operator std::string() const { return reason; } + auto& what() const { return reason; } + auto& why() const { return reason; } + bool operator==(bool b) const { return eq == b; } + bool const eq = true; + std::string const reason = "=="; + std::size_t idx = 0; +}; + +} // namespace PHARE::core + +#endif /* PHARE_CORE_UTILITIES_EQUALITY_HPP */ diff --git a/src/core/utilities/types.hpp b/src/core/utilities/types.hpp index b320a00c2..4d032edf2 100644 --- a/src/core/utilities/types.hpp +++ b/src/core/utilities/types.hpp @@ -447,6 +447,15 @@ auto inline float_equals(double const& a, double const& b, double diff = 1e-12) return std::abs(a - b) < diff; } + +template +bool inline float_equals(std::array const& a, std::array const& b, + double diff = 1e-15) +{ + return for_N_all([&](auto i) { return float_equals(a[i], b[i], diff); }); +} + + template struct Apply { @@ -581,6 +590,13 @@ struct SetMax D& d; }; + +template +auto static constexpr as_tuple(std::array const& arr) +{ + return for_N([&](auto i) -> auto& { return arr[i]; }); +}; + } // namespace PHARE::core diff --git a/tests/core/data/field/test_field_fixtures.hpp b/tests/core/data/field/test_field_fixtures.hpp index c63c1ce97..726721478 100644 --- a/tests/core/data/field/test_field_fixtures.hpp +++ b/tests/core/data/field/test_field_fixtures.hpp @@ -1,7 +1,10 @@ #ifndef PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_HPP #define PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_HPP +#include "core/data/grid/grid.hpp" #include "core/data/field/field.hpp" +#include "core/utilities/equality.hpp" +#include "core/hybrid/hybrid_quantities.hpp" namespace PHARE::core { @@ -9,6 +12,106 @@ namespace PHARE::core template using Field_t = Field; + +template +struct FieldComparator +{ + auto float_eq(auto const a, auto const b) const + { + if constexpr (binary_eq) + return a == b; + else + return any_float_eq(a, b, diff); + }; + + template + auto operator()(F0 const& ref, F1 const& cmp) + { + auto const& ref_dat = ref.data(); + auto const& cmp_dat = cmp.data(); + for (std::size_t i = 0; i < ref.size(); ++i) + { + ref0 = ref_dat[i] == 0 ? ref0 + 1 : ref0; + cmp0 = cmp_dat[i] == 0 ? cmp0 + 1 : cmp0; + + nan0 = std::isnan(ref_dat[i]) ? nan0 + 1 : nan0; + nan1 = std::isnan(cmp_dat[i]) ? nan1 + 1 : nan1; + + if (!(std::isnan(ref_dat[i]) || std::isnan(cmp_dat[i]))) + { + auto ret = std::abs(ref_dat[i] - cmp_dat[i]); + if (ret < diff) + { + ++eqvals; + if (ref_dat[i] != 0 and cmp_dat[i] != 0) + ++eqnot0; + } + else + max_diff = ret > max_diff ? ret : max_diff; + } + } + ok = eqvals == ref.size() and nan0 == 0 and nan1 == 0; + return std::make_tuple(eqvals, eqnot0, ref0, cmp0); + } + + operator bool() const { return ok; } + + double const diff = 1e-15; + std::size_t eqvals = 0, eqnot0 = 0, ref0 = 0, cmp0 = 0, nan0 = 0, nan1 = 0; + bool ok = true; + double max_diff = 0; +}; + +using FloatFieldComparator_t = FieldComparator; + + +template +EqualityReport compare_fields(Field const& ref, Field const& cmp, + double const diff = 1e-15) +{ + auto const same_sizes = ref.size() == cmp.size(); + + if (!same_sizes) + return EqualityReport{false, "Tensorfield shape/size mismatch"}; + + std::stringstream log; + + FloatFieldComparator_t eq{diff}; + auto const [eqvals, eqnot0, ref0, cmp0] = eq(ref, cmp); + + std::string const names + = ref.name() == cmp.name() ? ref.name() : ref.name() + std::string{"/"} + cmp.name(); + log << "Fields compare (" << names << ") "; + + if (!eq) + { + auto const bad = ref.size() - eqvals; + log << "value mismatch: \n"; + log << "ok(" << eqvals << ") - "; + log << "ok!=0(" << eqnot0 << ") - "; + log << "bad(" << bad << ") - "; + log << "ref0(" << ref0 << ") - "; + log << "cmp0(" << cmp0 << ") - "; + log << "diff(" << eq.max_diff << ") - "; + log << "nan0(" << eq.nan0 << ") - "; + log << "nan1(" << eq.nan1 << ")\n"; + return EqualityReport{false, log.str()}; + } + + log << "are == with "; + log << "ok(" << eqvals << ") - "; + log << "ok!=0(" << eqnot0 << ") "; + + return EqualityReport{true, log.str()}; +} + +template +EqualityReport compare_fields(Grid const& ref, Grid const& cmp, + double const diff = 1e-15) +{ + return compare_fields(*ref, *cmp, diff); +} + } // namespace PHARE::core diff --git a/tests/core/data/gridlayout/test_gridlayout.hpp b/tests/core/data/gridlayout/test_gridlayout.hpp index 4137b9249..90cbb5e60 100644 --- a/tests/core/data/gridlayout/test_gridlayout.hpp +++ b/tests/core/data/gridlayout/test_gridlayout.hpp @@ -1,27 +1,52 @@ + #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" +#include "core/utilities/box/box.hpp" +#include "core/utilities/point/point.hpp" template -class TestGridLayout : public GridLayout -{ // to expose a default constructor +class TestGridLayout : public GridLayout // to expose a default constructor +{ + auto constexpr static dimension = GridLayout::dimension; + double constexpr static level_0_dx = .1; + public: + using Super = GridLayout; + auto static constexpr dim = GridLayout::dimension; TestGridLayout() = default; - TestGridLayout(std::uint32_t const cells) - : GridLayout{PHARE::core::ConstArray(1.0 / cells), + TestGridLayout(double const dl, std::uint32_t const cells) + : GridLayout{PHARE::core::ConstArray(dl), PHARE::core::ConstArray(cells), PHARE::core::Point{PHARE::core::ConstArray(0)}} { } - auto static make(std::uint32_t cells) { return TestGridLayout{cells}; } + TestGridLayout(std::uint32_t const cells, int const level = 0) + : GridLayout{PHARE::core::ConstArray(level_0_dx / (level + 1)), + PHARE::core::ConstArray(cells), + PHARE::core::Point{PHARE::core::ConstArray(0)}, + {}, + level} + { + } + + TestGridLayout(PHARE::core::Box const& box, auto const& origin, + int const level = 0) + : GridLayout{PHARE::core::ConstArray(level_0_dx / (level + 1)), box.shape(), + origin, box, level} + { + } + + auto static make(std::uint32_t const cells) { return TestGridLayout{cells}; } + + Super& operator*() { return *this; } + Super const& operator*() const { return *this; } }; #endif /*TESTS_CORE_DATA_GRIDLAYOUT_TEST_GRIDLAYOUT_HPP*/ diff --git a/tests/core/data/ion_population/test_ion_population_fixtures.hpp b/tests/core/data/ion_population/test_ion_population_fixtures.hpp index 54b0083f3..f1900cc97 100644 --- a/tests/core/data/ion_population/test_ion_population_fixtures.hpp +++ b/tests/core/data/ion_population/test_ion_population_fixtures.hpp @@ -1,8 +1,6 @@ #ifndef PHARE_TEST_CORE_DATA_ION_POPULATIONS_ION_POPULATION_FIXTURES_HPP #define PHARE_TEST_CORE_DATA_ION_POPULATIONS_ION_POPULATION_FIXTURES_HPP -#include -#include #include "phare_core.hpp" #include "core/data/ions/ions.hpp" @@ -14,19 +12,23 @@ #include "tests/core/data/vecfield/test_vecfield_fixtures.hpp" #include "tests/core/data/particles/test_particles_fixtures.hpp" + +#include + namespace PHARE::core { -template + +template struct UsableIonsDefaultTypes { public: - auto static constexpr dim = ParticleArray_::dimension; - auto static constexpr interp = interp_; + auto static constexpr dim = ParticleArray_::dimension; SimOpts static constexpr opts{dim, interp_}; using PHARE_Types = PHARE::core::PHARE_Types; using ParticleArray_t = ParticleArray_; + using Quantity = Quantity_; using Array_t = NdArrayVector; using Grid_t = Grid; using UsableVecField_t = UsableVecField; @@ -39,43 +41,87 @@ struct UsableIonsDefaultTypes }; +auto inline pop_dict(std::string const& name, std::size_t const ppc = 0) +{ + initializer::PHAREDict popdict; + popdict["name"] = name; + popdict["mass"] = 1.0; + auto particle_initializer = initializer::PHAREDict{}; + particle_initializer["nbr_part_per_cell"] = static_cast(ppc); + particle_initializer["charge"] = 1.; + particle_initializer["basis"] = std::string{"cartesian"}; + popdict["particle_initializer"] = particle_initializer; + return popdict; +} + template class UsableIonsPopulation_ : public _defaults::IonPopulation_t { + using Quantity = _defaults::Quantity; using GridLayout_t = _defaults::GridLayout_t; using VecField_t = _defaults::VecField_t; using TensorField_t = _defaults::TensorField_t; using ParticleArray_t = _defaults::ParticleArray_t; using Super = IonPopulation; -public: - UsableIonsPopulation_(initializer::PHAREDict const& dict, GridLayout_t const& layout) - : Super{dict} - , particleDensity{this->name() + "_particleDensity", layout, HybridQuantity::Scalar::rho} - , chargeDensity{this->name() + "_chargeDensity", layout, HybridQuantity::Scalar::rho} - , F{this->name() + "_flux", layout, HybridQuantity::Vector::V} - , M{this->name() + "_momentumTensor", layout, HybridQuantity::Tensor::M} - , particles{this->name(), layout.AMRBox()} + + + void set() { auto&& [_F, _M, _pd, _cd, _particles] = Super::getCompileTimeResourcesViewList(); F.set_on(_F); M.set_on(_M); - _pd.setBuffer(&particleDensity); - _cd.setBuffer(&chargeDensity); + _pd.setBuffer(&rhoP); + _cd.setBuffer(&rhoC); _particles.setBuffer(&particles.pack()); + + // PHARE_LOG_LINE_SS(_pd.data()[0] << " == " << _pd.data()[_pd.size() - 1]); + assert(compare_fields(_pd, *rhoP)); } - Super& view() { return *this; } - Super const& view() const { return *this; } - auto& operator*() { return view(); } - auto& operator*() const { return view(); } +public: + UsableIonsPopulation_(initializer::PHAREDict const& dict, GridLayout_t const& layout) + : Super{dict} + , layout_{layout} + , rhoP{this->name() + "_particleDensity", layout_, Quantity::Scalar::rho, 1} + , rhoC{this->name() + "_chargeDensity", layout_, Quantity::Scalar::rho, 1} + , F{this->name() + "_flux", layout, Quantity::Vector::V, 1} + , M{this->name() + "_momentumTensor", layout, Quantity::Tensor::M, 1} + , particles{this->name(), grow(layout.AMRBox(), GridLayout_t::nbrParticleGhosts() + 1)} + { + set(); + } + + UsableIonsPopulation_(UsableIonsPopulation_ const& that) + : Super{pop_dict(that.name())} + , layout_{that.layout_} + , rhoP{that.rhoP} + , rhoC{that.rhoC} + , F{that.F} + , M{that.M} + , particles{that.particles} + { + set(); + } + + UsableIonsPopulation_(UsableIonsPopulation_&&) = default; + UsableIonsPopulation_& operator=(UsableIonsPopulation_ const&) = delete; + UsableIonsPopulation_& operator=(UsableIonsPopulation_&&) = default; - _defaults::Grid_t particleDensity, chargeDensity; + + Super& super() { return *this; } + Super const& super() const { return *this; } + auto& operator*() { return super(); } + auto& operator*() const { return super(); } + + GridLayout_t layout_; + _defaults::Grid_t rhoP, rhoC; _defaults::UsableVecField_t F; _defaults::UsableTensorField_t M; UsableParticlesPopulation particles; }; + template using UsableIonsPopulation = UsableIonsPopulation_>; @@ -84,68 +130,122 @@ template class UsableIons_ : public Ions { + using Quantity = _defaults::Quantity; using GridLayout_t = _defaults::GridLayout_t; using Super = Ions; - auto static pop_dict(std::string name) - { - initializer::PHAREDict popdict; - popdict["name"] = name; - popdict["mass"] = 1.0; - popdict["particle_initializer"] = initializer::PHAREDict{}; - return popdict; - } + template - auto static super(PopNames const& pop_names) + auto static super(PopNames const& pop_names, std::size_t const ppc) { initializer::PHAREDict dict; dict["nbrPopulations"] = pop_names.size(); for (std::size_t i = 0; i < pop_names.size(); ++i) - dict["pop" + std::to_string(i)] = pop_dict(pop_names[i]); + dict["pop" + std::to_string(i)] = pop_dict(pop_names[i], ppc); return dict; } -public: - UsableIons_(GridLayout_t const& layout, std::vector const& pop_names) - : Super{super(pop_names)} - , massDensity{"massDensity", layout, HybridQuantity::Scalar::rho} - , chargeDensity{"chargeDensity", layout, HybridQuantity::Scalar::rho} - , Vi{"bulkVel", layout, HybridQuantity::Vector::V} - , M{"momentumTensor", layout, HybridQuantity::Tensor::M} + + + auto static super(Super const supe) + { + initializer::PHAREDict dict; + dict["nbrPopulations"] = supe.size(); + for (std::size_t i = 0; i < supe.size(); ++i) + dict["pop" + std::to_string(i)] = pop_dict(supe[i].name()); + return dict; + } + + void set() { auto&& [_bV, _M, _cd, _md] = Super::getCompileTimeResourcesViewList(); Vi.set_on(_bV); M.set_on(_M); - _cd.setBuffer(&chargeDensity); - _md.setBuffer(&massDensity); + _cd.setBuffer(&rhoC); + _md.setBuffer(&rhoM); - for (std::size_t i = 0; i < pop_names.size(); ++i) - populations.emplace_back(pop_dict(pop_names[i]), layout); - Super::getRunTimeResourcesViewList().clear(); + auto& super_pops = Super::getRunTimeResourcesViewList(); + super_pops.clear(); for (auto& pop : populations) - Super::getRunTimeResourcesViewList().emplace_back(*pop); + super_pops.emplace_back(*pop); + } + + +public: + using ParticleArray_t = typename _defaults::ParticleArray_t; + + UsableIons_(GridLayout_t const& layout, initializer::PHAREDict const& dict) + : Super{dict} + , layout_{layout} + , rhoM{"massDensity", layout_, Quantity::Scalar::rho, 1} + , rhoC{"chargeDensity", layout_, Quantity::Scalar::rho, 1} + , Vi{"bulkVel", layout, Quantity::Vector::V, 1} + , M{"momentumTensor", layout, Quantity::Tensor::M, 0} + { + auto& super_pops = Super::getRunTimeResourcesViewList(); + populations.reserve(super_pops.size()); + for (std::size_t i = 0; i < super_pops.size(); ++i) + populations.emplace_back(dict["pop" + std::to_string(i)], layout); + set(); + } + + UsableIons_(GridLayout_t const& layout, std::vector const& pop_names, + std::size_t const ppc = 0) + : UsableIons_{layout, super(pop_names, ppc)} + { + } + + UsableIons_(GridLayout_t const& layout, std::string const& pop = "protons") + : UsableIons_{layout, std::vector{pop}} + { } - UsableIons_(GridLayout_t const& layout, std::string const& pop_name = "protons") - : UsableIons_{layout, std::vector{pop_name}} + UsableIons_(UsableIons_&& that) + : Super(super(*that)) + , layout_{that.layout_} + , rhoM{std::move(that.rhoM)} + , rhoC{std::move(that.rhoC)} + , Vi{std::move(that.Vi)} + , M{std::move(that.M)} + , populations{std::move(that.populations)} { + set(); } - Super& view() { return *this; } - Super const& view() const { return *this; } - auto& operator*() { return view(); } - auto& operator*() const { return view(); } + UsableIons_(UsableIons_ const& that) + : Super(super(*that)) + , layout_{that.layout_} + , rhoM{that.rhoM} + , rhoC{that.rhoC} + , Vi{that.Vi} + , M{that.M} + , populations{that.populations} + { + set(); + } + + UsableIons_& operator=(UsableIons_ const&) = delete; + UsableIons_& operator=(UsableIons_&&) = default; + + Super& super() { return *this; } + Super const& super() const { return *this; } + auto& operator*() { return super(); } + auto& operator*() const { return super(); } - _defaults::Grid_t massDensity, chargeDensity; + GridLayout_t layout_; + _defaults::Grid_t rhoM, rhoC; _defaults::UsableVecField_t Vi; _defaults::UsableTensorField_t M; std::vector> populations; }; + + template using UsableIons = UsableIons_>; + } // namespace PHARE::core diff --git a/tests/core/data/particles/test_particles_fixtures.hpp b/tests/core/data/particles/test_particles_fixtures.hpp index cba9887e7..456bde42a 100644 --- a/tests/core/data/particles/test_particles_fixtures.hpp +++ b/tests/core/data/particles/test_particles_fixtures.hpp @@ -1,6 +1,12 @@ #ifndef PHARE_TEST_CORE_DATA_PARTICLES_TEST_PARTICLES_FIXTURES_HPP #define PHARE_TEST_CORE_DATA_PARTICLES_TEST_PARTICLES_FIXTURES_HPP +#include "core/utilities/types.hpp" +#include "core/utilities/equality.hpp" +#include "core/data/particles/particle.hpp" +#include "core/data/ions/ion_population/particle_pack.hpp" + +#include #include namespace PHARE::core @@ -23,14 +29,146 @@ struct UsableParticlesPopulation ParticleArray_t domain_particles; ParticleArray_t patch_ghost_particles = domain_particles; ParticleArray_t level_ghost_particles = domain_particles; - core::ParticlesPack particles_pack{name, - &domain_particles, // - &patch_ghost_particles, - &level_ghost_particles, - /*levelGhostParticlesOld=*/nullptr, - /*levelGhostParticlesNew=*/nullptr}; + ParticlesPack particles_pack{name, + &domain_particles, // + &patch_ghost_particles, + &level_ghost_particles, + /*levelGhostParticlesOld=*/nullptr, + /*levelGhostParticlesNew=*/nullptr}; +}; + +template> +Particle_t particle(std::array const& icell) +{ + return { + /*.weight = */ .001, + /*.charge = */ 1, + /*.iCell = */ icell, + /*.delta = */ ConstArray(.5), + /*.v = */ {{.002002002002, .003003003003, .004004004004}} // + }; +} + +template +void add_particles_in(Particles& particles, Box const& box, std::size_t const ppc) +{ + for (auto const& bix : box) + for (std::size_t i = 0; i < ppc; ++i) + particles.emplace_back(particle(*bix)); +} + + +auto inline rando(std::optional seed = std::nullopt) +{ + if (seed.has_value()) + return std::mt19937_64(*seed); + std::random_device rd; + std::seed_seq seed_seq{rd(), rd(), rd(), rd(), rd()}; + return std::mt19937_64(seed_seq); +} + +template +void delta_disperse(Particles& particles, std::optional seed = std::nullopt) +{ + auto gen = rando(seed); + ParticleDeltaDistribution deltaDistrib; + for (auto& p : particles) + p.delta = for_N_make_array([&](auto) { return deltaDistrib(gen); }); +} + +template +void vary_velocity(Particles& particles, double const min, double const max, + std::optional seed = std::nullopt) +{ + auto gen = rando(seed); + std::uniform_real_distribution dist{min, max}; + for (auto& p : particles) + p.v = for_N_make_array<3>([&](auto) { return dist(gen); }); +} + +template +struct CellFlattener +{ + template + NO_DISCARD RValue operator()(Icell const& icell) const + { + if constexpr (Box_t::dimension == 2) + return icell[1] + icell[0] * shape[1] * shape[0]; + if constexpr (Box_t::dimension == 3) + return icell[2] + icell[1] * shape[2] + icell[0] * shape[1] * shape[2]; + return icell[0]; + } + + Box_t const box; + std::array shape = box.shape().toArray(); }; +template // support const vs non-const iterators +auto it_dist(I0&& i0, I1&& i1) +{ + return std::distance(i0, i1); +} + + +template +auto& sort_particles(GridLayout const& layout, auto& particles) +{ + auto constexpr cmp_deltas = [](auto const& a, auto const& b) -> bool { + return as_tuple(a.delta) < as_tuple(b.delta); + }; + auto const by_deltas = [&](std::uint64_t const& l, std::uint64_t const& r) { + std::sort(particles.begin() + l, particles.begin() + r, cmp_deltas); + }; + + auto const ghostBox = grow(layout.AMRBox(), GridLayout::nbrParticleGhosts()); + CellFlattener cell_flattener{ghostBox}; + + std::sort( + particles.begin(), particles.end(), + [cf = cell_flattener](auto const& a, auto const& b) { return cf(a.iCell) < cf(b.iCell); }); + + auto const end = particles.end(); + auto beg = particles.begin(); + auto lst = particles.begin(); + + auto const check = [&]() { return lst != end and lst->iCell == beg->iCell; }; + + while (lst != end) + { + lst = beg + 1; + while (check()) + ++lst; + auto const s = it_dist(particles.begin(), beg); + auto const e = it_dist(particles.begin(), lst); + by_deltas(s, e); + beg = lst; + } + + return particles; +} + +EqualityReport compare_particles(auto const& ps0, auto const& ps1, double const atol) +{ + if (ps0.size() != ps1.size()) + return EqualityReport{false, "different sizes: " + std::to_string(ps0.size()) + " vs " + + std::to_string(ps1.size())}; + + for (std::size_t i = 0; i < ps0.size(); ++i) + { + std::string const idx = std::to_string(i); + if (ps0[i].iCell != ps1[i].iCell) + return EqualityReport{false, "icell mismatch at index: " + idx, i}; + + if (!float_equals(ps0[i].v, ps1[i].v, atol)) + return EqualityReport{false, "v mismatch at index: " + idx, i}; + + if (!float_equals(ps0[i].delta, ps1[i].delta, atol)) + return EqualityReport{false, "delta mismatch at index: " + idx, i}; + } + + return EqualityReport{true}; +} + } // namespace PHARE::core diff --git a/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp b/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp index baf09ef6c..928b97d6a 100644 --- a/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp +++ b/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp @@ -3,6 +3,7 @@ #include "core/data/grid/grid.hpp" #include "core/data/field/field.hpp" +#include "core/utilities/equality.hpp" #include "core/hybrid/hybrid_quantities.hpp" #include "core/data/tensorfield/tensorfield.hpp" @@ -23,6 +24,11 @@ template class UsableTensorField : public TensorField, HybridQuantity, rank_> { auto constexpr static N_elements = detail::tensor_field_dim_from_rank(); + void _set() + { + for (std::size_t i = 0; i < N_elements; ++i) + super()[i].setBuffer(&xyz[i]); + } public: auto static constexpr dimension = dim; @@ -31,12 +37,29 @@ class UsableTensorField : public TensorField, HybridQuantity, rank_ using tensor_t = typename Super::tensor_t; template - UsableTensorField(std::string const& name, GridLayout const& layout, tensor_t qty) + UsableTensorField(std::string const& name, GridLayout const& layout, tensor_t qty, + std::optional v = std::nullopt) : Super{name, qty} , xyz{make_grids(Super::componentNames(), layout, qty)} { - for (std::size_t i = 0; i < N_elements; ++i) - super()[i].setBuffer(&xyz[i]); + if (v) + for (std::size_t i = 0; i < N_elements; ++i) + xyz[i].fill(*v); + _set(); + } + + UsableTensorField(UsableTensorField&& that) + : Super{std::forward(that)} + , xyz{std::move(that.xyz)} + { + _set(); + } + + UsableTensorField(UsableTensorField const& that) + : Super{that} + , xyz{that.xyz} + { + _set(); } void set_on(Super& tensorfield) @@ -62,6 +85,38 @@ class UsableTensorField : public TensorField, HybridQuantity, rank_ }; + + +template +EqualityReport compare_tensor_fields(TensorField const& ref, + TensorField const& cmp, + double const diff) +{ + auto constexpr static N_elements = detail::tensor_field_dim_from_rank(); + + if (ref.componentNames() != cmp.componentNames()) + return EqualityReport{false, "Tensorfield component mismatch"}; + + auto const same_sizes = [&]() { + return core::for_N_all([&](auto i) { return ref[i].size() == cmp[i].size(); }); + }(); + + if (!same_sizes) + return EqualityReport{false, "Tensorfield shape/size mismatch"}; + + std::stringstream log; + log << std::endl; + for (std::size_t ci = 0; ci < N_elements; ++ci) + if (auto eqr = compare_fields(ref[ci], cmp[ci], diff); !eqr) + return eqr; + else + log << eqr.what() << std::endl; + + return EqualityReport{true, log.str()}; +} + + + } // namespace PHARE::core diff --git a/tests/core/numerics/ion_updater/CMakeLists.txt b/tests/core/numerics/ion_updater/CMakeLists.txt index e206a3e58..97030b1bc 100644 --- a/tests/core/numerics/ion_updater/CMakeLists.txt +++ b/tests/core/numerics/ion_updater/CMakeLists.txt @@ -2,19 +2,18 @@ cmake_minimum_required (VERSION 3.20.1) project(test-updater) -set(SOURCES test_updater.cpp) +function(_updater_test test_name) -add_executable(${PROJECT_NAME} ${SOURCES}) + add_executable(${test_name} ${test_name}.cpp) -target_include_directories(${PROJECT_NAME} PRIVATE - ${GTEST_INCLUDE_DIRS} - ) + target_include_directories(${test_name} PRIVATE ${GTEST_INCLUDE_DIRS}) -target_link_libraries(${PROJECT_NAME} PRIVATE - phare_core - phare_simulator - ${GTEST_LIBS}) + target_link_libraries(${test_name} PRIVATE phare_core phare_initializer ${GTEST_LIBS}) -add_no_mpi_phare_test(${PROJECT_NAME} ${CMAKE_CURRENT_BINARY_DIR}) + add_no_mpi_phare_test(${test_name} ${CMAKE_CURRENT_BINARY_DIR}) +endfunction(_updater_test) + +_updater_test(test_updater) +_updater_test(test_updaters) diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index 4066b332b..9bbe08dfa 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -1,4 +1,5 @@ #include "gtest/gtest.h" +#include #include "phare_core.hpp" @@ -21,17 +22,17 @@ Return density(Param x) Return vx(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>(x.size(), .1); } Return vy(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>(x.size(), .1); } Return vz(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>(x.size(), .1); } Return vthx(Param x) @@ -77,64 +78,37 @@ PHARE::initializer::PHAREDict createDict() dict["simulation"]["algo"]["ion_updater"]["pusher"]["name"] = std::string{"modified_boris"}; - dict["ions"]["nbrPopulations"] = std::size_t{2}; - dict["ions"]["pop0"]["name"] = std::string{"protons"}; - dict["ions"]["pop0"]["mass"] = 1.; - dict["ions"]["pop0"]["particle_initializer"]["name"] = std::string{"maxwellian"}; - dict["ions"]["pop0"]["particle_initializer"]["density"] = static_cast(density); - - dict["ions"]["pop0"]["particle_initializer"]["bulk_velocity_x"] - = static_cast(vx); - - dict["ions"]["pop0"]["particle_initializer"]["bulk_velocity_y"] - = static_cast(vy); - - dict["ions"]["pop0"]["particle_initializer"]["bulk_velocity_z"] - = static_cast(vz); - - - dict["ions"]["pop0"]["particle_initializer"]["thermal_velocity_x"] - = static_cast(vthx); - - dict["ions"]["pop0"]["particle_initializer"]["thermal_velocity_y"] - = static_cast(vthy); - - dict["ions"]["pop0"]["particle_initializer"]["thermal_velocity_z"] - = static_cast(vthz); - - - dict["ions"]["pop0"]["particle_initializer"]["nbr_part_per_cell"] = int{nbrPartPerCell}; - dict["ions"]["pop0"]["particle_initializer"]["charge"] = 1.; - dict["ions"]["pop0"]["particle_initializer"]["basis"] = std::string{"cartesian"}; - - dict["ions"]["pop1"]["name"] = std::string{"alpha"}; - dict["ions"]["pop1"]["mass"] = 1.; - dict["ions"]["pop1"]["particle_initializer"]["name"] = std::string{"maxwellian"}; - dict["ions"]["pop1"]["particle_initializer"]["density"] = static_cast(density); - - dict["ions"]["pop1"]["particle_initializer"]["bulk_velocity_x"] - = static_cast(vx); - - dict["ions"]["pop1"]["particle_initializer"]["bulk_velocity_y"] - = static_cast(vy); - - dict["ions"]["pop1"]["particle_initializer"]["bulk_velocity_z"] - = static_cast(vz); - - - dict["ions"]["pop1"]["particle_initializer"]["thermal_velocity_x"] - = static_cast(vthx); - - dict["ions"]["pop1"]["particle_initializer"]["thermal_velocity_y"] - = static_cast(vthy); - - dict["ions"]["pop1"]["particle_initializer"]["thermal_velocity_z"] - = static_cast(vthz); - - - dict["ions"]["pop1"]["particle_initializer"]["nbr_part_per_cell"] = int{nbrPartPerCell}; - dict["ions"]["pop1"]["particle_initializer"]["charge"] = 1.; - dict["ions"]["pop1"]["particle_initializer"]["basis"] = std::string{"cartesian"}; + dict["ions"]["nbrPopulations"] = std::size_t{2}; + + dict["ions"]["pop0"]["name"] = std::string{"protons"}; + dict["ions"]["pop0"]["mass"] = 1.; + auto& pop0_initializer = dict["ions"]["pop0"]["particle_initializer"]; + pop0_initializer["name"] = std::string{"maxwellian"}; + pop0_initializer["density"] = static_cast(density); + pop0_initializer["bulk_velocity_x"] = static_cast(vx); + pop0_initializer["bulk_velocity_y"] = static_cast(vy); + pop0_initializer["bulk_velocity_z"] = static_cast(vz); + pop0_initializer["thermal_velocity_x"] = static_cast(vthx); + pop0_initializer["thermal_velocity_y"] = static_cast(vthy); + pop0_initializer["thermal_velocity_z"] = static_cast(vthz); + pop0_initializer["nbr_part_per_cell"] = int{nbrPartPerCell}; + pop0_initializer["charge"] = 1.; + pop0_initializer["basis"] = std::string{"cartesian"}; + + dict["ions"]["pop1"]["name"] = std::string{"alpha"}; + dict["ions"]["pop1"]["mass"] = 1.; + auto& pop1_initializer = dict["ions"]["pop1"]["particle_initializer"]; + pop1_initializer["name"] = std::string{"maxwellian"}; + pop1_initializer["density"] = static_cast(density); + pop1_initializer["bulk_velocity_x"] = static_cast(vx); + pop1_initializer["bulk_velocity_y"] = static_cast(vy); + pop1_initializer["bulk_velocity_z"] = static_cast(vz); + pop1_initializer["thermal_velocity_x"] = static_cast(vthx); + pop1_initializer["thermal_velocity_y"] = static_cast(vthy); + pop1_initializer["thermal_velocity_z"] = static_cast(vthz); + pop1_initializer["nbr_part_per_cell"] = int{nbrPartPerCell}; + pop1_initializer["charge"] = 1.; + pop1_initializer["basis"] = std::string{"cartesian"}; dict["electromag"]["name"] = std::string{"EM"}; dict["electromag"]["electric"]["name"] = std::string{"E"}; @@ -150,14 +124,6 @@ static auto init_dict = createDict(); -template -struct DimInterp -{ - static constexpr auto dimension = dim; - static constexpr auto interp_order = interporder; -}; - - // the Electromag and Ions used in this test // need their resources pointers (Fields and ParticleArrays) to set manually @@ -170,9 +136,9 @@ struct ElectromagBuffers { constexpr static PHARE::SimOpts opts{dim, interp_order}; using PHARETypes = PHARE::core::PHARE_Types; - using Grid = typename PHARETypes::Grid_t; - using GridLayout = typename PHARETypes::GridLayout_t; - using Electromag = typename PHARETypes::Electromag_t; + using Grid = PHARETypes::Grid_t; + using GridLayout = PHARETypes::GridLayout_t; + using Electromag = PHARETypes::Electromag_t; using UsableVecFieldND = UsableVecField; UsableVecFieldND B, E; @@ -356,23 +322,62 @@ struct IonsBuffers }; +template +struct TestParam +{ + auto static constexpr opts = opts_; + auto static constexpr impl = impl_; + static constexpr auto dimension = opts.dimension; + static constexpr auto interp_order = opts.interp_order; +}; + +template +struct updater_test_bits; -template -struct IonUpdaterTest : public ::testing::Test +template +struct updater_test_bits { - static constexpr auto dim = DimInterpT::dimension; - static constexpr auto interp_order = DimInterpT::interp_order; - constexpr static PHARE::SimOpts opts{dim, interp_order}; - using PHARETypes = PHARE::core::PHARE_Types; - using Ions = typename PHARETypes::Ions_t; - using Electromag = typename PHARETypes::Electromag_t; - using GridLayout = typename PHARE::core::GridLayout>; - using ParticleArray = typename PHARETypes::ParticleArray_t; - using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory_t; + using PHARETypes = PHARE::core::PHARE_Types; + using Ions = PHARETypes::Ions_t; + using IonUpdater = IonUpdaterProxy>; +}; + +template +struct updater_test_bits +{ + using PHARETypes = PHARE::core::PHARE_Types; + using Ions = PHARETypes::Ions_t; + using IonUpdater = IonUpdaterProxy>; +}; + +template +struct updater_test_bits +{ + using PHARETypes = PHARE::core::PHARE_Types; + using Ions = PHARETypes::Ions_t; + using IonUpdater = IonUpdaterProxy>; +}; + - using IonUpdater = typename PHARE::core::IonUpdater; - using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + +template +struct IonUpdaterTest : public ::testing::Test, + public updater_test_bits +{ + auto static constexpr opts = TestParam_t::opts; + static constexpr auto dim = opts.dimension; + static constexpr auto interp_order = opts.interp_order; + using PHARETypes = PHARE::core::PHARE_Types; + using Ions = PHARETypes::Ions_t; + using Electromag = PHARETypes::Electromag_t; + using ParticleArray = PHARETypes::ParticleArray_t; + using ParticleInitializerFactory = PHARETypes::ParticleInitializerFactory_t; + using GridLayout = PHARETypes::GridLayout_t; + + using basics = updater_test_bits; + using IonUpdater = basics::IonUpdater; + using Boxing_t = UpdaterSelectionBoxing; double dt{0.01}; @@ -380,8 +385,10 @@ struct IonUpdaterTest : public ::testing::Test // grid configuration std::array ncells; GridLayout layout; - // assumes no level ghost cells - Boxing_t const boxing{layout, {grow(layout.AMRBox(), GridLayout::nbrParticleGhosts())}}; + // assumes only level ghost cells + Box const domainBox = layout.AMRBox(); + Box const ghostBox = grow(domainBox, GridLayout::nbrParticleGhosts()); + Boxing_t const boxing{layout, remove(ghostBox, domainBox)}; // data for electromagnetic fields @@ -529,7 +536,7 @@ struct IonUpdaterTest : public ::testing::Test void fillIonsMomentsGhosts() { - using Interpolator = typename IonUpdater::Interpolator; + using Interpolator = IonUpdater::Interpolator_t; Interpolator interpolate; for (auto& pop : this->ions) @@ -656,11 +663,20 @@ struct IonUpdaterTest : public ::testing::Test -using DimInterps = ::testing::Types, DimInterp<1, 2>, DimInterp<1, 3>>; - - -TYPED_TEST_SUITE(IonUpdaterTest, DimInterps, ); - +// clang-format off +using Permutations = ::testing::Types< + TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam +>; +// clang-format on +TYPED_TEST_SUITE(IonUpdaterTest, Permutations, ); @@ -741,8 +757,8 @@ TYPED_TEST(IonUpdaterTest, particlesUntouchedInMomentOnlyMode) IonsBuffers ionsBufferCpy{this->ionsBuffers, this->layout}; - ionUpdater.updatePopulations(this->ions, this->EM, this->boxing, this->dt, - UpdaterMode::domain_only); + ionUpdater.updatePopulations(UpdaterMode::domain_only, this->ions, this->EM, this->boxing, + this->dt); this->fillIonsMomentsGhosts(); @@ -812,7 +828,7 @@ TYPED_TEST(IonUpdaterTest, momentsAreChangedInParticlesAndMomentsMode) IonsBuffers ionsBufferCpy{this->ionsBuffers, this->layout}; - ionUpdater.updatePopulations(this->ions, this->EM, this->boxing, this->dt, UpdaterMode::all); + ionUpdater.updatePopulations(UpdaterMode::all, this->ions, this->EM, this->boxing, this->dt); this->fillIonsMomentsGhosts(); @@ -832,8 +848,8 @@ TYPED_TEST(IonUpdaterTest, momentsAreChangedInMomentsOnlyMode) IonsBuffers ionsBufferCpy{this->ionsBuffers, this->layout}; - ionUpdater.updatePopulations(this->ions, this->EM, this->boxing, this->dt, - UpdaterMode::domain_only); + ionUpdater.updatePopulations(UpdaterMode::domain_only, this->ions, this->EM, this->boxing, + this->dt); this->fillIonsMomentsGhosts(); @@ -850,8 +866,8 @@ TYPED_TEST(IonUpdaterTest, thatNoNaNsExistOnPhysicalNodesMoments) typename IonUpdaterTest::IonUpdater ionUpdater{ init_dict["simulation"]["algo"]["ion_updater"]}; - ionUpdater.updatePopulations(this->ions, this->EM, this->boxing, this->dt, - UpdaterMode::domain_only); + ionUpdater.updatePopulations(UpdaterMode::domain_only, this->ions, this->EM, this->boxing, + this->dt); this->fillIonsMomentsGhosts(); diff --git a/tests/core/numerics/ion_updater/test_updaters.cpp b/tests/core/numerics/ion_updater/test_updaters.cpp new file mode 100644 index 000000000..4b50725a2 --- /dev/null +++ b/tests/core/numerics/ion_updater/test_updaters.cpp @@ -0,0 +1,287 @@ + + +#include "core/data/particles/particle_array.hpp" +#include "core/numerics/ion_updater/ion_updater.hpp" // IWYU pragma: keep + +#include "tests/core/data/electromag/test_electromag_fixtures.hpp" +#include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" + + +#include "gtest/gtest.h" +#include +#include + + +namespace PHARE::core +{ +// RUNTIME ENV VAR OVERRIDES +auto static const cells = get_env_as("PHARE_CELLS", std::uint32_t{10}); +auto static const ppc = get_env_as("PHARE_PPC", std::size_t{500}); +auto static const seeder = get_env_as("PHARE_SEED", std::size_t{0}); // 0 == NO SEED! +auto static const n_patches = get_env_as("PHARE_PATCHES", std::size_t{1}); +auto static const dt = get_env_as("PHARE_TIMESTEP", double{.001}); +auto static const do_cmp = get_env_as("PHARE_COMPARE", std::size_t{1}); +auto static const cmp_only = get_env_as("PHARE_CMP_ONLY", std::size_t{0}); +auto static const ref_only = get_env_as("PHARE_REF_ONLY", std::size_t{0}); + +bool static const premain = []() { return true; }(); + + +template +auto get_updater_for(Ions& /*ions*/, EM const& /*em*/, GridLayout_t const& /*layout*/) +{ + if constexpr (impl == 0) + return IonUpdaterProxy>(); + if constexpr (impl == 1) + return IonUpdaterProxy>(); + if constexpr (impl == 2) + return IonUpdaterProxy>(); +} + + +template +void ref_update(UpdaterMode mode, Patches& patches) +{ + using GridLayout_t = Patches::value_type::GridLayout_t; + using Particles = Patches::value_type::ParticleArray_t; + using Boxing_t = UpdaterSelectionBoxing; + + for (auto& [layout, ions, _, electromag] : patches) + { + auto const domainBox = layout.AMRBox(); + auto const ghostBox = grow(domainBox, GridLayout_t::nbrParticleGhosts()); + auto updater = get_updater_for<0>(*ions, electromag, layout); + + Boxing_t const boxing{layout, remove(ghostBox, domainBox)}; + updater.updatePopulations(mode, *ions, electromag, boxing, dt); + } +} + +template +void cmp_update(UpdaterMode mode, Patches& patches) +{ + using GridLayout_t = Patches::value_type::GridLayout_t; + using Particles = Patches::value_type::ParticleArray_t; + using Boxing_t = UpdaterSelectionBoxing; + + for (auto& [layout, ions, _, electromag] : patches) + { + auto const domainBox = layout.AMRBox(); + auto const ghostBox = grow(domainBox, GridLayout_t::nbrParticleGhosts()); + auto updater = get_updater_for(*ions, electromag, layout); + + Boxing_t const boxing{layout, remove(ghostBox, domainBox)}; + updater.updatePopulations(mode, *ions, electromag, boxing, dt); + } +} + + +template +auto make_ions(GridLayout_t const& layout) +{ + auto constexpr static interp = GridLayout_t::interp_order; + + UsableIons ions{layout}; + std::optional seed = std::nullopt; + if (seeder > 0) + seed = seeder; + + EXPECT_EQ(ions.populations[0].particles.domain_particles.size(), 0ull); + + auto const disperse = [&](auto& particles) { + delta_disperse(particles.domain_particles, seed); + vary_velocity(particles.domain_particles, -6, 6, seed); + + delta_disperse(particles.level_ghost_particles, seed); + vary_velocity(particles.level_ghost_particles, -6, 6, seed); + }; + + auto const particle_box = layout.AMRBox(); + + auto add_particles = [&](auto& particles) { + particles.domain_particles.reserve(particle_box.size() * ppc); + add_particles_in(particles.domain_particles, particle_box, ppc); + }; + + add_particles(ions.populations[0].particles); + disperse(ions.populations[0].particles); + + + EXPECT_EQ(ions.populations[0].particles.domain_particles.size(), particle_box.size() * ppc); + + return ions; +} + + + +template +auto from_ions(GridLayout_t const& layout, Ions const& from) +{ + auto constexpr static interp = GridLayout_t::interp_order; + + UsableIons ions{layout, "protons"}; + EXPECT_EQ(ions.populations[0].particles.domain_particles.size(), 0ull); + + auto _add_particles_from = [&](auto& src, auto& dst) { dst = src; }; + + _add_particles_from.template operator()(from.populations[0].particles.domain_particles, + ions.populations[0].particles.domain_particles); + + _add_particles_from.template operator()(from.populations[0].particles.level_ghost_particles, + ions.populations[0].particles.level_ghost_particles); + + auto const particle_box = layout.AMRBox(); + EXPECT_EQ(ions.populations[0].particles.domain_particles.size(), particle_box.size() * ppc); + return std::move(ions); +} + + +template +void check_particles(GridLayout_t const& layout, P0& ref, P1& cmp, double const atol) +{ + auto const box = layout.AMRBox(); + sort_particles(layout, ref); + sort_particles(layout, cmp); + + EXPECT_EQ(ref.size(), cmp.size()); + + auto const report = compare_particles(ref, cmp, atol); + PHARE_LOG_LINE_STR("results: " << report.why()); + PHARE_LOG_LINE_STR("eg: " << ref[0]); + EXPECT_TRUE(report); +} + +template +void compare(GridLayout_t const& layout, R& ref, C& cmp) +{ + double diff = 1e-14; + + check_particles(layout, ref.populations[0].particles.domain_particles, + cmp.populations[0].particles.domain_particles, diff); + + auto const freport = compare_tensor_fields(ref.populations[0].F, cmp.populations[0].F, diff); + PHARE_LOG_LINE_STR("results: " << freport.why()); + EXPECT_TRUE(freport); + + auto const rhoCport = compare_fields(ref.populations[0].rhoC, cmp.populations[0].rhoC, diff); + PHARE_LOG_LINE_STR("results: " << rhoCport.why()); + EXPECT_TRUE(rhoCport); + + auto const rhoPport = compare_fields(ref.populations[0].rhoP, cmp.populations[0].rhoP, diff); + PHARE_LOG_LINE_STR("results: " << rhoPport.why()); + EXPECT_TRUE(rhoPport); +} + +template +struct TestParam +{ + auto static constexpr opts = opts_; + auto static constexpr impl = impl_; + static constexpr auto dimension = opts.dimension; + static constexpr auto interp_order = opts.interp_order; +}; + + + + +template +struct UpdatersComparisonTest : public ::testing::Test +{ + auto constexpr static opts = Param::opts; + auto constexpr static dim = opts.dimension; + auto constexpr static interp = opts.interp_order; + + using PHARE_Types = core::PHARE_Types; + using GridLayout_t = TestGridLayout; + using RefParticleArray_t = ParticleArray; + using CmpParticleArray_t = RefParticleArray_t; + using UsableElectromag_t = UsableElectromag; + using RefIons_t = UsableIons; + using RefEM_t = UsableElectromag_t; + using CmpIons_t = UsableIons; + using CmpEM_t = UsableElectromag_t; + + GridLayout_t const layout{cells}; + UpdaterMode const updater_mode = UpdaterMode::all; + + UpdatersComparisonTest() {} + + void run() + { + cmp_patches.reserve(n_patches); + auto& ref = ref_patches.emplace_back(layout, make_ions(layout)); + + if (!ref_only) + for (std::size_t i = 0; i < n_patches; i++) + cmp_patches.emplace_back(layout, from_ions(layout, ref.ions)); + + if (!cmp_only) + ref_update(updater_mode, ref_patches); + + if (!ref_only) + cmp_update(updater_mode, cmp_patches); + + if (do_cmp) + for (auto& cmp : cmp_patches) + compare(*layout, ref_patches[0].ions, cmp.ions); + } + + + template + struct Patch + { + using ParticleArray_t = Ions_t::particle_array_type; + using GridLayout_t = UpdatersComparisonTest::GridLayout_t; + using Electromag_t = EM::Super; + + GridLayout_t layout; + Ions_t ions; + EM em{layout}; + Electromag_t electromag = *em; + + std::string patchID() const { return "patch_id"; } + }; + + + using RefPatch = Patch; + using CmpPatch = Patch; + + std::vector> ref_patches{}; + std::vector> cmp_patches{}; +}; + +// clang-format off +using Permutations_t = ::testing::Types< + TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam + , TestParam +>; +// clang-format on + +TYPED_TEST_SUITE(UpdatersComparisonTest, Permutations_t, ); + +TYPED_TEST(UpdatersComparisonTest, updater_domain_only) +{ + this->run(); +} + +template // used by gtest +void PrintTo(ParticleArray const& arr, std::ostream* os) +{ + *os << arr; +} + +} // namespace PHARE::core + + +int main(int argc, char** argv) +{ + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tests/core/numerics/pusher/test_pusher.cpp b/tests/core/numerics/pusher/test_pusher.cpp index 485126572..87520c456 100644 --- a/tests/core/numerics/pusher/test_pusher.cpp +++ b/tests/core/numerics/pusher/test_pusher.cpp @@ -1,32 +1,29 @@ + + +#include "core/utilities/types.hpp" +#include "core/utilities/box/box.hpp" +#include "core/utilities/range/range.hpp" +#include "core/data/particles/particle_array.hpp" +#include "core/numerics/boundary_condition/boundary_condition.hpp" +#include "core/numerics/ion_updater/ion_updater/ion_updater_impl0.hpp" + +#include "phare_core.hpp" + +#include "tests/core/data/gridlayout/test_gridlayout.hpp" + #include "gmock/gmock.h" #include "gtest/gtest.h" #include -#include -#include -#include -#include #include #include #include -#include #include #include #include -#include "core/data/particles/particle_array.hpp" -#include "core/numerics/boundary_condition/boundary_condition.hpp" -#include "core/numerics/pusher/boris.hpp" -#include "core/utilities/range/range.hpp" -#include "core/utilities/box/box.hpp" - -#include "tests/core/data/gridlayout/test_gridlayout.hpp" -#include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" - using namespace PHARE::core; - - struct Trajectory { std::vector x, y, z; @@ -97,8 +94,6 @@ class TestInterpolator return eb_interop; } - - void particleToMesh(auto&&... args) {} }; @@ -108,28 +103,6 @@ class TestElectromag { }; -template -struct TestIons // Some tests have large domains but no need for fields -{ - using particle_array_type = ParticleArray_; - - struct TestIonPop - { - auto& domainParticles() { return domain; } - auto& levelGhostParticles() { return levelGhost; } - auto& patchGhostParticles() { return patchGhost; } - double mass() const { return 1; } - - GridLayout const& layout; - ParticleArray_ domain{layout.AMRBox()}; - ParticleArray_ patchGhost{layout.AMRBox()}; - ParticleArray_ levelGhost{layout.AMRBox()}; - }; - - GridLayout layout; - TestIonPop pop{layout}; -}; - // with this mock, all particles are found inside class DummySelector @@ -147,76 +120,40 @@ class DummySelector template class APusher : public ::testing::Test { - static constexpr auto dimension = dim; - static constexpr auto interp_order = 1; - constexpr static PHARE::SimOpts opts{dimension, interp_order}; - - using PHARE_TYPES = PHARE::core::PHARE_Types; - using Interpolator = TestInterpolator; - using Electromag = TestElectromag; - using GridLayout_t = TestGridLayout; - using Ions_t = TestIons>; - using IonUpdater = typename PHARE::core::IonUpdater; - using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + using Particle = typename ParticleArray::Particle_t; + using PhareTypes = PHARE_Types; + using GridLayout_t = TestGridLayout; + using Ions_t = PhareTypes::Ions_t; public: - using Pusher_ = BorisPusher, Electromag, Interpolator, - BoundaryCondition, GridLayout_t>; + using Pusher = IonUpdater0::Pusher; APusher() - : expectedTrajectory{readExpectedTrajectory()} - , layout{30} - , pusher{std::make_unique()} - , mass{1} - , dt{0.0001} - , tstart{0} - , tend{10} - , nt{static_cast((tend - tstart) / dt + 1)} - { - particles.emplace_back( + particlesIn.emplace_back( + Particle{1., 1., ConstArray(5), ConstArray(0.), {0., 10., 0.}}); + particlesOut.emplace_back( Particle{1., 1., ConstArray(5), ConstArray(0.), {0., 10., 0.}}); - dxyz.fill(0.05); for (std::size_t i = 0; i < dim; i++) actual[i].resize(nt, 0.05); - pusher->setMeshAndTimeStep(dxyz, dt); - - std::transform(std::begin(dxyz), std::end(dxyz), std::begin(halfDtOverDl), - [dt = this->dt](double& x) { return 0.5 * dt / x; }); } protected: - using Particle = typename ParticleArray::Particle_t; - Trajectory expectedTrajectory; - GridLayout_t layout; - - std::unique_ptr pusher; - double mass; - double dt; - double tstart, tend; - std::size_t nt; - Electromag em; - Interpolator interpolator; + Trajectory expectedTrajectory{readExpectedTrajectory()}; + GridLayout_t layout{.05, 0ul}; + ParticleArray particlesIn{}; + ParticleArray particlesOut{}; + double mass = 1; + double dt = 0.0001; + double tstart = 0, tend = 10; + std::size_t nt = static_cast((tend - tstart) / dt + 1); + TestElectromag em; + TestInterpolator interpolator; + DummySelector selector; // BoundaryCondition bc; std::array, dim> actual; - std::array dxyz; - Ions_t ions{layout}; - - ParticleArray& particles = ions.pop.domainParticles(); - - std::array halfDtOverDl; - double const dto2m = 0.5 * dt / mass; - - void move() - { - for (auto& particle : particles) - { - particle.iCell = boris::advance(particle, halfDtOverDl); - boris::accelerate(particle, interpolator(particle, em, layout), dto2m); - particle.iCell = boris::advance(particle, halfDtOverDl); - } - } + std::array dxyz = layout.meshSize(); }; @@ -226,13 +163,21 @@ using APusher3D = APusher<3>; TEST_F(APusher3D, trajectoryIsOk) { + auto rangeIn = makeIndexRange(particlesIn); + auto rangeOut = makeIndexRange(particlesOut); + std::copy(rangeIn.begin(), rangeIn.end(), rangeOut.begin()); + + Pusher pusher{dt, mass, layout}; + for (decltype(nt) i = 0; i < nt; ++i) { - actual[0][i] = (particles[0].iCell[0] + particles[0].delta[0]) * dxyz[0]; - actual[1][i] = (particles[0].iCell[1] + particles[0].delta[1]) * dxyz[1]; - actual[2][i] = (particles[0].iCell[2] + particles[0].delta[2]) * dxyz[2]; + 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]; + actual[2][i] = (particlesOut[0].iCell[2] + particlesOut[0].delta[2]) * dxyz[2]; - move(); + pusher.move(rangeIn, rangeOut, em, interpolator, selector, selector); + + std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); } EXPECT_THAT(actual[0], ::testing::Pointwise(::testing::DoubleNear(1e-5), expectedTrajectory.x)); @@ -245,12 +190,19 @@ TEST_F(APusher3D, trajectoryIsOk) TEST_F(APusher2D, trajectoryIsOk) { + auto rangeIn = makeIndexRange(particlesIn); + auto rangeOut = makeIndexRange(particlesOut); + std::copy(rangeIn.begin(), rangeIn.end(), rangeOut.begin()); + + Pusher pusher{dt, mass, layout}; for (decltype(nt) i = 0; i < nt; ++i) { - actual[0][i] = (particles[0].iCell[0] + particles[0].delta[0]) * dxyz[0]; - actual[1][i] = (particles[0].iCell[1] + particles[0].delta[1]) * dxyz[1]; + 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]; - move(); + pusher.move(rangeIn, rangeOut, em, interpolator, selector, selector); + + std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); } EXPECT_THAT(actual[0], ::testing::Pointwise(::testing::DoubleNear(1e-5), expectedTrajectory.x)); @@ -261,11 +213,18 @@ TEST_F(APusher2D, trajectoryIsOk) TEST_F(APusher1D, trajectoryIsOk) { + auto rangeIn = makeIndexRange(particlesIn); + auto rangeOut = makeIndexRange(particlesOut); + std::copy(rangeIn.begin(), rangeIn.end(), rangeOut.begin()); + Pusher pusher{dt, mass, layout}; + for (decltype(nt) i = 0; i < nt; ++i) { - actual[0][i] = (particles[0].iCell[0] + particles[0].delta[0]) * dxyz[0]; + actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0]; - move(); + pusher.move(rangeIn, rangeOut, em, interpolator, selector, selector); + + std::copy(rangeOut.begin(), rangeOut.end(), rangeIn.begin()); } EXPECT_THAT(actual[0], ::testing::Pointwise(::testing::DoubleNear(1e-5), expectedTrajectory.x)); @@ -279,32 +238,14 @@ TEST_F(APusher1D, trajectoryIsOk) // and those that stay. class APusherWithLeavingParticles : public ::testing::Test { - using Interpolator = TestInterpolator; - using Electromag = TestElectromag; - - static constexpr auto dimension = 1; - static constexpr auto interp_order = 1; - constexpr static PHARE::SimOpts opts{dimension, interp_order}; - - using PHARE_TYPES = PHARE::core::PHARE_Types; - using ParticleArray_t = ParticleArray<1>; - using GridLayout_t = TestGridLayout; - using Ions_t = PHARE_TYPES::Ions_t; - using IonUpdater = typename PHARE::core::IonUpdater; - using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + using PhareTypes = PHARE_Types; + using GridLayout_t = TestGridLayout; + using Ions_t = PhareTypes::Ions_t; public: + using Pusher = IonUpdater0::Pusher; + APusherWithLeavingParticles() - : pusher{std::make_unique, TestElectromag, TestInterpolator, - BoundaryCondition<1, 1>, GridLayout_t>>()} - , mass{1} - , dt{0.001} - , tstart{0} - , tend{10} - , nt{static_cast((tend - tstart) / dt + 1)} - , layout{10} - // , particlesOut1{grow(domain, 1), 1000} - // , particlesOut2{grow(domain, 1), 1000} { std::random_device rd; std::mt19937 gen(rd()); @@ -312,38 +253,27 @@ class APusherWithLeavingParticles : public ::testing::Test std::uniform_real_distribution delta(0, 1); for (std::size_t iPart = 0; iPart < 1000; ++iPart) - { particlesIn.emplace_back( Particle<1>{1., 1., {{dis(gen)}}, {{delta(gen)}}, {{5., 0., 0.}}}); - } - pusher->setMeshAndTimeStep({{dx}}, dt); } protected: - std::unique_ptr, Electromag, Interpolator, - BoundaryCondition<1, 1>, GridLayout_t>> - pusher; - double mass; - double dt; - double tstart; - double tend; - std::size_t nt; - Electromag em; - Interpolator interpolator; + double mass = 1; + double dt = .001; + double tstart = 0; + double tend = 10; + std::size_t nt = static_cast((tend - tstart) / dt + 1); + TestElectromag em; + TestInterpolator interpolator; double dx = 0.1; - GridLayout_t layout; - Box domain = layout.AMRBox(); - // BoundaryCondition<1, 1> bc; - // ParticleArray<1> particlesOut1; - // ParticleArray<1> particlesOut2; - - // TestIons> ions{layout}; - UsableIons ions{layout}; - - ParticleArray& particlesIn = ions[0].domainParticles(); - - Boxing_t const boxing{layout, {layout.AMRBox()}}; + Box domain{Point{0.}, Point{1.}}; + Box cells{Point{0}, Point{9}}; + BoundaryCondition<1, 1> bc; + ParticleArray<1> particlesIn{grow(cells, 1)}; + ParticleArray<1> particlesOut1{{grow(cells, 1)}, 1000}; + ParticleArray<1> particlesOut2{{grow(cells, 1)}, 1000}; + GridLayout_t layout{10ul}; }; @@ -351,16 +281,33 @@ class APusherWithLeavingParticles : public ::testing::Test TEST_F(APusherWithLeavingParticles, splitLeavingFromNonLeavingParticles) { + auto rangeIn = makeIndexRange(particlesIn); + auto inDomain = rangeIn; + + auto selector = [this](auto& particleRange) // + { + auto& box = this->cells; + return particleRange.array().partition( + [&](auto const& cell) { return PHARE::core::isIn(Point{cell}, box); }); + }; + + Pusher pusher{dt, mass, layout}; for (decltype(nt) i = 0; i < nt; ++i) - pusher->move(ions[0], em, interpolator, boxing); - - auto& patchGhost = ions[0].patchGhostParticles(); - EXPECT_TRUE( - std::none_of(std::begin(patchGhost), std::end(patchGhost), - [this](auto const& particle) { return PHARE::core::isIn(particle, domain); })); - EXPECT_TRUE( - std::all_of(std::begin(particlesIn), std::end(particlesIn), - [this](auto const& particle) { return PHARE::core::isIn(particle, domain); })); + { + inDomain = pusher.move(rangeIn, rangeIn, em, interpolator, selector, selector); + + if (inDomain.end() != std::end(particlesIn)) + { + std::cout << inDomain.iend() << " and " << particlesIn.size() << "\n"; + break; + } + } + EXPECT_TRUE(std::none_of(inDomain.end(), std::end(particlesIn), [this](auto const& particle) { + return PHARE::core::isIn(Point{particle.iCell}, cells); + })); + EXPECT_TRUE(std::all_of(std::begin(inDomain), std::end(inDomain), [this](auto const& particle) { + return PHARE::core::isIn(Point{particle.iCell}, cells); + })); } @@ -386,9 +333,9 @@ TEST_F(APusherWithLeavingParticles, pusherWithOrWithoutBCReturnsSameNbrOfStaying { auto layout = DummyLayout<1>{}; newEndWithBC - = pusher->move(rangeIn, rangeOut1, em, mass, interpolator, selector, bc, layout); + = pusher.move(rangeIn, rangeOut1, em, mass, interpolator, selector, bc, layout); newEndWithoutBC - = pusher->move(rangeIn, rangeOut2, em, mass, interpolator, selector, layout); + = pusher.move(rangeIn, rangeOut2, em, mass, interpolator, selector, layout); if (newEndWithBC != std::end(particlesOut1) || newEndWithoutBC != std::end(particlesOut2)) { diff --git a/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp b/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp index b837af3fd..5e9e5424e 100644 --- a/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp +++ b/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp @@ -3,33 +3,36 @@ #include "tests/core/data/gridlayout/test_gridlayout.hpp" #include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" -using namespace PHARE; + template void updater_routine(benchmark::State& state) { + using namespace PHARE; + using namespace PHARE::core; + constexpr std::uint32_t cells = 30; constexpr std::uint32_t n_parts = 1e7; - auto static constexpr opts = PHARE::SimOpts{dim, interp}; + auto static constexpr opts = SimOpts{dim, interp}; - using PHARE_Types = core::PHARE_Types; + using PHARE_Types = PHARE_Types; using GridLayout_t = TestGridLayout; - using Electromag_t = core::UsableElectromag; + using Electromag_t = UsableElectromag; using ParticleArray = PHARE_Types::ParticleArray_t; using Particle_t = ParticleArray::value_type; - using Ions = PHARE::core::UsableIons; - using IonUpdater = core::IonUpdater; - using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + using Ions = UsableIons; + using IonUpdater = IonUpdaterProxy>; + using Boxing_t = UpdaterSelectionBoxing; GridLayout_t layout{cells}; Electromag_t em{layout}; - Ions ions{layout, "protons"}; + Ions ions{layout}; Boxing_t const boxing{layout, {grow(layout.AMRBox(), GridLayout_t::nbrParticleGhosts())}}; auto& patch_particles = ions.populations[0].particles; patch_particles.domain_particles.vector() - = std::vector(n_parts, core::bench::particle()); - core::bench::disperse(patch_particles.domain_particles, 0, cells - 1); + = std::vector(n_parts, bench::particle()); + bench::disperse(patch_particles.domain_particles, 0, cells - 1); std::sort(patch_particles.domain_particles); auto particles_copy = patch_particles.domain_particles; // tmp storage between update modes @@ -37,12 +40,12 @@ void updater_routine(benchmark::State& state) dict["pusher"]["name"] = std::string{"modified_boris"}; IonUpdater ionUpdater_{dict}; - double current_time = 1.0; - double new_time = 1.005; - auto dt = new_time - current_time; + double const current_time = 1.0; + double const new_time = 1.005; + auto const dt = new_time - current_time; while (state.KeepRunningBatch(1)) // while (state.KeepRunning()) { - ionUpdater_.updatePopulations(ions, em, boxing, dt, core::UpdaterMode::domain_only); + ionUpdater_.updatePopulations(UpdaterMode::domain_only, ions, em, boxing, dt); ionUpdater_.updateIons(ions); patch_particles.domain_particles = particles_copy; @@ -50,7 +53,8 @@ void updater_routine(benchmark::State& state) = std::get<4>(ions.getRunTimeResourcesViewList()[0].getCompileTimeResourcesViewList()); pack.setBuffer(&patch_particles.pack()); - ionUpdater_.updatePopulations(ions, em, boxing, dt, core::UpdaterMode::all); + + ionUpdater_.updatePopulations(UpdaterMode::all, ions, em, boxing, dt); ionUpdater_.updateIons(ions); } } diff --git a/tools/bench/core/numerics/pusher/push_raw_use.cpp b/tools/bench/core/numerics/pusher/push_raw_use.cpp index ae812f457..b0da0d4e5 100644 --- a/tools/bench/core/numerics/pusher/push_raw_use.cpp +++ b/tools/bench/core/numerics/pusher/push_raw_use.cpp @@ -1,45 +1,105 @@ +#include "core/numerics/ion_updater/ion_updater/ion_updater_impl0.hpp" +#define PHARE_UPDATER_IMPL 0 + #include "push_bench.hpp" + #include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" // Does not include google benchmark as we want to see only PHARE operations/instructions/etc + + +#if PHARE_UPDATER_IMPL == 0 + +template +void push() +{ + using namespace PHARE::core; + + auto static constexpr opts = PHARE::SimOpts{dim, interp}; + constexpr std::uint32_t cells = 65; + // constexpr std::uint32_t n_parts = 1e7; + + using PHARE_Types = PHARE_Types; + using GridLayout_t = TestGridLayout; + using Interpolator = Interpolator; + using Electromag_t = UsableElectromag; + using Ions_t = PHARE_Types::Ions_t; + using ParticleArray = Ions_t::particle_array_type; + using Pusher = IonUpdater0::Pusher; + // using ParticleRange = IndexRange; + // using BoundaryCondition = BoundaryCondition; + + GridLayout_t layout{cells}; + Electromag_t em{layout}; + + ParticleArray domainParticles{layout.AMRBox()}; + + std::stringstream ss; + ss << "unsorted_particles_" << dim << ".raw"; + bench::read_raw_from_file(domainParticles, ss.str()); + + // std::sort(domainParticles); + + ParticleArray tmpDomain{layout.AMRBox()}; + tmpDomain.vector() = domainParticles.vector(); + + auto rangeIn = makeIndexRange(domainParticles); + auto rangeOut = makeIndexRange(tmpDomain); + + Interpolator interpolator; + auto const no_op = [](auto& particleRange) { return particleRange; }; + + Pusher const pusher{/*dt=*/.001, /*mass=*/1, layout}; + pusher.move( + /*ParticleRange const&*/ rangeIn, // + /*ParticleRange& */ rangeOut, // + /*Electromag const& */ em, // + /*Interpolator&*/ interpolator, no_op, no_op); +} + +#endif + + +#if PHARE_UPDATER_IMPL == 1 + template void push() { + using namespace PHARE::core; + auto static constexpr opts = PHARE::SimOpts{dim, interp}; bool static constexpr copy_particle = true; 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 = PHARE::core::UsableElectromag; - using Ions_t = PHARE_Types::Ions_t; - using ParticleArray = Ions_t::particle_array_type; - using IonUpdater = typename PHARE::core::IonUpdater; - using Boxing_t = PHARE::core::UpdaterSelectionBoxing; - - using BorisPusher_t = PHARE::core::BorisPusher; + using PHARE_Types = PHARE_Types; + using GridLayout_t = TestGridLayout; + using Electromag_t = UsableElectromag; + using Interpolator = Interpolator; + using Ions_t = PHARE_Types::Ions_t; + using ParticleArray = Ions_t::particle_array_type; + using Pusher = IonUpdater1::Pusher; + using UsableIons_t = UsableIons; + using Boxing_t = PHARE::core::UpdaterSelectionBoxing; GridLayout_t layout{cells}; Electromag_t em{layout}; - Boxing_t boxing{layout, {layout.AMRBox()}}; - PHARE::core::UsableIons ions{layout}; + UsableIons_t ions{layout}; + Boxing_t const boxing{layout, {grow(layout.AMRBox(), GridLayout_t::nbrParticleGhosts())}}; std::stringstream ss; ss << "unsorted_particles_" << dim << ".raw"; - PHARE::core::bench::read_raw_from_file(ions[0].domainParticles(), ss.str()); - - BorisPusher_t pusher; - pusher.setMeshAndTimeStep(layout.meshSize(), .001); + bench::read_raw_from_file(ions[0].domainParticles(), ss.str()); Interpolator interpolator; - pusher.template move(ions[0], em, interpolator, boxing); + Pusher const pusher{/*dt=*/.001, /*mass=*/1, layout}; + pusher.move(ions[0].domainParticles(), em, boxing, interpolator); } +#endif + + + int main(int /*argc*/, char** /*argv*/) { push<3, 3>(); diff --git a/tools/bench/core/numerics/pusher/pusher.cpp b/tools/bench/core/numerics/pusher/pusher.cpp index 243c9321f..6d91e19ee 100644 --- a/tools/bench/core/numerics/pusher/pusher.cpp +++ b/tools/bench/core/numerics/pusher/pusher.cpp @@ -1,13 +1,17 @@ +#define PHARE_UPDATER_IMPL 1 + #include "push_bench.hpp" #include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" + +#if PHARE_UPDATER_IMPL == 0 + template void push(benchmark::State& state) { - auto static constexpr opts = PHARE::SimOpts{dim, interp}; - bool static constexpr copy_particle = true; - constexpr std::uint32_t cells = 65; - constexpr std::uint32_t n_parts = 1e7; + auto static constexpr opts = PHARE::SimOpts{dim, interp}; + constexpr std::uint32_t cells = 65; + constexpr std::uint32_t n_parts = 1e7; using PHARE_Types = PHARE::core::PHARE_Types; using GridLayout_t = TestGridLayout; @@ -17,37 +21,86 @@ void push(benchmark::State& state) using Ions_t = PHARE_Types::Ions_t; using ParticleArray = Ions_t::particle_array_type; using Particle_t = ParticleArray::value_type; + using ParticleRange = PHARE::core::IndexRange; - using IonUpdater = typename PHARE::core::IonUpdater; - using Boxing_t = PHARE::core::UpdaterSelectionBoxing; - - using BorisPusher_t = PHARE::core::BorisPusher; GridLayout_t layout{cells}; - Boxing_t boxing{layout, {layout.AMRBox()}}; - Electromag_t em{layout}; - PHARE::core::UsableIons ions{layout}; - auto& domainParticles = ions[0].domainParticles(); + ParticleArray domainParticles{layout.AMRBox()}; domainParticles.vector() = std::vector(n_parts, PHARE::core::bench::particle()); PHARE::core::bench::disperse(domainParticles, 0, cells - 1, 13337); + // std::sort(domainParticles); + ParticleArray tmpDomain{layout.AMRBox()}; + tmpDomain.vector() = domainParticles.vector(); + auto rangeIn = PHARE::core::makeIndexRange(domainParticles); + auto rangeOut = PHARE::core::makeIndexRange(tmpDomain); BorisPusher_t pusher; pusher.setMeshAndTimeStep(layout.meshSize(), .001); Interpolator interpolator; auto const no_op = [](auto& particleRange) { return particleRange; }; + while (state.KeepRunningBatch(1)) + { + pusher.move( + /*ParticleRange const&*/ rangeIn, // + /*ParticleRange& */ rangeOut, // + /*Electromag const& */ em, // + /*double mass */ 1, + /*Interpolator&*/ interpolator, + /*GridLayout const&*/ layout, // + no_op, no_op); + } +} + +#endif +#if PHARE_UPDATER_IMPL == 1 + +template +void push(benchmark::State& state) +{ + using namespace PHARE::core; + + auto static constexpr opts = PHARE::SimOpts{dim, interp}; + bool static constexpr copy_particle = true; + constexpr std::uint32_t cells = 65; + constexpr std::uint32_t n_parts = 1e7; + + using PHARE_Types = PHARE_Types; + using GridLayout_t = TestGridLayout; + using Electromag_t = UsableElectromag; + using Interpolator = PHARE::core::Interpolator; + using ParticleArray = PHARE_Types::ParticleArray_t; + using Particle_t = ParticleArray::value_type; + using Ions_t = PHARE_Types::Ions_t; + using Pusher = IonUpdater1::Pusher; + using Boxing_t = PHARE::core::UpdaterSelectionBoxing; + + GridLayout_t layout{cells}; + Electromag_t em{layout}; + UsableIons ions{layout}; + + auto& domainParticles = ions[0].domainParticles(); + domainParticles.vector() = std::vector(n_parts, bench::particle()); + bench::disperse(domainParticles, 0, cells - 1, 13337); + + Boxing_t const boxing{layout, {grow(layout.AMRBox(), GridLayout_t::nbrParticleGhosts())}}; + Interpolator interpolator; + Pusher pusher{.001, 1.0, layout}; while (state.KeepRunningBatch(1)) - pusher.template move(ions[0], em, interpolator, boxing); + pusher.template move_interpolate_and_sort(ions[0], em, boxing, interpolator); } +#endif + BENCHMARK_TEMPLATE(push, /*dim=*/1, /*interp=*/1)->Unit(benchmark::kMicrosecond); BENCHMARK_TEMPLATE(push, /*dim=*/1, /*interp=*/2)->Unit(benchmark::kMicrosecond); diff --git a/tools/bench/core/utilities/CMakeLists.txt b/tools/bench/core/utilities/CMakeLists.txt new file mode 100644 index 000000000..fcb3fb4e4 --- /dev/null +++ b/tools/bench/core/utilities/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required (VERSION 3.20.1) + +project(phare_bench_boxes) + +add_phare_cpp_benchmark(10 ${PROJECT_NAME} box_iterator ${CMAKE_CURRENT_BINARY_DIR}) +add_phare_cpp_benchmark(10 ${PROJECT_NAME} box_span ${CMAKE_CURRENT_BINARY_DIR}) +add_phare_cpp_benchmark(10 ${PROJECT_NAME} box_span_N_loop ${CMAKE_CURRENT_BINARY_DIR}) + diff --git a/tools/bench/core/utilities/bench_box.hpp b/tools/bench/core/utilities/bench_box.hpp new file mode 100644 index 000000000..d4f3ebd0b --- /dev/null +++ b/tools/bench/core/utilities/bench_box.hpp @@ -0,0 +1,39 @@ + +#include "core/utilities/types.hpp" +#include "core/utilities/box/box.hpp" +#include "core/utilities/box/box_span.hpp" +#include "core/data/field/field_box_span.hpp" + +#include "phare_core.hpp" +#include "phare_simulator_options.hpp" + +#include "tests/core/data/gridlayout/test_gridlayout.hpp" + +#include +#include + +#define PRINT(x) std::cout << __FILE__ << ":" << __LINE__ << " " << x << std::endl; + +namespace PHARE::bench::core +{ +static constexpr PHARE::SimOpts opts{3, 1}; +using PHARE_Types = PHARE::core::PHARE_Types; +using GridLayout_t = TestGridLayout; +using Grid_t = PHARE_Types::Grid_t; + +std::uint64_t static now() +{ + return std::chrono::duration_cast( + std::chrono::steady_clock::now().time_since_epoch()) + .count(); +} + +struct ScopeTimer +{ + ~ScopeTimer() { PRINT("Time: " << (now() - start) / div / 1e6 << " ms"); } + + double div = 1; + std::uint64_t const start = now(); +}; + +} // namespace PHARE::bench::core diff --git a/tools/bench/core/utilities/box_iterator.cpp b/tools/bench/core/utilities/box_iterator.cpp new file mode 100644 index 000000000..0b7aaeea7 --- /dev/null +++ b/tools/bench/core/utilities/box_iterator.cpp @@ -0,0 +1,45 @@ +#include "bench_box.hpp" +#include "core/def.hpp" + + + +namespace PHARE::bench::core +{ + +int test() +{ + std::size_t constexpr static FOR = 10; + + GridLayout_t layout{666}; + Grid_t dst{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 11}; + Grid_t src{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 11}; + auto const box = layout.ghostBoxFor(dst); + + { + ScopeTimer timer{FOR}; + for (std::size_t L = 0; L < FOR; ++L) + { + auto src_it = box.begin(); + auto dst_it = box.begin(); + for (; dst_it != box.end(); ++src_it, ++dst_it) + dst(*dst_it) += src(*src_it) + L; + } + } + + PHARE_DEBUG_DO({ + for (auto const bix : layout.ghostBoxFor(dst)) + if (dst(bix) != 11 * 11 + 45) + return 1; + }) + + return dst.data()[0] != dst.data()[dst.size() - 1]; +} + +} // namespace PHARE::bench::core + +int main() +{ + auto const ret = PHARE::bench::core::test(); + PRINT(ret); + return ret; +} diff --git a/tools/bench/core/utilities/box_span.cpp b/tools/bench/core/utilities/box_span.cpp new file mode 100644 index 000000000..d03652ff2 --- /dev/null +++ b/tools/bench/core/utilities/box_span.cpp @@ -0,0 +1,65 @@ +#include "bench_box.hpp" + + + +namespace PHARE::bench::core +{ + +int test() +{ + std::size_t constexpr static FOR = 10; + + GridLayout_t layout{666}; + Grid_t dst{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 11}; + Grid_t const src{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 11}; + auto const box = layout.ghostBoxFor(dst); + + { + ScopeTimer timer{FOR}; + for (std::size_t L = 0; L < FOR; ++L) + { + auto d_span = make_field_box_span(box, dst); + auto const s_span = make_field_box_span(box, src); + + auto d_slabs = d_span.begin(); + auto s_slabs = s_span.begin(); + auto const send = s_span.end(); + for (; s_slabs != send; ++s_slabs, ++d_slabs) + { + auto& s_slab = *s_slabs; + auto& d_slab = *d_slabs; + + auto d_rows = d_slab.begin(); + auto s_rows = s_slab.begin(); + auto const srend = s_slab.end(); + + for (; s_rows != srend; ++s_rows, ++d_rows) + { + auto& s_row = *s_rows; + auto& d_row = *d_rows; + + for (std::size_t i = 0; i < s_row.size(); ++i) + d_row[i] += s_row[i] + L; + } + } + } + } + + PHARE_DEBUG_DO({ + for (auto const bix : box) + if (dst(bix) != 11 * 11 + 45) + return 1; + }) + + return dst.data()[0] != dst.data()[dst.size() - 1]; +} + +} // namespace PHARE::bench::core + +int main() +{ + std::ios::sync_with_stdio(false); + auto const ret = PHARE::bench::core::test(); + PRINT(ret); + return ret; +} diff --git a/tools/bench/core/utilities/box_span_N_loop.cpp b/tools/bench/core/utilities/box_span_N_loop.cpp new file mode 100644 index 000000000..fa2546589 --- /dev/null +++ b/tools/bench/core/utilities/box_span_N_loop.cpp @@ -0,0 +1,85 @@ +#include "bench_box.hpp" +#include "core/def.hpp" + + + +namespace PHARE::bench::core +{ + +struct Opper +{ + std::size_t L, off; + double* dst; + double* src; + + void operator()(auto i) + { + auto const idx = off + i; + dst[idx] += src[idx] + L; + } +}; + +int test() +{ + std::size_t constexpr static FOR = 10; + std::size_t constexpr static N = 8; + + GridLayout_t layout{666}; // + Grid_t dst{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 11}; + Grid_t src{"rho", layout, PHARE::core::HybridQuantity::Scalar::rho, 11}; + + auto const box = layout.ghostBoxFor(dst); + auto const shape = box.shape(); + auto const size = shape[2]; + auto const mod = shape[2] & 100; + auto const times = shape[2] / 100; + + { + ScopeTimer timer{FOR}; + for (std::size_t L = 0; L < FOR; ++L) + { + auto tslabs = make_field_box_span(box, dst); + auto rslabs = make_field_box_span(box, src); + auto tslabit = tslabs.begin(); + + for (auto rslabit = rslabs.begin(); rslabit != rslabs.end(); ++rslabit, ++tslabit) + { + auto rrow = *rslabit; + auto trow = *tslabit; + auto trowit = trow.begin(); + + for (auto rrowit = rrow.begin(); rrowit != rrow.end(); ++rrowit, ++trowit) + { + std::size_t offset = 0; + + for (std::size_t i = 0; i < times; ++i) + { + Opper opper{L, offset, &(*trowit)[0], &(*rrowit)[0]}; + PHARE::core::for_N(opper); + offset += N; + } + + for (std::size_t i = times * N; i < size; ++i) + (*trowit)[i] += (*rrowit)[i] + L; + } + } + } + } + + PHARE_DEBUG_DO({ + for (auto const bix : layout.ghostBoxFor(dst)) + if (dst(bix) != 11 * 11 + 45) + return 1; + }) + + return dst.data()[0] != dst.data()[dst.size() - 1]; +} + +} // namespace PHARE::bench::core + +int main() +{ + auto const ret = PHARE::bench::core::test(); + PRINT(ret); + return ret; +} From ac332e473dc39a240e1a90773230463fe5872706 Mon Sep 17 00:00:00 2001 From: deegan Date: Mon, 30 Mar 2026 13:38:03 +0200 Subject: [PATCH 3/3] bench --- .../numerics/ion_updater/test_updaters.cpp | 10 +- tools/bench/core/bench.hpp | 1 - .../core/numerics/ion_updater/CMakeLists.txt | 9 +- .../ion_updater/bench_ion_updater.cpp | 76 +++---------- .../numerics/ion_updater/bench_updater.hpp | 100 ++++++++++++++++++ 5 files changed, 125 insertions(+), 71 deletions(-) create mode 100644 tools/bench/core/numerics/ion_updater/bench_updater.hpp diff --git a/tests/core/numerics/ion_updater/test_updaters.cpp b/tests/core/numerics/ion_updater/test_updaters.cpp index 4b50725a2..73be7ef18 100644 --- a/tests/core/numerics/ion_updater/test_updaters.cpp +++ b/tests/core/numerics/ion_updater/test_updaters.cpp @@ -15,9 +15,9 @@ namespace PHARE::core { // RUNTIME ENV VAR OVERRIDES +auto static const seeder = get_env_as("PHARE_SEED", std::size_t{0}); // 0 == NO SEED! auto static const cells = get_env_as("PHARE_CELLS", std::uint32_t{10}); auto static const ppc = get_env_as("PHARE_PPC", std::size_t{500}); -auto static const seeder = get_env_as("PHARE_SEED", std::size_t{0}); // 0 == NO SEED! auto static const n_patches = get_env_as("PHARE_PATCHES", std::size_t{1}); auto static const dt = get_env_as("PHARE_TIMESTEP", double{.001}); auto static const do_cmp = get_env_as("PHARE_COMPARE", std::size_t{1}); @@ -27,8 +27,8 @@ auto static const ref_only = get_env_as("PHARE_REF_ONLY", std::size_t{0}); bool static const premain = []() { return true; }(); -template -auto get_updater_for(Ions& /*ions*/, EM const& /*em*/, GridLayout_t const& /*layout*/) +template +auto get_updater_for(Ions const&) { if constexpr (impl == 0) return IonUpdaterProxy>(); @@ -50,7 +50,7 @@ void ref_update(UpdaterMode mode, Patches& patches) { auto const domainBox = layout.AMRBox(); auto const ghostBox = grow(domainBox, GridLayout_t::nbrParticleGhosts()); - auto updater = get_updater_for<0>(*ions, electromag, layout); + auto updater = get_updater_for<0>(*ions); Boxing_t const boxing{layout, remove(ghostBox, domainBox)}; updater.updatePopulations(mode, *ions, electromag, boxing, dt); @@ -68,7 +68,7 @@ void cmp_update(UpdaterMode mode, Patches& patches) { auto const domainBox = layout.AMRBox(); auto const ghostBox = grow(domainBox, GridLayout_t::nbrParticleGhosts()); - auto updater = get_updater_for(*ions, electromag, layout); + auto updater = get_updater_for(*ions); Boxing_t const boxing{layout, remove(ghostBox, domainBox)}; updater.updatePopulations(mode, *ions, electromag, boxing, dt); diff --git a/tools/bench/core/bench.hpp b/tools/bench/core/bench.hpp index 9a3419ba3..8480d7b47 100644 --- a/tools/bench/core/bench.hpp +++ b/tools/bench/core/bench.hpp @@ -8,7 +8,6 @@ #include "tests/core/data/vecfield/test_vecfield_fixtures.hpp" #include "tests/core/data/electromag/test_electromag_fixtures.hpp" -#include "benchmark/benchmark.h" #include diff --git a/tools/bench/core/numerics/ion_updater/CMakeLists.txt b/tools/bench/core/numerics/ion_updater/CMakeLists.txt index c033de426..dce74495a 100644 --- a/tools/bench/core/numerics/ion_updater/CMakeLists.txt +++ b/tools/bench/core/numerics/ion_updater/CMakeLists.txt @@ -2,4 +2,11 @@ cmake_minimum_required (VERSION 3.20.1) project(phare_bench_ion_updater) -add_phare_cpp_benchmark(10 ${PROJECT_NAME} bench_ion_updater ${CMAKE_CURRENT_BINARY_DIR}) +add_phare_cpp_benchmark(10 bench_updater_impl0 bench_ion_updater ${CMAKE_CURRENT_BINARY_DIR}) +target_compile_definitions(bench_updater_impl0_bench_ion_updater PRIVATE PHARE_UPDATER_IMPL=0) + +add_phare_cpp_benchmark(10 bench_updater_impl1 bench_ion_updater ${CMAKE_CURRENT_BINARY_DIR}) +target_compile_definitions(bench_updater_impl1_bench_ion_updater PRIVATE PHARE_UPDATER_IMPL=1) + +add_phare_cpp_benchmark(10 bench_updater_impl2 bench_ion_updater ${CMAKE_CURRENT_BINARY_DIR}) +target_compile_definitions(bench_updater_impl2_bench_ion_updater PRIVATE PHARE_UPDATER_IMPL=2) diff --git a/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp b/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp index 5e9e5424e..dc5051a33 100644 --- a/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp +++ b/tools/bench/core/numerics/ion_updater/bench_ion_updater.cpp @@ -1,75 +1,23 @@ -#include "tools/bench/core/bench.hpp" -#include "core/numerics/ion_updater/ion_updater.hpp" -#include "tests/core/data/gridlayout/test_gridlayout.hpp" -#include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" +#ifndef PHARE_UPDATER_IMPL +#error // no +#endif -template -void updater_routine(benchmark::State& state) -{ - using namespace PHARE; - using namespace PHARE::core; - - constexpr std::uint32_t cells = 30; - constexpr std::uint32_t n_parts = 1e7; - auto static constexpr opts = SimOpts{dim, interp}; - - using PHARE_Types = PHARE_Types; - using GridLayout_t = TestGridLayout; - using Electromag_t = UsableElectromag; - using ParticleArray = PHARE_Types::ParticleArray_t; - using Particle_t = ParticleArray::value_type; - using Ions = UsableIons; - using IonUpdater = IonUpdaterProxy>; - using Boxing_t = UpdaterSelectionBoxing; - - GridLayout_t layout{cells}; - Electromag_t em{layout}; - Ions ions{layout}; - Boxing_t const boxing{layout, {grow(layout.AMRBox(), GridLayout_t::nbrParticleGhosts())}}; - - auto& patch_particles = ions.populations[0].particles; - patch_particles.domain_particles.vector() - = std::vector(n_parts, bench::particle()); - bench::disperse(patch_particles.domain_particles, 0, cells - 1); - std::sort(patch_particles.domain_particles); - auto particles_copy = patch_particles.domain_particles; // tmp storage between update modes - - initializer::PHAREDict dict; - dict["pusher"]["name"] = std::string{"modified_boris"}; - IonUpdater ionUpdater_{dict}; - - double const current_time = 1.0; - double const new_time = 1.005; - auto const dt = new_time - current_time; - while (state.KeepRunningBatch(1)) // while (state.KeepRunning()) - { - ionUpdater_.updatePopulations(UpdaterMode::domain_only, ions, em, boxing, dt); - ionUpdater_.updateIons(ions); - - patch_particles.domain_particles = particles_copy; - auto& pack - = std::get<4>(ions.getRunTimeResourcesViewList()[0].getCompileTimeResourcesViewList()); - pack.setBuffer(&patch_particles.pack()); +#include "bench_updater.hpp" +namespace PHARE::core +{ - ionUpdater_.updatePopulations(UpdaterMode::all, ions, em, boxing, dt); - ionUpdater_.updateIons(ions); - } +BENCHMARK_TEMPLATE_DEFINE_F(UpdaterBencher, test_3_1_x, + Params<3, 1, PHARE_UPDATER_IMPL>)(benchmark::State& st) +{ + (*this)(st); } -BENCHMARK_TEMPLATE(updater_routine, 1, 1)->Unit(benchmark::kMicrosecond); -BENCHMARK_TEMPLATE(updater_routine, 1, 2)->Unit(benchmark::kMicrosecond); -BENCHMARK_TEMPLATE(updater_routine, 1, 3)->Unit(benchmark::kMicrosecond); - -BENCHMARK_TEMPLATE(updater_routine, 2, 1)->Unit(benchmark::kMicrosecond); -BENCHMARK_TEMPLATE(updater_routine, 2, 2)->Unit(benchmark::kMicrosecond); -BENCHMARK_TEMPLATE(updater_routine, 2, 3)->Unit(benchmark::kMicrosecond); +BENCHMARK_REGISTER_F(UpdaterBencher, test_3_1_x)->Unit(benchmark::kMicrosecond); -BENCHMARK_TEMPLATE(updater_routine, 3, 1)->Unit(benchmark::kMicrosecond); -BENCHMARK_TEMPLATE(updater_routine, 3, 2)->Unit(benchmark::kMicrosecond); -BENCHMARK_TEMPLATE(updater_routine, 3, 3)->Unit(benchmark::kMicrosecond); +} // namespace PHARE::core int main(int argc, char** argv) { diff --git a/tools/bench/core/numerics/ion_updater/bench_updater.hpp b/tools/bench/core/numerics/ion_updater/bench_updater.hpp new file mode 100644 index 000000000..5e9f27f30 --- /dev/null +++ b/tools/bench/core/numerics/ion_updater/bench_updater.hpp @@ -0,0 +1,100 @@ + +#include "core/numerics/ion_updater/ion_updater.hpp" + +#include "phare_core.hpp" +#include "phare_simulator_options.hpp" + +#include "tests/core/data/gridlayout/test_gridlayout.hpp" +#include "tests/core/data/ion_population/test_ion_population_fixtures.hpp" + +#include "tools/bench/core/bench.hpp" + +#include "benchmark/benchmark.h" + +namespace PHARE::core +{ + +template +auto get_updater_for(Ions const&) +{ + if constexpr (impl == 0) + return IonUpdaterProxy>(); + if constexpr (impl == 1) + return IonUpdaterProxy>(); + if constexpr (impl == 2) + return IonUpdaterProxy>(); +} + +template +struct Params +{ + static constexpr auto dim = dim_; + static constexpr auto interp = interp_; + static constexpr auto impl = impl_; +}; + +template +struct UpdaterBencher : public benchmark::Fixture +{ + static constexpr auto dim = Params::dim; + static constexpr auto interp = Params::interp; + static constexpr auto impl = Params::impl; + static constexpr std::uint32_t cells = 30; + static constexpr std::uint32_t n_parts = 1e7; + static constexpr auto opts = PHARE::SimOpts{dim, interp}; + + using PHARE_Types = core::PHARE_Types; + using GridLayout_t = TestGridLayout; + using Electromag_t = UsableElectromag; + using ParticleArray = PHARE_Types::ParticleArray_t; + using Particle_t = ParticleArray::value_type; + using Ions = UsableIons; + // using IonUpdater = IonUpdaterProxy>; + using Boxing_t = UpdaterSelectionBoxing; + + GridLayout_t layout{cells}; + Electromag_t em{layout}; + Ions ions{layout}; + Boxing_t const boxing{layout, {grow(layout.AMRBox(), GridLayout_t::nbrParticleGhosts())}}; + ParticleArray particles_copy; + + void SetUp(::benchmark::State& state) + { + auto& patch_particles = ions.populations[0].particles; + patch_particles.domain_particles.vector() + = std::vector(n_parts, bench::particle()); + bench::disperse(patch_particles.domain_particles, 0, cells - 1); + std::sort(patch_particles.domain_particles); + particles_copy = patch_particles.domain_particles; // tmp storage between update modes + } + + void TearDown(::benchmark::State& state) {} + + void operator()(benchmark::State& state) + { + initializer::PHAREDict dict; + dict["pusher"]["name"] = std::string{"modified_boris"}; + auto ionUpdater_ = get_updater_for(ions); + + auto& patch_particles = ions.populations[0].particles; + double const current_time = 1.0; + double const new_time = 1.005; + auto const dt = new_time - current_time; + while (state.KeepRunningBatch(1)) // while (state.KeepRunning()) + { + ionUpdater_.updatePopulations(UpdaterMode::domain_only, ions, em, boxing, dt); + ionUpdater_.updateIons(ions); + + patch_particles.domain_particles = particles_copy; + auto& pack = std::get<4>( + ions.getRunTimeResourcesViewList()[0].getCompileTimeResourcesViewList()); + pack.setBuffer(&patch_particles.pack()); + + + ionUpdater_.updatePopulations(UpdaterMode::all, ions, em, boxing, dt); + ionUpdater_.updateIons(ions); + } + } +}; + +} // namespace PHARE::core