CD_ItoKMCPhysicsImplem.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_ITOKMCPHYSICSIMPLEM_H
#define CD_ITOKMCPHYSICSIMPLEM_H
 
// Std includes
#include <algorithm>
 
// Chombo includes
#include <ParmParse.H>
 
// Our includes
#include <CD_ItoKMCPhysics.H>
#include <CD_ParticleManagement.H>
#include <CD_Random.H>
#include <CD_Units.H>
#include <CD_DataOps.H>
#include <CD_NamespaceHeader.H>
 
using namespace Physics::ItoKMC;
 
inline ItoKMCPhysics::ItoKMCPhysics() noexcept
{
  CH_TIME("ItoKMCPhysics::ItoKMCPhysics");
 
  m_className = "ItoKMCPhysics";
 
  m_kmcReactions.clear();
  m_photoReactions.clear();
 
  // Some default settings (mostly in case user forgets to call the parsing algorithms).
  m_isDefined         = false;
  m_debug             = true;
  m_hasKMCSolver      = false;
  m_downstreamSpecies = -1;
  m_maxNewParticles   = 32;
  m_maxNewPhotons     = 32;
  m_Ncrit             = 5;
  m_eps               = 2.0;
  m_NSSA              = 5;
  m_maxIter           = 10;
  m_SSAlim            = 5.0;
  m_exitTol           = 1.E-6;
  m_algorithm         = Algorithm::ExplicitEuler;
  m_particlePlacement = ParticlePlacement::Random;
}
 
inline ItoKMCPhysics::~ItoKMCPhysics() noexcept
{
  CH_TIME("ItoKMCPhysics::~ItoKMCPhysics");
}
 
inline void
ItoKMCPhysics::define() noexcept
{
  CH_TIME("ItoKMCPhysics::define");
 
  this->defineSpeciesMap();
  this->definePhotoPathways();
 
  // Safety hook -- make sure no one defines reactions using an out-of-range species index.
#ifndef NDEBUG
  for (const auto& R : m_kmcReactions) {
    const auto& lhsReactants = R.getReactants();
    const auto& rhsReactants = R.getReactiveProducts();
    const auto& rhsPhotons   = R.getNonReactiveProducts();
 
    for (const auto& r : lhsReactants) {
      CH_assert(r < m_itoSpecies.size() + m_cdrSpecies.size());
    }
    for (const auto& r : rhsReactants) {
      CH_assert(r < m_itoSpecies.size() + m_cdrSpecies.size());
    }
    for (const auto& r : rhsPhotons) {
      CH_assert(r < m_rtSpecies.size());
    }
  }
#endif
 
  m_isDefined = true;
}
 
inline void
ItoKMCPhysics::defineSpeciesMap() noexcept
{
  CH_TIME("ItoKMCPhysics::defineSpeciesMap");
 
  const int numItoSpecies = this->getNumItoSpecies();
  const int numCdrSpecies = this->getNumCdrSpecies();
 
  int species = 0;
  for (int i = 0; i < numItoSpecies; i++, species++) {
    m_speciesMap.emplace(species, std::make_pair(SpeciesType::Ito, i));
  }
 
  for (int i = 0; i < numCdrSpecies; i++, species++) {
    m_speciesMap.emplace(species, std::make_pair(SpeciesType::CDR, i));
  }
}
 
inline bool
ItoKMCPhysics::needGradient(const int a_plasmaIndex) const noexcept
{
  CH_assert(a_plasmaIndex >= 0);
  CH_assert(a_plasmaIndex < m_itoSpecies.size() + m_cdrSpecies.size());
 
  return this->needGradients();
}
 
inline void
ItoKMCPhysics::defineKMC() const noexcept
{
  CH_TIME("ItoKMCPhysics::defineKMC");
 
  CH_assert(!m_hasKMCSolver);
 
  // Deep copy of reaction rates
  m_kmcReactionsThreadLocal.resize(0);
  for (const auto& r : m_kmcReactions) {
    m_kmcReactionsThreadLocal.emplace_back(std::make_shared<const KMCReaction>(r));
  }
 
  m_kmcSolver.define(m_kmcReactionsThreadLocal);
  m_kmcSolver.setSolverParameters(m_Ncrit, m_NSSA, m_maxIter, m_eps, m_SSAlim, m_exitTol);
  m_kmcState.define(m_itoSpecies.size() + m_cdrSpecies.size(), m_rtSpecies.size());
  m_kmcStateScratch.define(m_itoSpecies.size() + m_cdrSpecies.size(), m_rtSpecies.size());
 
  m_kmcPropensityScratch.resize(m_kmcReactionsThreadLocal.size());
 
  // Compact the reactions that actually take part in the physics time step calculation. The dt probe at the
  // end of advanceKMC scales every propensity by m_reactiveDtFactors, which the user sets to zero for every
  // reaction excluded from that calculation -- so those propensities were evaluated and then multiplied away,
  // and KMCSolver::computeDt still reduced over the whole reaction set. The probe's cost therefore grew with
  // the size of the chemistry while only the included reactions could move the answer, and a chemistry that
  // includes two reactions out of forty paid for all forty, per grid cell, twice per probe.
  //
  // The compaction is exact, not an approximation. An excluded reaction contributes muIJ * 0 to every
  // reactant sum in computeDt, and a reactant whose accumulated mu is zero is skipped there by its own
  // 'mu > numeric_limits<Real>::min()' guard -- so dropping the excluded reactions and the reactants only
  // they reach leaves the returned time step bit-identical.
  //
  // A model that never sized m_reactiveDtFactors is read as "every reaction counts", matching the parser
  // default for a reaction with no 'include_dt_calc' entry.
  m_kmcReactionsDt.resize(0);
  m_reactiveDtFactorsDt.resize(0);
 
  for (size_t i = 0; i < m_kmcReactionsThreadLocal.size(); i++) {
    const Real factor = (i < m_reactiveDtFactors.size()) ? m_reactiveDtFactors[i] : 1.0;
 
    if (factor != 0.0) {
      m_kmcReactionsDt.emplace_back(m_kmcReactionsThreadLocal[i]);
      m_reactiveDtFactorsDt.emplace_back(factor);
    }
  }
 
  m_kmcPropensityScratchDt.resize(m_kmcReactionsDt.size());
 
  m_hasKMCSolver = true;
}
 
