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, 129 insertions, 0 deletions

diff --git a/Source/Particles/Collision/MCCProcess.cpp b/Source/Particles/Collision/MCCProcess.cpp
new file mode 100644
index 000000000..75778ce2f
--- /dev/null
+++ b/Source/Particles/Collision/MCCProcess.cpp
@@ -0,0 +1,129 @@
+/* Copyright 2021 Modern Electron
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "MCCProcess.H"
+#include "WarpX.H"
+
+MCCProcess::MCCProcess (
+                        const std::string& scattering_process,
+                        const std::string& cross_section_file,
+                        const amrex::Real energy )
+    : m_type(parseProcessType(scattering_process))
+    , m_energy_penalty(energy)
+{
+    amrex::Print() << "Reading file " << cross_section_file << " for "
+                   << scattering_process << " scattering cross-sections.\n";
+
+    // read the cross-section data file into memory
+    readCrossSectionFile(cross_section_file, m_energies, m_sigmas);
+
+    init();
+}
+
+template <typename InputVector>
+MCCProcess::MCCProcess (
+                        const std::string& scattering_process,
+                        const InputVector&& energies,
+                        const InputVector&& sigmas,
+                        const amrex::Real energy )
+    : m_type(parseProcessType(scattering_process))
+    , m_energy_penalty(energy)
+{
+    using std::begin;
+    using std::end;
+    m_energies.insert(m_energies.begin(), begin(energies), end(energies));
+    m_sigmas  .insert(m_sigmas.begin(),   begin(sigmas),   end(sigmas));
+
+    init();
+}
+
+void*
+MCCProcess::operator new ( std::size_t count )
+{
+    return amrex::The_Managed_Arena()->alloc(count);
+}
+
+void
+MCCProcess::operator delete ( void* ptr )
+{
+    amrex::The_Managed_Arena()->free(ptr);
+}
+
+void
+MCCProcess::init()
+{
+    m_energies_data = m_energies.data();
+    m_sigmas_data = m_sigmas.data();
+
+    // save energy grid parameters for easy use
+    m_grid_size = m_energies.size();
+    m_energy_lo = m_energies[0];
+    m_energy_hi = m_energies[m_grid_size-1];
+    m_sigma_lo = m_sigmas[0];
+    m_sigma_hi = m_sigmas[m_grid_size-1];
+    m_dE = (m_energy_hi - m_energy_lo)/(m_grid_size - 1.);
+
+    // sanity check cross-section energy grid
+    sanityCheckEnergyGrid(m_energies, m_dE);
+
+    // check that the cross-section is 0 at the energy cost if the energy
+    // cost is > 0 - this is to prevent the possibility of negative left
+    // over energy after a collision event
+    if (m_energy_penalty > 0) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            (getCrossSection(m_energy_penalty) == 0),
+            "Cross-section > 0 at energy cost for collision."
+        );
+    }
+}
+
+MCCProcessType
+MCCProcess::parseProcessType(const std::string& scattering_process)
+{
+    if (scattering_process == "elastic") {
+        return MCCProcessType::ELASTIC;
+    } else if (scattering_process == "back") {
+        return MCCProcessType::BACK;
+    } else if (scattering_process == "charge_exchange") {
+        return MCCProcessType::CHARGE_EXCHANGE;
+    } else if (scattering_process == "ionization") {
+        return MCCProcessType::IONIZATION;
+    } else if (scattering_process.find("excitation") != std::string::npos) {
+        return MCCProcessType::EXCITATION;
+    } else {
+        return MCCProcessType::INVALID;
+    }
+}
+
+void
+MCCProcess::readCrossSectionFile (
+                                  const std::string cross_section_file,
+                                  MCCProcess::VectorType<amrex::Real>& energies,
+                                  MCCProcess::VectorType<amrex::Real>& sigmas )
+{
+    std::ifstream infile(cross_section_file);
+    double energy, sigma;
+    while (infile >> energy >> sigma) {
+        energies.push_back(energy);
+        sigmas.push_back(sigma);
+    }
+}
+
+void
+MCCProcess::sanityCheckEnergyGrid (
+                                   const VectorType<amrex::Real>& energies,
+                                   amrex::Real dE
+                                   )
+{
+    // confirm that the input data for the cross-section was provided with
+    // equal energy steps, otherwise the linear interpolation will fail
+    for (unsigned i = 1; i < energies.size(); i++) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+                                         (std::abs(energies[i] - energies[i-1] - dE) < dE / 100.0),
+                                         "Energy grid not evenly spaced."
+                                         );
+    }
+}