CD_KMCDualStateReactionImplem.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_KMCDUALSTATEREACTIONIMPLEM_H
#define CD_KMCDUALSTATEREACTIONIMPLEM_H
 
// Std includes
#include <limits>
 
// Chombo includes
#include <CH_assert.H>
 
// Our includes
#include <CD_KMCDualStateReaction.H>
#include <CD_NamespaceHeader.H>
 
template <typename State, typename T>
inline KMCDualStateReaction<State, T>::KMCDualStateReaction(const std::list<size_t>& a_lhsReactives,
                                                            const std::list<size_t>& a_rhsReactives,
                                                            const std::list<size_t>& a_rhsNonReactives) noexcept
{
  m_lhsReactives    = a_lhsReactives;
  m_rhsReactives    = a_rhsReactives;
  m_rhsNonReactives = a_rhsNonReactives;
 
  CH_assert(a_lhsReactives.size() > 0);
 
  this->computeStateChanges();
}
 
template <typename State, typename T>
inline KMCDualStateReaction<State, T>::~KMCDualStateReaction()
{}
 
template <typename State, typename T>
inline void
KMCDualStateReaction<State, T>::computeStateChanges() noexcept
{
  // Consumed species.
  for (const auto& r : m_lhsReactives) {
    if (m_reactiveStateChange.find(r) == m_reactiveStateChange.end()) {
      m_reactiveStateChange.emplace(r, -1);
    }
    else {
      m_reactiveStateChange[r]--;
    }
  }
 
  // Produced species.
  for (const auto& p : m_rhsReactives) {
    if (m_reactiveStateChange.find(p) == m_reactiveStateChange.end()) {
      m_reactiveStateChange.emplace(p, +1);
    }
    else {
      m_reactiveStateChange[p]++;
    }
  }
 
  // Produced photons.
  for (const auto& p : m_rhsNonReactives) {
    if (m_nonReactiveStateChange.find(p) == m_nonReactiveStateChange.end()) {
      m_nonReactiveStateChange.emplace(p, +1);
    }
    else {
      m_nonReactiveStateChange[p]++;
    }
  }
 
  // Compute the propensity factor that we need when there are bi- or tri-particle reactions involving the same species.
  // We need to do this because for biparticle reactions of N particles of the same type, there are 0.5 * N * (N-1)
  // unique pairs of particles. Or in general when we have a k-order reaction we have 'N choose k' = N!/(k! * (N-k)!)
  // combinations. In the propensity function we compute the propensity as
  //
  //    a = rate * factor * N * (N-1) * (N-2) ...
  //
  // For this to make sense we have
  //
  //    a = rate * N * (N-1) * (N-2) ... (N-k+1) * (N-k)!/(k! * (N-k)!)
  //      = rate * N * (N-1) * (N-2) ... * 1/k!
  //
  // so the factor is just 1/k! (for each species).
 
  m_propensityFactor = 1.0;
 
  std::map<size_t, size_t> reactantNumbers;
  for (const auto& r : m_lhsReactives) {
    if (reactantNumbers.find(r) == reactantNumbers.end()) {
      reactantNumbers.emplace(r, 1);
    }
    else {
      reactantNumbers[r]++;
    }
  }
 
  // Factorial function.
  auto factorial = [](const T& N) -> T {
    T fac = (T)1;
    for (T i = 2; i <= N; i++) {
      fac *= i;
    }
 
    return fac;
  };
 
  for (const auto& rn : reactantNumbers) {
    m_propensityFactor *= 1.0 / factorial(rn.second);
  }
 
  // Compact the reactant tally into the form propensity() consumes, so that the falling factorial
  // can be formed straight off the incoming state rather than off a private copy of it.
  m_lhsOrders.clear();
  m_lhsOrders.reserve(reactantNumbers.size());
 
  for (const auto& rn : reactantNumbers) {
    m_lhsOrders.emplace_back(rn.first, rn.second);
  }
 
  // Dense mirror of m_reactiveStateChange for the per-cell lookups in KMCSolver. The map is keyed by
  // species index, so the highest key fixes the length; species the reaction does not touch keep a
  // zero state change.
  m_reactiveStateChangeDense.clear();
 
  if (!m_reactiveStateChange.empty()) {
    const size_t highestSpecies = m_reactiveStateChange.rbegin()->first;
 
    m_reactiveStateChangeDense.assign(highestSpecies + 1, (T)0);
 
    for (const auto& s : m_reactiveStateChange) {
      m_reactiveStateChangeDense[s.first] = s.second;
    }
  }
}
 