inline void
ItoKMCPhysics::killKMC() const noexcept
{
  CH_TIME("ItoKMCPhysics::killKMC");
 
  CH_assert(m_hasKMCSolver);
 
  m_kmcReactionsThreadLocal.resize(0);
  m_kmcSolver.define(m_kmcReactionsThreadLocal);
  m_kmcState.define(0, 0);
  m_kmcStateScratch.define(0, 0);
 
  m_kmcPropensityScratch.resize(0);
 
  m_kmcReactionsDt.resize(0);
  m_reactiveDtFactorsDt.resize(0);
  m_kmcPropensityScratchDt.resize(0);
 
  m_hasKMCSolver = false;
}
 
inline void
ItoKMCPhysics::definePhotoPathways() noexcept
{
  CH_TIME("ItoKMCPhysics::definePhotoPathways");
 
  // Build a temporary list of pathways. I.e. restructure the list of reactions
  //
  // Y1 -> A
  // Y1 -> B
  // Y2 -> C
  // Y2 -> D
  //
  // into separate lists for Y1, Y2, ....
  //
  std::map<int, std::vector<std::pair<int, Real>>> pathways;
 
  for (int i = 0; i < m_photoReactions.size(); i++) {
    const ItoKMCPhotoReaction& r = m_photoReactions[i];
 
    const size_t& src        = r.getSourcePhoton();
    const Real    efficiency = r.getEfficiency();
 
    pathways[src].emplace_back(std::make_pair(i, efficiency));
  }
 
  // Go through the temporary pathways list and compute the relative efficiencies of one of the
  // photons triggering a reaction. The relative efficiencies are given by
  //
  // p(i) = R(i)/sum_j R(j).
  //
  for (const auto& p : pathways) {
    const int                               photoSpecies             = p.first;
    const std::vector<std::pair<int, Real>> reactionsAndEfficiencies = p.second;
 
    std::map<int, int> localToGlobalMap;
    std::list<double>  efficiencies;
    double             sumEfficiencies = 0.0;
 
    for (int i = 0; i < reactionsAndEfficiencies.size(); i++) {
      sumEfficiencies += (double)reactionsAndEfficiencies[i].second;
    }
 
    for (int i = 0; i < reactionsAndEfficiencies.size(); i++) {
      localToGlobalMap.emplace(i, reactionsAndEfficiencies[i].first);
      efficiencies.emplace_back((double)reactionsAndEfficiencies[i].second / sumEfficiencies);
    }
 
    std::discrete_distribution<int> distribution(efficiencies.begin(), efficiencies.end());
 
    m_photoPathways.insert(std::make_pair((int)photoSpecies, std::make_pair(distribution, localToGlobalMap)));
  }
}
 
inline const std::map<int, std::pair<SpeciesType, int>>&
ItoKMCPhysics::getSpeciesMap() const noexcept
{
  CH_TIME("ItoKMCPhysics::getSpeciesMap");
 
  return m_speciesMap;
}
 
inline void
ItoKMCPhysics::parseRuntimeOptions() noexcept
{
  CH_TIME("ItoKMCPhysics::parseRuntimeOptions");
 
  this->parsePPC();
  this->parseDebug();
  this->parseAlgorithm();
}
 
inline void
ItoKMCPhysics::parsePPC() noexcept
{
  CH_TIME("ItoKMCPhysics::parsePPC");
 
  ParmParse pp(m_className.c_str());
 
  pp.get("max_new_particles", m_maxNewParticles);
  pp.get("max_new_photons", m_maxNewPhotons);
 
  if (m_maxNewParticles < 1) {
    MayDay::Abort(("ItoKMCPhysics::parsePPC - '" + m_className + ".max_new_particles' must be >= 1").c_str());
  }
}
 
inline void
ItoKMCPhysics::parseDebug() noexcept
{
  CH_TIME("ItoKMCPhysics::parseDebug");
 
  ParmParse pp(m_className.c_str());
 
  pp.get("debug", m_debug);
}
 
