Cleanup plasma injection (#4191)

* Extensive clean up of PlasmaInjection

* Small cleanup of parseMomentum

* Fix string argument reference

* Remove unneeded surface_flux

3 files changed, 426 insertions, 373 deletions

diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H
index d0c59f8d6..89ea4a00a 100644
--- a/Source/Initialization/PlasmaInjector.H
+++ b/Source/Initialization/PlasmaInjector.H
@@ -67,9 +67,11 @@ public:
     PhysicalSpecies getPhysicalSpecies() const {return physical_species;}
 
     // bool: whether the initial injection of particles should be done
-    // This routine is called during initialization of the plasma. When injecting
-    // a surface flux, no injection is done doing initialization so return false.
-    bool doInjection () const noexcept { return h_inj_pos != nullptr && !surface_flux;}
+    // This routine is called during initialization of the plasma.
+    bool doInjection () const noexcept { return h_inj_pos != nullptr;}
+
+    // bool: whether the flux injection of particles should be done.
+    bool doFluxInjection () const noexcept { return h_flux_pos != nullptr;}
 
     bool add_single_particle = false;
     amrex::Vector<amrex::ParticleReal> single_particle_pos;
@@ -107,7 +109,6 @@ public:
     std::unique_ptr<openPMD::Series> m_openpmd_input_series;
 #endif
 
-    bool surface_flux = false; // inject from a surface
     amrex::Real surface_flux_pos; // surface location
     amrex::Real flux_tmin = -1.; // Time after which we start injecting particles
     amrex::Real flux_tmax = -1.; // Time after which we stop injecting particles
@@ -132,6 +133,7 @@ public:
     amrex::Real density_max = std::numeric_limits<amrex::Real>::max();
 
     InjectorPosition* getInjectorPosition ();
+    InjectorPosition* getInjectorFluxPosition ();
     InjectorDensity*  getInjectorDensity ();
 
     InjectorFlux*  getInjectorFlux ();
@@ -153,6 +155,9 @@ protected:
     std::unique_ptr<InjectorPosition> h_inj_pos;
     InjectorPosition* d_inj_pos = nullptr;
 
+    std::unique_ptr<InjectorPosition> h_flux_pos;
+    InjectorPosition* d_flux_pos = nullptr;
+
     std::unique_ptr<InjectorDensity,InjectorDensityDeleter> h_inj_rho;
     InjectorDensity* d_inj_rho = nullptr;
     std::unique_ptr<amrex::Parser> density_parser;
@@ -172,9 +177,17 @@ protected:
     std::unique_ptr<TemperatureProperties> h_mom_temp;
     std::unique_ptr<VelocityProperties> h_mom_vel;
 
-    void parseDensity (const amrex::ParmParse& pp);
-    void parseFlux (const amrex::ParmParse& pp);
-    void parseMomentum (const amrex::ParmParse& pp);
+    void setupSingleParticle (const amrex::ParmParse& pp_species_name);
+    void setupMultipleParticles (const amrex::ParmParse& pp_species_name);
+    void setupGaussianBeam (const amrex::ParmParse& pp_species_name);
+    void setupNRandomPerCell (const amrex::ParmParse& pp_species_name);
+    void setupNFluxPerCell (const amrex::ParmParse& pp_species_name);
+    void setupNuniformPerCell (const amrex::ParmParse& pp_species_name);
+    void setupExternalFile (const amrex::ParmParse& pp_species_name);
+
+    void parseDensity (const amrex::ParmParse& pp_species_name);
+    void parseFlux (const amrex::ParmParse& pp_species_name);
+    void parseMomentum (const amrex::ParmParse& pp_species_name, const std::string& style);
 };
 
 #endif
diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp
index 03cee4c5e..d47517df3 100644
--- a/Source/Initialization/PlasmaInjector.cpp
+++ b/Source/Initialization/PlasmaInjector.cpp
@@ -180,68 +180,113 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name,
 
     num_particles_per_cell_each_dim.assign(3, 0);
 
-    if (injection_style == "none") {
-        return;
-    } else if (injection_style == "singleparticle") {
-        utils::parser::getArrWithParser(
-            pp_species_name, "single_particle_pos", single_particle_pos, 0, 3);
-        utils::parser::getArrWithParser(
-            pp_species_name, "single_particle_u", single_particle_u, 0, 3);
-        for (auto& x : single_particle_u) {
-            x *= PhysConst::c;
-        }
-        utils::parser::getWithParser(
-            pp_species_name, "single_particle_weight", single_particle_weight);
-        add_single_particle = true;
+    if (injection_style == "singleparticle") {
+        setupSingleParticle(pp_species_name);
         return;
     } else if (injection_style == "multipleparticles") {
-        utils::parser::getArrWithParser(
-            pp_species_name, "multiple_particles_pos_x", multiple_particles_pos_x);
-        utils::parser::getArrWithParser(
-            pp_species_name, "multiple_particles_pos_y", multiple_particles_pos_y);
-        utils::parser::getArrWithParser(
-            pp_species_name, "multiple_particles_pos_z", multiple_particles_pos_z);
-        utils::parser::getArrWithParser(
-            pp_species_name, "multiple_particles_ux", multiple_particles_ux);
-        utils::parser::getArrWithParser(
-            pp_species_name, "multiple_particles_uy", multiple_particles_uy);
-        utils::parser::getArrWithParser(
-            pp_species_name, "multiple_particles_uz", multiple_particles_uz);
-        utils::parser::getArrWithParser(
-            pp_species_name, "multiple_particles_weight", multiple_particles_weight);
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            ((multiple_particles_pos_x.size() == multiple_particles_pos_y.size()) &&
-             (multiple_particles_pos_x.size() == multiple_particles_pos_z.size()) &&
-             (multiple_particles_pos_x.size() == multiple_particles_ux.size()) &&
-             (multiple_particles_pos_x.size() == multiple_particles_uy.size()) &&
-             (multiple_particles_pos_x.size() == multiple_particles_uz.size()) &&
-             (multiple_particles_pos_x.size() == multiple_particles_weight.size())),
-            "Error: The multiple particles source quantities must all have the same number of elements");
-        for (auto& vx : multiple_particles_ux) { vx *= PhysConst::c; }
-        for (auto& vy : multiple_particles_uy) { vy *= PhysConst::c; }
-        for (auto& vz : multiple_particles_uz) { vz *= PhysConst::c; }
-        add_multiple_particles = true;
+        setupMultipleParticles(pp_species_name);
         return;
     } else if (injection_style == "gaussian_beam") {
+        setupGaussianBeam(pp_species_name);
+    } else if (injection_style == "nrandompercell") {
+        setupNRandomPerCell(pp_species_name);
+    } else if (injection_style == "nfluxpercell") {
+        setupNFluxPerCell(pp_species_name);
+    } else if (injection_style == "nuniformpercell") {
+        setupNuniformPerCell(pp_species_name);
+    } else if (injection_style == "external_file") {
+        setupExternalFile(pp_species_name);
+    } else if (injection_style != "none") {
+        StringParseAbortMessage("Injection style", injection_style);
+    }
+
+    amrex::Gpu::synchronize();
+}
 
