Skip to content

Commit

Permalink
ebless particles
Browse files Browse the repository at this point in the history
  • Loading branch information
@PhilipDeegan
PhilipDeegan committed Nov 18, 2023
1 parent be0eef0 commit fe2cb2a
Show file tree
Hide file tree
Showing 11 changed files with 192 additions and 255 deletions.
21 changes: 5 additions & 16 deletions src/core/data/particles/particle.hpp
View file
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;
}
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
79 changes: 44 additions & 35 deletions src/core/numerics/interpolator/interpolator.hpp
View file
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
Expand Down Expand Up @@ -457,26 +461,22 @@ class Interpolator : private Weighter<interpOrder>
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>>;
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 @@ class Interpolator : private Weighter<interpOrder>
// 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();
auto const& [Bx, By, Bz] = Em.B();

pEx = meshToParticle_.template operator()<GridLayout, Scalar::Ex>(Ex, indexWeights);
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 @@ class Interpolator : private Weighter<interpOrder>
* 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();
Expand All @@ -548,6 +551,12 @@ class Interpolator : private Weighter<interpOrder>
}
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);
}


/**
Expand Down
110 changes: 56 additions & 54 deletions src/core/numerics/pusher/boris.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@

namespace PHARE::core
{


template<std::size_t dim, typename ParticleRange, typename Electromag, typename Interpolator,
typename BoundaryCondition, typename GridLayout>
class BorisPusher
Expand Down Expand Up @@ -91,10 +93,14 @@ class BorisPusher

// get electromagnetic fields interpolated on the particles of rangeOut
// stop at newEnd.
interpolator(rangeOut, emFields, layout);
dto2m = 0.5 * dt_ / mass;
for (auto& currPart : rangeOut)
{
auto const& particle_EB = interpolator(currPart, emFields, layout);

// get the particle velocity from t=n to t=n+1
accelerate_(rangeOut, rangeOut, mass);
// get the particle velocity from t=n to t=n+1
accelerate_(currPart, particle_EB);
}

// now advance the particles from t=n+1/2 to t=n+1 using v_{n+1} just calculated
// and get a pointer to the first leaving particle
Expand Down Expand Up @@ -148,7 +154,7 @@ class BorisPusher
* @return the function returns and iterator on the first leaving particle, as
* detected by the ParticleSelector
*/
void pushStep_(ParticleRange const& rangeIn, ParticleRange& rangeOut, PushStep step)
void pushStep_(ParticleRange const& rangeIn, ParticleRange& rangeOut, PushStep const step)
{
auto& inParticles = rangeIn.array();
auto& outParticles = rangeOut.array();
Expand Down Expand Up @@ -180,79 +186,75 @@ class BorisPusher

/** Accelerate the particles in rangeIn and put the new velocity in rangeOut
*/
void accelerate_(ParticleRange rangeIn, ParticleRange rangeOut, double mass)
template<typename Particle_t, typename ParticleEB>
void accelerate_(Particle_t& part, ParticleEB const& particleEB)
{
double dto2m = 0.5 * dt_ / mass;
auto& [pE, pB] = particleEB;
auto& [pEx, pEy, pEz] = pE;
auto& [pBx, pBy, pBz] = pB;

auto& inParticles = rangeIn.array();
auto& outParticles = rangeOut.array();

for (auto inIdx = rangeIn.ibegin(), outIdx = rangeOut.ibegin(); inIdx < rangeIn.iend();
++inIdx, ++outIdx)
{
auto& inPart = inParticles[inIdx];
auto& outPart = outParticles[inIdx];
double coef1 = inPart.charge * dto2m;
double const coef1 = part.charge * dto2m;

// We now apply the 3 steps of the BORIS PUSHER
// We now apply the 3 steps of the BORIS PUSHER

// 1st half push of the electric field
double velx1 = inPart.v[0] + coef1 * inPart.Ex;
double vely1 = inPart.v[1] + coef1 * inPart.Ey;
double velz1 = inPart.v[2] + coef1 * inPart.Ez;
// 1st half push of the electric field
double velx1 = part.v[0] + coef1 * pEx;
double vely1 = part.v[1] + coef1 * pEy;
double velz1 = part.v[2] + coef1 * pEz;


// preparing variables for magnetic rotation
double const rx = coef1 * inPart.Bx;
double const ry = coef1 * inPart.By;
double const rz = coef1 * inPart.Bz;
// preparing variables for magnetic rotation
double const rx = coef1 * pBx;
double const ry = coef1 * pBy;
double const rz = coef1 * pBz;

double const rx2 = rx * rx;
double const ry2 = ry * ry;
double const rz2 = rz * rz;
double const rxry = rx * ry;
double const rxrz = rx * rz;
double const ryrz = ry * rz;
double const rx2 = rx * rx;
double const ry2 = ry * ry;
double const rz2 = rz * rz;
double const rxry = rx * ry;
double const rxrz = rx * rz;
double const ryrz = ry * rz;

double const invDet = 1. / (1. + rx2 + ry2 + rz2);
double const invDet = 1. / (1. + rx2 + ry2 + rz2);

// preparing rotation matrix due to the magnetic field
// m = invDet*(I + r*r - r x I) - I where x denotes the cross product
double const mxx = 1. + rx2 - ry2 - rz2;
double const mxy = 2. * (rxry + rz);
double const mxz = 2. * (rxrz - ry);
// preparing rotation matrix due to the magnetic field
// m = invDet*(I + r*r - r x I) - I where x denotes the cross product
double const mxx = 1. + rx2 - ry2 - rz2;
double const mxy = 2. * (rxry + rz);
double const mxz = 2. * (rxrz - ry);

double const myx = 2. * (rxry - rz);
double const myy = 1. + ry2 - rx2 - rz2;
double const myz = 2. * (ryrz + rx);
double const myx = 2. * (rxry - rz);
double const myy = 1. + ry2 - rx2 - rz2;
double const myz = 2. * (ryrz + rx);

double const mzx = 2. * (rxrz + ry);
double const mzy = 2. * (ryrz - rx);
double const mzz = 1. + rz2 - rx2 - ry2;
double const mzx = 2. * (rxrz + ry);
double const mzy = 2. * (ryrz - rx);
double const mzz = 1. + rz2 - rx2 - ry2;

// magnetic rotation
double const velx2 = (mxx * velx1 + mxy * vely1 + mxz * velz1) * invDet;
double const vely2 = (myx * velx1 + myy * vely1 + myz * velz1) * invDet;
double const velz2 = (mzx * velx1 + mzy * vely1 + mzz * velz1) * invDet;
// magnetic rotation
double const velx2 = (mxx * velx1 + mxy * vely1 + mxz * velz1) * invDet;
double const vely2 = (myx * velx1 + myy * vely1 + myz * velz1) * invDet;
double const velz2 = (mzx * velx1 + mzy * vely1 + mzz * velz1) * invDet;


// 2nd half push of the electric field
velx1 = velx2 + coef1 * inPart.Ex;
vely1 = vely2 + coef1 * inPart.Ey;
velz1 = velz2 + coef1 * inPart.Ez;
// 2nd half push of the electric field
velx1 = velx2 + coef1 * pEx;
vely1 = vely2 + coef1 * pEy;
velz1 = velz2 + coef1 * pEz;

// Update particle velocity
outPart.v[0] = velx1;
outPart.v[1] = vely1;
outPart.v[2] = velz1;
}
// Update particle velocity
part.v[0] = velx1;
part.v[1] = vely1;
part.v[2] = velz1;
}




std::array<double, dim> halfDtOverDl_;
double dt_;
double dto2m;
};

} // namespace PHARE::core
Expand Down
12 changes: 0 additions & 12 deletions tests/amr/data/particles/copy/test_particledata_copyNd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,12 +117,6 @@ TYPED_TEST(AParticlesDataND, PreservesAllParticleAttributesAfterCopy)
Pointwise(DoubleEq(), this->particle.delta));
EXPECT_THAT(this->destData.domainParticles[0].weight, DoubleEq(this->particle.weight));
EXPECT_THAT(this->destData.domainParticles[0].charge, DoubleEq(this->particle.charge));
EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Ex, this->particle.Ex);
EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Ey, this->particle.Ey);
EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Ez, this->particle.Ez);
EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Bx, this->particle.Bx);
EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].By, this->particle.By);
EXPECT_DOUBLE_EQ(this->destData.domainParticles[0].Bz, this->particle.Bz);