inline void
ItoKMCPhysics::parseAlgorithm() noexcept
{
  CH_TIME("ItoKMCPhysics::parseAlgorithm");
 
  ParmParse pp(m_className.c_str());
 
  std::string str;
 
  pp.get("algorithm", str);
  pp.get("crit_num", m_Ncrit);
  pp.get("SSA_num", m_NSSA);
  pp.get("prop_eps", m_eps);
  pp.get("SSA_lim", m_SSAlim);
  pp.get("max_iter", m_maxIter);
  pp.get("exit_tolerance", m_exitTol);
 
  if (str == "ssa") {
    m_algorithm = Algorithm::SSA;
  }
  else if (str == "explicit_euler") {
    m_algorithm = Algorithm::ExplicitEuler;
  }
  else if (str == "midpoint") {
    m_algorithm = Algorithm::Midpoint;
  }
  else if (str == "prc") {
    m_algorithm = Algorithm::PRC;
  }
  else if (str == "implicit_euler") {
    m_algorithm = Algorithm::ImplicitEuler;
  }
  else if (str == "hybrid_explicit_euler") {
    m_algorithm = Algorithm::HybridExplicitEuler;
  }
  else if (str == "hybrid_midpoint") {
    m_algorithm = Algorithm::HybridMidpoint;
  }
  else if (str == "hybrid_prc") {
    m_algorithm = Algorithm::HybridPRC;
  }
  else if (str == "hybrid_implicit_euler") {
    m_algorithm = Algorithm::HybridImplicitEuler;
  }
  else {
    MayDay::Error("ItoKMCPhysics::parseAlgorithm - unknown algorithm requested");
  }
}
 
inline const Vector<RefCountedPtr<ItoSpecies>>&
ItoKMCPhysics::getItoSpecies() const
{
  return m_itoSpecies;
}
 
inline const Vector<RefCountedPtr<CdrSpecies>>&
ItoKMCPhysics::getCdrSpecies() const
{
  return m_cdrSpecies;
}
 
inline const Vector<RefCountedPtr<RtSpecies>>&
ItoKMCPhysics::getRtSpecies() const
{
  return m_rtSpecies;
}
 
inline const Vector<DiffusionFunction>&
ItoKMCPhysics::getItoDiffusionFunctions() const noexcept
{
  return m_itoDiffusionFunctions;
}
 
inline int
ItoKMCPhysics::getNumItoSpecies() const
{
  return m_itoSpecies.size();
}
 
inline int
ItoKMCPhysics::getNumCdrSpecies() const
{
  return m_cdrSpecies.size();
}
 
inline int
ItoKMCPhysics::getNumPlasmaSpecies() const
{
  return m_itoSpecies.size() + m_cdrSpecies.size();
}
 
inline int
ItoKMCPhysics::getNumPhotonSpecies() const
{
  return m_rtSpecies.size();
}
 
inline Real
ItoKMCPhysics::initialSigma(const Real a_time, const RealVect& a_pos) const
{
  return 0.0;
}
 