-        utils::parser::getWithParser(pp_species_name, "x_m", x_m);
-        utils::parser::getWithParser(pp_species_name, "y_m", y_m);
-        utils::parser::getWithParser(pp_species_name, "z_m", z_m);
-        utils::parser::getWithParser(pp_species_name, "x_rms", x_rms);
-        utils::parser::getWithParser(pp_species_name, "y_rms", y_rms);
-        utils::parser::getWithParser(pp_species_name, "z_rms", z_rms);
-        utils::parser::queryWithParser(pp_species_name, "x_cut", x_cut);
-        utils::parser::queryWithParser(pp_species_name, "y_cut", y_cut);
-        utils::parser::queryWithParser(pp_species_name, "z_cut", z_cut);
-        utils::parser::getWithParser(pp_species_name, "q_tot", q_tot);
-        utils::parser::getWithParser(pp_species_name, "npart", npart);
-        pp_species_name.query("do_symmetrize", do_symmetrize);
-        pp_species_name.query("symmetrization_order", symmetrization_order);
-        const std::set<int> valid_symmetries = {4,8};
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE( valid_symmetries.count(symmetrization_order),
-            "Error: Symmetrization only supported to orders 4 or 8 ");
-        gaussian_beam = true;
-        parseMomentum(pp_species_name);
+#ifdef AMREX_USE_GPU
+PlasmaInjector::~PlasmaInjector ()
+{
+    if (d_inj_pos) {
+        amrex::The_Arena()->free(d_inj_pos);
+    }
+    if (d_flux_pos) {
+        amrex::The_Arena()->free(d_flux_pos);
+    }
+    if (d_inj_rho) {
+        amrex::The_Arena()->free(d_inj_rho);
+    }
+    if (d_inj_mom) {
+        amrex::The_Arena()->free(d_inj_mom);
+    }
+}
+#else
+PlasmaInjector::~PlasmaInjector () = default;
+#endif
+
+void PlasmaInjector::setupSingleParticle (const amrex::ParmParse& pp_species_name)
+{
+    utils::parser::getArrWithParser(
+        pp_species_name, "single_particle_pos", single_particle_pos, 0, 3);
+    utils::parser::getArrWithParser(
+        pp_species_name, "single_particle_u", single_particle_u, 0, 3);
+    for (auto& x : single_particle_u) {
+        x *= PhysConst::c;
+    }
+    utils::parser::getWithParser(
+        pp_species_name, "single_particle_weight", single_particle_weight);
+    add_single_particle = true;
+}
+
+void PlasmaInjector::setupMultipleParticles (const amrex::ParmParse& pp_species_name)
+{
+    utils::parser::getArrWithParser(
+        pp_species_name, "multiple_particles_pos_x", multiple_particles_pos_x);
+    utils::parser::getArrWithParser(
+        pp_species_name, "multiple_particles_pos_y", multiple_particles_pos_y);
+    utils::parser::getArrWithParser(
+        pp_species_name, "multiple_particles_pos_z", multiple_particles_pos_z);
+    utils::parser::getArrWithParser(
+        pp_species_name, "multiple_particles_ux", multiple_particles_ux);
+    utils::parser::getArrWithParser(
+        pp_species_name, "multiple_particles_uy", multiple_particles_uy);
+    utils::parser::getArrWithParser(
+        pp_species_name, "multiple_particles_uz", multiple_particles_uz);
+    utils::parser::getArrWithParser(
+        pp_species_name, "multiple_particles_weight", multiple_particles_weight);
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        ((multiple_particles_pos_x.size() == multiple_particles_pos_y.size()) &&
+         (multiple_particles_pos_x.size() == multiple_particles_pos_z.size()) &&
+         (multiple_particles_pos_x.size() == multiple_particles_ux.size()) &&
+         (multiple_particles_pos_x.size() == multiple_particles_uy.size()) &&
+         (multiple_particles_pos_x.size() == multiple_particles_uz.size()) &&
+         (multiple_particles_pos_x.size() == multiple_particles_weight.size())),
+        "Error: The multiple particles source quantities must all have the same number of elements");
+    for (auto& vx : multiple_particles_ux) { vx *= PhysConst::c; }
+    for (auto& vy : multiple_particles_uy) { vy *= PhysConst::c; }
+    for (auto& vz : multiple_particles_uz) { vz *= PhysConst::c; }
+    add_multiple_particles = true;
+}
+
+void PlasmaInjector::setupGaussianBeam (const amrex::ParmParse& pp_species_name)
+{
+    utils::parser::getWithParser(pp_species_name, "x_m", x_m);
+    utils::parser::getWithParser(pp_species_name, "y_m", y_m);
+    utils::parser::getWithParser(pp_species_name, "z_m", z_m);
+    utils::parser::getWithParser(pp_species_name, "x_rms", x_rms);
+    utils::parser::getWithParser(pp_species_name, "y_rms", y_rms);
+    utils::parser::getWithParser(pp_species_name, "z_rms", z_rms);
+    utils::parser::queryWithParser(pp_species_name, "x_cut", x_cut);
+    utils::parser::queryWithParser(pp_species_name, "y_cut", y_cut);
+    utils::parser::queryWithParser(pp_species_name, "z_cut", z_cut);
+    utils::parser::getWithParser(pp_species_name, "q_tot", q_tot);
+    utils::parser::getWithParser(pp_species_name, "npart", npart);
+    pp_species_name.query("do_symmetrize", do_symmetrize);
+    pp_species_name.query("symmetrization_order", symmetrization_order);
+    const std::set<int> valid_symmetries = {4,8};
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE( valid_symmetries.count(symmetrization_order),
+        "Error: Symmetrization only supported to orders 4 or 8 ");
+    gaussian_beam = true;
+    parseMomentum(pp_species_name, "gaussian_beam");
 #if defined(WARPX_DIM_XZ)
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE( y_rms > 0._rt,
         "Error: Gaussian beam y_rms must be strictly greater than 0 in 2D "
@@ -254,308 +299,278 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name,
         "Error: Gaussian beam y_rms must be strictly greater than 0 in 1D "
         "(it is used when computing the particles' weights from the total beam charge)");
 #endif
+}
 