// particle is in the domain of the source patchdata
// and in last ghost of the destination patchdata
Expand All @@ -138,12 +132,6 @@ TYPED_TEST(AParticlesDataND, PreservesAllParticleAttributesAfterCopy)
Pointwise(DoubleEq(), this->particle.delta));
EXPECT_THAT(this->destData.patchGhostParticles[0].weight, DoubleEq(this->particle.weight));
EXPECT_THAT(this->destData.patchGhostParticles[0].charge, DoubleEq(this->particle.charge));
EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Ex, this->particle.Ex);
EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Ey, this->particle.Ey);
EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Ez, this->particle.Ez);
EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Bx, this->particle.Bx);
EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].By, this->particle.By);
EXPECT_DOUBLE_EQ(this->destData.patchGhostParticles[0].Bz, this->particle.Bz);
}


Expand Down
6 changes: 0 additions & 6 deletions tests/amr/data/particles/copy_overlap/test_particledata_copy_periodicNd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,12 +135,6 @@ TYPED_TEST(twoParticlesDataNDTouchingPeriodicBorders, preserveParticleAttributes
EXPECT_THAT(this->destPdat.patchGhostParticles[0].delta, Eq(this->particle.delta));
EXPECT_THAT(this->destPdat.patchGhostParticles[0].weight, Eq(this->particle.weight));
EXPECT_THAT(this->destPdat.patchGhostParticles[0].charge, Eq(this->particle.charge));
EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Ex, this->particle.Ex);
EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Ey, this->particle.Ey);
EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Ez, this->particle.Ez);
EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Bx, this->particle.Bx);
EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].By, this->particle.By);
EXPECT_DOUBLE_EQ(this->destPdat.patchGhostParticles[0].Bz, this->particle.Bz);
}


Expand Down
Loading

0 comments on commit fe2cb2a

Please sign in to comment.