template <typename State, typename T>
inline Real&
KMCDualStateReaction<State, T>::rate() const noexcept
{
  return m_rate;
}
 
template <typename State, typename T>
inline T
KMCDualStateReaction<State, T>::population(const size_t& a_reactant, const State& a_state) noexcept
{
  const auto& reactiveState = a_state.getReactiveState();
 
  CH_assert(reactiveState.size() > a_reactant);
 
  return reactiveState[a_reactant];
}
 
template <typename State, typename T>
inline Real
KMCDualStateReaction<State, T>::propensity(const State& a_state) const noexcept
{
#ifndef NDEBUG
  this->sanityCheck(a_state);
#endif
 
  Real A = m_rate * m_propensityFactor;
 
  const auto& reactiveState = a_state.getReactiveState();
 
  // For a species entering the reaction k times the contribution is the falling factorial
  // X * (X-1) * ... * (X-k+1) -- the same product the previous implementation formed by walking
  // m_lhsReactives and decrementing a private copy of the state after each factor. Taking the order
  // from m_lhsOrders gives the identical product without copying the state.
  for (const auto& lhs : m_lhsOrders) {
    T X = reactiveState[lhs.first];
 
    for (size_t k = 0; k < lhs.second; k++) {
      A *= X;
 
      X -= (T)1;
    }
  }
 
  return A;
}
 
template <typename State, typename T>
inline T
KMCDualStateReaction<State, T>::computeCriticalNumberOfReactions(const State& a_state) const noexcept
{
#ifndef NDEBUG
  this->sanityCheck(a_state);
#endif
 
  T Lj = std::numeric_limits<T>::max();
 
  const auto& reactiveState = a_state.getReactiveState();
 
  for (const auto& s : m_reactiveStateChange) {
 
    const size_t rI   = s.first;
    const T      nuIJ = s.second;
 
    if (nuIJ < 0) {
      Lj = std::min(Lj, reactiveState[rI] / std::abs(nuIJ));
    }
  }
 
  return Lj;
}
 
template <typename State, typename T>
inline const std::list<size_t>&
KMCDualStateReaction<State, T>::getReactants() const noexcept
{
  return m_lhsReactives;
}
 
template <typename State, typename T>
inline std::list<size_t>
KMCDualStateReaction<State, T>::getReactiveProducts() const noexcept
{
  return m_rhsReactives;
}
 
template <typename State, typename T>
inline std::list<size_t>
KMCDualStateReaction<State, T>::getNonReactiveProducts() const noexcept
{
  return m_rhsNonReactives;
}
 
template <typename State, typename T>
inline T
KMCDualStateReaction<State, T>::getStateChange(const size_t a_particleReactant) const noexcept
{
  // Species beyond the dense vector are not part of this reaction, so their state change is zero --
  // the same answer the map lookup gave for an absent key.
  T nuIJ = 0;
 
  if (a_particleReactant < m_reactiveStateChangeDense.size()) {
    nuIJ = m_reactiveStateChangeDense[a_particleReactant];
  }
 
  return nuIJ;
}
 
template <typename State, typename T>
inline void
KMCDualStateReaction<State, T>::advanceState(State& a_state, const T& a_numReactions) const noexcept
{
#ifndef NDEBUG
  this->sanityCheck(a_state);
#endif
  CH_assert(a_state.isValidState());
 
  auto& reactiveState = a_state.getReactiveState();
  auto& photonState   = a_state.getNonReactiveState();
 
  for (const auto& s : m_reactiveStateChange) {
    reactiveState[s.first] += a_numReactions * s.second;
  }
 
  for (const auto& s : m_nonReactiveStateChange) {
    photonState[s.first] += a_numReactions * s.second;
  }
 
  CH_assert(a_state.isValidState());
}
 
template <typename State, typename T>
inline void
KMCDualStateReaction<State, T>::sanityCheck(const State& a_state) const noexcept
{
  const auto& reactiveState = a_state.getReactiveState();
  const auto& photonState   = a_state.getNonReactiveState();
 
  for (const auto& idx : m_lhsReactives) {
    CH_assert(reactiveState.size() > idx);
  }
 
  for (const auto& idx : m_rhsReactives) {
    CH_assert(reactiveState.size() > idx);
  }
 
  for (const auto& idx : m_rhsNonReactives) {
    CH_assert(photonState.size() > idx);
  }
}
 
#include <CD_NamespaceFooter.H>
 
#endif