-    }
-    // Depending on injection type at runtime, initialize inj_pos
-    // so that inj_pos->getPositionUnitBox calls
-    // InjectorPosition[Random or Regular].getPositionUnitBox.
-    else if (injection_style == "nrandompercell") {
-        utils::parser::getWithParser(
-            pp_species_name, "num_particles_per_cell", num_particles_per_cell);
+void PlasmaInjector::setupNRandomPerCell (const amrex::ParmParse& pp_species_name)
+{
+    utils::parser::getWithParser(
+        pp_species_name, "num_particles_per_cell", num_particles_per_cell);
 #if WARPX_DIM_RZ
-        if (WarpX::n_rz_azimuthal_modes > 1) {
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            num_particles_per_cell>=2*WarpX::n_rz_azimuthal_modes,
-            "Error: For accurate use of WarpX cylindrical geometry the number "
-            "of particles should be at least two times n_rz_azimuthal_modes "
-            "(Please visit PR#765 for more information.)");
-        }
+    if (WarpX::n_rz_azimuthal_modes > 1) {
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        num_particles_per_cell>=2*WarpX::n_rz_azimuthal_modes,
+        "Error: For accurate use of WarpX cylindrical geometry the number "
+        "of particles should be at least two times n_rz_azimuthal_modes "
+        "(Please visit PR#765 for more information.)");
+    }
 #endif
-        // Construct InjectorPosition with InjectorPositionRandom.
-        h_inj_pos = std::make_unique<InjectorPosition>(
-            (InjectorPositionRandom*)nullptr,
-            xmin, xmax, ymin, ymax, zmin, zmax);
-        parseDensity(pp_species_name);
-        parseMomentum(pp_species_name);
-    } else if (injection_style == "nfluxpercell") {
-        surface_flux = true;
-        utils::parser::getWithParser(
-            pp_species_name, "num_particles_per_cell", num_particles_per_cell_real);
+    // Construct InjectorPosition with InjectorPositionRandom.
+    h_inj_pos = std::make_unique<InjectorPosition>(
+        (InjectorPositionRandom*)nullptr,
+        xmin, xmax, ymin, ymax, zmin, zmax);
+#ifdef AMREX_USE_GPU
+    d_inj_pos = static_cast<InjectorPosition*>
+        (amrex::The_Arena()->alloc(sizeof(InjectorPosition)));
+    amrex::Gpu::htod_memcpy_async(d_inj_pos, h_inj_pos.get(), sizeof(InjectorPosition));
+#else
+    d_inj_pos = h_inj_pos.get();
+#endif
+
+    parseDensity(pp_species_name);
+    parseMomentum(pp_species_name, "nrandompercell");
+}
+
+void PlasmaInjector::setupNFluxPerCell (const amrex::ParmParse& pp_species_name)
+{
+    utils::parser::getWithParser(
+        pp_species_name, "num_particles_per_cell", num_particles_per_cell_real);
 #ifdef WARPX_DIM_RZ
-        if (WarpX::n_rz_azimuthal_modes > 1) {
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            num_particles_per_cell_real>=2*WarpX::n_rz_azimuthal_modes,
-            "Error: For accurate use of WarpX cylindrical geometry the number "
-            "of particles should be at least two times n_rz_azimuthal_modes "
-            "(Please visit PR#765 for more information.)");
-        }
+    if (WarpX::n_rz_azimuthal_modes > 1) {
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        num_particles_per_cell_real>=2*WarpX::n_rz_azimuthal_modes,
+        "Error: For accurate use of WarpX cylindrical geometry the number "
+        "of particles should be at least two times n_rz_azimuthal_modes "
+        "(Please visit PR#765 for more information.)");
+    }
 #endif
-        utils::parser::getWithParser(
-            pp_species_name, "surface_flux_pos", surface_flux_pos);
-        utils::parser::queryWithParser(
-            pp_species_name, "flux_tmin", flux_tmin);
-        utils::parser::queryWithParser(
-            pp_species_name, "flux_tmax", flux_tmax);
-        std::string flux_normal_axis_string;
-        pp_species_name.get("flux_normal_axis", flux_normal_axis_string);
-        flux_normal_axis = -1;
+    utils::parser::getWithParser(
+        pp_species_name, "surface_flux_pos", surface_flux_pos);
+    utils::parser::queryWithParser(
+        pp_species_name, "flux_tmin", flux_tmin);
+    utils::parser::queryWithParser(
+        pp_species_name, "flux_tmax", flux_tmax);
+    std::string flux_normal_axis_string;
+    pp_species_name.get("flux_normal_axis", flux_normal_axis_string);
+    flux_normal_axis = -1;
 #ifdef WARPX_DIM_RZ
-        if      (flux_normal_axis_string == "r" || flux_normal_axis_string == "R") {
-            flux_normal_axis = 0;
-        }
-        if      (flux_normal_axis_string == "t" || flux_normal_axis_string == "T") {
-            flux_normal_axis = 1;
-        }
+    if      (flux_normal_axis_string == "r" || flux_normal_axis_string == "R") {
+        flux_normal_axis = 0;
+    }
+    if      (flux_normal_axis_string == "t" || flux_normal_axis_string == "T") {
+        flux_normal_axis = 1;
+    }
 #else
 #    ifndef WARPX_DIM_1D_Z
-        if      (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") {
-            flux_normal_axis = 0;
-        }
+    if      (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") {
+        flux_normal_axis = 0;
+    }
 #    endif
 #endif
 #ifdef WARPX_DIM_3D
-        if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") {
-            flux_normal_axis = 1;
-        }
+    if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") {
+        flux_normal_axis = 1;
+    }
 #endif
-        if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") {
-            flux_normal_axis = 2;
-        }
+    if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") {
+        flux_normal_axis = 2;
+    }
 #ifdef WARPX_DIM_3D
-        const std::string flux_normal_axis_help = "'x', 'y', or 'z'.";
+    const std::string flux_normal_axis_help = "'x', 'y', or 'z'.";
 #else
 #    ifdef WARPX_DIM_RZ
-        const std::string flux_normal_axis_help = "'r' or 'z'.";
+    const std::string flux_normal_axis_help = "'r' or 'z'.";
 #    elif WARPX_DIM_XZ
-        const std::string flux_normal_axis_help = "'x' or 'z'.";
+    const std::string flux_normal_axis_help = "'x' or 'z'.";
 #    else
