From 0dd2d6f892f6c7886c6fd03f50bfa101b6494355 Mon Sep 17 00:00:00 2001 From: "L. Diana Amorim" Date: Fri, 24 Apr 2020 19:24:51 -0700 Subject: Read species distribution from OPMD - part 3 (#883) * Added .profile=external_file and .profile_file * Added description of input parameters to Docs * Changed from profile to injection option for external file * Fix typo in amrex abort message (due to copy paste) * Added the OpenPMD use amrex abort message * Minor fix - not sure how to remove EOL issue * Tried to add AddExternalFileBeam functon to PhysicalParticleContainer * Trued to fix EOL white space issue * Added read/print species name from OPMD file * Fixed OpenPMD charge and mass read * Added number of particles in species * Added nparts and converted charge/mass units to SI * Fix to nr of particles * Added q_tot parameter to determine part weight * Added macroparticle's weight * Fixed std::int typo - use only int * Added read x,y and z + ifdef 3D * Added velocity of particles and stored in container * Converted velocities to SI * No need to specfy momentum distribution if external file is used * Update Source/Initialization/PlasmaInjector.cpp Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.cpp Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.cpp Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.H Co-Authored-By: Axel Huebl * Update Source/Particles/PhysicalParticleContainer.cpp Co-Authored-By: Axel Huebl * No need to include openPMD header yet * Fixed const in AddExternalFileBeam * Fix compatibility with read_opmd and read_opmd_2 * Corrected position and charge units/sign * Added Doc note and abort message for RZ * Fix typo and EOL * Minor fixes * Fixed details. Added fix to EOL - testing * Changed to physical_charge and explained it in Docs * Fix header extra openPMD include files * Removed additional debugging comments * Try doxygen again * Fix plasmainjector physical_q_tot * Change apart to auto (not long as in gaussian beam injection style) Co-Authored-By: Axel Huebl * Correct for new type of npart Co-Authored-By: Axel Huebl * Removed Geant4 renormalization * Read openPMD file and checked its units * Trying to correct momentum information in plotfiles * Compilation error - > SegFault 11 * Added unitSI() to each position/momentum direction. Compilation SegFault * Re-structured code to use only once series.flush() * Commented out amrex::Print() s * Path to fix issue * Update after review * Chanegd to and made it optional * Fixed documentation typo and re-organized it * Added const to npart as per reviewer recommendation * Fixed issue with velocity - which became momentum * Fixed 2D and 3D options + documentation * Implemented fixes from past reviewsand fixed duplicate entry in Docs * Fix auto in iterator in for loop * Added fixes according to reviewer * Revert "Added fixes according to reviewer" This reverts commit 0a485d83c28014c8b6f53a30f26719f23f89bad5. * Fixed IO block and reordered logic with Axel * Fix spaces Co-authored-by: Axel Huebl --- Source/Particles/PhysicalParticleContainer.cpp | 121 +++++++++++++++++++------ 1 file changed, 93 insertions(+), 28 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0bfbeea16..59b6c5ca3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -254,7 +254,7 @@ PhysicalParticleContainer::AddGaussianBeam ( std::normal_distribution disty(y_m, y_rms); std::normal_distribution distz(z_m, z_rms); - // Allocate temporary vectors on the CPU + // Declare temporary vectors on the CPU Gpu::HostVector particle_x; Gpu::HostVector particle_y; Gpu::HostVector particle_z; @@ -326,36 +326,101 @@ PhysicalParticleContainer::AddGaussianBeam ( } void -PhysicalParticleContainer::AddPlasmaFromFile (const std::string s_f, amrex::Real q_tot) +PhysicalParticleContainer::AddPlasmaFromFile(const std::string s_f, + amrex::Real q_tot) { -#ifdef WARPX_USE_OPENPMD - openPMD::Series series = openPMD::Series(s_f, openPMD::AccessType::READ_ONLY); - amrex::Print() << "openPMD standard version " << series.openPMD() << "\n"; - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(series.iterations.size() == 1u, "External " - "file should contain only one iteration\n"); - openPMD::Iteration& i = series.iterations[1]; - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(i.particles.size() == 1u, "External file " - "should contain only one species\n"); - std::pair ps = *i.particles.begin(); - - //TODO: In future PRs will add AMREX_ALWAYS_ASSERT_WITH_MESSAGE to test if mass and charge are both const - amrex::Real p_m = ps.second["mass"][openPMD::RecordComponent::SCALAR].loadChunk().get()[0]; - amrex::Real p_q = ps.second["charge"][openPMD::RecordComponent::SCALAR].loadChunk().get()[0]; - int npart = ps.second["position"]["x"].getExtent()[0]; - series.flush(); - - mass = p_m*PhysConst::mevpc2_kg; - charge = p_q*PhysConst::q_e; - Real const weight = q_tot/(charge*amrex::Real(npart)); + // Declare temporary vectors on the CPU + Gpu::HostVector particle_x; + Gpu::HostVector particle_z; + Gpu::HostVector particle_ux; + Gpu::HostVector particle_uz; + Gpu::HostVector particle_w; + Gpu::HostVector particle_y; + Gpu::HostVector particle_uy; - amrex::Print() << npart << " parts of species " << ps.first << "\nWith" - << " mass = " << mass << " and charge = " << charge << "\nTo initialize " - << npart << " macroparticles with weight " << weight << "\n"; +#ifdef WARPX_USE_OPENPMD + //TODO: Make changes for read/write in multiple MPI ranks + if (ParallelDescriptor::IOProcessor()) { + openPMD::Series series = openPMD::Series(s_f, + openPMD::AccessType::READ_ONLY); + amrex::Print() << "openPMD standard version " << series.openPMD() << "\n"; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(series.iterations.size() == 1u, "External " + "file should contain only 1 iteration\n"); + openPMD::Iteration& i = series.iterations[1]; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(i.particles.size() == 1u, "External file " + "should contain only 1 species\n"); + std::pair ps = *i.particles.begin(); + + //TODO: Add ASSERT_WITH_MESSAGE to test if mass and charge are both const + amrex::ParticleReal p_m = ps.second["mass"][openPMD::RecordComponent::SCALAR].loadChunk().get()[0]; + double const mass_unit = ps.second["mass"][openPMD::RecordComponent::SCALAR].unitSI(); + amrex::ParticleReal p_q = ps.second["charge"][openPMD::RecordComponent::SCALAR].loadChunk().get()[0]; + double const charge_unit = ps.second["charge"][openPMD::RecordComponent::SCALAR].unitSI(); +# if (defined WARPX_DIM_3D) || (defined WARPX_DIM_2D) + auto const npart = ps.second["position"]["x"].getExtent()[0]; + std::shared_ptr ptr_x = ps.second["position"]["x"].loadChunk(); + double const position_unit_x = ps.second["position"]["x"].unitSI(); + std::shared_ptr ptr_z = ps.second["position"]["z"].loadChunk(); + double const position_unit_z = ps.second["position"]["z"].unitSI(); + std::shared_ptr ptr_ux = ps.second["momentum"]["x"].loadChunk(); + double const momentum_unit_x = ps.second["momentum"]["x"].unitSI(); + std::shared_ptr ptr_uz = ps.second["momentum"]["z"].loadChunk(); + double const momentum_unit_z = ps.second["momentum"]["z"].unitSI(); +# else + amrex::Abort("AddPlasmaFromFile is only implemented for 2D and 3D\n") +# endif +# if (defined WARPX_DIM_3D) + std::shared_ptr ptr_y = ps.second["position"]["y"].loadChunk(); + double const position_unit_y = ps.second["position"]["y"].unitSI(); + std::shared_ptr ptr_uy = ps.second["momentum"]["y"].loadChunk(); + double const momentum_unit_y = ps.second["momentum"]["y"].unitSI(); +# endif + series.flush(); + + mass=p_m*mass_unit; + charge=p_q*charge_unit; + + amrex::Real weight; + if (q_tot != 0.0){ + weight = q_tot/(p_q*amrex::Real(npart)); + } + else { + weight = charge; + } - amrex::Print()<<"WARNING: this is WIP, no particle has been injected!!"; -#endif + for (auto i = decltype(npart){0}; iinsideBounds(x, y, z)) { + amrex::ParticleReal const ux = ptr_ux.get()[i]*momentum_unit_x/PhysConst::m_e; + amrex::ParticleReal const uz = ptr_uz.get()[i]*momentum_unit_z/PhysConst::m_e; +# if (defined WARPX_DIM_2D) + amrex::ParticleReal const uy = 0.0; +# elif (defined WARPX_DIM_3D) + amrex::ParticleReal const uy = ptr_uy.get()[i]*momentum_unit_y/PhysConst::m_e; +# endif + CheckAndAddParticle(x, y, z, { ux, uy, uz}, weight, + particle_x, particle_y, particle_z, + particle_ux, particle_uy, particle_uz, + particle_w); + } + } + auto const np = particle_z.size(); + if (np < npart) { + amrex::Print()<<"WARNING: Simulation box doesn't cover all particles\n"; + } + } //IO Processor + auto const np = particle_z.size(); + AddNParticles(0,np, + particle_x.dataPtr(), particle_y.dataPtr(), particle_z.dataPtr(), + particle_ux.dataPtr(), particle_uy.dataPtr(), particle_uz.dataPtr(), + 1, particle_w.dataPtr(),1); +#endif //OPENPMD return; } -- cgit v1.2.3