Reconfigured the collision classes to allow for generalization (#1583)

* Reconfigured the collision classes to allow for generalization
* Added literals to PairWiseCoulombCollision
* Minor cleanup of PairWiseCoulombCollision.cpp, removing query of ndt
* Add boilerplate to CollisionBase class
* Fixed white space in CollisionBase.H

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

10 files changed, 289 insertions, 113 deletions

diff --git a/Source/Particles/Collision/CMakeLists.txt b/Source/Particles/Collision/CMakeLists.txt
index f77a13b2c..6fbcaa029 100644
--- a/Source/Particles/Collision/CMakeLists.txt
+++ b/Source/Particles/Collision/CMakeLists.txt
@@ -1,4 +1,7 @@
 target_sources(WarpX
   PRIVATE
-    CollisionType.cpp
+    CollisionHandler.cpp
+    CollisionBase.cpp
+    PairWiseCoulombCollision.cpp
 )
+
diff --git a/Source/Particles/Collision/CollisionBase.H b/Source/Particles/Collision/CollisionBase.H
new file mode 100644
index 000000000..f40476847
--- /dev/null
+++ b/Source/Particles/Collision/CollisionBase.H
@@ -0,0 +1,36 @@
+/* Copyright 2020 David Grote
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_COLLISIONBASE_H_
+#define WARPX_PARTICLES_COLLISION_COLLISIONBASE_H_
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_Vector.H>
+
+// Forward declaration needed
+class MultiParticleContainer;
+
+class CollisionBase
+{
+public:
+
+    CollisionBase (std::string collision_name);
+
+    virtual void doCollisions (amrex::Real /*cur_time*/, MultiParticleContainer* /*mypc*/ ){}
+
+    CollisionBase(CollisionBase const &) = delete;
+    CollisionBase(CollisionBase &&) = delete;
+    CollisionBase & operator=(CollisionBase const &) = delete;
+
+    virtual ~CollisionBase() = default;
+protected:
+
+    amrex::Vector<std::string> m_species_names;
+    int m_ndt;
+
+};
+
+#endif // WARPX_PARTICLES_COLLISION_COLLISIONBASE_H_
diff --git a/Source/Particles/Collision/CollisionBase.cpp b/Source/Particles/Collision/CollisionBase.cpp
new file mode 100644
index 000000000..c56586af6
--- /dev/null
+++ b/Source/Particles/Collision/CollisionBase.cpp
@@ -0,0 +1,21 @@
+/* Copyright 2020 David Grote
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "CollisionBase.H"
+
+CollisionBase::CollisionBase (std::string collision_name)
+{
+
+    // read collision species
+    amrex::ParmParse pp(collision_name);
+    pp.getarr("species", m_species_names);
+
+    // number of time steps between collisions
+    m_ndt = 1;
+    pp.query("ndt", m_ndt);
+
+}
+
diff --git a/Source/Particles/Collision/CollisionHandler.H b/Source/Particles/Collision/CollisionHandler.H
new file mode 100644
index 000000000..f4338e3bf
--- /dev/null
+++ b/Source/Particles/Collision/CollisionHandler.H
@@ -0,0 +1,35 @@
+/* Copyright 2020 David Grote
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_COLLISIONHANDLER_H_
+#define WARPX_PARTICLES_COLLISION_COLLISIONHANDLER_H_
+
+#include "CollisionBase.H"
+#include <AMReX_Vector.H>
+
+// Forward declaration needed
+class MultiParticleContainer;
+
+/* \brief CollisionHandler is a light weight class that contains the
+ * list of collisions to be done.
+ */
+class CollisionHandler
+{
+public:
+    CollisionHandler ();
+
+    /* Perform all of the collisions */
+    void doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc);
+
+private:
+
+    amrex::Vector<std::string> collision_names;
+    amrex::Vector<std::string> collision_types;
+    amrex::Vector< std::unique_ptr<CollisionBase> > allcollisions;
+
+};
+
+#endif // WARPX_PARTICLES_COLLISION_COLLISIONHANDLER_H_
diff --git a/Source/Particles/Collision/CollisionHandler.cpp b/Source/Particles/Collision/CollisionHandler.cpp
new file mode 100644
index 000000000..588e8d008
--- /dev/null
+++ b/Source/Particles/Collision/CollisionHandler.cpp
@@ -0,0 +1,52 @@
+/* Copyright 2020 David Grote
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "CollisionHandler.H"
+#include "PairWiseCoulombCollision.H"
+#include <AMReX_ParmParse.H>
+
+CollisionHandler::CollisionHandler()
+{
+
+    // Read in collision input
+    amrex::ParmParse pp_collisions("collisions");
+    pp_collisions.queryarr("collision_names", collision_names);
+
+    // Create instances based on the collision type
+    auto const ncollisions = collision_names.size();
+    collision_types.resize(ncollisions);
+    allcollisions.resize(ncollisions);
+    for (int i = 0; i < static_cast<int>(ncollisions); ++i) {
+        amrex::ParmParse pp_collision_name(collision_names[i]);
+
+        // For legacy, pairwisecoulomb is the default
+        std::string type = "pairwisecoulomb";
+
+        pp_collision_name.query("type", type);
+        collision_types[i] = type;
+
+        if (type == "pairwisecoulomb") {
+            allcollisions[i] = std::make_unique<PairWiseCoulombCollision>(collision_names[i]);
+        }
+
+    }
+
+}
+
+/** Perform all collisions
+ *
+ * @param cur_time Current time
+ * @param mypc MultiParticleContainer calling this method
+ *
+ */
+void CollisionHandler::doCollisions ( amrex::Real cur_time, MultiParticleContainer* mypc)
+{
+
+    for (auto& collision : allcollisions) {
+        collision->doCollisions( cur_time, mypc);
+    }
+
+}
diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H
deleted file mode 100644
index f1ffbe439..000000000
--- a/Source/Particles/Collision/CollisionType.H
+++ /dev/null
@@ -1,47 +0,0 @@
-/* Copyright 2019 Yinjian Zhao
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-#ifndef WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_
-#define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_
-
-#include "Particles/WarpXParticleContainer.H"
-#include <AMReX_DenseBins.H>
-#include <AMReX_REAL.H>
-#include <AMReX_ParmParse.H>
-
-class CollisionType
-{
-public:
-    int  m_species1_index;
-    int  m_species2_index;
-    int  m_ndt;
-    bool m_isSameSpecies;
-    amrex::Real m_CoulombLog;
-
-    CollisionType(
-        const std::vector<std::string>& species_names,
-        std::string const collision_name);
-
-    /** Perform all binary collisions within a tile
-     *
-     * @param lev AMR level of the tile
-     * @param mfi iterator for multifab
-     * @param species1/2 pointer to species container
-     * @param isSameSpecies true if collision is between same species
-     * @param CoulombLog user input Coulomb logrithm
-     *
-     */
-
-    static void doCoulombCollisionsWithinTile (
-        int const lev, amrex::MFIter const& mfi,
-        std::unique_ptr<WarpXParticleContainer>& species1,
-        std::unique_ptr<WarpXParticleContainer>& species2,
-        bool const isSameSpecies, amrex::Real const CoulombLog,
-        int const ndt );
-
-};
-
-#endif // WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_
diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package
index ffaa5fe92..9e6f5b7b6 100644
--- a/Source/Particles/Collision/Make.package
+++ b/Source/Particles/Collision/Make.package
@@ -1,3 +1,5 @@
-CEXE_sources += CollisionType.cpp
+CEXE_sources += CollisionHandler.cpp
+CEXE_sources += CollisionBase.cpp
+CEXE_sources += PairWiseCoulombCollision.cpp
 
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/Collision
diff --git a/Source/Particles/Collision/PairWiseCoulombCollision.H b/Source/Particles/Collision/PairWiseCoulombCollision.H
new file mode 100644
index 000000000..eea91d7f9
--- /dev/null
+++ b/Source/Particles/Collision/PairWiseCoulombCollision.H
@@ -0,0 +1,48 @@
+/* Copyright 2020 David Grote
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_PAIRWISECOULOMBCOLLISION_H_
+#define WARPX_PARTICLES_COLLISION_PAIRWISECOULOMBCOLLISION_H_
+
+#include "Particles/MultiParticleContainer.H"
+#include "CollisionBase.H"
+#include <AMReX_ParmParse.H>
+
+class PairWiseCoulombCollision
+    : public CollisionBase
+{
+public:
+
+    PairWiseCoulombCollision (std::string collision_name);
+
+    /** Perform the collisions
+     *
+     * @param cur_time Current time
+     * @param mypc Container of species involved
+     *
+     */
+    void doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc) override;
+
+    /** Perform all binary collisions within a tile
+     *
+     * @param lev AMR level of the tile
+     * @param mfi iterator for multifab
+     * @param species1/2 pointer to species container
+     *
+     */
+    void doCoulombCollisionsWithinTile (
+        int const lev, amrex::MFIter const& mfi,
+        WarpXParticleContainer& species1,
+        WarpXParticleContainer& species2);
+
+private:
+
+    amrex::Real m_CoulombLog;
+    bool m_isSameSpecies;
+
+};
+
+#endif // WARPX_PARTICLES_COLLISION_PAIRWISECOULOMBCOLLISION_H_
diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/PairWiseCoulombCollision.cpp
index 9ed31b483..e27b51603 100644
--- a/Source/Particles/Collision/CollisionType.cpp
+++ b/Source/Particles/Collision/PairWiseCoulombCollision.cpp
@@ -1,56 +1,71 @@
-/* Copyright 2019 Yinjian Zhao
+/* Copyright 2020 Yinjian Zhao, David Grote
  *
  * This file is part of WarpX.
  *
  * License: BSD-3-Clause-LBNL
  */
