add plasma subcycling

Hi-PACE · Aug 15, 2023 · 50c10c8 · 50c10c8
1 parent b587430
commit 50c10c8
Show file tree

Hide file tree

Showing 3 changed files with 6 additions and 1 deletion.
diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H
@@ -173,6 +173,7 @@ public:
     amrex::Real m_mass = 0; /**< mass of each particle of this species */
     amrex::Real m_charge = 0; /**< charge of each particle of this species, per Ion level */
     int m_init_ion_lev = -1; /**< initial Ion level of each particle */
+    int m_n_subcycles = 1; /**< number of subcycles in the plasma particle push */
     bool m_can_ionize = false; /**< whether this plasma can ionize */
     std::string m_product_name = ""; /**< name of Ionization product plasma */
     PlasmaParticleContainer* m_product_pc = nullptr; /**< Ionization product plasma */

diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp
@@ -49,6 +49,7 @@ PlasmaParticleContainer::ReadParameters ()
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mass_Da != 0, "Unknown Element");
     }
 
+    queryWithParserAlt(pp, "n_subcycles", m_n_subcycles, pp_alt);
     queryWithParser(pp, "mass_Da", mass_Da);
     if(mass_Da != 0) {
         m_mass = phys_const.m_p * mass_Da / 1.007276466621;

diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp
@@ -66,10 +66,11 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields,
         const auto ptd = pti.GetParticleTile().getParticleTileData();
 
         const bool can_ionize = plasma.m_can_ionize;
+        const int n_subcycles = plasma.m_n_subcycles;
 
         using PTileType = PlasmaParticleContainer::ParticleTileType;
         const auto setPositionEnforceBC = EnforceBCandSetPos<PTileType>(gm[0]);
-        const amrex::Real dz = gm[0].CellSize(2);
+        const amrex::Real dz = gm[0].CellSize(2) / n_subcycles;
 
         const amrex::Real me_clight_mass_ratio = phys_const.c * phys_const.m_e/plasma.m_mass;
         const amrex::Real clight = phys_const.c;
@@ -110,6 +111,7 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields,
                     amrex::Real Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt;
                     amrex::Real Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt;
 
+                    for (int i = 0; i < n_subcycles; i++) {
                     doGatherShapeN<depos_order.value>(xp, yp, ExmByp, EypBxp, Ezp, Bxp, Byp,
                             Bzp, slice_arr, psi_comp, ez_comp, bx_comp, by_comp,
                             bz_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset);
@@ -209,6 +211,7 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields,
                     ptd.rdata(PlasmaIdx::ux)[ip] = ux;
                     ptd.rdata(PlasmaIdx::uy)[ip] = uy;
                     ptd.rdata(PlasmaIdx::psi)[ip] = psi;
+                  } // for (int isub=0; isub<nsub; ++isub) {
 #else
                     amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip];
                     amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip];