aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles/PhysicalParticleContainer.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles/PhysicalParticleContainer.cpp')
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp123
1 files changed, 108 insertions, 15 deletions
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index ec994cbb7..6cfe20dd9 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -17,6 +17,8 @@
#include <UpdatePosition.H>
#include <UpdateMomentumBoris.H>
#include <UpdateMomentumVay.H>
+#include <UpdateMomentumBorisWithRadiationReaction.H>
+#include <UpdateMomentumHigueraCary.H>
using namespace amrex;
@@ -41,9 +43,26 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
pp.query("do_continuous_injection", do_continuous_injection);
// Whether to plot back-transformed (lab-frame) diagnostics
// for this species.
- pp.query("do_boosted_frame_diags", do_boosted_frame_diags);
+ pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics);
pp.query("do_field_ionization", do_field_ionization);
+
+ //check if Radiation Reaction is enabled and do consistency checks
+ pp.query("do_classical_radiation_reaction", do_classical_radiation_reaction);
+ //if the species is not a lepton, do_classical_radiation_reaction
+ //should be false
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+ !(do_classical_radiation_reaction && !AmIALepton()),
+ "Can't enable classical radiation reaction for non lepton species. " );
+
+ //Only Boris pusher is compatible with radiation reaction
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+ !(do_classical_radiation_reaction &&
+ WarpX::particle_pusher_algo != ParticlePusherAlgo::Boris),
+ "Radiation reaction can be enabled only if Boris pusher is used");
+ //_____________________________
+
+
#ifdef AMREX_USE_GPU
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
do_field_ionization == 0,
@@ -67,7 +86,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
//variable to set plot_flags size
int plot_flag_size = PIdx::nattribs;
- if(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
+ if(WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics)
plot_flag_size += 6;
#ifdef WARPX_QED
@@ -422,13 +441,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) {
Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] );
- overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, Real(0.)) * dx[dir]);
+ overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]);
} else {
no_overlap = true; break;
}
if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) {
Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] );
- overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, Real(0.)) * dx[dir]);
+ overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]);
} else {
no_overlap = true; break;
}
@@ -1059,7 +1078,7 @@ PhysicalParticleContainer::Evolve (int lev,
bool has_buffer = cEx || cjx;
- if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
+ if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics)
{
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
@@ -1408,9 +1427,9 @@ PhysicalParticleContainer::SplitParticles(int lev)
// before splitting results in a uniform distribution after splitting
const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim;
const std::array<Real,3>& dx = WarpX::CellSize(lev);
- amrex::Vector<amrex::Real> split_offset = {dx[0]/2./ppc_nd[0],
- dx[1]/2./ppc_nd[1],
- dx[2]/2./ppc_nd[2]};
+ amrex::Vector<Real> split_offset = {dx[0]/2._rt/ppc_nd[0],
+ dx[1]/2._rt/ppc_nd[1],
+ dx[2]/2._rt/ppc_nd[2]};
// particle Array Of Structs data
auto& particles = pti.GetArrayOfStructs();
@@ -1566,7 +1585,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
- if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf))
+ if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf))
{
copy_attribs(pti, x, y, z);
}
@@ -1579,7 +1598,23 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
// Loop over the particles and update their momentum
const Real q = this->charge;
const Real m = this-> mass;
- if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
+
+
+ //Assumes that all consistency checks have been done at initialization
+ if(do_classical_radiation_reaction){
+ amrex::ParallelFor(
+ pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i],
+ Ex[i], Ey[i], Ez[i], Bx[i],
+ By[i], Bz[i], qp, m, dt);
+ UpdatePosition( x[i], y[i], z[i],
+ ux[i], uy[i], uz[i], dt );
+ }
+ );
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
amrex::ParallelFor(
pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
@@ -1605,6 +1640,19 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
ux[i], uy[i], uz[i], dt );
}
);
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
+ amrex::ParallelFor(
+ pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumHigueraCary( ux[i], uy[i], uz[i],
+ Ex[i], Ey[i], Ez[i], Bx[i],
+ By[i], Bz[i], qp, m, dt);
+ UpdatePosition( x[i], y[i], z[i],
+ ux[i], uy[i], uz[i], dt );
+ }
+ );
} else {
amrex::Abort("Unknown particle pusher");
};
@@ -1687,17 +1735,55 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
// Loop over the particles and update their momentum
const Real q = this->charge;
const Real m = this-> mass;
- if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
+
+ int* AMREX_RESTRICT ion_lev = nullptr;
+ if (do_field_ionization){
+ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
+ }
+
+
+ //Assumes that all consistency checks have been done at initialization
+ if(do_classical_radiation_reaction){
+ amrex::ParallelFor(
+ pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumBorisWithRadiationReaction(
+ ux[i], uy[i], uz[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ qp, m, dt);
+ }
+ );
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
amrex::ParallelFor( pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
- UpdateMomentumBoris( ux[i], uy[i], uz[i],
- Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumBoris(
+ ux[i], uy[i], uz[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ qp, m, dt);
}
);
} else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) {
amrex::ParallelFor( pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
- UpdateMomentumVay( ux[i], uy[i], uz[i],
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumVay(
+ ux[i], uy[i], uz[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ qp, m, dt);
+ }
+ );
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
+ amrex::ParallelFor( pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ UpdateMomentumHigueraCary( ux[i], uy[i], uz[i],
Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
}
);
@@ -1756,7 +1842,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
// Note the the slice should always move in the negative boost direction.
AMREX_ALWAYS_ASSERT(z_new < z_old);
- AMREX_ALWAYS_ASSERT(do_boosted_frame_diags == 1);
+ AMREX_ALWAYS_ASSERT(do_back_transformed_diagnostics == 1);
const int nlevs = std::max(0, finestLevel()+1);
@@ -2142,6 +2228,13 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const
);
}
+//This function return true if the PhysicalParticleContainer contains electrons
+//or positrons, false otherwise
+bool
+PhysicalParticleContainer::AmIALepton(){
+ return (this-> mass == PhysConst::m_e);
+}
+
#ifdef WARPX_QED
bool PhysicalParticleContainer::has_quantum_sync()
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-25 18:34:12 +0000