-        const std::string flux_normal_axis_help = "'z'.";
+    const std::string flux_normal_axis_help = "'z'.";
 #    endif
 #endif
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0,
-            "Error: Invalid value for flux_normal_axis. It must be " + flux_normal_axis_help);
-        pp_species_name.get("flux_direction", flux_direction);
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_direction == +1 || flux_direction == -1,
-            "Error: flux_direction must be -1 or +1.");
-        // Construct InjectorPosition with InjectorPositionRandom.
-        h_inj_pos = std::make_unique<InjectorPosition>(
-            (InjectorPositionRandomPlane*)nullptr,
-            xmin, xmax, ymin, ymax, zmin, zmax,
-            flux_normal_axis);
-        parseFlux(pp_species_name);
-        parseMomentum(pp_species_name);
-    } else if (injection_style == "nuniformpercell") {
-        // Note that for RZ, three numbers are expected, r, theta, and z.
-        // For 2D, only two are expected. The third is overwritten with 1.
-        // For 1D, only one is expected. The second and third are overwritten with 1.
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0,
+        "Error: Invalid value for flux_normal_axis. It must be " + flux_normal_axis_help);
+    pp_species_name.get("flux_direction", flux_direction);
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(flux_direction == +1 || flux_direction == -1,
+        "Error: flux_direction must be -1 or +1.");
+    // Construct InjectorPosition with InjectorPositionRandom.
+    h_flux_pos = std::make_unique<InjectorPosition>(
+        (InjectorPositionRandomPlane*)nullptr,
+        xmin, xmax, ymin, ymax, zmin, zmax,
+        flux_normal_axis);
+#ifdef AMREX_USE_GPU
+    d_flux_pos = static_cast<InjectorPosition*>
+        (amrex::The_Arena()->alloc(sizeof(InjectorPosition)));
+    amrex::Gpu::htod_memcpy_async(d_flux_pos, h_flux_pos.get(), sizeof(InjectorPosition));
+#else
+    d_flux_pos = h_flux_pos.get();
+#endif
+
+    parseFlux(pp_species_name);
+    parseMomentum(pp_species_name, "nfluxpercell");
+}
+
+void PlasmaInjector::setupNuniformPerCell (const amrex::ParmParse& pp_species_name)
+{
+    // Note that for RZ, three numbers are expected, r, theta, and z.
+    // For 2D, only two are expected. The third is overwritten with 1.
+    // For 1D, only one is expected. The second and third are overwritten with 1.
 #if defined(WARPX_DIM_1D_Z)
-        constexpr int num_required_ppc_each_dim = 1;
+    constexpr int num_required_ppc_each_dim = 1;
 #elif defined(WARPX_DIM_XZ)
-        constexpr int num_required_ppc_each_dim = 2;
+    constexpr int num_required_ppc_each_dim = 2;
 #else
-        constexpr int num_required_ppc_each_dim = 3;
+    constexpr int num_required_ppc_each_dim = 3;
 #endif
-        utils::parser::getArrWithParser(
-            pp_species_name, "num_particles_per_cell_each_dim",
-            num_particles_per_cell_each_dim, 0, num_required_ppc_each_dim);
+    utils::parser::getArrWithParser(
+        pp_species_name, "num_particles_per_cell_each_dim",
+        num_particles_per_cell_each_dim, 0, num_required_ppc_each_dim);
 #if WARPX_DIM_XZ
-        num_particles_per_cell_each_dim.push_back(1);
+    num_particles_per_cell_each_dim.push_back(1);
 #endif
 #if WARPX_DIM_1D_Z
-        num_particles_per_cell_each_dim.push_back(1); // overwrite 2nd number with 1
-        num_particles_per_cell_each_dim.push_back(1); // overwrite 3rd number with 1
+    num_particles_per_cell_each_dim.push_back(1); // overwrite 2nd number with 1
+    num_particles_per_cell_each_dim.push_back(1); // overwrite 3rd number with 1
 #endif
 #if WARPX_DIM_RZ
-        if (WarpX::n_rz_azimuthal_modes > 1) {
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            num_particles_per_cell_each_dim[1]>=2*WarpX::n_rz_azimuthal_modes,
-            "Error: For accurate use of WarpX cylindrical geometry the number "
-            "of particles in the theta direction should be at least two times "
-            "n_rz_azimuthal_modes (Please visit PR#765 for more information.)");
-        }
-#endif
-        // Construct InjectorPosition from InjectorPositionRegular.
-        h_inj_pos = std::make_unique<InjectorPosition>(
-            (InjectorPositionRegular*)nullptr,
-            xmin, xmax, ymin, ymax, zmin, zmax,
-            amrex::Dim3{num_particles_per_cell_each_dim[0],
-                num_particles_per_cell_each_dim[1],
-                num_particles_per_cell_each_dim[2]});
-        num_particles_per_cell = num_particles_per_cell_each_dim[0] *
-                                 num_particles_per_cell_each_dim[1] *
-                                 num_particles_per_cell_each_dim[2];
-        parseDensity(pp_species_name);
-        parseMomentum(pp_species_name);
-    } else if (injection_style == "external_file") {
-#ifndef WARPX_USE_OPENPMD
-        WARPX_ABORT_WITH_MESSAGE(
-            "WarpX has to be compiled with USE_OPENPMD=TRUE to be able"
-            " to read the external openPMD file with species data");
-#endif
-        external_file = true;
-        std::string str_injection_file;
-        pp_species_name.get("injection_file", str_injection_file);
-        // optional parameters
-        utils::parser::queryWithParser(pp_species_name, "q_tot", q_tot);
-        utils::parser::queryWithParser(pp_species_name, "z_shift",z_shift);
-
-#ifdef WARPX_USE_OPENPMD
-        if (amrex::ParallelDescriptor::IOProcessor()) {
-            m_openpmd_input_series = std::make_unique<openPMD::Series>(
-                str_injection_file, openPMD::Access::READ_ONLY);
-
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                m_openpmd_input_series->iterations.size() == 1u,
-                "External file should contain only 1 iteration\n");
-            openPMD::Iteration it = m_openpmd_input_series->iterations.begin()->second;
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                it.particles.size() == 1u,
-                "External file should contain only 1 species\n");
-            std::string const ps_name = it.particles.begin()->first;
-            openPMD::ParticleSpecies ps = it.particles.begin()->second;
-
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                ps.contains("charge") || charge_is_specified || species_is_specified,
-                std::string("'") + ps_name +
-                ".injection_file' does not contain a 'charge' species record. "
-                "Please specify '" + ps_name + ".charge' or "
-                "'" + ps_name + ".species_type' in your input file!\n");
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                ps.contains("mass") || mass_is_specified || species_is_specified,
-                std::string("'") + ps_name +
-                ".injection_file' does not contain a 'mass' species record. "
-                "Please specify '" + ps_name + ".mass' or "
-                "'" + ps_name + ".species_type' in your input file!\n");
-
-            if (charge_is_specified) {
-                ablastr::warn_manager::WMRecordWarning("Species",
-                    "Both '" + ps_name + ".charge' and '" +
-                        ps_name + ".injection_file' specify a charge.\n'" +
-                        ps_name + ".charge' will take precedence.\n");
-            }
-            else if (species_is_specified) {
-                ablastr::warn_manager::WMRecordWarning("Species",
-                    "Both '" + ps_name + ".species_type' and '" +
-                        ps_name + ".injection_file' specify a charge.\n'" +
-                        ps_name + ".species_type' will take precedence.\n");
-            }
-            else {
-                // TODO: Add ASSERT_WITH_MESSAGE to test if charge is a constant record
-                amrex::ParticleReal const p_q =
-                    ps["charge"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::ParticleReal>().get()[0];
-                double const charge_unit = ps["charge"][openPMD::RecordComponent::SCALAR].unitSI();
-                charge = p_q * charge_unit;
-            }
-            if (mass_is_specified) {
-                ablastr::warn_manager::WMRecordWarning("Species",
-                    "Both '" + ps_name + ".mass' and '" +
-                        ps_name + ".injection_file' specify a charge.\n'" +
-                        ps_name + ".mass' will take precedence.\n");
-            }
-            else if (species_is_specified) {
-                ablastr::warn_manager::WMRecordWarning("Species",
-                    "Both '" + ps_name + ".species_type' and '" +
-                        ps_name + ".injection_file' specify a mass.\n'" +
-                        ps_name + ".species_type' will take precedence.\n");
-            }
-            else {
-                // TODO: Add ASSERT_WITH_MESSAGE to test if mass is a constant record
-                amrex::ParticleReal const p_m =
-                    ps["mass"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::ParticleReal>().get()[0];
-                double const mass_unit = ps["mass"][openPMD::RecordComponent::SCALAR].unitSI();
-                mass = p_m * mass_unit;
-            }
-        } // IOProcessor
-
-        // Broadcast charge and mass to non-IO processors
-        if (!charge_is_specified && !species_is_specified)
-            amrex::ParallelDescriptor::Bcast(&charge, 1,
-                amrex::ParallelDescriptor::IOProcessorNumber());
-        if (!mass_is_specified && !species_is_specified)
-            amrex::ParallelDescriptor::Bcast(&mass, 1,
-                amrex::ParallelDescriptor::IOProcessorNumber());
-#else
-        WARPX_ABORT_WITH_MESSAGE(
-            "Plasma injection via external_file requires openPMD support: "
-            "Add USE_OPENPMD=TRUE when compiling WarpX.");
-#endif  // WARPX_USE_OPENPMD
-
-    } else {
-        StringParseAbortMessage("Injection style", injection_style);
+    if (WarpX::n_rz_azimuthal_modes > 1) {
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        num_particles_per_cell_each_dim[1]>=2*WarpX::n_rz_azimuthal_modes,
+        "Error: For accurate use of WarpX cylindrical geometry the number "
+        "of particles in the theta direction should be at least two times "
+        "n_rz_azimuthal_modes (Please visit PR#765 for more information.)");
     }