-#include "CollisionType.H"
+#include "PairWiseCoulombCollision.H"
 #include "ShuffleFisherYates.H"
 #include "ElasticCollisionPerez.H"
 #include "Utils/ParticleUtils.H"
 #include "Utils/WarpXUtil.H"
 #include "WarpX.H"
 
-CollisionType::CollisionType(
-    const std::vector<std::string>& species_names,
-    std::string const collision_name)
+using namespace amrex::literals;
+
+PairWiseCoulombCollision::PairWiseCoulombCollision (std::string const collision_name)
+    : CollisionBase(collision_name)
 {
 
-    // read collision species
-    std::vector<std::string> collision_species;
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_species_names.size() == 2,
+                                     "Pair wise Coulomb must have exactly two species.");
+
     amrex::ParmParse pp(collision_name);
-    pp.getarr("species", collision_species);
-    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(collision_species.size() == 2,
-    "Collision species must name exactly two species.");
 
     // default Coulomb log, if < 0, will be computed automatically
-    m_CoulombLog = -1.0;
+    m_CoulombLog = -1.0_rt;
     queryWithParser(pp, "CoulombLog", m_CoulombLog);
 
-    // number of time steps between collisions
-    m_ndt = 1;
-    pp.query("ndt", m_ndt);
-
-    for (int i=0; i<static_cast<int>(species_names.size()); i++)
-    {
-        if (species_names[i] == collision_species[0])
-        { m_species1_index = i; }
-        if (species_names[i] == collision_species[1])
-        { m_species2_index = i; }
-    }
-
-    if (collision_species[0] == collision_species[1])
+    if (m_species_names[0] == m_species_names[1])
         m_isSameSpecies = true;
     else
         m_isSameSpecies = false;
 
 }
 
