[DO NOT MERGE] Test PR 1343 Prevent plasma ionization past the atomic number

.github/workflows/linux.yml

@@ -39,14 +39,23 @@ jobs:
     - name: Build & Install
       run: |
         pip install -U -e ./tools
+        export CXXFLAGS="-fsanitize=undefined,address -fno-sanitize-recover=all"
         cmake -S . -B build                        \
+            -DCMAKE_BUILD_TYPE=RelWithDebInfo      \
             -DHiPACE_COMPUTE=OMP                   \
             -DCMAKE_CXX_STANDARD=17                \
             -DCMAKE_INSTALL_PREFIX=/tmp/my-hipace  \
             -DCMAKE_VERBOSE_MAKEFILE=ON
         cmake --build build -j 2
     - name: Run Tests
-      run: ctest --test-dir build --output-on-failure
+      run: |
+        #MPI implementations often leak memory
+        export "ASAN_OPTIONS=detect_leaks=0"
+
+        # Print stacktrace for easier debugging
+        export UBSAN_OPTIONS=print_stacktrace=1
+
+        ctest --test-dir build --output-on-failure
 
   linux_gcc_cxx17:
     name: GNU@7.5 C++17 Serial

src/particles/plasma/PlasmaParticleContainer.H

@@ -279,6 +279,7 @@ public:
     bool m_can_field_ionize = false; /**< whether this plasma can ionize from the field */
     bool m_can_laser_ionize = false; /**< whether this plasma can ionize from a laser */
     int m_init_ion_lev = -1; /**< initial Ion level of each particle */
+    int m_max_ion_lev = 0; /**< maximum ionization level */
     std::string m_product_name = ""; /**< name of Ionization product plasma */
     PlasmaParticleContainer* m_product_pc = nullptr; /**< Ionization product plasma */
     /** to calculate Ionization probability with ADK formula */

src/particles/plasma/PlasmaParticleContainer.cpp

@@ -418,6 +418,7 @@ IonizationModule (const int lev,
         amrex::Real* AMREX_RESTRICT adk_prefactor = m_adk_prefactor.data();
         amrex::Real* AMREX_RESTRICT adk_exp_prefactor = m_adk_exp_prefactor.data();
         amrex::Real* AMREX_RESTRICT adk_power = m_adk_power.data();
+        const int max_ion_lev = m_max_ion_lev;
 
         long num_ions = ptile_ion.numParticles();
 
@@ -466,6 +467,9 @@ IonizationModule (const int lev,
             const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi,
                                                            /* Assumes Aabssq == 0 */ 0._rt);
             const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip];
+            if (ion_lev_loc >= max_ion_lev) {
+                return;
+            }
             // gamma / (psi + 1) to complete dt for QSA
             amrex::Real w_dtau = gamma_psi * adk_prefactor[ion_lev_loc] *
                 std::pow(Ep, adk_power[ion_lev_loc]) *
@@ -611,6 +615,7 @@ LaserIonization (const int islice,
         amrex::Real* AMREX_RESTRICT laser_adk_prefactor = m_laser_adk_prefactor.data();
         amrex::Real* AMREX_RESTRICT laser_dp_prefactor = m_laser_dp_prefactor.data();
         amrex::Real* AMREX_RESTRICT laser_dp_second_prefactor = m_laser_dp_second_prefactor.data();
+        const int max_ion_lev = m_max_ion_lev;
 
         long num_ions = ptile_ion.numParticles();
 
@@ -660,6 +665,9 @@ LaserIonization (const int islice,
             const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi,
                                                            /* Assumes Aabssq == 0 */ 0._rt);
             const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip];
+            if (ion_lev_loc >= max_ion_lev) {
+                return;
+            }
             // gamma / (psi + 1) to complete dt for QSA
             amrex::Real w_dtau_dc = gamma_psi * adk_prefactor[ion_lev_loc] *
                 std::pow(Ep, adk_power[ion_lev_loc]) *

src/particles/plasma/PlasmaParticleContainerInit.cpp

@@ -426,6 +426,7 @@ InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_
     // Get atomic number and ionization energies from file
     const int ion_element_id = ion_map_ids[physical_element];
     const int ion_atomic_number = ion_atomic_numbers[ion_element_id];
