diff options
| -rw-r--r-- | Docs/source/running_cpp/parameters.rst | 12 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.H | 2 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 15 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 121 | ||||
| -rw-r--r-- | Source/Utils/WarpXConst.H | 1 |
5 files changed, 114 insertions, 37 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 904741ebf..287f534b0 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -296,7 +296,7 @@ Particle initialization * ``gaussian_beam``: Inject particle beam with gaussian distribution in space in all directions. This requires additional parameters: - ``<species_name>.q_tot`` (beam charge), + ``<species_name>.q_tot`` (beam charge) optional (default is ``q_tot=0``), ``<species_name>.npart`` (number of particles in the beam), ``<species_name>.x/y/z_m`` (average position in `x/y/z`), ``<species_name>.x/y/z_rms`` (standard deviation in `x/y/z`), @@ -305,9 +305,13 @@ Particle initialization and optional argument ``<species_name>.do_symmetrize`` (whether to symmetrize the beam in the x and y directions). - * ``external_file``: inject macroparticles with properties (charge, mass, position, and momentum) according to data in external file. - It requires the additional arguments ``<species_name>.injection_file`` and ``<species_name>.q_tot``, which are the string corresponding to the openPMD file name and the beam charge. - When using this style, it is not necessary to add other ``<species_name>.(...)`` paramters, because they will be read directly from the file. + * ``external_file``: Inject macroparticles with properties (mass, charge, position, and momentum - :math:`\gamma \beta m c`) read from an external openPMD file. + It requires the additional arguments: + ``<species_name>.injection_file`` (`string`) openPMD file name and + ``<species_name>.q_tot`` (`double`) optional (default is ``q_tot=0`` and no re-scaling is done, ``weight=q_p``) when specified it is used to re-scale the weight of externally loaded ``N`` physical particles, each of charge ``q_p``, to inject macroparticles of ``weight=<species_name>.q_tot/q_p/N``. + The external file should include the species ``openPMD::Record`` s labeled ``mass`` and ``charge`` (`double` scalars) and also the ``position`` and ``momentum`` (`double` arrays), with dimensionality and units set via ``openPMD::setUnitDimension`` and ``setUnitSI``. + The ``external_file`` option is currently implemented for 2D and 3D geometries, with record components ``x``, ``z`` and ``y`` for 3D. + For more information on the `openPMD format <https://github.com/openPMD>`__ and how to build WarpX with it, please visit :doc:`../building/openpmd`. * ``<species_name>.num_particles_per_cell_each_dim`` (`3 integers in 3D and RZ, 2 integers in 2D`) With the NUniformPerCell injection style, this specifies the number of particles along each axis diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index 278680739..0681bdf4f 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -67,7 +67,7 @@ public: amrex::Real x_cut = std::numeric_limits<amrex::Real>::max(); amrex::Real y_cut = std::numeric_limits<amrex::Real>::max(); amrex::Real z_cut = std::numeric_limits<amrex::Real>::max(); - amrex::Real q_tot; + amrex::Real q_tot = 0.0; long npart; int do_symmetrize = 0; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index ec2c454e2..e34447411 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -229,10 +229,14 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) parseDensity(pp); parseMomentum(pp); } else if (part_pos_s == "external_file") { +#ifdef WARPX_DIM_RZ + amrex::Abort("The option of reading particle data from an external " + "file has not been implemented nor tested in RZ geometry"); +#endif #ifdef WARPX_USE_OPENPMD external_file = true; - pp.get("injection_file",str_injection_file); - pp.get("q_tot",q_tot); + pp.get("injection_file", str_injection_file); + pp.query("q_tot", q_tot); #else amrex::Abort("WarpX has to be compiled with USE_OPENPMD=TRUE to be able" " to read the external openPMD file with species data"); @@ -398,7 +402,12 @@ void PlasmaInjector::parseMomentum (ParmParse& pp) makeParser(str_momentum_function_uy,{"x","y","z"}), makeParser(str_momentum_function_uz,{"x","y","z"}))); } else { - StringParseAbortMessage("Momentum distribution type", mom_dist_s); + //No need for momentum definition if external file is used + std::string s_inj_style; + pp.query("injection_style", s_inj_style); + if (s_inj_style != "external_file") { + StringParseAbortMessage("Momentum distribution type", mom_dist_s); + } } } 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; } diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index ae6cd808e..9e214576f 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -34,7 +34,6 @@ namespace PhysConst (45.*m_e*m_e*m_e*m_e*c*c*c*c*c); // SI value is 1.3050122.e-52 static constexpr amrex::Real xi_c2 = xi * c * c; // This should be usable for single precision, though // very close to smallest number possible: smallest number = 1.2e-38, xi_c2 = 1.1e-35 - static constexpr amrex::Real mevpc2_kg = 1.7826619216279e-30; // to convert MeV in kg } #endif |