-
-    if (h_inj_pos) {
-#ifdef AMREX_USE_GPU
-        d_inj_pos = static_cast<InjectorPosition*>
-            (amrex::The_Arena()->alloc(sizeof(InjectorPosition)));
-        amrex::Gpu::htod_memcpy_async(d_inj_pos, h_inj_pos.get(), sizeof(InjectorPosition));
-#else
-        d_inj_pos = h_inj_pos.get();
 #endif
-    }
-
-    if (h_inj_rho) {
+    // Construct InjectorPosition from InjectorPositionRegular.
+    h_inj_pos = std::make_unique<InjectorPosition>(
+        (InjectorPositionRegular*)nullptr,
+        xmin, xmax, ymin, ymax, zmin, zmax,
+        amrex::Dim3{num_particles_per_cell_each_dim[0],
+            num_particles_per_cell_each_dim[1],
+            num_particles_per_cell_each_dim[2]});
 #ifdef AMREX_USE_GPU
-        d_inj_rho = static_cast<InjectorDensity*>
-            (amrex::The_Arena()->alloc(sizeof(InjectorDensity)));
-        amrex::Gpu::htod_memcpy_async(d_inj_rho, h_inj_rho.get(), sizeof(InjectorDensity));
+    d_inj_pos = static_cast<InjectorPosition*>
+        (amrex::The_Arena()->alloc(sizeof(InjectorPosition)));
+    amrex::Gpu::htod_memcpy_async(d_inj_pos, h_inj_pos.get(), sizeof(InjectorPosition));
 #else
-        d_inj_rho = h_inj_rho.get();
+    d_inj_pos = h_inj_pos.get();
 #endif
-    }
+    num_particles_per_cell = num_particles_per_cell_each_dim[0] *
+                             num_particles_per_cell_each_dim[1] *
+                             num_particles_per_cell_each_dim[2];
+    parseDensity(pp_species_name);
+    parseMomentum(pp_species_name, "nuniformpercell");
+}
 
-    if (h_inj_flux) {
-#ifdef AMREX_USE_GPU
-        d_inj_flux = static_cast<InjectorFlux*>
-            (amrex::The_Arena()->alloc(sizeof(InjectorFlux)));
-        amrex::Gpu::htod_memcpy_async(d_inj_flux, h_inj_flux.get(), sizeof(InjectorFlux));
-#else
-        d_inj_flux = h_inj_flux.get();
+void PlasmaInjector::setupExternalFile (const amrex::ParmParse& pp_species_name)
+{
+#ifndef WARPX_USE_OPENPMD
+    WARPX_ABORT_WITH_MESSAGE(
+        "WarpX has to be compiled with USE_OPENPMD=TRUE to be able"
+        " to read the external openPMD file with species data");
 #endif
-    }
+    external_file = true;
+    std::string str_injection_file;
+    pp_species_name.get("injection_file", str_injection_file);
+    // optional parameters
+    utils::parser::queryWithParser(pp_species_name, "q_tot", q_tot);
+    utils::parser::queryWithParser(pp_species_name, "z_shift",z_shift);
 
-    if (h_inj_mom) {
-#ifdef AMREX_USE_GPU
-        d_inj_mom = static_cast<InjectorMomentum*>
-            (amrex::The_Arena()->alloc(sizeof(InjectorMomentum)));
-        amrex::Gpu::htod_memcpy_async(d_inj_mom, h_inj_mom.get(), sizeof(InjectorMomentum));
-#else
-        d_inj_mom = h_inj_mom.get();
-#endif
-    }
+    const bool charge_is_specified = pp_species_name.contains("charge");
+    const bool mass_is_specified = pp_species_name.contains("mass");
+    const bool species_is_specified = pp_species_name.contains("species_type");
 
-    amrex::Gpu::synchronize();
-}
+#ifdef WARPX_USE_OPENPMD
+    if (amrex::ParallelDescriptor::IOProcessor()) {
+        m_openpmd_input_series = std::make_unique<openPMD::Series>(
+            str_injection_file, openPMD::Access::READ_ONLY);
 
-#ifdef AMREX_USE_GPU
-PlasmaInjector::~PlasmaInjector ()
-{
-    if (d_inj_pos) {
-        amrex::The_Arena()->free(d_inj_pos);
-    }
-    if (d_inj_rho) {
-        amrex::The_Arena()->free(d_inj_rho);
-    }
-    if (d_inj_mom) {
-        amrex::The_Arena()->free(d_inj_mom);
-    }
-}
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+            m_openpmd_input_series->iterations.size() == 1u,
+            "External file should contain only 1 iteration\n");
+        openPMD::Iteration it = m_openpmd_input_series->iterations.begin()->second;
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+            it.particles.size() == 1u,
+            "External file should contain only 1 species\n");
+        std::string const ps_name = it.particles.begin()->first;
+        openPMD::ParticleSpecies ps = it.particles.begin()->second;
+
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+            ps.contains("charge") || charge_is_specified || species_is_specified,
+            std::string("'") + ps_name +
+            ".injection_file' does not contain a 'charge' species record. "
+            "Please specify '" + ps_name + ".charge' or "
+            "'" + ps_name + ".species_type' in your input file!\n");
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+            ps.contains("mass") || mass_is_specified || species_is_specified,
+            std::string("'") + ps_name +
+            ".injection_file' does not contain a 'mass' species record. "
+            "Please specify '" + ps_name + ".mass' or "
+            "'" + ps_name + ".species_type' in your input file!\n");
+
+        if (charge_is_specified) {
+            ablastr::warn_manager::WMRecordWarning("Species",
+                "Both '" + ps_name + ".charge' and '" +
+                    ps_name + ".injection_file' specify a charge.\n'" +
+                    ps_name + ".charge' will take precedence.\n");
+        }
+        else if (species_is_specified) {
+            ablastr::warn_manager::WMRecordWarning("Species",
+                "Both '" + ps_name + ".species_type' and '" +
+                    ps_name + ".injection_file' specify a charge.\n'" +
+                    ps_name + ".species_type' will take precedence.\n");
+        }
+        else {
+            // TODO: Add ASSERT_WITH_MESSAGE to test if charge is a constant record
+            amrex::ParticleReal const p_q =
+                ps["charge"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::ParticleReal>().get()[0];
+            double const charge_unit = ps["charge"][openPMD::RecordComponent::SCALAR].unitSI();
+            charge = p_q * charge_unit;
+        }
+        if (mass_is_specified) {
+            ablastr::warn_manager::WMRecordWarning("Species",
+                "Both '" + ps_name + ".mass' and '" +
+                    ps_name + ".injection_file' specify a charge.\n'" +
+                    ps_name + ".mass' will take precedence.\n");
+        }
+        else if (species_is_specified) {
+            ablastr::warn_manager::WMRecordWarning("Species",
+                "Both '" + ps_name + ".species_type' and '" +
+                    ps_name + ".injection_file' specify a mass.\n'" +
+                    ps_name + ".species_type' will take precedence.\n");
+        }
+        else {
+            // TODO: Add ASSERT_WITH_MESSAGE to test if mass is a constant record
+            amrex::ParticleReal const p_m =
+                ps["mass"][openPMD::RecordComponent::SCALAR].loadChunk<amrex::ParticleReal>().get()[0];
+            double const mass_unit = ps["mass"][openPMD::RecordComponent::SCALAR].unitSI();
+            mass = p_m * mass_unit;
+        }
+    } // IOProcessor
+
+    // Broadcast charge and mass to non-IO processors
+    if (!charge_is_specified && !species_is_specified)
+        amrex::ParallelDescriptor::Bcast(&charge, 1,
+            amrex::ParallelDescriptor::IOProcessorNumber());
+    if (!mass_is_specified && !species_is_specified)
+        amrex::ParallelDescriptor::Bcast(&mass, 1,
+            amrex::ParallelDescriptor::IOProcessorNumber());
 #else
