try reserve refinement

Author: PhilipDeegan

.github/workflows/cmake_ubuntu.yml

@@ -78,11 +78,11 @@ jobs:
         make -j2 && sudo make install && cd ../.. && rm -rf samrai
         cd ${{runner.workspace}}/build && rm -rf *
         cmake $GITHUB_WORKSPACE -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON --fresh \
-              -DCMAKE_BUILD_TYPE=RelWithDebInfo -Dasan=OFF \
+              -DCMAKE_BUILD_TYPE=Debug -Dasan=OFF \
               -DCMAKE_C_COMPILER_LAUNCHER=ccache \
               -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \
               -DlowResourceTests=ON -DdevMode=ON -Dbench=ON \
-              -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 " -Dphare_configurator=ON
+              -DCMAKE_CXX_FLAGS="-O3 -DPHARE_DIAG_DOUBLES=1 " -Dphare_configurator=ON
 
     - name: Build
       working-directory: ${{runner.workspace}}/build

src/amr/data/particles/refine/particles_data_split.hpp

@@ -2,22 +2,22 @@
 #define PHARE_PARTICLES_DATA_SPLIT_HPP
 
 
+#include "core/def.hpp"
+#include "phare_core.hpp"
 #include "core/def/phare_mpi.hpp"
+#include "core/utilities/constants.hpp"
 
-#include "core/def.hpp"
-#include "amr/data/particles/particles_data.hpp"
-#include "amr/resources_manager/amr_utils.hpp"
 #include "split.hpp"
-#include "core/utilities/constants.hpp"
-#include "phare_core.hpp"
 #include "amr/amr_constants.hpp"
+#include "amr/utilities/box/amr_box.hpp"
+#include "amr/resources_manager/amr_utils.hpp"
+#include "amr/data/particles/particles_data.hpp"
 
 #include <SAMRAI/geom/CartesianPatchGeometry.h>
 #include <SAMRAI/hier/Box.h>
 #include <SAMRAI/hier/RefineOperator.h>
 #include <SAMRAI/pdat/CellOverlap.h>
 
-#include <functional>
 
 
 namespace PHARE
@@ -53,6 +53,12 @@ namespace amr
     template<typename ParticleArray, ParticlesDataSplitType splitType, typename Splitter>
     class ParticlesRefineOperator : public SAMRAI::hier::RefineOperator
     {
+        using Particle_t = typename ParticleArray::value_type;
+        static constexpr bool putParticlesInCoarseBoundary
+            = splitType == ParticlesDataSplitType::coarseBoundary
+              || splitType == ParticlesDataSplitType::coarseBoundaryOld
+              || splitType == ParticlesDataSplitType::coarseBoundaryNew;
+
     public:
         static constexpr auto dim           = Splitter::dimension;
         static constexpr auto interpOrder   = Splitter::interp_order;
@@ -99,7 +105,7 @@ namespace amr
                 = std::dynamic_pointer_cast<ParticlesData<ParticleArray>>(
                     source.getPatchData(sourceComponent));
 
-            // Finnaly we need the cartesion geometry of both patch.
+            // Finaly we need the cartesion geometry of both patch.
             auto patchGeomDestination
                 = std::dynamic_pointer_cast<SAMRAI::geom::CartesianPatchGeometry>(
                     destination.getPatchGeometry());
@@ -133,9 +139,14 @@ namespace amr
             }
         }
 
