From 20c785a67312b88e497a94e1a27a343918cdc9ce Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 8 Jul 2019 10:24:40 -0700 Subject: Cleaning code for slices in boosted frame diag --- 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 7e7c9534e..27155e1c6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1897,7 +1897,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real for (int lev = 0; lev < nlevs; ++lev) { - const Real* dx = Geom(lev).CellSize(); + const Real* dx = Geom(lev).CellSize(); const Real* plo = Geom(lev).ProbLo(); // first we touch each map entry in serial @@ -1951,7 +1951,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real // timestep, skip it. if ( not (((zp_new[i] >= z_new) && (zp_old[i] <= z_old)) || ((zp_new[i] <= z_new) && (zp_old[i] >= z_old))) ) continue; - + // Lorentz transform particles to lab frame Real gamma_new_p = std::sqrt(1.0 + inv_c2*(uxp_new[i]*uxp_new[i] + uyp_new[i]*uyp_new[i] + uzp_new[i]*uzp_new[i])); Real t_new_p = WarpX::gamma_boost*t_boost - uzfrm*zp_new[i]*inv_c2; -- cgit v1.2.3 From 045f22b541752b9940aae91b5b6031870c13e7c3 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 17 Sep 2019 11:45:40 +0200 Subject: Enabled selection of Boris+RR, bugfixing --- Source/Particles/PhysicalParticleContainer.cpp | 14 ++++++++++++++ .../Pusher/UpdateMomentumBorisWithRadiationReaction.H | 6 +++--- Source/Particles/RigidInjectedParticleContainer.cpp | 12 ++++++++++++ Source/Utils/WarpXAlgorithmSelection.H | 3 ++- Source/Utils/WarpXAlgorithmSelection.cpp | 1 + 5 files changed, 32 insertions(+), 4 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4bc0ee16e..95f2243b3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -15,6 +15,7 @@ #include #include #include +#include using namespace amrex; @@ -1584,6 +1585,19 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::BorisRR) { + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } + UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); + UpdatePosition( x[i], y[i], z[i], + ux[i], uy[i], uz[i], dt ); + } + ); } else { amrex::Abort("Unknown particle pusher"); }; diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H index dd7c41a94..dbb5c02da 100644 --- a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H +++ b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H @@ -41,9 +41,9 @@ void UpdateMomentumBorisWithRadiationReaction( const amrex::Real inv_gamma_n = 1.0/gamma_n; //Estimation of the velocity at intermediate (integer) time - const amrex::Real vx_n = ux_n*inv_gamma; - const amrex::Real vy_n = uy_n*inv_gamma; - const amrex::Real vz_n = uz_n*inv_gamma; + const amrex::Real vx_n = ux_n*inv_gamma_n; + const amrex::Real vy_n = uy_n*inv_gamma_n; + const amrex::Real vz_n = uz_n*inv_gamma_n; //Lorentz force const amrex::Real flx = q*(Ex + vy_n*Bz - vz_n*By); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index cd5e34770..fd002a9af 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -13,6 +13,7 @@ #include #include #include +#include using namespace amrex; @@ -443,6 +444,17 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::BorisRR) { + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumBorisWithRadiationReaction( + uxpp[i], uypp[i], uzpp[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); + } + ); } else { amrex::Abort("Unknown particle pusher"); }; diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 54c721abf..35a15c6fd 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -14,7 +14,8 @@ struct MaxwellSolverAlgo { struct ParticlePusherAlgo { enum { Boris = 0, - Vay = 1 + Vay = 1, + BorisRR = 2 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 8a6ff6dbf..45f94f52f 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -18,6 +18,7 @@ const std::map maxwell_solver_algo_to_int = { const std::map particle_pusher_algo_to_int = { {"boris", ParticlePusherAlgo::Boris }, {"vay", ParticlePusherAlgo::Vay }, + {"borisRR", ParticlePusherAlgo::BorisRR }, {"default", ParticlePusherAlgo::Boris } }; -- cgit v1.2.3 From 7dbb4c297420bb9360ec6667e081b8e56c7efcd3 Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Fri, 20 Sep 2019 08:53:22 -0700 Subject: fixing EOL space --- Source/Diagnostics/BoostedFrameDiagnostic.H | 80 ++++----- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 224 ++++++++++++------------- Source/Initialization/WarpXInitData.cpp | 6 +- Source/Particles/PhysicalParticleContainer.cpp | 4 +- 4 files changed, 157 insertions(+), 157 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index e2e2478aa..397c53c0e 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -22,17 +22,17 @@ /// /** \brief - * The capability for full 3D snapshot and slices in the lab-frame is implemented. - * LabFrameDiag class defines the parameters required to copy data from boost-frame at (z_boost t_boost) + * The capability for full 3D snapshot and slices in the lab-frame is implemented. + * LabFrameDiag class defines the parameters required to copy data from boost-frame at (z_boost t_boost) * to lab-frame at (z_lab, t_lab) - * Two derived classes, namely, LabFrameSnapShot and LabFrameSlice are defined to allow for + * Two derived classes, namely, LabFrameSnapShot and LabFrameSlice are defined to allow for * different functions required to copy boosted frame data to lab-frame slice. - * The approach here is to define an array of LabFrameDiag which could include both, snapshots and slices, - * sorted based on the t_lab of each of the diags. This is done to re-use the z_lab data for the diags - * that have the same t_lab, instead of re-generating the slice at z_lab for that t_lab. + * The approach here is to define an array of LabFrameDiag which could include both, snapshots and slices, + * sorted based on the t_lab of each of the diags. This is done to re-use the z_lab data for the diags + * that have the same t_lab, instead of re-generating the slice at z_lab for that t_lab. */ - -class LabFrameDiag { + +class LabFrameDiag { public: std::string file_name; amrex::Real t_lab; @@ -40,7 +40,7 @@ class LabFrameDiag { amrex::IntVect prob_ncells_lab_; amrex::RealBox diag_domain_lab_; amrex::Box buff_box_; - + amrex::Real current_z_lab; amrex::Real current_z_boost; amrex::Real inv_gamma_boost_; @@ -48,20 +48,20 @@ class LabFrameDiag { amrex::Real dz_lab_; amrex::Real dx_; amrex::Real dy_; - + int ncomp_to_dump_; std::vector mesh_field_names_; - + int file_num; int initial_i; - + // For back-transformed diagnostics of grid fields, data_buffer_ // stores a buffer of the fields in the lab frame (in a MultiFab, i.e. // with all box data etc.). When the buffer if full, dump to file. std::unique_ptr data_buffer_; // particles_buffer_ is currently blind to refinement level. // particles_buffer_[j] is a WarpXParticleContainer::DiagnosticParticleData where - // - j is the species number for the current diag + // - j is the species number for the current diag amrex::Vector particles_buffer_; // buff_counter_ is the number of z slices in data_buffer_ int buff_counter_; @@ -70,66 +70,66 @@ class LabFrameDiag { /// /// This snapshot is at time t_lab, and the simulation is at time t_boost. /// The Lorentz transformation picks out one slice corresponding to both - /// of those times, at position current_z_boost and current_z_lab in the - /// boosted and lab frames, respectively. + /// of those times, at position current_z_boost and current_z_lab in the + /// boosted and lab frames, respectively. /// void updateCurrentZPositions(amrex::Real t_boost, amrex::Real inv_gamma, amrex::Real inv_beta); - + void createLabFrameDirectories(); /// /// Write some useful metadata about this snapshot. /// void writeLabFrameHeader(); - + virtual int getId() {return 0;} virtual void AddDataToBuffer(amrex::MultiFab& tmp_slice_ptr, int i_lab, - amrex::Gpu::ManagedDeviceVector map_actual_fields_to_dump){}; + amrex::Gpu::ManagedDeviceVector map_actual_fields_to_dump){}; virtual void AddPartDataToParticleBuffer( amrex::Vector tmp_particle_buffer, int nSpeciesBoostedFrame) {}; }; /// LabFrameSnapShot stores metadata corresponding to a single time -/// snapshot in the lab frame. The snapshot is written to disk +/// snapshot in the lab frame. The snapshot is written to disk /// in lab_frame_data/snapshots. zmin_lab, zmax_lab, and t_lab -/// are all constant for a given snapshot. current_z_lab and -/// current_z_boost for each snapshot are updated as the +/// are all constant for a given snapshot. current_z_lab and +/// current_z_boost for each snapshot are updated as the /// simulation time in the boosted frame advances. /// class LabFrameSnapShot : public LabFrameDiag { public: int snapshot_id; int getId() {return snapshot_id;} - LabFrameSnapShot(amrex::Real t_lab_in, amrex::Real t_boost, - amrex::Real inv_gamma_boost_, amrex::Real inv_beta_boost_, - amrex::Real dz_lab_, amrex::RealBox prob_domain_lab, - amrex::IntVect prob_ncells_lab, int ncomp_to_dump, - std::vector mesh_field_names, - amrex::RealBox diag_domain_lab, + LabFrameSnapShot(amrex::Real t_lab_in, amrex::Real t_boost, + amrex::Real inv_gamma_boost_, amrex::Real inv_beta_boost_, + amrex::Real dz_lab_, amrex::RealBox prob_domain_lab, + amrex::IntVect prob_ncells_lab, int ncomp_to_dump, + std::vector mesh_field_names, + amrex::RealBox diag_domain_lab, amrex::Box diag_box, int file_num_in); void AddDataToBuffer( amrex::MultiFab& tmp_slice, int k_lab, - amrex::Gpu::ManagedDeviceVector map_actual_fields_to_dump); + amrex::Gpu::ManagedDeviceVector map_actual_fields_to_dump); void AddPartDataToParticleBuffer( amrex::Vector tmp_particle_buffer, int nSpeciesBoostedFrame); }; /** \brief - * LabFrameSlice stores the data corresponding to a single time at the user-defined slice location + * LabFrameSlice stores the data corresponding to a single time at the user-defined slice location * The slice is written to disk in the lab_frame_data/slices - * Similar to snapshots, zmin_lab, zmax_lab, and t_lab are constant for a given slice. + * Similar to snapshots, zmin_lab, zmax_lab, and t_lab are constant for a given slice. * current_z_lab and current_z_boost for each snapshot are updated as the sim time in boosted frame advances. */ class LabFrameSlice : public LabFrameDiag { public: int slice_id; int getId() {return slice_id;} - LabFrameSlice(amrex::Real t_lab_in, amrex::Real t_boost, - amrex::Real inv_gamma_boost_, amrex::Real inv_beta_boost_, - amrex::Real dz_lab_, amrex::RealBox prob_domain_lab, - amrex::IntVect prob_ncells_lab, int ncomp_to_dump, - std::vector mesh_field_names, - amrex::RealBox diag_domain_lab, - amrex::Box diag_box, int file_num_in, + LabFrameSlice(amrex::Real t_lab_in, amrex::Real t_boost, + amrex::Real inv_gamma_boost_, amrex::Real inv_beta_boost_, + amrex::Real dz_lab_, amrex::RealBox prob_domain_lab, + amrex::IntVect prob_ncells_lab, int ncomp_to_dump, + std::vector mesh_field_names, + amrex::RealBox diag_domain_lab, + amrex::Box diag_box, int file_num_in, amrex::Real cell_dx, amrex::Real cell_dy); void AddDataToBuffer( amrex::MultiFab& tmp_slice_ptr, int i_lab, amrex::Gpu::ManagedDeviceVector map_actual_fields_to_dump); @@ -155,7 +155,7 @@ class BoostedFrameDiagnostic { int max_box_size_ = 256; std::vector > LabFrameDiags_; - + void writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const int i_lab); @@ -167,11 +167,11 @@ public: BoostedFrameDiagnostic(amrex::Real zmin_lab, amrex::Real zmax_lab, amrex::Real v_window_lab, amrex::Real dt_snapshots_lab, - int N_snapshots, amrex::Real dt_slice_snapshots_lab, + int N_snapshots, amrex::Real dt_slice_snapshots_lab, int N_slice_snapshots, amrex::Real gamma_boost, amrex::Real t_boost, amrex::Real dt_boost, int boost_direction, const amrex::Geometry& geom, amrex::RealBox& slice_realbox); - + void Flush(const amrex::Geometry& geom); void writeLabFrameData(const amrex::MultiFab* cell_centered_data, diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 92dcad417..db40f8682 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -447,7 +447,7 @@ namespace } #endif -bool compare_tlab_uptr(const std::unique_ptr&a, +bool compare_tlab_uptr(const std::unique_ptr&a, const std::unique_ptr&b) { return a->t_lab < b->t_lab; @@ -531,10 +531,10 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) BoostedFrameDiagnostic:: BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, - Real dt_snapshots_lab, int N_snapshots, + Real dt_snapshots_lab, int N_snapshots, Real dt_slice_snapshots_lab, int N_slice_snapshots, - Real gamma_boost, Real t_boost, Real dt_boost, - int boost_direction, const Geometry& geom, + Real gamma_boost, Real t_boost, Real dt_boost, + int boost_direction, const Geometry& geom, amrex::RealBox& slice_realbox) : gamma_boost_(gamma_boost), dt_snapshots_lab_(dt_snapshots_lab), @@ -601,35 +601,35 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, // x and y bounds are the same for lab frame and boosted frame prob_domain_lab.setLo(AMREX_SPACEDIM-1, zmin_lab + v_window_lab * t_lab); prob_domain_lab.setHi(AMREX_SPACEDIM-1, zmax_lab + v_window_lab * t_lab); - Box diag_box = geom.Domain(); - LabFrameDiags_[i].reset(new LabFrameSnapShot(t_lab, t_boost, + Box diag_box = geom.Domain(); + LabFrameDiags_[i].reset(new LabFrameSnapShot(t_lab, t_boost, inv_gamma_boost_, inv_beta_boost_, dz_lab_, - prob_domain_lab, prob_ncells_lab, + prob_domain_lab, prob_ncells_lab, ncomp_to_dump, mesh_field_names, prob_domain_lab, diag_box, i)); } - + amrex::Real zmin_slice_lab; amrex::Real zmax_slice_lab; amrex::Real cell_dx = 0; amrex::Real cell_dy = 0; IntVect slice_ncells_lab ; Box slicediag_box; - + if ( N_slice_snapshots_ > 0) { const amrex::Real* current_slice_lo = slice_realbox.lo(); - const amrex::Real* current_slice_hi = slice_realbox.hi(); + const amrex::Real* current_slice_hi = slice_realbox.hi(); zmin_slice_lab = current_slice_lo[AMREX_SPACEDIM-1] / ( (1.+beta_boost_)*gamma_boost_); zmax_slice_lab = current_slice_hi[AMREX_SPACEDIM-1] / - ( (1.+beta_boost_)*gamma_boost_); + ( (1.+beta_boost_)*gamma_boost_); int Nz_slice_lab = static_cast((zmax_slice_lab - zmin_slice_lab) * inv_dz_lab_); int Nx_slice_lab = ( current_slice_hi[0] - current_slice_lo[0] ) / geom.CellSize(0); if (Nx_slice_lab == 0 ) {Nx_slice_lab = 1;} - cell_dx = geom.CellSize(0); + cell_dx = geom.CellSize(0); #if (AMREX_SPACEDIM == 3) - int Ny_slice_lab = ( current_slice_hi[1] - current_slice_lo[1]) / + int Ny_slice_lab = ( current_slice_hi[1] - current_slice_lo[1]) / geom.CellSize(1); if (Ny_slice_lab == 0 ) {Ny_slice_lab = 1;} slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab}; @@ -647,13 +647,13 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, slice_lo[i_dim] = (slice_realbox.lo(i_dim) - geom.ProbLo(i_dim) - 0.5*geom.CellSize(i_dim))/geom.CellSize(i_dim); slice_hi[i_dim] = (slice_realbox.hi(i_dim) - geom.ProbLo(i_dim) - 0.5*geom.CellSize(i_dim))/geom.CellSize(i_dim); } - Box stmp(slice_lo,slice_hi); + Box stmp(slice_lo,slice_hi); slicediag_box = stmp; } for (int i = 0; i < N_slice_snapshots; ++i) { Real t_slice_lab = i * dt_slice_snapshots_lab_ ; RealBox prob_domain_lab = geom.ProbDomain(); - // replace z bounds by lab-frame coordinates + // replace z bounds by lab-frame coordinates prob_domain_lab.setLo(AMREX_SPACEDIM-1, zmin_lab + v_window_lab * t_slice_lab); prob_domain_lab.setHi(AMREX_SPACEDIM-1, zmax_lab + v_window_lab * t_slice_lab); RealBox slice_dom_lab = slice_realbox; @@ -661,17 +661,17 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, // note : x and y bounds are the same for lab and boosted frames // initial lab slice extent // slice_dom_lab.setLo(AMREX_SPACEDIM-1, zmin_slice_lab + v_window_lab * t_slice_lab ); - slice_dom_lab.setHi(AMREX_SPACEDIM-1, zmax_slice_lab + + slice_dom_lab.setHi(AMREX_SPACEDIM-1, zmax_slice_lab + v_window_lab * t_slice_lab ); - + // construct labframeslice - LabFrameDiags_[i+N_snapshots].reset(new LabFrameSlice(t_slice_lab, t_boost, + LabFrameDiags_[i+N_snapshots].reset(new LabFrameSlice(t_slice_lab, t_boost, inv_gamma_boost_, inv_beta_boost_, dz_lab_, - prob_domain_lab, slice_ncells_lab, + prob_domain_lab, slice_ncells_lab, ncomp_to_dump, mesh_field_names, slice_dom_lab, slicediag_box, i, cell_dx, cell_dy)); } - // sort diags based on their respective t_lab + // sort diags based on their respective t_lab std::stable_sort(LabFrameDiags_.begin(), LabFrameDiags_.end(), compare_tlab_uptr); AMREX_ALWAYS_ASSERT(max_box_size_ >= num_buffer_); @@ -698,7 +698,7 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) const BoxArray& ba = LabFrameDiags_[i]->data_buffer_->boxArray(); const int hi = ba[0].bigEnd(boost_direction_); const int lo = hi - LabFrameDiags_[i]->buff_counter_ + 1; - + //Box buff_box = geom.Domain(); Box buff_box = LabFrameDiags_[i]->buff_box_; buff_box.setSmall(boost_direction_, lo); @@ -707,17 +707,17 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) BoxArray buff_ba(buff_box); buff_ba.maxSize(max_box_size_); DistributionMapping buff_dm(buff_ba); - + const int ncomp = LabFrameDiags_[i]->data_buffer_->nComp(); - + MultiFab tmp(buff_ba, buff_dm, ncomp, 0); - + tmp.copy(*LabFrameDiags_[i]->data_buffer_, 0, 0, ncomp); #ifdef WARPX_USE_HDF5 for (int comp = 0; comp < ncomp; ++comp) output_write_field(LabFrameDiags_[i]->file_name, mesh_field_names[comp], tmp, comp); -#else +#else std::stringstream ss; ss << LabFrameDiags_[i]->file_name << "/Level_0/" << Concatenate("buffer", i_lab, 5); VisMF::Write(tmp, ss.str()); @@ -780,7 +780,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, Real diag_zmin_lab = LabFrameDiags_[i]->diag_domain_lab_.lo(AMREX_SPACEDIM-1); Real diag_zmax_lab = LabFrameDiags_[i]->diag_domain_lab_.hi(AMREX_SPACEDIM-1); - + if ( ( LabFrameDiags_[i]->current_z_boost < zlo_boost) or ( LabFrameDiags_[i]->current_z_boost > zhi_boost) or ( LabFrameDiags_[i]->current_z_lab < diag_zmin_lab) or @@ -806,7 +806,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, buff_dm, ncomp_to_dump, 0) ); } // ... reset particle buffer particles_buffer_[i] - if (WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_particles) LabFrameDiags_[i]->particles_buffer_.resize(mypc.nSpeciesBoostedFrameDiags()); } @@ -815,9 +815,9 @@ writeLabFrameData(const MultiFab* cell_centered_data, const int start_comp = 0; const bool interpolate = true; // tmp_slice_ptr containing slice data is generated only if t_lab != prev_t_lab - if (LabFrameDiags_[i]->t_lab != prev_t_lab ) { + if (LabFrameDiags_[i]->t_lab != prev_t_lab ) { // Get slice in the boosted frame - if (tmp_slice_ptr) + if (tmp_slice_ptr) { tmp_slice_ptr.reset(new MultiFab()); tmp_slice_ptr.reset(nullptr); @@ -826,7 +826,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, LabFrameDiags_[i]->current_z_boost, *cell_centered_data, geom, start_comp, ncomp, interpolate); - + // transform it to the lab frame LorentzTransformZ(*slice, gamma_boost_, beta_boost_, ncomp); // Create a 2D box for the slice in the boosted frame @@ -838,46 +838,46 @@ writeLabFrameData(const MultiFab* cell_centered_data, // Make it a BoxArray slice_ba BoxArray slice_ba(slice_box); slice_ba.maxSize(max_box_size_); - // Create MultiFab tmp on slice_ba witih data from slice, that can be potentially re-used. - tmp_slice_ptr = std::unique_ptr(new MultiFab(slice_ba, - LabFrameDiags_[i]->data_buffer_->DistributionMap(), + // Create MultiFab tmp on slice_ba witih data from slice, that can be potentially re-used. + tmp_slice_ptr = std::unique_ptr(new MultiFab(slice_ba, + LabFrameDiags_[i]->data_buffer_->DistributionMap(), ncomp, 0)); - + tmp_slice_ptr->copy(*slice, 0, 0, ncomp); } // tmp_slice_ptr is re-used if the t_lab of a diag is equal to that of the previous diag LabFrameDiags_[i]->AddDataToBuffer(*tmp_slice_ptr, i_lab, map_actual_fields_to_dump); - //// if t_lab of the next diag is different from current t_lab, - //// deallocate and nullify tmp_slice_ptr. - //if ( (i==LabFrameDiags_.size()-1) || - // ((i+1)t_lab != prev_t_lab ) { + //// if t_lab of the next diag is different from current t_lab, + //// deallocate and nullify tmp_slice_ptr. + //if ( (i==LabFrameDiags_.size()-1) || + // ((i+1)t_lab != prev_t_lab ) { // tmp_slice_ptr.reset(new MultiFab()); // tmp_slice_ptr.reset(nullptr); //} } - if (WarpX::do_boosted_frame_particles) { - - if (LabFrameDiags_[i]->t_lab != prev_t_lab ) { - if (tmp_particle_buffer.size()>0) + if (WarpX::do_boosted_frame_particles) { + + if (LabFrameDiags_[i]->t_lab != prev_t_lab ) { + if (tmp_particle_buffer.size()>0) { tmp_particle_buffer.clear(); tmp_particle_buffer.shrink_to_fit(); } - tmp_particle_buffer.resize(mypc.nSpeciesBoostedFrameDiags()); + tmp_particle_buffer.resize(mypc.nSpeciesBoostedFrameDiags()); mypc.GetLabFrameData( LabFrameDiags_[i]->file_name, i_lab, boost_direction_, old_z_boost, LabFrameDiags_[i]->current_z_boost, t_boost, LabFrameDiags_[i]->t_lab, dt, tmp_particle_buffer); - } - LabFrameDiags_[i]->AddPartDataToParticleBuffer(tmp_particle_buffer, - mypc.nSpeciesBoostedFrameDiags()); + } + LabFrameDiags_[i]->AddPartDataToParticleBuffer(tmp_particle_buffer, + mypc.nSpeciesBoostedFrameDiags()); } ++LabFrameDiags_[i]->buff_counter_; prev_t_lab = LabFrameDiags_[i]->t_lab; - + // If buffer full, write to disk. if ( LabFrameDiags_[i]->buff_counter_ == num_buffer_) { @@ -911,7 +911,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // Write data to disk (custom) writeParticleData(LabFrameDiags_[i]->particles_buffer_[j], part_ss.str(), i_lab); #endif - } + } LabFrameDiags_[i]->particles_buffer_.clear(); } LabFrameDiags_[i]->buff_counter_ = 0; @@ -1057,16 +1057,16 @@ writeMetaData () std::ofstream HeaderFile_slice; HeaderFile_slice.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); std::string HeaderFileName_slice(WarpX::lab_data_directory+"/slices/Header"); - HeaderFile_slice.open(HeaderFileName_slice.c_str(), + HeaderFile_slice.open(HeaderFileName_slice.c_str(), std::ofstream::out | std::ofstream::trunc | std::ofstream::binary); if (!HeaderFile_slice.good()) FileOpenFailed(HeaderFileName_slice); - + HeaderFile_slice.precision(17); - + HeaderFile_slice << N_slice_snapshots_ << "\n"; HeaderFile_slice << dt_slice_snapshots_lab_ << "\n"; HeaderFile_slice << gamma_boost_ << "\n"; @@ -1076,14 +1076,14 @@ writeMetaData () } - + } LabFrameSnapShot:: -LabFrameSnapShot(Real t_lab_in, Real t_boost, Real inv_gamma_boost_, - Real inv_beta_boost_, Real dz_lab_, RealBox prob_domain_lab, - IntVect prob_ncells_lab, int ncomp_to_dump, - std::vector mesh_field_names, +LabFrameSnapShot(Real t_lab_in, Real t_boost, Real inv_gamma_boost_, + Real inv_beta_boost_, Real dz_lab_, RealBox prob_domain_lab, + IntVect prob_ncells_lab, int ncomp_to_dump, + std::vector mesh_field_names, amrex::RealBox diag_domain_lab, Box diag_box, int file_num_in) { t_lab = t_lab_in; @@ -1110,13 +1110,13 @@ LabFrameSnapShot(Real t_lab_in, Real t_boost, Real inv_gamma_boost_, void LabFrameDiag:: -updateCurrentZPositions(Real t_boost, Real inv_gamma, Real inv_beta) +updateCurrentZPositions(Real t_boost, Real inv_gamma, Real inv_beta) { current_z_boost = (t_lab*inv_gamma - t_boost)*PhysConst::c*inv_beta; current_z_lab = (t_lab - t_boost*inv_gamma)*PhysConst::c*inv_beta; } -void +void LabFrameDiag:: createLabFrameDirectories() { #ifdef WARPX_USE_HDF5 @@ -1197,7 +1197,7 @@ createLabFrameDirectories() { } -void +void LabFrameDiag:: writeLabFrameHeader() { #ifndef WARPX_USE_HDF5 @@ -1247,11 +1247,11 @@ writeLabFrameHeader() { LabFrameSlice:: -LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_, - Real inv_beta_boost_, Real dz_lab_, RealBox prob_domain_lab, - IntVect prob_ncells_lab, int ncomp_to_dump, - std::vector mesh_field_names, - RealBox diag_domain_lab, Box diag_box, int file_num_in, +LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_, + Real inv_beta_boost_, Real dz_lab_, RealBox prob_domain_lab, + IntVect prob_ncells_lab, int ncomp_to_dump, + std::vector mesh_field_names, + RealBox diag_domain_lab, Box diag_box, int file_num_in, amrex::Real cell_dx, amrex::Real cell_dy) { t_lab = t_lab_in; @@ -1272,13 +1272,13 @@ LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_, buff_counter_ = 0; dx_ = cell_dx; dy_ = cell_dy; - + if (WarpX::do_boosted_frame_fields) data_buffer_.reset(nullptr); } -void +void LabFrameSnapShot:: -AddDataToBuffer( MultiFab& tmp, int k_lab, +AddDataToBuffer( MultiFab& tmp, int k_lab, amrex::Gpu::ManagedDeviceVector map_actual_fields_to_dump) { const int ncomp_to_dump = map_actual_fields_to_dump.size(); @@ -1286,7 +1286,7 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) { Array4 tmp_arr = tmp[mfi].array(); Array4 buf_arr = buf[mfi].array(); - // For 3D runs, rmp is a 2D (x,y) multifab that contains only + // For 3D runs, rmp is a 2D (x,y) multifab that contains only // slice to write to file const Box& bx = mfi.tilebox(); const auto field_map_ptr = map_actual_fields_to_dump.dataPtr(); @@ -1305,22 +1305,22 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, } -void +void LabFrameSlice:: -AddDataToBuffer( MultiFab& tmp, int k_lab, +AddDataToBuffer( MultiFab& tmp, int k_lab, amrex::Gpu::ManagedDeviceVector map_actual_fields_to_dump) { const int ncomp_to_dump = map_actual_fields_to_dump.size(); MultiFab& buf = *data_buffer_; - - for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) + + for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) { Array4 tmp_arr = tmp[mfi].array(); Array4 buf_arr = buf[mfi].array(); Box& bx = buff_box_; const Box& bx_bf = mfi.tilebox(); bx.setSmall(AMREX_SPACEDIM-1,bx_bf.smallEnd(AMREX_SPACEDIM-1)); - bx.setBig(AMREX_SPACEDIM-1,bx_bf.bigEnd(AMREX_SPACEDIM-1)); + bx.setBig(AMREX_SPACEDIM-1,bx_bf.bigEnd(AMREX_SPACEDIM-1)); if (bx_bf.contains(bx)) { const auto field_map_ptr = map_actual_fields_to_dump.dataPtr(); ParallelFor(bx, ncomp_to_dump, @@ -1331,16 +1331,16 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp); #else buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp); -#endif +#endif } - ); + ); } } - + } -void +void LabFrameSnapShot:: AddPartDataToParticleBuffer(Vector tmp_particle_buffer, int nspeciesBoostedFrame) { for (int isp = 0; isp < nspeciesBoostedFrame; ++isp) { @@ -1348,42 +1348,42 @@ AddPartDataToParticleBuffer(Vector tmp_particle_buffer, int nSpeciesBoostedFrame) { for (int isp = 0; isp < nSpeciesBoostedFrame; ++isp) { auto np = tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size(); - + if (np == 0) return; - + auto const& wpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::w); auto const& xpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::x); auto const& ypc = tmp_particle_buffer[isp].GetRealData(DiagIdx::y); @@ -1409,7 +1409,7 @@ AddPartDataToParticleBuffer(Vector= (diag_domain_lab_.lo(0)-dx_) && xpc[i] <= (diag_domain_lab_.hi(0)+dx_) ) { #if (AMREX_SPACEDIM == 3) - if( ypc[i] >= (diag_domain_lab_.lo(1)-dy_) && ypc[i] <= (diag_domain_lab_.hi(1) + dy_)) + if( ypc[i] >= (diag_domain_lab_.lo(1)-dy_) && ypc[i] <= (diag_domain_lab_.hi(1) + dy_)) #endif { wpc_buff[partcounter] = wpc[i]; xpc_buff[partcounter] = xpc[i]; ypc_buff[partcounter] = ypc[i]; - zpc_buff[partcounter] = zpc[i]; + zpc_buff[partcounter] = zpc[i]; uxpc_buff[partcounter] = uxpc[i]; uypc_buff[partcounter] = uypc[i]; uzpc_buff[partcounter] = uzpc[i]; @@ -1439,8 +1439,8 @@ AddPartDataToParticleBuffer(Vector= z_new) && (zp_old[i] <= z_old)) || ((zp_new[i] <= z_new) && (zp_old[i] >= z_old))) ) continue; - + // Lorentz transform particles to lab frame Real gamma_new_p = std::sqrt(1.0 + inv_c2*(uxp_new[i]*uxp_new[i] + uyp_new[i]*uyp_new[i] + uzp_new[i]*uzp_new[i])); Real t_new_p = WarpX::gamma_boost*t_boost - uzfrm*zp_new[i]*inv_c2; -- cgit v1.2.3 From 389d96a9225dfeb03d351e6726b7b52c72476626 Mon Sep 17 00:00:00 2001 From: Luca Date: Tue, 24 Sep 2019 20:11:13 +0200 Subject: bugfix --- Examples/Tests/RadiationReaction/test_RR.3d | 1 + Source/Particles/PhysicalParticleContainer.cpp | 8 ++++++++ Source/Utils/WarpXAlgorithmSelection.cpp | 2 +- 3 files changed, 10 insertions(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Tests/RadiationReaction/test_RR.3d b/Examples/Tests/RadiationReaction/test_RR.3d index 40cd651b4..9bf25198e 100644 --- a/Examples/Tests/RadiationReaction/test_RR.3d +++ b/Examples/Tests/RadiationReaction/test_RR.3d @@ -25,6 +25,7 @@ warpx.B_external = 0.0 0.0 1.0e5 warpx.verbose = 1 # Algorithms +algo.particle_pusher = "boris_rr" # CFL warpx.cfl = 1.0 diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 8d0c988b4..3d727af84 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1567,6 +1567,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; + if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), @@ -1701,6 +1702,13 @@ PhysicalParticleContainer::PushP (int lev, Real dt, Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } ); + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::BorisRR) { + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + UpdateMomentumBorisWithRadiationReaction( 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"); }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 45f94f52f..8e65cb8d0 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -18,7 +18,7 @@ const std::map maxwell_solver_algo_to_int = { const std::map particle_pusher_algo_to_int = { {"boris", ParticlePusherAlgo::Boris }, {"vay", ParticlePusherAlgo::Vay }, - {"borisRR", ParticlePusherAlgo::BorisRR }, + {"boris_rr", ParticlePusherAlgo::BorisRR }, {"default", ParticlePusherAlgo::Boris } }; -- cgit v1.2.3 From 8f375c0c07f17b27115b61b95d21cf9bb8f554e7 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Mon, 30 Sep 2019 09:50:03 +0200 Subject: merged with dev --- Source/Particles/PhysicalParticleContainer.H | 6 +++++ Source/Particles/PhysicalParticleContainer.cpp | 31 +++++++++++++++++++++----- 2 files changed, 31 insertions(+), 6 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ace1ec7f8..d789c7a79 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -145,6 +145,12 @@ protected: // Inject particles during the whole simulation void ContinuousInjection (const amrex::RealBox& injection_box) override; +private: + + //This function return true if the PhysicalParticleContainer contains electrons + //or positrons, false otherwise + bool AmIALepton(); + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3d727af84..f3b1bdce7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1568,7 +1568,13 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real q = this->charge; const Real m = this-> mass; - if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + + // Boris + RR is enabled only for leptons + auto algo = WarpX::particle_pusher_algo; + if (!AmIALepton() && algo == ParticlePusherAlgo::BorisRR) + algo = ParticlePusherAlgo::Boris; + + if (algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { @@ -1581,7 +1587,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + } else if (algo == ParticlePusherAlgo::Vay) { amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { @@ -1594,7 +1600,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, ux[i], uy[i], uz[i], dt ); } ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::BorisRR) { + } else if (algo == ParticlePusherAlgo::BorisRR) { amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { @@ -1688,21 +1694,27 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; - if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ + + // Boris + RR is enabled only for leptons + auto algo = WarpX::particle_pusher_algo; + if (!AmIALepton() && algo == ParticlePusherAlgo::BorisRR) + algo = ParticlePusherAlgo::Boris; + + if (algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { UpdateMomentumBoris( ux[i], uy[i], uz[i], Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); } ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + } else if (algo == ParticlePusherAlgo::Vay) { amrex::ParallelFor( pti.numParticles(), [=] 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::BorisRR) { + } else if (algo == ParticlePusherAlgo::BorisRR) { amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], @@ -2149,3 +2161,10 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const } ); } + +//This function return true if the PhysicalParticleContainer contains electrons +//or positrons, false otherwise +bool +PhysicalParticleContainer::AmIALepton(){ + return (this-> mass == PhysConst::m_e); +} -- cgit v1.2.3 From 4c62803f740dc7f9eebe1861103096e73d0a1eaf Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Mon, 30 Sep 2019 13:46:39 +0200 Subject: Do radiation reaction is now a species property --- Examples/Tests/RadiationReaction/test_RR.3d | 61 -------------- Source/Particles/PhysicalParticleContainer.H | 6 +- Source/Particles/PhysicalParticleContainer.cpp | 92 ++++++++++++++-------- .../Particles/RigidInjectedParticleContainer.cpp | 30 ++++--- Source/Utils/WarpXAlgorithmSelection.H | 3 +- Source/Utils/WarpXAlgorithmSelection.cpp | 1 - 6 files changed, 85 insertions(+), 108 deletions(-) delete mode 100644 Examples/Tests/RadiationReaction/test_RR.3d (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Tests/RadiationReaction/test_RR.3d b/Examples/Tests/RadiationReaction/test_RR.3d deleted file mode 100644 index 9bf25198e..000000000 --- a/Examples/Tests/RadiationReaction/test_RR.3d +++ /dev/null @@ -1,61 +0,0 @@ -max_step = 200 -amr.n_cell = 64 64 64 -amr.blocking_factor = 16 -amr.max_grid_size = 64 -amr.max_level = 0 - -warpx.fine_tag_lo = -0.8 -0.8 -0.8 -warpx.fine_tag_hi = 0.8 0.8 0.8 - -amr.plot_int = 1 # How often to write plotfiles. "<= 0" means no plotfiles. - -warpx.plot_raw_fields = 0 -warpx.plot_finepatch = 0 -warpx.plot_crsepatch = 0 - -# Geometry -geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 1 # Is periodic? -geometry.prob_lo = -5.0e-6 -5.0e-6 -5.0e-6 # physical domain -geometry.prob_hi = 5.0e-6 5.0e-6 5.0e-6 # physical domain - -warpx.B_external = 0.0 0.0 1.0e5 - -# Verbosity -warpx.verbose = 1 - -# Algorithms -algo.particle_pusher = "boris_rr" - -# CFL -warpx.cfl = 1.0 - -# particles -particles.nspecies = 2 -particles.species_names = electron positron - - -electron.charge = -q_e -electron.mass = m_e -electron.injection_style = "SingleParticle" -electron.single_particle_pos = 0.0 0.0 0.0 -electron.single_particle_vel = 100.0 0.0 0.0 # gamma*beta - -positron.charge = q_e -positron.mass = m_e -positron.injection_style = "SingleParticle" -positron.single_particle_pos = 0.0 0.0 0.0 -positron.single_particle_vel = 100.0 0.0 0.0 # gamma*beta - -electron.single_particle_weight = 1.6e-19 -positron.single_particle_weight = 1.6e-19 - -# interpolation -interpolation.nox = 3 -interpolation.noy = 3 -interpolation.noz = 3 - -# Moving window -warpx.do_moving_window = 0 - -warpx.do_dive_cleaning = 1 diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ed95cffb9..9541a3652 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -186,12 +186,14 @@ protected: // Inject particles during the whole simulation void ContinuousInjection (const amrex::RealBox& injection_box) override; -private: - //This function return true if the PhysicalParticleContainer contains electrons //or positrons, false otherwise bool AmIALepton(); + //When true PhysicalParticleContainer tries to use a pusher including + //radiation reaction + bool do_classical_radiation_reaction = false; + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 89953120d..e5553bdf3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -43,6 +43,22 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("do_field_ionization", do_field_ionization); + + //check if Radiation Reaction is enabled and do consistency checks + pp.query("do_classical_radiation_reaction", do_classical_radiation_reaction); + //if the species is not a lepton, do_classical_radiation_reaction + //should be false + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + AmIALepton() || !do_classical_radiation_reaction, + "Can't enable classical radiation reaction for non lepton species. " ); + + //Only Boris pusher is compatible with radiation reaction + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris, + "Radiation reaction can be enabled only if Boris pusher is used"); + //_____________________________ + + #ifdef AMREX_USE_GPU AMREX_ALWAYS_ASSERT_WITH_MESSAGE( do_field_ionization == 0, @@ -1489,44 +1505,40 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const Real m = this-> mass; - // Boris + RR is enabled only for leptons - auto algo = WarpX::particle_pusher_algo; - if (!AmIALepton() && algo == ParticlePusherAlgo::BorisRR) - algo = ParticlePusherAlgo::Boris; - - if (algo == ParticlePusherAlgo::Boris){ + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; if (ion_lev){ qp *= ion_lev[i]; } - UpdateMomentumBoris( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], qp, m, dt); + UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ux[i], uy[i], uz[i], dt ); } ); - } else if (algo == ParticlePusherAlgo::Vay) { + } 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]; } - UpdateMomentumVay( ux[i], uy[i], uz[i], - Ex[i], Ey[i], Ez[i], Bx[i], - By[i], Bz[i], qp, m, dt); + UpdateMomentumBoris( ux[i], uy[i], uz[i], + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ux[i], uy[i], uz[i], dt ); } ); - } else if (algo == ParticlePusherAlgo::BorisRR) { + } 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]; } - UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], + UpdateMomentumVay( ux[i], uy[i], uz[i], Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], qp, m, dt); UpdatePosition( x[i], y[i], z[i], @@ -1615,30 +1627,48 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const Real q = this->charge; const Real m = this-> mass; - // Boris + RR is enabled only for leptons - auto algo = WarpX::particle_pusher_algo; - if (!AmIALepton() && algo == ParticlePusherAlgo::BorisRR) - algo = ParticlePusherAlgo::Boris; + int* AMREX_RESTRICT ion_lev = nullptr; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } - if (algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( pti.numParticles(), + + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumBoris( ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } + UpdateMomentumBorisWithRadiationReaction( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + qp, m, dt); } ); - } else if (algo == ParticlePusherAlgo::Vay) { + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), [=] 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); + 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 (algo == ParticlePusherAlgo::BorisRR) { + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } + UpdateMomentumVay( + ux[i], uy[i], uz[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + qp, m, dt); } ); } else { diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index a3ae5f4e8..24d4ac24c 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -430,25 +430,33 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this->mass; - if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( pti.numParticles(), + + //Assumes that all consistency checks have been done at initialization + if(do_classical_radiation_reaction){ + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumBoris( uxpp[i], uypp[i], uzpp[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + UpdateMomentumBorisWithRadiationReaction( + uxpp[i], uypp[i], uzpp[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); } ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumVay( uxpp[i], uypp[i], uzpp[i], - Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt); + UpdateMomentumBoris( + uxpp[i], uypp[i], uzpp[i], + Expp[i], Eypp[i], Ezpp[i], + Bxpp[i], Bypp[i], Bzpp[i], + q, m, dt); } ); - } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::BorisRR) { - amrex::ParallelFor( - pti.numParticles(), + } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - UpdateMomentumBorisWithRadiationReaction( + UpdateMomentumVay( uxpp[i], uypp[i], uzpp[i], Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 35a15c6fd..54c721abf 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -14,8 +14,7 @@ struct MaxwellSolverAlgo { struct ParticlePusherAlgo { enum { Boris = 0, - Vay = 1, - BorisRR = 2 + Vay = 1 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 8e65cb8d0..8a6ff6dbf 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -18,7 +18,6 @@ const std::map maxwell_solver_algo_to_int = { const std::map particle_pusher_algo_to_int = { {"boris", ParticlePusherAlgo::Boris }, {"vay", ParticlePusherAlgo::Vay }, - {"boris_rr", ParticlePusherAlgo::BorisRR }, {"default", ParticlePusherAlgo::Boris } }; -- cgit v1.2.3 From 43bb052cc31a019e69d3cde5676f59a2c25806a2 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 10 Oct 2019 16:58:29 +0200 Subject: Fixed bug --- Source/Particles/PhysicalParticleContainer.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 283be4e1c..07cef4277 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -48,14 +48,15 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_classical_radiation_reaction", do_classical_radiation_reaction); //if the species is not a lepton, do_classical_radiation_reaction //should be false - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - AmIALepton() || !do_classical_radiation_reaction, + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + !(do_classical_radiation_reaction && !AmIALepton()), "Can't enable classical radiation reaction for non lepton species. " ); //Only Boris pusher is compatible with radiation reaction AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris, - "Radiation reaction can be enabled only if Boris pusher is used"); + !(do_classical_radiation_reaction && + WarpX::particle_pusher_algo != ParticlePusherAlgo::Boris), + "Radiation reaction can be enabled only if Boris pusher is used"); //_____________________________ -- cgit v1.2.3 From 12768748673bb89375b97566f29949ab1c9a2d85 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 10 Oct 2019 17:00:27 +0200 Subject: Style change --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 07cef4277..db969e58f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -54,7 +54,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp //Only Boris pusher is compatible with radiation reaction AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - !(do_classical_radiation_reaction && + !(do_classical_radiation_reaction && WarpX::particle_pusher_algo != ParticlePusherAlgo::Boris), "Radiation reaction can be enabled only if Boris pusher is used"); //_____________________________ -- cgit v1.2.3 From df30827050883f879ab09560ec1d82a04ac3e0f3 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 10 Oct 2019 17:02:15 +0200 Subject: Style change --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index db969e58f..22ece03f7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -48,7 +48,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_classical_radiation_reaction", do_classical_radiation_reaction); //if the species is not a lepton, do_classical_radiation_reaction //should be false - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( !(do_classical_radiation_reaction && !AmIALepton()), "Can't enable classical radiation reaction for non lepton species. " ); -- cgit v1.2.3 From 3e3670639e3f57cdc9f1a43a6c69d7e8308f100d Mon Sep 17 00:00:00 2001 From: lucafedeli88 Date: Fri, 11 Oct 2019 19:14:03 -0500 Subject: Removed merge conflict from file --- Source/Particles/PhysicalParticleContainer.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7f6b2daa9..c3ea3b763 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2205,14 +2205,13 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const ); } -<<<<<<< HEAD //This function return true if the PhysicalParticleContainer contains electrons //or positrons, false otherwise bool PhysicalParticleContainer::AmIALepton(){ return (this-> mass == PhysConst::m_e); } -======= + #ifdef WARPX_QED bool PhysicalParticleContainer::has_quantum_sync() @@ -2239,4 +2238,3 @@ set_quantum_sync_engine_ptr(std::shared_ptr ptr) shr_ptr_qs_engine = ptr; } #endif ->>>>>>> upstream/dev -- cgit v1.2.3