-PlasmaInjector::~PlasmaInjector () = default;
-#endif
+    WARPX_ABORT_WITH_MESSAGE(
+        "Plasma injection via external_file requires openPMD support: "
+        "Add USE_OPENPMD=TRUE when compiling WarpX.");
+#endif  // WARPX_USE_OPENPMD
+}
 
 // Depending on injection type at runtime, initialize inj_rho
 // so that inj_rho->getDensity calls
 // InjectorPosition[Constant or Predefined or etc.].getDensity.
-void PlasmaInjector::parseDensity (const amrex::ParmParse& pp)
+void PlasmaInjector::parseDensity (const amrex::ParmParse& pp_species_name)
 {
     // parse density information
     std::string rho_prof_s;
-    pp.get("profile", rho_prof_s);
+    pp_species_name.get("profile", rho_prof_s);
     std::transform(rho_prof_s.begin(), rho_prof_s.end(),
                    rho_prof_s.begin(), ::tolower);
     if (rho_prof_s == "constant") {
-        utils::parser::getWithParser(pp, "density", density);
+        utils::parser::getWithParser(pp_species_name, "density", density);
         // Construct InjectorDensity with InjectorDensityConstant.
         h_inj_rho.reset(new InjectorDensity((InjectorDensityConstant*)nullptr, density));
     } else if (rho_prof_s == "predefined") {
@@ -563,39 +578,44 @@ void PlasmaInjector::parseDensity (const amrex::ParmParse& pp)
         h_inj_rho.reset(new InjectorDensity((InjectorDensityPredefined*)nullptr,species_name));
     } else if (rho_prof_s == "parse_density_function") {
         utils::parser::Store_parserString(
-            pp, "density_function(x,y,z)", str_density_function);
+            pp_species_name, "density_function(x,y,z)", str_density_function);
         // Construct InjectorDensity with InjectorDensityParser.
         density_parser = std::make_unique<amrex::Parser>(
             utils::parser::makeParser(str_density_function,{"x","y","z"}));
         h_inj_rho.reset(new InjectorDensity((InjectorDensityParser*)nullptr,
             density_parser->compile<3>()));
     } else {
-        //No need for profile definition if external file is used
-        std::string injection_style = "none";
-        pp.query("injection_style", injection_style);
-        if (injection_style != "external_file") {
-            StringParseAbortMessage("Density profile type", rho_prof_s);
-        }
+        StringParseAbortMessage("Density profile type", rho_prof_s);
+    }
+
+    if (h_inj_rho) {
+#ifdef AMREX_USE_GPU
+        d_inj_rho = static_cast<InjectorDensity*>
+            (amrex::The_Arena()->alloc(sizeof(InjectorDensity)));
+        amrex::Gpu::htod_memcpy_async(d_inj_rho, h_inj_rho.get(), sizeof(InjectorDensity));
+#else
+        d_inj_rho = h_inj_rho.get();
+#endif
     }
 }
 
 // Depending on injection type at runtime, initialize inj_flux
 // so that inj_flux->getFlux calls
 // InjectorFlux[Constant or Parser or etc.].getFlux.
-void PlasmaInjector::parseFlux (const amrex::ParmParse& pp)
+void PlasmaInjector::parseFlux (const amrex::ParmParse& pp_species_name)
 {
     // parse flux information
     std::string flux_prof_s;
-    pp.get("flux_profile", flux_prof_s);
+    pp_species_name.get("flux_profile", flux_prof_s);
     std::transform(flux_prof_s.begin(), flux_prof_s.end(),
                    flux_prof_s.begin(), ::tolower);
     if (flux_prof_s == "constant") {
-        utils::parser::getWithParser(pp, "flux", flux);
+        utils::parser::getWithParser(pp_species_name, "flux", flux);
         // Construct InjectorFlux with InjectorFluxConstant.
         h_inj_flux.reset(new InjectorFlux((InjectorFluxConstant*)nullptr, flux));
     } else if (flux_prof_s == "parse_flux_function") {
         utils::parser::Store_parserString(
-            pp, "flux_function(x,y,z,t)", str_flux_function);
+            pp_species_name, "flux_function(x,y,z,t)", str_flux_function);
         // Construct InjectorFlux with InjectorFluxParser.
         flux_parser = std::make_unique<amrex::Parser>(
             utils::parser::makeParser(str_flux_function,{"x","y","z","t"}));
@@ -604,18 +624,28 @@ void PlasmaInjector::parseFlux (const amrex::ParmParse& pp)
     } else {
         StringParseAbortMessage("Flux profile type", flux_prof_s);
     }
