Feature - Monte Carlo Collisions with static background neutrals (#1857)

* Update copyright notices

* allow specification of boundary potentials at runtime when using Dirichlet boundary conditions in the electrostatic solver (labframe)

* added parsing to boundary potentials specified at runtime to allow time dependence through a mathematical expression with t (time)

* updated to picmistandard 0.0.14 in order to set the electrostatic solver convergence threshold

* update docs

* various changes requested during PR review

* fixed issue causing old tests to break and added a new test for time varying boundary potentials

* possibly a fix for the failed time varying boundary condition test

* changed permission on the analysis file for the time varying BCs test

* switched to using yt for data analysis since h5py is not available

* made changes compatible with PR#1730; changed potential boundary setting routine to use the ParallelFor construct and set all boundaries in a single call

* fixed typo in computePhiRZ

* updated docs and fixed other minor typos

* fixed bug in returning from loop over dimensions when setting boundary potentials rather than continuing

* changed to setting potentials on domain boundaries rather than tilebox boundaries and changed picmi.py to accept boundary potentials

* now using domain.surroundingNodes() to get the proper boundary cells for the physical domain

* fixed typo in variable name specifying z-boundary potential

* Initial commit of MCC development. Collision type background_mcc handles collisions with a neutral background species

* added back scattering and started expanding the multiple scattering processes functionality some

* added charge exchange collision handling

* added CrossSectionHandler class to install collision process cross-section calculators

* added file reading for cross-section data

* added input parameter for energy lost during inelastic collisions and changed how secondary species are passed for ionization events

* added ionization - requires work to add to the amrex::ParallelForRNG loop

* switched the MCC ionization handling to use the same workflow as other particle creation processes i.e. using the FilterCopyTransform functionality

* updated the docs with the input parameters needed to include MCC in a run

* added test for MCC and a function to ensure that cross-section data is provided with equal energy steps

* fixed issue with build failing when USE_OMP=TRUE and some of the naming issues in Examples/Physics_applications/capacitive_discharge but I am not sure what to do about the other files in that directory

* Improve file name construction to be strictly C++ compliant

* WIP GPU Support

* Fix QED Build (CUDA 10.0)

Replace capture of a host-side array with unnamed members for E & B
field transport with a nicely named struct that transports the
Array4's as members.

This is harder to mix up and thus more self-documenting and solves an
issue with NVCC 10.0 of the form:
```
nvcc_internal_extended_lambda_implementation: In instantiation of '__nv_dl_wrapper_t<Tag, F1, F2, F3, F4, F5>::__nv_dl_wrapper_t(F1, F2, F3, F4, F5) [with Tag = __nv_dl_tag<int (*)(amrex::ParticleTile<0, 0, 4, 0, amrex::ArenaAllocator>&, amrex::ParticleTile<0, 0, 4, 0, amrex::ArenaAllocator>&, amrex::Box, const amrex::Array4<const double> (&)[6], int, int, const SchwingerFilterFunc&, const SmartCreate&, const SmartCreate&, const SchwingerTransformFunc&), filterCreateTransformFromFAB<1, amrex::ParticleTile<0, 0, 4, 0, amrex::ArenaAllocator>, amrex::Array4<const double> [6], int, const SchwingerFilterFunc&, const SmartCreate&, const SmartCreate&, const SchwingerTransformFunc&>, 1>; F1 = amrex::Array4<double>; F2 = const SchwingerFilterFunc; F3 = const amrex::Array4<const double> [6]; F4 = const amrex::Box; F5 = int*]':
/home/ubuntu/repos/WarpX/Source/Particles/ParticleCreation/FilterCreateTransformFromFAB.H:174:28:   required from 'Index filterCreateTransformFromFAB(DstTile&, DstTile&, amrex::Box, const FABs&, Index, Index, FilterFunc&&, CreateFunc1&&, CreateFunc2&&, TransFunc&&) [with int N = 1; DstTile = amrex::ParticleTile<0, 0, 4, 0, amrex::ArenaAllocator>; FABs = amrex::Array4<const double> [6]; Index = int; FilterFunc = const SchwingerFilterFunc&; CreateFunc1 = const SmartCreate&; CreateFunc2 = const SmartCreate&; TransFunc = const SchwingerTransformFunc&]'
/home/ubuntu/repos/WarpX/Source/Particles/MultiParticleContainer.cpp:1169:167:   required from here
nvcc_internal_extended_lambda_implementation:70:103: error: invalid initializer for array member 'const amrex::Array4<const double> __nv_dl_wrapper_t<__nv_dl_tag<int (*)(amrex::ParticleTile<0, 0, 4, 0, amrex::ArenaAllocator>&, amrex::ParticleTile<0, 0, 4, 0, amrex::ArenaAllocator>&, amrex::Box, const amrex::Array4<const double> (&)[6], int, int, const SchwingerFilterFunc&, const SmartCreate&, const SmartCreate&, const SchwingerTransformFunc&), filterCreateTransformFromFAB<1, amrex::ParticleTile<0, 0, 4, 0, amrex::ArenaAllocator>, amrex::Array4<const double> [6], int, const SchwingerFilterFunc&, const SmartCreate&, const SmartCreate&, const SchwingerTransformFunc&>, 1>, amrex::Array4<double>, const SchwingerFilterFunc, const amrex::Array4<const double> [6], const amrex::Box, int*>::f3 [6]'
```

* CUDA: Quiet numerous warnings about unused-variable-warning suppressions being unreachable statements

* Compiles on GPU; may even run as intended

* Delete overwrought attempt at polymorphic implementation

* Fix compilation error from nvcc being stupid

* fixed improper input file for MCC test and updated reference results - a statistical test of the MCC routine would be better so that reference results should not change with changes in the RNG

* Runs on CPU and GPU now

* Clean up GPU-related memory/allocation management and function usage

* Try inlining MCCProcess::getCrossSection to appease HIP and SYCL compilers

* Fix up style/formatting issues

* Typedef to make stuff cleaner and simpler

* MCC: Make helper functions static

* MCC: Pull parsing out to a helper

* MCC: Name member variables according to convention

* MCC: Pull out part of constructor

* MCC: Add constructor that will take any iterable source for energies/cross sections

* MCC: Overload operator new/delete to allocate in managed memory, to make later use more straightforward

* MCC: Add process type for ionization

* MCC: Expose a method for adding processes programmatically

* MCC: Follow convention of all types being 'class', which keeps grep easy

* MCC: Fix a formatting silliness

* added a check that the collision cross-section is zero at the energy penalty for the collision to ensure that no collision will happen with a particle with insufficient energy to pay the energy cost

* updated MCC input files to new standard inputs

* reverted incorrect changes that was messed up during various upstream and branch merges

* moved the MCC benchmark results to the Examples section in the documentation, which allows us to meet the style requirements - the tests are ongoing and the results will be provided in a following commit

* Add GPU synchronization after collisions

* added benchmark results and updated test results with the refined cross-sections needed to accurately calculate the benchmark cases

* removed example input files for benchmarks since the style requirement prohibits input files not included in an automated test; also updated the reference results for the MCC test which changed slightly after merging upstream development and updated amrex

* CLean up indentation and bit of commented out code

* Inline addProcess method and refactor

* Apply suggestions from code review

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* switched MCC files copyright to Modern Electron

* Remove sync calls, which are unnecessary on support modern hardware

* removed He collision cross-sections and instead the new warpx-data repository to access those files; also added a call in run_test.sh to clone the new repository during tests

* Apply suggestions from code review

Co-authored-by: David Grote <dpgrote@lbl.gov>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* cleaned up the MCC documentation a bit

* added include statements now needed by the MCC (after recent PR merges) and updated the MCC test reference values which changed slightly due to changing the value of Boltzmann's constant

* added plot results for 3rd benchmark case from literature and changed documentation to reference uploaded image rather than local image in repo

* updated MCC test file to match earlier execution which changed due to the new warpx use_filter default value being 1

* Apply suggestions from code review

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Apply suggestions from code review

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* added warp-data repository clone command to docs

* fix breaking change from earlier commit

Co-authored-by: Peter Scherpelz <peter.scherpelz@modernelectron.com>
Co-authored-by: Phil Miller <phil@intensecomputing.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: Phil Miller <phil.miller@intensecomputing.com>
Co-authored-by: Ubuntu <ubuntu@ip-172-31-31-245.us-east-2.compute.internal>
Co-authored-by: David Grote <dpgrote@lbl.gov>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

1 files changed, 382 insertions, 0 deletions

diff --git a/Source/Particles/Collision/BackgroundMCCCollision.cpp b/Source/Particles/Collision/BackgroundMCCCollision.cpp
new file mode 100644
index 000000000..c723a2800
--- /dev/null
+++ b/Source/Particles/Collision/BackgroundMCCCollision.cpp
@@ -0,0 +1,382 @@
+/* Copyright 2021 Modern Electron
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "BackgroundMCCCollision.H"
+#include "MCCScattering.H"
+#include "Particles/ParticleCreation/FilterCopyTransform.H"
+#include "Particles/ParticleCreation/SmartCopy.H"
+#include "Utils/ParticleUtils.H"
+#include "Utils/WarpXUtil.H"
+#include "Utils/WarpXProfilerWrapper.H"
+#include "WarpX.H"
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_REAL.H>
+#include <AMReX_Vector.H>
+
+#include <string>
+
+BackgroundMCCCollision::BackgroundMCCCollision (std::string const collision_name)
+    : CollisionBase(collision_name)
+{
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_species_names.size() == 1,
+                                     "Background MCC must have exactly one species.");
+
+    amrex::ParmParse pp(collision_name);
+
+    pp.query("background_density", m_background_density);
+    pp.query("background_temperature", m_background_temperature);
+
+    // if the neutral mass is specified use it, but if ionization is
+    // included the mass of the secondary species of that interaction
+    // will be used. If no neutral mass is specified and ionization is not
+    // included the mass of the colliding species will be used
+    m_background_mass = -1;
+    pp.query("background_mass", m_background_mass);
+
+    // query for a list of collision processes
+    // these could be elastic, excitation, charge_exchange, back, etc.
+    amrex::Vector<std::string> scattering_process_names;
+    pp.queryarr("scattering_processes", scattering_process_names);
+
+    // create a vector of MCCProcess objects from each scattering
+    // process name
+    for (auto scattering_process : scattering_process_names) {
+        std::string kw_cross_section = scattering_process + "_cross_section";
+        std::string cross_section_file;
+        pp.query(kw_cross_section.c_str(), cross_section_file);
+
+        amrex::Real energy = 0.0;
+        // if the scattering process is excitation or ionization get the
+        // energy associated with that process
+        if (scattering_process.find("excitation") != std::string::npos ||
+            scattering_process.find("ionization") != std::string::npos) {
+            std::string kw_energy = scattering_process + "_energy";
+            pp.get(kw_energy.c_str(), energy);
+        }
+
+        auto process = new MCCProcess(scattering_process, cross_section_file, energy);
+
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(process->m_type != MCCProcessType::INVALID,
+                                         "Cannot add an unknown MCC process type");
+
+        // if the scattering process is ionization get the secondary species
+        // only one ionization process is supported, the vector
+        // m_ionization_processes is only used to make it simple to calculate
+        // the maximum collision frequency with the same function used for
+        // particle conserving processes
+        if (process->m_type == MCCProcessType::IONIZATION) {
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!ionization_flag,
+                                             "Background MCC only supports a single ionization process");
+            ionization_flag = true;
+
+            std::string secondary_species;
+            pp.get("ionization_species", secondary_species);
+            m_species_names.push_back(secondary_species);
+
+            m_ionization_processes.push_back(process);
+        } else {
+            m_scattering_processes.push_back(process);
+        }
+
+        amrex::Gpu::synchronize();
+    }
+}
+
+/** Calculate the maximum collision frequency using a fixed energy grid that
+ *  ranges from 1e-4 to 5000 eV in 0.2 eV increments
+ */
+amrex::Real
+BackgroundMCCCollision::get_nu_max(amrex::Gpu::ManagedVector<MCCProcess*> const& mcc_processes)
+{
+    using namespace amrex::literals;
+    amrex::Real nu, nu_max = 0.0;
+
+    for (double E = 1e-4; E < 5000; E+=0.2) {
+        amrex::Real sigma_E = 0.0;
+
+        // loop through all collision pathways
+        for (const auto &scattering_process : mcc_processes) {
+            // get collision cross-section
+            sigma_E += scattering_process->getCrossSection(E);
+        }
+
+        // calculate collision frequency
+        nu = (
+              m_background_density * std::sqrt(2.0_rt / m_mass1 * PhysConst::q_e)
+              * sigma_E * std::sqrt(E)
+              );
+        if (nu > nu_max) {
+            nu_max = nu;
+        }
+    }
+    return nu_max;
+}
+
+void
+BackgroundMCCCollision::doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc)
+{
+    WARPX_PROFILE("BackgroundMCCCollision::doCollisions()");
+    using namespace amrex::literals;
+
+    const amrex::Real dt = WarpX::GetInstance().getdt(0);
+    if ( int(std::floor(cur_time/dt)) % m_ndt != 0 ) return;
+
+    auto& species1 = mypc->GetParticleContainerFromName(m_species_names[0]);
+    // this is a very ugly hack to have species2 be a reference and be
+    // defined in the scope of doCollisions
+    auto& species2 = (
+                      (m_species_names.size() == 2) ?
+                      mypc->GetParticleContainerFromName(m_species_names[1]) :
+                      mypc->GetParticleContainerFromName(m_species_names[0])
+                      );
+
+    if (!init_flag) {
+        m_mass1 = species1.getMass();
+
+        // calculate maximum collision frequency without ionization
+        m_nu_max = get_nu_max(m_scattering_processes);
+
+        // calculate total collision probability
+        auto coll_n = m_nu_max * dt;
+        m_total_collision_prob = 1.0_rt - std::exp(-coll_n);
+
+        // dt has to be small enough that a linear expansion of the collision
+        // probability is sufficiently accurately, otherwise the MCC results
+        // will be very heavily affected by small changes in the timestep
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n < 0.1_rt,
+            "dt is too large to ensure accurate MCC results"
+        );
+
+        if (ionization_flag) {
+            // calculate maximum collision frequency for ionization
+            m_nu_max_ioniz = get_nu_max(m_ionization_processes);
+
+            // calculate total ionization probability
+            auto coll_n_ioniz = m_nu_max_ioniz * dt;
+            m_total_collision_prob_ioniz = 1.0_rt - std::exp(-coll_n_ioniz);
+
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n_ioniz < 0.1_rt,
+                "dt is too large to ensure accurate MCC results"
+            );
+
+            // if an ionization process is included the secondary species mass
+            // is taken as the background mass
+            m_background_mass = species2.getMass();
+        }
+        // if no neutral species mass was specified and ionization is not
+        // included assume that the collisions will be with neutrals of the
+        // same mass as the colliding species (like with ion-neutral collisions)
+        else if (m_background_mass == -1) {
+            m_background_mass = species1.getMass();
+        }
+
+        amrex::Print() <<
+            "Setting up collisions for " << m_species_names[0] << " with total "
+            "collision probability: " <<
+            m_total_collision_prob + m_total_collision_prob_ioniz << "\n";
+
+        init_flag = true;
+    }
+
+    // Loop over refinement levels
+    auto const flvl = species1.finestLevel();
+    for (int lev = 0; lev <= flvl; ++lev) {
+
+        // firstly loop over particles box by box and do all particle conserving
+        // scattering
+#ifdef _OPENMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (WarpXParIter pti(species1, lev); pti.isValid(); ++pti) {
+            doBackgroundCollisionsWithinTile(pti);
+        }
+
+        // secondly perform ionization through the SmartCopyFactory if needed
+        if (ionization_flag) {
+            doBackgroundIonization(lev, species1, species2);
+        }
+    }
+}
+
+
+/** Perform all particle conserving MCC collisions within a tile
+ *
+ * @param pti particle iterator
+ *
+ */
+void BackgroundMCCCollision::doBackgroundCollisionsWithinTile
+( WarpXParIter& pti )
+{
+    using namespace amrex::literals;
+
+    // So that CUDA code gets its intrinsic, not the host-only C++ library version
+    using std::sqrt;
+
+    // get particle count
+    const long np = pti.numParticles();
+
+    // get collider properties
+    amrex::Real mass1 = m_mass1;
+
+    // get neutral properties
+    amrex::Real n_a = m_background_density;
+    amrex::Real T_a = m_background_temperature;
+    amrex::Real mass_a = m_background_mass;
+    amrex::Real vel_std = sqrt(PhysConst::kb * T_a / mass_a);
+
+    // get collision parameters
+    auto scattering_processes = m_scattering_processes.data();
+    auto process_count = m_scattering_processes.size();
+
+    amrex::Real total_collision_prob = m_total_collision_prob;
+    amrex::Real nu_max = m_nu_max;
+
+    // get Struct-Of-Array particle data, also called attribs
+    auto& attribs = pti.GetAttribs();
+    amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+    amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+    amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+
+    amrex::ParallelForRNG(np,
+                          [=] AMREX_GPU_HOST_DEVICE (long ip, amrex::RandomEngine const& engine)
+                          {
+                              // determine if this particle should collide
+                              if (amrex::Random(engine) > total_collision_prob) return;
+
+                              amrex::Real v_coll, v_coll2, E_coll, sigma_E, nu_i = 0;
+                              amrex::Real col_select = amrex::Random(engine);
+                              amrex::ParticleReal ua_x, ua_y, ua_z;
+                              amrex::ParticleReal uCOM_x, uCOM_y, uCOM_z;
+
+                              // get velocities of gas particles from a Maxwellian distribution
+                              ua_x = vel_std * amrex::RandomNormal(0_rt, 1.0_rt, engine);
+                              ua_y = vel_std * amrex::RandomNormal(0_rt, 1.0_rt, engine);
+                              ua_z = vel_std * amrex::RandomNormal(0_rt, 1.0_rt, engine);
+
+                              // calculate the center of momentum velocity
+                              uCOM_x = (mass1 * ux[ip] + mass_a * ua_x) / (mass1 + mass_a);
+                              uCOM_y = (mass1 * uy[ip] + mass_a * ua_y) / (mass1 + mass_a);
+                              uCOM_z = (mass1 * uz[ip] + mass_a * ua_z) / (mass1 + mass_a);
+
+                              // calculate relative velocity of collision and collision energy if
+                              // the colliding particle is an ion. For electron collisions we
+                              // cannot use the relative velocity since that allows the
+                              // possibility where the electron kinetic energy in the lab frame
+                              // is insufficient to cause excitation but not in the COM frame -
+                              // for energy to balance this situation requires the neutral to
+                              // lose energy during the collision which we don't currently
+                              // account for.
+                              if (mass_a / mass1 > 1e3) {
+                                  v_coll2 = ux[ip]*ux[ip] + uy[ip]*uy[ip] + uz[ip]*uz[ip];
+                                  E_coll = 0.5_rt * mass1 * v_coll2 / PhysConst::q_e;
+                              }
+                              else {
+                                  v_coll2 = (
+                                             (ux[ip] - ua_x)*(ux[ip] - ua_x)
+                                             + (uy[ip] - ua_y)*(uy[ip] - ua_y)
+                                             + (uz[ip] - ua_z)*(uz[ip] - ua_z)
+                                             );
+                                  E_coll = (
+                                            0.5_rt * mass1 * mass_a / (mass1 + mass_a) * v_coll2
+                                            / PhysConst::q_e
+                                            );
+                              }
+                              v_coll = sqrt(v_coll2);
+
+                              // loop through all collision pathways
+                              for (size_t i = 0; i < process_count; i++) {
+                                  auto const& scattering_process = **(scattering_processes + i);
+
+                                  // get collision cross-section
+                                  sigma_E = scattering_process.getCrossSection(E_coll);
+
+                                  // calculate normalized collision frequency
+                                  nu_i += n_a * sigma_E * v_coll / nu_max;
+
+                                  // check if this collision should be performed and call
+                                  // the appropriate scattering function
+                                  if (col_select > nu_i) continue;
+
+                                  if (scattering_process.m_type == MCCProcessType::ELASTIC) {
+                                      ElasticScattering(
+                                                        ux[ip], uy[ip], uz[ip], uCOM_x, uCOM_y, uCOM_z, engine
+                                                        );
+                                  }
+                                  else if (scattering_process.m_type == MCCProcessType::BACK) {
+                                      BackScattering(
+                                                     ux[ip], uy[ip], uz[ip], uCOM_x, uCOM_y, uCOM_z
+                                                     );
+                                  }
+                                  else if (scattering_process.m_type == MCCProcessType::CHARGE_EXCHANGE) {
+                                      ChargeExchange(ux[ip], uy[ip], uz[ip], ua_x, ua_y, ua_z);
+                                  }
+                                  else if (scattering_process.m_type == MCCProcessType::EXCITATION) {
+                                      // get the new velocity magnitude
+                                      amrex::Real vp = sqrt(
+                                                            2.0_rt / mass1 * PhysConst::q_e
+                                                            * (E_coll - scattering_process.m_energy_penalty)
+                                                            );
+                                      RandomizeVelocity(ux[ip], uy[ip], uz[ip], vp, engine);
+                                  }
+                                  break;
+                              }
+                          }
+                          );
+}
+
+
+/** Perform MCC ionization interactions
+ *
+ * @param pti particle iterator
+ * @param species1/2 reference to species container used to inject new
+ particles from ionization events
+ *
+ */
+void BackgroundMCCCollision::doBackgroundIonization
+( int lev, WarpXParticleContainer& species1,
+  WarpXParticleContainer& species2)
+{
+    WARPX_PROFILE("BackgroundMCCCollision::doBackgroundIonization()");
+
+    SmartCopyFactory copy_factory_elec(species1, species1);
+    SmartCopyFactory copy_factory_ion(species1, species2);
+    const auto CopyElec = copy_factory_elec.getSmartCopy();
+    const auto CopyIon = copy_factory_ion.getSmartCopy();
+
+    const auto Filter = ImpactIonizationFilterFunc(
+                                                   *m_ionization_processes[0],
+                                                   m_mass1, m_total_collision_prob_ioniz,
+                                                   m_nu_max_ioniz / m_background_density
+                                                   );
+
+    amrex::Real vel_std = std::sqrt(
+                                    PhysConst::kb * m_background_temperature / m_background_mass
+                                    );
+
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+    for (WarpXParIter pti(species1, lev); pti.isValid(); ++pti) {
+        auto& elec_tile = species1.ParticlesAt(lev, pti);
+        auto& ion_tile = species2.ParticlesAt(lev, pti);
+
+        const auto np_elec = elec_tile.numParticles();
+        const auto np_ion = ion_tile.numParticles();
+
+        auto Transform = ImpactIonizationTransformFunc(
+                                                       m_ionization_processes[0]->m_energy_penalty, m_mass1, vel_std
+                                                       );
+
+        const auto num_added = filterCopyTransformParticles<1>(
+                                                               elec_tile, ion_tile, elec_tile, np_elec, np_ion,
+                                                               Filter, CopyElec, CopyIon, Transform
+                                                               );
+
+        setNewParticleIDs(elec_tile, np_elec, num_added);
+        setNewParticleIDs(ion_tile, np_ion, num_added);
+    }
+}