CD_ItoKMCPhysics.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_ITOKMCPHYSICS_H
#define CD_ITOKMCPHYSICS_H
 
// Std includes
#include <memory>
#include <vector>
 
// Chombo includes
#include <RealVect.H>
#include <RefCountedPtr.H>
#include <List.H>
 
// Our includes
#include <CD_ItoSpecies.H>
#include <CD_CdrSpecies.H>
#include <CD_RtSpecies.H>
#include <CD_Photon.H>
#include <CD_ItoParticle.H>
#include <CD_ParticleSoA.H>
#include <CD_ItoKMCPhotoReaction.H>
#include <CD_ItoKMCSurfaceReactionSet.H>
#include <CD_KMCDualState.H>
#include <CD_KMCDualStateReaction.H>
#include <CD_KMCSolver.H>
#include <CD_NamespaceHeader.H>
 
namespace Physics {
namespace ItoKMC {
 
using FPR = Real;
 
using KMCState = KMCDualState<FPR>;
 
using KMCReaction = KMCDualStateReaction<KMCState, FPR>;
 
using KMCSolverType = KMCSolver<KMCReaction, KMCState, FPR>;
 
using DiffusionFunction = std::function<RealVect(const ItoParticle& p, const Real& a_dt)>;
 
enum SpeciesType
{
  Ito, 
  CDR  
};
 
class ItoKMCPhysics
{
public:
  inline ItoKMCPhysics() noexcept;
 
  inline virtual ~ItoKMCPhysics() noexcept;
 
  inline void
  defineKMC() const noexcept;
 
  inline void
  killKMC() const noexcept;
 
  virtual Real
  getNeutralDensity(const RealVect& a_pos) const noexcept = 0;
 
  virtual Real
  computeAlpha(const Real a_E, const RealVect& a_x) const = 0;
 
  virtual Real
  computeEta(const Real a_E, const RealVect& a_x) const = 0;
 
  const Vector<RefCountedPtr<ItoSpecies>>&
  getItoSpecies() const;
 
  const Vector<RefCountedPtr<CdrSpecies>>&
  getCdrSpecies() const;
 
  const Vector<RefCountedPtr<RtSpecies>>&
  getRtSpecies() const;
 
  const Vector<DiffusionFunction>&
  getItoDiffusionFunctions() const noexcept;
 
  virtual Vector<std::string>
  getPlotVariableNames() const noexcept;
 
  virtual Vector<Real>
  getPlotVariables(const RealVect&         a_E,
                   const RealVect&         a_pos,
                   const Vector<Real>&     a_phi,
                   const Vector<RealVect>& a_gradPhi,
                   const Real              a_dx,
                   const Real              a_kappa) const noexcept;
 
  virtual int
  getNumberOfPlotVariables() const noexcept;
 
  inline int
  getNumItoSpecies() const;
 
  inline int
  getNumCdrSpecies() const;
 
  inline int
  getNumPlasmaSpecies() const;
 
  inline int
  getNumPhotonSpecies() const;
 
  virtual bool
  needGradients() const noexcept;
 
  virtual bool
  needGradient(const int a_plasmaIndex) const noexcept;
 
  inline const std::map<int, std::pair<SpeciesType, int>>&
  getSpeciesMap() const noexcept;
 
  inline virtual void
  parseRuntimeOptions() noexcept;
 
  inline virtual Real
  initialSigma(const Real a_time, const RealVect& a_pos) const;
 
  virtual void
  computeMobilities(Vector<Real>&   a_mobilities,
                    const Real      a_time,
                    const RealVect& a_pos,
                    const RealVect& a_E) const noexcept = 0;
 
  virtual void
  computeDiffusionCoefficients(Vector<Real>&   a_diffusionCoefficients,
                               const Real      a_time,
                               const RealVect& a_pos,
                               const RealVect& a_E) const noexcept = 0;
 
  virtual void
  secondaryEmissionEB(Vector<ParticleSoA<ItoParticle>>&       a_secondaryParticles,
                      Vector<Real>&                           a_cdrFluxes,
                      Vector<ParticleSoA<Photon>>&            a_secondaryPhotons,
                      const Vector<ParticleSoA<ItoParticle>>& a_primaryParticles,
                      const Vector<Real>&                     a_cdrFluxesExtrap,
                      const Vector<ParticleSoA<Photon>>&      a_primaryPhotons,
                      const RealVect&                         a_E,
                      const RealVect&                         a_physicalCellCenter,
                      const RealVect&                         a_cellCentroid,
                      const RealVect&                         a_bndryCentroid,
                      const RealVect&                         a_bndryNormal,
                      const Real                              a_bndryArea,
                      const Real                              a_dx,
                      const Real                              a_dt,
                      const bool                              a_isDielectric,
                      const int                               a_matIndex) const noexcept = 0;
 
  inline void
  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;
 
  inline void
  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;
 