+    if (h_inj_flux) {
+#ifdef AMREX_USE_GPU
+        d_inj_flux = static_cast<InjectorFlux*>
+            (amrex::The_Arena()->alloc(sizeof(InjectorFlux)));
+        amrex::Gpu::htod_memcpy_async(d_inj_flux, h_inj_flux.get(), sizeof(InjectorFlux));
+#else
+        d_inj_flux = h_inj_flux.get();
+#endif
+    }
+
 }
 
 // Depending on injection type at runtime, initialize inj_mom
 // so that inj_mom->getMomentum calls
 // InjectorMomentum[Constant or Gaussian or etc.].getMomentum.
-void PlasmaInjector::parseMomentum (const amrex::ParmParse& pp)
+void PlasmaInjector::parseMomentum (const amrex::ParmParse& pp_species_name, const std::string& style)
 {
     using namespace amrex::literals;
 
     // parse momentum information
     std::string mom_dist_s;
-    pp.get("momentum_distribution_type", mom_dist_s);
+    pp_species_name.get("momentum_distribution_type", mom_dist_s);
     std::transform(mom_dist_s.begin(),
                    mom_dist_s.end(),
                    mom_dist_s.begin(),
@@ -630,9 +660,9 @@ void PlasmaInjector::parseMomentum (const amrex::ParmParse& pp)
         amrex::Real ux = 0._rt;
         amrex::Real uy = 0._rt;
         amrex::Real uz = 0._rt;
-        utils::parser::queryWithParser(pp, "ux", ux);
-        utils::parser::queryWithParser(pp, "uy", uy);
-        utils::parser::queryWithParser(pp, "uz", uz);
+        utils::parser::queryWithParser(pp_species_name, "ux", ux);
+        utils::parser::queryWithParser(pp_species_name, "uy", uy);
+        utils::parser::queryWithParser(pp_species_name, "uz", uz);
         // Construct InjectorMomentum with InjectorMomentumConstant.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumConstant*)nullptr, ux, uy, uz));
     } else if (mom_dist_s == "gaussian") {
@@ -642,17 +672,17 @@ void PlasmaInjector::parseMomentum (const amrex::ParmParse& pp)
         amrex::Real ux_th = 0._rt;
         amrex::Real uy_th = 0._rt;
         amrex::Real uz_th = 0._rt;
-        utils::parser::queryWithParser(pp, "ux_m", ux_m);
-        utils::parser::queryWithParser(pp, "uy_m", uy_m);
-        utils::parser::queryWithParser(pp, "uz_m", uz_m);
-        utils::parser::queryWithParser(pp, "ux_th", ux_th);
-        utils::parser::queryWithParser(pp, "uy_th", uy_th);
-        utils::parser::queryWithParser(pp, "uz_th", uz_th);
+        utils::parser::queryWithParser(pp_species_name, "ux_m", ux_m);
+        utils::parser::queryWithParser(pp_species_name, "uy_m", uy_m);
+        utils::parser::queryWithParser(pp_species_name, "uz_m", uz_m);
+        utils::parser::queryWithParser(pp_species_name, "ux_th", ux_th);
+        utils::parser::queryWithParser(pp_species_name, "uy_th", uy_th);
+        utils::parser::queryWithParser(pp_species_name, "uz_th", uz_th);
         // Construct InjectorMomentum with InjectorMomentumGaussian.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumGaussian*)nullptr,
                                              ux_m, uy_m, uz_m, ux_th, uy_th, uz_th));
     } else if (mom_dist_s == "gaussianflux") {
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(surface_flux,
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(style == "nfluxpercell",
             "Error: gaussianflux can only be used with injection_style = NFluxPerCell");
         amrex::Real ux_m = 0._rt;
         amrex::Real uy_m = 0._rt;
@@ -660,12 +690,12 @@ void PlasmaInjector::parseMomentum (const amrex::ParmParse& pp)
         amrex::Real ux_th = 0._rt;
         amrex::Real uy_th = 0._rt;
         amrex::Real uz_th = 0._rt;
-        utils::parser::queryWithParser(pp, "ux_m", ux_m);
-        utils::parser::queryWithParser(pp, "uy_m", uy_m);
-        utils::parser::queryWithParser(pp, "uz_m", uz_m);
-        utils::parser::queryWithParser(pp, "ux_th", ux_th);
-        utils::parser::queryWithParser(pp, "uy_th", uy_th);
-        utils::parser::queryWithParser(pp, "uz_th", uz_th);
+        utils::parser::queryWithParser(pp_species_name, "ux_m", ux_m);
+        utils::parser::queryWithParser(pp_species_name, "uy_m", uy_m);
+        utils::parser::queryWithParser(pp_species_name, "uz_m", uz_m);
+        utils::parser::queryWithParser(pp_species_name, "ux_th", ux_th);
+        utils::parser::queryWithParser(pp_species_name, "uy_th", uy_th);
+        utils::parser::queryWithParser(pp_species_name, "uz_th", uz_th);
         // Construct InjectorMomentum with InjectorMomentumGaussianFlux.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumGaussianFlux*)nullptr,
                                              ux_m, uy_m, uz_m, ux_th, uy_th, uz_th,
@@ -677,41 +707,41 @@ void PlasmaInjector::parseMomentum (const amrex::ParmParse& pp)
         amrex::Real ux_max = 0._rt;
         amrex::Real uy_max = 0._rt;
         amrex::Real uz_max = 0._rt;
-        utils::parser::queryWithParser(pp, "ux_min", ux_min);
-        utils::parser::queryWithParser(pp, "uy_min", uy_min);
-        utils::parser::queryWithParser(pp, "uz_min", uz_min);
-        utils::parser::queryWithParser(pp, "ux_max", ux_max);
-        utils::parser::queryWithParser(pp, "uy_max", uy_max);
-        utils::parser::queryWithParser(pp, "uz_max", uz_max);
+        utils::parser::queryWithParser(pp_species_name, "ux_min", ux_min);
+        utils::parser::queryWithParser(pp_species_name, "uy_min", uy_min);
+        utils::parser::queryWithParser(pp_species_name, "uz_min", uz_min);
+        utils::parser::queryWithParser(pp_species_name, "ux_max", ux_max);
+        utils::parser::queryWithParser(pp_species_name, "uy_max", uy_max);
+        utils::parser::queryWithParser(pp_species_name, "uz_max", uz_max);
         // Construct InjectorMomentum with InjectorMomentumUniform.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumUniform*)nullptr,
                                              ux_min, uy_min, uz_min, ux_max, uy_max, uz_max));
     } else if (mom_dist_s == "maxwell_boltzmann"){
-        h_mom_temp = std::make_unique<TemperatureProperties>(pp);
+        h_mom_temp = std::make_unique<TemperatureProperties>(pp_species_name);
         const GetTemperature getTemp(*h_mom_temp);
-        h_mom_vel = std::make_unique<VelocityProperties>(pp);
+        h_mom_vel = std::make_unique<VelocityProperties>(pp_species_name);
         const GetVelocity getVel(*h_mom_vel);
         // Construct InjectorMomentum with InjectorMomentumBoltzmann.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumBoltzmann*)nullptr, getTemp, getVel));
     } else if (mom_dist_s == "maxwell_juttner"){
-        h_mom_temp = std::make_unique<TemperatureProperties>(pp);
+        h_mom_temp = std::make_unique<TemperatureProperties>(pp_species_name);
         const GetTemperature getTemp(*h_mom_temp);
-        h_mom_vel = std::make_unique<VelocityProperties>(pp);
+        h_mom_vel = std::make_unique<VelocityProperties>(pp_species_name);
         const GetVelocity getVel(*h_mom_vel);
         // Construct InjectorMomentum with InjectorMomentumJuttner.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumJuttner*)nullptr, getTemp, getVel));
     } else if (mom_dist_s == "radial_expansion") {
         amrex::Real u_over_r = 0._rt;
-        utils::parser::queryWithParser(pp, "u_over_r", u_over_r);
+        utils::parser::queryWithParser(pp_species_name, "u_over_r", u_over_r);
         // Construct InjectorMomentum with InjectorMomentumRadialExpansion.
         h_inj_mom.reset(new InjectorMomentum
                         ((InjectorMomentumRadialExpansion*)nullptr, u_over_r));
     } else if (mom_dist_s == "parse_momentum_function") {
-        utils::parser::Store_parserString(pp, "momentum_function_ux(x,y,z)",
+        utils::parser::Store_parserString(pp_species_name, "momentum_function_ux(x,y,z)",
             str_momentum_function_ux);
-        utils::parser::Store_parserString(pp, "momentum_function_uy(x,y,z)",
+        utils::parser::Store_parserString(pp_species_name, "momentum_function_uy(x,y,z)",
             str_momentum_function_uy);
-        utils::parser::Store_parserString(pp, "momentum_function_uz(x,y,z)",
+        utils::parser::Store_parserString(pp_species_name, "momentum_function_uz(x,y,z)",
             str_momentum_function_uz);
         // Construct InjectorMomentum with InjectorMomentumParser.
         ux_parser = std::make_unique<amrex::Parser>(
@@ -725,12 +755,16 @@ void PlasmaInjector::parseMomentum (const amrex::ParmParse& pp)
                                              uy_parser->compile<3>(),
                                              uz_parser->compile<3>()));
     } else {
-        //No need for momentum definition if external file is used
-        std::string injection_style = "none";
-        pp.query("injection_style", injection_style);
-        if (injection_style != "external_file") {
-            StringParseAbortMessage("Momentum distribution type", mom_dist_s);
-        }
+        StringParseAbortMessage("Momentum distribution type", mom_dist_s);
+    }
+    if (h_inj_mom) {
+#ifdef AMREX_USE_GPU
+        d_inj_mom = static_cast<InjectorMomentum*>
+            (amrex::The_Arena()->alloc(sizeof(InjectorMomentum)));
+        amrex::Gpu::htod_memcpy_async(d_inj_mom, h_inj_mom.get(), sizeof(InjectorMomentum));
+#else
+        d_inj_mom = h_inj_mom.get();
+#endif
     }
 }
 
@@ -762,6 +796,12 @@ PlasmaInjector::getInjectorPosition ()
     return d_inj_pos;
 }
 
