diff options
Diffstat (limited to 'Source')
45 files changed, 2485 insertions, 537 deletions
diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BackTransformedDiagnostic.H index b09f089c1..6300f5dfb 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BackTransformedDiagnostic.H @@ -1,5 +1,5 @@ -#ifndef WARPX_BoostedFrameDiagnostic_H_ -#define WARPX_BoostedFrameDiagnostic_H_ +#ifndef WARPX_BackTransformedDiagnostic_H_ +#define WARPX_BackTransformedDiagnostic_H_ #include <vector> #include <map> @@ -148,7 +148,7 @@ class LabFrameSlice : public LabFrameDiag { }; /** \brief - * BoostedFrameDiagnostic class handles the back-transformation of data when + * BackTransformedDiagnostic class handles the back-transformation of data when * running simulations in a boosted frame of reference to the lab-frame. * Because of the relativity of simultaneity, events that are synchronized * in the simulation boosted frame are not @@ -161,13 +161,13 @@ class LabFrameSlice : public LabFrameDiag { * to the output directory. The functions Flush() and writeLabFrameData() * are called at the end of the simulation and when the * the buffer for data storage is full, respectively. The particle data - * is collected and written only if particle.do_boosted_frame_diagnostic = 1. + * is collected and written only if particle.do_back_transformed_diagnostics = 1. */ -class BoostedFrameDiagnostic { +class BackTransformedDiagnostic { public: - BoostedFrameDiagnostic(amrex::Real zmin_lab, amrex::Real zmax_lab, + BackTransformedDiagnostic(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_slice_snapshots, amrex::Real gamma_boost, diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BackTransformedDiagnostic.cpp index d2a9a8ad7..4881c7b24 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BackTransformedDiagnostic.cpp @@ -1,7 +1,7 @@ #include <AMReX_MultiFabUtil.H> #include <AMReX_MultiFabUtil_C.H> -#include "BoostedFrameDiagnostic.H" +#include "BackTransformedDiagnostic.H" #include "SliceDiagnostic.H" #include "WarpX_f.H" #include "WarpX.H" @@ -468,7 +468,7 @@ void LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) { // Loop over tiles/boxes and in-place convert each slice from boosted - // frame to lab frame. + // frame to back-transformed lab frame. #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -514,8 +514,8 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) } } -BoostedFrameDiagnostic:: -BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, +BackTransformedDiagnostic:: +BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, 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, @@ -533,10 +533,10 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, { - BL_PROFILE("BoostedFrameDiagnostic::BoostedFrameDiagnostic"); + BL_PROFILE("BackTransformedDiagnostic::BackTransformedDiagnostic"); - AMREX_ALWAYS_ASSERT(WarpX::do_boosted_frame_fields or - WarpX::do_boosted_frame_particles); + AMREX_ALWAYS_ASSERT(WarpX::do_back_transformed_fields or + WarpX::do_back_transformed_particles); m_inv_gamma_boost_ = 1.0 / m_gamma_boost_; m_beta_boost_ = std::sqrt(1.0 - m_inv_gamma_boost_*m_inv_gamma_boost_); @@ -585,7 +585,7 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, // Get simulation domain physical coordinates (in boosted frame). RealBox prob_domain_lab = geom.ProbDomain(); // Replace z bounds by lab-frame coordinates - // x and y bounds are the same for lab frame and boosted frame + // x and y bounds are the same for back-transformed 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(); @@ -677,9 +677,9 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, AMREX_ALWAYS_ASSERT(m_max_box_size_ >= m_num_buffer_); } -void BoostedFrameDiagnostic::Flush(const Geometry& geom) +void BackTransformedDiagnostic::Flush(const Geometry& geom) { - BL_PROFILE("BoostedFrameDiagnostic::Flush"); + BL_PROFILE("BackTransformedDiagnostic::Flush"); VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); @@ -694,7 +694,7 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) int i_lab = (m_LabFrameDiags_[i]->m_current_z_lab - zmin_lab) / m_dz_lab_; if (m_LabFrameDiags_[i]->m_buff_counter_ != 0) { - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { const BoxArray& ba = m_LabFrameDiags_[i]->m_data_buffer_->boxArray(); const int hi = ba[0].bigEnd(m_boost_direction_); const int lo = hi - m_LabFrameDiags_[i]->m_buff_counter_ + 1; @@ -729,12 +729,12 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) #endif } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_back_transformed_particles) { // Loop over species to be dumped to BFD - for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) { // Get species name std::string species_name = - species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; + species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)]; #ifdef WARPX_USE_HDF5 // Dump species data writeParticleDataHDF5(m_LabFrameDiags_[i]->m_particles_buffer_[j], @@ -763,12 +763,12 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeLabFrameData(const MultiFab* cell_centered_data, const MultiParticleContainer& mypc, const Geometry& geom, const Real t_boost, const Real dt) { - BL_PROFILE("BoostedFrameDiagnostic::writeLabFrameData"); + BL_PROFILE("BackTransformedDiagnostic::writeLabFrameData"); VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); @@ -806,7 +806,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // If buffer of snapshot i is empty... if ( m_LabFrameDiags_[i]->m_buff_counter_ == 0) { // ... reset fields buffer data_buffer_ - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { m_LabFrameDiags_[i]->m_buff_box_.setSmall(m_boost_direction_, i_lab - m_num_buffer_ + 1); m_LabFrameDiags_[i]->m_buff_box_.setBig(m_boost_direction_, i_lab); @@ -818,12 +818,12 @@ writeLabFrameData(const MultiFab* cell_centered_data, buff_dm, m_ncomp_to_dump, 0) ); } // ... reset particle buffer particles_buffer_[i] - if (WarpX::do_boosted_frame_particles) + if (WarpX::do_back_transformed_particles) m_LabFrameDiags_[i]->m_particles_buffer_.resize( - mypc.nSpeciesBoostedFrameDiags()); + mypc.nSpeciesBackTransformedDiagnostics()); } - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { const int ncomp = cell_centered_data->nComp(); const int start_comp = 0; const bool interpolate = true; @@ -871,7 +871,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, tmp_slice_ptr.reset(nullptr); } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_back_transformed_particles) { if (m_LabFrameDiags_[i]->m_t_lab != prev_t_lab ) { if (tmp_particle_buffer.size()>0) @@ -879,7 +879,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, tmp_particle_buffer.clear(); tmp_particle_buffer.shrink_to_fit(); } - tmp_particle_buffer.resize(mypc.nSpeciesBoostedFrameDiags()); + tmp_particle_buffer.resize(mypc.nSpeciesBackTransformedDiagnostics()); mypc.GetLabFrameData(m_LabFrameDiags_[i]->m_file_name, i_lab, m_boost_direction_, old_z_boost, m_LabFrameDiags_[i]->m_current_z_boost, @@ -887,7 +887,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, tmp_particle_buffer); } m_LabFrameDiags_[i]->AddPartDataToParticleBuffer(tmp_particle_buffer, - mypc.nSpeciesBoostedFrameDiags()); + mypc.nSpeciesBackTransformedDiagnostics()); } ++m_LabFrameDiags_[i]->m_buff_counter_; @@ -895,7 +895,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // If buffer full, write to disk. if (m_LabFrameDiags_[i]->m_buff_counter_ == m_num_buffer_) { - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { #ifdef WARPX_USE_HDF5 Box buff_box = m_LabFrameDiags_[i]->m_buff_box_; @@ -913,12 +913,12 @@ writeLabFrameData(const MultiFab* cell_centered_data, #endif } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_back_transformed_particles) { // Loop over species to be dumped to BFD - for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) { // Get species name const std::string species_name = species_names[ - mypc.mapSpeciesBoostedFrameDiags(j)]; + mypc.mapSpeciesBackTransformedDiagnostics(j)]; #ifdef WARPX_USE_HDF5 // Write data to disk (HDF5) writeParticleDataHDF5(m_LabFrameDiags_[i]->m_particles_buffer_[j], @@ -946,7 +946,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, #ifdef WARPX_USE_HDF5 void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const std::string& species_name) { @@ -994,11 +994,11 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat #endif void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const int i_lab) { - BL_PROFILE("BoostedFrameDiagnostic::writeParticleData"); + BL_PROFILE("BackTransformedDiagnostic::writeParticleData"); std::string field_name; std::ofstream ofs; @@ -1044,10 +1044,10 @@ writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, } void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeMetaData () { - BL_PROFILE("BoostedFrameDiagnostic::writeMetaData"); + BL_PROFILE("BackTransformedDiagnostic::writeMetaData"); if (ParallelDescriptor::IOProcessor()) { const std::string fullpath = WarpX::lab_data_directory + "/snapshots"; @@ -1131,7 +1131,7 @@ LabFrameSnapShot(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in, m_file_num, 5); createLabFrameDirectories(); m_buff_counter_ = 0; - if (WarpX::do_boosted_frame_fields) m_data_buffer_.reset(nullptr); + if (WarpX::do_back_transformed_fields) m_data_buffer_.reset(nullptr); } void @@ -1155,7 +1155,7 @@ createLabFrameDirectories() { if (ParallelDescriptor::IOProcessor()) { - if (WarpX::do_boosted_frame_fields) + if (WarpX::do_back_transformed_fields) { const auto lo = lbound(m_buff_box_); for (int comp = 0; comp < m_ncomp_to_dump_; ++comp) { @@ -1173,15 +1173,15 @@ createLabFrameDirectories() { ParallelDescriptor::Barrier(); - if (WarpX::do_boosted_frame_particles){ + if (WarpX::do_back_transformed_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); // Loop over species to be dumped to BFD - for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) + for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) { // Loop over species to be dumped to BFD std::string species_name = - species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; + species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)]; output_create_species_group(m_file_name, species_name); for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k) { @@ -1208,10 +1208,10 @@ createLabFrameDirectories() { const std::string particles_prefix = "particle"; // Loop over species to be dumped to BFD - for(int i = 0; i < mypc.nSpeciesBoostedFrameDiags(); ++i) { + for(int i = 0; i < mypc.nSpeciesBackTransformedDiagnostics(); ++i) { // Get species name std::string species_name = - species_names[mypc.mapSpeciesBoostedFrameDiags(i)]; + species_names[mypc.mapSpeciesBackTransformedDiagnostics(i)]; const std::string fullpath = m_file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); @@ -1301,7 +1301,7 @@ LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in, m_buff_counter_ = 0; m_particle_slice_dx_lab_ = particle_slice_dx_lab; - if (WarpX::do_boosted_frame_fields) m_data_buffer_.reset(nullptr); + if (WarpX::do_back_transformed_fields) m_data_buffer_.reset(nullptr); } void @@ -1430,10 +1430,10 @@ void LabFrameSlice:: AddPartDataToParticleBuffer( Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer, - int nSpeciesBoostedFrame) { + int nSpeciesBackTransformedDiagnostics) { - for (int isp = 0; isp < nSpeciesBoostedFrame; ++isp) { + for (int isp = 0; isp < nSpeciesBackTransformedDiagnostics; ++isp) { auto np = tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size(); if (np == 0) return; diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package index dfd947d53..b624d6ebe 100644 --- a/Source/Diagnostics/Make.package +++ b/Source/Diagnostics/Make.package @@ -1,9 +1,9 @@ CEXE_sources += WarpXIO.cpp -CEXE_sources += BoostedFrameDiagnostic.cpp +CEXE_sources += BackTransformedDiagnostic.cpp CEXE_sources += ParticleIO.cpp CEXE_sources += FieldIO.cpp CEXE_headers += FieldIO.H -CEXE_headers += BoostedFrameDiagnostic.H +CEXE_headers += BackTransformedDiagnostic.H CEXE_headers += ElectrostaticIO.cpp CEXE_headers += SliceDiagnostic.H CEXE_sources += SliceDiagnostic.cpp diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index d55660b39..3b7481b8c 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -113,7 +113,7 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const } #ifdef WARPX_QED - if(pc->do_qed){ + if(pc->m_do_qed){ real_names.push_back("tau"); } #endif diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 7a3262703..b5fd52bdc 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -141,9 +141,9 @@ WarpX::EvolveEM (int numsteps) bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { std::unique_ptr<MultiFab> cell_centered_data = nullptr; - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { cell_centered_data = GetCellCenteredData(); } myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); @@ -261,7 +261,7 @@ WarpX::EvolveEM (int numsteps) WriteCheckPointFile(); } - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { myBFD->Flush(geom[0]); } diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index 6ecae93e0..4ab2fa022 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -29,21 +29,25 @@ struct InjectorPositionRegular // is a_ppc*(ref_fac**AMREX_SPACEDIM). AMREX_GPU_HOST_DEVICE amrex::XDim3 - getPositionUnitBox (int i_part, int ref_fac=1) const noexcept + getPositionUnitBox (int const i_part, int const ref_fac=1) const noexcept { - int nx = ref_fac*ppc.x; - int ny = ref_fac*ppc.y; + using namespace amrex; + + int const nx = ref_fac*ppc.x; + int const ny = ref_fac*ppc.y; #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) - int nz = ref_fac*ppc.z; + int const nz = ref_fac*ppc.z; #else - int nz = 1; + int const nz = 1; #endif - int ix_part = i_part/(ny*nz); // written this way backward compatibility - int iz_part = (i_part-ix_part*(ny*nz)) / ny; - int iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part; - return amrex::XDim3{(amrex::Real(0.5)+ix_part)/nx, - (amrex::Real(0.5)+iy_part)/ny, - (amrex::Real(0.5)+iz_part) / nz}; + int const ix_part = i_part / (ny*nz); // written this way backward compatibility + int const iz_part = (i_part-ix_part*(ny*nz)) / ny; + int const iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part; + return XDim3{ + (0.5_rt + ix_part) / nx, + (0.5_rt + iy_part) / ny, + (0.5_rt + iz_part) / nz + }; } private: amrex::Dim3 ppc; @@ -100,7 +104,7 @@ struct InjectorPosition // (the union is called Object, and the instance is called object). AMREX_GPU_HOST_DEVICE amrex::XDim3 - getPositionUnitBox (int i_part, int ref_fac=1) const noexcept + getPositionUnitBox (int const i_part, int const ref_fac=1) const noexcept { switch (type) { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 97870ee67..609399ece 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -82,14 +82,14 @@ WarpX::InitData () void WarpX::InitDiagnostics () { - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { const Real* current_lo = geom[0].ProbLo(); const Real* current_hi = geom[0].ProbHi(); Real dt_boost = dt[0]; // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0 Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); - myBFD.reset(new BoostedFrameDiagnostic(zmin_lab, + myBFD.reset(new BackTransformedDiagnostic(zmin_lab, zmax_lab, moving_window_v, dt_snapshots_lab, num_snapshots_lab, @@ -305,18 +305,18 @@ WarpX::InitLevelData (int lev, Real time) { for (int i = 0; i < 3; ++i) { current_fp[lev][i]->setVal(0.0); - Efield_fp[lev][i]->setVal(0.0); - Bfield_fp[lev][i]->setVal(0.0); + Efield_fp[lev][i]->setVal(E_external_grid[i]); + Bfield_fp[lev][i]->setVal(B_external_grid[i]); } if (lev > 0) { for (int i = 0; i < 3; ++i) { - Efield_aux[lev][i]->setVal(0.0); - Bfield_aux[lev][i]->setVal(0.0); + Efield_aux[lev][i]->setVal(E_external_grid[i]); + Bfield_aux[lev][i]->setVal(B_external_grid[i]); current_cp[lev][i]->setVal(0.0); - Efield_cp[lev][i]->setVal(0.0); - Bfield_cp[lev][i]->setVal(0.0); + Efield_cp[lev][i]->setVal(E_external_grid[i]); + Bfield_cp[lev][i]->setVal(B_external_grid[i]); } } diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 8571c74ad..9493672e0 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -26,7 +26,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, { charge = 1.0; mass = std::numeric_limits<Real>::max(); - do_boosted_frame_diags = 0; + do_back_transformed_diagnostics = 0; ParmParse pp(laser_name); @@ -100,7 +100,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, } // Plane normal - Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]); + Real s = 1.0_rt / std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]); nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s }; if (WarpX::gamma_boost > 1.) { @@ -119,19 +119,19 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, } // The first polarization vector - s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]); + s = 1.0_rt / std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]); p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s }; - Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0); + Real const dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, "Laser plane vector is not perpendicular to the main polarization vector"); p_Y = CrossProduct(nvec, p_X); // The second polarization vector - s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]); + s = 1.0_rt / std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]); stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s }; - dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, + Real const dp2 = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp2) < 1.0e-14, "stc_direction is not perpendicular to the laser plane vector"); // Get angle between p_X and stc_direction @@ -266,20 +266,20 @@ LaserParticleContainer::InitData (int lev) position = updated_position; } - auto Transform = [&](int i, int j) -> Vector<Real>{ + auto Transform = [&](int const i, int const j) -> Vector<Real>{ #if (AMREX_SPACEDIM == 3) - return { position[0] + (S_X*(i+0.5))*u_X[0] + (S_Y*(j+0.5))*u_Y[0], - position[1] + (S_X*(i+0.5))*u_X[1] + (S_Y*(j+0.5))*u_Y[1], - position[2] + (S_X*(i+0.5))*u_X[2] + (S_Y*(j+0.5))*u_Y[2] }; + return { position[0] + (S_X*(Real(i)+0.5_rt))*u_X[0] + (S_Y*(Real(j)+0.5_rt))*u_Y[0], + position[1] + (S_X*(Real(i)+0.5_rt))*u_X[1] + (S_Y*(Real(j)+0.5_rt))*u_Y[1], + position[2] + (S_X*(Real(i)+0.5_rt))*u_X[2] + (S_Y*(Real(j)+0.5_rt))*u_Y[2] }; #else # if (defined WARPX_DIM_RZ) - return { position[0] + (S_X*(i+0.5)), + return { position[0] + (S_X*(Real(i)+0.5)), 0.0, position[2]}; # else - return { position[0] + (S_X*(i+0.5))*u_X[0], + return { position[0] + (S_X*(Real(i)+0.5))*u_X[0], 0.0, - position[2] + (S_X*(i+0.5))*u_X[2] }; + position[2] + (S_X*(Real(i)+0.5))*u_X[2] }; # endif #endif }; @@ -449,9 +449,9 @@ LaserParticleContainer::Evolve (int lev, #endif { #ifdef _OPENMP - int thread_num = omp_get_thread_num(); + int const thread_num = omp_get_thread_num(); #else - int thread_num = 0; + int const thread_num = 0; #endif Cuda::ManagedDeviceVector<Real> plane_Xp, plane_Yp, amplitude_E; @@ -610,7 +610,7 @@ void LaserParticleContainer::ComputeWeightMobility (Real Sx, Real Sy) { constexpr Real eps = 0.01; - constexpr Real fac = 1.0/(2.0*MathConst::pi*PhysConst::mu0*PhysConst::c*PhysConst::c*eps); + constexpr Real fac = 1.0_rt / (2.0_rt * MathConst::pi * PhysConst::mu0 * PhysConst::c * PhysConst::c * eps); weight = fac * wavelength * Sx * Sy / std::min(Sx,Sy) * e_max; // The mobility is the constant of proportionality between the field to diff --git a/Source/Laser/LaserProfiles.cpp b/Source/Laser/LaserProfiles.cpp index 281ab2101..44411cedf 100644 --- a/Source/Laser/LaserProfiles.cpp +++ b/Source/Laser/LaserProfiles.cpp @@ -28,16 +28,16 @@ LaserParticleContainer::gaussian_laser_profile ( const Real oscillation_phase = k0 * PhysConst::c * ( t - profile_t_peak ); // The coefficients below contain info about Gouy phase, // laser diffraction, and phase front curvature - const Complex diffract_factor = Real(1.) + I * profile_focal_distance - * Real(2.)/( k0 * profile_waist * profile_waist ); - const Complex inv_complex_waist_2 = Real(1.)/( profile_waist*profile_waist * diffract_factor ); + const Complex diffract_factor = 1._rt + I * profile_focal_distance + * 2._rt/( k0 * profile_waist * profile_waist ); + const Complex inv_complex_waist_2 = 1._rt / ( profile_waist*profile_waist * diffract_factor ); // Time stretching due to STCs and phi2 complex envelope // (1 if zeta=0, beta=0, phi2=0) - const Complex stretch_factor = Real(1.) + Real(4.) * + const Complex stretch_factor = 1._rt + 4._rt * (zeta+beta*profile_focal_distance) * (zeta+beta*profile_focal_distance) * (inv_tau2*inv_complex_waist_2) + - Real(2.)*I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; + 2._rt *I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; // Amplitude and monochromatic oscillations Complex prefactor = e_max * MathFunc::exp( I * oscillation_phase ); @@ -61,10 +61,10 @@ LaserParticleContainer::gaussian_laser_profile ( amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) { - const Complex stc_exponent = Real(1.)/stretch_factor * inv_tau2 * + const Complex stc_exponent = 1._rt / stretch_factor * inv_tau2 * MathFunc::pow((t - tmp_profile_t_peak - tmp_beta*k0*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) - - Real(2.)*I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) + 2._rt *I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) *( tmp_zeta - tmp_beta*tmp_profile_focal_distance ) * inv_complex_waist_2),2); // stcfactor = everything but complex transverse envelope const Complex stcfactor = prefactor * MathFunc::exp( - stc_exponent ); diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 54baecbf6..b81fed147 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -83,11 +83,6 @@ ifeq ($(USE_SENSEI_INSITU),TRUE) endif ifeq ($(QED),TRUE) - BOOST_ROOT ?= NOT_SET - ifneq ($(BOOST_ROOT),NOT_SET) - VPATH_LOCATIONS += $(BOOST_ROOT) - INCLUDE_LOCATIONS += $(BOOST_ROOT) - endif INCLUDE_LOCATIONS += $(PICSAR_HOME)/src/multi_physics/QED/src CXXFLAGS += -DWARPX_QED CFLAGS += -DWARPX_QED @@ -95,8 +90,20 @@ ifeq ($(QED),TRUE) F90FLAGS += -DWARPX_QED include $(WARPX_HOME)/Source/QED/Make.package USERSuffix := $(USERSuffix).QED -endif + ifeq ($(QED_TABLE_GEN),TRUE) + BOOST_ROOT ?= NOT_SET + ifneq ($(BOOST_ROOT),NOT_SET) + VPATH_LOCATIONS += $(BOOST_ROOT) + INCLUDE_LOCATIONS += $(BOOST_ROOT) + endif + CXXFLAGS += -DWARPX_QED_TABLE_GEN + CFLAGS += -DWARPX_QED_TABLE_GEN + FFLAGS += -DWARPX_QED_TABLE_GEN + F90FLAGS += -DWARPX_QED_TABLE_GEN + USERSuffix := $(USERSuffix).GENTABLES + endif +endif include $(PICSAR_HOME)/src/Make.package diff --git a/Source/Parallelization/InterpolateCurrentFineToCoarse.H b/Source/Parallelization/InterpolateCurrentFineToCoarse.H index 148b725d0..cbbcdfab5 100644 --- a/Source/Parallelization/InterpolateCurrentFineToCoarse.H +++ b/Source/Parallelization/InterpolateCurrentFineToCoarse.H @@ -56,6 +56,8 @@ public: int const k ) const noexcept { + using namespace amrex; + auto const & fine_unsafe = m_fine; // out-of-bounds access not secured with zero-values yet auto const & coarse = m_coarse; // out-of-bounds access not secured but will also not occur @@ -71,29 +73,29 @@ public: int const kk = k * m_refinement_ratio; #if AMREX_SPACEDIM == 2 if (IDim == 0) { - coarse(i, j, k) = 0.25 * ( + coarse(i, j, k) = 0.25_rt * ( fine(ii, jj, kk) + fine(ii + 1, jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + fine(ii, jj + 1, kk) + fine(ii + 1, jj + 1, kk) ) ); } else if (IDim == 2) { - coarse(i, j, k) = 0.25 * ( + coarse(i, j, k) = 0.25_rt * ( fine(ii, jj, kk) + fine(ii, jj + 1, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii - 1, jj, kk) + fine(ii - 1, jj + 1, kk) + fine(ii + 1, jj, kk) + fine(ii + 1, jj + 1, kk) ) ); } else { - coarse(i, j, k) = 0.25 * ( + coarse(i, j, k) = 0.25_rt * ( fine(ii, jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii - 1, jj , kk) + fine(ii + 1, jj , kk) + fine(ii , jj - 1, kk) + fine(ii , jj + 1, kk) ) + - 0.25 * ( + 0.25_rt * ( fine(ii - 1, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + fine(ii - 1, jj + 1, kk) + fine(ii + 1, jj + 1, kk) ) @@ -101,64 +103,64 @@ public: } #elif AMREX_SPACEDIM == 3 if (IDim == 0) { - coarse(i,j,k) = 0.125 * ( + coarse(i,j,k) = 0.125_rt * ( fine(ii , jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) + fine(ii , jj , kk-1) + fine(ii , jj , kk+1) ) + - 0.25 * ( + 0.25_rt * ( fine(ii , jj-1, kk-1) + fine(ii , jj+1, kk-1) + fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) ) + fine(ii+1, jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii+1, jj-1, kk ) + fine(ii+1, jj+1, kk ) + fine(ii+1, jj , kk-1) + fine(ii+1, jj , kk+1) ) + - 0.25 * ( + 0.25_rt * ( fine(ii+1, jj-1, kk-1) + fine(ii+1, jj+1, kk-1) + fine(ii+1, jj-1, kk+1) + fine(ii+1, jj+1, kk+1) ) ); } else if (IDim == 1) { - coarse(i, j, k) = 0.125 * ( + coarse(i, j, k) = 0.125_rt * ( fine(ii, jj , kk) + - 0.5 * ( + 0.5_rt * ( fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + fine(ii , jj , kk-1) + fine(ii , jj , kk+1) - ) + - 0.25 * ( + ) + + 0.25_rt * ( fine(ii-1, jj , kk-1) + fine(ii+1, jj , kk-1) + fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) - ) + - fine(ii, jj+1, kk) + - 0.5 * ( + ) + + fine(ii, jj+1, kk) + + 0.5_rt * ( fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) + fine(ii , jj+1, kk-1) + fine(ii , jj+1, kk+1) - ) + - 0.25 * ( + ) + + 0.25_rt * ( fine(ii-1, jj+1, kk-1) + fine(ii+1, jj+1, kk-1) + fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) - ) + ) ); } else { - coarse(i, j, k) = 0.125 * ( + coarse(i, j, k) = 0.125_rt * ( fine(ii, jj, kk ) + - 0.5 * ( + 0.5_rt * ( fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) ) + - 0.25 * ( + 0.25_rt * ( fine(ii-1, jj-1, kk ) + fine(ii+1, jj-1, kk ) + fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) ) + fine(ii, jj, kk+1) + - 0.5 * ( + 0.5_rt * ( fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) + fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) ) + - 0.25 * ( + 0.25_rt * ( fine(ii-1, jj-1, kk+1) + fine(ii+1, jj-1, kk+1) + fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) ) diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index 5da867c9f..169cd0ee1 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -5,38 +5,38 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_bfield_x (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bxa, + amrex::Array4<amrex::Real > const& Bxa, amrex::Array4<amrex::Real const> const& Bxf, amrex::Array4<amrex::Real const> const& Bxc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wx = (j == jg*2) ? 0.0 : 0.5; - Real owx = 1.0-wx; + Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt; + Real const owx = 1.0_rt - wx; Bxa(j,k,l) = owx * Bxc(jg,kg,lg) + wx * Bxc(jg+1,kg,lg) + Bxf(j,k,l); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_bfield_y (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bya, + amrex::Array4<amrex::Real > const& Bya, amrex::Array4<amrex::Real const> const& Byf, amrex::Array4<amrex::Real const> const& Byc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); // Note that for 2d, l=0, because the amrex convention is used here. #if (AMREX_SPACEDIM == 3) - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; Bya(j,k,l) = owy * Byc(jg,kg,lg) + wy * Byc(jg,kg+1,lg) + Byf(j,k,l); #else Bya(j,k,l) = Byc(jg,kg,lg) + Byf(j,k,l); @@ -45,47 +45,47 @@ void warpx_interp_bfield_y (int j, int k, int l, AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_bfield_z (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bza, + amrex::Array4<amrex::Real > const& Bza, amrex::Array4<amrex::Real const> const& Bzf, amrex::Array4<amrex::Real const> const& Bzc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); // Note that for 2d, l=0, because the amrex convention is used here. #if (AMREX_SPACEDIM == 3) - Real wz = (l == lg*2) ? 0.0 : 0.5; - Real owz = 1.0-wz; + Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt; + Real const owz = 1.0_rt - wz; Bza(j,k,l) = owz * Bzc(jg,kg,lg) + owz * Bzc(jg,kg,lg+1) + Bzf(j,k,l); #else - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; Bza(j,k,l) = owy * Bzc(jg,kg,lg) + owy * Bzc(jg,kg+1,lg) + Bzf(j,k,l); #endif } AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_efield_x (int j, int k, int l, - amrex::Array4<amrex::Real> const& Exa, + amrex::Array4<amrex::Real > const& Exa, amrex::Array4<amrex::Real const> const& Exf, amrex::Array4<amrex::Real const> const& Exc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; #if (AMREX_SPACEDIM == 3) - Real wz = (l == lg*2) ? 0.0 : 0.5; - Real owz = 1.0-wz; + Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt; + Real const owz = 1.0_rt - wz; Exa(j,k,l) = owy * owz * Exc(jg ,kg ,lg ) + wy * owz * Exc(jg ,kg+1,lg ) + owy * wz * Exc(jg ,kg ,lg+1) @@ -98,30 +98,30 @@ void warpx_interp_efield_x (int j, int k, int l, AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_efield_y (int j, int k, int l, - amrex::Array4<amrex::Real> const& Eya, + amrex::Array4<amrex::Real > const& Eya, amrex::Array4<amrex::Real const> const& Eyf, amrex::Array4<amrex::Real const> const& Eyc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wx = (j == jg*2) ? 0.0 : 0.5; - Real owx = 1.0-wx; + Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt; + Real const owx = 1.0_rt - wx; #if (AMREX_SPACEDIM == 3) - Real wz = (l == lg*2) ? 0.0 : 0.5; - Real owz = 1.0-wz; + Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt; + Real const owz = 1.0_rt - wz; Eya(j,k,l) = owx * owz * Eyc(jg ,kg ,lg ) + wx * owz * Eyc(jg+1,kg ,lg ) + owx * wz * Eyc(jg ,kg ,lg+1) + wx * wz * Eyc(jg+1,kg ,lg+1) + Eyf(j,k,l); #else - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; Eya(j,k,l) = owx * owy * Eyc(jg ,kg ,lg) + wx * owy * Eyc(jg+1,kg ,lg) + owx * wy * Eyc(jg ,kg+1,lg) @@ -132,22 +132,22 @@ void warpx_interp_efield_y (int j, int k, int l, AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_efield_z (int j, int k, int l, - amrex::Array4<amrex::Real> const& Eza, + amrex::Array4<amrex::Real > const& Eza, amrex::Array4<amrex::Real const> const& Ezf, amrex::Array4<amrex::Real const> const& Ezc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wx = (j == jg*2) ? 0.0 : 0.5; - Real owx = 1.0-wx; + Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt; + Real const owx = 1.0_rt - wx; #if (AMREX_SPACEDIM == 3) - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real owy = 1.0_rt - wy; Eza(j,k,l) = owx * owy * Ezc(jg ,kg ,lg ) + wx * owy * Ezc(jg+1,kg ,lg ) + owx * wy * Ezc(jg ,kg+1,lg ) diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 6da0f1155..2737eb008 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -4,6 +4,9 @@ #include "ShapeFactors.H" #include <WarpX_Complex.H> +#include <AMReX_Array4.H> +#include <AMReX_REAL.H> + /* \brief Current Deposition for thread thread_num * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. @@ -208,69 +211,71 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, const amrex::Real q, const long n_rz_azimuthal_modes) { + using namespace amrex; + // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) - const bool do_ionization = ion_lev; - const amrex::Real dxi = 1.0/dx[0]; - const amrex::Real dtsdx0 = dt*dxi; - const amrex::Real xmin = xyzmin[0]; + bool const do_ionization = ion_lev; + Real const dxi = 1.0_rt / dx[0]; + Real const dtsdx0 = dt*dxi; + Real const xmin = xyzmin[0]; #if (defined WARPX_DIM_3D) - const amrex::Real dyi = 1.0/dx[1]; - const amrex::Real dtsdy0 = dt*dyi; - const amrex::Real ymin = xyzmin[1]; + Real const dyi = 1.0_rt / dx[1]; + Real const dtsdy0 = dt*dyi; + Real const ymin = xyzmin[1]; #endif - const amrex::Real dzi = 1.0/dx[2]; - const amrex::Real dtsdz0 = dt*dzi; - const amrex::Real zmin = xyzmin[2]; + Real const dzi = 1.0_rt / dx[2]; + Real const dtsdz0 = dt*dzi; + Real const zmin = xyzmin[2]; #if (defined WARPX_DIM_3D) - const amrex::Real invdtdx = 1.0/(dt*dx[1]*dx[2]); - const amrex::Real invdtdy = 1.0/(dt*dx[0]*dx[2]); - const amrex::Real invdtdz = 1.0/(dt*dx[0]*dx[1]); + Real const invdtdx = 1.0_rt / (dt*dx[1]*dx[2]); + Real const invdtdy = 1.0_rt / (dt*dx[0]*dx[2]); + Real const invdtdz = 1.0_rt / (dt*dx[0]*dx[1]); #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) - const amrex::Real invdtdx = 1.0/(dt*dx[2]); - const amrex::Real invdtdz = 1.0/(dt*dx[0]); - const amrex::Real invvol = 1.0/(dx[0]*dx[2]); + Real const invdtdx = 1.0_rt / (dt*dx[2]); + Real const invdtdz = 1.0_rt / (dt*dx[0]); + Real const invvol = 1.0_rt / (dx[0]*dx[2]); #endif #if (defined WARPX_DIM_RZ) - const Complex I = Complex{0., 1.}; + Complex const I = Complex{0., 1.}; #endif - const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c; + Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c ); // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr amrex::ParallelFor( np_to_depose, - [=] AMREX_GPU_DEVICE (long ip) { + [=] AMREX_GPU_DEVICE (long const ip) { // --- Get particle quantities - const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + Real const gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + uyp[ip]*uyp[ip]*clightsq + uzp[ip]*uzp[ip]*clightsq); // wqx, wqy wqz are particle current in each direction - amrex::Real wq = q*wp[ip]; + Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; } - const amrex::Real wqx = wq*invdtdx; + Real const wqx = wq*invdtdx; #if (defined WARPX_DIM_3D) - const amrex::Real wqy = wq*invdtdy; + Real const wqy = wq*invdtdy; #endif - const amrex::Real wqz = wq*invdtdz; + Real const wqz = wq*invdtdz; // computes current and old position in grid units #if (defined WARPX_DIM_RZ) - const amrex::Real xp_mid = xp[ip] - 0.5*dt*uxp[ip]*gaminv; - const amrex::Real yp_mid = yp[ip] - 0.5*dt*uyp[ip]*gaminv; - const amrex::Real xp_old = xp[ip] - dt*uxp[ip]*gaminv; - const amrex::Real yp_old = yp[ip] - dt*uyp[ip]*gaminv; - const amrex::Real rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); - const amrex::Real rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid); - const amrex::Real rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old); - amrex::Real costheta_new, sintheta_new; - if (rp_new > 0.) { + Real const xp_mid = xp[ip] - 0.5_rt * dt*uxp[ip]*gaminv; + Real const yp_mid = yp[ip] - 0.5_rt * dt*uyp[ip]*gaminv; + Real const xp_old = xp[ip] - dt*uxp[ip]*gaminv; + Real const yp_old = yp[ip] - dt*uyp[ip]*gaminv; + Real const rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid); + Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old); + Real costheta_new, sintheta_new; + if (rp_new > 0._rt) { costheta_new = xp[ip]/rp_new; sintheta_new = yp[ip]/rp_new; } else { @@ -278,7 +283,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, sintheta_new = 0.; } amrex::Real costheta_mid, sintheta_mid; - if (rp_mid > 0.) { + if (rp_mid > 0._rt) { costheta_mid = xp_mid/rp_mid; sintheta_mid = yp_mid/rp_mid; } else { @@ -286,7 +291,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, sintheta_mid = 0.; } amrex::Real costheta_old, sintheta_old; - if (rp_old > 0.) { + if (rp_old > 0._rt) { costheta_old = xp_old/rp_old; sintheta_old = yp_old/rp_old; } else { @@ -296,37 +301,37 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, const Complex xy_new0 = Complex{costheta_new, sintheta_new}; const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid}; const Complex xy_old0 = Complex{costheta_old, sintheta_old}; - const amrex::Real x_new = (rp_new - xmin)*dxi; - const amrex::Real x_old = (rp_old - xmin)*dxi; + Real const x_new = (rp_new - xmin)*dxi; + Real const x_old = (rp_old - xmin)*dxi; #else - const amrex::Real x_new = (xp[ip] - xmin)*dxi; - const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv; + Real const x_new = (xp[ip] - xmin)*dxi; + Real const x_old = x_new - dtsdx0*uxp[ip]*gaminv; #endif #if (defined WARPX_DIM_3D) - const amrex::Real y_new = (yp[ip] - ymin)*dyi; - const amrex::Real y_old = y_new - dtsdy0*uyp[ip]*gaminv; + Real const y_new = (yp[ip] - ymin)*dyi; + Real const y_old = y_new - dtsdy0*uyp[ip]*gaminv; #endif - const amrex::Real z_new = (zp[ip] - zmin)*dzi; - const amrex::Real z_old = z_new - dtsdz0*uzp[ip]*gaminv; + Real const z_new = (zp[ip] - zmin)*dzi; + Real const z_old = z_new - dtsdz0*uzp[ip]*gaminv; #if (defined WARPX_DIM_RZ) - const amrex::Real vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; + Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; #elif (defined WARPX_DIM_XZ) - const amrex::Real vy = uyp[ip]*gaminv; + Real const vy = uyp[ip]*gaminv; #endif // Shape factor arrays // Note that there are extra values above and below // to possibly hold the factor for the old particle // which can be at a different grid location. - amrex::Real sx_new[depos_order + 3] = {0.}; - amrex::Real sx_old[depos_order + 3] = {0.}; + Real sx_new[depos_order + 3] = {0.}; + Real sx_old[depos_order + 3] = {0.}; #if (defined WARPX_DIM_3D) - amrex::Real sy_new[depos_order + 3] = {0.}; - amrex::Real sy_old[depos_order + 3] = {0.}; + Real sy_new[depos_order + 3] = {0.}; + Real sy_old[depos_order + 3] = {0.}; #endif - amrex::Real sz_new[depos_order + 3] = {0.}; - amrex::Real sz_old[depos_order + 3] = {0.}; + Real sz_new[depos_order + 3] = {0.}; + Real sz_old[depos_order + 3] = {0.}; // --- Compute shape factors // Compute shape factors for position as they are now and at old positions @@ -397,7 +402,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djr_cmplx = amrex::Real(2.)*sdxi*xy_mid; + const Complex djr_cmplx = 2._rt *sdxi*xy_mid; amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real()); amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag()); xy_mid = xy_mid*xy_mid0; @@ -407,8 +412,8 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, } for (int k=dkl; k<=depos_order+2-dku; k++) { for (int i=dil; i<=depos_order+2-diu; i++) { - const amrex::Real sdyj = wq*vy*invvol*((sz_new[k] + 0.5*(sz_old[k] - sz_new[k]))*sx_new[i] + - (0.5*sz_new[k] + 1./3.*(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i])); + Real const sdyj = wq*vy*invvol*((sz_new[k] + 0.5_rt * (sz_old[k] - sz_new[k]))*sx_new[i] + + (0.5_rt * sz_new[k] + 1._rt / 3._rt *(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i])); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj); #if (defined WARPX_DIM_RZ) Complex xy_new = xy_new0; @@ -418,7 +423,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes // The minus sign comes from the different convention with respect to Davidson et al. - const Complex djt_cmplx = -amrex::Real(2.)*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode* + const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode* (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old)); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag()); @@ -430,15 +435,15 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, } } for (int i=dil; i<=depos_order+2-diu; i++) { - amrex::Real sdzk = 0.; + Real sdzk = 0.; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5*(sx_old[i] - sx_new[i])); + sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5_rt * (sx_old[i] - sx_new[i])); amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk); #if (defined WARPX_DIM_RZ) Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djz_cmplx = amrex::Real(2.)*sdzk*xy_mid; + const Complex djz_cmplx = 2._rt * sdzk * xy_mid; amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real()); amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag()); xy_mid = xy_mid*xy_mid0; diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index c9d520873..a6c124ddd 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -19,6 +19,7 @@ include $(WARPX_HOME)/Source/Particles/Pusher/Make.package include $(WARPX_HOME)/Source/Particles/Deposition/Make.package include $(WARPX_HOME)/Source/Particles/Gather/Make.package include $(WARPX_HOME)/Source/Particles/Sorting/Make.package +include $(WARPX_HOME)/Source/Particles/ParticleCreation/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 58546a106..f9a0e51d7 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -1,14 +1,15 @@ - #ifndef WARPX_ParticleContainer_H_ #define WARPX_ParticleContainer_H_ -#include <AMReX_Particles.H> +#include "ElementaryProcess.H" + #include <WarpXParticleContainer.H> #include <PhysicalParticleContainer.H> #include <RigidInjectedParticleContainer.H> #include <PhotonParticleContainer.H> #include <LaserParticleContainer.H> +#include <AMReX_Particles.H> #ifdef WARPX_QED #include <BreitWheelerEngineWrapper.H> #include <QuantumSyncEngineWrapper.H> @@ -18,6 +19,8 @@ #include <map> #include <string> #include <algorithm> +#include <utility> +#include <tuple> // // MultiParticleContainer holds multiple (nspecies or npsecies+1 when @@ -162,9 +165,9 @@ public: int nSpecies() const {return nspecies;} - int nSpeciesBoostedFrameDiags() const {return nspecies_boosted_frame_diags;} - int mapSpeciesBoostedFrameDiags(int i) const {return map_species_boosted_frame_diags[i];} - int doBoostedFrameDiags() const {return do_boosted_frame_diags;} + int nSpeciesBackTransformedDiagnostics() const {return nspecies_back_transformed_diagnostics;} + int mapSpeciesBackTransformedDiagnostics(int i) const {return map_species_back_transformed_diagnostics[i];} + int doBackTransformedDiagnostics() const {return do_back_transformed_diagnostics;} int nSpeciesDepositOnMainGrid () const { bool const onMainGrid = true; @@ -196,6 +199,8 @@ public: PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } + IonizationProcess ionization_process; + protected: // Particle container types @@ -215,13 +220,58 @@ protected: #ifdef WARPX_QED // The QED engines - BreitWheelerEngine bw_engine; - QuantumSynchrotronEngine qs_engine; + std::shared_ptr<BreitWheelerEngine> m_shr_p_bw_engine; + std::shared_ptr<QuantumSynchrotronEngine> m_shr_p_qs_engine; //_______________________________ - //Initialize QED engines and provides smart pointers - //to species who need QED processes + /** + * Initialize QED engines and provides smart pointers + * to species who need QED processes + */ void InitQED (); + + //Variables to store how many species need a QED process + int m_nspecies_quantum_sync = 0; + int m_nspecies_breit_wheeler = 0; + //________ + + /** + * Returns the number of species having Quantum Synchrotron process enabled + */ + int NSpeciesQuantumSync() const { return m_nspecies_quantum_sync;} + + /** + * Returns the number of species having Breit Wheeler process enabled + */ + int NSpeciesBreitWheeler() const { return m_nspecies_breit_wheeler;} + + /** + * Initializes the Quantum Synchrotron engine + */ + void InitQuantumSync (); + + /** + * Initializes the Quantum Synchrotron engine + */ + void InitBreitWheeler (); + + /** + * Parses inputfile parameters for Quantum Synchrotron engine + * @return {a tuple containing a flag which is true if tables + * have to be generate, a filename (where tables should be stored + * or read from) and control parameters.} + */ + std::tuple<bool, std::string, PicsarQuantumSynchrotronCtrl> + ParseQuantumSyncParams (); + + /** + * Parses inputfile parameters for Breit Wheeler engine + * @return {a tuple containing a flag which is true if tables + * have to be generate, a filename (where tables should be stored + * or read from) and control parameters.} + */ + std::tuple<bool, std::string, PicsarBreitWheelerCtrl> + ParseBreitWheelerParams (); #endif private: @@ -236,12 +286,12 @@ private: void mapSpeciesProduct (); int getSpeciesID (std::string product_str); - // Number of species dumped in BoostedFrameDiagnostics - int nspecies_boosted_frame_diags = 0; - // map_species_boosted_frame_diags[i] is the species ID in - // MultiParticleContainer for 0<i<nspecies_boosted_frame_diags - std::vector<int> map_species_boosted_frame_diags; - int do_boosted_frame_diags = 0; + // Number of species dumped in BackTransformedDiagnostics + int nspecies_back_transformed_diagnostics = 0; + // map_species_back_transformed_diagnostics[i] is the species ID in + // MultiParticleContainer for 0<i<nspecies_back_transformed_diagnostics + std::vector<int> map_species_back_transformed_diagnostics; + int do_back_transformed_diagnostics = 0; // runtime parameters int nlasers = 0; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index c860d21f5..fca22daa2 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,11 +1,16 @@ #include <MultiParticleContainer.H> + +#include <AMReX_Vector.H> + #include <WarpX_f.H> #include <WarpX.H> +//This is now needed for writing a binary file on disk. +#include <WarpXUtil.H> + #include <limits> #include <algorithm> #include <string> -#include <memory> using namespace amrex; @@ -38,16 +43,17 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) // Compute the number of species for which lab-frame data is dumped // nspecies_lab_frame_diags, and map their ID to MultiParticleContainer // particle IDs in map_species_lab_diags. - map_species_boosted_frame_diags.resize(nspecies); - nspecies_boosted_frame_diags = 0; + map_species_back_transformed_diagnostics.resize(nspecies); + nspecies_back_transformed_diagnostics = 0; for (int i=0; i<nspecies; i++){ auto& pc = allcontainers[i]; - if (pc->do_boosted_frame_diags){ - map_species_boosted_frame_diags[nspecies_boosted_frame_diags] = i; - do_boosted_frame_diags = 1; - nspecies_boosted_frame_diags += 1; + if (pc->do_back_transformed_diagnostics){ + map_species_back_transformed_diagnostics[nspecies_back_transformed_diagnostics] = i; + do_back_transformed_diagnostics = 1; + nspecies_back_transformed_diagnostics += 1; } } + ionization_process = IonizationProcess(); } void @@ -315,7 +321,11 @@ void MultiParticleContainer::Redistribute () { for (auto& pc : allcontainers) { - pc->Redistribute(); + if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { + pc->RedistributeCPU(); + } else { + pc->Redistribute(); + } } } @@ -323,7 +333,11 @@ void MultiParticleContainer::RedistributeLocal (const int num_ghost) { for (auto& pc : allcontainers) { - pc->Redistribute(0, 0, 0, num_ghost); + if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { + pc->RedistributeCPU(0, 0, 0, num_ghost); + } else { + pc->Redistribute(0, 0, 0, num_ghost); + } } } @@ -387,8 +401,8 @@ MultiParticleContainer BL_PROFILE("MultiParticleContainer::GetLabFrameData"); // Loop over particle species - for (int i = 0; i < nspecies_boosted_frame_diags; ++i){ - int isp = map_species_boosted_frame_diags[i]; + for (int i = 0; i < nspecies_back_transformed_diagnostics; ++i){ + int isp = map_species_back_transformed_diagnostics[i]; WarpXParticleContainer* pc = allcontainers[isp].get(); WarpXParticleContainer::DiagnosticParticles diagnostic_particles; pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); @@ -526,146 +540,6 @@ MultiParticleContainer::getSpeciesID (std::string product_str) return i_product; } -namespace -{ - // For particle i in mfi, if is_ionized[i]=1, copy particle - // particle i from container pc_source into pc_product - void createIonizedParticles ( - int lev, const MFIter& mfi, - std::unique_ptr< WarpXParticleContainer>& pc_source, - std::unique_ptr< WarpXParticleContainer>& pc_product, - amrex::Gpu::ManagedDeviceVector<int>& is_ionized) - { - BL_PROFILE("createIonizedParticles"); - - const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - // Get source particle data - auto& ptile_source = pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - const int np_source = ptile_source.GetArrayOfStructs().size(); - if (np_source == 0) return; - // --- source AoS particle data - WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); - // --- source SoA particle data - auto& soa_source = ptile_source.GetStructOfArrays(); - GpuArray<ParticleReal*,PIdx::nattribs> attribs_source; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - attribs_source[ia] = soa_source.GetRealData(ia).data(); - } - // --- source runtime attribs - GpuArray<ParticleReal*,3> runtime_uold_source; - // Prepare arrays for boosted frame diagnostics. - runtime_uold_source[0] = soa_source.GetRealData(PIdx::ux).data(); - runtime_uold_source[1] = soa_source.GetRealData(PIdx::uy).data(); - runtime_uold_source[2] = soa_source.GetRealData(PIdx::uz).data(); - - // Indices of product particle for each ionized source particle. - // i_product[i]-1 is the location in product tile of product particle - // from source particle i. - amrex::Gpu::ManagedDeviceVector<int> i_product; - i_product.resize(np_source); - // 0<i<np_source - // 0<i_product<np_ionized - // Strictly speaking, i_product should be an exclusive_scan of - // is_ionized. However, for indices where is_ionized is 1, the - // inclusive scan gives the same result with an offset of 1. - // The advantage of inclusive_scan is that the sum of is_ionized - // is in the last element, so no other reduction is required to get - // number of particles. - // Gpu::inclusive_scan runs on the current GPU stream, and synchronizes - // with the CPU, so that the next line (executed by the CPU) has the - // updated values of i_product - amrex::Gpu::inclusive_scan(is_ionized.begin(), is_ionized.end(), i_product.begin()); - int np_ionized = i_product[np_source-1]; - if (np_ionized == 0) return; - int* AMREX_RESTRICT p_i_product = i_product.dataPtr(); - - // Get product particle data - auto& ptile_product = pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - // old and new (i.e., including ionized particles) number of particles - // for product species - const int np_product_old = ptile_product.GetArrayOfStructs().size(); - const int np_product_new = np_product_old + np_ionized; - // Allocate extra space in product species for ionized particles. - ptile_product.resize(np_product_new); - // --- product AoS particle data - // First element is the first newly-created product particle - WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; - // --- product SoA particle data - auto& soa_product = ptile_product.GetStructOfArrays(); - GpuArray<ParticleReal*,PIdx::nattribs> attribs_product; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - // First element is the first newly-created product particle - attribs_product[ia] = soa_product.GetRealData(ia).data() + np_product_old; - } - // --- product runtime attribs - GpuArray<ParticleReal*,6> runtime_attribs_product; - bool do_boosted_product = WarpX::do_boosted_frame_diagnostic - && pc_product->DoBoostedFrameDiags(); - if (do_boosted_product) { - std::map<std::string, int> comps_product = pc_product->getParticleComps(); - runtime_attribs_product[0] = soa_product.GetRealData(comps_product[ "xold"]).data() + np_product_old; - runtime_attribs_product[1] = soa_product.GetRealData(comps_product[ "yold"]).data() + np_product_old; - runtime_attribs_product[2] = soa_product.GetRealData(comps_product[ "zold"]).data() + np_product_old; - runtime_attribs_product[3] = soa_product.GetRealData(comps_product["uxold"]).data() + np_product_old; - runtime_attribs_product[4] = soa_product.GetRealData(comps_product["uyold"]).data() + np_product_old; - runtime_attribs_product[5] = soa_product.GetRealData(comps_product["uzold"]).data() + np_product_old; - } - - int pid_product; -#pragma omp critical (doFieldIonization_nextid) - { - // ID of first newly-created product particle - pid_product = pc_product->NextID(); - // Update NextID to include particles created in this function - pc_product->setNextID(pid_product+np_ionized); - } - const int cpuid = ParallelDescriptor::MyProc(); - - // Loop over all source particles. If is_ionized, copy particle data - // to corresponding product particle. - amrex::For( - np_source, [=] AMREX_GPU_DEVICE (int is) noexcept - { - if(p_is_ionized[is]){ - // offset of 1 due to inclusive scan - int ip = p_i_product[is]-1; - // is: index of ionized particle in source species - // ip: index of corresponding new particle in product species - WarpXParticleContainer::ParticleType& p_product = particles_product[ip]; - WarpXParticleContainer::ParticleType& p_source = particles_source[is]; - // Copy particle from source to product: AoS - p_product.id() = pid_product + ip; - p_product.cpu() = cpuid; - p_product.pos(0) = p_source.pos(0); - p_product.pos(1) = p_source.pos(1); -#if (AMREX_SPACEDIM == 3) - p_product.pos(2) = p_source.pos(2); -#endif - // Copy particle from source to product: SoA - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - attribs_product[ia][ip] = attribs_source[ia][is]; - } - // Update xold etc. if boosted frame diagnostics required - // for product species. Fill runtime attribs with a copy of - // current properties (xold = x etc.). - if (do_boosted_product) { - runtime_attribs_product[0][ip] = p_source.pos(0); - runtime_attribs_product[1][ip] = p_source.pos(1); - runtime_attribs_product[2][ip] = p_source.pos(2); - runtime_attribs_product[3][ip] = runtime_uold_source[0][ip]; - runtime_attribs_product[4][ip] = runtime_uold_source[1][ip]; - runtime_attribs_product[5][ip] = runtime_uold_source[2][ip]; - } - } - } - ); - } -} - void MultiParticleContainer::doFieldIonization () { @@ -727,7 +601,17 @@ MultiParticleContainer::doFieldIonization () amrex::Gpu::ManagedDeviceVector<int> is_ionized; pc_source->buildIonizationMask(mfi, lev, is_ionized); // Create particles in pc_product - createIonizedParticles(lev, mfi, pc_source, pc_product, is_ionized); + int do_boost = WarpX::do_back_transformed_diagnostics + && pc_product->doBackTransformedDiagnostics(); + amrex::Gpu::ManagedDeviceVector<int> v_do_back_transformed_product{do_boost}; + const amrex::Vector<WarpXParticleContainer*> v_pc_product {pc_product.get()}; + // Copy source to product particles, and increase ionization + // level of source particle + ionization_process.createParticles(lev, mfi, pc_source, v_pc_product, + is_ionized, v_do_back_transformed_product); + // Synchronize to prevent the destruction of temporary arrays (at the + // end of the function call) before the kernel executes. + Gpu::streamSynchronize(); } } // lev } // pc_source @@ -736,15 +620,295 @@ MultiParticleContainer::doFieldIonization () #ifdef WARPX_QED void MultiParticleContainer::InitQED () { + m_shr_p_qs_engine = std::make_shared<QuantumSynchrotronEngine>(); + m_shr_p_bw_engine = std::make_shared<BreitWheelerEngine>(); + + m_nspecies_quantum_sync = 0; + m_nspecies_breit_wheeler = 0; + for (auto& pc : allcontainers) { if(pc->has_quantum_sync()){ pc->set_quantum_sync_engine_ptr - (std::make_shared<QuantumSynchrotronEngine>(qs_engine)); + (m_shr_p_qs_engine); + m_nspecies_quantum_sync++; } if(pc->has_breit_wheeler()){ pc->set_breit_wheeler_engine_ptr - (std::make_shared<BreitWheelerEngine>(bw_engine)); + (m_shr_p_bw_engine); + m_nspecies_breit_wheeler++; } } + + if(m_nspecies_quantum_sync != 0) + InitQuantumSync(); + + if(m_nspecies_breit_wheeler !=0) + InitBreitWheeler(); + +} + +void MultiParticleContainer::InitQuantumSync () +{ + bool generate_table; + PicsarQuantumSynchrotronCtrl ctrl; + std::string filename; + std::tie(generate_table, filename, ctrl) = ParseQuantumSyncParams(); + + //Only temporary for test purposes, will be removed + ParmParse pp("qed_qs"); + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) return; + //_________________________________________________ + + + if(generate_table && ParallelDescriptor::IOProcessor()){ + m_shr_p_qs_engine->compute_lookup_tables(ctrl); + Vector<char> all_data = m_shr_p_qs_engine->export_lookup_tables_data(); + WarpXUtilIO::WriteBinaryDataOnFile(filename, all_data); + } + ParallelDescriptor::Barrier(); + + Vector<char> table_data; + ParallelDescriptor::ReadAndBcastFile(filename, table_data); + ParallelDescriptor::Barrier(); + + //No need to initialize from raw data for the processor that + //has just generated the table + if(!generate_table || !ParallelDescriptor::IOProcessor()){ + m_shr_p_qs_engine->init_lookup_tables_from_raw_data(table_data); + } + + if(!m_shr_p_qs_engine->are_lookup_tables_initialized()) + amrex::Error("Table initialization has failed!\n"); +} + +void MultiParticleContainer::InitBreitWheeler () +{ + bool generate_table; + PicsarBreitWheelerCtrl ctrl; + std::string filename; + std::tie(generate_table, filename, ctrl) = ParseBreitWheelerParams(); + + //Only temporary for test purposes, will be removed + ParmParse pp("qed_bw"); + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) return; + //_________________________________________________ + + if(generate_table && ParallelDescriptor::IOProcessor()){ + m_shr_p_bw_engine->compute_lookup_tables(ctrl); + Vector<char> all_data = m_shr_p_bw_engine->export_lookup_tables_data(); + WarpXUtilIO::WriteBinaryDataOnFile(filename, all_data); + } + ParallelDescriptor::Barrier(); + + Vector<char> table_data; + ParallelDescriptor::ReadAndBcastFile(filename, table_data); + ParallelDescriptor::Barrier(); + + //No need to initialize from raw data for the processor that + //has just generated the table + if(!generate_table || !ParallelDescriptor::IOProcessor()){ + m_shr_p_bw_engine->init_lookup_tables_from_raw_data(table_data); + } + + if(!m_shr_p_bw_engine->are_lookup_tables_initialized()) + amrex::Error("Table initialization has failed!\n"); +} + +std::tuple<bool,std::string,PicsarQuantumSynchrotronCtrl> +MultiParticleContainer::ParseQuantumSyncParams () +{ + PicsarQuantumSynchrotronCtrl ctrl = + m_shr_p_qs_engine->get_default_ctrl(); + bool generate_table{false}; + std::string table_name; + + ParmParse pp("qed_qs"); + + // Engine paramenter: chi_part_min is the minium chi parameter to be + // considered by the engine. If a lepton has chi < chi_part_min, + // the optical depth is not evolved and photon generation is ignored + pp.query("chi_min", ctrl.chi_part_min); + + //Only temporary for test purposes, will be removed + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) + return std::make_tuple(false, "__DUMMY__", ctrl); + //_________________________________________________ + + pp.query("generate_table", generate_table); + if(generate_table){ + int t_int = 0; + + //==Table parameters== + + //--- sub-table 1 (1D) + //These parameters are used to pre-compute a function + //which appears in the evolution of the optical depth + + //Minimun chi for the table. If a lepton has chi < chi_part_tdndt_min, + ///chi is considered as it were equal to chi_part_tdndt_min + pp.query("tab_dndt_chi_min", ctrl.chi_part_tdndt_min); + + //Maximum chi for the table. If a lepton has chi > chi_part_tdndt_max, + ///chi is considered as it were equal to chi_part_tdndt_max + pp.query("tab_dndt_chi_max", ctrl.chi_part_tdndt_max); + + //How many points should be used for chi in the table + pp.query("tab_dndt_how_many", t_int); + ctrl.chi_part_tdndt_how_many= t_int; + //------ + + //--- sub-table 2 (2D) + //These parameters are used to pre-compute a function + //which is used to extract the properties of the generated + //photons. + + //Minimun chi for the table. If a lepton has chi < chi_part_tem_min, + ///chi is considered as it were equal to chi_part_tem_min + pp.query("tab_em_chi_min", ctrl.chi_part_tem_min); + + //Maximum chi for the table. If a lepton has chi > chi_part_tem_max, + ///chi is considered as it were equal to chi_part_tem_max + pp.query("tab_em_chi_max", ctrl.chi_part_tem_max); + + //How many points should be used for chi in the table + pp.query("tab_em_chi_how_many", t_int); + ctrl.chi_part_tem_how_many = t_int; + + //The other axis of the table is a cumulative probability distribution + //(corresponding to different energies of the generated particles) + //This parameter is the number of different points to consider + pp.query("tab_em_prob_how_many", t_int); + ctrl.prob_tem_how_many = t_int; + //==================== + + pp.query("save_table_in", table_name); + + } + + std::string load_table_name; + pp.query("load_table_from", load_table_name); + if(!load_table_name.empty()){ + if(generate_table && ParallelDescriptor::IOProcessor()){ + amrex::Print() << "Warning, Quantum Synchrotron table will be loaded, not generated. \n"; + } + table_name = load_table_name; + generate_table = false; + } + +#ifndef WARPX_QED_TABLE_GEN + if(generate_table){ + amrex::Error("Error: Compile with QED_TABLE_GEN=TRUE to enable table generation!\n"); + } +#endif + + if(table_name.empty()){ + amrex::Error("Error: Quantum Synchrotron table has either to be generated or to be loaded.\n"); + } + + return std::make_tuple(generate_table, table_name, ctrl); +} + +std::tuple<bool,std::string,PicsarBreitWheelerCtrl> +MultiParticleContainer::ParseBreitWheelerParams () +{ + PicsarBreitWheelerCtrl ctrl = + m_shr_p_bw_engine->get_default_ctrl(); + bool generate_table{false}; + std::string table_name; + + ParmParse pp("qed_bw"); + + // Engine paramenter: chi_phot_min is the minium chi parameter to be + // considered by the engine. If a photon has chi < chi_phot_min, + // the optical depth is not evolved and pair generation is ignored + pp.query("chi_min", ctrl.chi_phot_min); + + //Only temporary for test purposes, will be removed + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) + return std::make_tuple(false, "__DUMMY__", ctrl); + //_________________________________________________ + + pp.query("generate_table", generate_table); + if(generate_table){ + + int t_int; + + //==Table parameters== + + //--- sub-table 1 (1D) + //These parameters are used to pre-compute a function + //which appears in the evolution of the optical depth + + //Minimun chi for the table. If a photon has chi < chi_phot_tdndt_min, + //an analytical approximation is used. + pp.query("tab_dndt_chi_min", ctrl.chi_phot_tdndt_min); + + //Maximum chi for the table. If a photon has chi > chi_phot_tdndt_min, + //an analytical approximation is used. + pp.query("tab_dndt_chi_max", ctrl.chi_phot_tdndt_max); + + //How many points should be used for chi in the table + pp.query("tab_dndt_how_many", t_int); + ctrl.chi_phot_tdndt_how_many = t_int; + //------ + + //--- sub-table 2 (2D) + //These parameters are used to pre-compute a function + //which is used to extract the properties of the generated + //particles. + + //Minimun chi for the table. If a photon has chi < chi_phot_tpair_min + //chi is considered as it were equal to chi_phot_tpair_min + pp.query("tab_pair_chi_min", ctrl.chi_phot_tpair_min); + + //Maximum chi for the table. If a photon has chi > chi_phot_tpair_max + //chi is considered as it were equal to chi_phot_tpair_max + pp.query("tab_pair_chi_max", ctrl.chi_phot_tpair_max); + + //How many points should be used for chi in the table + pp.query("tab_pair_chi_how_many", t_int); + ctrl.chi_phot_tpair_how_many = t_int; + + //The other axis of the table is the fraction of the initial energy + //'taken away' by the most energetic particle of the pair. + //This parameter is the number of different fractions to consider + pp.query("tab_pair_frac_how_many", t_int); + ctrl.chi_frac_tpair_how_many = t_int; + //==================== + + pp.query("save_table_in", table_name); + } + + std::string load_table_name; + pp.query("load_table_from", load_table_name); + if(!load_table_name.empty()){ + if(generate_table && ParallelDescriptor::IOProcessor()){ + amrex::Print() << "Warning, Breit Wheeler table will be loaded, not generated. \n"; + } + table_name = load_table_name; + generate_table = false; + } + +#ifndef WARPX_QED_TABLE_GEN + if(generate_table){ + if(ParallelDescriptor::IOProcessor()){ + amrex::Error("Error: Compile with QED_TABLE_GEN=TRUE to enable table generation!\n"); + } + } +#endif + + if(table_name.empty()){ + amrex::Error("Error: Breit Wheeler table has either to be generated or to be loaded.\n"); + } + + return std::make_tuple(generate_table, table_name, ctrl); } #endif diff --git a/Source/Particles/ParticleCreation/CopyParticle.H b/Source/Particles/ParticleCreation/CopyParticle.H new file mode 100644 index 000000000..bd2d530fb --- /dev/null +++ b/Source/Particles/ParticleCreation/CopyParticle.H @@ -0,0 +1,90 @@ +#ifndef COPYPARTICLE_H_ +#define COPYPARTICLE_H_ + +#include "WarpXParticleContainer.H" + +/** + * \brief Functor to copy one particle + * + * This is meant to be a very small class captured by value in kernel launches, + * that can be initialized on the host and copied to the device at each + * iteration. It should be general enough to be used by all processes. + * + * Pointers to SoA data are saved when constructor is called. + * AoS data is passed as argument of operator (). + */ +class copyParticle +{ + +public: + + // ID of MPI rank + int m_cpuid; + // If true, will copy old attribs for back-transforme diagnostics + bool m_do_back_transformed_product; + // Source old (runtime) attribs for back-transformed diagnostics + amrex::GpuArray<amrex::ParticleReal*,3> m_runtime_uold_source; + // Source attribs + amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> m_attribs_source; + // Product attribs + amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> m_attribs_product; + // Product runtime attribs + amrex::GpuArray<amrex::ParticleReal*,6> m_runtime_attribs_product; + + // Simple constructor + AMREX_GPU_HOST_DEVICE + copyParticle( + const int cpuid, const int do_back_transformed_product, + const amrex::GpuArray<amrex::ParticleReal*,3> runtime_uold_source, + const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_source, + const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_product, + const amrex::GpuArray<amrex::ParticleReal*,6> runtime_attribs_product) + : + m_cpuid(cpuid), + m_do_back_transformed_product(do_back_transformed_product), + m_runtime_uold_source(runtime_uold_source), + m_attribs_source(attribs_source), + m_attribs_product(attribs_product), + m_runtime_attribs_product(runtime_attribs_product){} + + /** + * \brief Overload operator () to do the copy of 1 particle + * + * \param is index of source particle + * \param ip index of product particle + * \param pid_product ID of first product particle + * \param p_source Struct with data for 1 source particle (position etc.) + * \param p_source Struct with data for 1 product particle (position etc.) + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator () (int is, int ip, int pid_product, + WarpXParticleContainer::ParticleType& p_source, + WarpXParticleContainer::ParticleType& p_product) const noexcept + { + // Copy particle from source to product: AoS + p_product.id() = pid_product + ip; + p_product.cpu() = m_cpuid; + p_product.pos(0) = p_source.pos(0); + p_product.pos(1) = p_source.pos(1); +#if (AMREX_SPACEDIM == 3) + p_product.pos(2) = p_source.pos(2); +#endif + // Copy particle from source to product: SoA + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + m_attribs_product[ia][ip] = m_attribs_source[ia][is]; + } + // Update xold etc. if boosted frame diagnostics required + // for product species. Fill runtime attribs with a copy of + // current properties (xold = x etc.). + if (m_do_back_transformed_product) { + m_runtime_attribs_product[0][ip] = p_source.pos(0); + m_runtime_attribs_product[1][ip] = p_source.pos(1); + m_runtime_attribs_product[2][ip] = p_source.pos(2); + m_runtime_attribs_product[3][ip] = m_runtime_uold_source[0][ip]; + m_runtime_attribs_product[4][ip] = m_runtime_uold_source[1][ip]; + m_runtime_attribs_product[5][ip] = m_runtime_uold_source[2][ip]; + } + } +}; + +#endif // COPYPARTICLE_H_ diff --git a/Source/Particles/ParticleCreation/ElementaryProcess.H b/Source/Particles/ParticleCreation/ElementaryProcess.H new file mode 100644 index 000000000..fa791c244 --- /dev/null +++ b/Source/Particles/ParticleCreation/ElementaryProcess.H @@ -0,0 +1,257 @@ +#ifndef ELEMENTARYPROCESS_H_ +#define ELEMENTARYPROCESS_H_ + +#include "WarpXParticleContainer.H" +#include "CopyParticle.H" +#include "TransformParticle.H" + +/** + * \brief Base class for particle creation processes + * + * Particles in a source particle container are copied to product particle + * containers if they are flagged. Both source and product particles can be + * modified. + * + * - initialize_copy initializes a general copy functor + * - createParticles formats the data to prepare for the copy, then + * calls copyAndTransformParticles + * - copyAndTransformParticles loops over particles, performs the copy and + * transform source and product particles if needed. + * + * The class is templated on the process type. This gives flexibility + * on source and product particle transformations. + */ +template<elementaryProcessType ProcessT> +class elementaryProcess +{ +public: + + /** + * \brief initialize and return functor for particle copy + * + * \param cpuid id of MPI rank + * \param do_back_transformed_product whether to copy old attribs + * \param runtime_uold_source momentum of source particles + * \param attribs_source attribs of source particles + * \param attribs_product attribs of product particles + * \param runtime_attribs_product runtime attribs for product, to store old attribs + */ + copyParticle initialize_copy( + const int cpuid, const int do_back_transformed_product, + const amrex::GpuArray<amrex::ParticleReal*,3> runtime_uold_source, + const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_source, + const amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_product, + const amrex::GpuArray<amrex::ParticleReal*,6> runtime_attribs_product) const noexcept + { + return copyParticle ( + cpuid, do_back_transformed_product, runtime_uold_source, + attribs_source, attribs_product, runtime_attribs_product); + }; + + /** + * \brief create particles in product particle containers + * + * For particle i in mfi, if is_flagged[i]=1, copy particle + * particle i from container pc_source into each container in v_pc_product + * + * \param lev MR level + * \param mfi MFIter where source particles are located + * \param pc_source source particle container + * \param v_pc_product vector of product particle containers + * \param is_flagged particles for which is_flagged is 1 are copied + * \param v_do_back_transformed_product vector: whether to copy old attribs for + * each product container + */ + void createParticles ( + int lev, const amrex::MFIter& mfi, + std::unique_ptr< WarpXParticleContainer>& pc_source, + amrex::Vector<WarpXParticleContainer*> v_pc_product, + amrex::Gpu::ManagedDeviceVector<int>& is_flagged, + amrex::Gpu::ManagedDeviceVector<int> v_do_back_transformed_product) + { + + BL_PROFILE("createParticles"); + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + // Get source particle data + auto& ptile_source = pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + const int np_source = ptile_source.GetArrayOfStructs().size(); + if (np_source == 0) return; + // --- source AoS particle data + WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); + // --- source SoA particle data + auto& soa_source = ptile_source.GetStructOfArrays(); + amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_source; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + attribs_source[ia] = soa_source.GetRealData(ia).data(); + } + // --- source runtime attribs + amrex::GpuArray<amrex::ParticleReal*,3> runtime_uold_source; + // Prepare arrays for boosted frame diagnostics. + runtime_uold_source[0] = soa_source.GetRealData(PIdx::ux).data(); + runtime_uold_source[1] = soa_source.GetRealData(PIdx::uy).data(); + runtime_uold_source[2] = soa_source.GetRealData(PIdx::uz).data(); + + // --- source runtime integer attribs + amrex::GpuArray<int*,1> runtime_iattribs_source; + std::map<std::string, int> icomps_source = pc_source->getParticleiComps(); + runtime_iattribs_source[0] = soa_source.GetIntData(icomps_source["ionization_level"]).data(); + + int nproducts = v_pc_product.size(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + v_do_back_transformed_product.size() == nproducts, + "v_pc_product and v_do_back_transformed_product must have the same size."); + + // Indices of product particle for each source particle. + // i_product[i]-1 is the location in product tile of product particle + // from source particle i. + amrex::Gpu::ManagedDeviceVector<int> i_product; + i_product.resize(np_source); + // 0<i<np_source + // 0<i_product<np_flagged + // Strictly speaking, i_product should be an exclusive_scan of + // is_flagged. However, for indices where is_flagged is 1, the + // inclusive scan gives the same result with an offset of 1. + // The advantage of inclusive_scan is that the sum of is_flagged + // is in the last element, so no other reduction is required to get + // number of particles. + // Gpu::inclusive_scan runs on the current GPU stream, and synchronizes + // with the CPU, so that the next line (executed by the CPU) has the + // updated values of i_product + amrex::Gpu::inclusive_scan(is_flagged.begin(), is_flagged.end(), i_product.begin()); + int np_flagged = i_product[np_source-1]; + if (np_flagged == 0) return; + + amrex::Vector<copyParticle> v_copy_functor; + amrex::Gpu::ManagedDeviceVector<int> v_pid_product(nproducts); + amrex::Gpu::ManagedDeviceVector<WarpXParticleContainer::ParticleType*> v_particles_product(nproducts); + for (int iproduct=0; iproduct<nproducts; iproduct++){ + WarpXParticleContainer*& pc_product = v_pc_product[iproduct]; + // Get product particle data + auto& ptile_product = pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + // old and new (i.e., including flagged particles) number of particles + // for product species + const int np_product_old = ptile_product.GetArrayOfStructs().size(); + const int np_product_new = np_product_old + np_flagged; + // Allocate extra space in product species for flagged particles. + ptile_product.resize(np_product_new); + // --- product AoS particle data + // WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; + v_particles_product[iproduct] = ptile_product.GetArrayOfStructs()().data() + np_product_old; + // First element is the first newly-created product particle + // --- product SoA particle data + auto& soa_product = ptile_product.GetStructOfArrays(); + amrex::GpuArray<amrex::ParticleReal*,PIdx::nattribs> attribs_product; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + // First element is the first newly-created product particle + attribs_product[ia] = soa_product.GetRealData(ia).data() + np_product_old; + } + // --- product runtime attribs + amrex::GpuArray<amrex::ParticleReal*,6> runtime_attribs_product; + const int do_boost = v_do_back_transformed_product[iproduct]; + if (do_boost) { + std::map<std::string, int> comps_product = pc_product->getParticleComps(); + runtime_attribs_product[0] = soa_product.GetRealData(comps_product[ "xold"]).data() + np_product_old; + runtime_attribs_product[1] = soa_product.GetRealData(comps_product[ "yold"]).data() + np_product_old; + runtime_attribs_product[2] = soa_product.GetRealData(comps_product[ "zold"]).data() + np_product_old; + runtime_attribs_product[3] = soa_product.GetRealData(comps_product["uxold"]).data() + np_product_old; + runtime_attribs_product[4] = soa_product.GetRealData(comps_product["uyold"]).data() + np_product_old; + runtime_attribs_product[5] = soa_product.GetRealData(comps_product["uzold"]).data() + np_product_old; + } + + // --- product runtime integer attribs + int pid_product; +#pragma omp critical (doParticleCreation_nextid) + { + // ID of first newly-created product particle + pid_product = pc_product->NextID(); + // Update NextID to include particles created in this function + pc_product->setNextID(pid_product+np_flagged); + } + const int cpuid = amrex::ParallelDescriptor::MyProc(); + + // Create instance of copy functor, and add it to the vector + v_copy_functor.push_back (initialize_copy( + cpuid, v_do_back_transformed_product[iproduct], + runtime_uold_source, + attribs_source, + attribs_product, + runtime_attribs_product) ); + v_pid_product[iproduct] = pid_product; + + } + // Loop over source particles and create product particles + copyAndTransformParticles(is_flagged, i_product, np_source, v_pid_product, + v_particles_product, particles_source, v_copy_functor, + runtime_iattribs_source); + // Synchronize to prevent the destruction of temporary arrays (at the + // end of the function call) before the kernel executes. + amrex::Gpu::streamSynchronize(); + } + + /** + * \brief Loop over source particles and create product particles + * + * \param is_flagged particles for which is_flagged is 1 are copied + * \param i_product relative indices of product particle. This is the same + * for all product particle containers. + * \param np_source number of particles in source container + * \param v_pid_product for each product particle container: ID of the + * first product particle. Others IDs are incremented from this one + * \param v_particles_product for each product particle container: pointer + * to the AoS particle data + * \param particles_source pointter to the AoS source particle data + * \param v_copy_functor vector of copy functors, one for each product particle container + * \param runtime_iattribs_source pointer to source runtime integer attribs + */ + void copyAndTransformParticles( + amrex::Gpu::ManagedDeviceVector<int>& is_flagged, + amrex::Gpu::ManagedDeviceVector<int>& i_product, + int np_source, + amrex::Gpu::ManagedDeviceVector<int> v_pid_product, + amrex::Gpu::ManagedDeviceVector<WarpXParticleContainer::ParticleType*> v_particles_product, + WarpXParticleContainer::ParticleType* particles_source, + const amrex::Vector<copyParticle>& v_copy_functor, + amrex::GpuArray<int*,1> runtime_iattribs_source) + { + int const * const AMREX_RESTRICT p_is_flagged = is_flagged.dataPtr(); + int const * const AMREX_RESTRICT p_i_product = i_product.dataPtr(); + copyParticle const * const p_copy_functor = v_copy_functor.data(); + WarpXParticleContainer::ParticleType * * p_particles_product = v_particles_product.data(); + int const * const p_pid_product = v_pid_product.data(); + + // Loop over all source particles. If is_flagged, copy particle data + // to corresponding product particle. + int nproducts = v_particles_product.size(); + amrex::For( + np_source, [=] AMREX_GPU_DEVICE (int is) noexcept + { + if(p_is_flagged[is]){ + // offset of 1 due to inclusive scan + int ip = p_i_product[is]-1; + WarpXParticleContainer::ParticleType& p_source = particles_source[is]; + for (int iproduct=0; iproduct<nproducts; iproduct++){ + // is: index of flagged particle in source species + // ip: index of corresponding new particle in product species + WarpXParticleContainer::ParticleType* particles_product = p_particles_product[iproduct]; + WarpXParticleContainer::ParticleType& p_product = particles_product[ip]; + p_copy_functor[iproduct](is, ip, p_pid_product[iproduct], p_source, p_product); + } + // Transform source particle + transformSourceParticle<ProcessT>(is, p_source, runtime_iattribs_source); + // Transform product particles. The loop over product + // species is done inside the function to allow for + // more flexibility. + transformProductParticle<ProcessT>(ip, p_particles_product); + } + } + ); + } +}; + +// Derived class for ionization process +class IonizationProcess: public elementaryProcess<elementaryProcessType::Ionization> +{}; + +#endif // ELEMENTARYPROCESS_H_ diff --git a/Source/Particles/ParticleCreation/Make.package b/Source/Particles/ParticleCreation/Make.package new file mode 100644 index 000000000..6e32f4a77 --- /dev/null +++ b/Source/Particles/ParticleCreation/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += ElementaryProcess.H +CEXE_headers += CopyParticle.H +CEXE_headers += TransformParticle.H + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/ParticleCreation/ +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/ParticleCreation/ diff --git a/Source/Particles/ParticleCreation/TransformParticle.H b/Source/Particles/ParticleCreation/TransformParticle.H new file mode 100644 index 000000000..14e534d0e --- /dev/null +++ b/Source/Particles/ParticleCreation/TransformParticle.H @@ -0,0 +1,44 @@ +#ifndef TRANSFORMPARTICLE_H_ +#define TRANSFORMPARTICLE_H_ + +#include "WarpXParticleContainer.H" + +enum struct elementaryProcessType { Ionization }; + +/** + * \brief small modifications on source particle + * \param i index of particle + * \param particle particle struct + * \param runtime_iattribs integer attribs + */ +template < elementaryProcessType ProcessT > +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void transformSourceParticle( + int i, + WarpXParticleContainer::ParticleType& particle, + amrex::GpuArray<int*,1> runtime_iattribs){} + +/** + * \brief small modifications on product particle + * \param i index of particle + * \param particle particle struct + */ +template < elementaryProcessType ProcessT > +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void transformProductParticle( + int i, + WarpXParticleContainer::ParticleType* * v_particle){} + +// For ionization, increase ionization level of source particle +template <> +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void transformSourceParticle < elementaryProcessType::Ionization > ( + int i, + WarpXParticleContainer::ParticleType& particle, + amrex::GpuArray<int*,1> runtime_iattribs) +{ + // increment particle's ionization level + runtime_iattribs[0][i] += 1; +} + +#endif // TRANSFORMPARTICLE_H_ diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 3e678fb26..7e52b52e1 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -26,9 +26,9 @@ PhotonParticleContainer::PhotonParticleContainer (AmrCore* amr_core, int ispecie ParmParse pp(species_name); #ifdef WARPX_QED - //IF do_qed is enabled, find out if Breit Wheeler process is enabled - if(do_qed) - pp.query("do_qed_breit_wheeler", do_qed_breit_wheeler); + //IF m_do_qed is enabled, find out if Breit Wheeler process is enabled + if(m_do_qed) + pp.query("do_qed_breit_wheeler", m_do_qed_breit_wheeler); //Check for processes which do not make sense for photons bool test_quantum_sync = false; @@ -45,9 +45,8 @@ void PhotonParticleContainer::InitData() { AddParticles(0); // Note - add on level 0 - if (maxLevel() > 0) { - Redistribute(); // We then redistribute - } + Redistribute(); // We then redistribute + } void @@ -74,7 +73,7 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) { copy_attribs(pti, x, y, z); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b18c9b5f8..17a504719 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -10,6 +10,7 @@ #ifdef WARPX_QED #include <QuantumSyncEngineWrapper.H> #include <BreitWheelerEngineWrapper.H> + #include <QedChiFunctions.H> #endif #include <map> @@ -219,16 +220,16 @@ protected: #ifdef WARPX_QED // A flag to enable quantum_synchrotron process for leptons - bool do_qed_quantum_sync = false; + bool m_do_qed_quantum_sync = false; // A flag to enable breit_wheeler process [photons only!!] - bool do_qed_breit_wheeler = false; + bool m_do_qed_breit_wheeler = false; // A smart pointer to an instance of a Quantum Synchrotron engine - std::shared_ptr<QuantumSynchrotronEngine> shr_ptr_qs_engine; + std::shared_ptr<QuantumSynchrotronEngine> m_shr_p_qs_engine; // A smart pointer to an instance of a Breit Wheeler engine [photons only!] - std::shared_ptr<BreitWheelerEngine> shr_ptr_bw_engine; + std::shared_ptr<BreitWheelerEngine> m_shr_p_bw_engine; #endif }; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 785454e09..a1cd703b0 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -43,7 +43,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_continuous_injection", do_continuous_injection); // Whether to plot back-transformed (lab-frame) diagnostics // for this species. - pp.query("do_boosted_frame_diags", do_boosted_frame_diags); + pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics); pp.query("do_field_ionization", do_field_ionization); @@ -62,35 +62,36 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp "Radiation reaction can be enabled only if Boris pusher is used"); //_____________________________ - #ifdef AMREX_USE_GPU - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - do_field_ionization == 0, - "Field ionization does not work on GPU so far, because the current " - "version of Redistribute in AMReX does not work with runtime parameters"); + Print()<<"\n-----------------------------------------------------\n"; + Print()<<"WARNING: field ionization on GPU uses RedistributeCPU\n"; + Print()<<"-----------------------------------------------------\n\n"; + //AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + //do_field_ionization == 0, + //"Field ionization does not work on GPU so far, because the current " + //"version of Redistribute in AMReX does not work with runtime parameters"); #endif - #ifdef WARPX_QED //Add real component if QED is enabled - pp.query("do_qed", do_qed); - if(do_qed) + pp.query("do_qed", m_do_qed); + if(m_do_qed) AddRealComp("tau"); //IF do_qed is enabled, find out if Quantum Synchrotron process is enabled - if(do_qed) - pp.query("do_qed_quantum_sync", do_qed_quantum_sync); + if(m_do_qed) + pp.query("do_qed_quantum_sync", m_do_qed_quantum_sync); //TODO: SHOULD CHECK IF SPECIES IS EITHER ELECTRONS OR POSITRONS!! #endif //variable to set plot_flags size int plot_flag_size = PIdx::nattribs; - if(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if(WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) plot_flag_size += 6; #ifdef WARPX_QED - if(do_qed){ + if(m_do_qed){ // plot_flag will have an entry for the optical depth plot_flag_size++; } @@ -122,7 +123,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } #ifdef WARPX_QED - if(do_qed){ + if(m_do_qed){ //Optical depths is always plotted if QED is on plot_flags[plot_flag_size-1] = 1; } @@ -441,13 +442,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int dir=0; dir<AMREX_SPACEDIM; dir++) { if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); - overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, Real(0.)) * dx[dir]); + overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]); } else { no_overlap = true; break; } if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, Real(0.)) * dx[dir]); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); } else { no_overlap = true; break; } @@ -543,11 +544,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) BreitWheelerGetOpticalDepth breit_wheeler_get_opt; if(loc_has_quantum_sync){ quantum_sync_get_opt = - shr_ptr_qs_engine->build_optical_depth_functor(); + m_shr_p_qs_engine->build_optical_depth_functor(); } if(loc_has_breit_wheeler){ breit_wheeler_get_opt = - shr_ptr_bw_engine->build_optical_depth_functor(); + m_shr_p_bw_engine->build_optical_depth_functor(); } #endif @@ -1012,12 +1013,12 @@ PhysicalParticleContainer::FieldGather (int lev, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,WarpX::E_external[0]); - Eyp.assign(np,WarpX::E_external[1]); - Ezp.assign(np,WarpX::E_external[2]); - Bxp.assign(np,WarpX::B_external[0]); - Byp.assign(np,WarpX::B_external[1]); - Bzp.assign(np,WarpX::B_external[2]); + 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 @@ -1078,7 +1079,7 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) { for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -1148,13 +1149,13 @@ PhysicalParticleContainer::Evolve (int lev, exfab, eyfab, ezfab, bxfab, byfab, bzfab); } - Exp.assign(np,WarpX::E_external[0]); - Eyp.assign(np,WarpX::E_external[1]); - Ezp.assign(np,WarpX::E_external[2]); + 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[0]); - Byp.assign(np,WarpX::B_external[1]); - Bzp.assign(np,WarpX::B_external[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]); // Determine which particles deposit/gather in the buffer, and // which particles deposit/gather in the fine patch @@ -1427,9 +1428,9 @@ PhysicalParticleContainer::SplitParticles(int lev) // before splitting results in a uniform distribution after splitting const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); - amrex::Vector<amrex::Real> split_offset = {dx[0]/2./ppc_nd[0], - dx[1]/2./ppc_nd[1], - dx[2]/2./ppc_nd[2]}; + amrex::Vector<Real> split_offset = {dx[0]/2._rt/ppc_nd[0], + dx[1]/2._rt/ppc_nd[1], + dx[2]/2._rt/ppc_nd[2]}; // particle Array Of Structs data auto& particles = pti.GetArrayOfStructs(); @@ -1585,7 +1586,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf)) + if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf)) { copy_attribs(pti, x, y, z); } @@ -1701,13 +1702,13 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,WarpX::E_external[0]); - Eyp.assign(np,WarpX::E_external[1]); - Ezp.assign(np,WarpX::E_external[2]); + 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[0]); - Byp.assign(np,WarpX::B_external[1]); - Bzp.assign(np,WarpX::B_external[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 @@ -1842,7 +1843,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real // Note the the slice should always move in the negative boost direction. AMREX_ALWAYS_ASSERT(z_new < z_old); - AMREX_ALWAYS_ASSERT(do_boosted_frame_diags == 1); + AMREX_ALWAYS_ASSERT(do_back_transformed_diagnostics == 1); const int nlevs = std::max(0, finestLevel()+1); @@ -2294,8 +2295,6 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const Real p = 1. - std::exp( - w_dtau ); if (random_draw < p){ - // increment particle's ionization level - ion_lev[i] += 1; // update mask p_ionization_mask[i] = 1; } @@ -2315,25 +2314,25 @@ PhysicalParticleContainer::AmIALepton(){ bool PhysicalParticleContainer::has_quantum_sync() { - return do_qed_quantum_sync; + return m_do_qed_quantum_sync; } bool PhysicalParticleContainer::has_breit_wheeler() { - return do_qed_breit_wheeler; + return m_do_qed_breit_wheeler; } void PhysicalParticleContainer:: set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine> ptr) { - shr_ptr_bw_engine = ptr; + m_shr_p_bw_engine = ptr; } void PhysicalParticleContainer:: set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine> ptr) { - shr_ptr_qs_engine = ptr; + m_shr_p_qs_engine = ptr; } #endif diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 535ffec6f..507206041 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -392,12 +392,12 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,WarpX::E_external[0]); - Eyp.assign(np,WarpX::E_external[1]); - Ezp.assign(np,WarpX::E_external[2]); - Bxp.assign(np,WarpX::B_external[0]); - Byp.assign(np,WarpX::B_external[1]); - Bzp.assign(np,WarpX::B_external[2]); + 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 diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 879f4782e..dbd913c5b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -273,13 +273,14 @@ public: AddIntComp(comm); } - int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + int doBackTransformedDiagnostics () const { return do_back_transformed_diagnostics; } virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) {}; std::map<std::string, int> getParticleComps () { return particle_comps;} + std::map<std::string, int> getParticleiComps () { return particle_icomps;} protected: @@ -316,10 +317,10 @@ protected: amrex::Gpu::ManagedVector<amrex::Real> adk_exp_prefactor; std::string physical_element; - int do_boosted_frame_diags = 1; + int do_back_transformed_diagnostics = 1; #ifdef WARPX_QED - bool do_qed = false; + bool m_do_qed = false; virtual bool has_quantum_sync(){return false;}; virtual bool has_breit_wheeler(){return false;}; diff --git a/Source/QED/BreitWheelerEngineInnards.H b/Source/QED/BreitWheelerEngineInnards.H new file mode 100644 index 000000000..640cdfa94 --- /dev/null +++ b/Source/QED/BreitWheelerEngineInnards.H @@ -0,0 +1,46 @@ +#ifndef WARPX_breit_wheeler_engine_innards_h_ +#define WARPX_breit_wheeler_engine_innards_h_ + +#include "QedWrapperCommons.H" + +#include <AMReX_Gpu.H> + +//This includes only the definition of a simple datastructure +//used to control the Breit Wheeler engine. +#include <breit_wheeler_engine_ctrl.h> + +/** + * This structure holds all the parameters required to use the + * Breit Wheeler engine: a POD control structure and lookup + * tables data. + */ +struct BreitWheelerEngineInnards +{ + // Control parameters (a POD struct) + // ctrl contains several parameters: + // - chi_phot_min : the minium chi parameter to be + // considered by the engine + // - chi_phot_tdndt_min : minimun chi for sub-table 1 (1D) + // - chi_phot_tdndt_max : maximum chi for sub-table 1 (1D) + // - chi_phot_tdndt_how_many : how many points to use for sub-table 1 (1D) + // - chi_phot_tpair_min : minimun chi for sub-table 2 (1D) + // - chi_phot_tpair_max : maximum chi for sub-table 2 (1D) + // - chi_phot_tpair_how_many : how many points to use for chi for sub-table 2 (2D) + // - chi_frac_tpair_how_many : how many points to use for the second axis of sub-table 2 (2D) + picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real> ctrl; + + //Lookup table data + //---sub-table 1 (1D) + amrex::Gpu::ManagedVector<amrex::Real> TTfunc_coords; + amrex::Gpu::ManagedVector<amrex::Real> TTfunc_data; + //--- + + //---sub-table 2 (2D) + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_1; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_2; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_data; + //______ +}; +//========================================================== + +#endif //WARPX_breit_wheeler_engine_innards_h_ diff --git a/Source/QED/BreitWheelerEngineTableBuilder.H b/Source/QED/BreitWheelerEngineTableBuilder.H new file mode 100644 index 000000000..e30b82208 --- /dev/null +++ b/Source/QED/BreitWheelerEngineTableBuilder.H @@ -0,0 +1,26 @@ +#ifndef WARPX_breit_wheeler_engine_table_builder_h_ +#define WARPX_breit_wheeler_engine_table_builder_h_ + +#include "QedWrapperCommons.H" +#include "BreitWheelerEngineInnards.H" + +//This includes only the definition of a simple datastructure +//used to control the Breit Wheeler engine. +#include <breit_wheeler_engine_ctrl.h> + +/** + * A class which computes the lookup tables for the Breit Wheeler engine. + */ +class BreitWheelerEngineTableBuilder{ + public: + /** + * Computes the tables. + * @param[in] ctrl control parameters to generate the tables + * @param[out] innards structure holding both a copy of ctrl and lookup tables data + */ + void compute_table + (picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real> ctrl, + BreitWheelerEngineInnards& innards) const; +}; + +#endif //WARPX_breit_wheeler_engine_table_builder_h_ diff --git a/Source/QED/BreitWheelerEngineTableBuilder.cpp b/Source/QED/BreitWheelerEngineTableBuilder.cpp new file mode 100644 index 000000000..3326d5b59 --- /dev/null +++ b/Source/QED/BreitWheelerEngineTableBuilder.cpp @@ -0,0 +1,58 @@ +#include "BreitWheelerEngineTableBuilder.H" + +//Include the full Breit Wheeler engine with table generation support +//(after some consistency tests). This requires to have a recent version +// of the Boost library. +#ifdef PXRMP_CORE_ONLY + #error The Table Builder is incompatible with PXRMP_CORE_ONLY +#endif + +#ifdef __PICSAR_MULTIPHYSICS_BREIT_WHEELER_ENGINE__ + #warning breit_wheeler_engine.hpp should not have been included before reaching this point. +#endif +#include <breit_wheeler_engine.hpp> +//_______________________________________________ + +//Some handy aliases +using PicsarBreitWheelerEngine = picsar::multi_physics:: + breit_wheeler_engine<amrex::Real, QedUtils::DummyStruct>; + +using PicsarBreitWheelerCtrl = + picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real>; +//_______________________________________________ + +void +BreitWheelerEngineTableBuilder::compute_table + (PicsarBreitWheelerCtrl ctrl, + BreitWheelerEngineInnards& innards) const +{ + PicsarBreitWheelerEngine bw_engine( + std::move(QedUtils::DummyStruct()), 1.0, ctrl); + + bw_engine.compute_dN_dt_lookup_table(); + bw_engine.compute_cumulative_pair_table(); + + auto bw_innards_picsar = bw_engine.export_innards(); + + //Copy data in a GPU-friendly data-structure + innards.ctrl = bw_innards_picsar.bw_ctrl; + innards.TTfunc_coords.assign(bw_innards_picsar.TTfunc_table_coords_ptr, + bw_innards_picsar.TTfunc_table_coords_ptr + + bw_innards_picsar.TTfunc_table_coords_how_many); + innards.TTfunc_data.assign(bw_innards_picsar.TTfunc_table_data_ptr, + bw_innards_picsar.TTfunc_table_data_ptr + + bw_innards_picsar.TTfunc_table_data_how_many); + innards.cum_distrib_coords_1.assign( + bw_innards_picsar.cum_distrib_table_coords_1_ptr, + bw_innards_picsar.cum_distrib_table_coords_1_ptr + + bw_innards_picsar.cum_distrib_table_coords_1_how_many); + innards.cum_distrib_coords_2.assign( + bw_innards_picsar.cum_distrib_table_coords_2_ptr, + bw_innards_picsar.cum_distrib_table_coords_2_ptr + + bw_innards_picsar.cum_distrib_table_coords_2_how_many); + innards.cum_distrib_data.assign( + bw_innards_picsar.cum_distrib_table_data_ptr, + bw_innards_picsar.cum_distrib_table_data_ptr + + bw_innards_picsar.cum_distrib_table_data_how_many); + //____ +} diff --git a/Source/QED/BreitWheelerEngineWrapper.H b/Source/QED/BreitWheelerEngineWrapper.H index 7cdc66b44..1033ff7c9 100644 --- a/Source/QED/BreitWheelerEngineWrapper.H +++ b/Source/QED/BreitWheelerEngineWrapper.H @@ -2,53 +2,302 @@ #define WARPX_breit_wheeler_engine_wrapper_h_ #include "QedWrapperCommons.H" +#include "BreitWheelerEngineInnards.H" -//BW ENGINE from PICSAR +#include <AMReX_Array.H> +#include <AMReX_Vector.H> +#include <AMReX_Gpu.H> + +//#define PXRMP_CORE_ONLY allows importing only the 'core functions' of the +//Breit Wheeler engine of the QED PICSAR library. #define PXRMP_CORE_ONLY -#include "breit_wheeler_engine.hpp" +#include <breit_wheeler_engine.hpp> + +//Lookup table building function is in a dedicated (optional) class to +//avoid including heavy dependencies if they are not needed. +#ifdef WARPX_QED_TABLE_GEN + #include "BreitWheelerEngineTableBuilder.H" +#endif + +#include <string> + +//Some handy aliases + +// The engine has two templated arguments: the numerical type +// and a random number generator. However, random numbers are not +// used to generate the lookup tables and the static member +// functions which are called in the functors do not use +// random numbers as well. Therefore, an empty "DummyStruct" +// can be passed as a template parameter. +using PicsarBreitWheelerEngine = picsar::multi_physics:: + breit_wheeler_engine<amrex::Real, QedUtils::DummyStruct>; -using WarpXBreitWheelerWrapper = - picsar::multi_physics::breit_wheeler_engine<amrex::Real, DummyStruct>; +using PicsarBreitWheelerCtrl = + picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real>; +//__________ // Functors ================================== -// These functors provide the core elementary functions of the library -// Can be included in GPU kernels +// These functors allow using the core elementary functions of the library. +// They are generated by a factory class (BreitWheelerEngine, see below). +// They can be included in GPU kernels. /** - * \brief Functor to initialize the optical depth of photons for the + * Functor to initialize the optical depth of photons for the * Breit-Wheeler process */ class BreitWheelerGetOpticalDepth { public: + /** + * Constructor does nothing because optical depth initialization + * does not require control parameters or lookup tables. + */ BreitWheelerGetOpticalDepth () {}; + /** + * () operator is just a thin wrapper around a very simple function to + * generate the optical depth. It can be used on GPU. + */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE - amrex::Real operator() () const + amrex::Real operator() () const noexcept { - return WarpXBreitWheelerWrapper:: + //A random number in [0,1) should be provided as an argument. + return PicsarBreitWheelerEngine:: internal_get_optical_depth(amrex::Random()); } }; //____________________________________________ +/** + * Functor to evolve the optical depth of photons due to the + * Breit-Wheeler process + */ +class BreitWheelerEvolveOpticalDepth +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + BreitWheelerEvolveOpticalDepth(BreitWheelerEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_TTfunc_size{r_innards.TTfunc_coords.size()}, + m_p_TTfunc_coords{r_innards.TTfunc_coords.dataPtr()}, + m_p_TTfunc_data{r_innards.TTfunc_data.dataPtr()} + {}; + + /** + * Evolves the optical depth. It can be used on GPU. + * @param[in] px,py,pz momentum components of the photon (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] dt timestep (SI units) + * @param[in,out] opt_depth optical depth of the photon. It is modified by the method. + * @return a flag which is 1 if optical depth becomes negative (i.e. a pair has to be generated). + */ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + int operator()( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real dt, amrex::Real& opt_depth) const noexcept + { + bool has_event_happened{false}; + + //the library provides the time (< dt) at which the event occurs, but this + //feature won't be used in WarpX for now. + amrex::Real unused_event_time{0.0}; + + PicsarBreitWheelerEngine:: + internal_evolve_opt_depth_and_determine_event( + px, py, pz, + ex, ey, ez, + bx, by, bz, + dt, opt_depth, + has_event_happened, unused_event_time, + m_dummy_lambda, + picsar::multi_physics::lookup_1d<amrex::Real>{ + m_TTfunc_size, + m_p_TTfunc_coords, + m_p_TTfunc_data}, + m_ctrl); + + return has_event_happened; + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarBreitWheelerCtrl m_ctrl; + + //lookup table data + size_t m_TTfunc_size; + amrex::Real* m_p_TTfunc_coords; + amrex::Real* m_p_TTfunc_data; +}; + +/** + * Functor to generate a pair via the + * Breit-Wheeler process + */ +class BreitWheelerGeneratePairs +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + BreitWheelerGeneratePairs( + BreitWheelerEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_cum_distrib_coords_1_size{r_innards.cum_distrib_coords_1.size()}, + m_cum_distrib_coords_2_size{r_innards.cum_distrib_coords_2.size()}, + m_p_distrib_coords_1{r_innards.cum_distrib_coords_1.data()}, + m_p_distrib_coords_2{r_innards.cum_distrib_coords_2.data()}, + m_p_cum_distrib_data{r_innards.cum_distrib_data.data() + }{}; + + /** + * Generates sampling (template parameter) pairs according to Breit Wheeler process. + * It can be used on GPU. + * @param[in] px,py,pz momentum components of the photon (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] weight of the photon (code units) + * @param[out] e_px,e_py,e_pz momenta of generated electrons. Each array should have size=sampling (SI units) + * @param[out] p_px,p_py,p_pz momenta of generated positrons. Each array should have size=sampling (SI units) + * @param[out] e_weight,p_weight weight of the generated particles Each array should have size=sampling (code units). + */ + template <size_t sampling> + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void operator()( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real weight, + amrex::Real* e_px, amrex::Real* e_py, amrex::Real* e_pz, + amrex::Real* p_px, amrex::Real* p_py, amrex::Real* p_pz, + amrex::Real* e_weight, amrex::Real* p_weight) const noexcept + { + //[sampling] random numbers are needed + picsar::multi_physics::picsar_array<amrex::Real, sampling> + rand_zero_one_minus_epsi; + for(auto& el : rand_zero_one_minus_epsi) el = amrex::Random(); + + PicsarBreitWheelerEngine:: + internal_generate_breit_wheeler_pairs( + px, py, pz, + ex, ey, ez, + bx, by, bz, + weight, sampling, + e_px, e_py, e_pz, + p_px, p_py, p_pz, + e_weight, p_weight, + m_dummy_lambda, + picsar::multi_physics::lookup_2d<amrex::Real>{ + m_cum_distrib_coords_1_size, + m_p_distrib_coords_1, + m_cum_distrib_coords_2_size, + m_p_distrib_coords_2, + m_p_cum_distrib_data + }, + m_ctrl, + rand_zero_one_minus_epsi.data()); + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarBreitWheelerCtrl m_ctrl; + + //lookup table data + size_t m_cum_distrib_coords_1_size; + size_t m_cum_distrib_coords_2_size; + amrex::Real* m_p_distrib_coords_1; + amrex::Real* m_p_distrib_coords_2; + amrex::Real* m_p_cum_distrib_data; +}; + // Factory class ============================= /** - * \brief Wrapper for the Breit Wheeler engine of the PICSAR library + * Wrapper for the Breit Wheeler engine of the PICSAR library */ class BreitWheelerEngine { public: + /** + * Constructor requires no arguments. + */ BreitWheelerEngine (); /** - * \brief Builds the functor to initialize the optical depth + * Builds the functor to initialize the optical depth */ BreitWheelerGetOpticalDepth build_optical_depth_functor (); + + /** + * Builds the functor to evolve the optical depth + */ + BreitWheelerEvolveOpticalDepth build_evolve_functor (); + + /** + * Builds the functor to generate the pairs + */ + BreitWheelerGeneratePairs build_pair_functor (); + + /** + * Checks if the optical tables are properly initialized + */ + bool are_lookup_tables_initialized () const; + + /** + * Init lookup tables from raw binary data. + * @param[in] raw_data a Vector of char + * @return true if it succeeds, false if it cannot parse raw_data + */ + bool init_lookup_tables_from_raw_data (const amrex::Vector<char>& raw_data); + + /** + * Export lookup tables data into a raw binary Vector + * @return the data in binary format. The Vector is empty if tables were + * not previously initialized. + */ + amrex::Vector<char> export_lookup_tables_data () const; + + /** + * Computes the lookup tables. It does nothing unless WarpX is compiled with QED_TABLE_GEN=TRUE + * @param[in] ctrl control params to generate the tables + */ + void compute_lookup_tables (PicsarBreitWheelerCtrl ctrl); + + /** + * gets default (reasonable) values for the control parameters + * @return default control params to generate the tables + */ + PicsarBreitWheelerCtrl get_default_ctrl() const; + +private: + bool m_lookup_tables_initialized = false; + + BreitWheelerEngineInnards m_innards; + +//Table builing is available only if WarpX is compiled with QED_TABLE_GEN=TRUE +#ifdef WARPX_QED_TABLE_GEN + BreitWheelerEngineTableBuilder m_table_builder; +#endif + }; //============================================ diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp index 33e6d48d7..42953c97f 100644 --- a/Source/QED/BreitWheelerEngineWrapper.cpp +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -1,16 +1,188 @@ #include "BreitWheelerEngineWrapper.H" + +#include "QedTableParserHelperFunctions.H" + +#include <utility> + +using namespace std; +using namespace QedUtils; +using namespace amrex; + //This file provides a wrapper aroud the breit_wheeler engine //provided by the PICSAR library -using namespace picsar::multi_physics; - // Factory class ============================= BreitWheelerEngine::BreitWheelerEngine (){} -BreitWheelerGetOpticalDepth BreitWheelerEngine::build_optical_depth_functor () +BreitWheelerGetOpticalDepth +BreitWheelerEngine::build_optical_depth_functor () { return BreitWheelerGetOpticalDepth(); } +BreitWheelerEvolveOpticalDepth +BreitWheelerEngine::build_evolve_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return BreitWheelerEvolveOpticalDepth(m_innards); +} + +BreitWheelerGeneratePairs +BreitWheelerEngine::build_pair_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return BreitWheelerGeneratePairs(m_innards); +} + +bool BreitWheelerEngine::are_lookup_tables_initialized () const +{ + return m_lookup_tables_initialized; +} + +bool +BreitWheelerEngine::init_lookup_tables_from_raw_data ( + const Vector<char>& raw_data) +{ + const char* p_data = raw_data.data(); + const char* const p_last = &raw_data.back(); + bool is_ok; + + //Header (control parameters) + tie(is_ok, m_innards.ctrl.chi_phot_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tdndt_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tdndt_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tdndt_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tdndt_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tdndt_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tdndt_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tpair_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tpair_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tpair_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tpair_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tpair_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tpair_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_frac_tpair_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_frac_tpair_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + //___________________________ + + //Data + Vector<Real> tndt_coords(m_innards.ctrl.chi_phot_tdndt_how_many); + Vector<Real> tndt_data(m_innards.ctrl.chi_phot_tdndt_how_many); + Vector<Real> cum_tab_coords1(m_innards.ctrl.chi_phot_tpair_how_many); + Vector<Real> cum_tab_coords2(m_innards.ctrl.chi_frac_tpair_how_many); + Vector<Real> cum_tab_data(m_innards.ctrl.chi_phot_tpair_how_many* + m_innards.ctrl.chi_frac_tpair_how_many); + + tie(is_ok, tndt_coords, p_data) = + parse_raw_data_vec<Real>( + p_data, tndt_coords.size(), p_last); + if(!is_ok) return false; + m_innards.TTfunc_coords.assign(tndt_coords.begin(), tndt_coords.end()); + + tie(is_ok, tndt_data, p_data) = + parse_raw_data_vec<Real>( + p_data, tndt_data.size(), p_last); + if(!is_ok) return false; + m_innards.TTfunc_data.assign(tndt_data.begin(), tndt_data.end()); + + tie(is_ok, cum_tab_coords1, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_coords1.size(), p_last); + if(!is_ok) return false; + m_innards.cum_distrib_coords_1.assign( + cum_tab_coords1.begin(), cum_tab_coords1.end()); + + tie(is_ok, cum_tab_coords2, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_coords2.size(), p_last); + if(!is_ok) return false; + m_innards.cum_distrib_coords_2.assign( + cum_tab_coords2.begin(), cum_tab_coords2.end()); + + tie(is_ok, cum_tab_data, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_data.size(), p_last); + if(!is_ok) return false; + m_innards.cum_distrib_data.assign( + cum_tab_data.begin(), cum_tab_data.end()); + + //___________________________ + m_lookup_tables_initialized = true; + + return true; +} + +Vector<char> BreitWheelerEngine::export_lookup_tables_data () const +{ + Vector<char> res{}; + + if(!m_lookup_tables_initialized) + return res; + + add_data_to_vector_char(&m_innards.ctrl.chi_phot_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tdndt_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tdndt_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tdndt_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tpair_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tpair_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tpair_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_frac_tpair_how_many, 1, res); + + add_data_to_vector_char(m_innards.TTfunc_coords.data(), + m_innards.TTfunc_coords.size(), res); + add_data_to_vector_char(m_innards.TTfunc_data.data(), + m_innards.TTfunc_data.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_1.data(), + m_innards.cum_distrib_coords_1.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_2.data(), + m_innards.cum_distrib_coords_2.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_data.data(), + m_innards.cum_distrib_data.size(), res); + + return res; +} + +PicsarBreitWheelerCtrl +BreitWheelerEngine::get_default_ctrl() const +{ + return PicsarBreitWheelerCtrl(); +} + +void BreitWheelerEngine::compute_lookup_tables ( + PicsarBreitWheelerCtrl ctrl) +{ +#ifdef WARPX_QED_TABLE_GEN + m_table_builder.compute_table(ctrl, m_innards); + m_lookup_tables_initialized = true; +#endif +} + //============================================ diff --git a/Source/QED/Make.package b/Source/QED/Make.package index d4bad3f18..c9cd73cc2 100644 --- a/Source/QED/Make.package +++ b/Source/QED/Make.package @@ -1,8 +1,21 @@ CEXE_headers += QedWrapperCommons.H +CEXE_headers += QedChiFunctions.H +CEXE_headers += QedTableParserHelperFunctions.H +CEXE_headers += BreitWheelerEngineInnards.H +CEXE_headers += QuantumSyncEngineInnards.H CEXE_headers += BreitWheelerEngineWrapper.H -CEXE_headers += QuantumSyncsEngineWrapper.H +CEXE_headers += QuantumSyncEngineWrapper.H CEXE_sources += BreitWheelerEngineWrapper.cpp CEXE_sources += QuantumSyncEngineWrapper.cpp +#Table generation is enabled only if QED_TABLE_GEN is +#set to true +ifeq ($(QED_TABLE_GEN),TRUE) + CEXE_headers += BreitWheelerEngineTableBuilder.H + CEXE_headers += QuantumSyncEngineTableBuilder.H + CEXE_sources += BreitWheelerEngineTableBuilder.cpp + CEXE_sources += QuantumSyncEngineTableBuilder.cpp +endif + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/QED VPATH_LOCATIONS += $(WARPX_HOME)/Source/QED diff --git a/Source/QED/QedChiFunctions.H b/Source/QED/QedChiFunctions.H new file mode 100644 index 000000000..dd8ffac0e --- /dev/null +++ b/Source/QED/QedChiFunctions.H @@ -0,0 +1,62 @@ +#ifndef WARPX_amrex_qed_chi_functions_h_ +#define WARPX_amrex_qed_chi_functions_h_ + +/** + * This header contains wrappers around functions provided by + * the PICSAR QED library to calculate the 'chi' parameter + * for photons or electrons and positrons. + */ + +#include "QedWrapperCommons.H" + +//#define PXRMP_CORE_ONLY allows importing only the 'core functions' of the +//QED library. +#define PXRMP_CORE_ONLY +#include <chi_functions.hpp> + +namespace QedUtils{ + /** + * Function to calculate the 'chi' parameter for photons. + * Suitable for GPU kernels. + * @param[in] px,py,pz components of momentum (SI units) + * @param[in] ex,ey,ez components of electric field (SI units) + * @param[in] bx,by,bz components of magnetic field (SI units) + * @return chi parameter + */ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real chi_photon( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz) + { + //laser wavelength is unused if SI units are set + const amrex::Real dummy_lambda = 1.0; + return picsar::multi_physics::chi_photon( + px, py, pz, ex, ey, ez, bx, by, bz, dummy_lambda); + } + + /** + * Function to calculate the 'chi' parameter for electrons or positrons. + * Suitable for GPU kernels. + * @param[in] px,py,pz components of momentum (SI units) + * @param[in] ex,ey,ez components of electric field (SI units) + * @param[in] bx,by,bz components of magnetic field (SI units) + * @return chi parameter + */ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real chi_lepton( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz) + { + //laser wavelength is unused if SI units are set + const amrex::Real dummy_lambda = 1.0; + return picsar::multi_physics::chi_lepton( + px, py, pz, ex, ey, ez, bx, by, bz, dummy_lambda); + } + //_________ +}; + +#endif //WARPX_amrex_qed_chi_functions_h_ diff --git a/Source/QED/QedTableParserHelperFunctions.H b/Source/QED/QedTableParserHelperFunctions.H new file mode 100644 index 000000000..9f9f37017 --- /dev/null +++ b/Source/QED/QedTableParserHelperFunctions.H @@ -0,0 +1,90 @@ +#ifndef WARPX_amrex_qed_table_parser_helper_functions_h_ +#define WARPX_amrex_qed_table_parser_helper_functions_h_ + +/** + * This header contains helper functions to safely extract data + * (e.g. integers, floating point numbers) from raw binary data + * (i.e. a char*) and to convert arrays into raw binary data. + */ + +#include <AMReX_Vector.H> +#include <tuple> + +namespace QedUtils{ + /** + * This function safely extracts an amrex::Vector<T> from raw binary data. + * T must be a simple datatype (e.g. an int, a float, a double...). + * + * @param[in] p_char a pointer to the binary stream + * @param[in] how_many how many T should be read from stream + * @param[in] p_last a pointer to the last element of the char* array + * @return {a tuple containing + * 1) flag (which is false if p_last is exceeded) + * 2) a Vector of T + * 3) a pointer to a new location of the binary data (after having read how_many T)} + */ + template <class T> + std::tuple<bool, amrex::Vector<T>, const char*>parse_raw_data_vec( + const char* p_data, size_t how_many, const char* const p_last) + { + amrex::Vector<T> res; + if(p_data + sizeof(T)*how_many > p_last) + return std::make_tuple(false, res, nullptr); + + auto r_data = reinterpret_cast<const T*>(p_data); + + res.assign(r_data, r_data + how_many); + + p_data += sizeof(T)*how_many; + return std::make_tuple(true, res, p_data); + } + + /** + * This function safely extracts a T from raw binary data. + * T must be a simple datatype (e.g. an int, a float, a double...). + * + * @param[in] p_char a pointer to the binary stream + * @param[in] p_last a pointer to the last element of the char* array + * @return {a tuple containing + * 1) flag (which is false if p_last is exceeded) + * 2) a T + * 3) a pointer to a new location of the binary data (after having read 1 T)} + */ + template <class T> + std::tuple<bool, T, const char*> parse_raw_data( + const char* p_data, const char* const p_last) + { + T res; + if(p_data + sizeof(T) > p_last) + return std::make_tuple(false, res, nullptr); + + auto r_data = reinterpret_cast<const T*>(p_data); + + res = *r_data; + + p_data += sizeof(T); + return std::make_tuple(true, res, p_data); + } + + /** + * This function converts a C-style array of T into + * a Vector<char> (i.e. raw binary data) and adds it + * to an existing Vector<char> passed by reference + * @param[in] p_data a pointer to the beginning of the array + * @param[in] how_many number of elements of type T in the array + * @param[in,out] raw_data data will be appended to this vector + */ + template <class T> + void add_data_to_vector_char ( + const T* p_data, size_t how_many, amrex::Vector<char>& raw_data) + { + raw_data.insert( + raw_data.end(), + reinterpret_cast<const char*>(p_data), + reinterpret_cast<const char*>(p_data) + + sizeof(T)*how_many + ); + } +}; + +#endif //WARPX_amrex_qed_table_parser_helper_functions_h_ diff --git a/Source/QED/QedWrapperCommons.H b/Source/QED/QedWrapperCommons.H index 821034c06..210e7247e 100644 --- a/Source/QED/QedWrapperCommons.H +++ b/Source/QED/QedWrapperCommons.H @@ -1,16 +1,36 @@ #ifndef WARPX_amrex_qed_wrapper_commons_h_ #define WARPX_amrex_qed_wrapper_commons_h_ -//Common definitions for the QED library wrappers +/** + * This header contains some common #define directives and a + * 'dummy' class used by the QED library wrappers and related + * components. + */ #include <AMReX_AmrCore.H> +#include <AMReX_GpuQualifiers.H> -//Sets the decorator for GPU -#define PXRMP_GPU AMREX_GPU_DEVICE -//Sets SI units in the library +/** + * PICSAR uses PXRMP_GPU to decorate methods which should be + * compiled for GPU. The user has to set it to the right value + * (AMREX_GPU_DEVICE in this case). + * PXRMP_WITH_SI_UNITS sets the library to use International + * System units. + */ +#define PXRMP_GPU AMREX_GPU_HOST_DEVICE #define PXRMP_WITH_SI_UNITS +//_________________________ + +/** + * A namespace called 'QedUtils' is used to encapsulate + * free functions (defined elsewhere) and an + * empty datastructure (DummyStruct), which is re-used by several + * components. + */ +namespace QedUtils{ + struct DummyStruct{}; +}; +//_________________________ -//An empty data type -struct DummyStruct{}; #endif //WARPX_amrex_qed_wrapper_commons_h_ diff --git a/Source/QED/QuantumSyncEngineInnards.H b/Source/QED/QuantumSyncEngineInnards.H new file mode 100644 index 000000000..6206b103a --- /dev/null +++ b/Source/QED/QuantumSyncEngineInnards.H @@ -0,0 +1,46 @@ +#ifndef WARPX_quantum_sync_engine_innards_h_ +#define WARPX_quantum_sync_engine_innards_h_ + +#include "QedWrapperCommons.H" + +#include <AMReX_Gpu.H> + +//This includes only the definition of a simple datastructure +//used to control the Quantum Synchrotron engine. +#include <quantum_sync_engine_ctrl.h> + +/** + * This structure holds all the parameters required to use the + * Quantum Synchrotron engine: a POD control structure and lookup + * tables data. + */ +struct QuantumSynchrotronEngineInnards +{ + // Control parameters (a POD struct) + // ctrl contains several parameters: + // - chi_part_min : the minium chi parameter to be + // considered by the engine + // - chi_part_tdndt_min : minimun chi for sub-table 1 (1D) + // - chi_part_tdndt_max : maximum chi for sub-table 1 (1D) + // - chi_part_tdndt_how_many : how many points to use for sub-table 1 (1D) + // - chi_part_tem_min : minimun chi for sub-table 2 (1D) + // - chi_part_tem_max : maximum chi for sub-table 2 (1D) + // - chi_part_tem_how_many : how many points to use for chi for sub-table 2 (2D) + // - prob_tem_how_many : how many points to use for the second axis of sub-table 2 (2D) + picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real> ctrl; + + //Lookup table data + //---sub-table 1 (1D) + amrex::Gpu::ManagedDeviceVector<amrex::Real> KKfunc_coords; + amrex::Gpu::ManagedDeviceVector<amrex::Real> KKfunc_data; + //--- + + //---sub-table 2 (2D) + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_1; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_2; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_data; + //______ +}; +//========================================================== + +#endif //WARPX_quantum_sync_engine_innards_h_ diff --git a/Source/QED/QuantumSyncEngineTableBuilder.H b/Source/QED/QuantumSyncEngineTableBuilder.H new file mode 100644 index 000000000..e70f5d02f --- /dev/null +++ b/Source/QED/QuantumSyncEngineTableBuilder.H @@ -0,0 +1,26 @@ +#ifndef WARPX_quantum_sync_engine_table_builder_h_ +#define WARPX_quantum_sync_engine_table_builder_h_ + +#include "QedWrapperCommons.H" +#include "QuantumSyncEngineInnards.H" + +//This includes only the definition of a simple datastructure +//used to control the Quantum Synchrotron engine. +#include <quantum_sync_engine_ctrl.h> + +/** + * A class which computes the lookup tables for the Quantum Synchrotron engine. + */ +class QuantumSynchrotronEngineTableBuilder{ +public: + /** + * Computes the tables. + * @param[in] ctrl control parameters to generate the tables + * @param[out] innards structure holding both a copy of ctrl and lookup tables data + */ + void compute_table + (picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real> ctrl, + QuantumSynchrotronEngineInnards& innards) const; +}; + +#endif //WARPX_quantum_sync_engine_table_builder_h_ diff --git a/Source/QED/QuantumSyncEngineTableBuilder.cpp b/Source/QED/QuantumSyncEngineTableBuilder.cpp new file mode 100644 index 000000000..51c3720f2 --- /dev/null +++ b/Source/QED/QuantumSyncEngineTableBuilder.cpp @@ -0,0 +1,58 @@ +#include "QuantumSyncEngineTableBuilder.H" + +//Include the full Quantum Synchrotron engine with table generation support +//(after some consistency tests). This requires to have a recent version +// of the Boost library. +#ifdef PXRMP_CORE_ONLY + #error The Table Builder is incompatible with PXRMP_CORE_ONLY +#endif + +#ifdef __PICSAR_MULTIPHYSICS_BREIT_WHEELER_ENGINE__ + #warning quantum_sync_engine.hpp should not have been included before reaching this point. +#endif +#include <quantum_sync_engine.hpp> +//_______________________________________________ + +//Some handy aliases +using PicsarQuantumSynchrotronEngine = picsar::multi_physics:: + quantum_synchrotron_engine<amrex::Real, QedUtils::DummyStruct>; + +using PicsarQuantumSynchrotronCtrl = + picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real>; +//_______________________________________________ + +void +QuantumSynchrotronEngineTableBuilder::compute_table + (PicsarQuantumSynchrotronCtrl ctrl, + QuantumSynchrotronEngineInnards& innards) const +{ + PicsarQuantumSynchrotronEngine qs_engine( + std::move(QedUtils::DummyStruct()), 1.0, ctrl); + + qs_engine.compute_dN_dt_lookup_table(); + qs_engine.compute_cumulative_phot_em_table(); + + auto qs_innards_picsar = qs_engine.export_innards(); + + //Copy data in a GPU-friendly data-structure + innards.ctrl = qs_innards_picsar.qs_ctrl; + innards.KKfunc_coords.assign(qs_innards_picsar.KKfunc_table_coords_ptr, + qs_innards_picsar.KKfunc_table_coords_ptr + + qs_innards_picsar.KKfunc_table_coords_how_many); + innards.KKfunc_data.assign(qs_innards_picsar.KKfunc_table_data_ptr, + qs_innards_picsar.KKfunc_table_data_ptr + + qs_innards_picsar.KKfunc_table_data_how_many); + innards.cum_distrib_coords_1.assign( + qs_innards_picsar.cum_distrib_table_coords_1_ptr, + qs_innards_picsar.cum_distrib_table_coords_1_ptr + + qs_innards_picsar.cum_distrib_table_coords_1_how_many); + innards.cum_distrib_coords_2.assign( + qs_innards_picsar.cum_distrib_table_coords_2_ptr, + qs_innards_picsar.cum_distrib_table_coords_2_ptr + + qs_innards_picsar.cum_distrib_table_coords_2_how_many); + innards.cum_distrib_data.assign( + qs_innards_picsar.cum_distrib_table_data_ptr, + qs_innards_picsar.cum_distrib_table_data_ptr + + qs_innards_picsar.cum_distrib_table_data_how_many); + //____ +} diff --git a/Source/QED/QuantumSyncEngineWrapper.H b/Source/QED/QuantumSyncEngineWrapper.H index a285dc8b6..1a6ffe4f3 100644 --- a/Source/QED/QuantumSyncEngineWrapper.H +++ b/Source/QED/QuantumSyncEngineWrapper.H @@ -2,18 +2,45 @@ #define WARPX_quantum_sync_engine_wrapper_h_ #include "QedWrapperCommons.H" +#include "QuantumSyncEngineInnards.H" -//QS ENGINE from PICSAR +#include <AMReX_Array.H> +#include <AMReX_Vector.H> +#include <AMReX_Gpu.H> + +//#define PXRMP_CORE_ONLY allows importing only the 'core functions' of the +//Quantum Synchrotron engine of the QED PICSAR library. #define PXRMP_CORE_ONLY -#include "quantum_sync_engine.hpp" +#include <quantum_sync_engine.hpp> + +//Lookup table building function is in a dedicated (optional) class to +//avoid including heavy dependencies if they are not needed. +#ifdef WARPX_QED_TABLE_GEN + #include "QuantumSyncEngineTableBuilder.H" +#endif + +#include <string> + +//Some handy aliases + +// The engine has two templated arguments: the numerical type +// and a random number generator. However, random numbers are not +// used to generate the lookup tables and the static member +// functions which are called in the functors do not use +// random numbers as well. Therefore, an empty "DummyStruct" +// can be passed as a template parameter. +using PicsarQuantumSynchrotronEngine = picsar::multi_physics:: + quantum_synchrotron_engine<amrex::Real, QedUtils::DummyStruct>; -using WarpXQuantumSynchrotronWrapper = - picsar::multi_physics::quantum_synchrotron_engine<amrex::Real, DummyStruct>; +using PicsarQuantumSynchrotronCtrl = + picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real>; +//__________ // Functors ================================== -// These functors provide the core elementary functions of the library -// Can be included in GPU kernels +// These functors allow using the core elementary functions of the library. +// They are generated by a factory class (QuantumSynchrotronEngine, see below). +// They can be included in GPU kernels. /** * Functor to initialize the optical depth of leptons for the @@ -22,33 +49,252 @@ using WarpXQuantumSynchrotronWrapper = class QuantumSynchrotronGetOpticalDepth { public: + /** + * Constructor does nothing because optical depth initialization + * does not require control parameters or lookup tables. + */ QuantumSynchrotronGetOpticalDepth () {}; + /** + * () operator is just a thin wrapper around a very simple function to + * generate the optical depth. It can be used on GPU. + */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE - amrex::Real operator() () const + amrex::Real operator() () const noexcept { - return WarpXQuantumSynchrotronWrapper:: + //A random number in [0,1) should be provided as an argument. + return PicsarQuantumSynchrotronEngine:: internal_get_optical_depth(amrex::Random()); } }; //____________________________________________ +/** + * Functor to evolve the optical depth of leptons due to the + * Quantum Synchrotron process + */ +class QuantumSynchrotronEvolveOpticalDepth +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + QuantumSynchrotronEvolveOpticalDepth( + QuantumSynchrotronEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_KKfunc_size{r_innards.KKfunc_coords.size()}, + m_p_KKfunc_coords{r_innards.KKfunc_coords.dataPtr()}, + m_p_KKfunc_data{r_innards.KKfunc_data.dataPtr()} + {}; + + /** + * Evolves the optical depth. It can be used on GPU. + * @param[in] px,py,pz momentum components of the lepton (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] dt timestep (SI units) + * @param[in,out] opt_depth optical depth of the lepton. It is modified by the method. + * @return a flag which is 1 if optical depth becomes negative (i.e. a photon has to be generated). + */ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + int operator()( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real dt, amrex::Real& opt_depth) const noexcept + { + bool has_event_happened{false}; + + //the library provides the time (< dt) at which the event occurs, but this + //feature won't be used in WarpX for now. + amrex::Real unused_event_time{0.0}; + + PicsarQuantumSynchrotronEngine:: + internal_evolve_opt_depth_and_determine_event( + px, py, pz, + ex, ey, ez, + bx, by, bz, + dt, opt_depth, + has_event_happened, unused_event_time, + m_dummy_lambda, + picsar::multi_physics::lookup_1d<amrex::Real>{ + m_KKfunc_size, + m_p_KKfunc_coords, + m_p_KKfunc_data}, + m_ctrl); + + return has_event_happened; + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarQuantumSynchrotronCtrl m_ctrl; + + //lookup table data + size_t m_KKfunc_size; + amrex::Real* m_p_KKfunc_coords; + amrex::Real* m_p_KKfunc_data; +}; + +/** + * Functor to generate a photon via the Quantum Synchrotron process + * and to update momentum accordingly + */ +class QuantumSynchrotronGeneratePhotonAndUpdateMomentum +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + QuantumSynchrotronGeneratePhotonAndUpdateMomentum( + QuantumSynchrotronEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_cum_distrib_coords_1_size{r_innards.cum_distrib_coords_1.size()}, + m_cum_distrib_coords_2_size{r_innards.cum_distrib_coords_2.size()}, + m_p_distrib_coords_1{r_innards.cum_distrib_coords_1.data()}, + m_p_distrib_coords_2{r_innards.cum_distrib_coords_2.data()}, + m_p_cum_distrib_data{r_innards.cum_distrib_data.data()} + {}; + + /** + * Generates sampling (template parameter) photons according to Quantum Synchrotron process. + * It can be used on GPU. + * @param[in,out] px,py,pz momentum components of the lepton. They are modified (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] weight of the lepton (code units) + * @param[out] g_px,g_py,g_pz momenta of generated photons. Each array should have size=sampling (SI units) + * @param[out] g_weight weight of the generated photons. Array should have size=sampling (code units) + */ + template <size_t sampling> + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void operator()( + amrex::Real* px, amrex::Real* py, amrex::Real* pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real weight, + amrex::Real* g_px, amrex::Real* g_py, amrex::Real* g_pz, + amrex::Real* g_weight) const noexcept + { + //[sampling] random numbers are needed + amrex::GpuArray<amrex::Real, sampling> + rand_zero_one_minus_epsi; + for(auto& el : rand_zero_one_minus_epsi) el = amrex::Random(); + + PicsarQuantumSynchrotronEngine:: + internal_generate_photons_and_update_momentum( + *px, *py, *pz, + ex, ey, ez, + bx, by, bz, + weight, sampling, + g_px, g_py, g_pz, + g_weight, + m_dummy_lambda, + picsar::multi_physics::lookup_2d<amrex::Real>{ + m_cum_distrib_coords_1_size, + m_p_distrib_coords_1, + m_cum_distrib_coords_2_size, + m_p_distrib_coords_2, + m_p_cum_distrib_data}, + m_ctrl, + rand_zero_one_minus_epsi.data()); + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarQuantumSynchrotronCtrl m_ctrl; + + //lookup table data + size_t m_cum_distrib_coords_1_size; + size_t m_cum_distrib_coords_2_size; + amrex::Real* m_p_distrib_coords_1; + amrex::Real* m_p_distrib_coords_2; + amrex::Real* m_p_cum_distrib_data; +}; + // Factory class ============================= /** - * \brief Wrapper for the Quantum Synchrotron engine of the PICSAR library + * Wrapper for the Quantum Synchrotron engine of the PICSAR library */ class QuantumSynchrotronEngine { public: + /** + * Constructor requires no arguments. + */ QuantumSynchrotronEngine (); /** - * \brief Builds the functor to initialize the optical depth + * Builds the functor to initialize the optical depth */ QuantumSynchrotronGetOpticalDepth build_optical_depth_functor (); + + /** + * Builds the functor to evolve the optical depth + */ + QuantumSynchrotronEvolveOpticalDepth build_evolve_functor (); + + /** + * Builds the functor to generate photons + */ + QuantumSynchrotronGeneratePhotonAndUpdateMomentum build_phot_em_functor (); + + /** + * Checks if the optical tables are properly initialized + */ + bool are_lookup_tables_initialized () const; + + /** + * Init lookup tables from raw binary data. + * @param[in] raw_data a Vector of char + * @return true if it succeeds, false if it cannot parse raw_data + */ + bool init_lookup_tables_from_raw_data (const amrex::Vector<char>& raw_data); + + /** + * Export lookup tables data into a raw binary Vector + * @return the data in binary format. The Vector is empty if tables were + * not previously initialized. + */ + amrex::Vector<char> export_lookup_tables_data () const; + + /** + * Computes the lookup tables. It does nothing unless WarpX is compiled with QED_TABLE_GEN=TRUE + * @param[in] ctrl control params to generate the tables + */ + void compute_lookup_tables (PicsarQuantumSynchrotronCtrl ctrl); + + /** + * gets default (reasonable) values for the control parameters + * @return default control params to generate the tables + */ + PicsarQuantumSynchrotronCtrl get_default_ctrl() const; + +private: + bool m_lookup_tables_initialized = false; + + QuantumSynchrotronEngineInnards m_innards; + +//Table builing is available only if the libray is compiled with QED_TABLE_GEN=TRUE +#ifdef WARPX_QED_TABLE_GEN + QuantumSynchrotronEngineTableBuilder m_table_builder; +#endif + }; //============================================ diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp index b55149187..b2630dc4d 100644 --- a/Source/QED/QuantumSyncEngineWrapper.cpp +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -1,9 +1,16 @@ #include "QuantumSyncEngineWrapper.H" + +#include "QedTableParserHelperFunctions.H" + +#include <utility> + +using namespace std; +using namespace QedUtils; +using namespace amrex; + //This file provides a wrapper aroud the quantum_sync engine //provided by the PICSAR library -using namespace picsar::multi_physics; - // Factory class ============================= QuantumSynchrotronEngine::QuantumSynchrotronEngine (){} @@ -14,4 +21,167 @@ QuantumSynchrotronEngine::build_optical_depth_functor () return QuantumSynchrotronGetOpticalDepth(); } +QuantumSynchrotronEvolveOpticalDepth QuantumSynchrotronEngine::build_evolve_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return QuantumSynchrotronEvolveOpticalDepth(m_innards); +} + +QuantumSynchrotronGeneratePhotonAndUpdateMomentum QuantumSynchrotronEngine::build_phot_em_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return QuantumSynchrotronGeneratePhotonAndUpdateMomentum(m_innards); + +} + +bool QuantumSynchrotronEngine::are_lookup_tables_initialized () const +{ + return m_lookup_tables_initialized; +} + +bool +QuantumSynchrotronEngine::init_lookup_tables_from_raw_data ( + const Vector<char>& raw_data) +{ + const char* p_data = raw_data.data(); + const char* const p_last = &raw_data.back(); + bool is_ok; + + //Header (control parameters) + tie(is_ok, m_innards.ctrl.chi_part_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tdndt_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tdndt_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tdndt_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tdndt_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tdndt_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tdndt_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tem_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tem_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tem_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tem_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tem_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tem_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.prob_tem_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.prob_tem_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + //___________________________ + + //Data + Vector<Real> tndt_coords(m_innards.ctrl.chi_part_tdndt_how_many); + Vector<Real> tndt_data(m_innards.ctrl.chi_part_tdndt_how_many); + Vector<Real> cum_tab_coords1(m_innards.ctrl.chi_part_tem_how_many); + Vector<Real> cum_tab_coords2(m_innards.ctrl.prob_tem_how_many); + Vector<Real> cum_tab_data(m_innards.ctrl.chi_part_tem_how_many* + m_innards.ctrl.prob_tem_how_many); + + tie(is_ok, tndt_coords, p_data) = + parse_raw_data_vec<Real>( + p_data, tndt_coords.size(), p_last); + if(!is_ok) return false; + m_innards.KKfunc_coords.assign(tndt_coords.begin(), tndt_coords.end()); + + tie(is_ok, tndt_data, p_data) = + parse_raw_data_vec<Real>( + p_data, tndt_data.size(), p_last); + if(!is_ok) return false; + m_innards.KKfunc_data.assign(tndt_data.begin(), tndt_data.end()); + + tie(is_ok, cum_tab_coords1, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_coords1.size(), p_last); + if(!is_ok) return false; + m_innards.cum_distrib_coords_1.assign( + cum_tab_coords1.begin(), cum_tab_coords1.end()); + + tie(is_ok, cum_tab_coords2, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_coords2.size(), p_last); + if(!is_ok) return false; + m_innards.cum_distrib_coords_2.assign( + cum_tab_coords2.begin(), cum_tab_coords2.end()); + + tie(is_ok, cum_tab_data, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_data.size(), p_last); + if(!is_ok) return false; + m_innards.cum_distrib_data.assign( + cum_tab_data.begin(), cum_tab_data.end()); + + //___________________________ + m_lookup_tables_initialized = true; + + return true; +} + +Vector<char> QuantumSynchrotronEngine::export_lookup_tables_data () const +{ + Vector<char> res{}; + + if(!m_lookup_tables_initialized) + return res; + + add_data_to_vector_char(&m_innards.ctrl.chi_part_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tdndt_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tdndt_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tdndt_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tem_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tem_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tem_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.prob_tem_how_many, 1, res); + + add_data_to_vector_char(m_innards.KKfunc_coords.data(), + m_innards.KKfunc_coords.size(), res); + add_data_to_vector_char(m_innards.KKfunc_data.data(), + m_innards.KKfunc_data.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_1.data(), + m_innards.cum_distrib_coords_1.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_2.data(), + m_innards.cum_distrib_coords_2.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_data.data(), + m_innards.cum_distrib_data.size(), res); + + return res; +} + +PicsarQuantumSynchrotronCtrl +QuantumSynchrotronEngine::get_default_ctrl() const +{ + return PicsarQuantumSynchrotronCtrl(); +} + +void QuantumSynchrotronEngine::compute_lookup_tables ( + PicsarQuantumSynchrotronCtrl ctrl) +{ +#ifdef WARPX_QED_TABLE_GEN + m_table_builder.compute_table(ctrl, m_innards); + m_lookup_tables_initialized = true; +#endif +} + //============================================ diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index f84701a02..c577da7f3 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -99,8 +99,8 @@ WarpX::MoveWindow (bool move_j) for (int dim = 0; dim < 3; ++dim) { // Fine grid - shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir); - shiftMF(*Efield_fp[lev][dim], geom[lev], num_shift, dir); + shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir, B_external_grid[dim]); + shiftMF(*Efield_fp[lev][dim], geom[lev], num_shift, dir, E_external_grid[dim]); if (move_j) { shiftMF(*current_fp[lev][dim], geom[lev], num_shift, dir); } @@ -113,8 +113,8 @@ WarpX::MoveWindow (bool move_j) if (lev > 0) { // Coarse grid - shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir); - shiftMF(*Efield_cp[lev][dim], geom[lev-1], num_shift_crse, dir); + shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, B_external_grid[dim]); + shiftMF(*Efield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, E_external_grid[dim]); shiftMF(*Bfield_aux[lev][dim], geom[lev], num_shift, dir); shiftMF(*Efield_aux[lev][dim], geom[lev], num_shift, dir); if (move_j) { @@ -203,7 +203,8 @@ WarpX::MoveWindow (bool move_j) } void -WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir) +WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir, + amrex::Real external_field) { BL_PROFILE("WarpX::shiftMF()"); const BoxArray& ba = mf.boxArray(); @@ -257,7 +258,7 @@ WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir) if (outbox.ok()) { AMREX_PARALLEL_FOR_4D ( outbox, nc, i, j, k, n, { - srcfab(i,j,k,n) = 0.0; + srcfab(i,j,k,n) = external_field; }); } diff --git a/Source/Utils/WarpXTagging.cpp b/Source/Utils/WarpXTagging.cpp index 8ea3211a3..91bb802e8 100644 --- a/Source/Utils/WarpXTagging.cpp +++ b/Source/Utils/WarpXTagging.cpp @@ -22,9 +22,9 @@ WarpX::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/) for (BoxIterator bi(bx); bi.ok(); ++bi) { const IntVect& cell = bi(); - RealVect pos {AMREX_D_DECL((cell[0]+0.5)*dx[0]+problo[0], - (cell[1]+0.5)*dx[1]+problo[1], - (cell[2]+0.5)*dx[2]+problo[2])}; + RealVect pos {AMREX_D_DECL((cell[0]+0.5_rt)*dx[0]+problo[0], + (cell[1]+0.5_rt)*dx[1]+problo[1], + (cell[2]+0.5_rt)*dx[2]+problo[2])}; if (pos > fine_tag_lo && pos < fine_tag_hi) { fab(cell) = TagBox::SET; } diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index fd6e72dc6..ab28c5446 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -2,6 +2,8 @@ #include <AMReX_Vector.H> #include <AMReX_MultiFab.H> +#include <string> + void ReadBoostedFrameParameters(amrex::Real& gamma_boost, amrex::Real& beta_boost, amrex::Vector<int>& boost_direction); @@ -9,3 +11,13 @@ void ConvertLabParamsToBoost(); void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax); + +namespace WarpXUtilIO{ + /** + * A helper function to write binary data on disk. + * @param[in] filename where to write + * @param[in] data Vector containing binary data to write on disk + * return true if it succeeds, false otherwise + */ + bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data); +} diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 4b11eb69d..8764a09c6 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -1,10 +1,11 @@ -#include <cmath> - #include <WarpXUtil.H> #include <WarpXConst.H> #include <AMReX_ParmParse.H> #include <WarpX.H> +#include <cmath> +#include <fstream> + using namespace amrex; void ReadBoostedFrameParameters(Real& gamma_boost, Real& beta_boost, @@ -152,3 +153,14 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) } } } + +namespace WarpXUtilIO{ + bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data) + { + std::ofstream of{filename, std::ios::binary}; + of.write(data.data(), data.size()); + of.close(); + return of.good(); + } +} + diff --git a/Source/WarpX.H b/Source/WarpX.H index ead9c6102..04231b3be 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -4,7 +4,7 @@ #include "WarpXDtType.H" #include <MultiParticleContainer.H> #include <PML.H> -#include <BoostedFrameDiagnostic.H> +#include <BackTransformedDiagnostic.H> #include <BilinearFilter.H> #include <NCIGodfreyFilter.H> @@ -69,13 +69,18 @@ public: MultiParticleContainer& GetPartContainer () { return *mypc; } - static void shiftMF(amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, int dir); + static void shiftMF(amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, int dir, + amrex::Real external_field = 0.); static void GotoNextLine (std::istream& is); - // External fields - static amrex::Vector<amrex::Real> B_external; - static amrex::Vector<amrex::Real> E_external; + // External fields added to particle fields. + static amrex::Vector<amrex::Real> B_external_particle; + static amrex::Vector<amrex::Real> E_external_particle; + + // Initial field on the grid. + static amrex::Vector<amrex::Real> E_external_grid; + static amrex::Vector<amrex::Real> B_external_grid; // Algorithms static long current_deposition_algo; @@ -100,12 +105,12 @@ public: static bool serialize_ics; // Back transformation diagnostic - static bool do_boosted_frame_diagnostic; + static bool do_back_transformed_diagnostics; static std::string lab_data_directory; static int num_snapshots_lab; static amrex::Real dt_snapshots_lab; - static bool do_boosted_frame_fields; - static bool do_boosted_frame_particles; + static bool do_back_transformed_fields; + static bool do_back_transformed_particles; // Boosted frame parameters static amrex::Real gamma_boost; @@ -479,7 +484,7 @@ private: std::unique_ptr<MultiParticleContainer> mypc; // Boosted Frame Diagnostics - std::unique_ptr<BoostedFrameDiagnostic> myBFD; + std::unique_ptr<BackTransformedDiagnostic> myBFD; // // Fields: First array for level, second for direction diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 00d1cb4ac..25b33ea82 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -24,8 +24,11 @@ using namespace amrex; -Vector<Real> WarpX::B_external(3, 0.0); -Vector<Real> WarpX::E_external(3, 0.0); +Vector<Real> WarpX::B_external_particle(3, 0.0); +Vector<Real> WarpX::E_external_particle(3, 0.0); + +Vector<Real> WarpX::E_external_grid(3, 0.0); +Vector<Real> WarpX::B_external_grid(3, 0.0); int WarpX::do_moving_window = 0; int WarpX::moving_window_dir = -1; @@ -61,12 +64,12 @@ int WarpX::num_mirrors = 0; int WarpX::sort_int = -1; -bool WarpX::do_boosted_frame_diagnostic = false; +bool WarpX::do_back_transformed_diagnostics = false; std::string WarpX::lab_data_directory = "lab_frame_data"; int WarpX::num_snapshots_lab = std::numeric_limits<int>::lowest(); Real WarpX::dt_snapshots_lab = std::numeric_limits<Real>::lowest(); -bool WarpX::do_boosted_frame_fields = true; -bool WarpX::do_boosted_frame_particles = true; +bool WarpX::do_back_transformed_fields = true; +bool WarpX::do_back_transformed_particles = true; int WarpX::num_slice_snapshots_lab = 0; Real WarpX::dt_slice_snapshots_lab; @@ -169,7 +172,7 @@ WarpX::WarpX () current_injection_position = geom[0].ProbLo(moving_window_dir); } } - do_boosted_frame_particles = mypc->doBoostedFrameDiags(); + do_back_transformed_particles = mypc->doBackTransformedDiagnostics(); Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); @@ -295,8 +298,11 @@ WarpX::ReadParameters () pp.query("zmax_plasma_to_compute_max_step", zmax_plasma_to_compute_max_step); - pp.queryarr("B_external", B_external); - pp.queryarr("E_external", E_external); + pp.queryarr("B_external_particle", B_external_particle); + pp.queryarr("E_external_particle", E_external_particle); + + pp.queryarr("E_external_grid", E_external_grid); + pp.queryarr("B_external_grid", B_external_grid); pp.query("do_moving_window", do_moving_window); if (do_moving_window) @@ -325,8 +331,8 @@ WarpX::ReadParameters () moving_window_v *= PhysConst::c; } - pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); - if (do_boosted_frame_diagnostic) { + pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics); + if (do_back_transformed_diagnostics) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, "gamma_boost must be > 1 to use the boosted frame diagnostic."); @@ -354,7 +360,7 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); - pp.query("do_boosted_frame_fields", do_boosted_frame_fields); + pp.query("do_back_transformed_fields", do_back_transformed_fields); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, "The moving window should be on if using the boosted frame diagnostic."); @@ -601,7 +607,7 @@ WarpX::ReadParameters () } } - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, "gamma_boost must be > 1 to use the boost frame diagnostic"); pp.query("num_slice_snapshots_lab", num_slice_snapshots_lab); |