  inline void
  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;
 
  inline void
  reconcilePhotoionization(Vector<ParticleSoA<ItoParticle>*>&  a_itoParticles,
                           Vector<ParticleSoA<NoPayload>*>&    a_cdrParticles,
                           const Vector<ParticleSoA<Photon>*>& a_absorbedPhotons) const noexcept;
 
  inline void
  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;
 
protected:
  enum class Algorithm
  {
    SSA,                 
    ExplicitEuler,       
    Midpoint,            
    PRC,                 
    ImplicitEuler,       
    HybridExplicitEuler, 
    HybridMidpoint,      
    HybridPRC,           
    HybridImplicitEuler  
  };
 
  enum class ParticlePlacement
  {
    Random,     
    Centroid,   
    Downstream, 
    Parent      
  };
 
  Algorithm m_algorithm;
 
  ParticlePlacement m_particlePlacement;
 
  std::map<int, std::pair<SpeciesType, int>> m_speciesMap;
 
  std::string m_className;
 
  bool m_debug;
 
  bool m_isDefined;
 
  static thread_local bool m_hasKMCSolver;
 
  static thread_local KMCSolverType m_kmcSolver;
 
  static thread_local KMCState m_kmcState;
 
  static thread_local KMCState m_kmcStateScratch;
 
  static thread_local std::vector<Real> m_kmcPropensityScratch;
 
  static thread_local std::vector<long long> m_weightScratch;
 
  static thread_local std::vector<Real> m_parentWeightScratch;
 
  static thread_local std::vector<std::shared_ptr<const KMCReaction>> m_kmcReactionsThreadLocal;
 
  static thread_local std::vector<std::shared_ptr<const KMCReaction>> m_kmcReactionsDt;
 
  static thread_local std::vector<Real> m_reactiveDtFactorsDt;
 
  static thread_local std::vector<Real> m_kmcPropensityScratchDt;
 
  std::vector<KMCReaction> m_kmcReactions;
 
  std::vector<ItoKMCPhotoReaction> m_photoReactions;
 
  std::vector<Real> m_reactiveDtFactors;
 
  std::map<int, std::pair<std::discrete_distribution<int>, std::map<int, int>>> m_photoPathways;
 
  ItoKMCSurfaceReactionSet m_surfaceReactions;
 
  Vector<DiffusionFunction> m_itoDiffusionFunctions;
 
  Vector<RefCountedPtr<ItoSpecies>> m_itoSpecies;
 
  Vector<RefCountedPtr<CdrSpecies>> m_cdrSpecies;
 
  Vector<RefCountedPtr<RtSpecies>> m_rtSpecies;
 
  int m_downstreamSpecies;
 
  int m_maxNewParticles;
 
  int m_maxNewPhotons;
 
  int m_Ncrit;
 
  int m_NSSA;
 
  int m_maxIter;
 
  Real m_SSAlim;
 
  Real m_eps;
 
  Real m_exitTol;
 
  inline void
  define() noexcept;
 
  inline void
  defineSpeciesMap() noexcept;
 
  inline void
  definePhotoPathways() noexcept;
 
  inline void
  parsePPC() noexcept;
 
  inline void
  parseDebug() noexcept;
 
  inline void
  parseAlgorithm() noexcept;
 
  virtual void
  updateReactionRates(std::vector<std::shared_ptr<const KMCReaction>>& a_kmcReactions,
                      const RealVect                                   a_E,
                      const RealVect                                   a_pos,
                      const Vector<Real>&                              a_phi,
                      const Vector<RealVect>&                          a_gradPhi,
                      const Real                                       a_dt,
                      const Real                                       a_dx,
                      const Real                                       a_kappa) const noexcept = 0;
 
  inline void
  removeParticles(ParticleSoA<ItoParticle>& a_particles, const long long a_numToRemove) const;
 
  inline RealVect
  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;
 
  inline bool
  sampleParentPosition(RealVect&                       a_pos,
                       const ParticleSoA<ItoParticle>& a_particles,
                       const std::vector<Real>&        a_cumulativeWeights) const noexcept;
 
  inline void
  buildParentWeights(std::vector<Real>&              a_cumulativeWeights,
                     const ParticleSoA<ItoParticle>& a_particles,
                     const std::size_t               a_numEligible) const noexcept;
 
  inline bool
  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;
 
  inline RealVect
  noDiffusion(const ItoParticle& a_particle, const Real a_dt) const noexcept;
 
  inline RealVect
  isotropicDiffusion(const ItoParticle& a_particle, const Real a_dt) const noexcept;
 
  inline RealVect
  forwardIsotropicDiffusion(const ItoParticle& a_particle, const Real a_dt) const noexcept;
};
} // namespace ItoKMC
} // namespace Physics
 
#include <CD_NamespaceFooter.H>
 
#include <CD_ItoKMCPhysicsImplem.H>
 
#endif