Skip to content

Updater does not copy vector #1166

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
PhilipDeegan wants to merge 3 commits into
base: master
Choose a base branch
from
Draft

Updater does not copy vector #1166

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1a763a0
Updater does not copy vector
PhilipDeegan Feb 27, 2026
3a6856b
MORE!
PhilipDeegan Mar 24, 2026
ac332e4
bench
PhilipDeegan Mar 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 28 additions & 5 deletions src/amr/resources_manager/amr_utils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,29 +227,39 @@ namespace amr
noDomainOverlapsOn(hierarchy, iLevel);
}


// potentially to replace with SAMRAI coarse to fine boundary stuff
template<typename GridLayoutT> // 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<dimension>(domain);
auto const particleGhostBox = grow(domBox, GridLayoutT::nbrParticleGhosts());

auto const neighbors = getSameLevelNeighbors(patch, hierarchy);
auto const neighbors = getSameLevelNeighbors(patch, hierarchy);
std::vector<core::Box<int, GridLayoutT::dimension>> patchGhostLayerBoxes;
patchGhostLayerBoxes.reserve(neighbors.size() + 1);
patchGhostLayerBoxes.emplace_back(domBox);
for (auto const& neighbox : neighbors)
patchGhostLayerBoxes.emplace_back(
*(particleGhostBox * phare_box_from<dimension>(neighbox)));

return patchGhostLayerBoxes;
}

template<typename GridLayoutT>
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<dimension>(patch.getBox());
auto const particleGhostBox = grow(domBox, GridLayoutT::nbrParticleGhosts());
return core::generate(
[&](auto& neighbox) {
return *(particleGhostBox * phare_box_from<dimension>(neighbox));
},
getSameLevelNeighbors(patch, hierarchy));
}

inline auto to_string(auto const& id)
{
std::stringstream patchID;
Expand Down Expand Up @@ -340,7 +350,20 @@ namespace amr
orMissing(ilvl);
}

template<std::size_t dim>
auto removeIntersections(core::Box<int, dim> const& box, core::Box<int, dim> const& remove)
{
return generate([&](auto const& rem) { return phare_box_from<dim>(rem); },
samrai_box_from(box).removeIntersections(samrai_box_from(remove)));
}

template<std::size_t dim>
auto removeIntersections(core::Box<std::uint32_t, dim> const& box,
core::Box<uint32_t, dim> const& remove)
{
return generate([&](auto const& rem) { return as_unsigned_phare_box<dim>(rem); },
samrai_box_from(box).removeIntersections(samrai_box_from(remove)));
}

} // namespace amr
} // namespace PHARE
Expand Down
10 changes: 5 additions & 5 deletions src/amr/solvers/solver_ppc.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ class SolverPPC : public ISolver<AMR_Types>
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<Ions, Electromag, GridLayout>;
using IonUpdater_t = PHARE::core::IonUpdater<Ions>;

Electromag electromagPred_{"EMPred"};
Electromag electromagAvg_{"EMAvg"};
Expand Down Expand Up @@ -184,7 +184,7 @@ class SolverPPC : public ISolver<AMR_Types>
if (auto [it, suc] = levelBoxing.try_emplace(
amr::to_string(patch->getGlobalId()),
Boxing_t{amr::layoutFromPatch<GridLayout>(*patch),
amr::makeNonLevelGhostBoxFor<GridLayout>(*patch, hierarchy)});
amr::patchGhostBoxOverlaps<GridLayout>(*patch, hierarchy)});
!suc)
throw std::runtime_error("boxing map insertion failure");
}
Expand All @@ -198,7 +198,7 @@ class SolverPPC : public ISolver<AMR_Types>
}


using Boxing_t = core::UpdaterSelectionBoxing<IonUpdater_t, GridLayout>;
using Boxing_t = core::UpdaterSelectionBoxing<GridLayout, ParticleArray>;
std::unordered_map<int /*level*/, std::unordered_map<std::string /*patchid*/, Boxing_t>> boxing;