-using namespace amrex;
+void
+PairWiseCoulombCollision::doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc)
+{
+
+    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]);
+    auto& species2 = mypc->GetParticleContainerFromName(m_species_names[1]);
+
+    // Enable tiling
+    amrex::MFItInfo info;
+    if (amrex::Gpu::notInLaunchRegion()) info.EnableTiling(species1.tile_size);
+
+    // Loop over refinement levels
+    for (int lev = 0; lev <= species1.finestLevel(); ++lev){
+
+        // Loop over all grids/tiles at this level
+#ifdef _OPENMP
+        info.SetDynamic(true);
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (amrex::MFIter mfi = species1.MakeMFIter(lev, info); mfi.isValid(); ++mfi){
+            doCoulombCollisionsWithinTile( lev, mfi, species1, species2 );
+        }
+    }
+}
+
+
 // Define shortcuts for frequently-used type names
 using ParticleType = WarpXParticleContainer::ParticleType;
 using ParticleTileType = WarpXParticleContainer::ParticleTileType;
-using ParticleBins = DenseBins<ParticleType>;
+using ParticleBins = amrex::DenseBins<ParticleType>;
 using index_type = ParticleBins::index_type;
 
 using namespace ParticleUtils;
@@ -60,22 +75,22 @@ using namespace ParticleUtils;
  * @param lev AMR level of the tile
  * @param mfi iterator for multifab
  * @param species1/2 pointer to species container