inline void
ItoKMCPhysics::advanceKMC(Vector<FPR>&            a_numParticles,
                          Vector<FPR>&            a_numNewPhotons,
                          Real&                   a_physicsDt,
                          const Vector<Real>&     a_phi,
                          const Vector<RealVect>& a_gradPhi,
                          const Real              a_dt,
                          const RealVect          a_E,
                          const RealVect          a_pos,
                          const Real              a_dx,
                          const Real              a_kappa) const
{
  // Note: This is called PER GRID CELL, i.e. within OpenMP parallel regions. For this reason the KMC solver
  //       must be defined through defineKMC() (which must be later killed).
  CH_assert(m_isDefined);
  CH_assert(m_hasKMCSolver);
 
  std::vector<FPR>& kmcParticles = m_kmcState.getReactiveState();
  std::vector<FPR>& kmcPhotons   = m_kmcState.getNonReactiveState();
 
  for (size_t i = 0; i < a_numParticles.size(); i++) {
    kmcParticles[i] = a_numParticles[i];
  }
 
  for (auto& p : kmcPhotons) {
    p = 0LL;
  }
 
  // Lambda function used for computing charge before and after reactions. Used only in debug mode
  // for ensuring that nothing goes wrong with charge conservation in the chemistry integration.
  auto computeCharge = [&]() -> long long {
    long long Q = 0.0;
    for (int i = 0; i < kmcParticles.size(); i++) {
      const SpeciesType& speciesType  = m_speciesMap.at(i).first;
      const int&         speciesIndex = m_speciesMap.at(i).second;
 
      int Z = 0;
 
      switch (speciesType) {
      case SpeciesType::Ito: {
        Z = m_itoSpecies[speciesIndex]->getChargeNumber();
 
        break;
      }
      case SpeciesType::CDR: {
        Z = m_cdrSpecies[speciesIndex]->getChargeNumber();
 
        break;
      }
      default: {
        MayDay::Abort("ItoKMCPhysics::advanceKMC -- logic bust in computeCharge()");
 
        break;
      }
      }
 
      Q += llround(kmcParticles[i]) * Z;
    }
 
    return Q;
  };
 
  // In debug mode, compute the total charge.
  const long long chargeBefore = m_debug ? computeCharge() : 0LL;
 
  // Update the reaction rates to be used by the KMC solver.
  this->updateReactionRates(m_kmcReactionsThreadLocal, a_E, a_pos, a_phi, a_gradPhi, a_dt, a_dx, a_kappa);
 
  // Run the KMC solver.
  switch (m_algorithm) {
  case Algorithm::SSA: {
    m_kmcSolver.advanceSSA(m_kmcState, a_dt);
 
    break;
  }
  case Algorithm::ExplicitEuler: {
    m_kmcSolver.advanceTau(m_kmcState, a_dt, KMCLeapPropagator::ExplicitEuler);
 
    break;
  }
  case Algorithm::Midpoint: {
    m_kmcSolver.advanceTau(m_kmcState, a_dt, KMCLeapPropagator::Midpoint);
 
    break;
  }
  case Algorithm::PRC: {
    m_kmcSolver.advanceTau(m_kmcState, a_dt, KMCLeapPropagator::PRC);
 
    break;
  }
  case Algorithm::ImplicitEuler: {
    m_kmcSolver.advanceTau(m_kmcState, a_dt, KMCLeapPropagator::ImplicitEuler);
 
    break;
  }
  case Algorithm::HybridExplicitEuler: {
    m_kmcSolver.advanceHybrid(m_kmcState, a_dt, KMCLeapPropagator::ExplicitEuler);
 
    break;
  }
  case Algorithm::HybridMidpoint: {
    m_kmcSolver.advanceHybrid(m_kmcState, a_dt, KMCLeapPropagator::Midpoint);
 
    break;
  }
  case Algorithm::HybridPRC: {
    m_kmcSolver.advanceHybrid(m_kmcState, a_dt, KMCLeapPropagator::PRC);
 
    break;
  }
  case Algorithm::HybridImplicitEuler: {
    m_kmcSolver.advanceHybrid(m_kmcState, a_dt, KMCLeapPropagator::ImplicitEuler);
 
    break;
  }
  default: {
    MayDay::Error("ItoKMCPhysics::advanceKMC - logic bust");
  }
  }
 
  // Put KMC back into ItoKMC
  for (size_t i = 0; i < a_numParticles.size(); i++) {
    a_numParticles[i] = (FPR)kmcParticles[i];
  }
  for (size_t i = 0; i < a_numNewPhotons.size(); i++) {
    a_numNewPhotons[i] = (FPR)kmcPhotons[i];
  }
 
  const long long chargeAfter = m_debug ? computeCharge() : 0LL;
 
  if (chargeAfter != chargeBefore) {
    MayDay::Warning("ItoKMCPhysics::advanceKMC -- charge not conserved!");
  }
 
  // This loop is for isolating reactions that the user will explicitly ask to fire before computing the physics-based
  // time step. It exists because if there are no electrons but lots of ions, one may get a time step that is too large
  // because X/|sum mu| is zero for the electrons. Similarly, one may get a time step that is too small away from
  // ionization regions because X/|sum mu| may be tiny if there is, say, 1 electron but lots of detachment.
  //
  // Everything below runs on the compacted reaction set built by defineKMC(), not on the full one: only the
  // reactions the user included can move the answer, and the excluded ones cost a propensity evaluation and a
  // reduction apiece, per probe, per grid cell. With nothing included there is nothing to limit, so the whole
  // block drops out.
  if (m_kmcReactionsDt.empty()) {
    return;
  }
 
  // Fills the reusable propensity buffer rather than returning one, so that the probes below do not
  // allocate once per reaction per grid cell.
  auto fillPropensitiesDt = [&](const KMCState& a_state) -> void {
    for (size_t i = 0; i < m_kmcReactionsDt.size(); i++) {
      m_kmcPropensityScratchDt[i] = m_kmcReactionsDt[i]->propensity(a_state) * m_reactiveDtFactorsDt[i];
    }
  };
 
  // Do a time step limitation on the complete state, using the user-specified reactions.
  fillPropensitiesDt(m_kmcState);
 
  a_physicsDt = std::min(a_physicsDt,
                         m_kmcSolver.computeDt(m_kmcState, m_kmcReactionsDt, m_kmcPropensityScratchDt, 1.0));
 
  // Go through reactions that are potentially detaching species and do the calculation again. This probe walks
  // the FULL reaction set even though the limit is evaluated on the compacted one: an excluded reaction can
  // still deplete a species that an included reaction depends on, and that depletion is exactly what the probe
  // is looking for.
  for (const auto& reaction : m_kmcReactionsThreadLocal) {
    const auto N = reaction->computeCriticalNumberOfReactions(m_kmcState);
 
    // Trigger on both the reactions and the reactions with completely consumed reactants. The state
    // is only copied inside the branch because that is the only place it is read -- copying it above
    // the test meant paying for a copy per reaction per grid cell, including for the reactions that
    // never enter here.
    if (N < std::numeric_limits<FPR>::max()) {
      m_kmcStateScratch = m_kmcState;
 
      reaction->advanceState(m_kmcStateScratch, N);
 
      fillPropensitiesDt(m_kmcStateScratch);
 
      a_physicsDt = std::min(a_physicsDt,
                             m_kmcSolver.computeDt(m_kmcStateScratch, m_kmcReactionsDt, m_kmcPropensityScratchDt, 1.0));
    }
  }
}
 
