aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles/PhysicalParticleContainer.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles/PhysicalParticleContainer.cpp')
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp121
1 files changed, 93 insertions, 28 deletions
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<double> disty(y_m, y_rms);
std::normal_distribution<double> distz(z_m, z_rms);
- // Allocate temporary vectors on the CPU
+ // Declare temporary vectors on the CPU
Gpu::HostVector<ParticleReal> particle_x;
Gpu::HostVector<ParticleReal> particle_y;
Gpu::HostVector<ParticleReal> 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<std::string,openPMD::ParticleSpecies> 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<amrex::Real>().get()[0];
- amrex::Real p_q = ps.second["charge"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::Real>().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<ParticleReal> particle_x;
+ Gpu::HostVector<ParticleReal> particle_z;
+ Gpu::HostVector<ParticleReal> particle_ux;
+ Gpu::HostVector<ParticleReal> particle_uz;
+ Gpu::HostVector<ParticleReal> particle_w;
+ Gpu::HostVector<ParticleReal> particle_y;
+ Gpu::HostVector<ParticleReal> 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<std::string,openPMD::ParticleSpecies> 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<amrex::ParticleReal>().get()[0];
+ double const mass_unit = ps.second["mass"][openPMD::RecordComponent::SCALAR].unitSI();
+ amrex::ParticleReal p_q = ps.second["charge"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::ParticleReal>().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<amrex::ParticleReal> ptr_x = ps.second["position"]["x"].loadChunk<amrex::ParticleReal>();
+ double const position_unit_x = ps.second["position"]["x"].unitSI();
+ std::shared_ptr<amrex::ParticleReal> ptr_z = ps.second["position"]["z"].loadChunk<amrex::ParticleReal>();
+ double const position_unit_z = ps.second["position"]["z"].unitSI();
+ std::shared_ptr<amrex::ParticleReal> ptr_ux = ps.second["momentum"]["x"].loadChunk<amrex::ParticleReal>();
+ double const momentum_unit_x = ps.second["momentum"]["x"].unitSI();
+ std::shared_ptr<amrex::ParticleReal> ptr_uz = ps.second["momentum"]["z"].loadChunk<amrex::ParticleReal>();
+ 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<amrex::ParticleReal> ptr_y = ps.second["position"]["y"].loadChunk<amrex::ParticleReal>();
+ double const position_unit_y = ps.second["position"]["y"].unitSI();
+ std::shared_ptr<amrex::ParticleReal> ptr_uy = ps.second["momentum"]["y"].loadChunk<amrex::ParticleReal>();
+ 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}; i<npart; ++i){
+ amrex::ParticleReal const x = ptr_x.get()[i]*position_unit_x;
+ amrex::ParticleReal const z = ptr_z.get()[i]*position_unit_z;
+# if (defined WARPX_DIM_2D)
+ amrex::Real const y = 0.0;
+# elif (defined WARPX_DIM_3D)
+ amrex::ParticleReal const y = ptr_y.get()[i]*position_unit_y;
+# endif
+ if (plasma_injector->insideBounds(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;
}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-13 14:25:12 +0000