+InjectorPosition*
+PlasmaInjector::getInjectorFluxPosition ()
+{
+    return d_flux_pos;
+}
+
 InjectorDensity*
 PlasmaInjector::getInjectorDensity ()
 {
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 7fd85833e..958c800ff 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -1499,7 +1499,7 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
         fine_injection_box.coarsen(rrfac);
     }
 
-    InjectorPosition* inj_pos = plasma_injector->getInjectorPosition();
+    InjectorPosition* flux_pos = plasma_injector->getInjectorFluxPosition();
     InjectorFlux*  inj_flux = plasma_injector->getInjectorFlux();
     InjectorMomentum* inj_mom = plasma_injector->getInjectorMomentumDevice();
     constexpr int level_zero = 0;
@@ -1619,7 +1619,7 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
 
             const int num_ppc_int = static_cast<int>(num_ppc_real + amrex::Random(engine));
 
-            if (inj_pos->overlapsWith(lo, hi))
+            if (flux_pos->overlapsWith(lo, hi))
             {
                 auto index = overlap_box.index(iv);
                 int r;
@@ -1767,12 +1767,12 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
                 p.id() = pid+ip;
                 p.cpu() = cpuid;
 
-                // This assumes the inj_pos is of type InjectorPositionRandomPlane
+                // This assumes the flux_pos is of type InjectorPositionRandomPlane
                 const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ?
                   // In the refined injection region: use refinement ratio `lrrfac`
-                  inj_pos->getPositionUnitBox(i_part, lrrfac, engine) :
+                  flux_pos->getPositionUnitBox(i_part, lrrfac, engine) :
                   // Otherwise: use 1 as the refinement ratio
-                  inj_pos->getPositionUnitBox(i_part, amrex::IntVect::TheUnitVector(), engine);
+                  flux_pos->getPositionUnitBox(i_part, amrex::IntVect::TheUnitVector(), engine);
                 auto pos = getCellCoords(overlap_corner, dx, r, iv);
                 auto ppos = PDim3(pos);
 
@@ -1809,7 +1809,7 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
                 // Lab-frame simulation
                 // If the particle's initial position is not within or on the species's
                 // xmin, xmax, ymin, ymax, zmin, zmax, go to the next generated particle.
-                if (!inj_pos->insideBoundsInclusive(ppos.x, ppos.y, ppos.z)) {
+                if (!flux_pos->insideBoundsInclusive(ppos.x, ppos.y, ppos.z)) {
                     p.id() = -1;
                     continue;
                 }
@@ -2620,7 +2620,7 @@ PhysicalParticleContainer::ContinuousInjection (const RealBox& injection_box)
 void
 PhysicalParticleContainer::ContinuousFluxInjection (amrex::Real t, amrex::Real dt)
 {
-    if (plasma_injector->surface_flux){
+    if (plasma_injector->doFluxInjection()){
         // Check the optional parameters for start and stop of injection
         if ( ((plasma_injector->flux_tmin<0) || (t>=plasma_injector->flux_tmin)) &&
              ((plasma_injector->flux_tmax<0) || (t< plasma_injector->flux_tmax)) ){