inline void
ItoKMCPhysics::reconcileParticles(Vector<ParticleSoA<ItoParticle>*>& a_particles,
                                  const Vector<FPR>&                 a_newNumParticles,
                                  const Vector<FPR>&                 a_oldNumParticles,
                                  const RealVect                     a_electricField,
                                  const RealVect                     a_cellPos,
                                  const RealVect                     a_centroidPos,
                                  const RealVect                     a_lo,
                                  const RealVect                     a_hi,
                                  const RealVect                     a_bndryCentroid,
                                  const RealVect                     a_bndryNormal,
                                  const Real                         a_dx,
                                  const Real                         a_kappa) const noexcept
{
  CH_assert(m_isDefined);
  CH_assert(a_particles.size() == a_newNumParticles.size());
  CH_assert(a_oldNumParticles.size() == a_newNumParticles.size());
 
  if (m_debug) {
    for (int i = 0; i < a_particles.size(); i++) {
      const FPR& numNew = a_newNumParticles[i];
      const FPR& numOld = a_oldNumParticles[i];
 
      if (numNew < (FPR)0) {
        MayDay::Warning("ItoKMCPhysics::reconcileParticles - new number of particles is < 0 (overflow issue?)");
      }
      else if (static_cast<long long>(numNew) < 0LL) {
        MayDay::Warning("ItoKMCPhysics::reconcileParticles - integer overflow!");
      }
 
      if (numOld < 0) {
        MayDay::Warning("ItoKMCPhysics::reconcileParticles - old number of particles is < 0");
      }
      else if (static_cast<long long>(numOld) < 0LL) {
        MayDay::Warning("ItoKMCPhysics::reconcileParticles - integer overflow for old particles!");
      }
    }
  }
 
  // Compute the upstream position of the particles (which is usually the electrons).
  bool     hasDownstream    = false;
  RealVect upstreamPosition = RealVect::Zero;
  RealVect upstreamLo       = -0.5 * RealVect::Unit;
  RealVect upstreamHi       = +0.5 * RealVect::Unit;
  RealVect v                = RealVect::Zero;
  int      Z                = 0;
 
  if (m_particlePlacement == ParticlePlacement::Downstream) {
    CH_assert(m_downstreamSpecies >= 0);
 
    Z = m_itoSpecies[m_downstreamSpecies]->getChargeNumber();
 
    if (Z != 0) {
      v = Z * a_electricField;
      v = v / v.vectorLength();
    }
 
    hasDownstream = this->computeUpstreamPosition(upstreamPosition,
                                                  upstreamLo,
                                                  upstreamHi,
                                                  Z,
                                                  *a_particles[m_downstreamSpecies],
                                                  a_electricField,
                                                  a_cellPos,
                                                  a_dx);
  }
 
  for (int i = 0; i < a_particles.size(); i++) {
    const long long diff = static_cast<long long>(a_newNumParticles[i] - a_oldNumParticles[i]);
 
    if (diff > 0LL) {
      const long long numParticles = static_cast<long long>(a_particles[i]->size());
 
      // Choose weights for the new particles and go with one of the placement algorithms.
      ParticleManagement::partitionParticleWeights(m_weightScratch, diff, static_cast<long long>(m_maxNewParticles));
 
      const std::vector<long long>& particleWeights = m_weightScratch;
 
      if (m_particlePlacement == ParticlePlacement::Parent) {
        this->buildParentWeights(m_parentWeightScratch, *a_particles[i], static_cast<std::size_t>(numParticles));
      }
 
      for (const auto& w : particleWeights) {
        RealVect parentPos = RealVect::Zero;
 
        const bool hasParent = (m_particlePlacement == ParticlePlacement::Parent) &&
                               this->sampleParentPosition(parentPos, *a_particles[i], m_parentWeightScratch);
 
        const RealVect x = this->drawNewParticlePosition(hasParent,
                                                         hasDownstream,
                                                         parentPos,
                                                         upstreamPosition,
                                                         upstreamLo,
                                                         upstreamHi,
                                                         v,
                                                         a_cellPos,
                                                         a_centroidPos,
                                                         a_lo,
                                                         a_hi,
                                                         a_bndryCentroid,
                                                         a_bndryNormal,
                                                         a_dx,
                                                         a_kappa);
 
        a_particles[i]->append(x, 1.0 * w, ItoParticle{});
      }
    }
    else if (diff < 0LL) {
      // Removing particles is a bit more difficult because we need to manipulate weights.
      this->removeParticles(*a_particles[i], -diff);
    }
  }
}
 