Expand Down Expand Up @@ -559,8 +559,8 @@ void SolverPPC<HybridModel, AMR_Types>::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)
{
Expand Down
7 changes: 2 additions & 5 deletions src/amr/utilities/box/amr_box.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
{
Expand Down
2 changes: 0 additions & 2 deletions src/core/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/core/data/grid/grid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"};
Expand Down
6 changes: 3 additions & 3 deletions src/core/data/grid/gridlayout.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ namespace core
static constexpr std::size_t interp_order = GridLayoutImpl::interp_order;
using This = GridLayout<GridLayoutImpl>;
using implT = GridLayoutImpl;

using AMRBox_t = Box<int, dimension>;

/**
* @brief Constructor of a GridLayout
Expand All @@ -109,8 +109,8 @@ namespace core
*/
GridLayout(std::array<double, dimension> const& meshSize,
std::array<std::uint32_t, dimension> const& nbrCells,
Point<double, dimension> const& origin,
Box<int, dimension> AMRBox = Box<int, dimension>{}, int level_number = 0)
Point<double, dimension> const& origin, AMRBox_t const AMRBox = {},
int level_number = 0)
: meshSize_{meshSize}
, origin_{origin}
, nbrPhysicalCells_{nbrCells}
Expand Down
10 changes: 5 additions & 5 deletions src/core/data/ions/ions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
6 changes: 6 additions & 0 deletions src/core/data/ndarray/ndarray_vector.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,12 @@ class NdArrayVector
return MaskedView{*this, std::forward<Mask>(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_; }
Expand Down
55 changes: 43 additions & 12 deletions src/core/data/particles/particle_array.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <vector>
#include <cstddef>
#include <utility>
#include <vector>

#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
{
Expand Down Expand Up @@ -41,7 +42,7 @@ class ParticleArray


public:
ParticleArray(box_t box)
ParticleArray(box_t box = {})
: box_{box}
, cellMap_{box_}
{
Expand Down Expand Up @@ -146,6 +147,27 @@ class ParticleArray
void empty_map() { cellMap_.empty(); }


template<typename Cell>
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;
Comment on lines +151 to +167
Copy link
Copy Markdown

@coderabbitai coderabbitai bot Feb 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Add explicit precondition guards in swap_last_reduce_by_one.

If this helper is ever called with an empty array or an out-of-range index, size() - 1 underflows and indexing is UB.

🔧 Proposed guard 
 template<typename Cell>
 bool swap_last_reduce_by_one(Cell const& oldCell, std::size_t const idx)
 {
+    assert(size() > 0);
+    assert(idx < size());
+
     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];
🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@src/core/data/particles/particle_array.hpp` around lines 151 - 167, Add
explicit precondition checks at the start of swap_last_reduce_by_one: verify
size() > 0 and that the provided idx is < size() before using size() - 1; if
either check fails, return false or handle as appropriate. Locate the checks
inside the swap_last_reduce_by_one function (which uses size(), particles_,
cellMap_.swap, cellMap_.erase and resize) so you avoid underflowing size() - 1
and prevent out-of-range access when calling particles_[size() - 1] and
performing cellMap_.swap/erase.

}


NO_DISCARD auto nbr_particles_in(box_t const& box) const { return cellMap_.size(box); }

using cell_t = std::array<int, dim>;
Expand Down Expand Up @@ -180,9 +202,18 @@ class ParticleArray


template<typename Cell>
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<typename Cell>
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())
{
Expand Down
1 change: 1 addition & 0 deletions src/core/errors.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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); }
Copy link
Copy Markdown

@coderabbitai coderabbitai bot Feb 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

🧩 Analysis chain

🏁 Script executed:

cat -n src/core/errors.hpp

Repository: PHAREHUB/PHARE

Length of output: 3070

🏁 Script executed:

#!/bin/bash
# Search for DictionaryException constructor call sites
rg -n 'DictionaryException\s*\(' --type=cpp --type=c --type=hpp --type=h -A 2

Repository: PHAREHUB/PHARE

Length of output: 85

🏁 Script executed:

#!/bin/bash
# Search for DictionaryException constructor call sites
rg -n 'DictionaryException\s*\(' --type=cpp --type=c --type=h -A 2

Repository: PHAREHUB/PHARE

Length of output: 402

🏁 Script executed:

#!/bin/bash
# Search for instantiation/throw patterns
rg -n 'throw\s+.*DictionaryException|DictionaryException\s*{' --type=cpp --type=c --type=h -B 1 -A 3

Repository: PHAREHUB/PHARE

Length of output: 227

🏁 Script executed:

#!/bin/bash
# Search for throw statements with DictionaryException
rg -n 'throw.*DictionaryException' --type=cpp --type=c --type=h -A 2

Repository: PHAREHUB/PHARE

Length of output: 1027

🏁 Script executed:

#!/bin/bash
# Search for usage of the two-arg constructor pattern
rg -n 'DictionaryException\s*\([^)]*,\s*[^)]*\)' --type=cpp --type=c --type=h -B 1 -A 1

Repository: PHAREHUB/PHARE

Length of output: 308

🏁 Script executed:

#!/bin/bash
# Get more context from the boris.hpp usage
sed -n '155,170p' src/core/numerics/pusher/boris.hpp

Repository: PHAREHUB/PHARE

Length of output: 646

Make the single-argument constructor explicit and type-aligned.

Line 23 uses a generic constructor signature (auto const&) but forwards to operator()(std::string const, std::string const), which only accepts strings. This mismatch obscures the actual API contract. Additionally, the non-explicit constructor allows unintended implicit conversions.

Suggested change 
-    DictionaryException(auto const& cause) { (*this)("cause", cause); }
+    explicit DictionaryException(std::string const& cause) { (*this)("cause", cause); }
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
DictionaryException(auto const& cause) { (*this)("cause", cause); }
explicit DictionaryException(std::string const& cause) { (*this)("cause", cause); }
🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@src/core/errors.hpp` at line 23, Make the single-argument constructor
explicit and match the parameter type expected by operator(): change
DictionaryException(auto const& cause) to an explicit constructor that takes a
string type (e.g., explicit DictionaryException(std::string const& cause) or
explicit DictionaryException(std::string_view cause) and convert to std::string)
and then call (*this)("cause", cause); to ensure the API contract and prevent
implicit conversions; update the constructor signature in the
DictionaryException declaration accordingly.


using Dict_t = cppdict::Dict<std::string>;

Expand Down
61 changes: 33 additions & 28 deletions src/core/numerics/interpolator/interpolator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -459,7 +459,6 @@
}



/**\brief interpolate electromagnetic fields on all particles in the range
*
* For each particle :
Expand All @@ -468,41 +467,47 @@
* - then it uses Interpol<> to calculate the interpolation of E and B components
* onto the particle.
*/
template<typename Particle_t, typename GridLayout, typename VecField, typename Field>
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();

Check notice

Code scanning / CodeQL

Unused local variable

Variable (unnamed local variable) is not used.

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable (unnamed local variable) is not used.

indexAndWeights_<QtyCentering, QtyCentering::primal>(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<typename ParticleRange, typename VecField, typename GridLayout, typename Field>
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_<QtyCentering, QtyCentering::primal>(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<typename ParticleRange, typename VecField, typename GridLayout, typename Field>
Expand Down
Loading
Loading