+        static void _reserve(ParticleArray& particles, std::size_t const& size)
+        {
+            if constexpr (!putParticlesInCoarseBoundary)
+                particles.reserve(size);
+        }
 
         /** @brief given two ParticlesData (destination and source),
-         * an overlap , a ratio and the geometry of both patches, perform the
+         * an overlap, a ratio and the geometry of both patches, perform the
          * splitting of coarse particles onto the destination patch
          */
         void refine_(ParticlesData<ParticleArray>& destParticlesData,
@@ -166,43 +177,43 @@ namespace amr
             // same for destinationGhostBox and destinationDomainBox the later will allow to get an
             // index relative to the interior
 
+            std::array const particlesArrays{&srcInteriorParticles, &srcGhostParticles};
+
+            {
+                std::size_t incoming_estimate = 0;
+                for (auto const& destinationBox : destBoxes)
+                    for (auto const& sourceParticlesArray : particlesArrays)
+                        incoming_estimate += sourceParticlesArray->nbr_particles_in(
+                            phare_box_from<dim>(destinationBox));
+                _reserve(destDomainParticles, (incoming_estimate * nbRefinedPart));
+            }
+
             Splitter split;
 
             // The PatchLevelFillPattern had compute boxes that correspond to the expected filling.
             // In case of a coarseBoundary it will most likely give multiple boxes
             // in case of interior, this will be just one box usually
             for (auto const& destinationBox : destBoxes)
             {
-                std::array particlesArrays{&srcInteriorParticles, &srcGhostParticles};
-                auto splitBox = getSplitBox(destinationBox);
+                auto const splitBox = getSplitBox(destinationBox);
 
-                auto isInDest = [&destinationBox](auto const& particle) //
+                auto const isInDest = [&destinationBox](auto const& particle) //
                 { return isInBox(destinationBox, particle); };
 
-
                 for (auto const& sourceParticlesArray : particlesArrays)
                 {
                     for (auto const& particle : *sourceParticlesArray)
                     {
-                        std::array<typename ParticleArray::value_type, nbRefinedPart>
-                            refinedParticles;
-                        auto particleRefinedPos = toFineGrid<interpOrder>(particle);
+                        auto const particleRefinedPos = toFineGrid<interpOrder>(particle);
 
                         if (isInBox(splitBox, particleRefinedPos))
                         {
+                            std::array<Particle_t, nbRefinedPart> refinedParticles;
                             split(particleRefinedPos, refinedParticles);
 
-
                             // we need to know in which of interior or levelGhostParticlesXXXX
                             // arrays we must put particles
 
-                            bool constexpr putParticlesInCoarseBoundary
-                                = splitType == ParticlesDataSplitType::coarseBoundary
-                                  || splitType == ParticlesDataSplitType::coarseBoundaryOld
-                                  || splitType == ParticlesDataSplitType::coarseBoundaryNew;
-
-
-
                             if constexpr (putParticlesInCoarseBoundary)
                             {
                                 if constexpr (splitType == ParticlesDataSplitType::coarseBoundary)
@@ -243,9 +254,9 @@ namespace amr
                                              std::back_inserter(destDomainParticles), isInDest);
                             }
                         } // end is candidate for split
-                    }     // end loop on particles
-                }         // end loop on source particle arrays
-            }             // loop on destination box
+                    } // end loop on particles
+                } // end loop on source particle arrays
+            } // loop on destination box
         }
 
 

src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp

@@ -44,7 +44,7 @@ class MaxwellianParticleInitializer : public ParticleInitializer<ParticleArray,
         std::uint32_t const& nbrParticlesPerCell, std::optional<std::size_t> seed = {},
         Basis const basis                                 = Basis::Cartesian,
         std::array<InputFunction, 3> const& magneticField = {nullptr, nullptr, nullptr},
-        double densityCutOff                              = 1e-5)
+        double const densityCutOff = 1e-5, double const over_alloc_factor = .1)
         : density_{density}
         , bulkVelocity_{bulkVelocity}
         , thermalVelocity_{thermalVelocity}
@@ -54,6 +54,7 @@ class MaxwellianParticleInitializer : public ParticleInitializer<ParticleArray,
         , nbrParticlePerCell_{nbrParticlesPerCell}
         , basis_{basis}
         , rngSeed_{seed}
+        , over_alloc_factor_{over_alloc_factor}
     {
     }
 
@@ -92,6 +93,7 @@ class MaxwellianParticleInitializer : public ParticleInitializer<ParticleArray,
     std::uint32_t nbrParticlePerCell_;
     Basis basis_;
     std::optional<std::size_t> rngSeed_;
+    double over_alloc_factor_ = .1;
 };
 
 
@@ -178,20 +180,16 @@ void MaxwellianParticleInitializer<ParticleArray, GridLayout>::loadParticles(
     auto randGen           = getRNG(rngSeed_);
     ParticleDeltaDistribution<double> deltaDistrib;
 
-    std::size_t const expected_size
-        = [&](auto const& _n /* or error: reference to local binding 'n' declared*/) {
-              std::size_t size = 0;
-              for (std::size_t flatCellIdx = 0; flatCellIdx < ndCellIndices.size(); ++flatCellIdx)
-                  if (_n[flatCellIdx] >= densityCutOff_)
-                      ++size;
-              return size * nbrParticlePerCell_;
-          }(n);
-
-    double const over_alloc_factor = .1; // todo from dict
-    std::size_t const allocate_size
-        = expected_size
-          + (layout.AMRBox().surface_cell_count() * (nbrParticlePerCell_ * over_alloc_factor));
-    particles.reserve(allocate_size);
+    auto const expected_size = std::accumulate(n, n + ndCellIndices.size(), std::size_t{0},
+                                               [&](auto const& sum, auto const& density_value) {
+                                                   return (density_value > densityCutOff_)
+                                                              ? sum + nbrParticlePerCell_
+                                                              : sum;
+                                               });
+
+    auto const incoming_estimate
+        = (layout.AMRBox().surface_cell_count() * (nbrParticlePerCell_ * over_alloc_factor_));
+    particles.reserve(incoming_estimate);
 
     for (std::size_t flatCellIdx = 0; flatCellIdx < ndCellIndices.size(); ++flatCellIdx)
     {