inline RealVect
ItoKMCPhysics::drawNewParticlePosition(const bool     a_hasParent,
                                       const bool     a_hasDownstream,
                                       const RealVect a_parentPos,
                                       const RealVect a_upstreamPos,
                                       const RealVect a_upstreamLo,
                                       const RealVect a_upstreamHi,
                                       const RealVect a_downstreamDirection,
                                       const RealVect a_cellPos,
                                       const RealVect a_centroidPos,
                                       const RealVect a_lo,
                                       const RealVect a_hi,
                                       const RealVect a_bndryCentroid,
                                       const RealVect a_bndryNormal,
                                       const Real     a_dx,
                                       const Real     a_kappa) const noexcept
{
  RealVect x = RealVect::Zero;
 
  switch (m_particlePlacement) {
  case ParticlePlacement::Centroid: {
    x = a_cellPos + a_centroidPos * a_dx;
 
    break;
  }
  case ParticlePlacement::Random: {
    x = Random::randomPosition(a_cellPos, a_lo, a_hi, a_bndryCentroid, a_bndryNormal, a_dx, a_kappa);
 
    break;
  }
  case ParticlePlacement::Parent: {
    // Ionization happens at the electrons, so put the new particle on top of one of the particles that were
    // already there, drawn with probability proportional to the parent weight. This is the placement that
    // introduces no sub-grid transport of its own: it neither scatters the new weight across the cell (which
    // acts as a numerical diffusion) nor pulls it towards the cell centre (which acts as a drag on the front).
    // The inherited position is by construction on the fluid side of the EB, so no cut-cell test is needed.
    //
    // Falls back to a random position when there is nothing to inherit from, e.g. photoionization products
    // landing in a cell that holds no particles of this species yet.
    if (a_hasParent) {
      x = a_parentPos;
    }
    else {
      x = Random::randomPosition(a_cellPos, a_lo, a_hi, a_bndryCentroid, a_bndryNormal, a_dx, a_kappa);
    }
 
    break;
  }
  case ParticlePlacement::Downstream: {
    if (a_hasDownstream) {
      x = Random::randomPosition(a_upstreamLo, a_upstreamHi, a_upstreamPos, a_downstreamDirection);
 
      if ((x - a_bndryCentroid).dotProduct(a_bndryNormal) < 0.0) {
        x = a_upstreamPos;
      }
 
      x = a_cellPos + x * a_dx;
    }
    else {
 
      x = Random::randomPosition(a_cellPos, a_lo, a_hi, a_bndryCentroid, a_bndryNormal, a_dx, a_kappa);
    }
 
    break;
  }
  default: {
    MayDay::Error("ItoKMCPhysics::drawNewParticlePosition - logic bust");
 
    break;
  }
  }
 
  return x;
}
 
inline void
ItoKMCPhysics::buildParentWeights(std::vector<Real>&              a_cumulativeWeights,
                                  const ParticleSoA<ItoParticle>& a_particles,
                                  const std::size_t               a_numEligible) const noexcept
{
  a_cumulativeWeights.clear();
  a_cumulativeWeights.reserve(a_numEligible);
 
  Real runningSum = 0.0;
 
  for (std::size_t j = 0; j < a_numEligible; j++) {
    runningSum += a_particles.weight(j);
 
    a_cumulativeWeights.push_back(runningSum);
  }
}
 
inline bool
ItoKMCPhysics::sampleParentPosition(RealVect&                       a_pos,
                                    const ParticleSoA<ItoParticle>& a_particles,
                                    const std::vector<Real>&        a_cumulativeWeights) const noexcept
{
  if (a_cumulativeWeights.empty() || a_cumulativeWeights.back() <= 0.0) {
    return false;
  }
 
  const Real r = Random::getUniformReal01() * a_cumulativeWeights.back();
 
  // First entry whose running sum reaches r. The end() case is round-off leaving the last sum a hair below r.
  const auto it = std::lower_bound(a_cumulativeWeights.cbegin(), a_cumulativeWeights.cend(), r);
 
  const std::size_t parent = (it == a_cumulativeWeights.cend())
                               ? a_cumulativeWeights.size() - 1
                               : static_cast<std::size_t>(it - a_cumulativeWeights.cbegin());
 
  a_pos = a_particles.position(parent);
 
  return true;
}
 
inline void
ItoKMCPhysics::reconcileCdrParticles(Vector<ParticleSoA<NoPayload>*>&         a_cdrParticles,
                                     const Vector<ParticleSoA<ItoParticle>*>& a_itoParticles,
                                     const Vector<long long>&                 a_numNewParticles,
                                     const RealVect                           a_electricField,
                                     const RealVect                           a_cellPos,
                                     const RealVect                           a_centroidPos,
                                     const RealVect                           a_lo,
                                     const RealVect                           a_hi,
                                     const RealVect                           a_bndryCentroid,
                                     const RealVect                           a_bndryNormal,
                                     const Real                               a_dx,
                                     const Real                               a_kappa) const noexcept
{
  CH_assert(m_isDefined);
  CH_assert(a_cdrParticles.size() == a_numNewParticles.size());
  CH_assert(a_itoParticles.size() == m_itoSpecies.size());
 
  // Most cells produce nothing for most species, and everything below is per-cell setup that such a cell must not
  // pay for. Note that the productions are >= 0 by contract -- the caller has already split off the removal.
  bool hasProduction = false;
 
  for (int i = 0; i < a_numNewParticles.size(); i++) {
    CH_assert(a_numNewParticles[i] >= 0LL);
 
    hasProduction = hasProduction || (a_numNewParticles[i] > 0LL);
  }
 
  if (!hasProduction) {
    return;
  }
 
  // Downstream placement needs the drift direction and the downstream region, exactly as in reconcileParticles.
  bool     hasDownstream    = false;
  RealVect upstreamPosition = RealVect::Zero;
  RealVect upstreamLo       = -0.5 * RealVect::Unit;
  RealVect upstreamHi       = +0.5 * RealVect::Unit;
  RealVect v                = RealVect::Zero;
 
  if (m_particlePlacement == ParticlePlacement::Downstream) {
    CH_assert(m_downstreamSpecies >= 0);
 
    const int Z = m_itoSpecies[m_downstreamSpecies]->getChargeNumber();
 
    if (Z != 0) {
      v = Z * a_electricField;
      v = v / v.vectorLength();
    }
 
    hasDownstream = this->computeUpstreamPosition(upstreamPosition,
                                                  upstreamLo,
                                                  upstreamHi,
                                                  Z,
                                                  *a_itoParticles[m_downstreamSpecies],
                                                  a_electricField,
                                                  a_cellPos,
                                                  a_dx);
  }
 
  // ParticlePlacement::Parent has no meaning here: a CDR product has no parents in the cell, its container holding
  // only what the current step created. drawNewParticlePosition falls back to a uniformly random position.
  constexpr bool hasParent = false;
  const RealVect parentPos = RealVect::Zero;
 
  for (int i = 0; i < a_cdrParticles.size(); i++) {
    if (a_numNewParticles[i] > 0LL) {
 
      // Same partition as the Ito products: at most m_maxNewParticles computational particles, each carrying an
      // integer number of physical ones.
      ParticleManagement::partitionParticleWeights(m_weightScratch,
                                                   a_numNewParticles[i],
                                                   static_cast<long long>(m_maxNewParticles));
 
      for (const auto& w : m_weightScratch) {
        const RealVect x = this->drawNewParticlePosition(hasParent,
                                                         hasDownstream,
                                                         parentPos,
                                                         upstreamPosition,
                                                         upstreamLo,
                                                         upstreamHi,
                                                         v,
                                                         a_cellPos,
                                                         a_centroidPos,
                                                         a_lo,
                                                         a_hi,
                                                         a_bndryCentroid,
                                                         a_bndryNormal,
                                                         a_dx,
                                                         a_kappa);
 
        a_cdrParticles[i]->append(x, 1.0 * w);
      }
    }
  }
}
 