- * @param isSameSpecies true if collision is between same species
- * @param CoulombLog user input Coulomb logrithm
  * @param ndt user input number of time stpes between collisions
  *
  */
-void CollisionType::doCoulombCollisionsWithinTile
-    ( int const lev, MFIter const& mfi,
-    std::unique_ptr<WarpXParticleContainer>& species_1,
-    std::unique_ptr<WarpXParticleContainer>& species_2,
-    bool const isSameSpecies, Real const CoulombLog, int const ndt )
+void PairWiseCoulombCollision::doCoulombCollisionsWithinTile
+    ( int const lev, amrex::MFIter const& mfi,
+    WarpXParticleContainer& species_1,
+    WarpXParticleContainer& species_2)
 {
 
-    if ( isSameSpecies ) // species_1 == species_2
+    int const ndt = m_ndt;
+    amrex::Real CoulombLog = m_CoulombLog;
+
+    if ( m_isSameSpecies ) // species_1 == species_2
     {
         // Extract particles in the tile that `mfi` points to
-        ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi);
+        ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
 
         // Find the particles that are in each cell of this tile
         ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
@@ -86,22 +101,22 @@ void CollisionType::doCoulombCollisionsWithinTile
         int const n_cells = bins_1.numBins();
         // - Species 1
         auto& soa_1 = ptile_1.GetStructOfArrays();
-        ParticleReal * const AMREX_RESTRICT ux_1 =
+        amrex::ParticleReal * const AMREX_RESTRICT ux_1 =
             soa_1.GetRealData(PIdx::ux).data();
-        ParticleReal * const AMREX_RESTRICT uy_1 =
+        amrex::ParticleReal * const AMREX_RESTRICT uy_1 =
             soa_1.GetRealData(PIdx::uy).data();
-        ParticleReal * const AMREX_RESTRICT uz_1  =
+        amrex::ParticleReal * const AMREX_RESTRICT uz_1  =
             soa_1.GetRealData(PIdx::uz).data();
-        ParticleReal const * const AMREX_RESTRICT w_1 =
+        amrex::ParticleReal const * const AMREX_RESTRICT w_1 =
             soa_1.GetRealData(PIdx::w).data();
         index_type* indices_1 = bins_1.permutationPtr();
         index_type const* cell_offsets_1 = bins_1.offsetsPtr();
-        Real q1 = species_1->getCharge();
-        Real m1 = species_1->getMass();
+        amrex::Real q1 = species_1.getCharge();
+        amrex::Real m1 = species_1.getMass();
 
-        const Real dt = WarpX::GetInstance().getdt(lev);
-        Geometry const& geom = WarpX::GetInstance().Geom(lev);
-        Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box
+        const amrex::Real dt = WarpX::GetInstance().getdt(lev);
+        amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
+        amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
         const auto lo = lbound(cbx);
         const auto hi = ubound(cbx);
         int nz = hi.y-lo.y+1;
@@ -144,7 +159,7 @@ void CollisionType::doCoulombCollisionsWithinTile
                         cell_half_1, cell_stop_1,
                         indices_1, indices_1,
                         ux_1, uy_1, uz_1, ux_1, uy_1, uz_1, w_1, w_1,
-                        q1, q1, m1, m1, Real(-1.0), Real(-1.0),
+                        q1, q1, m1, m1, amrex::Real(-1.0), amrex::Real(-1.0),
                         dt*ndt, CoulombLog, dV, engine );
                 }
             }
@@ -153,8 +168,8 @@ void CollisionType::doCoulombCollisionsWithinTile
     else // species_1 != species_2
     {
         // Extract particles in the tile that `mfi` points to
-        ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi);
-        ParticleTileType& ptile_2 = species_2->ParticlesAt(lev, mfi);
+        ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
+        ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi);
 
         // Find the particles that are in each cell of this tile
         ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
