Skip to content
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

eb-less particles - mesh to particle eb handled directory in boris #781

Merged
merged 4 commits into from
Nov 25, 2023
Merged

eb-less particles - mesh to particle eb handled directory in boris #781

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3e92405
push benchmark updates
PhilipDeegan Nov 18, 2023
523a5dd
ebless particles
PhilipDeegan Nov 19, 2023
bb34fe5
expect_true(abs)) -> expect_near
PhilipDeegan Nov 21, 2023
09746ff
pr comments /cleanup
PhilipDeegan Nov 23, 2023
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
21 changes: 5 additions & 16 deletions src/core/data/particles/particle.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,26 +61,17 @@ struct Particle
std::array<double, dim> delta = ConstArray<double, dim>();
std::array<double, 3> v = ConstArray<double, 3>();

double Ex = 0, Ey = 0, Ez = 0;
double Bx = 0, By = 0, Bz = 0;

NO_DISCARD bool operator==(Particle<dim> const& that) const
{
return (this->weight == that.weight) && //
(this->charge == that.charge) && //
(this->iCell == that.iCell) && //
(this->delta == that.delta) && //
(this->v == that.v) && //
(this->Ex == that.Ex) && //
(this->Ey == that.Ey) && //
(this->Ez == that.Ez) && //
(this->Bx == that.Bx) && //
(this->By == that.By) && //
(this->Bz == that.Bz);
(this->v == that.v);
}

template<std::size_t dimension>
friend std::ostream& operator<<(std::ostream& out, const Particle<dimension>& particle);
friend std::ostream& operator<<(std::ostream& out, Particle<dimension> const& particle);
};

template<std::size_t dim>
Expand All @@ -102,8 +93,6 @@ std::ostream& operator<<(std::ostream& out, Particle<dim> const& particle)
out << v << ",";
}
out << "), charge : " << particle.charge << ", weight : " << particle.weight;
out << ", Exyz : " << particle.Ex << "," << particle.Ey << "," << particle.Ez;
out << ", Bxyz : " << particle.Bx << "," << particle.By << "," << particle.Bz;
out << '\n';
return out;
}
nicolasaunai marked this conversation as resolved.
Show resolved Hide resolved
Expand Down Expand Up @@ -132,9 +121,9 @@ inline constexpr auto is_phare_particle_type