inline bool
ItoKMCPhysics::computeUpstreamPosition(RealVect&                       a_pos,
                                       RealVect&                       a_lo,
                                       RealVect&                       a_hi,
                                       const int&                      a_Z,
                                       const ParticleSoA<ItoParticle>& a_particles,
                                       const RealVect&                 a_electricField,
                                       const RealVect&                 a_cellPos,
                                       const Real&                     a_dx) const noexcept
{
  CH_assert(a_dx > 0.0);
 
  a_pos = RealVect::Zero;
  a_lo  = -0.5 * RealVect::Unit;
  a_hi  = +0.5 * RealVect::Unit;
 
  const int      Z = (a_Z > 0) ? 1 : (a_Z < 0) ? -1 : 0;
  const RealVect E = a_electricField / a_electricField.vectorLength();
  const RealVect v = Z * E;
 
  // Upstream can only exist if there is a velocity direction and there are particles.
  const bool hasDownstream = (a_particles.size() > 0) && (v.vectorLength() > 0.0);
 
  if (hasDownstream) {
 
    Real D = std::numeric_limits<Real>::max();
 
    for (std::size_t i = 0; i < a_particles.size(); i++) {
      const RealVect x = (a_particles.position(i) - a_cellPos) / a_dx;
      const Real     d = x.dotProduct(v);
 
      if (d < D) {
        D     = d;
        a_pos = x;
      }
    }
 
    DataOps::computeMinValidBox(a_lo, a_hi, v, a_pos);
  }
 
  if (m_debug) {
    for (int dir = 0; dir < SpaceDim; dir++) {
      if (a_pos[dir] > 0.5 || a_pos[dir] < -0.5) {
        MayDay::Abort("ItoKMCPhysics::computeUpstreamPosition - logic bust");
      }
    }
  }
 
  return hasDownstream;
}
 
inline void
ItoKMCPhysics::removeParticles(ParticleSoA<ItoParticle>& a_particles, const long long a_numParticlesToRemove) const
{
  constexpr long long zero = 0LL;
 
  CH_assert(m_isDefined);
  CH_assert(a_numParticlesToRemove >= zero);
 
  // Quick lambda for getting total particle weight. Used for debugging.
  auto getTotalWeight = [&]() -> long long {
    long long W = zero;
 
    for (std::size_t i = 0; i < a_particles.size(); i++) {
      W += llround(a_particles.weight(i));
 
      if (a_particles.weight(i) < 1.0) {
        MayDay::Error("ItoKMCPhysics::removeParticles -- bad particle mass!");
      }
    }
 
    return W;
  };
 
  if (a_numParticlesToRemove > zero) {
 
    // For debugging only.
    long long totalWeightBefore = 0;
    long long totalWeightAfter  = 0;
 
    // Debug hook, compute the total particle weight before we start removing weights.
    if (m_debug) {
      totalWeightBefore = getTotalWeight();
 
      if (totalWeightBefore < a_numParticlesToRemove) {
        MayDay::Error("ItoKMCPhysics::removeParticles: logic bust (trying to remove too many particles)");
      }
    }
 
    // Remove physical particles.
    ParticleManagement::removePhysicalParticles(a_particles, a_numParticlesToRemove);
 
    // Remove particles with too low weight.
    ParticleManagement::deleteParticles(a_particles, std::numeric_limits<Real>::min());
 
    // Debug hook, make sure that particle weights are > 0 AND we've removed the desired
    // particle weight.
    if (m_debug) {
      totalWeightAfter = getTotalWeight();
 
      const long long errDiff = std::abs(totalWeightBefore - totalWeightAfter) - a_numParticlesToRemove;
      if (std::abs(errDiff) != zero) {
 
        pout() << "ItoKMCPhysics::removeParticles: Total weight before = " << totalWeightBefore << endl;
        pout() << "ItoKMCPhysics::removeParticles: Total weight after  = " << totalWeightAfter << endl;
        pout() << "ItoKMCPhysics::removeParticles: Should have removed = " << a_numParticlesToRemove << endl;
        pout() << "ItoKMCPhysics::removeParticles: Error               = " << errDiff << endl;
 
        MayDay::Abort("ItoKMCPhysics::removeParticles - incorrect mass removed");
      }
    }
  }
}
 
