From 7b9aa6dd3108f10ce088614abccb239645a5dfe0 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sat, 3 Aug 2019 17:27:40 -0700 Subject: add ionization_level component --- Source/Particles/PhysicalParticleContainer.cpp | 56 ++++++++++++++++++++------ 1 file changed, 43 insertions(+), 13 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7d63bd8e7..d0db0fbdc 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -37,6 +37,12 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); + pp.query("do_field_ionization", do_field_ionization); + if (do_field_ionization){ + AddRealComp("ionization_level"); + plot_flags.resize(PIdx::nattribs + 1, 1); + } + pp.query("plot_species", plot_species); int do_user_plot_vars; do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); @@ -73,6 +79,10 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) : WarpXParticleContainer(amr_core, 0) { plasma_injector.reset(new PlasmaInjector()); + if (do_field_ionization){ + Print()<<"AddRealComp\n"; + AddRealComp("ionization_level"); + } } void PhysicalParticleContainer::InitData() @@ -196,17 +206,20 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { // need to create old values auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } } // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); @@ -289,7 +302,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); GetParticles(lev)[std::make_pair(grid_id, tile_id)]; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { DefineAndReturnParticleTile(lev, grid_id, tile_id); } } @@ -416,10 +429,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; bool do_boosted = false; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - do_boosted = true; + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { DefineAndReturnParticleTile(lev, grid_id, tile_id); } + do_boosted = WarpX::do_boosted_frame_diagnostic + && do_boosted_frame_diags; + auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); @@ -439,7 +454,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; } - + Real* pi; + if (do_field_ionization) { + pi = soa.GetRealData(particle_comps[ "ionization_level"]).data() + old_size; + } + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -568,6 +587,10 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) u.z = gamma_boost * ( u.z -beta_boost*gamma_lab ); } + if (do_field_ionization) { + pi[ip] = 2.; // species_ionization_level; + } + u.x *= PhysConst::c; u.y *= PhysConst::c; u.z *= PhysConst::c; @@ -1351,13 +1374,20 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt); } } else { + Real* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetAttribs(particle_comps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + // Deposit inside domains - DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, lev, lev, dt); if (has_buffer){ // Deposit in buffers - DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, np_current, np-np_current, thread_num, lev, lev-1, dt); } -- cgit v1.2.3 From 58d728b3160c69ce8e9f117b255e7b057cc419fe Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sat, 3 Aug 2019 17:45:23 -0700 Subject: add ionization initialization --- Source/Evolve/WarpXEvolveEM.cpp | 2 +- Source/Particles/MultiParticleContainer.H | 5 ++ Source/Particles/MultiParticleContainer.cpp | 51 ++++++++++++++++++++ Source/Particles/PhysicalParticleContainer.H | 22 +++++++++ Source/Particles/PhysicalParticleContainer.cpp | 66 ++++++++++++++++++++++---- Source/Particles/WarpXParticleContainer.H | 14 +++++- 6 files changed, 149 insertions(+), 11 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index aa180cc95..4bd9d4566 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -278,7 +278,7 @@ WarpX::OneStep_nosub (Real cur_time) { // Loop over species. For each ionizable species, create particles in // product species. - // mypc->doFieldIonization(); + mypc->doFieldIonization(); // Push particle from x^{n} to x^{n+1} // from p^{n-1/2} to p^{n+1/2} // Deposit current j^{n+1/2} diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 7c9ede411..4ea9c2763 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -128,6 +128,8 @@ public: /// std::unique_ptr GetChargeDensity(int lev, bool local = false); + void doFieldIonization(); + void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; @@ -207,6 +209,9 @@ private: void ReadParameters (); + 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 diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 982e04e39..16c9b9762 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,10 +1,12 @@ #include #include #include +#include #include #include #include +#include using namespace amrex; @@ -29,6 +31,10 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); + // For each species, get the ID of its product species. + // This is used for ionization and pair creation processes. + mapSpeciesProduct(); + // 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. @@ -461,3 +467,48 @@ MultiParticleContainer::doContinuousInjection() const } return warpx_do_continuous_injection; } + +/* \brief Get ID of product species of each species. + * The users specifies the name of the product species, + * this routine get its ID. + */ +void +MultiParticleContainer::mapSpeciesProduct() +{ + for (int i=0; iionization_product_name and store its ID into + // pc->ionization_product. + if (pc->do_field_ionization){ + int i_product = getSpeciesID(pc->ionization_product_name); + pc->ionization_product = i_product; + } + } +} + +/* \brief Given a species name, return its ID. + */ +int +MultiParticleContainer::getSpeciesID(std::string product_str) +{ + int i_product; + bool found = 0; + // Loop over species + for (int i=0; i, 3> >& E, const amrex::Vector > > >& masks) override; @@ -105,6 +109,24 @@ public: void SplitParticles(int lev); + virtual void copyParticles(int lev) override; + + virtual int doFieldIonization( + const int lev, + amrex::Cuda::ManagedDeviceVector& x_buf, + amrex::Cuda::ManagedDeviceVector& y_buf, + amrex::Cuda::ManagedDeviceVector& z_buf, + amrex::Cuda::ManagedDeviceVector& ux_buf, + amrex::Cuda::ManagedDeviceVector& uy_buf, + amrex::Cuda::ManagedDeviceVector& uz_buf, + amrex::Cuda::ManagedDeviceVector& ex_buf, + amrex::Cuda::ManagedDeviceVector& ey_buf, + amrex::Cuda::ManagedDeviceVector& ez_buf, + amrex::Cuda::ManagedDeviceVector& bx_buf, + amrex::Cuda::ManagedDeviceVector& by_buf, + amrex::Cuda::ManagedDeviceVector& bz_buf, + amrex::Cuda::ManagedDeviceVector& w_buf) override; + // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d0db0fbdc..d2fcc93e0 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -38,10 +39,10 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("do_field_ionization", do_field_ionization); - if (do_field_ionization){ - AddRealComp("ionization_level"); - plot_flags.resize(PIdx::nattribs + 1, 1); - } + // If do_field_ionization, read initialization data from input file and + // read ionization energies from table. + if (do_field_ionization) + InitIonizationModule(); pp.query("plot_species", plot_species); int do_user_plot_vars; @@ -79,10 +80,12 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) : WarpXParticleContainer(amr_core, 0) { plasma_injector.reset(new PlasmaInjector()); +/* if (do_field_ionization){ Print()<<"AddRealComp\n"; AddRealComp("ionization_level"); } +*/ } void PhysicalParticleContainer::InitData() @@ -1604,25 +1607,38 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, copy_attribs(pti, x, y, z); } + Real* AMREX_RESTRICT ion_lev = nullptr; + if (do_field_ionization){ + ion_lev = pti.GetAttribs(particle_comps["ionization_level"]).dataPtr(); + } + // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( pti.numParticles(), + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } UpdateMomentumBoris( ux[i], uy[i], uz[i], gi[i], - Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); UpdatePosition( x[i], y[i], z[i], ux[i], uy[i], uz[i], dt ); } ); } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) { - amrex::ParallelFor( pti.numParticles(), + amrex::ParallelFor( + pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { + Real qp = q; + if (ion_lev){ qp *= ion_lev[i]; } UpdateMomentumVay( ux[i], uy[i], uz[i], gi[i], - Ex[i], Ey[i], Ez[i], Bx[i], By[i], Bz[i], q, m, dt); + Ex[i], Ey[i], Ez[i], Bx[i], + By[i], Bz[i], qp, m, dt); UpdatePosition( x[i], y[i], z[i], - ux[i], uy[i], uz[i], dt ); + ux[i], uy[i], uz[i], dt ); } ); } else { @@ -2058,3 +2074,35 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } } } + + +void PhysicalParticleContainer::InitIonizationModule() +{ + ParmParse pp(species_name); + pp.query("species_ionization_level", species_ionization_level); + pp.get("ionization_product", ionization_product_name); + pp.get("physical_element", physical_element); + // Add Real component for ionization level + AddRealComp("ionization_level"); + plot_flags.resize(PIdx::nattribs + 1, 1); + // Get atomic number and ionization energies from file + int ion_element_id = ion_map_ids[physical_element]; + ion_atomic_number = ion_atomic_numbers[ion_element_id]; + ionization_energies.resize(ion_atomic_number); + int offset = ion_energy_offsets[ion_element_id]; + for(int i=0; i particle_comps; @@ -290,8 +292,18 @@ protected: // This is currently required because continuous injection does not // support all features allowed by direct injection. int do_continuous_injection = 0; - int do_field_ionization = 0; + int do_field_ionization = 0; + int ionization_product; + std::string ionization_product_name; + int ion_atomic_number; + int species_ionization_level = 0; + amrex::Gpu::ManagedVector ionization_energies; + amrex::Gpu::ManagedVector adk_power; + amrex::Gpu::ManagedVector adk_prefactor; + amrex::Gpu::ManagedVector adk_exp_prefactor; + std::string physical_element; + int do_boosted_frame_diags = 1; amrex::Vector local_rho; -- cgit v1.2.3 From 43c06604530ac2c532d1c731fd60bbb5570446e1 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sat, 3 Aug 2019 17:56:59 -0700 Subject: add ionization table data --- Source/Particles/PhysicalParticleContainer.H | 18 -- Source/Particles/PhysicalParticleContainer.cpp | 2 +- Source/Utils/IonizationEnergiesTable.H | 136 ++++++++ Source/Utils/Make.package | 1 + Source/Utils/WarpXConst.H | 25 +- Source/Utils/atomic_data.txt | 427 +++++++++++++++++++++++++ Source/Utils/write_atomic_data_cpp.py | 67 ++++ 7 files changed, 647 insertions(+), 29 deletions(-) create mode 100644 Source/Utils/IonizationEnergiesTable.H create mode 100644 Source/Utils/atomic_data.txt create mode 100644 Source/Utils/write_atomic_data_cpp.py (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ddf7d1253..a6cf9d81f 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -108,25 +108,7 @@ public: virtual void PostRestart () final {} void SplitParticles(int lev); - - virtual void copyParticles(int lev) override; - virtual int doFieldIonization( - const int lev, - amrex::Cuda::ManagedDeviceVector& x_buf, - amrex::Cuda::ManagedDeviceVector& y_buf, - amrex::Cuda::ManagedDeviceVector& z_buf, - amrex::Cuda::ManagedDeviceVector& ux_buf, - amrex::Cuda::ManagedDeviceVector& uy_buf, - amrex::Cuda::ManagedDeviceVector& uz_buf, - amrex::Cuda::ManagedDeviceVector& ex_buf, - amrex::Cuda::ManagedDeviceVector& ey_buf, - amrex::Cuda::ManagedDeviceVector& ez_buf, - amrex::Cuda::ManagedDeviceVector& bx_buf, - amrex::Cuda::ManagedDeviceVector& by_buf, - amrex::Cuda::ManagedDeviceVector& bz_buf, - amrex::Cuda::ManagedDeviceVector& w_buf) override; - // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d2fcc93e0..0183a44b9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2079,7 +2079,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, void PhysicalParticleContainer::InitIonizationModule() { ParmParse pp(species_name); - pp.query("species_ionization_level", species_ionization_level); + pp.query("ionization_level", species_ionization_level); pp.get("ionization_product", ionization_product_name); pp.get("physical_element", physical_element); // Add Real component for ionization level diff --git a/Source/Utils/IonizationEnergiesTable.H b/Source/Utils/IonizationEnergiesTable.H new file mode 100644 index 000000000..475236abd --- /dev/null +++ b/Source/Utils/IonizationEnergiesTable.H @@ -0,0 +1,136 @@ +#include +#include + +#ifndef WARPX_IONIZATION_TABLE_H_ +#define WARPX_IONIZATION_TABLE_H_ + +std::map ion_map_ids = { + {"H", 0}, + {"He", 1}, + {"Li", 2}, + {"Be", 3}, + {"B", 4}, + {"C", 5}, + {"N", 6}, + {"O", 7}, + {"F", 8}, + {"Ne", 9}, + {"Na", 10}, + {"Mg", 11}, + {"Al", 12}, + {"Si", 13}, + {"P", 14}, + {"S", 15}, + {"Cl", 16}, + {"Ar", 17}, + {"Kr", 18}, + {"Rb", 19}, + {"Xe", 20}, + {"Rn", 21} }; + +const int nelements = 22; + +const int ion_atomic_numbers[nelements] = { + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, + 11, 12, 13, 14, 15, 16, 17, 18, 36, 37, + 54, 86}; + +const int ion_energy_offsets[nelements] = { + 0, 1, 3, 6, 10, 15, 21, 28, 36, 45, + 55, 66, 78, 91, 105, 120, 136, 153, 171, 207, + 244, 298}; + +const int energies_tab_length = 384; + +const amrex::Real table_ionization_energies[energies_tab_length]{ + // H + 13.59843449, + // He + 24.58738880, 54.4177650, + // Li + 5.39171495, 75.6400964, 122.4543581, + // Be + 9.322699, 18.21115, 153.896203, 217.7185843, + // B + 8.298019, 25.15483, 37.93058, 259.3715, 340.226020, + // C + 11.2602880, 24.383154, 47.88778, 64.49352, 392.090515, 489.993194, + // N + 14.53413, 29.60125, 47.4453, 77.4735, 97.8901, 552.06732, 667.046116, + // O + 13.618055, 35.12112, 54.93554, 77.41350, 113.8990, 138.1189, 739.32682, + 871.40988, + // F + 17.42282, 34.97081, 62.70798, 87.175, 114.249, 157.16311, 185.1868, + 953.89804, 1103.11747, + // Ne + 21.564540, 40.96297, 63.4233, 97.1900, 126.247, 157.934, 207.271, + 239.0970, 1195.80783, 1362.19915, + // Na + 5.1390769, 47.28636, 71.6200, 98.936, 138.404, 172.23, 208.504, + 264.192, 299.856, 1465.13449, 1648.70218, + // Mg + 7.646236, 15.035271, 80.1436, 109.2654, 141.33, 186.76, 225.02, + 265.924, 327.99, 367.489, 1761.80487, 1962.66365, + // Al + 5.985769, 18.82855, 28.447642, 119.9924, 153.8252, 190.49, 241.76, + 284.64, 330.21, 398.65, 442.005, 2085.97700, 2304.14005, + // Si + 8.15168, 16.34585, 33.49300, 45.14179, 166.767, 205.279, 246.57, + 303.59, 351.28, 401.38, 476.273, 523.415, 2437.65813, 2673.17753, + // P + 10.486686, 19.76949, 30.20264, 51.44387, 65.02511, 220.430, 263.57, + 309.60, 372.31, 424.40, 479.44, 560.62, 611.741, 2816.90876, + 3069.8415, + // S + 10.36001, 23.33788, 34.86, 47.222, 72.5945, 88.0529, 280.954, + 328.794, 379.84, 447.7, 504.55, 564.41, 651.96, 706.994, + 3223.7807, 3494.1879, + // Cl + 12.967632, 23.81364, 39.80, 53.24, 67.68, 96.94, 114.2013, + 348.306, 400.851, 456.7, 530.0, 591.58, 656.30, 750.23, + 809.198, 3658.3437, 3946.2909, + // Ar + 15.7596117, 27.62967, 40.735, 59.58, 74.84, 91.290, 124.41, + 143.4567, 422.60, 479.76, 540.4, 619.0, 685.5, 755.13, + 855.5, 918.375, 4120.6656, 4426.2228, + // Kr + 13.9996053, 24.35984, 35.838, 50.85, 64.69, 78.49, 109.13, + 125.802, 233.0, 268, 308, 350, 391, 446, + 492, 540, 591, 640, 785, 831.6, 882.8, + 945, 999.0, 1042, 1155.0, 1205.23, 2928.9, 3072, + 3228, 3380, 3584, 3752.0, 3971, 4109.083, 17296.420, + 17936.209, + // Rb + 4.1771280, 27.28954, 39.247, 52.20, 68.44, 82.9, 98.67, + 132.79, 150.628, 277.12, 313.1, 356.0, 400, 443, + 502, 550, 601, 654, 706.0, 857, 905.3, + 958.9, 1024, 1080, 1125, 1242.5, 1294.57, 3133.3, + 3281, 3443, 3600, 3815, 3988, 4214, 4356.865, + 18305.884, 18965.516, + // Xe + 12.1298436, 20.975, 31.05, 42.20, 54.1, 66.703, 91.6, + 105.9778, 179.84, 202.0, 229.02, 255.0, 281, 314, + 343, 374, 404, 434, 549, 582, 616, + 650, 700, 736, 818, 857.0, 1493, 1571, + 1653, 1742, 1826, 1919, 2023, 2113, 2209, + 2300, 2556, 2637, 2726, 2811, 2975, 3068, + 3243, 3333.8, 7660, 7889, 8144, 8382, 8971, + 9243, 9581, 9810.37, 40271.724, 41299.71, + // Rn + 10.74850, 21.4, 29.4, 36.9, 52.9, 64.0, 88.0, + 102.0, 154.0, 173.9, 195.0, 218.0, 240, 264, + 293, 317, 342, 367, 488, 520, 550, + 580, 640, 680, 760, 800, 850, 920, + 980, 1050, 1110, 1180, 1250, 1310, 1390, + 1460, 1520, 1590, 1660, 1720, 2033, 2094, + 2158, 2227, 2293, 2357, 2467, 2535, 2606, + 2674, 2944, 3010, 3082, 3149, 3433, 3510, + 3699, 3777, 6169, 6318, 6476, 6646, 6807, + 6964, 7283, 7450, 7630, 7800, 8260, 8410, + 8570, 8710, 9610, 9780, 10120, 10290, 21770, + 22160, 22600, 22990, 26310, 26830, 27490, 27903.1, + 110842.0, 112843.7 }; + +#endif // #ifndef WARPX_IONIZATION_TABLE_H_ + diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index cd335dbcf..389b3670a 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -7,6 +7,7 @@ CEXE_headers += WarpXUtil.H CEXE_headers += WarpXAlgorithmSelection.H CEXE_sources += WarpXAlgorithmSelection.cpp CEXE_headers += NCIGodfreyTables.H +CEXE_headers += IonizationEnergiesTable.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index 4e226e55f..bfb67b3d5 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -3,20 +3,25 @@ #include -// Physical constants -namespace PhysConst +// Math constants +namespace MathConst { - static constexpr amrex::Real c = 299792458.; - static constexpr amrex::Real ep0 = 8.854187817e-12; - static constexpr amrex::Real mu0 = 1.2566370614359173e-06; - static constexpr amrex::Real q_e = 1.6021764620000001e-19; - static constexpr amrex::Real m_e = 9.10938291e-31; - static constexpr amrex::Real m_p = 1.6726231000000001e-27; + static constexpr amrex::Real pi = 3.14159265358979323846; } -namespace MathConst +// Physical constants +namespace PhysConst { - static constexpr amrex::Real pi = 3.14159265358979323846; + static constexpr amrex::Real c = 299792458.; + static constexpr amrex::Real ep0 = 8.854187817e-12; + static constexpr amrex::Real mu0 = 1.2566370614359173e-06; + static constexpr amrex::Real q_e = 1.6021764620000001e-19; + static constexpr amrex::Real m_e = 9.10938291e-31; + static constexpr amrex::Real m_p = 1.6726231000000001e-27; + static constexpr amrex::Real hbar = 1.05457180010e-34; + static constexpr amrex::Real alpha = mu0/(4*MathConst::pi)*q_e*q_e*c/hbar; + static constexpr amrex::Real r_e = 1./(4*MathConst::pi*ep0) * q_e*q_e/(m_e*c*c); + // static constexpr amrex::Real alpha = 0.0072973525693; } #endif diff --git a/Source/Utils/atomic_data.txt b/Source/Utils/atomic_data.txt new file mode 100644 index 000000000..140f5a26a --- /dev/null +++ b/Source/Utils/atomic_data.txt @@ -0,0 +1,427 @@ +# Reference: +# Kramida, A., Ralchenko, Yu., Reader, J., and NIST ASD Team (2014). +# NIST Atomic Spectra Database (ver. 5.2), [Online]. +# Available: http://physics.nist.gov/asd [2017, March 3]. + +# License (from https://www.nist.gov/director/licensing): +# This data was developed by employees of the National Institute of Standards +# and Technology (NIST), an agency of the Federal Government. Pursuant to title +# 17 United States Code Section 105, works of NIST employees are not subject to +# copyright protection in the United States and are considered to be in the +# public domain. +# The data is provided by NIST as a public service and is expressly provided +# "AS IS." NIST MAKES NO WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR STATUTORY, +# INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTY OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE, NON-INFRINGEMENT AND DATA ACCURACY. NIST +# does not warrant or make any representations regarding the use of the data or +# the results thereof, including but not limited to the correctness, accuracy, +# reliability or usefulness of the data. NIST SHALL NOT BE LIABLE AND YOU +# HEREBY RELEASE NIST FROM LIABILITY FOR ANY INDIRECT, CONSEQUENTIAL, SPECIAL, +# OR INCIDENTAL DAMAGES (INCLUDING DAMAGES FOR LOSS OF BUSINESS PROFITS, +# BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, AND THE LIKE), WHETHER +# ARISING IN TORT, CONTRACT, OR OTHERWISE, ARISING FROM OR RELATING TO THE DATA +# (OR THE USE OF OR INABILITY TO USE THIS DATA), EVEN IF NIST HAS BEEN ADVISED +# OF THE POSSIBILITY OF SUCH DAMAGES. +# To the extent that NIST may hold copyright in countries other than the United +# States, you are hereby granted the non-exclusive irrevocable and unconditional +# right to print, publish, prepare derivative works and distribute the NIST +# data, in any medium, or authorize others to do so on your behalf, on a +# royalty-free basis throughout the world. +# You may improve, modify, and create derivative works of the data or any +# portion of the data, and you may copy and distribute such modifications or +# works. Modified works should carry a notice stating that you changed the data +# and should note the date and nature of any such change. Please explicitly +# acknowledge the National Institute of Standards and Technology as the source +# of the data: Data citation recommendations are provided below. +# Permission to use this data is contingent upon your acceptance of the terms +# of this agreement and upon your providing appropriate acknowledgments of +# NIST's creation of the data. + +----------------------------------------------------------------------------- +At. num | Sp. Name | Ion Charge | Ionization Energy (eV) | +--------|---------------|------------|--------------------------------------| + 1 | H I | 0 | (13.59843449) | + 2 | He I | 0 | 24.58738880 | + 2 | He II | +1 | (54.4177650) | + 3 | Li I | 0 | 5.39171495 | + 3 | Li II | +1 | [75.6400964] | + 3 | Li III | +2 | (122.4543581) | + 4 | Be I | 0 | 9.322699 | + 4 | Be II | +1 | 18.21115 | + 4 | Be III | +2 | [153.896203] | + 4 | Be IV | +3 | (217.7185843) | + 5 | B I | 0 | 8.298019 | + 5 | B II | +1 | 25.15483 | + 5 | B III | +2 | 37.93058 | + 5 | B IV | +3 | 259.3715 | + 5 | B V | +4 | (340.226020) | + 6 | C I | 0 | 11.2602880 | + 6 | C II | +1 | 24.383154 | + 6 | C III | +2 | 47.88778 | + 6 | C IV | +3 | 64.49352 | + 6 | C V | +4 | [392.090515] | + 6 | C VI | +5 | (489.993194) | + 7 | N I | 0 | 14.53413 | + 7 | N II | +1 | [29.60125] | + 7 | N III | +2 | [47.4453] | + 7 | N IV | +3 | 77.4735 | + 7 | N V | +4 | 97.8901 | + 7 | N VI | +5 | [552.06732] | + 7 | N VII | +6 | (667.046116) | + 8 | O I | 0 | 13.618055 | + 8 | O II | +1 | 35.12112 | + 8 | O III | +2 | [54.93554] | + 8 | O IV | +3 | 77.41350 | + 8 | O V | +4 | 113.8990 | + 8 | O VI | +5 | [138.1189] | + 8 | O VII | +6 | [739.32682] | + 8 | O VIII | +7 | (871.40988) | + 9 | F I | 0 | 17.42282 | + 9 | F II | +1 | 34.97081 | + 9 | F III | +2 | [62.70798] | + 9 | F IV | +3 | [87.175] | + 9 | F V | +4 | [114.249] | + 9 | F VI | +5 | 157.16311 | + 9 | F VII | +6 | [185.1868] | + 9 | F VIII | +7 | [953.89804] | + 9 | F IX | +8 | (1103.11747) | + 10 | Ne I | 0 | 21.564540 | + 10 | Ne II | +1 | 40.96297 | + 10 | Ne III | +2 | 63.4233 | + 10 | Ne IV | +3 | [97.1900] | + 10 | Ne V | +4 | 126.247 | + 10 | Ne VI | +5 | 157.934 | + 10 | Ne VII | +6 | 207.271 | + 10 | Ne VIII | +7 | [239.0970] | + 10 | Ne IX | +8 | [1195.80783] | + 10 | Ne X | +9 | (1362.19915) | + 11 | Na I | 0 | 5.1390769 | + 11 | Na II | +1 | 47.28636 | + 11 | Na III | +2 | [71.6200] | + 11 | Na IV | +3 | [98.936] | + 11 | Na V | +4 | 138.404 | + 11 | Na VI | +5 | [172.23] | + 11 | Na VII | +6 | [208.504] | + 11 | Na VIII | +7 | [264.192] | + 11 | Na IX | +8 | (299.856) | + 11 | Na X | +9 | (1465.13449) | + 11 | Na XI | +10 | (1648.70218) | + 12 | Mg I | 0 | 7.646236 | + 12 | Mg II | +1 | 15.035271 | + 12 | Mg III | +2 | 80.1436 | + 12 | Mg IV | +3 | [109.2654] | + 12 | Mg V | +4 | [141.33] | + 12 | Mg VI | +5 | [186.76] | + 12 | Mg VII | +6 | [225.02] | + 12 | Mg VIII | +7 | [265.924] | + 12 | Mg IX | +8 | [327.99] | + 12 | Mg X | +9 | (367.489) | + 12 | Mg XI | +10 | (1761.80487) | + 12 | Mg XII | +11 | (1962.66365) | + 13 | Al I | 0 | 5.985769 | + 13 | Al II | +1 | 18.82855 | + 13 | Al III | +2 | 28.447642 | + 13 | Al IV | +3 | 119.9924 | + 13 | Al V | +4 | [153.8252] | + 13 | Al VI | +5 | [190.49] | + 13 | Al VII | +6 | [241.76] | + 13 | Al VIII | +7 | [284.64] | + 13 | Al IX | +8 | [330.21] | + 13 | Al X | +9 | [398.65] | + 13 | Al XI | +10 | (442.005) | + 13 | Al XII | +11 | (2085.97700) | + 13 | Al XIII | +12 | (2304.14005) | + 14 | Si I | 0 | 8.15168 | + 14 | Si II | +1 | 16.34585 | + 14 | Si III | +2 | 33.49300 | + 14 | Si IV | +3 | 45.14179 | + 14 | Si V | +4 | 166.767 | + 14 | Si VI | +5 | [205.279] | + 14 | Si VII | +6 | [246.57] | + 14 | Si VIII | +7 | [303.59] | + 14 | Si IX | +8 | [351.28] | + 14 | Si X | +9 | [401.38] | + 14 | Si XI | +10 | [476.273] | + 14 | Si XII | +11 | (523.415) | + 14 | Si XIII | +12 | (2437.65813) | + 14 | Si XIV | +13 | (2673.17753) | + 15 | P I | 0 | 10.486686 | + 15 | P II | +1 | [19.76949] | + 15 | P III | +2 | 30.20264 | + 15 | P IV | +3 | 51.44387 | + 15 | P V | +4 | 65.02511 | + 15 | P VI | +5 | [220.430] | + 15 | P VII | +6 | [263.57] | + 15 | P VIII | +7 | [309.60] | + 15 | P IX | +8 | [372.31] | + 15 | P X | +9 | [424.40] | + 15 | P XI | +10 | [479.44] | + 15 | P XII | +11 | [560.62] | + 15 | P XIII | +12 | (611.741) | + 15 | P XIV | +13 | (2816.90876) | + 15 | P XV | +14 | (3069.8415) | + 16 | S I | 0 | 10.36001 | + 16 | S II | +1 | 23.33788 | + 16 | S III | +2 | [34.86] | + 16 | S IV | +3 | [47.222] | + 16 | S V | +4 | 72.5945 | + 16 | S VI | +5 | 88.0529 | + 16 | S VII | +6 | 280.954 | + 16 | S VIII | +7 | 328.794 | + 16 | S IX | +8 | [379.84] | + 16 | S X | +9 | [447.7] | + 16 | S XI | +10 | [504.55] | + 16 | S XII | +11 | [564.41] | + 16 | S XIII | +12 | [651.96] | + 16 | S XIV | +13 | (706.994) | + 16 | S XV | +14 | (3223.7807) | + 16 | S XVI | +15 | (3494.1879) | + 17 | Cl I | 0 | 12.967632 | + 17 | Cl II | +1 | 23.81364 | + 17 | Cl III | +2 | [39.80] | + 17 | Cl IV | +3 | [53.24] | + 17 | Cl V | +4 | [67.68] | + 17 | Cl VI | +5 | [96.94] | + 17 | Cl VII | +6 | 114.2013 | + 17 | Cl VIII | +7 | 348.306 | + 17 | Cl IX | +8 | 400.851 | + 17 | Cl X | +9 | [456.7] | + 17 | Cl XI | +10 | [530.0] | + 17 | Cl XII | +11 | [591.58] | + 17 | Cl XIII | +12 | [656.30] | + 17 | Cl XIV | +13 | [750.23] | + 17 | Cl XV | +14 | (809.198) | + 17 | Cl XVI | +15 | (3658.3437) | + 17 | Cl XVII | +16 | (3946.2909) | + 18 | Ar I | 0 | 15.7596117 | + 18 | Ar II | +1 | 27.62967 | + 18 | Ar III | +2 | [40.735] | + 18 | Ar IV | +3 | [59.58] | + 18 | Ar V | +4 | [74.84] | + 18 | Ar VI | +5 | 91.290 | + 18 | Ar VII | +6 | [124.41] | + 18 | Ar VIII | +7 | [143.4567] | + 18 | Ar IX | +8 | [422.60] | + 18 | Ar X | +9 | [479.76] | + 18 | Ar XI | +10 | [540.4] | + 18 | Ar XII | +11 | [619.0] | + 18 | Ar XIII | +12 | [685.5] | + 18 | Ar XIV | +13 | [755.13] | + 18 | Ar XV | +14 | [855.5] | + 18 | Ar XVI | +15 | (918.375) | + 18 | Ar XVII | +16 | (4120.6656) | + 18 | Ar XVIII | +17 | (4426.2228) | + 36 | Kr I | 0 | 13.9996053 | + 36 | Kr II | +1 | 24.35984 | + 36 | Kr III | +2 | 35.838 | + 36 | Kr IV | +3 | 50.85 | + 36 | Kr V | +4 | [64.69] | + 36 | Kr VI | +5 | [78.49] | + 36 | Kr VII | +6 | (109.13) | + 36 | Kr VIII | +7 | 125.802 | + 36 | Kr IX | +8 | [233.0] | + 36 | Kr X | +9 | (268) | + 36 | Kr XI | +10 | (308) | + 36 | Kr XII | +11 | (350) | + 36 | Kr XIII | +12 | (391) | + 36 | Kr XIV | +13 | (446) | + 36 | Kr XV | +14 | (492) | + 36 | Kr XVI | +15 | (540) | + 36 | Kr XVII | +16 | (591) | + 36 | Kr XVIII | +17 | (640) | + 36 | Kr XIX | +18 | [785] | + 36 | Kr XX | +19 | [831.6] | + 36 | Kr XXI | +20 | [882.8] | + 36 | Kr XXII | +21 | [945] | + 36 | Kr XXIII | +22 | [999.0] | + 36 | Kr XXIV | +23 | [1042] | + 36 | Kr XXV | +24 | [1155.0] | + 36 | Kr XXVI | +25 | [1205.23] | + 36 | Kr XXVII | +26 | [2928.9] | + 36 | Kr XXVIII | +27 | [3072] | + 36 | Kr XXIX | +28 | [3228] | + 36 | Kr XXX | +29 | [3380] | + 36 | Kr XXXI | +30 | [3584] | + 36 | Kr XXXII | +31 | [3752.0] | + 36 | Kr XXXIII | +32 | [3971] | + 36 | Kr XXXIV | +33 | (4109.083) | + 36 | Kr XXXV | +34 | (17296.420) | + 36 | Kr XXXVI | +35 | (17936.209) | + 37 | Rb I | 0 | 4.1771280 | + 37 | Rb II | +1 | 27.28954 | + 37 | Rb III | +2 | [39.247] | + 37 | Rb IV | +3 | [52.20] | + 37 | Rb V | +4 | (68.44) | + 37 | Rb VI | +5 | (82.9) | + 37 | Rb VII | +6 | [98.67] | + 37 | Rb VIII | +7 | [132.79] | + 37 | Rb IX | +8 | 150.628 | + 37 | Rb X | +9 | 277.12 | + 37 | Rb XI | +10 | (313.1) | + 37 | Rb XII | +11 | (356.0) | + 37 | Rb XIII | +12 | (400) | + 37 | Rb XIV | +13 | (443) | + 37 | Rb XV | +14 | (502) | + 37 | Rb XVI | +15 | (550) | + 37 | Rb XVII | +16 | (601) | + 37 | Rb XVIII | +17 | (654) | + 37 | Rb XIX | +18 | (706.0) | + 37 | Rb XX | +19 | [857] | + 37 | Rb XXI | +20 | [905.3] | + 37 | Rb XXII | +21 | [958.9] | + 37 | Rb XXIII | +22 | [1024] | + 37 | Rb XXIV | +23 | [1080] | + 37 | Rb XXV | +24 | [1125] | + 37 | Rb XXVI | +25 | [1242.5] | + 37 | Rb XXVII | +26 | [1294.57] | + 37 | Rb XXVIII | +27 | [3133.3] | + 37 | Rb XXIX | +28 | [3281] | + 37 | Rb XXX | +29 | [3443] | + 37 | Rb XXXI | +30 | [3600] | + 37 | Rb XXXII | +31 | [3815] | + 37 | Rb XXXIII | +32 | [3988] | + 37 | Rb XXXIV | +33 | [4214] | + 37 | Rb XXXV | +34 | (4356.865) | + 37 | Rb XXXVI | +35 | (18305.884) | + 37 | Rb XXXVII | +36 | (18965.516) | + 54 | Xe I | 0 | 12.1298436 | + 54 | Xe II | +1 | [20.975] | + 54 | Xe III | +2 | [31.05] | + 54 | Xe IV | +3 | 42.20 | + 54 | Xe V | +4 | [54.1] | + 54 | Xe VI | +5 | 66.703 | + 54 | Xe VII | +6 | [91.6] | + 54 | Xe VIII | +7 | 105.9778 | + 54 | Xe IX | +8 | 179.84 | + 54 | Xe X | +9 | (202.0) | + 54 | Xe XI | +10 | [229.02] | + 54 | Xe XII | +11 | (255.0) | + 54 | Xe XIII | +12 | (281) | + 54 | Xe XIV | +13 | (314) | + 54 | Xe XV | +14 | (343) | + 54 | Xe XVI | +15 | (374) | + 54 | Xe XVII | +16 | (404) | + 54 | Xe XVIII | +17 | (434) | + 54 | Xe XIX | +18 | (549) | + 54 | Xe XX | +19 | (582) | + 54 | Xe XXI | +20 | (616) | + 54 | Xe XXII | +21 | (650) | + 54 | Xe XXIII | +22 | (700) | + 54 | Xe XXIV | +23 | (736) | + 54 | Xe XXV | +24 | (818) | + 54 | Xe XXVI | +25 | [857.0] | + 54 | Xe XXVII | +26 | (1493) | + 54 | Xe XXVIII | +27 | (1571) | + 54 | Xe XXIX | +28 | (1653) | + 54 | Xe XXX | +29 | (1742) | + 54 | Xe XXXI | +30 | (1826) | + 54 | Xe XXXII | +31 | (1919) | + 54 | Xe XXXIII | +32 | (2023) | + 54 | Xe XXXIV | +33 | (2113) | + 54 | Xe XXXV | +34 | (2209) | + 54 | Xe XXXVI | +35 | (2300) | + 54 | Xe XXXVII | +36 | (2556) | + 54 | Xe XXXVIII | +37 | (2637) | + 54 | Xe XXXIX | +38 | (2726) | + 54 | Xe XL | +39 | (2811) | + 54 | Xe XLI | +40 | (2975) | + 54 | Xe XLII | +41 | (3068) | + 54 | Xe XLIII | +42 | (3243) | + 54 | Xe XLIV | +43 | [3333.8] | + 54 | Xe XLV | +44 | (7660) | + 54 | Xe XLVI | +45 | (7889) | + 54 | Xe XLVII | +46 | (8144) | + 54 | Xe XLVIII | +47 | (8382) | + 54 | Xe XLIX | +48 | (8971) | + 54 | Xe L | +49 | (9243) | + 54 | Xe LI | +50 | (9581) | + 54 | Xe LII | +51 | (9810.37) | + 54 | Xe LIII | +52 | (40271.724) | + 54 | Xe LIV | +53 | (41299.71) | + 86 | Rn I | 0 | 10.74850 | + 86 | Rn II | +1 | [21.4] | + 86 | Rn III | +2 | [29.4] | + 86 | Rn IV | +3 | (36.9) | + 86 | Rn V | +4 | (52.9) | + 86 | Rn VI | +5 | (64.0) | + 86 | Rn VII | +6 | (88.0) | + 86 | Rn VIII | +7 | (102.0) | + 86 | Rn IX | +8 | (154.0) | + 86 | Rn X | +9 | (173.9) | + 86 | Rn XI | +10 | (195.0) | + 86 | Rn XII | +11 | (218.0) | + 86 | Rn XIII | +12 | (240) | + 86 | Rn XIV | +13 | (264) | + 86 | Rn XV | +14 | (293) | + 86 | Rn XVI | +15 | (317) | + 86 | Rn XVII | +16 | (342) | + 86 | Rn XVIII | +17 | (367) | + 86 | Rn XIX | +18 | (488) | + 86 | Rn XX | +19 | (520) | + 86 | Rn XXI | +20 | (550) | + 86 | Rn XXII | +21 | (580) | + 86 | Rn XXIII | +22 | (640) | + 86 | Rn XXIV | +23 | (680) | + 86 | Rn XXV | +24 | (760) | + 86 | Rn XXVI | +25 | (800) | + 86 | Rn XXVII | +26 | (850) | + 86 | Rn XXVIII | +27 | (920) | + 86 | Rn XXIX | +28 | (980) | + 86 | Rn XXX | +29 | (1050) | + 86 | Rn XXXI | +30 | (1110) | + 86 | Rn XXXII | +31 | (1180) | + 86 | Rn XXXIII | +32 | (1250) | + 86 | Rn XXXIV | +33 | (1310) | + 86 | Rn XXXV | +34 | (1390) | + 86 | Rn XXXVI | +35 | (1460) | + 86 | Rn XXXVII | +36 | (1520) | + 86 | Rn XXXVIII | +37 | (1590) | + 86 | Rn XXXIX | +38 | (1660) | + 86 | Rn XL | +39 | (1720) | + 86 | Rn XLI | +40 | (2033) | + 86 | Rn XLII | +41 | (2094) | + 86 | Rn XLIII | +42 | (2158) | + 86 | Rn XLIV | +43 | (2227) | + 86 | Rn XLV | +44 | (2293) | + 86 | Rn XLVI | +45 | (2357) | + 86 | Rn XLVII | +46 | (2467) | + 86 | Rn XLVIII | +47 | (2535) | + 86 | Rn XLIX | +48 | (2606) | + 86 | Rn L | +49 | (2674) | + 86 | Rn LI | +50 | (2944) | + 86 | Rn LII | +51 | (3010) | + 86 | Rn LIII | +52 | (3082) | + 86 | Rn LIV | +53 | (3149) | + 86 | Rn LV | +54 | (3433) | + 86 | Rn LVI | +55 | (3510) | + 86 | Rn LVII | +56 | (3699) | + 86 | Rn LVIII | +57 | [3777] | + 86 | Rn LIX | +58 | (6169) | + 86 | Rn LX | +59 | (6318) | + 86 | Rn LXI | +60 | (6476) | + 86 | Rn LXII | +61 | (6646) | + 86 | Rn LXIII | +62 | (6807) | + 86 | Rn LXIV | +63 | (6964) | + 86 | Rn LXV | +64 | (7283) | + 86 | Rn LXVI | +65 | (7450) | + 86 | Rn LXVII | +66 | (7630) | + 86 | Rn LXVIII | +67 | (7800) | + 86 | Rn LXIX | +68 | (8260) | + 86 | Rn LXX | +69 | (8410) | + 86 | Rn LXXI | +70 | (8570) | + 86 | Rn LXXII | +71 | (8710) | + 86 | Rn LXXIII | +72 | (9610) | + 86 | Rn LXXIV | +73 | (9780) | + 86 | Rn LXXV | +74 | (10120) | + 86 | Rn LXXVI | +75 | (10290) | + 86 | Rn LXXVII | +76 | (21770) | + 86 | Rn LXXVIII | +77 | (22160) | + 86 | Rn LXXIX | +78 | (22600) | + 86 | Rn LXXX | +79 | (22990) | + 86 | Rn LXXXI | +80 | (26310) | + 86 | Rn LXXXII | +81 | (26830) | + 86 | Rn LXXXIII | +82 | (27490) | + 86 | Rn LXXXIV | +83 | (27903.1) | + 86 | Rn LXXXV | +84 | (110842.0) | + 86 | Rn LXXXVI | +85 | (112843.7) | +----------------------------------------------------------------------------- diff --git a/Source/Utils/write_atomic_data_cpp.py b/Source/Utils/write_atomic_data_cpp.py new file mode 100644 index 000000000..21d61a075 --- /dev/null +++ b/Source/Utils/write_atomic_data_cpp.py @@ -0,0 +1,67 @@ +import re, os +import numpy as np +from scipy.constants import e + +filename = os.path.join( '.', 'atomic_data.txt' ) +with open(filename) as f: + text_data = f.read() + +# Read full table from file and get names, atomic numbers and offsets +# position in table of ionization energies for all species +regex_command = '\n\s+(\d+)\s+\|\s+([A-Z]+[a-z]*)\s+\w+\s+\|\s+\+*(\d+)\s+\|\s+\(*\[*(\d+\.*\d*)' +list_of_tuples = re.findall( regex_command, text_data ) +ion_atom_numbers = [int(i) for i in list(dict.fromkeys( [x[0] for x in list_of_tuples] ))] +ion_names = list(dict.fromkeys( [x[1] for x in list_of_tuples] )) +ion_offsets = np.concatenate(([0], np.cumsum(np.array(ion_atom_numbers)[:-1])), axis=0) + +# Head of CPP file +cpp_string = '' +cpp_string += '#include \n' +cpp_string += '#include \n\n' +cpp_string += '#ifndef WARPX_IONIZATION_TABLE_H_\n' +cpp_string += '#define WARPX_IONIZATION_TABLE_H_\n\n' + +# Map each element to ID in table +cpp_string += 'std::map ion_map_ids = {' +for count, name in enumerate(ion_names): + cpp_string += '\n {"' + name + '", ' + str(count) + '},' +cpp_string = cpp_string[:-1] +cpp_string += ' };\n\n' + +# Atomic number of each species +cpp_string += 'const int nelements = ' + str(len(ion_names)) + ';\n\n' +cpp_string += 'const int ion_atomic_numbers[nelements] = {' +for count, atom_num in enumerate(ion_atom_numbers): + if count%10==0: cpp_string += '\n ' + cpp_string += str(atom_num) + ', ' +cpp_string = cpp_string[:-2] +cpp_string += '};\n\n' + +# Offset of each element in table of ionization energies +cpp_string += 'const int ion_energy_offsets[nelements] = {' +for count, offset in enumerate(ion_offsets): + if count%10==0: cpp_string += '\n ' + cpp_string += str(offset) + ', ' +cpp_string = cpp_string[:-2] +cpp_string += '};\n\n' + +# Table of ionization energies +cpp_string += 'const int energies_tab_length = ' + str(len(list_of_tuples)) + ';\n\n' +cpp_string += 'const amrex::Real table_ionization_energies[energies_tab_length]{' +for element in ion_names: + cpp_string += '\n // ' + element + regex_command = \ + '\n\s+(\d+)\s+\|\s+%s\s+\w+\s+\|\s+\+*(\d+)\s+\|\s+\(*\[*(\d+\.*\d*)' \ + %element + list_of_tuples = re.findall( regex_command, text_data ) + for count, energy in enumerate([x[2] for x in list_of_tuples]): + if count%7==0: cpp_string += '\n ' + cpp_string += energy + ', ' +cpp_string = cpp_string[:-2] +cpp_string += ' };\n\n' + +# Write the string to file +cpp_string += '#endif // #ifndef WARPX_IONIZATION_TABLE_H_\n' +f= open("IonizationEnergiesTable.H","w") +f.write(cpp_string) +f.close() -- cgit v1.2.3 From a0f579f27969c3f86e54874b5a3c6a732dff9f53 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 4 Aug 2019 15:47:12 -0700 Subject: add target species handling --- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 155 ++++++++++++++++++++++++- Source/Particles/PhysicalParticleContainer.H | 2 + Source/Particles/PhysicalParticleContainer.cpp | 21 +++- Source/Particles/WarpXParticleContainer.H | 5 +- 5 files changed, 179 insertions(+), 6 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 4ea9c2763..8f7fd56e3 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -7,7 +7,7 @@ #include #include - +#include #include #include #include diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 16c9b9762..736b906df 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -3,6 +3,9 @@ #include #include +#include +#include +#include #include #include #include @@ -508,7 +511,155 @@ MultiParticleContainer::getSpeciesID(std::string product_str) return i_product; } - void MultiParticleContainer::doFieldIonization() -{} +{ + /* + amrex::Gpu::SharedMemory grid_ids; + amrex::Gpu::SharedMemory tile_ids; + amrex::Gpu::SharedMemory is_ionized; + */ + + for (auto& pc : allcontainers){ + if (!pc->do_field_ionization){ continue; } + const Real * const AMREX_RESTRICT p_ionization_energies = pc->ionization_energies.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_prefactor = pc->adk_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_exp_prefactor = pc->adk_exp_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_power = pc->adk_power.dataPtr(); + for (int lev = 0; lev <= pc->finestLevel(); ++lev){ + +#ifdef _OPENMP + // First touch all tiles in the map in serial + for (MFIter mfi = pc->MakeMFIter(lev); mfi.isValid(); ++mfi) { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { + pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + } +#endif + MFItInfo info; + if (pc->do_tiling && Gpu::notInLaunchRegion()) { + info.EnableTiling(pc->tile_size); + } +#ifdef _OPENMP + info.SetDynamic(true); +#pragma omp parallel +#endif + for (MFIter mfi = pc->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + auto& particle_tile = pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + auto& soa = particle_tile.GetStructOfArrays(); + const int np = particle_tile.GetArrayOfStructs().size(); + if (np == 0) break; + amrex::Gpu::ManagedVector is_ionized; + // const int np_lev = pc->NumberOfParticlesAtLevel(lev); + is_ionized.resize(np); + int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); + const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); + const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); + const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); + const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); + const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); + const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); + const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); + const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); + Real * const AMREX_RESTRICT ilev_real = soa.GetRealData(pc->particle_comps["ionization_level"]).data(); + Real c = PhysConst::c; + Real c2_inv = 1./c/c; + + ParallelFor( + np, + [=] AMREX_GPU_DEVICE (long i) { + Real random_draw = Random(); + Real ga = std::sqrt(1. + + (ux[i]*ux[i] + + uy[i]*uy[i] + + uz[i]*uz[i]) * c2_inv); + Real E = std::sqrt( + - ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * c2_inv + + ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) * ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) + + ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) * ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) + + ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) * ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) + ); + int ilev = (int) round(ilev_real[i]); + // int ilev = static_cast(round(ilev_real[i])); + Real p; + Real w_dtau; + Print()< is_ionized_cumsum_vector; + is_ionized_cumsum_vector.resize(np); + int np_ionized = p_is_ionized[0]; + for(int i=1; iionization_product]; + auto& prod_particle_tile = prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + const int np_old = prod_particle_tile.GetArrayOfStructs().size(); + const int np_new = np_old + np_ionized; + prod_particle_tile.resize(np_new); + auto& prod_soa = prod_particle_tile.GetStructOfArrays(); + + GpuArray prod_pa; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + prod_pa[ia] = prod_soa.GetRealData(ia).data() + np_old; + } + WarpXParticleContainer::ParticleType* prod_pp = prod_particle_tile.GetArrayOfStructs()().data() + np_old; + WarpXParticleContainer::ParticleType* pp = particle_tile.GetArrayOfStructs()().data(); + + GpuArray pa; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + pa[ia] = soa.GetRealData(ia).data(); + } + + const int cpuid = ParallelDescriptor::MyProc(); + int* AMREX_RESTRICT is_ionized_cumsum = is_ionized_cumsum_vector.dataPtr(); + amrex::For( + np, [=] AMREX_GPU_DEVICE (int ip) noexcept + { + if(is_ionized[ip]){ + int i = is_ionized_cumsum[ip]; + WarpXParticleContainer::ParticleType& prod_p = prod_pp[i]; + WarpXParticleContainer::ParticleType& p = pp[ip]; + prod_p.id() = 9999; + prod_p.cpu() = cpuid; + prod_p.pos(0) = p.pos(0); + prod_p.pos(1) = p.pos(1); +#if (AMREX_SPACEDIM == 3) + prod_p.pos(20 = p.pos(2); +#endif + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + prod_pa[ia][i] = pa[ia][ip]; + } + } + } + ); + } // MFIter + } + } +} diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index a6cf9d81f..6f9fdc485 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -108,6 +108,8 @@ public: virtual void PostRestart () final {} void SplitParticles(int lev); + + virtual int doFieldIonization (const int lev) override; // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0183a44b9..8232b621f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -215,7 +215,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); particle_tile.push_back_real(particle_comps["zold"], z); @@ -591,7 +591,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) } if (do_field_ionization) { - pi[ip] = 2.; // species_ionization_level; + pi[ip] = species_ionization_level; } u.x *= PhysConst::c; @@ -2105,4 +2105,21 @@ void PhysicalParticleContainer::InitIonizationModule() adk_power.resize(ion_atomic_number); adk_prefactor.resize(ion_atomic_number); adk_exp_prefactor.resize(ion_atomic_number); + for (int i=0; i Date: Mon, 5 Aug 2019 17:48:59 -0700 Subject: get dt from WarpX instance, to initialize ionization prefactors --- Source/Particles/MultiParticleContainer.cpp | 45 ++++++++++++++------------ Source/Particles/PhysicalParticleContainer.cpp | 18 ++++++++--- 2 files changed, 38 insertions(+), 25 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 329a43ca7..a79670582 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -34,10 +34,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); - // For each species, get the ID of its product species. - // This is used for ionization and pair creation processes. - mapSpeciesProduct(); - // 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. @@ -144,6 +140,9 @@ MultiParticleContainer::InitData () pc->InitData(); } pc_tmp->InitData(); + // For each species, get the ID of its product species. + // This is used for ionization and pair creation processes. + mapSpeciesProduct(); } @@ -511,21 +510,23 @@ MultiParticleContainer::getSpeciesID(std::string product_str) return i_product; } +//void +//MultiParticleContainer::InitIonizationModule() + void MultiParticleContainer::doFieldIonization() { - /* - amrex::Gpu::SharedMemory grid_ids; - amrex::Gpu::SharedMemory tile_ids; - amrex::Gpu::SharedMemory is_ionized; - */ - for (auto& pc : allcontainers){ + if (!pc->do_field_ionization){ continue; } + + auto& prod_pc = allcontainers[pc->ionization_product]; + const Real * const AMREX_RESTRICT p_ionization_energies = pc->ionization_energies.dataPtr(); const Real * const AMREX_RESTRICT p_adk_prefactor = pc->adk_prefactor.dataPtr(); const Real * const AMREX_RESTRICT p_adk_exp_prefactor = pc->adk_exp_prefactor.dataPtr(); const Real * const AMREX_RESTRICT p_adk_power = pc->adk_power.dataPtr(); + for (int lev = 0; lev <= pc->finestLevel(); ++lev){ #ifdef _OPENMP @@ -537,6 +538,10 @@ MultiParticleContainer::doFieldIonization() if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); } + prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + if ( (prod_pc->NumRuntimeRealComps()>0) || (prod_pc->NumRuntimeIntComps()>0) ) { + prod_pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } } #endif MFItInfo info; @@ -575,7 +580,7 @@ MultiParticleContainer::doFieldIonization() ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { - Real random_draw = Random(); + Real random_draw = amrex::Random(); Real ga = std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + @@ -589,15 +594,14 @@ MultiParticleContainer::doFieldIonization() int ilev = (int) round(ilev_real[i]); // int ilev = static_cast(round(ilev_real[i])); Real p; - Real w_dtau; - if (E<1.e-100*(p_ionization_energies[0])){ - p = 0.; - } else { - w_dtau = 1./ ga * p_adk_prefactor[ilev] * - std::pow(E,p_adk_power[ilev]) * - std::exp( p_adk_exp_prefactor[ilev]/E ); - p = 1. - std::exp( - w_dtau ); - } + //if (E<1.e-20*(p_ionization_energies[0])){ + // p = 0.; + //} else { + Real w_dtau = 1./ ga * p_adk_prefactor[ilev] * + std::pow(E,p_adk_power[ilev]) * + std::exp( p_adk_exp_prefactor[ilev]/E ); + p = 1. - std::exp( - w_dtau ); + //} p_is_ionized[i] = 0; if (random_draw < p){ ilev_real[i] += 1.; @@ -617,7 +621,6 @@ MultiParticleContainer::doFieldIonization() break; } - auto& prod_pc = allcontainers[pc->ionization_product]; auto& prod_particle_tile = prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; const int np_old = prod_particle_tile.GetArrayOfStructs().size(); const int np_new = np_old + np_ionized; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 8232b621f..fa01f472d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -41,8 +41,8 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_field_ionization", do_field_ionization); // If do_field_ionization, read initialization data from input file and // read ionization energies from table. - if (do_field_ionization) - InitIonizationModule(); + //if (do_field_ionization) + // InitIonizationModule(); pp.query("plot_species", plot_species); int do_user_plot_vars; @@ -90,6 +90,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) void PhysicalParticleContainer::InitData() { + if (do_field_ionization) {InitIonizationModule();} AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute } @@ -2078,6 +2079,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, void PhysicalParticleContainer::InitIonizationModule() { + if (!do_field_ionization) return; ParmParse pp(species_name); pp.query("ionization_level", species_ionization_level); pp.get("ionization_product", ionization_product_name); @@ -2092,16 +2094,22 @@ void PhysicalParticleContainer::InitIonizationModule() int offset = ion_energy_offsets[ion_element_id]; for(int i=0; i Date: Tue, 6 Aug 2019 13:07:57 -0700 Subject: cleaning, and add capability for boosted frame runtime particle components --- Source/Particles/MultiParticleContainer.cpp | 318 +++++++++++++++---------- Source/Particles/PhysicalParticleContainer.H | 3 + Source/Particles/PhysicalParticleContainer.cpp | 82 +++++++ Source/Particles/WarpXParticleContainer.H | 8 + 4 files changed, 283 insertions(+), 128 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index a79670582..6c6a03231 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -513,155 +513,217 @@ MultiParticleContainer::getSpeciesID(std::string product_str) //void //MultiParticleContainer::InitIonizationModule() +namespace +{ + static void createIonizedParticles( + int lev, const MFIter& mfi, + std::unique_ptr< WarpXParticleContainer>& pc_source, + std::unique_ptr< WarpXParticleContainer>& pc_product, + const int * const p_is_ionized) + { + + 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(); + // --- source AoS particle data + WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); + // --- source SoA particle data + auto& soa_source = ptile_source.GetStructOfArrays(); + GpuArray attribs_source; + for (int ia = 0; ia < PIdx::nattribs; ++ia) { + attribs_source[ia] = soa_source.GetRealData(ia).data(); + } + // --- source runtime attribs + bool do_boosted_product = WarpX::do_boosted_frame_diagnostic + && pc_product->DoBoostedFrameDiags(); + bool do_boosted_source = WarpX::do_boosted_frame_diagnostic + && pc_source->DoBoostedFrameDiags(); + GpuArray runtime_attribs_source; + // + if (do_boosted_product && do_boosted_source) { + // If boosted frame diagnostics for source species, store them + std::map comps_source = pc_source->getParticleComps(); + runtime_attribs_source[0] = soa_source.GetRealData(comps_source[ "xold"]).data(); + runtime_attribs_source[1] = soa_source.GetRealData(comps_source[ "yold"]).data(); + runtime_attribs_source[2] = soa_source.GetRealData(comps_source[ "zold"]).data(); + runtime_attribs_source[3] = soa_source.GetRealData(comps_source["uxold"]).data(); + runtime_attribs_source[4] = soa_source.GetRealData(comps_source["uyold"]).data(); + runtime_attribs_source[5] = soa_source.GetRealData(comps_source["uzold"]).data(); + } else if (do_boosted_product && !do_boosted_source){ + // Otherwise, store current particle momenta. + // Positions are copied from AoS data. + runtime_attribs_source[3] = soa_source.GetRealData(PIdx::ux).data(); + runtime_attribs_source[4] = soa_source.GetRealData(PIdx::uy).data(); + runtime_attribs_source[5] = soa_source.GetRealData(PIdx::uz).data(); + } + + // Indices of product particle for each ionized source particle. + // i_product[i] is the location in product tile of product particle + // from source particle i. + amrex::Gpu::ManagedVector i_product; + i_product.resize(np_source); + // 0GetParticles(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 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 runtime_attribs_product; + if (do_boosted_product) { + std::map 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]){ + int ip = p_i_product[is]; + // is: index of ionized particle in source species + // ip: index of corresponding created 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. For position, we need a different + // handling depending on do_boosted_source. For momentum, + // runtime_attribs_source[3-5] contains appropriate data. + if (do_boosted_product) { + if (do_boosted_source) { + runtime_attribs_product[0][ip] = runtime_attribs_source[0][ip]; + runtime_attribs_product[1][ip] = runtime_attribs_source[1][ip]; + runtime_attribs_product[2][ip] = runtime_attribs_source[2][ip]; + } else { + 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_attribs_source[3][ip]; + runtime_attribs_product[4][ip] = runtime_attribs_source[4][ip]; + runtime_attribs_product[5][ip] = runtime_attribs_source[5][ip]; + } + } + } + ); + } +} + void MultiParticleContainer::doFieldIonization() { - for (auto& pc : allcontainers){ - - if (!pc->do_field_ionization){ continue; } - auto& prod_pc = allcontainers[pc->ionization_product]; + // Loop over all species. + // Ionized particles in pc_source create particles in pc_product + for (auto& pc_source : allcontainers){ + + // Skip if not ionizable + if (!pc_source->do_field_ionization){ continue; } - const Real * const AMREX_RESTRICT p_ionization_energies = pc->ionization_energies.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_prefactor = pc->adk_prefactor.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_exp_prefactor = pc->adk_exp_prefactor.dataPtr(); - const Real * const AMREX_RESTRICT p_adk_power = pc->adk_power.dataPtr(); + // Get product species + auto& pc_product = allcontainers[pc_source->ionization_product]; - for (int lev = 0; lev <= pc->finestLevel(); ++lev){ + for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ #ifdef _OPENMP - // First touch all tiles in the map in serial - for (MFIter mfi = pc->MakeMFIter(lev); mfi.isValid(); ++mfi) { + // Touch all tiles of source species in serial if additional arguments + for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - if ( (pc->NumRuntimeRealComps()>0) || (pc->NumRuntimeIntComps()>0) ) { - pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); - } - prod_pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - if ( (prod_pc->NumRuntimeRealComps()>0) || (prod_pc->NumRuntimeIntComps()>0) ) { - prod_pc->DefineAndReturnParticleTile(lev, grid_id, tile_id); + pc_source->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + if ( (pc_source->NumRuntimeRealComps()>0) || (pc_source->NumRuntimeIntComps()>0) ) { + pc_source->DefineAndReturnParticleTile(lev, grid_id, tile_id); } } #endif + // Touch all tiles of product species in serial + for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + pc_product->DefineAndReturnParticleTile(lev, grid_id, tile_id); + } + + // Enable tiling MFItInfo info; - if (pc->do_tiling && Gpu::notInLaunchRegion()) { - info.EnableTiling(pc->tile_size); + if (pc_source->do_tiling && Gpu::notInLaunchRegion()) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + pc_product->do_tiling, + "For ionization, either all or none of the " + "particle species must use tiling."); + info.EnableTiling(pc_source->tile_size); } + #ifdef _OPENMP info.SetDynamic(true); #pragma omp parallel #endif - for (MFIter mfi = pc->MakeMFIter(lev, info); mfi.isValid(); ++mfi) + // Loop over all grids (if not tiling) or grids and tiles (if tiling) + for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi) { - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - auto& particle_tile = pc->GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - auto& soa = particle_tile.GetStructOfArrays(); - const int np = particle_tile.GetArrayOfStructs().size(); - if (np == 0) break; + // Ionization mask: one element per source particles. + // 0 if not ionized, 1 if ionized. amrex::Gpu::ManagedVector is_ionized; - // const int np_lev = pc->NumberOfParticlesAtLevel(lev); - is_ionized.resize(np); - int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); - const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); - const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); - const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); - const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); - const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); - const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); - const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); - const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); - Real * const AMREX_RESTRICT ilev_real = soa.GetRealData(pc->particle_comps["ionization_level"]).data(); - Real c = PhysConst::c; - Real c2_inv = 1./c/c; - - ParallelFor( - np, - [=] AMREX_GPU_DEVICE (long i) { - Real random_draw = amrex::Random(); - Real ga = std::sqrt(1. + - (ux[i]*ux[i] + - uy[i]*uy[i] + - uz[i]*uz[i]) * c2_inv); - Real E = std::sqrt( - - ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * c2_inv - + ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) * ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) - + ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) * ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) - + ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) * ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) - ); - int ilev = (int) round(ilev_real[i]); - // int ilev = static_cast(round(ilev_real[i])); - Real p; - //if (E<1.e-20*(p_ionization_energies[0])){ - // p = 0.; - //} else { - Real w_dtau = 1./ ga * p_adk_prefactor[ilev] * - std::pow(E,p_adk_power[ilev]) * - std::exp( p_adk_exp_prefactor[ilev]/E ); - p = 1. - std::exp( - w_dtau ); - //} - p_is_ionized[i] = 0; - if (random_draw < p){ - ilev_real[i] += 1.; - p_is_ionized[i] = 1; - } - } - ); // ParallelFor - - amrex::Gpu::ManagedVector is_ionized_cumsum_vector; - is_ionized_cumsum_vector.resize(np); - int np_ionized = p_is_ionized[0]; - for(int i=1; iGetParticles(lev)[std::make_pair(grid_id,tile_id)]; - const int np_old = prod_particle_tile.GetArrayOfStructs().size(); - const int np_new = np_old + np_ionized; - prod_particle_tile.resize(np_new); - auto& prod_soa = prod_particle_tile.GetStructOfArrays(); - - GpuArray prod_pa; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - prod_pa[ia] = prod_soa.GetRealData(ia).data() + np_old; - } - WarpXParticleContainer::ParticleType* prod_pp = prod_particle_tile.GetArrayOfStructs()().data() + np_old; - WarpXParticleContainer::ParticleType* pp = particle_tile.GetArrayOfStructs()().data(); - - GpuArray pa; - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - pa[ia] = soa.GetRealData(ia).data(); - } + pc_source->buildIonizationMask(mfi, lev, is_ionized); + const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - const int cpuid = ParallelDescriptor::MyProc(); - int* AMREX_RESTRICT is_ionized_cumsum = is_ionized_cumsum_vector.dataPtr(); - amrex::For( - np, [=] AMREX_GPU_DEVICE (int ip) noexcept - { - if(is_ionized[ip]){ - int i = is_ionized_cumsum[ip]; - WarpXParticleContainer::ParticleType& prod_p = prod_pp[i]; - WarpXParticleContainer::ParticleType& p = pp[ip]; - prod_p.id() = 9999; - prod_p.cpu() = cpuid; - prod_p.pos(0) = p.pos(0); - prod_p.pos(1) = p.pos(1); -#if (AMREX_SPACEDIM == 3) - prod_p.pos(2) = p.pos(2); -#endif - for (int ia = 0; ia < PIdx::nattribs; ++ia) { - prod_pa[ia][i] = pa[ia][ip]; - } - } - } - ); + createIonizedParticles(lev, mfi, pc_source, pc_product, p_is_ionized); } // MFIter - } - } + } // lev + } // pc_source } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 6f9fdc485..8dab1a2f7 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -111,6 +111,9 @@ public: virtual int doFieldIonization (const int lev) override; + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedVector& ionization_mask) override; + // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index fa01f472d..7870e683c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2133,3 +2133,85 @@ PhysicalParticleContainer::doFieldIonization(const int lev) Print()<<"in PhysicalParticleContainer::doFieldIonization\n"; return 0; } + +/* \brief create mask of ionized particles (1 if ionized, 0 otherwise) + * + * \param mfi: tile or grid + * \param lev: MR level + * \param ionization_mask: Array with as many elements as particles in mfi. + * This function initialized the array, and set each element to 1 or 0 + * depending on whether the particle is ionized or not. + */ +void +PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedVector& ionization_mask) +{ + // Get pointers to ionization data from pc_source + const Real * const AMREX_RESTRICT p_ionization_energies = ionization_energies.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_prefactor = adk_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_exp_prefactor = adk_exp_prefactor.dataPtr(); + const Real * const AMREX_RESTRICT p_adk_power = adk_power.dataPtr(); + + // Current tile info + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + // Get GPU-friendly arrays of particle data + auto& ptile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + // Only need attribs (i.e., SoA data) + auto& soa = ptile.GetStructOfArrays(); + const int np = ptile.GetArrayOfStructs().size(); + + // If no particle, nothing to do. + if (np == 0) return; + // Otherwise, resize ionization_mask, and get poiters to attribs arrays. + ionization_mask.resize(np); + int * const AMREX_RESTRICT p_ionization_mask = ionization_mask.data(); + const Real * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + const Real * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data(); + const Real * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data(); + const Real * const AMREX_RESTRICT ex = soa.GetRealData(PIdx::Ex).data(); + const Real * const AMREX_RESTRICT ey = soa.GetRealData(PIdx::Ey).data(); + const Real * const AMREX_RESTRICT ez = soa.GetRealData(PIdx::Ez).data(); + const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); + const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); + const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); + Real* ilev_real = soa.GetRealData(particle_comps["ionization_level"]).data(); + + Real c = PhysConst::c; + Real c2_inv = 1./c/c; + + // Loop over all particles in grid/tile. If ionized, set mask to 1 + // and increment ionization level. + ParallelFor( + np, + [=] AMREX_GPU_DEVICE (long i) { + Real random_draw = amrex::Random(); + // Compute electric field amplitude in the particle's frame of + // reference (particularly important when in boosted frame). + Real ga = std::sqrt(1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i]) * c2_inv); + Real E = std::sqrt( + - ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * ( ux[i]*ex[i] + uy[i]*ey[i] + uz[i]*ez[i] ) * c2_inv + + ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) * ( ga *ex[i] + uy[i]*bz[i] - uz[i]*by[i] ) + + ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) * ( ga *ey[i] + uz[i]*bx[i] - ux[i]*bz[i] ) + + ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) * ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) + ); + // Get index of ionization_level + int ilev = (int) round(ilev_real[i]); + // Compute probability of ionization p + Real p; + Real w_dtau = 1./ ga * p_adk_prefactor[ilev] * + std::pow(E,p_adk_power[ilev]) * + std::exp( p_adk_exp_prefactor[ilev]/E ); + p = 1. - std::exp( - w_dtau ); + p_ionization_mask[i] = 0; + + if (random_draw < p){ + // increment particle's ionization level + ilev_real[i] += 1.; + // update mask + p_ionization_mask[i] = 1; + } + } + ); +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 1c1b3fd22..017a857a7 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -250,6 +250,8 @@ public: static int NextID () { return ParticleType::NextID(); } + void setNextID(int next_id) { ParticleType::NextID(next_id); } + bool do_splitting = false; // split along axes (0) or diagonals (1) @@ -277,6 +279,12 @@ public: virtual int doFieldIonization ( const int lev) { return 0; }; + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedVector& ionization_mask) + {}; + + std::map getParticleComps() { return particle_comps;} + protected: std::map particle_comps; -- cgit v1.2.3 From dab4a0081eac6929db26d16b79ef57c5882df0b4 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 7 Aug 2019 15:06:49 -0700 Subject: use c++ deposition by default. includes old attribs --- Source/Particles/MultiParticleContainer.cpp | 28 +++----------- Source/Particles/PhysicalParticleContainer.cpp | 53 ++++++++++++-------------- Source/WarpX.cpp | 2 +- 3 files changed, 31 insertions(+), 52 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 84cd107f8..36b41e0aa 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -516,12 +516,11 @@ MultiParticleContainer::getSpeciesID(std::string product_str) namespace { static void createIonizedParticles( - int lev, const MFIter& mfi, + int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, const int * const p_is_ionized) { - // Print()<<"in createIonizedParticles\n"; const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); @@ -530,7 +529,6 @@ namespace const int np_source = ptile_source.GetArrayOfStructs().size(); if (np_source == 0) return; // --- source AoS particle data - //Print()<<"la 1\n"; WarpXParticleContainer::ParticleType* particles_source = ptile_source.GetArrayOfStructs()().data(); // --- source SoA particle data auto& soa_source = ptile_source.GetStructOfArrays(); @@ -544,10 +542,10 @@ namespace bool do_boosted_source = WarpX::do_boosted_frame_diagnostic && pc_source->DoBoostedFrameDiags(); GpuArray runtime_attribs_source; - // - //Print()<<"la 2\n"; + // Prepare arrays for boosted frame diagnostics. + // If do_boosted_product, need different treatment + // depending on do_boosted_source if (do_boosted_product && do_boosted_source) { - //Print()<<"do_boosted_product && do_boosted_source\n"; // If boosted frame diagnostics for source species, store them std::map comps_source = pc_source->getParticleComps(); runtime_attribs_source[0] = soa_source.GetRealData(comps_source[ "xold"]).data(); @@ -557,7 +555,6 @@ namespace runtime_attribs_source[4] = soa_source.GetRealData(comps_source["uyold"]).data(); runtime_attribs_source[5] = soa_source.GetRealData(comps_source["uzold"]).data(); } else if (do_boosted_product && !do_boosted_source){ - //Print()<<"do_boosted_product && !do_boosted_source\n"; // Otherwise, store current particle momenta. // Positions are copied from AoS data. runtime_attribs_source[3] = soa_source.GetRealData(PIdx::ux).data(); @@ -568,27 +565,20 @@ namespace // Indices of product particle for each ionized source particle. // i_product[i] is the location in product tile of product particle // from source particle i. - //Print()<<"la 3\n"; amrex::Gpu::ManagedVector i_product; i_product.resize(np_source); // 0GetParticles(lev)[std::make_pair(grid_id,tile_id)]; // old and new (i.e., including ionized particles) number of particles @@ -601,18 +591,15 @@ namespace // First element is the first newly-created product particle WarpXParticleContainer::ParticleType* particles_product = ptile_product.GetArrayOfStructs()().data() + np_product_old; // --- product SoA particle data - //Print()<<"la 8\n"; auto& soa_product = ptile_product.GetStructOfArrays(); GpuArray 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; } - //Print()<<"la 9\n"; // --- product runtime attribs GpuArray runtime_attribs_product; if (do_boosted_product) { - //Print()<<"do_boosted_product 1\n"; std::map 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; @@ -640,7 +627,7 @@ namespace if(p_is_ionized[is]){ int ip = p_i_product[is]; // is: index of ionized particle in source species - // ip: index of corresponding created particle in product 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 @@ -660,23 +647,18 @@ namespace // handling depending on do_boosted_source. For momentum, // runtime_attribs_source[3-5] contains appropriate data. if (do_boosted_product) { - //Print()<<"do_boosted_product\n"; if (do_boosted_source) { - //Print()<<"do_boosted_product do_boosted_source\n"; runtime_attribs_product[0][ip] = runtime_attribs_source[0][ip]; runtime_attribs_product[1][ip] = runtime_attribs_source[1][ip]; runtime_attribs_product[2][ip] = runtime_attribs_source[2][ip]; } else { - //Print()<<"do_boosted_product NO do_boosted_source\n"; 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); } - //Print()<<"momenta\n"; runtime_attribs_product[3][ip] = runtime_attribs_source[3][ip]; runtime_attribs_product[4][ip] = runtime_attribs_source[4][ip]; runtime_attribs_product[5][ip] = runtime_attribs_source[5][ip]; - //Print()<<"done\n"; } } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7870e683c..75e975131 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2094,14 +2094,11 @@ void PhysicalParticleContainer::InitIonizationModule() int offset = ion_energy_offsets[ion_element_id]; for(int i=0; i WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::use_picsar_deposition = 1; +long WarpX::use_picsar_deposition = 0; long WarpX::current_deposition_algo; long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; -- cgit v1.2.3 From 4be45ed13bcd1dfab7455e09b2a9b96cd2064fca Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 7 Aug 2019 15:18:36 -0700 Subject: cleaning, remote outdated functions --- Source/Particles/MultiParticleContainer.cpp | 9 +-------- Source/Particles/PhysicalParticleContainer.H | 2 -- Source/Particles/PhysicalParticleContainer.cpp | 13 ------------- Source/Particles/WarpXParticleContainer.H | 3 --- 4 files changed, 1 insertion(+), 26 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 136d23f6c..48f685e8f 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -646,8 +646,6 @@ void MultiParticleContainer::doFieldIonization() { - Print()<<"in MultiParticleContainer::doFieldIonization\n"; - // Loop over all species. // Ionized particles in pc_source create particles in pc_product for (auto& pc_source : allcontainers){ @@ -698,16 +696,11 @@ MultiParticleContainer::doFieldIonization() { // Ionization mask: one element per source particles. // 0 if not ionized, 1 if ionized. - // Print()<<"here 1\n"; amrex::Gpu::ManagedVector is_ionized; - // Print()<<"here 2\n"; pc_source->buildIonizationMask(mfi, lev, is_ionized); - // Print()<<"here 3\n"; const int * const AMREX_RESTRICT p_is_ionized = is_ionized.dataPtr(); - // Print()<<"here 4\n"; createIonizedParticles(lev, mfi, pc_source, pc_product, p_is_ionized); - // Print()<<"here 5\n"; - } // MFIter + } } // lev } // pc_source } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 8dab1a2f7..c10dcc61b 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -108,8 +108,6 @@ public: virtual void PostRestart () final {} void SplitParticles(int lev); - - virtual int doFieldIonization (const int lev) override; virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedVector& ionization_mask) override; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 75e975131..02602b844 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -80,12 +80,6 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) : WarpXParticleContainer(amr_core, 0) { plasma_injector.reset(new PlasmaInjector()); -/* - if (do_field_ionization){ - Print()<<"AddRealComp\n"; - AddRealComp("ionization_level"); - } -*/ } void PhysicalParticleContainer::InitData() @@ -2121,13 +2115,6 @@ void PhysicalParticleContainer::InitIonizationModule() } } -int -PhysicalParticleContainer::doFieldIonization(const int lev) -{ - Print()<<"in PhysicalParticleContainer::doFieldIonization\n"; - return 0; -} - /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) * * \param mfi: tile or grid diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 017a857a7..82190b307 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -276,9 +276,6 @@ public: virtual void setIonizationProduct (int i_product) {}; - virtual int doFieldIonization ( - const int lev) { return 0; }; - virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedVector& ionization_mask) {}; -- cgit v1.2.3 From 30fa1002b9d106acbbbbfd0f01c7e80f47c72747 Mon Sep 17 00:00:00 2001 From: Maxence Thevenet Date: Wed, 7 Aug 2019 18:57:59 -0400 Subject: add profile data, remove typo in example --- Examples/Modules/ionization/inputs.2d | 1 - Source/Particles/MultiParticleContainer.cpp | 4 +++- Source/Particles/PhysicalParticleContainer.cpp | 1 + 3 files changed, 4 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Modules/ionization/inputs.2d b/Examples/Modules/ionization/inputs.2d index 1f2b77a7e..ddf30763a 100644 --- a/Examples/Modules/ionization/inputs.2d +++ b/Examples/Modules/ionization/inputs.2d @@ -4,7 +4,6 @@ max_step = 500 warpx.verbose = 1 warpx.cfl = 1 -warpx.plot_rho = 1 warpx.do_pml = 1 warpx.use_filter = 1 diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 48f685e8f..dafe8fe93 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -521,6 +521,8 @@ namespace std::unique_ptr< WarpXParticleContainer>& pc_product, const int * const p_is_ionized) { + BL_PROFILE("createIonizedParticles"); + const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); @@ -645,7 +647,7 @@ namespace void MultiParticleContainer::doFieldIonization() { - + BL_PROFILE("MPC::doFieldIonization"); // Loop over all species. // Ionized particles in pc_source create particles in pc_product for (auto& pc_source : allcontainers){ diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 02602b844..22b7df88f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2127,6 +2127,7 @@ void PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedVector& ionization_mask) { + BL_PROFILE("PPC::buildIonizationMask"); // Get pointers to ionization data from pc_source const Real * const AMREX_RESTRICT p_ionization_energies = ionization_energies.dataPtr(); const Real * const AMREX_RESTRICT p_adk_prefactor = adk_prefactor.dataPtr(); -- cgit v1.2.3 From b2754b63e9a25c550cc849f84078503bc94a417f Mon Sep 17 00:00:00 2001 From: Maxence Thevenet Date: Wed, 7 Aug 2019 21:17:42 -0400 Subject: exclusive_scan requires ManagedDeviceVector --- Source/Particles/MultiParticleContainer.cpp | 6 +++--- Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.H | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 9eb978d05..6d977e16c 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -519,7 +519,7 @@ namespace int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, - amrex::Gpu::ManagedVector& is_ionized) + amrex::Gpu::ManagedDeviceVector& is_ionized) { BL_PROFILE("createIonizedParticles"); @@ -550,7 +550,7 @@ namespace // 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::ManagedVector i_product; + amrex::Gpu::ManagedDeviceVector i_product; i_product.resize(np_source); // 0 is_ionized; + amrex::Gpu::ManagedDeviceVector is_ionized; pc_source->buildIonizationMask(mfi, lev, is_ionized); // Create particles in pc_product createIonizedParticles(lev, mfi, pc_source, pc_product, is_ionized); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index c10dcc61b..44c5030fa 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -110,7 +110,7 @@ public: void SplitParticles(int lev); virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedVector& ionization_mask) override; + amrex::Gpu::ManagedDeviceVector& ionization_mask) override; // Inject particles in Box 'part_box' virtual void AddParticles (int lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 22b7df88f..f1544d90c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2125,7 +2125,7 @@ void PhysicalParticleContainer::InitIonizationModule() */ void PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedVector& ionization_mask) + amrex::Gpu::ManagedDeviceVector& ionization_mask) { BL_PROFILE("PPC::buildIonizationMask"); // Get pointers to ionization data from pc_source diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 82190b307..334eb3282 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -277,7 +277,7 @@ public: virtual void setIonizationProduct (int i_product) {}; virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedVector& ionization_mask) + amrex::Gpu::ManagedDeviceVector& ionization_mask) {}; std::map getParticleComps() { return particle_comps;} -- cgit v1.2.3 From c7a56a07ee33f5d1cce29719a946036e36b5897c Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 11:00:08 -0700 Subject: change ionization level from Real to Int --- Source/Diagnostics/ParticleIO.cpp | 3 ++- Source/Laser/LaserParticleContainer.cpp | 2 +- Source/Particles/Deposition/CurrentDeposition.H | 4 ++-- Source/Particles/PhysicalParticleContainer.cpp | 29 +++++++++++++------------ Source/Particles/WarpXParticleContainer.H | 9 ++++++-- Source/Particles/WarpXParticleContainer.cpp | 2 +- 6 files changed, 28 insertions(+), 21 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 2b70c9379..618ecc6c5 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -114,7 +114,8 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const } if(pc->do_field_ionization){ - real_names.push_back("ionization_level"); + int_names.push_back("ionization_level"); + int_flags.resize(1, 1); } // Convert momentum to SI diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 37e8f8bf2..e2aebdd5e 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -504,7 +504,7 @@ LaserParticleContainer::Evolve (int lev, // Current Deposition // // Deposit inside domains - Real* ion_lev = nullptr; + int* ion_lev = nullptr; DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, lev, lev, dt); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index deff65424..86cf89cd3 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -26,7 +26,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, - const amrex::Real * const ion_lev, + const int * const ion_lev, const amrex::Array4& jx_arr, const amrex::Array4& jy_arr, const amrex::Array4& jz_arr, @@ -184,7 +184,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real * const uxp, const amrex::Real * const uyp, const amrex::Real * const uzp, - const amrex::Real * ion_lev, + const int * ion_lev, const amrex::Array4& Jx_arr, const amrex::Array4& Jy_arr, const amrex::Array4& Jz_arr, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7fdb783b7..98d9e5d18 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -452,9 +452,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; } - Real* pi; + int* pi; if (do_field_ionization) { - pi = soa.GetRealData(particle_comps[ "ionization_level"]).data() + old_size; + pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } const GpuArray overlap_corner @@ -1293,9 +1293,9 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt); } } else { - Real* AMREX_RESTRICT ion_lev; + int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ - ion_lev = pti.GetAttribs(particle_comps["ionization_level"]).dataPtr(); + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } else { ion_lev = nullptr; } @@ -1523,9 +1523,9 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, copy_attribs(pti, x, y, z); } - Real* AMREX_RESTRICT ion_lev = nullptr; + int* AMREX_RESTRICT ion_lev = nullptr; if (do_field_ionization){ - ion_lev = pti.GetAttribs(particle_comps["ionization_level"]).dataPtr(); + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } // Loop over the particles and update their momentum @@ -1972,7 +1972,7 @@ void PhysicalParticleContainer::InitIonizationModule() pp.get("ionization_product", ionization_product_name); pp.get("physical_element", physical_element); // Add Real component for ionization level - AddRealComp("ionization_level"); + AddIntComp("ionization_level"); plot_flags.resize(PIdx::nattribs + 1, 1); // Get atomic number and ionization energies from file int ion_element_id = ion_map_ids[physical_element]; @@ -2051,7 +2051,8 @@ PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const i const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); - Real* ilev_real = soa.GetRealData(particle_comps["ionization_level"]).data(); + // Real* ilev_real = soa.GetIntData(particle_icomps["ionization_level"]).data(); + int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); Real c = PhysConst::c; Real c2_inv = 1./c/c; @@ -2063,9 +2064,9 @@ PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const i np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level - int ilev = (int) round(ilev_real[i]); + // int ilev = (int) round(ilev_real[i]); p_ionization_mask[i] = 0; - if ( ilev particle_comps; + std::map particle_icomps; int species_id; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 2e418aae0..a7c30764d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -428,7 +428,7 @@ void WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, RealVector& wp, RealVector& uxp, RealVector& uyp, RealVector& uzp, - Real* ion_lev, + int* ion_lev, MultiFab* jx, MultiFab* jy, MultiFab* jz, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, -- cgit v1.2.3 From 6123a53ecc4a75b617cd63b5ebde50ada940a78e Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 14:21:07 -0700 Subject: minor cleaning, remove commented lines --- Source/Particles/PhysicalParticleContainer.cpp | 2 -- 1 file changed, 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 98d9e5d18..a9b57aef6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2051,7 +2051,6 @@ PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const i const Real * const AMREX_RESTRICT bx = soa.GetRealData(PIdx::Bx).data(); const Real * const AMREX_RESTRICT by = soa.GetRealData(PIdx::By).data(); const Real * const AMREX_RESTRICT bz = soa.GetRealData(PIdx::Bz).data(); - // Real* ilev_real = soa.GetIntData(particle_icomps["ionization_level"]).data(); int* ion_lev = soa.GetIntData(particle_icomps["ionization_level"]).data(); Real c = PhysConst::c; @@ -2064,7 +2063,6 @@ PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const i np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level - // int ilev = (int) round(ilev_real[i]); p_ionization_mask[i] = 0; if ( ion_lev[i] Date: Thu, 8 Aug 2019 14:34:40 -0700 Subject: use Weiqun's convention: space before brackets for function declaration and definition --- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 10 +++++----- Source/Particles/PhysicalParticleContainer.H | 4 +--- Source/Particles/PhysicalParticleContainer.cpp | 12 +++++------- Source/Particles/WarpXParticleContainer.H | 2 +- 5 files changed, 13 insertions(+), 17 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index fd9a33e86..aaed96423 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -127,7 +127,7 @@ public: /// std::unique_ptr GetChargeDensity(int lev, bool local = false); - void doFieldIonization(); + void doFieldIonization (); void Checkpoint (const std::string& dir) const; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 24f9f2f11..04f052086 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -453,7 +453,7 @@ MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const } int -MultiParticleContainer::doContinuousInjection() const +MultiParticleContainer::doContinuousInjection () const { int warpx_do_continuous_injection = 0; for (int i=0; i& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, @@ -647,7 +647,7 @@ namespace } void -MultiParticleContainer::doFieldIonization() +MultiParticleContainer::doFieldIonization () { BL_PROFILE("MPC::doFieldIonization"); // Loop over all species. diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 44c5030fa..c7494fbdf 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -22,9 +22,7 @@ public: virtual void InitData () override; - void InitIonizationLevel (); - - void InitIonizationModule(); + void InitIonizationModule (); #ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector, 3> >& E, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a9b57aef6..b20ea8090 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -39,10 +39,6 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("do_field_ionization", do_field_ionization); - // If do_field_ionization, read initialization data from input file and - // read ionization energies from table. - //if (do_field_ionization) - // InitIonizationModule(); pp.query("plot_species", plot_species); int do_user_plot_vars; @@ -84,6 +80,8 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) void PhysicalParticleContainer::InitData() { + // Init ionization module here instead of in the PhysicalParticleContainer + // constructor to have access to dt if (do_field_ionization) {InitIonizationModule();} AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute @@ -1964,7 +1962,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, } -void PhysicalParticleContainer::InitIonizationModule() +void PhysicalParticleContainer::InitIonizationModule () { if (!do_field_ionization) return; ParmParse pp(species_name); @@ -2017,8 +2015,8 @@ void PhysicalParticleContainer::InitIonizationModule() * depending on whether the particle is ionized or not. */ void -PhysicalParticleContainer::buildIonizationMask(const amrex::MFIter& mfi, const int lev, - amrex::Gpu::ManagedDeviceVector& ionization_mask) +PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const int lev, + amrex::Gpu::ManagedDeviceVector& ionization_mask) { BL_PROFILE("PPC::buildIonizationMask"); // Get pointers to ionization data from pc_source diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index dd5a150c4..efa299c43 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -284,7 +284,7 @@ public: amrex::Gpu::ManagedDeviceVector& ionization_mask) {}; - std::map getParticleComps() { return particle_comps;} + std::map getParticleComps () { return particle_comps;} protected: -- cgit v1.2.3 From c87ee4b5e5d62f8e90a719a42ff164f99983a62c Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 14:44:06 -0700 Subject: final minor cleaning --- Source/Particles/MultiParticleContainer.cpp | 6 +++--- Source/Particles/PhysicalParticleContainer.cpp | 6 +++--- Source/Utils/WarpXConst.H | 1 - 3 files changed, 6 insertions(+), 7 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 04f052086..e39930216 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -511,8 +511,8 @@ MultiParticleContainer::getSpeciesID (std::string product_str) namespace { - // For particle i in mfi, if is_ionized[i]=1, copy particle from - // particle i container pc_source into pc_product + // For particle i in mfi, if is_ionized[i]=1, copy particle + // particle i from container pc_source into pc_product static void createIonizedParticles ( int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, @@ -663,7 +663,7 @@ MultiParticleContainer::doFieldIonization () for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ #ifdef _OPENMP - // Touch all tiles of source species in serial if additional arguments + // Touch all tiles of source species in serial if runtime attribs for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b20ea8090..466446157 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -81,7 +81,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) void PhysicalParticleContainer::InitData() { // Init ionization module here instead of in the PhysicalParticleContainer - // constructor to have access to dt + // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute @@ -208,10 +208,10 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); particle_tile.push_back_real(particle_comps["zold"], z); - + particle_tile.push_back_real(particle_comps["uxold"], u[0]); particle_tile.push_back_real(particle_comps["uyold"], u[1]); particle_tile.push_back_real(particle_comps["uzold"], u[2]); diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index bfb67b3d5..6901d4629 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -21,7 +21,6 @@ namespace PhysConst static constexpr amrex::Real hbar = 1.05457180010e-34; static constexpr amrex::Real alpha = mu0/(4*MathConst::pi)*q_e*q_e*c/hbar; static constexpr amrex::Real r_e = 1./(4*MathConst::pi*ep0) * q_e*q_e/(m_e*c*c); - // static constexpr amrex::Real alpha = 0.0072973525693; } #endif -- cgit v1.2.3 From 8b50555b58cb88a61f806193d34aead9ffbf1612 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 19:17:00 -0700 Subject: hidden typo '-.- --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 466446157..2c6deefef 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2001,7 +2001,7 @@ void PhysicalParticleContainer::InitIonizationModule () adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) - * std::pow(std::pow(2*(Uion/UH),3./2)*Ea,2*n_eff - 1); + * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } -- cgit v1.2.3 From 1b9682b5d2705369a8a3c7da6ea1a00b9eb824d2 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 12 Aug 2019 09:52:36 -0700 Subject: typos in doc, and rename ionization_product --- Docs/source/running_cpp/parameters.rst | 8 ++++---- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 6b567d5b1..6d7000cd8 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -340,16 +340,16 @@ Particle initialization Do field ionization for this species (using the ADK theory). * ``.physical_element`` (`string`) - Only read of `do_field_ionization = 1`. Symbol of chemical element for + Only read if `do_field_ionization = 1`. Symbol of chemical element for this species. Example: for Helium, use ``physical_element = He``. -* ``.ionization_product`` (`string`) - Only read of `do_field_ionization = 1`. Name of species in which ionized +* ``.ionization_product_species`` (`string`) + Only read if `do_field_ionization = 1`. Name of species in which ionized electrons are stored. This species must be created as a regular species in the input file (in particular, it must be in `particles.species_names`). * ``.ionization_level`` (`int`) optional (default `0`) - Only read of `do_field_ionization = 1`. Initial ionization level of the + Only read if `do_field_ionization = 1`. Initial ionization level of the species (must be smaller than the atomic number of chemical element given in `physical_element`). diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 2c6deefef..69d7fa686 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1967,7 +1967,7 @@ void PhysicalParticleContainer::InitIonizationModule () if (!do_field_ionization) return; ParmParse pp(species_name); pp.query("ionization_level", species_ionization_level); - pp.get("ionization_product", ionization_product_name); + pp.get("ionization_product_species", ionization_product_name); pp.get("physical_element", physical_element); // Add Real component for ionization level AddIntComp("ionization_level"); -- cgit v1.2.3 From 1cc426d477fd160c0b59bf777a8be715527ccfb8 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 12 Aug 2019 09:56:33 -0700 Subject: use ionization_initial_level --- Docs/source/running_cpp/parameters.rst | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- Source/Particles/WarpXParticleContainer.H | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 6d7000cd8..59628366b 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -348,7 +348,7 @@ Particle initialization electrons are stored. This species must be created as a regular species in the input file (in particular, it must be in `particles.species_names`). -* ``.ionization_level`` (`int`) optional (default `0`) +* ``.ionization_initial_level`` (`int`) optional (default `0`) Only read if `do_field_ionization = 1`. Initial ionization level of the species (must be smaller than the atomic number of chemical element given in `physical_element`). diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 69d7fa686..a92f36a3c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -584,7 +584,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) } if (do_field_ionization) { - pi[ip] = species_ionization_level; + pi[ip] = ionization_initial_level; } u.x *= PhysConst::c; @@ -1966,7 +1966,7 @@ void PhysicalParticleContainer::InitIonizationModule () { if (!do_field_ionization) return; ParmParse pp(species_name); - pp.query("ionization_level", species_ionization_level); + pp.query("ionization_level", ionization_initial_level); pp.get("ionization_product_species", ionization_product_name); pp.get("physical_element", physical_element); // Add Real component for ionization level diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 2a625b4e1..b89406b3e 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -308,7 +308,7 @@ protected: int ionization_product; std::string ionization_product_name; int ion_atomic_number; - int species_ionization_level = 0; + int ionization_initial_level = 0; amrex::Gpu::ManagedVector ionization_energies; amrex::Gpu::ManagedVector adk_power; amrex::Gpu::ManagedVector adk_prefactor; -- cgit v1.2.3 From fab45a507406dfd764dd793b0b4d612e86297fbb Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 12 Aug 2019 10:57:02 -0700 Subject: typo --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a92f36a3c..a221be01c 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1966,7 +1966,7 @@ void PhysicalParticleContainer::InitIonizationModule () { if (!do_field_ionization) return; ParmParse pp(species_name); - pp.query("ionization_level", ionization_initial_level); + pp.query("ionization_initial_level", ionization_initial_level); pp.get("ionization_product_species", ionization_product_name); pp.get("physical_element", physical_element); // Add Real component for ionization level -- cgit v1.2.3 From 25690770d463071483bc78c99f29bbe056ec9355 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 19 Aug 2019 11:04:28 -0700 Subject: must check the overlap box is not empty --- Source/Particles/PhysicalParticleContainer.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d10390204..3d321c495 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -390,14 +390,16 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) // then how many new particles will be injected is not that simple // We have to shift fine_injection_box because overlap_box has been shifted. Box fine_overlap_box = overlap_box & amrex::shift(fine_injection_box,shifted); - max_new_particles += fine_overlap_box.numPts() * num_ppc - * (AMREX_D_TERM(rrfac,*rrfac,*rrfac)-1); - for (int icell = 0, ncells = overlap_box.numPts(); icell < ncells; ++icell) { - IntVect iv = overlap_box.atOffset(icell); - int r = (fine_overlap_box.contains(iv)) ? AMREX_D_TERM(rrfac,*rrfac,*rrfac) : 1; - for (int ipart = 0; ipart < r; ++ipart) { - cellid_v.push_back(icell); - cellid_v.push_back(ipart); + if (fine_overlap_box.ok()) { + max_new_particles += fine_overlap_box.numPts() * num_ppc + * (AMREX_D_TERM(rrfac,*rrfac,*rrfac)-1); + for (int icell = 0, ncells = overlap_box.numPts(); icell < ncells; ++icell) { + IntVect iv = overlap_box.atOffset(icell); + int r = (fine_overlap_box.contains(iv)) ? AMREX_D_TERM(rrfac,*rrfac,*rrfac) : 1; + for (int ipart = 0; ipart < r; ++ipart) { + cellid_v.push_back(icell); + cellid_v.push_back(ipart); + } } } } -- cgit v1.2.3 From 83f9298a244c57749e4c4a1a159c779466996278 Mon Sep 17 00:00:00 2001 From: Maxence Thevenet Date: Tue, 20 Aug 2019 14:25:51 -0400 Subject: copy member variables to tmp local variables --- Source/Particles/PhysicalParticleContainer.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 27d8cc5ae..4096d1911 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -463,6 +463,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) std::size_t shared_mem_bytes = plasma_injector->sharedMemoryNeeded(); int lrrfac = rrfac; + bool loc_do_field_ionization = do_field_ionization; + int loc_ionization_initial_level = ionization_initial_level; + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). // The invalid ones are given negative ID and are deleted during the @@ -583,8 +586,8 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) u.z = gamma_boost * ( u.z -beta_boost*gamma_lab ); } - if (do_field_ionization) { - pi[ip] = ionization_initial_level; + if (loc_do_field_ionization) { + pi[ip] = loc_ionization_initial_level; } u.x *= PhysConst::c; @@ -2084,11 +2087,10 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const + ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) * ( ga *ez[i] + ux[i]*by[i] - uy[i]*bx[i] ) ); // Compute probability of ionization p - Real p; Real w_dtau = 1./ ga * p_adk_prefactor[ion_lev[i]] * std::pow(E,p_adk_power[ion_lev[i]]) * std::exp( p_adk_exp_prefactor[ion_lev[i]]/E ); - p = 1. - std::exp( - w_dtau ); + Real p = 1. - std::exp( - w_dtau ); if (random_draw < p){ // increment particle's ionization level -- cgit v1.2.3 From 9795e23c12b3f3553c748e2e335b716f3df62e22 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 08:45:25 -0700 Subject: add comments and make sure charge=q_e for ionizable species --- Source/Particles/Deposition/CurrentDeposition.H | 12 ++++++++++++ Source/Particles/MultiParticleContainer.cpp | 9 ++++++++- Source/Particles/PhysicalParticleContainer.cpp | 9 +++++++-- Source/Utils/write_atomic_data_cpp.py | 6 ++++++ 4 files changed, 33 insertions(+), 3 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 86cf89cd3..7dde90e96 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -7,6 +7,10 @@ * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx_arr : Array4 of current density, either full array or tile. * \param jy_arr : Array4 of current density, either full array or tile. * \param jz_arr : Array4 of current density, either full array or tile. @@ -37,6 +41,8 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real stagger_shift, const amrex::Real q) { + // 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 dzi = 1.0/dx[2]; @@ -166,6 +172,10 @@ void doDepositionShapeN(const amrex::Real * const xp, * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. * \param uxp uyp uzp : Pointer to arrays of particle momentum. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param Jx_arr : Array4 of current density, either full array or tile. * \param Jy_arr : Array4 of current density, either full array or tile. * \param Jz_arr : Array4 of current density, either full array or tile. @@ -195,6 +205,8 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Dim3 lo, const amrex::Real q) { + // 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; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index e39930216..9e2cd097f 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -513,7 +513,7 @@ namespace { // For particle i in mfi, if is_ionized[i]=1, copy particle // particle i from container pc_source into pc_product - static void createIonizedParticles ( + void createIonizedParticles ( int lev, const MFIter& mfi, std::unique_ptr< WarpXParticleContainer>& pc_source, std::unique_ptr< WarpXParticleContainer>& pc_product, @@ -558,6 +558,9 @@ namespace // 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; @@ -662,6 +665,10 @@ MultiParticleContainer::doFieldIonization () for (int lev = 0; lev <= pc_source->finestLevel(); ++lev){ + // When using runtime components, AMReX requires to touch all tiles + // in serial and create particles tiles with runtime components if + // they do not exist (or if they were defined by default, i.e., + // without runtime component). #ifdef _OPENMP // Touch all tiles of source species in serial if runtime attribs for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4096d1911..4562b5c0a 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1979,12 +1979,17 @@ void PhysicalParticleContainer::InitIonizationModule () { if (!do_field_ionization) return; ParmParse pp(species_name); + if (charge != PhysConst::q_e){ + amrex::Warning( + "charge != q_e for ionizable species: overriding user value and setting charge = q_e."); + charge = PhysConst::q_e; + } pp.query("ionization_initial_level", ionization_initial_level); pp.get("ionization_product_species", ionization_product_name); pp.get("physical_element", physical_element); - // Add Real component for ionization level + // Add runtime integer component for ionization level AddIntComp("ionization_level"); - plot_flags.resize(PIdx::nattribs + 1, 1); + // plot_flags.resize(PIdx::nattribs + 1, 1); // Get atomic number and ionization energies from file int ion_element_id = ion_map_ids[physical_element]; ion_atomic_number = ion_atomic_numbers[ion_element_id]; diff --git a/Source/Utils/write_atomic_data_cpp.py b/Source/Utils/write_atomic_data_cpp.py index 21d61a075..b085d50eb 100644 --- a/Source/Utils/write_atomic_data_cpp.py +++ b/Source/Utils/write_atomic_data_cpp.py @@ -1,3 +1,9 @@ +''' +This python script reads ionization tables in atomic_data.txt (generated from +the NIST website) and extracts ionization levels into C++ file +IonizationEnergiesTable.H, which contains tables + metadata. +''' + import re, os import numpy as np from scipy.constants import e -- cgit v1.2.3 From 16dd4d0491b32fb669a1ee07ed2866f68aafd6da Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 09:14:21 -0700 Subject: charge deposition depends on ionization level --- Source/Laser/LaserParticleContainer.cpp | 12 ++++++++---- Source/Particles/Deposition/ChargeDeposition.H | 13 ++++++++++++- Source/Particles/PhysicalParticleContainer.cpp | 25 +++++++++++++++++++++---- Source/Particles/WarpXParticleContainer.H | 1 + Source/Particles/WarpXParticleContainer.cpp | 23 ++++++++++++++++------- 5 files changed, 58 insertions(+), 16 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index e9a2d42de..9fac71894 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -454,9 +454,10 @@ LaserParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev = nullptr; if (crho) { - DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 0, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -529,9 +530,12 @@ LaserParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev = nullptr; + DepositCharge(pti, wp, ion_lev, rho, 1, 0, + np_current, thread_num, lev, lev); if (crho) { - DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 1, np_current, + np-np_current, thread_num, lev, lev-1); } } diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 26d42370e..4253f5e76 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -6,6 +6,10 @@ /* \brief Charge Deposition for thread thread_num * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param rho_arr : Array4 of charge density, either full array or tile. * \param np_to_depose : Number of particles for which current is deposited. * \param dx : 3D cell size @@ -18,6 +22,7 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Real * const yp, const amrex::Real * const zp, const amrex::Real * const wp, + const int * const ion_lev, const amrex::Array4& rho_arr, const long np_to_depose, const std::array& dx, @@ -25,6 +30,9 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const amrex::Dim3 lo, const amrex::Real q) { + // 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 dzi = 1.0/dx[2]; #if (AMREX_SPACEDIM == 2) @@ -43,7 +51,10 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { // --- Get particle quantities - const amrex::Real wq = q*wp[ip]*invvol; + amrex::Real wq = q*wp[ip]*invvol; + if (do_ionization){ + wq *= ion_lev[ip]; + } // --- Compute shape factors // x direction diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4562b5c0a..3db413104 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1179,9 +1179,18 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev); + + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 0, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 0, np_current, + np-np_current, thread_num, lev, lev-1); } } @@ -1327,9 +1336,17 @@ PhysicalParticleContainer::Evolve (int lev, } if (rho) { - DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + DepositCharge(pti, wp, ion_lev, rho, 1, 0, + np_current, thread_num, lev, lev); if (has_buffer){ - DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); + DepositCharge(pti, wp, ion_lev, crho, 1, np_current, + np-np_current, thread_num, lev, lev-1); } } diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 65eec3c2d..c252eee6b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -171,6 +171,7 @@ public: virtual void DepositCharge(WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, amrex::MultiFab* rho, int icomp, const long offset, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 87c2d5e90..20f8c3e00 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -562,6 +562,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, + const int * const ion_lev, MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) @@ -623,14 +624,14 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, BL_PROFILE_VAR_START(blp_ppc_chd); if (WarpX::nox == 1){ - doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<1>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 2){ - doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<2>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } else if (WarpX::nox == 3){ - doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, - np_to_depose, dx, xyzmin, lo, q); + doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, ion_lev, + rho_arr, np_to_depose, dx, xyzmin, lo, q); } BL_PROFILE_VAR_STOP(blp_ppc_chd); @@ -740,7 +741,15 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - DepositCharge(pti, wp, rho.get(), 0, 0, np, thread_num, lev, lev); + int* AMREX_RESTRICT ion_lev; + if (do_field_ionization){ + ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); + } else { + ion_lev = nullptr; + } + + DepositCharge(pti, wp, ion_lev, rho.get(), 0, 0, np, + thread_num, lev, lev); } #ifdef _OPENMP } -- cgit v1.2.3 From d871e6f89e33aa5e9d5b308d19f84e86c5383d9a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 09:23:45 -0700 Subject: add comments for charge deposition --- Source/Particles/PhysicalParticleContainer.cpp | 3 ++- Source/Particles/WarpXParticleContainer.cpp | 22 ++++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3db413104..a82c97945 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1179,7 +1179,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_STOP(blp_copy); if (rho) { - + // Deposit charge before particle push, in component 0 of MultiFab rho. int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); @@ -1336,6 +1336,7 @@ PhysicalParticleContainer::Evolve (int lev, } if (rho) { + // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 20f8c3e00..af56c58fb 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -416,6 +416,10 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti, * \param pti : Particle iterator * \param wp : Array of particle weights * \param uxp uyp uzp : Array of particle + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. * \param jx jy jz : Full array of current density * \param offset : Index of first particle for which current is deposited * \param np_to_depose: Number of particles for which current is deposited. @@ -560,6 +564,24 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #endif } +/* \brief Charge Deposition for thread thread_num + * \param pti : Particle iterator + * \param wp : Array of particle weights + * \param ion_lev : Pointer to array of particle ionization level. This is + required to have the charge of each macroparticle + since q is a scalar. For non-ionizable species, + ion_lev is a null pointer. + * \param rho : Full array of charge density + * \param icomp : Component of rho into which charge is deposited. + 0: old value (before particle push). + 1: new value (after particle push). + * \param offset : Index of first particle for which charge is deposited + * \param np_to_depose: Number of particles for which charge is deposited. + Particles [offset,offset+np_tp_depose] deposit charge + * \param thread_num : Thread number (if tiling) + * \param lev : Level of box that contains particles + * \param depos_lev : Level on which particles deposit (if buffers are used) + */ void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, const int * const ion_lev, -- cgit v1.2.3 From fe26d39a01c785b946ef42103555154c13674f9e Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 09:39:24 -0700 Subject: add back element to plot flags --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a82c97945..147f8fac6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2007,7 +2007,7 @@ void PhysicalParticleContainer::InitIonizationModule () pp.get("physical_element", physical_element); // Add runtime integer component for ionization level AddIntComp("ionization_level"); - // plot_flags.resize(PIdx::nattribs + 1, 1); + plot_flags.resize(PIdx::nattribs + 1, 1); // Get atomic number and ionization energies from file int ion_element_id = ion_map_ids[physical_element]; ion_atomic_number = ion_atomic_numbers[ion_element_id]; -- cgit v1.2.3 From b675832da2c1f96eb62c5fe1c2a1e71733296af2 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 10:17:28 -0700 Subject: fix bug in particle output --- Source/Diagnostics/ParticleIO.cpp | 6 +++--- Source/Particles/PhysicalParticleContainer.cpp | 1 - 2 files changed, 3 insertions(+), 4 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index be4d87746..743df4ce6 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -76,14 +76,15 @@ MultiParticleContainer::Checkpoint (const std::string& dir) const void MultiParticleContainer::WritePlotFile (const std::string& dir) const { - Vector int_names; - Vector int_flags; for (unsigned i = 0, n = species_names.size(); i < n; ++i) { auto& pc = allcontainers[i]; if (pc->plot_species) { Vector real_names; + Vector int_names; + Vector int_flags; + real_names.push_back("weight"); real_names.push_back("momentum_x"); @@ -121,7 +122,6 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const // when ionization is on. int_flags.resize(1, 1); } - // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 147f8fac6..18a2c9646 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2007,7 +2007,6 @@ void PhysicalParticleContainer::InitIonizationModule () pp.get("physical_element", physical_element); // Add runtime integer component for ionization level AddIntComp("ionization_level"); - plot_flags.resize(PIdx::nattribs + 1, 1); // Get atomic number and ionization energies from file int ion_element_id = ion_map_ids[physical_element]; ion_atomic_number = ion_atomic_numbers[ion_element_id]; -- cgit v1.2.3 From eb4080b05342edd63c09091f393638c9053b7b19 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Wed, 21 Aug 2019 11:31:10 -0700 Subject: reimplement the boosted frame diagnostics without the runtime components --- Source/Diagnostics/ParticleIO.cpp | 11 ---- Source/Evolve/WarpXEvolveEM.cpp | 16 +++--- Source/Particles/MultiParticleContainer.cpp | 20 ------- Source/Particles/PhysicalParticleContainer.cpp | 65 ++++++---------------- .../Particles/RigidInjectedParticleContainer.cpp | 43 +++++--------- Source/Particles/WarpXParticleContainer.H | 2 + Source/Particles/WarpXParticleContainer.cpp | 27 --------- 7 files changed, 41 insertions(+), 143 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index f159e5302..4ac5cb6a1 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -102,17 +102,6 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("theta"); #endif - if (WarpX::do_boosted_frame_diagnostic && pc->DoBoostedFrameDiags()) - { - real_names.push_back("xold"); - real_names.push_back("yold"); - real_names.push_back("zold"); - - real_names.push_back("uxold"); - real_names.push_back("uyold"); - real_names.push_back("uzold"); - } - // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 16b5905d1..ff55724f8 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -140,6 +140,14 @@ WarpX::EvolveEM (int numsteps) bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); + if (do_boosted_frame_diagnostic) { + std::unique_ptr cell_centered_data = nullptr; + if (WarpX::do_boosted_frame_fields) { + cell_centered_data = GetCellCenteredData(); + } + myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); + } + bool move_j = is_synchronized || to_make_plot || do_insitu; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. @@ -173,14 +181,6 @@ WarpX::EvolveEM (int numsteps) t_new[i] = cur_time; } - if (do_boosted_frame_diagnostic) { - std::unique_ptr cell_centered_data = nullptr; - if (WarpX::do_boosted_frame_fields) { - cell_centered_data = GetCellCenteredData(); - } - myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); - } - // slice gen // if (to_make_plot || do_insitu || to_make_slice_plot) { diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 982e04e39..4980bc62a 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -42,26 +42,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) nspecies_boosted_frame_diags += 1; } } - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - for (int i = 0; i < nspecies_boosted_frame_diags; ++i) - { - int is = map_species_boosted_frame_diags[i]; - allcontainers[is]->AddRealComp("xold"); - allcontainers[is]->AddRealComp("yold"); - allcontainers[is]->AddRealComp("zold"); - allcontainers[is]->AddRealComp("uxold"); - allcontainers[is]->AddRealComp("uyold"); - allcontainers[is]->AddRealComp("uzold"); - } - pc_tmp->AddRealComp("xold"); - pc_tmp->AddRealComp("yold"); - pc_tmp->AddRealComp("zold"); - pc_tmp->AddRealComp("uxold"); - pc_tmp->AddRealComp("uyold"); - pc_tmp->AddRealComp("uzold"); - } } void diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d10390204..14807d0d9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -196,18 +196,6 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - // need to create old values - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); - } // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); } @@ -430,15 +418,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } - GpuArray pb; - if (do_boosted) { - pb[0] = soa.GetRealData(particle_comps[ "xold"]).data() + old_size; - pb[1] = soa.GetRealData(particle_comps[ "yold"]).data() + old_size; - pb[2] = soa.GetRealData(particle_comps[ "zold"]).data() + old_size; - pb[3] = soa.GetRealData(particle_comps["uxold"]).data() + old_size; - pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; - pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; - } const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), @@ -589,15 +568,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pa[PIdx::uy][ip] = u.y; pa[PIdx::uz][ip] = u.z; - if (do_boosted) { - pb[0][ip] = x; - pb[1][ip] = y; - pb[2][ip] = z; - pb[3][ip] = u.x; - pb[4][ip] = u.y; - pb[5][ip] = u.z; - } - #if (AMREX_SPACEDIM == 3) p.pos(0) = x; p.pos(1) = y; @@ -1633,21 +1603,22 @@ PhysicalParticleContainer::PushP (int lev, Real dt, void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const Real* yp, const Real* zp) { - auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - Real* AMREX_RESTRICT xpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - Real* AMREX_RESTRICT uxpold = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - Real* AMREX_RESTRICT uypold = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - Real* AMREX_RESTRICT uzpold = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - - const long np = pti.numParticles(); + const auto np = pti.numParticles(); + const auto lev = pti.GetLevel(); + const auto index = pti.GetPairIndex(); + tmp_particle_data.resize(finestLevel()+1); + for (int i = 0; i < 6; ++i) tmp_particle_data[lev][index][i].resize(np); + Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][0].dataPtr(); + Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][1].dataPtr(); + Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][2].dataPtr(); + Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][3].dataPtr(); + Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][4].dataPtr(); + Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][5].dataPtr(); ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { @@ -1733,15 +1704,15 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uyp_new = attribs[PIdx::uy ]; auto& uzp_new = attribs[PIdx::uz ]; - auto& xp_old = pti.GetAttribs(particle_comps["xold"]); - auto& yp_old = pti.GetAttribs(particle_comps["yold"]); - auto& zp_old = pti.GetAttribs(particle_comps["zold"]); - auto& uxp_old = pti.GetAttribs(particle_comps["uxold"]); - auto& uyp_old = pti.GetAttribs(particle_comps["uyold"]); - auto& uzp_old = pti.GetAttribs(particle_comps["uzold"]); + auto& xp_old = tmp_particle_data[lev][index][0]; + auto& yp_old = tmp_particle_data[lev][index][1]; + auto& zp_old = tmp_particle_data[lev][index][2]; + auto& uxp_old = tmp_particle_data[lev][index][3]; + auto& uyp_old = tmp_particle_data[lev][index][4]; + auto& uzp_old = tmp_particle_data[lev][index][5]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 36cb9d224..0c94d1e33 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -239,15 +239,13 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Real* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); if (!done_injecting_lev) { - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - // If the old values are not already saved, create copies here. - xp_save = xp; - yp_save = yp; - zp_save = zp; - uxp_save = uxp; - uyp_save = uyp; - uzp_save = uzp; - } + // If the old values are not already saved, create copies here. + xp_save = xp; + yp_save = yp; + zp_save = zp; + uxp_save = uxp; + uyp_save = uyp; + uzp_save = uzp; // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles @@ -275,27 +273,12 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (!done_injecting_lev) { - Real* AMREX_RESTRICT x_save; - Real* AMREX_RESTRICT y_save; - Real* AMREX_RESTRICT z_save; - Real* AMREX_RESTRICT ux_save; - Real* AMREX_RESTRICT uy_save; - Real* AMREX_RESTRICT uz_save; - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - x_save = xp_save.dataPtr(); - y_save = yp_save.dataPtr(); - z_save = zp_save.dataPtr(); - ux_save = uxp_save.dataPtr(); - uy_save = uyp_save.dataPtr(); - uz_save = uzp_save.dataPtr(); - } else { - x_save = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - y_save = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - z_save = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - ux_save = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - uy_save = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - uz_save = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - } + Real* AMREX_RESTRICT x_save = xp_save.dataPtr(); + Real* AMREX_RESTRICT y_save = yp_save.dataPtr(); + Real* AMREX_RESTRICT z_save = zp_save.dataPtr(); + Real* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + Real* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + Real* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index ac5b47ada..cbdf0cdbc 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -308,6 +308,8 @@ protected: // list of names of attributes to dump. amrex::Vector plot_vars; + amrex::Vector, std::array, 6> > > tmp_particle_data; + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index befa5cfed..db0f683c7 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -85,17 +85,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - particle_comps["xold"] = PIdx::nattribs; - particle_comps["yold"] = PIdx::nattribs+1; - particle_comps["zold"] = PIdx::nattribs+2; - particle_comps["uxold"] = PIdx::nattribs+3; - particle_comps["uyold"] = PIdx::nattribs+4; - particle_comps["uzold"] = PIdx::nattribs+5; - - } - // Initialize temporary local arrays for charge/current deposition int num_threads = 1; #ifdef _OPENMP @@ -238,14 +227,6 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["xold"], x[i]); - ptile.push_back_real(particle_comps["yold"], y[i]); - ptile.push_back_real(particle_comps["zold"], z[i]); - } - particle_tile.push_back(p); } @@ -256,14 +237,6 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); - } - for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_DIM_RZ -- cgit v1.2.3 From d2a176f231788f1a04472f09bbe7134a09d14ac5 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Wed, 21 Aug 2019 11:43:28 -0700 Subject: some typedefs and an enum for clarity --- Source/Particles/PhysicalParticleContainer.cpp | 34 ++++++++++---------------- Source/Particles/WarpXParticleContainer.H | 18 +++++++++++--- 2 files changed, 28 insertions(+), 24 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 14807d0d9..c6d3cd9e7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -277,9 +277,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); GetParticles(lev)[std::make_pair(grid_id, tile_id)]; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - DefineAndReturnParticleTile(lev, grid_id, tile_id); - } } #endif @@ -403,11 +400,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int cpuid = ParallelDescriptor::MyProc(); auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - bool do_boosted = false; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - do_boosted = true; - DefineAndReturnParticleTile(lev, grid_id, tile_id); - } auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + max_new_particles; particle_tile.resize(new_size); @@ -1612,13 +1604,13 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); tmp_particle_data.resize(finestLevel()+1); - for (int i = 0; i < 6; ++i) tmp_particle_data[lev][index][i].resize(np); - Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][0].dataPtr(); - Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][1].dataPtr(); - Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][2].dataPtr(); - Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][3].dataPtr(); - Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][4].dataPtr(); - Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][5].dataPtr(); + for (int i = 0; i < TmpIdx::nattribs; ++i) tmp_particle_data[lev][index][i].resize(np); + Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); + Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); + Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); + Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); + Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); + Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { @@ -1704,12 +1696,12 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uyp_new = attribs[PIdx::uy ]; auto& uzp_new = attribs[PIdx::uz ]; - auto& xp_old = tmp_particle_data[lev][index][0]; - auto& yp_old = tmp_particle_data[lev][index][1]; - auto& zp_old = tmp_particle_data[lev][index][2]; - auto& uxp_old = tmp_particle_data[lev][index][3]; - auto& uyp_old = tmp_particle_data[lev][index][4]; - auto& uzp_old = tmp_particle_data[lev][index][5]; + auto& xp_old = tmp_particle_data[lev][index][TmpIdx::xold]; + auto& yp_old = tmp_particle_data[lev][index][TmpIdx::yold]; + auto& zp_old = tmp_particle_data[lev][index][TmpIdx::zold]; + auto& uxp_old = tmp_particle_data[lev][index][TmpIdx::uxold]; + auto& uyp_old = tmp_particle_data[lev][index][TmpIdx::uyold]; + auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index cbdf0cdbc..a21506ec3 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -29,6 +29,15 @@ struct DiagIdx }; }; +struct TmpIdx +{ + enum { + xold = 0, + yold, zold, uxold, uyold, uzold, + nattribs + }; +}; + namespace ParticleStringNames { const std::map to_index = { @@ -297,7 +306,10 @@ protected: amrex::Vector local_jy; amrex::Vector local_jz; - amrex::Vector > m_xp, m_yp, m_zp, m_giv; + using DataContainer = amrex::Gpu::ManagedDeviceVector; + using PairIndex = std::pair; + + amrex::Vector m_xp, m_yp, m_zp, m_giv; // Whether to dump particle quantities. // If true, particle position is always dumped. @@ -307,8 +319,8 @@ protected: amrex::Vector plot_flags; // list of names of attributes to dump. amrex::Vector plot_vars; - - amrex::Vector, std::array, 6> > > tmp_particle_data; + + amrex::Vector > > tmp_particle_data; private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, -- cgit v1.2.3 From d5fc8e40d0373bd4ba1e306f4ce5dc0a895ba6e7 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Wed, 21 Aug 2019 12:30:35 -0700 Subject: touch all map entries for tmp_particle_data --- Source/Particles/PhysicalParticleContainer.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c6d3cd9e7..95909feba 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -274,9 +274,10 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) #ifdef _OPENMP // First touch all tiles in the map in serial for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - GetParticles(lev)[std::make_pair(grid_id, tile_id)]; + auto index = std::make_pair(mfi.index(), mfi.LocalTileIndex()); + GetParticles(lev)[index]; + tmp_particle_data.resize(finestLevel()+1); + tmp_particle_data[lev][index]; } #endif @@ -1611,7 +1612,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); - + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { xpold[i]=xp[i]; -- cgit v1.2.3 From e71f1cc255f169652ac0d6c0138035ea732b8e65 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 21 Aug 2019 12:54:08 -0700 Subject: error if ionization on GPU --- Source/Particles/PhysicalParticleContainer.cpp | 6 ++++++ 1 file changed, 6 insertions(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 18a2c9646..1d1d6c1b7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -39,6 +39,12 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("do_field_ionization", do_field_ionization); +#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"); +#endif pp.query("plot_species", plot_species); int do_user_plot_vars; -- cgit v1.2.3 From 11d362d949107d9c51a6a9771b05d62bc8e1cd8e Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Mon, 26 Aug 2019 10:26:06 -0700 Subject: fix race condition --- Source/Particles/PhysicalParticleContainer.cpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 79a6f88f1..df8dcf836 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -949,6 +949,19 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + const auto np = pti.numParticles(); + const auto lev = pti.GetLevel(); + const auto index = pti.GetPairIndex(); + tmp_particle_data.resize(finestLevel()+1); + for (int i = 0; i < TmpIdx::nattribs; ++i) + tmp_particle_data[lev][index][i].resize(np); + } + } + #ifdef _OPENMP #pragma omp parallel #endif @@ -1680,8 +1693,6 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); - tmp_particle_data.resize(finestLevel()+1); - for (int i = 0; i < TmpIdx::nattribs; ++i) tmp_particle_data[lev][index][i].resize(np); Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); -- cgit v1.2.3 From 49dc1a600e30443c6d946689bfa80d3cd1091b1f Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 26 Aug 2019 11:06:56 -0700 Subject: add comment --- Source/Particles/PhysicalParticleContainer.cpp | 1 + 1 file changed, 1 insertion(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index df8dcf836..9169e6c1f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -295,6 +295,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) auto index = std::make_pair(mfi.index(), mfi.LocalTileIndex()); GetParticles(lev)[index]; tmp_particle_data.resize(finestLevel()+1); + // Create map entry if not there tmp_particle_data[lev][index]; if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { DefineAndReturnParticleTile(lev, mfi.index(), mfi.LocalTileIndex()); -- cgit v1.2.3