+    m_max_ion_lev = ion_atomic_number;
     amrex::Vector<amrex::Real> h_ionization_energies(ion_atomic_number);
     const int offset = ion_energy_offsets[ion_element_id];
     for(int i=0; i<ion_atomic_number; i++){

tests/checksum/benchmarks_json/ionization.2Rank.json

@@ -1,34 +1,34 @@
 {
   "lev=0": {
-    "Bx": 2682248.7024631,
-    "By": 2693948.2701062,
-    "Bz": 461542.25627809,
-    "Ex": 2027829034919000.0,
-    "ExmBy": 1615808683503300.0,
-    "Ey": 2028075445665300.0,
-    "EypBx": 1618198281458700.0,
-    "Ez": 661283577718300.0,
-    "Psi": 10769981625.527,
-    "Sx": 7.6153143258378e+17,
-    "Sy": 7.6749661894861e+17,
-    "chi": 2095164842379900.0,
-    "jx": 5.1803326128365e+17,
-    "jx_beam": 231581099135240.0,
-    "jy": 5.166188820914e+17,
-    "jy_beam": 230138000597450.0,
-    "jz_beam": 1.6907142961473e+17,
-    "rhomjz": 8996307093.0909
+    "Bx": 2620311.1374515,
+    "By": 2635709.7531506,
+    "Bz": 414833.05353295,
+    "Ex": 2012939150389200.0,
+    "ExmBy": 1596237863501100.0,
+    "Ey": 1985459967332800.0,
+    "EypBx": 1577025752807400.0,
+    "Ez": 679659182354100.0,
+    "Psi": 10548202179.363,
+    "Sx": 7.6385240482099e+17,
+    "Sy": 7.6031507997545e+17,
+    "chi": 2108090067779600.0,
+    "jx": 5.1771260837627e+17,
+    "jx_beam": 229049530186220.0,
+    "jy": 5.1922551204841e+17,
+    "jy_beam": 224659467158600.0,
+    "jz_beam": 1.6907142196764e+17,
+    "rhomjz": 9012917997.7387
   },
   "beam": {
     "charge": 1.127932350336e-16,
     "id": 248160,
     "mass": 6.4130061345856e-28,
-    "x": 0.00024861156909874,
-    "y": 0.00024733653838739,
-    "z": 0.014359549367759,
-    "ux": 3002.9712551376,
-    "uy": 3008.114106705,
-    "uz": 1386735.4183111,
+    "x": 0.00024635902135382,
+    "y": 0.00024286960511978,
+    "z": 0.014359553378628,
+    "ux": 3013.9898183181,
+    "uy": 3014.5223974211,
+    "uz": 1386709.9483869,
     "w": 825000000.0
   }
 }

tests/checksum/benchmarks_json/laser_ionization.1Rank.json

@@ -1,41 +1,41 @@
 {
   "lev=0": {
-    "Bx": 3.5819859018319,
-    "By": 3.7269958544827,
-    "Bz": 2204.1195988099,
-    "Ex": 1117325248.1712,
+    "Bx": 3.7149515661811,
+    "By": 3.7865224175932,
+    "Bz": 2314.4198911339,
+    "Ex": 1135170862.8424,
     "ExmBy": 0.0,
-    "Ey": 1073852358.0315,
+    "Ey": 1113714461.3764,
     "EypBx": 0.0,
-    "Ez": 775468649901.41,
+    "Ez": 707934921502.67,
     "Psi": 0.0,
-    "Sx": 308477798484.84,
-    "Sy": 326849048142.86,
+    "Sx": 316753118383.02,
+    "Sy": 324030818911.92,
     "aabs": 0.46877495720552,
-    "chi": 265988685160760.0,
-    "jx": 700828092230540.0,
+    "chi": 267414764966510.0,
+    "jx": 702652162565720.0,
     "jx_beam": 0.0,
     "jy": 0.0,
     "jy_beam": 0.0,
     "jz_beam": 0.0,
-    "laserChi": 907831282357500.0,
+    "laserChi": 909177572466960.0,
     "laserEnvelope": 74.898701643489,
-    "rho_elec": 1202810276.6959,
-    "rho_ion": 1202806069.8845,
+    "rho_elec": 1209259053.802,
+    "rho_ion": 1209254830.8364,
     "rhomjz": 0.0,
-    "ux^2_elec": 0.26699411903602,
+    "ux^2_elec": 0.26780158956573,
     "ux^2_ion": 0.0,
-    "ux_elec": 83.376249717662,
+    "ux_elec": 83.593256064089,
     "ux_ion": 0.0,
     "uy^2_elec": 0.0,
     "uy^2_ion": 0.0,
     "uy_elec": 0.0,
     "uy_ion": 0.0,
-    "uz^2_elec": 1.6555627362415e-06,
-    "uz^2_ion": 1.4839081226415e-20,
-    "uz_elec": 0.14984820708954,
-    "uz_ion": 4.8429931176486e-09,
-    "w_elec": 480470478299.92,
+    "uz^2_elec": 1.6568753834647e-06,
+    "uz^2_ion": 1.5025085926483e-20,
+    "uz_elec": 0.15044101093153,
+    "uz_ion": 4.8989912126984e-09,
+    "w_elec": 483046484939.32,
     "w_ion": 1400000000000.1,
     "|a^2|": 0.46877495720552
   }