CD_KMCSolver.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_KMCSOLVER_H
#define CD_KMCSOLVER_H
 
// Std includes
#include <functional>
#include <memory>
#include <vector>
 
// Chombo includes
#include <REAL.H>
 
// Our includes
#include <CD_NamespaceHeader.H>
 
enum class KMCLeapPropagator
{
  ExplicitEuler, 
  Midpoint,      
  PRC,           
  ImplicitEuler  
};
 
template <typename R, typename State, typename T = long long>
class KMCSolver
{
public:
  using ReactionList = std::vector<std::shared_ptr<const R>>;
 
  KMCSolver() noexcept;
 
  KMCSolver(const KMCSolver&) = default;
 
  KMCSolver(const KMCSolver&&) = delete;
 
  inline KMCSolver(const ReactionList& a_reactions) noexcept;
 
  virtual ~KMCSolver() noexcept;
 
  KMCSolver&
  operator=(const KMCSolver& a_other) = default;
 
  KMCSolver&
  operator=(const KMCSolver&&) = delete;
 
  inline void
  define(const ReactionList& a_reactions) noexcept;
 
  inline void
  setSolverParameters(T a_numCrit, T a_numSSA, T a_maxIter, Real a_eps, Real a_SSAlim, Real a_exitTol) noexcept;
 
  inline std::vector<std::vector<T>>
  getNu(const State& a_state, const ReactionList& a_reactions) const noexcept;
 
  inline std::vector<Real>
  propensities(const State& a_state) const noexcept;
 
  inline std::vector<Real>
  propensities(const State& a_state, const ReactionList& a_reactions) const noexcept;
 
  inline Real
  totalPropensity(const State& a_state) const noexcept;
 
  inline Real
  totalPropensity(const State& a_state, const ReactionList& a_reactions) const noexcept;
 
  inline std::pair<ReactionList, ReactionList>
  partitionReactions(const State& a_state) const noexcept;
 
  inline std::pair<ReactionList, ReactionList>
  partitionReactions(const State& a_state, const ReactionList& a_reactions) const noexcept;
 
  inline Real
  getCriticalTimeStep(const State& a_state) const noexcept;
 
  inline Real
  getCriticalTimeStep(const State& a_state, const ReactionList& a_criticalReactions) const noexcept;
 
  inline Real
  getCriticalTimeStep(const std::vector<Real>& a_propensities) const noexcept;
 
  inline Real
  getCriticalTimeStep(const Real& a_totalPropensity) const noexcept;
 
  inline Real
  getNonCriticalTimeStep(const State& a_state) const noexcept;
 
  inline Real
  getNonCriticalTimeStep(const State& a_state, const ReactionList& a_reactions) const noexcept;
 
  inline Real
  getNonCriticalTimeStep(const State&             a_state,
                         const ReactionList&      a_nonCriticalReactions,
                         const std::vector<Real>& a_nonCriticalPropensities) const noexcept;
 
  inline Real
  computeDt(const State&             a_state,
            const ReactionList&      a_reactions,
            const std::vector<Real>& a_propensities,
            Real                     a_epsilon) const noexcept;
 
  inline void
  stepSSA(State& a_state) const noexcept;
 
  inline void
  stepSSA(State& a_state, const ReactionList& a_reactions) const noexcept;
 
  inline void
  stepSSA(State& a_state, const ReactionList& a_reactions, const std::vector<Real>& a_propensities) const noexcept;
 
  inline void
  advanceSSA(State& a_state, Real a_dt) const noexcept;
 
  inline void
  advanceSSA(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
 
  inline void
  stepExplicitEuler(State& a_state, Real a_dt) const noexcept;
 
  inline void
  stepExplicitEuler(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
 
  inline void
  stepMidpoint(State& a_state, Real a_dt) const noexcept;
 
  inline void
  stepMidpoint(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
 
  inline void
  stepPRC(State& a_state, Real a_dt) const noexcept;
 
  inline void
  stepPRC(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
 
  inline void
  stepImplicitEuler(State& a_state, Real a_dt) const noexcept;
 
  inline void
  stepImplicitEuler(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
 
  inline void
  advanceTau(State&                   a_state,
             const Real&              a_dt,
             const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
 
  inline void
  advanceTau(State&                   a_state,
             const ReactionList&      a_reactions,
             const Real&              a_dt,
             const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
 
  inline void
  advanceHybrid(State&                   a_state,
                Real                     a_dt,
                const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
 
  inline void
  advanceHybrid(State&                   a_state,
                const ReactionList&      a_reactions,
                Real                     a_dt,
                const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
 
  inline void
  advanceHybrid(
    State&                                                                               a_state,
    const ReactionList&                                                                  a_reactions,
    Real                                                                                 a_dt,
    const std::function<void(State&, const ReactionList& a_reactions, const Real a_dt)>& a_propagator) const noexcept;
 
protected:
  ReactionList m_reactions;
 
  T m_Ncrit;
 
  T m_numSSA;
 
  T m_maxIter;
 
  Real m_eps;
 
  Real m_SSAlim;
 
  Real m_exitTol;
 
  std::vector<std::vector<int>> m_nu;
 
  mutable std::vector<Real> m_muScratch;
 
  mutable std::vector<Real> m_sigmaScratch;
 
  mutable std::vector<size_t> m_reactantScratch;
 
  mutable std::vector<char> m_seenScratch;
 
  inline void
  gatherDistinctReactants(const ReactionList& a_reactions) const noexcept;
};
 
#include <CD_NamespaceFooter.H>
 
#include <CD_KMCSolverImplem.H>
 
#endif