From 31b762139179b9c8dc74429c73f979636d56dc29 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 25 Oct 2019 14:28:34 +0200 Subject: Separate PushPX kernels in a dedicated function --- Source/Particles/PhysicalParticleContainer.cpp | 31 +++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d12a4dbff..723a35ac1 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1599,11 +1599,32 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real q = this->charge; const Real m = this-> mass; + PushPX_classical(x, y, z, ux, uy, uz, + Ex, Ey, Ez, Bx, By, Bz, + q, m, ion_lev, dt, pti.numParticles()); +} +void PhysicalParticleContainer::PushPX_classical( + ParticleReal* const AMREX_RESTRICT x, + ParticleReal* const AMREX_RESTRICT y, + ParticleReal* const AMREX_RESTRICT z, + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Ex, + const ParticleReal* const AMREX_RESTRICT Ey, + const ParticleReal* const AMREX_RESTRICT Ez, + const ParticleReal* const AMREX_RESTRICT Bx, + const ParticleReal* const AMREX_RESTRICT By, + const ParticleReal* const AMREX_RESTRICT Bz, + Real q, Real m, int* AMREX_RESTRICT ion_lev, + amrex::Real dt, long num_particles +) +{ //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1616,7 +1637,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1629,7 +1650,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1642,7 +1663,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { amrex::ParallelFor( - pti.numParticles(), + num_particles, [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1655,7 +1676,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ); } else { amrex::Abort("Unknown particle pusher"); - }; + } } void -- cgit v1.2.3 From 3a37ca0c138080a60213b7a3353769079d46589c Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 25 Oct 2019 15:18:16 +0200 Subject: added some of the kernels to evolve optical depth --- Source/Particles/PhysicalParticleContainer.H | 22 +++++ Source/Particles/PhysicalParticleContainer.cpp | 132 +++++++++++++++++++++++++ 2 files changed, 154 insertions(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 8f3a21b38..098acbcd5 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -252,6 +252,28 @@ private: ); //________________________ +#ifdef WARPX_QED + //PushPX + QED (Quantum Synchrotron) + void PushPX_QedQuantumSynchrotron( + amrex::ParticleReal* const AMREX_RESTRICT x, + amrex::ParticleReal* const AMREX_RESTRICT y, + amrex::ParticleReal* const AMREX_RESTRICT z, + amrex::ParticleReal* const AMREX_RESTRICT ux, + amrex::ParticleReal* const AMREX_RESTRICT uy, + amrex::ParticleReal* const AMREX_RESTRICT uz, + const amrex::ParticleReal* const AMREX_RESTRICT Ex, + const amrex::ParticleReal* const AMREX_RESTRICT Ey, + const amrex::ParticleReal* const AMREX_RESTRICT Ez, + const amrex::ParticleReal* const AMREX_RESTRICT Bx, + const amrex::ParticleReal* const AMREX_RESTRICT By, + const amrex::ParticleReal* const AMREX_RESTRICT Bz, + amrex::Real q, amrex::Real m, + amrex::ParticleReal* const AMREX_RESTRICT tau, + amrex::Real dt, long num_particles + ); + //________________________ +#endif + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 723a35ac1..d3212baf3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1599,9 +1599,25 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real q = this->charge; const Real m = this-> mass; +#ifdef WARPX_QED + if(do_qed_quantum_sync){ + amrex::Real* AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); + + PushPX_QedQuantumSynchrotron(x, y, z, ux, uy, uz, + Ex, Ey, Ez, Bx, By, Bz, + q, m, p_tau, dt, pti.numParticles()); + } + else { + PushPX_classical(x, y, z, ux, uy, uz, + Ex, Ey, Ez, Bx, By, Bz, + q, m, ion_lev, dt, pti.numParticles()); + } +#else PushPX_classical(x, y, z, ux, uy, uz, Ex, Ey, Ez, Bx, By, Bz, q, m, ion_lev, dt, pti.numParticles()); +#endif } void PhysicalParticleContainer::PushPX_classical( @@ -1679,6 +1695,122 @@ void PhysicalParticleContainer::PushPX_classical( } } +#ifdef WARPX_QED +void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( + amrex::ParticleReal* const AMREX_RESTRICT x, + amrex::ParticleReal* const AMREX_RESTRICT y, + amrex::ParticleReal* const AMREX_RESTRICT z, + amrex::ParticleReal* const AMREX_RESTRICT ux, + amrex::ParticleReal* const AMREX_RESTRICT uy, + amrex::ParticleReal* const AMREX_RESTRICT uz, + const amrex::ParticleReal* const AMREX_RESTRICT Ex, + const amrex::ParticleReal* const AMREX_RESTRICT Ey, + const amrex::ParticleReal* const AMREX_RESTRICT Ez, + const amrex::ParticleReal* const AMREX_RESTRICT Bx, + const amrex::ParticleReal* const AMREX_RESTRICT By, + const amrex::ParticleReal* const AMREX_RESTRICT Bz, + amrex::Real q, amrex::Real m, + amrex::ParticleReal* const AMREX_RESTRICT tau, + amrex::Real dt, long num_particles +) +{ + QuantumSynchrotronEvolveOpticalDepth evolve_opt = + shr_ptr_qs_engine->build_evolve_functor(); + + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::Abort("QED + Classical Radiation Reaction has still to be implemented!"); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + amrex::ParallelFor( + num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const amrex::ParticleReal ux_old = ux[i]; + const amrex::ParticleReal uy_old = uy[i]; + const amrex::ParticleReal uz_old = uz[i]; + + UpdateMomentumBoris( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + + const amrex::ParticleReal half_mass = 0.5*m; + + const amrex::ParticleReal px_n = half_mass * (ux[i]+ux_old); + const amrex::ParticleReal py_n = half_mass * (uy[i]+uy_old); + const amrex::ParticleReal pz_n = half_mass * (uz[i]+uz_old); + + bool has_event_happened = evolve_opt( + px_n, py_n, pz_n, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, tau[i]); + + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + amrex::ParallelFor( + num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const amrex::ParticleReal ux_old = ux[i]; + const amrex::ParticleReal uy_old = uy[i]; + const amrex::ParticleReal uz_old = uz[i]; + + UpdateMomentumVay( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + + const amrex::ParticleReal half_mass = 0.5*m; + + const amrex::ParticleReal px_n = half_mass * (ux[i]+ux_old); + const amrex::ParticleReal py_n = half_mass * (uy[i]+uy_old); + const amrex::ParticleReal pz_n = half_mass * (uz[i]+uz_old); + + bool has_event_happened = evolve_opt( + px_n, py_n, pz_n, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, tau[i]); + + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { + amrex::ParallelFor( + num_particles, + [=] AMREX_GPU_DEVICE (long i) { + + const amrex::ParticleReal ux_old = ux[i]; + const amrex::ParticleReal uy_old = uy[i]; + const amrex::ParticleReal uz_old = uz[i]; + + UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + + const amrex::ParticleReal half_mass = 0.5*m; + + const amrex::ParticleReal px_n = half_mass * (ux[i]+ux_old); + const amrex::ParticleReal py_n = half_mass * (uy[i]+uy_old); + const amrex::ParticleReal pz_n = half_mass * (uz[i]+uz_old); + + bool has_event_happened = evolve_opt( + px_n, py_n, pz_n, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, tau[i]); + + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + } else { + amrex::Abort("Unknown particle pusher"); + } +} +#endif + void PhysicalParticleContainer::PushP (int lev, Real dt, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, -- cgit v1.2.3 From 905295431bbcd2a8824ffe97e340bb1910589062 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 29 Oct 2019 11:37:00 +0100 Subject: adding tau --- Source/Particles/PhysicalParticleContainer.H | 30 +++ Source/Particles/PhysicalParticleContainer.cpp | 277 +++++++++++++++++-------- 2 files changed, 221 insertions(+), 86 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index c72312fee..c82d14ebf 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -251,6 +251,21 @@ public: amrex::Real q, amrex::Real m, int* AMREX_RESTRICT ion_lev, amrex::Real dt, long num_particles ); + + void PushP_classical( + amrex::ParticleReal* const AMREX_RESTRICT ux, + amrex::ParticleReal* const AMREX_RESTRICT uy, + amrex::ParticleReal* const AMREX_RESTRICT uz, + const amrex::ParticleReal* const AMREX_RESTRICT Expp, + const amrex::ParticleReal* const AMREX_RESTRICT Eypp, + const amrex::ParticleReal* const AMREX_RESTRICT Ezpp, + const amrex::ParticleReal* const AMREX_RESTRICT Bxpp, + const amrex::ParticleReal* const AMREX_RESTRICT Bypp, + const amrex::ParticleReal* const AMREX_RESTRICT Bzpp, + amrex::Real q, amrex::Real m, int* AMREX_RESTRICT ion_lev, + amrex::Real dt, long num_particles + ); + //________________________ #ifdef WARPX_QED @@ -272,6 +287,21 @@ public: amrex::ParticleReal* const AMREX_RESTRICT tau, amrex::Real dt, long num_particles ); + + void PushP_QedQuantumSynchrotron( + amrex::ParticleReal* const AMREX_RESTRICT ux, + amrex::ParticleReal* const AMREX_RESTRICT uy, + amrex::ParticleReal* const AMREX_RESTRICT uz, + const amrex::ParticleReal* const AMREX_RESTRICT Expp, + const amrex::ParticleReal* const AMREX_RESTRICT Eypp, + const amrex::ParticleReal* const AMREX_RESTRICT Ezpp, + const amrex::ParticleReal* const AMREX_RESTRICT Bxpp, + const amrex::ParticleReal* const AMREX_RESTRICT Bypp, + const amrex::ParticleReal* const AMREX_RESTRICT Bzpp, + amrex::Real q, amrex::Real m, + amrex::ParticleReal* const AMREX_RESTRICT tau, + amrex::Real dt, long num_particles + ); //________________________ #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7e6293a41..bc8c8c3c4 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1600,8 +1600,8 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real m = this-> mass; #ifdef WARPX_QED - if(do_qed_quantum_sync){ - amrex::Real* AMREX_RESTRICT p_tau = + if(has_quantum_sync()){ + Real* AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); PushPX_QedQuantumSynchrotron(x, y, z, ux, uy, uz, @@ -1697,21 +1697,21 @@ void PhysicalParticleContainer::PushPX_classical( #ifdef WARPX_QED void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( - amrex::ParticleReal* const AMREX_RESTRICT x, - amrex::ParticleReal* const AMREX_RESTRICT y, - amrex::ParticleReal* const AMREX_RESTRICT z, - amrex::ParticleReal* const AMREX_RESTRICT ux, - amrex::ParticleReal* const AMREX_RESTRICT uy, - amrex::ParticleReal* const AMREX_RESTRICT uz, - const amrex::ParticleReal* const AMREX_RESTRICT Ex, - const amrex::ParticleReal* const AMREX_RESTRICT Ey, - const amrex::ParticleReal* const AMREX_RESTRICT Ez, - const amrex::ParticleReal* const AMREX_RESTRICT Bx, - const amrex::ParticleReal* const AMREX_RESTRICT By, - const amrex::ParticleReal* const AMREX_RESTRICT Bz, - amrex::Real q, amrex::Real m, - amrex::ParticleReal* const AMREX_RESTRICT tau, - amrex::Real dt, long num_particles + ParticleReal* const AMREX_RESTRICT x, + ParticleReal* const AMREX_RESTRICT y, + ParticleReal* const AMREX_RESTRICT z, + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Ex, + const ParticleReal* const AMREX_RESTRICT Ey, + const ParticleReal* const AMREX_RESTRICT Ez, + const ParticleReal* const AMREX_RESTRICT Bx, + const ParticleReal* const AMREX_RESTRICT By, + const ParticleReal* const AMREX_RESTRICT Bz, + Real q, amrex::Real m, + ParticleReal* const AMREX_RESTRICT tau, + Real dt, long num_particles ) { QuantumSynchrotronEvolveOpticalDepth evolve_opt = @@ -1724,19 +1724,19 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( amrex::ParallelFor( num_particles, [=] AMREX_GPU_DEVICE (long i) { - const amrex::ParticleReal ux_old = ux[i]; - const amrex::ParticleReal uy_old = uy[i]; - const amrex::ParticleReal uz_old = uz[i]; + const ParticleReal ux_old = ux[i]; + const ParticleReal uy_old = uy[i]; + const ParticleReal uz_old = uz[i]; UpdateMomentumBoris( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - const amrex::ParticleReal half_mass = 0.5*m; + const ParticleReal half_mass = 0.5*m; - const amrex::ParticleReal px_n = half_mass * (ux[i]+ux_old); - const amrex::ParticleReal py_n = half_mass * (uy[i]+uy_old); - const amrex::ParticleReal pz_n = half_mass * (uz[i]+uz_old); + const ParticleReal px_n = half_mass * (ux[i]+ux_old); + const ParticleReal py_n = half_mass * (uy[i]+uy_old); + const ParticleReal pz_n = half_mass * (uz[i]+uz_old); bool has_event_happened = evolve_opt( px_n, py_n, pz_n, @@ -1752,19 +1752,19 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( amrex::ParallelFor( num_particles, [=] AMREX_GPU_DEVICE (long i) { - const amrex::ParticleReal ux_old = ux[i]; - const amrex::ParticleReal uy_old = uy[i]; - const amrex::ParticleReal uz_old = uz[i]; + const ParticleReal ux_old = ux[i]; + const ParticleReal uy_old = uy[i]; + const ParticleReal uz_old = uz[i]; UpdateMomentumVay( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - const amrex::ParticleReal half_mass = 0.5*m; + const ParticleReal half_mass = 0.5*m; - const amrex::ParticleReal px_n = half_mass * (ux[i]+ux_old); - const amrex::ParticleReal py_n = half_mass * (uy[i]+uy_old); - const amrex::ParticleReal pz_n = half_mass * (uz[i]+uz_old); + const ParticleReal px_n = half_mass * (ux[i]+ux_old); + const ParticleReal py_n = half_mass * (uy[i]+uy_old); + const ParticleReal pz_n = half_mass * (uz[i]+uz_old); bool has_event_happened = evolve_opt( px_n, py_n, pz_n, @@ -1781,19 +1781,19 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( num_particles, [=] AMREX_GPU_DEVICE (long i) { - const amrex::ParticleReal ux_old = ux[i]; - const amrex::ParticleReal uy_old = uy[i]; - const amrex::ParticleReal uz_old = uz[i]; + const ParticleReal ux_old = ux[i]; + const ParticleReal uy_old = uy[i]; + const ParticleReal uz_old = uz[i]; UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - const amrex::ParticleReal half_mass = 0.5*m; + const ParticleReal half_mass = 0.5*m; - const amrex::ParticleReal px_n = half_mass * (ux[i]+ux_old); - const amrex::ParticleReal py_n = half_mass * (uy[i]+uy_old); - const amrex::ParticleReal pz_n = half_mass * (uz[i]+uz_old); + const ParticleReal px_n = half_mass * (ux[i]+ux_old); + const ParticleReal py_n = half_mass * (uy[i]+uy_old); + const ParticleReal pz_n = half_mass * (uz[i]+uz_old); bool has_event_happened = evolve_opt( px_n, py_n, pz_n, @@ -1894,59 +1894,164 @@ PhysicalParticleContainer::PushP (int lev, Real dt, ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } +#ifndef WARPX_QED + PushP_classical( + ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, + q, m, ion_lev, dt, pti.numParticles()); +#else + if(has_quantum_sync()){ + Real* AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); + PushP_QedQuantumSynchrotron( + ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, + q, m, p_tau, dt, pti.numParticles()); + } + else{ + PushP_classical( + ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, + q, m, ion_lev, dt, pti.numParticles()); + } +#endif + - //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) { - 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) { - 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); - } - ); - } else { - amrex::Abort("Unknown particle pusher"); - }; } } } +void PhysicalParticleContainer::PushP_classical( + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Expp, + const ParticleReal* const AMREX_RESTRICT Eypp, + const ParticleReal* const AMREX_RESTRICT Ezpp, + const ParticleReal* const AMREX_RESTRICT Bxpp, + const ParticleReal* const AMREX_RESTRICT Bypp, + const ParticleReal* const AMREX_RESTRICT Bzpp, + Real q, Real m, int* AMREX_RESTRICT ion_lev, + Real dt, long num_particles +) +{ + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::ParallelFor(num_particles, + [=] 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(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + 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(num_particles, + [=] AMREX_GPU_DEVICE (long 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(num_particles, + [=] 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); + } + ); + } else { + amrex::Abort("Unknown particle pusher"); + } +} + +void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( + ParticleReal* const AMREX_RESTRICT ux, + ParticleReal* const AMREX_RESTRICT uy, + ParticleReal* const AMREX_RESTRICT uz, + const ParticleReal* const AMREX_RESTRICT Expp, + const ParticleReal* const AMREX_RESTRICT Eypp, + const ParticleReal* const AMREX_RESTRICT Ezpp, + const ParticleReal* const AMREX_RESTRICT Bxpp, + const ParticleReal* const AMREX_RESTRICT Bypp, + const ParticleReal* const AMREX_RESTRICT Bzpp, + Real q, Real m, + ParticleReal* const AMREX_RESTRICT tau, + Real dt, long num_particles +) +{ + QuantumSynchrotronEvolveOpticalDepth evolve_opt = + shr_ptr_qs_engine->build_evolve_functor(); + + + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::Abort("QED + Classical Radiation Reaction has still to be implemented!"); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal ux_old = ux[i]; + const ParticleReal uy_old = uy[i]; + const ParticleReal uz_old = uz[i]; + + UpdateMomentumBoris( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + + const ParticleReal half_mass = 0.5*m; + + const ParticleReal px_n = half_mass * (ux[i]+ux_old); + const ParticleReal py_n = half_mass * (uy[i]+uy_old); + const ParticleReal pz_n = half_mass * (uz[i]+uz_old); + + bool has_event_happened = evolve_opt( + px_n, py_n, pz_n, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); + + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumVay( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ + amrex::ParallelFor(num_particles, + [=] 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); + } + ); + } else { + amrex::Abort("Unknown particle pusher"); + } +} + void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, const ParticleReal* yp, const ParticleReal* zp) { -- cgit v1.2.3 From c98c94d8bd87d653bd76fb4e6fda47f70d596c8f Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 29 Oct 2019 12:09:18 +0100 Subject: simplification --- Source/Particles/PhotonParticleContainer.cpp | 15 +++- Source/Particles/PhysicalParticleContainer.cpp | 110 ++++++++++++------------- 2 files changed, 64 insertions(+), 61 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 9197be6e3..30d010e6b 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -224,8 +224,13 @@ PhotonParticleContainer::DoBreitWheeler(int lev, amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + + const ParticleReal px = me * ux[i]; + const ParticleReal py = me * uy[i]; + const ParticleReal pz = me * uz[i]; + evolve_opt( - me*ux[i], me*uy[i], me*uz[i], + px, py, pz, Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], dt, p_tau[i]); @@ -261,8 +266,12 @@ PhotonParticleContainer::DoBreitWheelerPti(WarpXParIter& pti, amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - evolve_opt( - me*ux[i], me*uy[i], me*uz[i], + const ParticleReal px = me * ux[i]; + const ParticleReal py = me * uy[i]; + const ParticleReal pz = me * uz[i]; + + bool has_event_happened = evolve_opt( + px, py, pz, Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], dt, p_tau[i]); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index bc8c8c3c4..410476da5 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1724,26 +1724,20 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( amrex::ParallelFor( num_particles, [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal ux_old = ux[i]; - const ParticleReal uy_old = uy[i]; - const ParticleReal uz_old = uz[i]; - - UpdateMomentumBoris( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); - - const ParticleReal half_mass = 0.5*m; - - const ParticleReal px_n = half_mass * (ux[i]+ux_old); - const ParticleReal py_n = half_mass * (uy[i]+uy_old); - const ParticleReal pz_n = half_mass * (uz[i]+uz_old); + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; bool has_event_happened = evolve_opt( - px_n, py_n, pz_n, + px, py, pz, Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], dt, tau[i]); + UpdateMomentumBoris( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + UpdatePosition( x[i], y[i], z[i], ux[i], uy[i], uz[i], dt ); } @@ -1752,26 +1746,20 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( amrex::ParallelFor( num_particles, [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal ux_old = ux[i]; - const ParticleReal uy_old = uy[i]; - const ParticleReal uz_old = uz[i]; - - UpdateMomentumVay( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); - - const ParticleReal half_mass = 0.5*m; - - const ParticleReal px_n = half_mass * (ux[i]+ux_old); - const ParticleReal py_n = half_mass * (uy[i]+uy_old); - const ParticleReal pz_n = half_mass * (uz[i]+uz_old); + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; bool has_event_happened = evolve_opt( - px_n, py_n, pz_n, + px, py, pz, Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], dt, tau[i]); + UpdateMomentumVay( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + UpdatePosition( x[i], y[i], z[i], ux[i], uy[i], uz[i], dt ); } @@ -1780,27 +1768,20 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( amrex::ParallelFor( num_particles, [=] AMREX_GPU_DEVICE (long i) { - - const ParticleReal ux_old = ux[i]; - const ParticleReal uy_old = uy[i]; - const ParticleReal uz_old = uz[i]; - - UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); - - const ParticleReal half_mass = 0.5*m; - - const ParticleReal px_n = half_mass * (ux[i]+ux_old); - const ParticleReal py_n = half_mass * (uy[i]+uy_old); - const ParticleReal pz_n = half_mass * (uz[i]+uz_old); + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; bool has_event_happened = evolve_opt( - px_n, py_n, pz_n, + px, py, pz, Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], dt, tau[i]); + UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + UpdatePosition( x[i], y[i], z[i], ux[i], uy[i], uz[i], dt ); } @@ -2006,33 +1987,36 @@ void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal ux_old = ux[i]; - const ParticleReal uy_old = uy[i]; - const ParticleReal uz_old = uz[i]; + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; + + bool has_event_happened = evolve_opt( + px, py, pz, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); UpdateMomentumBoris( ux[i], uy[i], uz[i], Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); - - const ParticleReal half_mass = 0.5*m; - - const ParticleReal px_n = half_mass * (ux[i]+ux_old); - const ParticleReal py_n = half_mass * (uy[i]+uy_old); - const ParticleReal pz_n = half_mass * (uz[i]+uz_old); + } + ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; bool has_event_happened = evolve_opt( - px_n, py_n, pz_n, + px, py, pz, Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], dt, tau[i]); - } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ - amrex::ParallelFor(num_particles, - [=] AMREX_GPU_DEVICE (long i) { UpdateMomentumVay( ux[i], uy[i], uz[i], Expp[i], Eypp[i], Ezpp[i], @@ -2043,6 +2027,16 @@ void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; + + bool has_event_happened = evolve_opt( + px, py, pz, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); + UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } -- cgit v1.2.3 From d951dafd9ed46dc6a0ee9a2821773b1e53bef970 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 29 Oct 2019 15:08:51 +0100 Subject: Added QS to PushPX --- Source/Particles/PhysicalParticleContainer.cpp | 110 +++++++++++++++++++++---- Source/QED/BreitWheelerEngineWrapper.H | 6 ++ Source/QED/BreitWheelerEngineWrapper.cpp | 6 ++ Source/QED/QuantumSyncEngineWrapper.H | 6 ++ Source/QED/QuantumSyncEngineWrapper.cpp | 6 ++ 5 files changed, 118 insertions(+), 16 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 410476da5..40ab1d0ca 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1717,12 +1717,54 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( QuantumSynchrotronEvolveOpticalDepth evolve_opt = shr_ptr_qs_engine->build_evolve_functor(); + const auto chi_min = shr_ptr_qs_engine->get_ref_ctrl().chi_part_min; + //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ - amrex::Abort("QED + Classical Radiation Reaction has still to be implemented!"); + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; + + const ParticleReal chi = QedUtils::chi_lepton(px, py, pz, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i]); + + std::cout << chi << " " << chi_min << " "; + + if(chi < chi_min){ + UpdateMomentumBorisWithRadiationReaction( + ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + q, m, dt); + + std::cout << "CRR \n"; + } + else + { + bool has_event_happened = evolve_opt( + px, py, pz, + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + dt, tau[i]); + + UpdateMomentumBoris( + ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i], + q, m, dt); + + std::cout << "QRR \n"; + } + + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( - num_particles, + amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (long i) { const ParticleReal px = m * ux[i]; const ParticleReal py = m * uy[i]; @@ -1743,8 +1785,7 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { - amrex::ParallelFor( - num_particles, + amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (long i) { const ParticleReal px = m * ux[i]; const ParticleReal py = m * uy[i]; @@ -1875,11 +1916,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } -#ifndef WARPX_QED - PushP_classical( - ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, - q, m, ion_lev, dt, pti.numParticles()); -#else +#ifdef WARPX_QED if(has_quantum_sync()){ Real* AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); @@ -1892,6 +1929,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt, ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, q, m, ion_lev, dt, pti.numParticles()); } +#else + PushP_classical( + ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, + q, m, ion_lev, dt, pti.numParticles()); #endif @@ -1962,6 +2003,8 @@ void PhysicalParticleContainer::PushP_classical( } } +#ifdef WARPX_QED + void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( ParticleReal* const AMREX_RESTRICT ux, ParticleReal* const AMREX_RESTRICT uy, @@ -1980,10 +2023,43 @@ void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( QuantumSynchrotronEvolveOpticalDepth evolve_opt = shr_ptr_qs_engine->build_evolve_functor(); + const auto chi_min = shr_ptr_qs_engine->get_ref_ctrl().chi_part_min; //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ - amrex::Abort("QED + Classical Radiation Reaction has still to be implemented!"); + amrex::ParallelFor(num_particles, + [=] AMREX_GPU_DEVICE (long i) { + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; + + const ParticleReal chi = QedUtils::chi_lepton(px, py, pz, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i]); + + if(chi < chi_min){ + UpdateMomentumBorisWithRadiationReaction( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + else + { + bool has_event_happened = evolve_opt( + px, py, pz, + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + dt, tau[i]); + + UpdateMomentumBoris( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + } + ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (long i) { @@ -2002,8 +2078,8 @@ void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); - } - ); + } + ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (long i) { @@ -2022,8 +2098,8 @@ void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); - } - ); + } + ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ amrex::ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (long i) { @@ -2039,13 +2115,15 @@ void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); - } + } ); } else { amrex::Abort("Unknown particle pusher"); } } +#endif + void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, const ParticleReal* yp, const ParticleReal* zp) { diff --git a/Source/QED/BreitWheelerEngineWrapper.H b/Source/QED/BreitWheelerEngineWrapper.H index e3324c3f3..7cb60a3c0 100644 --- a/Source/QED/BreitWheelerEngineWrapper.H +++ b/Source/QED/BreitWheelerEngineWrapper.H @@ -270,6 +270,12 @@ public: */ PicsarBreitWheelerCtrl get_default_ctrl() const; + /** + * returns a constant reference to the control parameters + * @return const reference to control parameters + */ + const PicsarBreitWheelerCtrl& get_ref_ctrl() const; + private: bool lookup_tables_initialized = false; diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp index ec2c65232..78ad0f9de 100644 --- a/Source/QED/BreitWheelerEngineWrapper.cpp +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -200,6 +200,12 @@ BreitWheelerEngine::get_default_ctrl() const return PicsarBreitWheelerCtrl(); } +const PicsarBreitWheelerCtrl& +BreitWheelerEngine::get_ref_ctrl() const +{ + return innards.ctrl; +} + void BreitWheelerEngine::compute_lookup_tables ( PicsarBreitWheelerCtrl ctrl) { diff --git a/Source/QED/QuantumSyncEngineWrapper.H b/Source/QED/QuantumSyncEngineWrapper.H index df21a1ea0..8750cff5e 100644 --- a/Source/QED/QuantumSyncEngineWrapper.H +++ b/Source/QED/QuantumSyncEngineWrapper.H @@ -270,6 +270,12 @@ public: */ PicsarQuantumSynchrotronCtrl get_default_ctrl() const; + /** + * returns a constant reference to the control parameters + * @return const reference to control parameters + */ + const PicsarQuantumSynchrotronCtrl& get_ref_ctrl() const; + private: bool lookup_tables_initialized = false; diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp index 50a3b675c..0abe7010d 100644 --- a/Source/QED/QuantumSyncEngineWrapper.cpp +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -200,6 +200,12 @@ QuantumSynchrotronEngine::get_default_ctrl() const return PicsarQuantumSynchrotronCtrl(); } +const PicsarQuantumSynchrotronCtrl& +QuantumSynchrotronEngine::get_ref_ctrl() const +{ + return innards.ctrl; +} + void QuantumSynchrotronEngine::compute_lookup_tables ( PicsarQuantumSynchrotronCtrl ctrl) { -- cgit v1.2.3 From 7c1a7f9fd8239fa1eeccb9fe5cd2b34e2cf589c3 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 29 Oct 2019 15:12:58 +0100 Subject: removed unwanted debug messages --- Source/Particles/PhysicalParticleContainer.cpp | 6 ------ 1 file changed, 6 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 40ab1d0ca..b7b2d5df4 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1731,16 +1731,12 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i]); - std::cout << chi << " " << chi_min << " "; - if(chi < chi_min){ UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - - std::cout << "CRR \n"; } else { @@ -1755,8 +1751,6 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); - - std::cout << "QRR \n"; } UpdatePosition( x[i], y[i], z[i], -- cgit v1.2.3 From 2b944018f27058262617358da04694d0e0f6808a Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 30 Oct 2019 16:34:49 +0100 Subject: fixed some bugs due to merging --- Source/Particles/PhotonParticleContainer.cpp | 4 ++-- Source/Particles/PhysicalParticleContainer.cpp | 8 ++++---- Source/QED/BreitWheelerEngineWrapper.cpp | 2 +- Source/QED/QuantumSyncEngineWrapper.cpp | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index eb0dfb373..354b9587c 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -213,7 +213,7 @@ PhotonParticleContainer::DoBreitWheeler(int lev, const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); BreitWheelerEvolveOpticalDepth evolve_opt = - shr_ptr_bw_engine->build_evolve_functor(); + m_shr_p_bw_engine->build_evolve_functor(); amrex::Real* AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); @@ -256,7 +256,7 @@ PhotonParticleContainer::DoBreitWheelerPti(WarpXParIter& pti, const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); BreitWheelerEvolveOpticalDepth evolve_opt = - shr_ptr_bw_engine->build_evolve_functor(); + m_shr_p_bw_engine->build_evolve_functor(); amrex::Real* AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a02ac8593..ee0711372 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1716,9 +1716,9 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( ) { QuantumSynchrotronEvolveOpticalDepth evolve_opt = - shr_ptr_qs_engine->build_evolve_functor(); + m_shr_p_qs_engine->build_evolve_functor(); - const auto chi_min = shr_ptr_qs_engine->get_ref_ctrl().chi_part_min; + const auto chi_min = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ @@ -2016,9 +2016,9 @@ void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( ) { QuantumSynchrotronEvolveOpticalDepth evolve_opt = - shr_ptr_qs_engine->build_evolve_functor(); + m_shr_p_qs_engine->build_evolve_functor(); - const auto chi_min = shr_ptr_qs_engine->get_ref_ctrl().chi_part_min; + const auto chi_min = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp index f9ab6e0aa..78ff13fc5 100644 --- a/Source/QED/BreitWheelerEngineWrapper.cpp +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -179,7 +179,7 @@ BreitWheelerEngine::get_default_ctrl() const const PicsarBreitWheelerCtrl& BreitWheelerEngine::get_ref_ctrl() const { - return innards.ctrl; + return m_innards.ctrl; } void BreitWheelerEngine::compute_lookup_tables ( diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp index 4c1389ce5..c2d9f0abf 100644 --- a/Source/QED/QuantumSyncEngineWrapper.cpp +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -178,7 +178,7 @@ QuantumSynchrotronEngine::get_default_ctrl() const const PicsarQuantumSynchrotronCtrl& QuantumSynchrotronEngine::get_ref_ctrl() const { - return innards.ctrl; + return m_innards.ctrl; } void QuantumSynchrotronEngine::compute_lookup_tables ( -- cgit v1.2.3 From ecaedd97cb868a266971c81545c68930a2802126 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 31 Oct 2019 11:19:20 +0100 Subject: removed QED from PushP --- Source/Particles/PhotonParticleContainer.H | 18 +- Source/Particles/PhotonParticleContainer.cpp | 124 +----------- Source/Particles/PhysicalParticleContainer.H | 57 +----- Source/Particles/PhysicalParticleContainer.cpp | 250 +++++-------------------- 4 files changed, 54 insertions(+), 395 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index b1cc283c3..fbbc1077e 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -46,6 +46,7 @@ public: amrex::Cuda::ManagedDeviceVector& zp, amrex::Real dt, DtType a_dt_type=DtType::Full) override; + // Do nothing virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, @@ -53,7 +54,7 @@ public: const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, - const amrex::MultiFab& Bz) override; + const amrex::MultiFab& Bz) override {}; // DepositCurrent should do nothing for photons @@ -79,24 +80,11 @@ public: }; #ifdef WARPX_QED - /** - * This functions (called by PushP) evolves optical depths - * of the photons according to the Breit Wheeler process - */ - void DoBreitWheeler(int lev, - amrex::Real dt, - const amrex::MultiFab& Ex, - const amrex::MultiFab& Ey, - const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, - const amrex::MultiFab& By, - const amrex::MultiFab& Bz); - /** * This functions (called by PushPX) evolves optical depths * of the photons according to the Breit Wheeler process */ - void DoBreitWheelerPti(WarpXParIter& pti, + void DoBreitWheeler(WarpXParIter& pti, amrex::Real dt); #endif diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 354b9587c..afe697b7e 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -88,29 +88,11 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, ); #ifdef WARPX_QED - if(has_breit_wheeler()) DoBreitWheelerPti(pti, dt); + if(has_breit_wheeler()) DoBreitWheeler(pti, dt); #endif } -void -PhotonParticleContainer::PushP (int lev, - amrex::Real dt, - const amrex::MultiFab& Ex, - const amrex::MultiFab& Ey, - const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, - const amrex::MultiFab& By, - const amrex::MultiFab& Bz) -{ - BL_PROFILE("PhotonParticleContainer::PushP"); - if (do_not_push) return; - -#ifdef WARPX_QED - if(has_breit_wheeler()) DoBreitWheeler(lev,dt, Ex,Ey,Ez,Bx,By,Bz); -#endif -} - void PhotonParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, @@ -139,109 +121,7 @@ PhotonParticleContainer::Evolve (int lev, #ifdef WARPX_QED void -PhotonParticleContainer::DoBreitWheeler(int lev, - amrex::Real dt, - const amrex::MultiFab& Ex, - const amrex::MultiFab& Ey, - const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, - const amrex::MultiFab& By, - const amrex::MultiFab& Bz) -{ -#ifdef _OPENMP - #pragma omp parallel -#endif - { -#ifdef _OPENMP - int thread_num = omp_get_thread_num(); -#else - int thread_num = 0; -#endif - - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) - { - const Box& box = pti.validbox(); - - auto& attribs = pti.GetAttribs(); - - auto& Exp = attribs[PIdx::Ex]; - auto& Eyp = attribs[PIdx::Ey]; - auto& Ezp = attribs[PIdx::Ez]; - auto& Bxp = attribs[PIdx::Bx]; - auto& Byp = attribs[PIdx::By]; - auto& Bzp = attribs[PIdx::Bz]; - - const long np = pti.numParticles(); - - // Data on the grid - const FArrayBox& exfab = Ex[pti]; - const FArrayBox& eyfab = Ey[pti]; - const FArrayBox& ezfab = Ez[pti]; - const FArrayBox& bxfab = Bx[pti]; - const FArrayBox& byfab = By[pti]; - const FArrayBox& bzfab = Bz[pti]; - - Exp.assign(np,WarpX::E_external_particle[0]); - Eyp.assign(np,WarpX::E_external_particle[1]); - Ezp.assign(np,WarpX::E_external_particle[2]); - - Bxp.assign(np,WarpX::B_external_particle[0]); - Byp.assign(np,WarpX::B_external_particle[1]); - Bzp.assign(np,WarpX::B_external_particle[2]); - - // - // copy data from particle container to temp arrays - // - pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - - int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, - Ex.nGrow(), e_is_nodal, - 0, np, thread_num, lev, lev); - - // This wraps the momentum advance so that inheritors can modify the call. - // Extract pointers to the different particle quantities - ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); - ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); - ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr(); - const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr(); - - BreitWheelerEvolveOpticalDepth evolve_opt = - m_shr_p_bw_engine->build_evolve_functor(); - - amrex::Real* AMREX_RESTRICT p_tau = - pti.GetAttribs(particle_comps["tau"]).dataPtr(); - - const auto me = PhysConst::m_e; - - // Loop over the particles and update their optical_depth - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - - const ParticleReal px = me * ux[i]; - const ParticleReal py = me * uy[i]; - const ParticleReal pz = me * uz[i]; - - evolve_opt( - px, py, pz, - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - dt, p_tau[i]); - } - ); - } - } -} - -void -PhotonParticleContainer::DoBreitWheelerPti(WarpXParIter& pti, +PhotonParticleContainer::DoBreitWheeler(WarpXParIter& pti, amrex::Real dt) { auto& attribs = pti.GetAttribs(); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 9e113a24b..fae4e0f4b 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -257,7 +257,7 @@ protected: std::shared_ptr m_shr_p_bw_engine; #endif -// PushPX and PushP has been split into two methods (classical and QED versions) +// PushPX has been split into two methods (classical and QED versions) // these methods must be public to be compatible with CUDA) public: //Classical versions @@ -291,33 +291,6 @@ public: amrex::Real q, amrex::Real m, int* AMREX_RESTRICT ion_lev, amrex::Real dt, long num_particles ); - - /** - * This function actually pushes momenta when QED effects - * are disabled. - * @param[in,out] x,y,z positions - * @param[in,out] ux,uy,uz momenta - * @param[in] Expp,Eypp,Ezpp electric fields - * @param[in] Bxpp,Bypp,Bzpp magnetic fields - * @param[in] q charge - * @param[in] m mass - * @param[in,out] ion_lev ionization level - * @param[in] dt temporal step - * @param[in] num_particles number of particles - */ - void PushP_classical( - amrex::ParticleReal* const AMREX_RESTRICT ux, - amrex::ParticleReal* const AMREX_RESTRICT uy, - amrex::ParticleReal* const AMREX_RESTRICT uz, - const amrex::ParticleReal* const AMREX_RESTRICT Expp, - const amrex::ParticleReal* const AMREX_RESTRICT Eypp, - const amrex::ParticleReal* const AMREX_RESTRICT Ezpp, - const amrex::ParticleReal* const AMREX_RESTRICT Bxpp, - const amrex::ParticleReal* const AMREX_RESTRICT Bypp, - const amrex::ParticleReal* const AMREX_RESTRICT Bzpp, - amrex::Real q, amrex::Real m, int* AMREX_RESTRICT ion_lev, - amrex::Real dt, long num_particles - ); //________________________ //QED versions @@ -353,34 +326,6 @@ public: amrex::ParticleReal* const AMREX_RESTRICT tau, amrex::Real dt, long num_particles ); - - /** - * This function actually pushes momenta when QED effects - * are enabled. - * @param[in,out] x,y,z positions - * @param[in,out] ux,uy,uz momenta - * @param[in] Expp,Eypp,Ezpp electric fields - * @param[in] Bxpp,Bypp,Bzpp magnetic fields - * @param[in] q charge - * @param[in] m mass - * @param[in,out] tau optical depth - * @param[in] dt temporal step - * @param[in] num_particles number of particles - */ - void PushP_QedQuantumSynchrotron( - amrex::ParticleReal* const AMREX_RESTRICT ux, - amrex::ParticleReal* const AMREX_RESTRICT uy, - amrex::ParticleReal* const AMREX_RESTRICT uz, - const amrex::ParticleReal* const AMREX_RESTRICT Expp, - const amrex::ParticleReal* const AMREX_RESTRICT Eypp, - const amrex::ParticleReal* const AMREX_RESTRICT Ezpp, - const amrex::ParticleReal* const AMREX_RESTRICT Bxpp, - const amrex::ParticleReal* const AMREX_RESTRICT Bypp, - const amrex::ParticleReal* const AMREX_RESTRICT Bzpp, - amrex::Real q, amrex::Real m, - amrex::ParticleReal* const AMREX_RESTRICT tau, - amrex::Real dt, long num_particles - ); //________________________ #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ee0711372..2fa80a7fc 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1911,214 +1911,60 @@ PhysicalParticleContainer::PushP (int lev, Real dt, ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } -#ifdef WARPX_QED - if(has_quantum_sync()){ - Real* AMREX_RESTRICT p_tau = - pti.GetAttribs(particle_comps["tau"]).dataPtr(); - PushP_QedQuantumSynchrotron( - ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, - q, m, p_tau, dt, pti.numParticles()); - } - else{ - PushP_classical( - ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, - q, m, ion_lev, dt, pti.numParticles()); - } -#else - PushP_classical( - ux, uy, uz, Expp, Eypp, Ezpp, Bxpp, Bypp, Bzpp, - q, m, ion_lev, dt, pti.numParticles()); -#endif - - - } - } -} - -void PhysicalParticleContainer::PushP_classical( - ParticleReal* const AMREX_RESTRICT ux, - ParticleReal* const AMREX_RESTRICT uy, - ParticleReal* const AMREX_RESTRICT uz, - const ParticleReal* const AMREX_RESTRICT Expp, - const ParticleReal* const AMREX_RESTRICT Eypp, - const ParticleReal* const AMREX_RESTRICT Ezpp, - const ParticleReal* const AMREX_RESTRICT Bxpp, - const ParticleReal* const AMREX_RESTRICT Bypp, - const ParticleReal* const AMREX_RESTRICT Bzpp, - Real q, Real m, int* AMREX_RESTRICT ion_lev, - Real dt, long num_particles -) -{ - //Assumes that all consistency checks have been done at initialization - if(do_classical_radiation_reaction){ - amrex::ParallelFor(num_particles, - [=] 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(num_particles, - [=] AMREX_GPU_DEVICE (long i) { - 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(num_particles, - [=] AMREX_GPU_DEVICE (long 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(num_particles, - [=] 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); - } - ); - } else { + //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) { + 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) { + 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); + } + ); + } else { amrex::Abort("Unknown particle pusher"); - } -} - -#ifdef WARPX_QED - -void PhysicalParticleContainer::PushP_QedQuantumSynchrotron( - ParticleReal* const AMREX_RESTRICT ux, - ParticleReal* const AMREX_RESTRICT uy, - ParticleReal* const AMREX_RESTRICT uz, - const ParticleReal* const AMREX_RESTRICT Expp, - const ParticleReal* const AMREX_RESTRICT Eypp, - const ParticleReal* const AMREX_RESTRICT Ezpp, - const ParticleReal* const AMREX_RESTRICT Bxpp, - const ParticleReal* const AMREX_RESTRICT Bypp, - const ParticleReal* const AMREX_RESTRICT Bzpp, - Real q, Real m, - ParticleReal* const AMREX_RESTRICT tau, - Real dt, long num_particles -) -{ - QuantumSynchrotronEvolveOpticalDepth evolve_opt = - m_shr_p_qs_engine->build_evolve_functor(); - - const auto chi_min = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; - - //Assumes that all consistency checks have been done at initialization - if(do_classical_radiation_reaction){ - amrex::ParallelFor(num_particles, - [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; - - const ParticleReal chi = QedUtils::chi_lepton(px, py, pz, - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i]); - - if(chi < chi_min){ - UpdateMomentumBorisWithRadiationReaction( - ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - q, m, dt); - } - else - { - bool has_event_happened = evolve_opt( - px, py, pz, - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - dt, tau[i]); - - UpdateMomentumBoris( - ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - q, m, dt); - } } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor(num_particles, - [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; - - bool has_event_happened = evolve_opt( - px, py, pz, - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - dt, tau[i]); - UpdateMomentumBoris( - ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - q, m, dt); - } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){ - amrex::ParallelFor(num_particles, - [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; - - bool has_event_happened = evolve_opt( - px, py, pz, - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - dt, tau[i]); - - UpdateMomentumVay( - ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - q, m, dt); - } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){ - amrex::ParallelFor(num_particles, - [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; - - bool has_event_happened = evolve_opt( - px, py, pz, - Expp[i], Eypp[i], Ezpp[i], - Bxpp[i], Bypp[i], Bzpp[i], - dt, tau[i]); - - UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); - } - ); - } else { - amrex::Abort("Unknown particle pusher"); + } } } -#endif - void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const ParticleReal* xp, const ParticleReal* yp, const ParticleReal* zp) { -- cgit v1.2.3 From f01e1f4a772644a0b81d054a6057c8cf56203f7d Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 31 Oct 2019 11:42:42 +0100 Subject: Allow initialization tests without QED --- Source/Particles/PhotonParticleContainer.cpp | 6 +++++- Source/Particles/PhysicalParticleContainer.cpp | 5 ++++- 2 files changed, 9 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index afe697b7e..7e05656ef 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -88,7 +88,11 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, ); #ifdef WARPX_QED - if(has_breit_wheeler()) DoBreitWheeler(pti, dt); + //m_shr_p_bw_engine->are_lookup_tables_initialized() is necessary here if we want + //to perform just initialization tests of the optical depth without actually + //enabling QED effects (this requires lookup tables). + if(has_breit_wheeler() && m_shr_p_bw_engine->are_lookup_tables_initialized()) + DoBreitWheeler(pti, dt); #endif } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 2fa80a7fc..123f4eb8e 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1601,7 +1601,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real m = this-> mass; #ifdef WARPX_QED - if(has_quantum_sync()){ + //m_shr_p_qs_engine->are_lookup_tables_initialized() is necessary here if we want + //to perform just initialization tests of the optical depth without actually + //enabling QED effects (this requires lookup tables). + if(has_quantum_sync()&& m_shr_p_qs_engine->are_lookup_tables_initialized()){ Real* AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); -- cgit v1.2.3 From d44b32225c22e870b909d4243168575ead339fc0 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 31 Oct 2019 12:44:40 +0100 Subject: Moved evolution of optical depth in a dedicated function --- Source/Particles/PhotonParticleContainer.H | 11 +- Source/Particles/PhotonParticleContainer.cpp | 20 ++- Source/Particles/PhysicalParticleContainer.H | 81 ++-------- Source/Particles/PhysicalParticleContainer.cpp | 195 +++++-------------------- 4 files changed, 64 insertions(+), 243 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index fbbc1077e..6a474a92a 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -81,11 +81,14 @@ public: #ifdef WARPX_QED /** - * This functions (called by PushPX) evolves optical depths - * of the photons according to the Breit Wheeler process + * This function evolves the optical depth of the photons if QED effects + * are enabled. + * @param[in,out] pti particle iterator (optical depth will be modified) + * @param[in] dt temporal step */ - void DoBreitWheeler(WarpXParIter& pti, - amrex::Real dt); + virtual void EvolveOpticalDepth(WarpXParIter& pti, + amrex::Real dt) override; + #endif }; diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 7e05656ef..7d63fff11 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -86,15 +86,6 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); - -#ifdef WARPX_QED - //m_shr_p_bw_engine->are_lookup_tables_initialized() is necessary here if we want - //to perform just initialization tests of the optical depth without actually - //enabling QED effects (this requires lookup tables). - if(has_breit_wheeler() && m_shr_p_bw_engine->are_lookup_tables_initialized()) - DoBreitWheeler(pti, dt); -#endif - } void @@ -124,10 +115,17 @@ PhotonParticleContainer::Evolve (int lev, } #ifdef WARPX_QED + void -PhotonParticleContainer::DoBreitWheeler(WarpXParIter& pti, - amrex::Real dt) +PhotonParticleContainer::EvolveOpticalDepth( + WarpXParIter& pti,amrex::Real dt) { + //m_shr_p_bw_engine->are_lookup_tables_initialized() is necessary here if we want + //to perform just initialization tests of the optical depth without actually + //enabling QED effects (this requires lookup tables). + if(!has_breit_wheeler() || !m_shr_p_bw_engine->are_lookup_tables_initialized()) + return; + auto& attribs = pti.GetAttribs(); ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index fae4e0f4b..4701c559d 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -219,6 +219,15 @@ public: void set_quantum_sync_engine_ptr (std::shared_ptr ptr) override; //__________ + + /** + * This function evolves the optical depth of the particles if QED effects + * are enabled. + * @param[in,out] pti particle iterator (optical depth will be modified) + * @param[in] dt temporal step + */ + virtual void EvolveOpticalDepth(WarpXParIter& pti, + amrex::Real dt); #endif protected: @@ -257,78 +266,6 @@ protected: std::shared_ptr m_shr_p_bw_engine; #endif -// PushPX has been split into two methods (classical and QED versions) -// these methods must be public to be compatible with CUDA) -public: - //Classical versions - - /** - * This function actually pushes momenta and positions when QED effects - * are disabled. - * @param[in,out] x,y,z positions - * @param[in,out] ux,uy,uz momenta - * @param[in] Ex,Ey,Ez electric fields - * @param[in] Bx,By,Bz magnetic fields - * @param[in] q charge - * @param[in] m mass - * @param[in,out] ion_lev ionization level - * @param[in] dt temporal step - * @param[in] num_particles number of particles - */ - void PushPX_classical( - amrex::ParticleReal* const AMREX_RESTRICT x, - amrex::ParticleReal* const AMREX_RESTRICT y, - amrex::ParticleReal* const AMREX_RESTRICT z, - amrex::ParticleReal* const AMREX_RESTRICT ux, - amrex::ParticleReal* const AMREX_RESTRICT uy, - amrex::ParticleReal* const AMREX_RESTRICT uz, - const amrex::ParticleReal* const AMREX_RESTRICT Ex, - const amrex::ParticleReal* const AMREX_RESTRICT Ey, - const amrex::ParticleReal* const AMREX_RESTRICT Ez, - const amrex::ParticleReal* const AMREX_RESTRICT Bx, - const amrex::ParticleReal* const AMREX_RESTRICT By, - const amrex::ParticleReal* const AMREX_RESTRICT Bz, - amrex::Real q, amrex::Real m, int* AMREX_RESTRICT ion_lev, - amrex::Real dt, long num_particles - ); - //________________________ - - //QED versions -#ifdef WARPX_QED - - /** - * This function actually pushes momenta and positions when QED effects - * are enabled. - * @param[in,out] x,y,z positions - * @param[in,out] ux,uy,uz momenta - * @param[in] Ex,Ey,Ez electric fields - * @param[in] Bx,By,Bz magnetic fields - * @param[in] q charge - * @param[in] m mass - * @param[in,out] tau optical depth - * @param[in] dt temporal step - * @param[in] num_particles number of particles - */ - void PushPX_QedQuantumSynchrotron( - amrex::ParticleReal* const AMREX_RESTRICT x, - amrex::ParticleReal* const AMREX_RESTRICT y, - amrex::ParticleReal* const AMREX_RESTRICT z, - amrex::ParticleReal* const AMREX_RESTRICT ux, - amrex::ParticleReal* const AMREX_RESTRICT uy, - amrex::ParticleReal* const AMREX_RESTRICT uz, - const amrex::ParticleReal* const AMREX_RESTRICT Ex, - const amrex::ParticleReal* const AMREX_RESTRICT Ey, - const amrex::ParticleReal* const AMREX_RESTRICT Ez, - const amrex::ParticleReal* const AMREX_RESTRICT Bx, - const amrex::ParticleReal* const AMREX_RESTRICT By, - const amrex::ParticleReal* const AMREX_RESTRICT Bz, - amrex::Real q, amrex::Real m, - amrex::ParticleReal* const AMREX_RESTRICT tau, - amrex::Real dt, long num_particles - ); - //________________________ -#endif - }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 123f4eb8e..2b3e0d94d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1062,6 +1062,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE("PPC::Evolve()"); BL_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy); BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); + BL_PROFILE_VAR_NS("PPC::EvolveOpticalDepth", blp_ppc_qed_ev); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); const std::array& dx = WarpX::CellSize(lev); @@ -1253,6 +1254,15 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_fg); +#ifdef WARPX_QED + // + //Evolve Optical Depth + // + BL_PROFILE_VAR_START(blp_ppc_qed_ev); + EvolveOpticalDepth(pti, dt); + BL_PROFILE_VAR_STOP(blp_ppc_qed_ev); +#endif + // // Particle Push // @@ -1600,51 +1610,10 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real q = this->charge; const Real m = this-> mass; -#ifdef WARPX_QED - //m_shr_p_qs_engine->are_lookup_tables_initialized() is necessary here if we want - //to perform just initialization tests of the optical depth without actually - //enabling QED effects (this requires lookup tables). - if(has_quantum_sync()&& m_shr_p_qs_engine->are_lookup_tables_initialized()){ - Real* AMREX_RESTRICT p_tau = - pti.GetAttribs(particle_comps["tau"]).dataPtr(); - - PushPX_QedQuantumSynchrotron(x, y, z, ux, uy, uz, - Ex, Ey, Ez, Bx, By, Bz, - q, m, p_tau, dt, pti.numParticles()); - } - else { - PushPX_classical(x, y, z, ux, uy, uz, - Ex, Ey, Ez, Bx, By, Bz, - q, m, ion_lev, dt, pti.numParticles()); - } -#else - PushPX_classical(x, y, z, ux, uy, uz, - Ex, Ey, Ez, Bx, By, Bz, - q, m, ion_lev, dt, pti.numParticles()); -#endif -} - -void PhysicalParticleContainer::PushPX_classical( - ParticleReal* const AMREX_RESTRICT x, - ParticleReal* const AMREX_RESTRICT y, - ParticleReal* const AMREX_RESTRICT z, - ParticleReal* const AMREX_RESTRICT ux, - ParticleReal* const AMREX_RESTRICT uy, - ParticleReal* const AMREX_RESTRICT uz, - const ParticleReal* const AMREX_RESTRICT Ex, - const ParticleReal* const AMREX_RESTRICT Ey, - const ParticleReal* const AMREX_RESTRICT Ez, - const ParticleReal* const AMREX_RESTRICT Bx, - const ParticleReal* const AMREX_RESTRICT By, - const ParticleReal* const AMREX_RESTRICT Bz, - Real q, Real m, int* AMREX_RESTRICT ion_lev, - amrex::Real dt, long num_particles -) -{ //Assumes that all consistency checks have been done at initialization if(do_classical_radiation_reaction){ amrex::ParallelFor( - num_particles, + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1657,7 +1626,7 @@ void PhysicalParticleContainer::PushPX_classical( ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( - num_particles, + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1670,7 +1639,7 @@ void PhysicalParticleContainer::PushPX_classical( ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { amrex::ParallelFor( - num_particles, + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1683,7 +1652,7 @@ void PhysicalParticleContainer::PushPX_classical( ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { amrex::ParallelFor( - num_particles, + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } @@ -1696,94 +1665,39 @@ void PhysicalParticleContainer::PushPX_classical( ); } else { amrex::Abort("Unknown particle pusher"); - } + }; } #ifdef WARPX_QED -void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( - ParticleReal* const AMREX_RESTRICT x, - ParticleReal* const AMREX_RESTRICT y, - ParticleReal* const AMREX_RESTRICT z, - ParticleReal* const AMREX_RESTRICT ux, - ParticleReal* const AMREX_RESTRICT uy, - ParticleReal* const AMREX_RESTRICT uz, - const ParticleReal* const AMREX_RESTRICT Ex, - const ParticleReal* const AMREX_RESTRICT Ey, - const ParticleReal* const AMREX_RESTRICT Ez, - const ParticleReal* const AMREX_RESTRICT Bx, - const ParticleReal* const AMREX_RESTRICT By, - const ParticleReal* const AMREX_RESTRICT Bz, - Real q, amrex::Real m, - ParticleReal* const AMREX_RESTRICT tau, - Real dt, long num_particles -) +void PhysicalParticleContainer::EvolveOpticalDepth( + WarpXParIter& pti, amrex::Real dt) { + //m_shr_p_qs_engine->are_lookup_tables_initialized() is necessary here if we want + //to perform just initialization tests of the optical depth without actually + //enabling QED effects (this requires lookup tables). + if(!has_quantum_sync() || !m_shr_p_qs_engine->are_lookup_tables_initialized()) + return; + QuantumSynchrotronEvolveOpticalDepth evolve_opt = m_shr_p_qs_engine->build_evolve_functor(); - const auto chi_min = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; - - //Assumes that all consistency checks have been done at initialization - if(do_classical_radiation_reaction){ - amrex::ParallelFor(num_particles, - [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; - - const ParticleReal chi = QedUtils::chi_lepton(px, py, pz, - Ex[i], Ey[i], Ez[i], - Bx[i], By[i], Bz[i]); - - if(chi < chi_min){ - UpdateMomentumBorisWithRadiationReaction( - ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], - Bx[i], By[i], Bz[i], - q, m, dt); - } - else - { - bool has_event_happened = evolve_opt( - px, py, pz, - Ex[i], Ey[i], Ez[i], - Bx[i], By[i], Bz[i], - dt, tau[i]); - - UpdateMomentumBoris( - ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], - Bx[i], By[i], Bz[i], - q, 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(num_particles, - [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; + auto& attribs = pti.GetAttribs(); + const ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + const ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + const ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr(); + const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); + const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); - bool has_event_happened = evolve_opt( - px, py, pz, - Ex[i], Ey[i], Ez[i], - Bx[i], By[i], Bz[i], - dt, tau[i]); + ParticleReal* const AMREX_RESTRICT p_tau = + pti.GetAttribs(particle_comps["tau"]).dataPtr(); - UpdateMomentumBoris( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); + const ParticleReal m = this->mass; - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); - } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { - amrex::ParallelFor(num_particles, + amrex::ParallelFor(pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { const ParticleReal px = m * ux[i]; const ParticleReal py = m * uy[i]; @@ -1793,41 +1707,10 @@ void PhysicalParticleContainer::PushPX_QedQuantumSynchrotron( px, py, pz, Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], - dt, tau[i]); - - UpdateMomentumVay( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); - - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + dt, p_tau[i]); } - ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) { - amrex::ParallelFor( - num_particles, - [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; - - bool has_event_happened = evolve_opt( - px, py, pz, - Ex[i], Ey[i], Ez[i], - Bx[i], By[i], Bz[i], - dt, tau[i]); - - UpdateMomentumHigueraCary( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], q, m, dt); + ); - UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); - } - ); - } else { - amrex::Abort("Unknown particle pusher"); - } } #endif -- cgit v1.2.3 From a1fab0ee7b8e864516d69b66edd1a26bf9f2aefa Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 5 Nov 2019 19:26:10 +0100 Subject: added missing commit --- Source/Particles/MultiParticleContainer.H | 1 + Source/Particles/PhysicalParticleContainer.cpp | 41 +++++++++++++++++++++++++- 2 files changed, 41 insertions(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index f9a0e51d7..06898a5c1 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -11,6 +11,7 @@ #include #ifdef WARPX_QED + #include #include #include #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 2b3e0d94d..54cc1a880 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1610,7 +1610,45 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real q = this->charge; const Real m = this-> mass; - //Assumes that all consistency checks have been done at initialization +#ifdef WARPX_QED + + auto t_chi_max = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min; + + if(do_classical_radiation_reaction){ + if(m_do_qed_quantum_sync){ + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + auto chi = QedUtils::chi_lepton(m*ux[i], m*uy[i], m*uz[i], + Ex[i], Ey[i], Ez[i], + Bx[i], By[i], Bz[i]); + if(chi < t_chi_max){ + UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + } + else{ + UpdateMomentumBoris( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + } + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + }else{ + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], q, m, dt); + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); + } +#else if(do_classical_radiation_reaction){ amrex::ParallelFor( pti.numParticles(), @@ -1624,6 +1662,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); +#endif } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), -- cgit v1.2.3 From 7220ab302e9a62116f52f0345018d48d5884bca4 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 6 Nov 2019 02:43:07 +0100 Subject: fixed bug --- Source/Particles/PhotonParticleContainer.cpp | 8 ++++---- Source/Particles/PhysicalParticleContainer.cpp | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 7d63fff11..001f748d0 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -143,14 +143,14 @@ PhotonParticleContainer::EvolveOpticalDepth( amrex::Real* AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); - const auto me = PhysConst::m_e; + const auto mec = PhysConst::m_e * PhysConst::c; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = me * ux[i]; - const ParticleReal py = me * uy[i]; - const ParticleReal pz = me * uz[i]; + const ParticleReal px = mec * ux[i]; + const ParticleReal py = mec * uy[i]; + const ParticleReal pz = mec * uz[i]; bool has_event_happened = evolve_opt( px, py, pz, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 54cc1a880..cf609308b 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1734,13 +1734,13 @@ void PhysicalParticleContainer::EvolveOpticalDepth( ParticleReal* const AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); - const ParticleReal m = this->mass; + const auto mc = PhysConst::c * this->mass; amrex::ParallelFor(pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = m * ux[i]; - const ParticleReal py = m * uy[i]; - const ParticleReal pz = m * uz[i]; + const ParticleReal px = mc * ux[i]; + const ParticleReal py = mc * uy[i]; + const ParticleReal pz = mc * uz[i]; bool has_event_happened = evolve_opt( px, py, pz, -- cgit v1.2.3 From 5b7ae4e641721564eaf6a29699f945998e4be5c0 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 6 Nov 2019 11:09:24 +0100 Subject: Revert "fixed bug" This reverts commit 7220ab302e9a62116f52f0345018d48d5884bca4. --- Source/Particles/PhotonParticleContainer.cpp | 8 ++++---- Source/Particles/PhysicalParticleContainer.cpp | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 001f748d0..7d63fff11 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -143,14 +143,14 @@ PhotonParticleContainer::EvolveOpticalDepth( amrex::Real* AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); - const auto mec = PhysConst::m_e * PhysConst::c; + const auto me = PhysConst::m_e; amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = mec * ux[i]; - const ParticleReal py = mec * uy[i]; - const ParticleReal pz = mec * uz[i]; + const ParticleReal px = me * ux[i]; + const ParticleReal py = me * uy[i]; + const ParticleReal pz = me * uz[i]; bool has_event_happened = evolve_opt( px, py, pz, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index cf609308b..54cc1a880 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1734,13 +1734,13 @@ void PhysicalParticleContainer::EvolveOpticalDepth( ParticleReal* const AMREX_RESTRICT p_tau = pti.GetAttribs(particle_comps["tau"]).dataPtr(); - const auto mc = PhysConst::c * this->mass; + const ParticleReal m = this->mass; amrex::ParallelFor(pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - const ParticleReal px = mc * ux[i]; - const ParticleReal py = mc * uy[i]; - const ParticleReal pz = mc * uz[i]; + const ParticleReal px = m * ux[i]; + const ParticleReal py = m * uy[i]; + const ParticleReal pz = m * uz[i]; bool has_event_happened = evolve_opt( px, py, pz, -- cgit v1.2.3 From cec596f7936e98e0b6a71ef4ed3137c0dfe13758 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Nov 2019 12:20:44 +0100 Subject: cleaning --- Source/Particles/PhotonParticleContainer.cpp | 5 +---- Source/Particles/PhysicalParticleContainer.cpp | 5 +---- 2 files changed, 2 insertions(+), 8 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 424569c9a..c03ed00c2 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -120,10 +120,7 @@ void PhotonParticleContainer::EvolveOpticalDepth( WarpXParIter& pti,amrex::Real dt) { - //m_shr_p_bw_engine->are_lookup_tables_initialized() is necessary here if we want - //to perform just initialization tests of the optical depth without actually - //enabling QED effects (this requires lookup tables). - if(!has_breit_wheeler() || !m_shr_p_bw_engine->are_lookup_tables_initialized()) + if(!has_breit_wheeler()) return; auto& attribs = pti.GetAttribs(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c67f9d27d..ef8115569 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1701,10 +1701,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, void PhysicalParticleContainer::EvolveOpticalDepth( WarpXParIter& pti, amrex::Real dt) { - //m_shr_p_qs_engine->are_lookup_tables_initialized() is necessary here if we want - //to perform just initialization tests of the optical depth without actually - //enabling QED effects (this requires lookup tables). - if(!has_quantum_sync() || !m_shr_p_qs_engine->are_lookup_tables_initialized()) + if(!has_quantum_sync()) return; QuantumSynchrotronEvolveOpticalDepth evolve_opt = -- cgit v1.2.3 From 6f74a7ac262c26a276fe4b818ceb1e8ca26e5575 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 15 Nov 2019 11:22:50 -0800 Subject: Add option to set relative precision --- Docs/source/running_cpp/parameters.rst | 8 ++++++++ Examples/Modules/relativistic_space_charge_initialization/inputs | 1 + Source/Initialization/InitSpaceChargeField.cpp | 8 ++++---- Source/Particles/PhysicalParticleContainer.cpp | 1 + Source/Particles/WarpXParticleContainer.H | 1 + Source/WarpX.H | 3 ++- 6 files changed, 17 insertions(+), 5 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 50803560d..2aee22344 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -228,6 +228,14 @@ Particle initialization Whether to calculate the space-charge fields associated with this species at the beginning of the simulation. +* ``.self_fields_required_precision`` (`float`, default: 1.e-11) + The relative precision with which the initial space-charge fields should + be calculated. More specifically, the initial space-charge fields are + computed with an iterative Multi-Level Multi-Grid (MLMG) solver. + For highly-relativistic beams, this solver can fail to reach the default + precision within a reasonable time ; in that case, users can set a + relaxed precision requirement through ``self_fields_required_precision``. + * ``.profile`` (`string`) Density profile for this species. The options are: diff --git a/Examples/Modules/relativistic_space_charge_initialization/inputs b/Examples/Modules/relativistic_space_charge_initialization/inputs index 6d936aa27..eca38e074 100644 --- a/Examples/Modules/relativistic_space_charge_initialization/inputs +++ b/Examples/Modules/relativistic_space_charge_initialization/inputs @@ -15,6 +15,7 @@ beam.charge = -q_e beam.mass = m_e beam.injection_style = "gaussian_beam" beam.initialize_self_fields = 1 +beam.self_fields_required_precision = 1.e-4 beam.x_rms = 2.e-6 beam.y_rms = 2.e-6 beam.z_rms = 2.e-6 diff --git a/Source/Initialization/InitSpaceChargeField.cpp b/Source/Initialization/InitSpaceChargeField.cpp index f0225e5d3..2ae3e181a 100644 --- a/Source/Initialization/InitSpaceChargeField.cpp +++ b/Source/Initialization/InitSpaceChargeField.cpp @@ -39,7 +39,7 @@ WarpX::InitSpaceChargeField (WarpXParticleContainer& pc) for (Real& beta_comp : beta) beta_comp /= PhysConst::c; // Normalize // Compute the potential phi, by solving the Poisson equation - computePhi( rho, phi, beta ); + computePhi( rho, phi, beta, pc.self_fields_required_precision ); // Compute the corresponding electric and magnetic field, from the potential phi computeE( Efield_fp, phi, beta ); @@ -63,7 +63,8 @@ WarpX::InitSpaceChargeField (WarpXParticleContainer& pc) void WarpX::computePhi (const amrex::Vector >& rho, amrex::Vector >& phi, - std::array const beta) const + std::array const beta, + Real const required_precision) const { // Define the boundary conditions Array lobc, hibc; @@ -94,8 +95,7 @@ WarpX::computePhi (const amrex::Vector >& rho, // Solve the Poisson equation MLMG mlmg(linop); mlmg.setVerbose(2); - const Real reltol = 1.e-11; - mlmg.solve( GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), reltol, 0.0); + mlmg.solve( GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), required_precision, 0.0); // Normalize by the correct physical constant for (int lev=0; lev < rho.size(); lev++){ diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index cb2a7c0e8..9dc9db2c5 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -42,6 +42,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_continuous_injection", do_continuous_injection); pp.query("initialize_self_fields", initialize_self_fields); + pp.query("self_fields_required_precision", self_fields_required_precision); // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 7da25b452..49fa09262 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -255,6 +255,7 @@ public: bool do_splitting = false; bool initialize_self_fields = false; + bool self_fields_required_precision = 1.e-11; // split along diagonals (0) or axes (1) int split_type = 0; diff --git a/Source/WarpX.H b/Source/WarpX.H index b04f56993..37e2427a5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -305,7 +305,8 @@ public: void InitSpaceChargeField (WarpXParticleContainer& pc); void computePhi (const amrex::Vector >& rho, amrex::Vector >& phi, - std::array const beta = {{0,0,0}} ) const; + std::array const beta = {{0,0,0}}, + amrex::Real const required_precision=1.e-11 ) const; void computeE (amrex::Vector, 3> >& E, const amrex::Vector >& phi, std::array const beta = {{0,0,0}} ) const; -- cgit v1.2.3 From 25a614fe9e3bdbae8008b408976e65232286b1c5 Mon Sep 17 00:00:00 2001 From: Yinjian Zhao Date: Mon, 25 Nov 2019 11:16:01 -0700 Subject: Change do_not_depoist to be species related. --- Docs/source/running_cpp/parameters.rst | 8 ++++---- Source/Particles/PhysicalParticleContainer.cpp | 1 + Source/Particles/WarpXParticleContainer.cpp | 1 - 3 files changed, 5 insertions(+), 5 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 59c4be5fd..a16e0bfec 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -363,15 +363,15 @@ Particle initialization Split particles of the species when crossing the boundary from a lower resolution domain to a higher resolution domain. -* ``.do_not_deposit`` (`0` or `1` optional; default `0`) - If `1` is given, both charge deposition and current deposition will - not be done, so the fields due to particles remain zero. - * ``.split_type`` (`int`) optional (default `0`) Splitting technique. When `0`, particles are split along the simulation axes (4 particles in 2D, 6 particles in 3D). When `1`, particles are split along the diagonals (4 particles in 2D, 8 particles in 3D). +* ``.do_not_deposit`` (`0` or `1` optional; default `0`) + If `1` is given, both charge deposition and current deposition will + not be done, thus that species does not contribute to the fields. + * ``.plot_species`` (`0` or `1` optional; default `1`) Whether to plot particle quantities for this species. diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4f6d032e9..55290cf1f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -39,6 +39,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // Initialize splitting pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); + pp.query("do_not_deposit", do_not_deposit); pp.query("do_continuous_injection", do_continuous_injection); pp.query("initialize_self_fields", initialize_self_fields); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c256df6aa..7b5bf16d1 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -122,7 +122,6 @@ WarpXParticleContainer::ReadParameters () #endif pp.query("do_tiling", do_tiling); pp.query("do_not_push", do_not_push); - pp.query("do_not_deposit", do_not_deposit); initialized = true; } -- cgit v1.2.3 From a2c70d824036f371c07e670cb0aea6dde1dc383f Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 4 Dec 2019 11:31:32 -0800 Subject: When momentum-conserving gather is used, the index type of Ex, Ey, etc. is not the same as WarpX::*_nodal_flag. --- Source/Particles/PhysicalParticleContainer.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 008eb9c89..3128868e4 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1362,40 +1362,41 @@ PhysicalParticleContainer::applyNCIFilter ( #endif // Filter Ex (Both 2D and 3D) - filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + filtered_Ex.resize(amrex::convert(tbox,Ex.box().ixType())); // Safeguard for GPU exeli = filtered_Ex.elixir(); // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex, filtered_Ex.box()); // Update ex_ptr reference ex_ptr = &filtered_Ex; // Filter Ez - filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); + filtered_Ez.resize(amrex::convert(tbox,Ez.box().ixType())); ezeli = filtered_Ez.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez, filtered_Ez.box()); ez_ptr = &filtered_Ez; // Filter By - filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); + filtered_By.resize(amrex::convert(tbox,By.box().ixType())); byeli = filtered_By.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By, filtered_By.box()); by_ptr = &filtered_By; #if (AMREX_SPACEDIM == 3) // Filter Ey - filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); + filtered_Ey.resize(amrex::convert(tbox,Ey.box().ixType())); eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey, filtered_Ey.box()); ey_ptr = &filtered_Ey; // Filter Bx - filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); + filtered_Bx.resize(amrex::convert(tbox,Bx.box().ixType())); bxeli = filtered_Bx.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx, filtered_Bx.box()); bx_ptr = &filtered_Bx; // Filter Bz - filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); + filtered_Bz.resize(amrex::convert(tbox,Bz.box().ixType())); bzeli = filtered_Bz.elixir(); nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz, filtered_Bz.box()); bz_ptr = &filtered_Bz; -- cgit v1.2.3 From 8a7cd073a469487d7c8e88a3c13c7ae0e2aed588 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 5 Dec 2019 09:24:25 -0800 Subject: Updated field gather for generalized centering --- Source/Particles/Gather/FieldGather.H | 238 ++++++++++++++++++------- Source/Particles/PhysicalParticleContainer.cpp | 33 ++-- 2 files changed, 182 insertions(+), 89 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index bdeab16e8..b4fc84908 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -25,17 +25,16 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, - const amrex::Array4& ex_arr, - const amrex::Array4& ey_arr, - const amrex::Array4& ez_arr, - const amrex::Array4& bx_arr, - const amrex::Array4& by_arr, - const amrex::Array4& bz_arr, + amrex::FArrayBox const * const exfab, + amrex::FArrayBox const * const eyfab, + amrex::FArrayBox const * const ezfab, + amrex::FArrayBox const * const bxfab, + amrex::FArrayBox const * const byfab, + amrex::FArrayBox const * const bzfab, const long np_to_gather, const std::array& dx, const std::array xyzmin, const amrex::Dim3 lo, - const amrex::Real stagger_shift, const long n_rz_azimuthal_modes) { const amrex::Real dxi = 1.0/dx[0]; @@ -50,6 +49,24 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, #endif const amrex::Real zmin = xyzmin[2]; + amrex::Array4 const& ex_arr = exfab->array(); + amrex::Array4 const& ey_arr = eyfab->array(); + amrex::Array4 const& ez_arr = ezfab->array(); + amrex::Array4 const& bx_arr = bxfab->array(); + amrex::Array4 const& by_arr = byfab->array(); + amrex::Array4 const& bz_arr = bzfab->array(); + + amrex::IntVect const ex_type = exfab->box().type(); + amrex::IntVect const ey_type = eyfab->box().type(); + amrex::IntVect const ez_type = ezfab->box().type(); + amrex::IntVect const bx_type = bxfab->box().type(); + amrex::IntVect const by_type = byfab->box().type(); + amrex::IntVect const bz_type = bzfab->box().type(); + + constexpr int zdir = (AMREX_SPACEDIM - 1); + constexpr int NODE = amrex::IndexType::NODE; + constexpr int CELL = amrex::IndexType::CELL; + // Loop over particles and gather fields from // {e,b}{x,y,z}_arr to {E,B}{xyz}p. amrex::ParallelFor( @@ -64,31 +81,116 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, #else const amrex::Real x = (xp[ip]-xmin)*dxi; #endif - // Compute shape factors for node-centered quantities - amrex::Real sx [depos_order + 1]; - // j: leftmost grid point (node-centered) that particle touches - const int j = compute_shape_factor(sx, x); - // Compute shape factors for cell-centered quantities - amrex::Real sx0[depos_order + 1 - lower_in_v]; - // j0: leftmost grid point (cell-centered) that particle touches - const int j0 = compute_shape_factor( - sx0, x-stagger_shift); + + // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current + // sx_[eb][xyz] shape factor along x for the centering of each current + // There are only two possible centerings, node or cell centered, so at most only two shape factor + // arrays will be needed. + amrex::Real sx_node[depos_order + 1]; + amrex::Real sx_cell[depos_order + 1]; + amrex::Real sx_node_v[depos_order + 1 - lower_in_v]; + amrex::Real sx_cell_v[depos_order + 1 - lower_in_v]; + int j_node; + int j_cell; + int j_node_v; + int j_cell_v; + if ((ey_type[0] == NODE) || (ez_type[0] == NODE) || (bx_type[0] == NODE)) { + j_node = compute_shape_factor(sx_node, x); + } + if ((ey_type[0] == CELL) || (ez_type[0] == CELL) || (bx_type[0] == CELL)) { + j_cell = compute_shape_factor(sx_cell, x - 0.5); + } + if ((ex_type[0] == NODE) || (by_type[0] == NODE) || (bz_type[0] == NODE)) { + j_node_v = compute_shape_factor(sx_node_v, x); + } + if ((ex_type[0] == CELL) || (by_type[0] == CELL) || (bz_type[0] == CELL)) { + j_cell_v = compute_shape_factor(sx_cell_v, x - 0.5); + } + const amrex::Real (&sx_ex)[depos_order + 1 - lower_in_v] = ((ex_type[0] == NODE) ? sx_node_v : sx_cell_v); + const amrex::Real (&sx_ey)[depos_order + 1 ] = ((ey_type[0] == NODE) ? sx_node : sx_cell ); + const amrex::Real (&sx_ez)[depos_order + 1 ] = ((ez_type[0] == NODE) ? sx_node : sx_cell ); + const amrex::Real (&sx_bx)[depos_order + 1 ] = ((bx_type[0] == NODE) ? sx_node : sx_cell ); + const amrex::Real (&sx_by)[depos_order + 1 - lower_in_v] = ((by_type[0] == NODE) ? sx_node_v : sx_cell_v); + const amrex::Real (&sx_bz)[depos_order + 1 - lower_in_v] = ((bz_type[0] == NODE) ? sx_node_v : sx_cell_v); + int const j_ex = ((ex_type[0] == NODE) ? j_node_v : j_cell_v); + int const j_ey = ((ey_type[0] == NODE) ? j_node : j_cell ); + int const j_ez = ((ez_type[0] == NODE) ? j_node : j_cell ); + int const j_bx = ((bx_type[0] == NODE) ? j_node : j_cell ); + int const j_by = ((by_type[0] == NODE) ? j_node_v : j_cell_v); + int const j_bz = ((bz_type[0] == NODE) ? j_node_v : j_cell_v); + #if (AMREX_SPACEDIM == 3) // y direction const amrex::Real y = (yp[ip]-ymin)*dyi; - amrex::Real sy [depos_order + 1]; - const int k = compute_shape_factor(sy, y); - amrex::Real sy0[depos_order + 1 - lower_in_v]; - const int k0 = compute_shape_factor( - sy0, y-stagger_shift); + amrex::Real sy_node[depos_order + 1]; + amrex::Real sy_cell[depos_order + 1]; + amrex::Real sy_node_v[depos_order + 1 - lower_in_v]; + amrex::Real sy_cell_v[depos_order + 1 - lower_in_v]; + int k_node; + int k_cell; + int k_node_v; + int k_cell_v; + if ((ex_type[1] == NODE) || (ez_type[1] == NODE) || (by_type[1] == NODE)) { + k_node = compute_shape_factor(sy_node, y); + } + if ((ex_type[1] == CELL) || (ez_type[1] == CELL) || (by_type[1] == CELL)) { + k_cell = compute_shape_factor(sy_cell, y - 0.5); + } + if ((ey_type[1] == NODE) || (bx_type[1] == NODE) || (bz_type[1] == NODE)) { + k_node_v = compute_shape_factor(sy_node_v, y); + } + if ((ey_type[1] == CELL) || (bx_type[1] == CELL) || (bz_type[1] == CELL)) { + k_cell_v = compute_shape_factor(sy_cell_v, y - 0.5); + } + const amrex::Real (&sy_ex)[depos_order + 1 ] = ((ex_type[1] == NODE) ? sy_node : sy_cell ); + const amrex::Real (&sy_ey)[depos_order + 1 - lower_in_v] = ((ey_type[1] == NODE) ? sy_node_v : sy_cell_v); + const amrex::Real (&sy_ez)[depos_order + 1 ] = ((ez_type[1] == NODE) ? sy_node : sy_cell ); + const amrex::Real (&sy_bx)[depos_order + 1 - lower_in_v] = ((bx_type[1] == NODE) ? sy_node_v : sy_cell_v); + const amrex::Real (&sy_by)[depos_order + 1 ] = ((by_type[1] == NODE) ? sy_node : sy_cell ); + const amrex::Real (&sy_bz)[depos_order + 1 - lower_in_v] = ((bz_type[1] == NODE) ? sy_node_v : sy_cell_v); + int const k_ex = ((ex_type[1] == NODE) ? k_node : k_cell ); + int const k_ey = ((ey_type[1] == NODE) ? k_node_v : k_cell_v); + int const k_ez = ((ez_type[1] == NODE) ? k_node : k_cell ); + int const k_bx = ((bx_type[1] == NODE) ? k_node_v : k_cell_v); + int const k_by = ((by_type[1] == NODE) ? k_node : k_cell ); + int const k_bz = ((bz_type[1] == NODE) ? k_node_v : k_cell_v); + #endif // z direction const amrex::Real z = (zp[ip]-zmin)*dzi; - amrex::Real sz [depos_order + 1]; - const int l = compute_shape_factor(sz, z); - amrex::Real sz0[depos_order + 1 - lower_in_v]; - const int l0 = compute_shape_factor( - sz0, z-stagger_shift); + amrex::Real sz_node[depos_order + 1]; + amrex::Real sz_cell[depos_order + 1]; + amrex::Real sz_node_v[depos_order + 1 - lower_in_v]; + amrex::Real sz_cell_v[depos_order + 1 - lower_in_v]; + int l_node; + int l_cell; + int l_node_v; + int l_cell_v; + if ((ex_type[zdir] == NODE) || (ey_type[zdir] == NODE) || (bz_type[zdir] == NODE)) { + l_node = compute_shape_factor(sz_node, z); + } + if ((ex_type[zdir] == CELL) || (ey_type[zdir] == CELL) || (bz_type[zdir] == CELL)) { + l_cell = compute_shape_factor(sz_cell, z - 0.5); + } + if ((ez_type[zdir] == NODE) || (bx_type[zdir] == NODE) || (by_type[zdir] == NODE)) { + l_node_v = compute_shape_factor(sz_node_v, z); + } + if ((ez_type[zdir] == CELL) || (bx_type[zdir] == CELL) || (by_type[zdir] == CELL)) { + l_cell_v = compute_shape_factor(sz_cell_v, z - 0.5); + } + const amrex::Real (&sz_ex)[depos_order + 1 ] = ((ex_type[zdir] == NODE) ? sz_node : sz_cell ); + const amrex::Real (&sz_ey)[depos_order + 1 ] = ((ey_type[zdir] == NODE) ? sz_node : sz_cell ); + const amrex::Real (&sz_ez)[depos_order + 1 - lower_in_v] = ((ez_type[zdir] == NODE) ? sz_node_v : sz_cell_v); + const amrex::Real (&sz_bx)[depos_order + 1 - lower_in_v] = ((bx_type[zdir] == NODE) ? sz_node_v : sz_cell_v); + const amrex::Real (&sz_by)[depos_order + 1 - lower_in_v] = ((by_type[zdir] == NODE) ? sz_node_v : sz_cell_v); + const amrex::Real (&sz_bz)[depos_order + 1 ] = ((bz_type[zdir] == NODE) ? sz_node : sz_cell ); + int const l_ex = ((ex_type[zdir] == NODE) ? l_node : l_cell ); + int const l_ey = ((ey_type[zdir] == NODE) ? l_node : l_cell ); + int const l_ez = ((ez_type[zdir] == NODE) ? l_node_v : l_cell_v); + int const l_bx = ((bx_type[zdir] == NODE) ? l_node_v : l_cell_v); + int const l_by = ((by_type[zdir] == NODE) ? l_node_v : l_cell_v); + int const l_bz = ((bz_type[zdir] == NODE) ? l_node : l_cell ); + // Each field is gathered in a separate block of // AMREX_SPACEDIM nested loops because the deposition @@ -98,35 +200,35 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, // Gather field on particle Eyp[i] from field on grid ey_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ - Eyp[ip] += sx[ix]*sz[iz]* - ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 0); + Eyp[ip] += sx_ey[ix]*sz_ey[iz]* + ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 0); } } // Gather field on particle Exp[i] from field on grid ex_arr // Gather field on particle Bzp[i] from field on grid bz_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - Exp[ip] += sx0[ix]*sz[iz]* - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 0); - Bzp[ip] += sx0[ix]*sz[iz]* - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 0); + Exp[ip] += sx_ex[ix]*sz_ex[iz]* + ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 0); + Bzp[ip] += sx_bz[ix]*sz_bz[iz]* + bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 0); } } // Gather field on particle Ezp[i] from field on grid ez_arr // Gather field on particle Bxp[i] from field on grid bx_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order; ix++){ - Ezp[ip] += sx[ix]*sz0[iz]* - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 0); - Bxp[ip] += sx[ix]*sz0[iz]* - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 0); + Ezp[ip] += sx_ez[ix]*sz_ez[iz]* + ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 0); + Bxp[ip] += sx_bx[ix]*sz_bx[iz]* + bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 0); } } // Gather field on particle Byp[i] from field on grid by_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - Byp[ip] += sx0[ix]*sz0[iz]* - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 0); + Byp[ip] += sx_by[ix]*sz_by[iz]* + by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 0); } } @@ -149,41 +251,41 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, // Gather field on particle Eyp[i] from field on grid ey_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ - const amrex::Real dEy = (+ ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode-1)*xy.real() - - ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode)*xy.imag()); - Eyp[ip] += sx[ix]*sz[iz]*dEy; + const amrex::Real dEy = (+ ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode-1)*xy.real() + - ey_arr(lo.x+j_ey+ix, lo.y+l_ey+iz, 0, 2*imode)*xy.imag()); + Eyp[ip] += sx_ey[ix]*sz_ey[iz]*dEy; } } // Gather field on particle Exp[i] from field on grid ex_arr // Gather field on particle Bzp[i] from field on grid bz_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - const amrex::Real dEx = (+ ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode-1)*xy.real() - - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.imag()); - Exp[ip] += sx0[ix]*sz[iz]*dEx; - const amrex::Real dBz = (+ bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode-1)*xy.real() - - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.imag()); - Bzp[ip] += sx0[ix]*sz[iz]*dBz; + const amrex::Real dEx = (+ ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode-1)*xy.real() + - ex_arr(lo.x+j_ex+ix, lo.y+l_ex+iz, 0, 2*imode)*xy.imag()); + Exp[ip] += sx_ex[ix]*sz_ex[iz]*dEx; + const amrex::Real dBz = (+ bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode-1)*xy.real() + - bz_arr(lo.x+j_bz+ix, lo.y+l_bz+iz, 0, 2*imode)*xy.imag()); + Bzp[ip] += sx_bz[ix]*sz_bz[iz]*dBz; } } // Gather field on particle Ezp[i] from field on grid ez_arr // Gather field on particle Bxp[i] from field on grid bx_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order; ix++){ - const amrex::Real dEz = (+ ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode-1)*xy.real() - - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.imag()); - Ezp[ip] += sx[ix]*sz0[iz]*dEz; - const amrex::Real dBx = (+ bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode-1)*xy.real() - - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.imag()); - Bxp[ip] += sx[ix]*sz0[iz]*dBx; + const amrex::Real dEz = (+ ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode-1)*xy.real() + - ez_arr(lo.x+j_ez+ix, lo.y+l_ez+iz, 0, 2*imode)*xy.imag()); + Ezp[ip] += sx_ez[ix]*sz_ez[iz]*dEz; + const amrex::Real dBx = (+ bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode-1)*xy.real() + - bx_arr(lo.x+j_bx+ix, lo.y+l_bx+iz, 0, 2*imode)*xy.imag()); + Bxp[ip] += sx_bx[ix]*sz_bx[iz]*dBx; } } // Gather field on particle Byp[i] from field on grid by_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - const amrex::Real dBy = (+ by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode-1)*xy.real() - - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode)*xy.imag()); - Byp[ip] += sx0[ix]*sz0[iz]*dBy; + const amrex::Real dBy = (+ by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode-1)*xy.real() + - by_arr(lo.x+j_by+ix, lo.y+l_by+iz, 0, 2*imode)*xy.imag()); + Byp[ip] += sx_by[ix]*sz_by[iz]*dBy; } } xy = xy*xy0; @@ -203,8 +305,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int iy=0; iy<=depos_order; iy++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - Exp[ip] += sx0[ix]*sy[iy]*sz[iz]* - ex_arr(lo.x+j0+ix, lo.y+k+iy, lo.z+l+iz); + Exp[ip] += sx_ex[ix]*sy_ex[iy]*sz_ex[iz]* + ex_arr(lo.x+j_ex+ix, lo.y+k_ex+iy, lo.z+l_ex+iz); } } } @@ -212,8 +314,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int iy=0; iy<=depos_order-lower_in_v; iy++){ for (int ix=0; ix<=depos_order; ix++){ - Eyp[ip] += sx[ix]*sy0[iy]*sz[iz]* - ey_arr(lo.x+j+ix, lo.y+k0+iy, lo.z+l+iz); + Eyp[ip] += sx_ey[ix]*sy_ey[iy]*sz_ey[iz]* + ey_arr(lo.x+j_ey+ix, lo.y+k_ey+iy, lo.z+l_ey+iz); } } } @@ -221,8 +323,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int iy=0; iy<=depos_order; iy++){ for (int ix=0; ix<=depos_order; ix++){ - Ezp[ip] += sx[ix]*sy[iy]*sz0[iz]* - ez_arr(lo.x+j+ix, lo.y+k+iy, lo.z+l0+iz); + Ezp[ip] += sx_ez[ix]*sy_ez[iy]*sz_ez[iz]* + ez_arr(lo.x+j_ez+ix, lo.y+k_ez+iy, lo.z+l_ez+iz); } } } @@ -230,8 +332,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int iy=0; iy<=depos_order-lower_in_v; iy++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - Bzp[ip] += sx0[ix]*sy0[iy]*sz[iz]* - bz_arr(lo.x+j0+ix, lo.y+k0+iy, lo.z+l+iz); + Bzp[ip] += sx_bz[ix]*sy_bz[iy]*sz_bz[iz]* + bz_arr(lo.x+j_bz+ix, lo.y+k_bz+iy, lo.z+l_bz+iz); } } } @@ -239,8 +341,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int iy=0; iy<=depos_order; iy++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - Byp[ip] += sx0[ix]*sy[iy]*sz0[iz]* - by_arr(lo.x+j0+ix, lo.y+k+iy, lo.z+l0+iz); + Byp[ip] += sx_by[ix]*sy_by[iy]*sz_by[iz]* + by_arr(lo.x+j_by+ix, lo.y+k_by+iy, lo.z+l_by+iz); } } } @@ -248,8 +350,8 @@ void doGatherShapeN(const amrex::ParticleReal * const xp, for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int iy=0; iy<=depos_order-lower_in_v; iy++){ for (int ix=0; ix<=depos_order; ix++){ - Bxp[ip] += sx[ix]*sy0[iy]*sz0[iz]* - bx_arr(lo.x+j+ix, lo.y+k0+iy, lo.z+l0+iz); + Bxp[ip] += sx_bx[ix]*sy_bx[iy]*sz_bx[iz]* + bx_arr(lo.x+j_bx+ix, lo.y+k_bx+iy, lo.z+l_bx+iz); } } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3128868e4..95e8ce4fe 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2169,8 +2169,6 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); - // Set staggering shift depending on e_is_nodal - const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. @@ -2185,13 +2183,6 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, // Add guard cells to the box. box.grow(ngE); - const Array4& ex_arr = exfab->array(); - const Array4& ey_arr = eyfab->array(); - const Array4& ez_arr = ezfab->array(); - const Array4& bx_arr = bxfab->array(); - const Array4& by_arr = byfab->array(); - const Array4& bz_arr = bzfab->array(); - const ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; @@ -2209,25 +2200,25 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, np_to_gather, dx, - xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); + xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doGatherShapeN<2,1>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, np_to_gather, dx, - xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); + xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doGatherShapeN<3,1>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, np_to_gather, dx, - xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); + xyzmin, lo, WarpX::n_rz_azimuthal_modes); } } else { if (WarpX::nox == 1){ @@ -2235,25 +2226,25 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, np_to_gather, dx, - xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); + xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doGatherShapeN<2,0>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, np_to_gather, dx, - xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); + xyzmin, lo, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doGatherShapeN<3,0>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, Ezp.dataPtr() + offset, Bxp.dataPtr() + offset, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, np_to_gather, dx, - xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); + xyzmin, lo, WarpX::n_rz_azimuthal_modes); } } } -- cgit v1.2.3 From b2aff9c35a19b3dd0fed550242cc444c45c5d5d3 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 11 Dec 2019 10:42:02 +0100 Subject: fixed local variable masking parameter --- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3128868e4..e7e6b0695 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1078,11 +1078,11 @@ PhysicalParticleContainer::Evolve (int lev, for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const auto np = pti.numParticles(); - const auto lev = pti.GetLevel(); + const auto t_lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); tmp_particle_data.resize(finestLevel()+1); for (int i = 0; i < TmpIdx::nattribs; ++i) - tmp_particle_data[lev][index][i].resize(np); + tmp_particle_data[t_lev][index][i].resize(np); } } -- cgit v1.2.3 From 49ebe085e0e05fab0df3c0a99498a302df2f9ef4 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 15 Dec 2019 09:11:01 -0800 Subject: Doxygen: add amrex:: prefix for function arguments in cpp files --- Source/BoundaryConditions/PML.cpp | 4 ++-- Source/Evolve/WarpXEvolveEM.cpp | 4 ++-- .../SpectralSolver/SpectralFieldData.cpp | 8 ++++---- Source/FieldSolver/WarpXPushFieldsEM.cpp | 14 ++++++------- Source/Filter/Filter.cpp | 6 +++--- Source/Particles/PhysicalParticleContainer.cpp | 20 ++++++++++-------- Source/Particles/WarpXParticleContainer.cpp | 8 ++++---- Source/WarpX.cpp | 24 +++++++++++----------- 8 files changed, 46 insertions(+), 42 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index f6d5e05ae..51439430d 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -752,14 +752,14 @@ PML::CopyJtoPMLs (const std::array& j_fp, void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) +PML::ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain) { ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) +PML::ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index fe0fb9157..470d12610 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -475,7 +475,7 @@ WarpX::OneStep_sub1 (Real curtime) } void -WarpX::PushParticlesandDepose (Real cur_time) +WarpX::PushParticlesandDepose (amrex::Real cur_time) { // Evolve particles to p^{n+1/2} and x^{n+1} // Depose current, j^{n+1/2} @@ -485,7 +485,7 @@ WarpX::PushParticlesandDepose (Real cur_time) } void -WarpX::PushParticlesandDepose (int lev, Real cur_time, DtType a_dt_type) +WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type) { mypc->Evolve(lev, *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 8f0853484..edd4df34d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -3,10 +3,10 @@ using namespace amrex; /* \brief Initialize fields in spectral space, and FFT plans */ -SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, - const SpectralKSpace& k_space, - const DistributionMapping& dm, - const int n_field_required ) +SpectralFieldData::SpectralFieldData( const amrex::BoxArray& realspace_ba, + const SpectralKSpace& k_space, + const amrex::DistributionMapping& dm, + const int n_field_required ) { const BoxArray& spectralspace_ba = k_space.spectralspace_ba; diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 9807665c6..4848b051e 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -90,7 +90,7 @@ WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) #endif void -WarpX::EvolveB (Real a_dt) +WarpX::EvolveB (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { EvolveB(lev, a_dt); @@ -98,7 +98,7 @@ WarpX::EvolveB (Real a_dt) } void -WarpX::EvolveB (int lev, Real a_dt) +WarpX::EvolveB (int lev, amrex::Real a_dt) { BL_PROFILE("WarpX::EvolveB()"); EvolveB(lev, PatchType::fine, a_dt); @@ -303,7 +303,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) } void -WarpX::EvolveE (Real a_dt) +WarpX::EvolveE (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -312,7 +312,7 @@ WarpX::EvolveE (Real a_dt) } void -WarpX::EvolveE (int lev, Real a_dt) +WarpX::EvolveE (int lev, amrex::Real a_dt) { BL_PROFILE("WarpX::EvolveE()"); EvolveE(lev, PatchType::fine, a_dt); @@ -611,7 +611,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) } void -WarpX::EvolveF (Real a_dt, DtType a_dt_type) +WarpX::EvolveF (amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -622,7 +622,7 @@ WarpX::EvolveF (Real a_dt, DtType a_dt_type) } void -WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) +WarpX::EvolveF (int lev, amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -631,7 +631,7 @@ WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) } void -WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) +WarpX::EvolveF (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 5d3c14c14..93729a0ae 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -144,7 +144,7 @@ void Filter::DoFilter (const Box& tbx, * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil()"); ncomp = std::min(ncomp, srcmf.nComp()); @@ -179,8 +179,8 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, - const Box& tbx, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (amrex::FArrayBox& dstfab, const amrex::FArrayBox& srcfab, + const amrex::Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.nComp()); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c77ddf83b..d5c08074a 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -963,8 +963,12 @@ PhysicalParticleContainer::EvolveES (const Vector& dx = WarpX::CellSize(lev); @@ -2148,12 +2152,12 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, RealVector& Bxp, RealVector& Byp, RealVector& Bzp, - FArrayBox const * exfab, - FArrayBox const * eyfab, - FArrayBox const * ezfab, - FArrayBox const * bxfab, - FArrayBox const * byfab, - FArrayBox const * bzfab, + amrex::FArrayBox const * exfab, + amrex::FArrayBox const * eyfab, + amrex::FArrayBox const * ezfab, + amrex::FArrayBox const * bxfab, + amrex::FArrayBox const * byfab, + amrex::FArrayBox const * bzfab, const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 8926301b2..e6c6d31bd 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -420,7 +420,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, const int * const ion_lev, - MultiFab* rho, int icomp, + amrex::MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) { @@ -507,7 +507,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, } void -WarpXParticleContainer::DepositCharge (Vector >& rho, +WarpXParticleContainer::DepositCharge (amrex::Vector >& rho, bool local, bool reset, bool do_rz_volume_scaling) { @@ -760,7 +760,7 @@ WarpXParticleContainer::PushXES (Real dt) } void -WarpXParticleContainer::PushX (Real dt) +WarpXParticleContainer::PushX (amrex::Real dt) { const int nLevels = finestLevel(); for (int lev = 0; lev <= nLevels; ++lev) { @@ -769,7 +769,7 @@ WarpXParticleContainer::PushX (Real dt) } void -WarpXParticleContainer::PushX (int lev, Real dt) +WarpXParticleContainer::PushX (int lev, amrex::Real dt) { BL_PROFILE("WPC::PushX()"); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d94e35f61..53728b54f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1116,9 +1116,9 @@ WarpX::Evolve (int numsteps) { } void -WarpX::ComputeDivB (MultiFab& divB, int dcomp, - const std::array& B, - const std::array& dx) +WarpX::ComputeDivB (amrex::MultiFab& divB, int dcomp, + const std::array& B, + const std::array& dx) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1150,9 +1150,9 @@ WarpX::ComputeDivB (MultiFab& divB, int dcomp, } void -WarpX::ComputeDivB (MultiFab& divB, int dcomp, - const std::array& B, - const std::array& dx, int ngrow) +WarpX::ComputeDivB (amrex::MultiFab& divB, int dcomp, + const std::array& B, + const std::array& dx, int ngrow) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1184,9 +1184,9 @@ WarpX::ComputeDivB (MultiFab& divB, int dcomp, } void -WarpX::ComputeDivE (MultiFab& divE, int dcomp, - const std::array& E, - const std::array& dx) +WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp, + const std::array& E, + const std::array& dx) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1218,9 +1218,9 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, } void -WarpX::ComputeDivE (MultiFab& divE, int dcomp, - const std::array& E, - const std::array& dx, int ngrow) +WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp, + const std::array& E, + const std::array& dx, int ngrow) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; -- cgit v1.2.3