template<std::size_t dim, template<std::size_t> typename ParticleA,
template<std::size_t> typename ParticleB>
NO_DISCARD typename std::enable_if_t<is_phare_particle_type<dim, ParticleA<dim>>
and is_phare_particle_type<dim, ParticleB<dim>>,
bool>
NO_DISCARD typename std::enable_if_t<
is_phare_particle_type<dim, ParticleA<dim>> and is_phare_particle_type<dim, ParticleB<dim>>,
bool>
operator==(ParticleA<dim> const& particleA, ParticleB<dim> const& particleB)
{
return particleA.weight == particleB.weight and //
Expand Down
2 changes: 2 additions & 0 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 @@ -231,6 +231,8 @@ class ParticleArray
NO_DISCARD auto& vector() { return particles_; }
NO_DISCARD auto& vector() const { return particles_; }

auto& box() const { return box_; }

private:
Vector particles_;
box_t box_;
Expand Down
79 changes: 44 additions & 35 deletions src/core/numerics/interpolator/interpolator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,12 +9,16 @@
#include "core/data/grid/gridlayout.hpp"
#include "core/data/vecfield/vecfield_component.hpp"
#include "core/utilities/point/point.hpp"
#include "core/def.hpp"

#include "core/def.hpp"
#include "core/logger.hpp"



namespace PHARE::core
{


//! return the size of the index and weights arrays for
//! interpolation at a given order
// Number of points where the interpolator of order interpOrder
nicolasaunai marked this conversation as resolved.
Show resolved Hide resolved
Expand Down Expand Up @@ -457,26 +461,22 @@
public:
auto static constexpr interp_order = interpOrder;
auto static constexpr dimension = dim;
/**\brief interpolate electromagnetic fields on all particles in the range

/**\brief interpolate electromagnetic fields on a particle and return particles EB
*
* For each particle :
* - The function first calculates the startIndex and weights for interpolation at
* order InterpOrder and in dimension dim for dual and primal nodes
* - then it uses Interpol<> to calculate the interpolation of E and B components
* onto the particle.
*/
template<typename ParticleRange, typename Electromag, typename GridLayout>
inline void operator()(ParticleRange&& particleRange, Electromag const& Em,
GridLayout const& layout)
template<typename Particle_t, typename Electromag, typename GridLayout>
inline auto operator()(Particle_t& currPart, Electromag const& Em, GridLayout const& layout)
{
PHARE_LOG_SCOPE("Interpolator::operator()");

auto begin = particleRange.begin();
auto end = particleRange.end();

using Scalar = HybridQuantity::Scalar;
auto const& [Ex, Ey, Ez] = Em.E();
auto const& [Bx, By, Bz] = Em.B();
using E_B_tuple = std::tuple<std::array<double, 3>, std::array<double, 3>>;
nicolasaunai marked this conversation as resolved.
Show resolved Hide resolved
using Scalar = HybridQuantity::Scalar;

// for each particle, first calculate the startIndex and weights for dual and
// primal quantities. then, knowing the centering (primal or dual) of each
Expand All @@ -485,33 +485,36 @@
// calculated twice, and not for each E,B component.

PHARE_LOG_START("MeshToParticle::operator()");
for (auto currPart = begin; currPart != end; ++currPart)
{
auto& iCell = currPart->iCell;
auto& delta = currPart->delta;
indexAndWeights_<QtyCentering, QtyCentering::dual>(layout, iCell, delta);
indexAndWeights_<QtyCentering, QtyCentering::primal>(layout, iCell, delta);

auto indexWeights = std::forward_as_tuple(dual_startIndex_, dual_weights_,
primal_startIndex_, primal_weights_);

currPart->Ex
= meshToParticle_.template operator()<GridLayout, Scalar::Ex>(Ex, indexWeights);
currPart->Ey
= meshToParticle_.template operator()<GridLayout, Scalar::Ey>(Ey, indexWeights);
currPart->Ez
= meshToParticle_.template operator()<GridLayout, Scalar::Ez>(Ez, indexWeights);
currPart->Bx
= meshToParticle_.template operator()<GridLayout, Scalar::Bx>(Bx, indexWeights);
currPart->By
= meshToParticle_.template operator()<GridLayout, Scalar::By>(By, indexWeights);
currPart->Bz
= meshToParticle_.template operator()<GridLayout, Scalar::Bz>(Bz, indexWeights);
}

auto& iCell = currPart.iCell;
auto& delta = currPart.delta;
indexAndWeights_<QtyCentering, QtyCentering::dual>(layout, iCell, delta);
indexAndWeights_<QtyCentering, QtyCentering::primal>(layout, iCell, delta);

auto indexWeights = std::forward_as_tuple(dual_startIndex_, dual_weights_,
primal_startIndex_, primal_weights_);

E_B_tuple particle_EB;
auto& [pE, pB] = particle_EB;
auto& [pEx, pEy, pEz] = pE;
auto& [pBx, pBy, pBz] = pB;

auto const& [Ex, Ey, Ez] = Em.E();

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable (unnamed local variable) is not used.
auto const& [Bx, By, Bz] = Em.B();

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable (unnamed local variable) is not used.

pEx = meshToParticle_.template operator()<GridLayout, Scalar::Ex>(Ex, indexWeights);
nicolasaunai marked this conversation as resolved.
Show resolved Hide resolved
pEy = meshToParticle_.template operator()<GridLayout, Scalar::Ey>(Ey, indexWeights);
pEz = meshToParticle_.template operator()<GridLayout, Scalar::Ez>(Ez, indexWeights);
pBx = meshToParticle_.template operator()<GridLayout, Scalar::Bx>(Bx, indexWeights);
pBy = meshToParticle_.template operator()<GridLayout, Scalar::By>(By, indexWeights);
pBz = meshToParticle_.template operator()<GridLayout, Scalar::Bz>(Bz, indexWeights);

PHARE_LOG_STOP("MeshToParticle::operator()");
return particle_EB;
}



/**\brief interpolate electromagnetic fields on all particles in the range
*
* For each particle :
Expand All @@ -521,7 +524,7 @@
* onto the particle.
*/
template<typename ParticleRange, typename VecField, typename GridLayout, typename Field>
inline void operator()(ParticleRange&& particleRange, Field& density, VecField& flux,
inline void operator()(ParticleRange& particleRange, Field& density, VecField& flux,
GridLayout const& layout, double coef = 1.)
{
auto begin = particleRange.begin();
nicolasaunai marked this conversation as resolved.
Show resolved Hide resolved
Expand All @@ -548,6 +551,12 @@
}
PHARE_LOG_STOP("ParticleToMesh::operator()");
}
template<typename ParticleRange, typename VecField, typename GridLayout, typename Field>
inline void operator()(ParticleRange&& range, Field& density, VecField& flux,
GridLayout const& layout, double coef = 1.)
{
(*this)(range, density, flux, layout, coef);
}


/**
nicolasaunai marked this conversation as resolved.
Show resolved Hide resolved
Comment on lines 561 to 562
Copy link

Choose a reason for hiding this comment

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

The comment block here is not updated to reflect the changes in the code. Since the Interpolator class now operates on individual particles rather than ranges, the comment should be updated to reflect this change. This is important for maintainability and to avoid confusion for future developers who will work with this code.

Expand Down
Loading
Loading