diff options
Diffstat (limited to 'Source/Particles/MultiParticleContainer.cpp')
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 484 |
1 files changed, 331 insertions, 153 deletions
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 715c97b99..55768c6fc 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,11 +1,17 @@ -#include <limits> -#include <algorithm> -#include <string> - #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> + using namespace amrex; MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) @@ -37,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 @@ -149,6 +156,11 @@ MultiParticleContainer::InitData () // For each species, get the ID of its product species. // This is used for ionization and pair creation processes. mapSpeciesProduct(); + +#ifdef WARPX_QED + InitQED(); +#endif + } @@ -381,8 +393,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); @@ -520,146 +532,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 () { @@ -721,8 +593,314 @@ 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 } + +#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 + (m_shr_p_qs_engine); + m_nspecies_quantum_sync++; + } + if(pc->has_breit_wheeler()){ + pc->set_breit_wheeler_engine_ptr + (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 |