inline void
ItoKMCPhysics::reconcilePhotons(Vector<ParticleSoA<Photon>*>& a_newPhotons,
                                const Vector<FPR>&            a_numNewPhotons,
                                const RealVect                a_cellPos,
                                const RealVect                a_centroidPos,
                                const RealVect                a_lo,
                                const RealVect                a_hi,
                                const RealVect                a_bndryCentroid,
                                const RealVect                a_bndryNormal,
                                const Real                    a_dx,
                                const Real                    a_kappa) const noexcept
{
  CH_assert(m_isDefined);
 
  for (int i = 0; i < a_newPhotons.size(); i++) {
    if (a_numNewPhotons[i] > 0LL) {
 
      ParticleManagement::partitionParticleWeights(m_weightScratch,
                                                   static_cast<long long>(a_numNewPhotons[i]),
                                                   static_cast<long long>(m_maxNewPhotons));
 
      const std::vector<long long>& photonWeights = m_weightScratch;
 
      for (const auto& w : photonWeights) {
        const RealVect x = Random::randomPosition(a_cellPos, a_lo, a_hi, a_bndryCentroid, a_bndryNormal, a_dx, a_kappa);
        const RealVect v = Units::c * Random::getDirection();
        const Real     kappa = m_rtSpecies[i]->getAbsorptionCoefficient(x);
 
        a_newPhotons[i]->append(
          x,
          1.0 * w,
          Photon{
            static_cast<ParticleReal>(kappa),
            D_DECL(static_cast<ParticleReal>(v[0]), static_cast<ParticleReal>(v[1]), static_cast<ParticleReal>(v[2]))});
      }
    }
  }
}
 
inline void
ItoKMCPhysics::reconcilePhotoionization(Vector<ParticleSoA<ItoParticle>*>&  a_itoParticles,
                                        Vector<ParticleSoA<NoPayload>*>&    a_cdrParticles,
                                        const Vector<ParticleSoA<Photon>*>& a_absorbedPhotons) const noexcept
{
  CH_assert(m_isDefined);
  CH_assert(a_itoParticles.size() == m_itoSpecies.size());
  CH_assert(a_cdrParticles.size() == m_cdrSpecies.size());
 
  for (int i = 0; i < a_absorbedPhotons.size(); i++) {
    if (m_photoPathways.find(i) != m_photoPathways.end()) {
      std::discrete_distribution<int> d                = m_photoPathways.at(i).first;
      const std::map<int, int>&       localToGlobalMap = m_photoPathways.at(i).second;
 
      const ParticleSoA<Photon>& absorbedPhotons = *a_absorbedPhotons[i];
 
      for (std::size_t p = 0; p < absorbedPhotons.size(); p++) {
        const RealVect x = absorbedPhotons.position(p);
        const Real     w = absorbedPhotons.weight(p);
 
        // Determine the photo-reaction type.
        const int localReaction  = Random::get(d);
        const int globalReaction = localToGlobalMap.at(localReaction);
 
        const ItoKMCPhotoReaction& photoReaction = m_photoReactions[globalReaction];
        const std::list<size_t>&   plasmaTargets = photoReaction.getTargetSpecies();
 
        for (const auto& t : plasmaTargets) {
          const SpeciesType& type       = m_speciesMap.at(t).first;
          const int&         localIndex = m_speciesMap.at(t).second;
 
          if (type == SpeciesType::Ito) {
            a_itoParticles[localIndex]->append(x, w, ItoParticle{});
          }
          else if (type == SpeciesType::CDR) {
            a_cdrParticles[localIndex]->append(x, w);
          }
          else {
            MayDay::Error("CD_ItoKMCPhysics.H - logic bust in reconcilePhotoionization");
          }
        }
      }
    }
  }
}
 
inline RealVect
ItoKMCPhysics::noDiffusion(const ItoParticle& a_particle, const Real a_dt) const noexcept
{
  return RealVect::Zero;
}
 
inline RealVect
ItoKMCPhysics::isotropicDiffusion(const ItoParticle& a_particle, const Real a_dt) const noexcept
{
  RealVect r = RealVect::Zero;
 
  for (int dir = 0; dir < SpaceDim; dir++) {
    r[dir] = Random::getNormal01();
  }
 
  return sqrt(2.0 * a_particle.diffusion * a_dt) * r;
}
 
inline RealVect
ItoKMCPhysics::forwardIsotropicDiffusion(const ItoParticle& a_particle, const Real a_dt) const noexcept
{
  RealVect hop = this->isotropicDiffusion(a_particle, a_dt);
 
  const RealVect v = RealVect(D_DECL(a_particle.vx, a_particle.vy, a_particle.vz));
 
  if (v != RealVect::Zero) {
    const RealVect u = v / v.vectorLength();
    const Real     d = u.dotProduct(hop);
 
    hop -= std::min(d, 0.0) * u;
  }
 
  return hop;
}
 
#include <CD_NamespaceFooter.H>
 
#endif