@@ -166,32 +181,32 @@ void CollisionType::doCoulombCollisionsWithinTile
         int const n_cells = bins_1.numBins();
         // - Species 1
         auto& soa_1 = ptile_1.GetStructOfArrays();
-        ParticleReal * const AMREX_RESTRICT ux_1 =
+        amrex::ParticleReal * const AMREX_RESTRICT ux_1 =
             soa_1.GetRealData(PIdx::ux).data();
-        ParticleReal * const AMREX_RESTRICT uy_1 =
+        amrex::ParticleReal * const AMREX_RESTRICT uy_1 =
             soa_1.GetRealData(PIdx::uy).data();
-        ParticleReal * const AMREX_RESTRICT uz_1 =
+        amrex::ParticleReal * const AMREX_RESTRICT uz_1 =
             soa_1.GetRealData(PIdx::uz).data();
-        ParticleReal const * const AMREX_RESTRICT w_1 =
+        amrex::ParticleReal const * const AMREX_RESTRICT w_1 =
             soa_1.GetRealData(PIdx::w).data();
         index_type* indices_1 = bins_1.permutationPtr();
         index_type const* cell_offsets_1 = bins_1.offsetsPtr();
-        Real q1 = species_1->getCharge();
-        Real m1 = species_1->getMass();
+        amrex::Real q1 = species_1.getCharge();
+        amrex::Real m1 = species_1.getMass();
         // - Species 2
         auto& soa_2 = ptile_2.GetStructOfArrays();
-        Real* ux_2  = soa_2.GetRealData(PIdx::ux).data();
-        Real* uy_2  = soa_2.GetRealData(PIdx::uy).data();
-        Real* uz_2  = soa_2.GetRealData(PIdx::uz).data();
-        Real* w_2   = soa_2.GetRealData(PIdx::w).data();
+        amrex::Real* ux_2  = soa_2.GetRealData(PIdx::ux).data();
+        amrex::Real* uy_2  = soa_2.GetRealData(PIdx::uy).data();
+        amrex::Real* uz_2  = soa_2.GetRealData(PIdx::uz).data();
+        amrex::Real* w_2   = soa_2.GetRealData(PIdx::w).data();
         index_type* indices_2 = bins_2.permutationPtr();
         index_type const* cell_offsets_2 = bins_2.offsetsPtr();
-        Real q2 = species_2->getCharge();
-        Real m2 = species_2->getMass();
+        amrex::Real q2 = species_2.getCharge();
+        amrex::Real m2 = species_2.getMass();
 
-        const Real dt = WarpX::GetInstance().getdt(lev);
-        Geometry const& geom = WarpX::GetInstance().Geom(lev);
-        Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box
+        const amrex::Real dt = WarpX::GetInstance().getdt(lev);
+        amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
+        amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
         const auto lo = lbound(cbx);
         const auto hi = ubound(cbx);
         int nz = hi.y-lo.y+1;
@@ -240,11 +255,11 @@ void CollisionType::doCoulombCollisionsWithinTile
                         cell_start_1, cell_stop_1, cell_start_2, cell_stop_2,
                         indices_1, indices_2,
                         ux_1, uy_1, uz_1, ux_2, uy_2, uz_2, w_1, w_2,
-                        q1, q2, m1, m2, Real(-1.0), Real(-1.0),
+                        q1, q2, m1, m2, amrex::Real(-1.0), amrex::Real(-1.0),
                         dt*ndt, CoulombLog, dV, engine );
                 }
             }
         );
-    } // end if ( isSameSpecies)
+    } // end if ( m_isSameSpecies)
 
 }
diff --git a/Source/Particles/Collision/README.md b/Source/Particles/Collision/README.md
new file mode 100644
index 000000000..9784e2f24
--- /dev/null
+++ b/Source/Particles/Collision/README.md
@@ -0,0 +1,11 @@
+How to add a new collision operator:
+
+Create a class that inherits from CollisionBase and include the method doCollisions
+that does the operation. The constructor should call the CollisionBase constructor with
+the collision name argument and read in any input parameters with a prefix of the
+collision name that are specific to the operator.
+
+Then modify CollisionHandler.cpp to include the header file of the new class.
+In its constructor, add an if block for the input collision type creating an
+instance of the new class in the allcollisions vector. Follow PairWiseCoulombCollision
+as an example.