chombo-discharge
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Protected Attributes | List of all members
KMCSolver< R, State, T > Class Template Reference

Class for running Kinetic Monte-Carlo simulations. More...

#include <CD_KMCSolver.H>

Public Types

using ReactionList = std::vector< std::shared_ptr< const R > >
 

Public Member Functions

 KMCSolver () noexcept
 Default constructor – must subsequently define the object.
 
 KMCSolver (const KMCSolver &)=default
 Disallowed copy constructor.
 
 KMCSolver (const KMCSolver &&)=delete
 Disallowed move constructor.
 
 KMCSolver (const ReactionList &a_reactions) noexcept
 Full constructor.
 
virtual ~KMCSolver () noexcept
 Destructor.
 
KMCSolver & operator= (const KMCSolver &)=default
 Copy assignment operator.
 
KMCSolver & operator= (const KMCSolver &&)=delete
 Disallowed move assignment operator.
 
void define (const ReactionList &a_reactions) noexcept
 Define function. Sets the reactions.
 
void setSolverParameters (const T a_numCrit, const T a_numSSA, const T a_maxIter, const Real a_eps, const Real a_SSAlim, const Real a_exitTol) noexcept
 Set solver parameters.
 
std::vector< std::vector< T > > getNu (const State &a_state, const ReactionList &a_reactions) const noexcept
 Compute the state vector changes for all reactions.
 
std::vector< Real > propensities (const State &a_state) const noexcept
 Compute propensities for ALL reactions.
 
std::vector< Real > propensities (const State &a_state, const ReactionList &a_reactions) const noexcept
 Compute propensities for a subset of reactions.
 
Real totalPropensity (const State &a_state) const noexcept
 Compute the total propensity for ALL reactions.
 
Real totalPropensity (const State &a_state, const ReactionList &a_reactions) const noexcept
 Compute the total propensity for a subset of reactions.
 
std::pair< ReactionList, ReactionList > partitionReactions (const State &a_state) const noexcept
 Partition reactions into critical and non-critical reactions.
 
std::pair< ReactionList, ReactionList > partitionReactions (const State &a_state, const ReactionList &a_reactions) const noexcept
 Partition reactions into critical and non-critical reactions.
 
Real getCriticalTimeStep (const State &a_state) const noexcept
 Get the time to the next critical reaction.
 
Real getCriticalTimeStep (const State &a_state, const ReactionList &a_criticalReactions) const noexcept
 Get the time to the next critical reaction.
 
Real getCriticalTimeStep (const std::vector< Real > &a_propensities) const noexcept
 Get the time to the next critical reaction.
 
Real getCriticalTimeStep (const Real &a_totalPropensity) const noexcept
 Get the time to the next critical reaction.
 
Real getNonCriticalTimeStep (const State &a_state) const noexcept
 Get the non-critical time step.
 
Real getNonCriticalTimeStep (const State &a_state, const ReactionList &a_reactions) const noexcept
 Get the non-critical time step.
 
Real getNonCriticalTimeStep (const State &a_state, const ReactionList &a_nonCriticalReactions, const std::vector< Real > &a_nonCriticalPropensities) const noexcept
 Get the non-critical time step.
 
Real computeDt (const State &a_state, const ReactionList &a_reactions, const std::vector< Real > &a_propensities, const Real a_epsilon) const noexcept
 Compute a time step, using the leap condition on the mean value with eps = 1.
 
void stepSSA (State &a_state) const noexcept
 Perform a single SSA step.
 
void stepSSA (State &a_state, const ReactionList &a_reactions) const noexcept
 Perform a single SSA step.
 
void stepSSA (State &a_state, const ReactionList &a_reactions, const std::vector< Real > &a_propensities) const noexcept
 Perform a single SSA step. This version has pre-computed propensities (for optimization reasons)
 
void advanceSSA (State &a_state, const Real a_dt) const noexcept
 Advance with the SSA over the input time. This can end up using substepping.
 
void advanceSSA (State &a_state, const ReactionList &a_reactions, const Real a_dt) const noexcept
 Advance with the SSA over the input time. This can end up using substepping.
 
void stepExplicitEuler (State &a_state, const Real a_dt) const noexcept
 Perform one plain tau-leaping step using ALL reactions.
 
void stepExplicitEuler (State &a_state, const ReactionList &a_reactions, const Real a_dt) const noexcept
 Perform one plain tau-leaping step over the input reactions using a time step a_dt.
 
void stepMidpoint (State &a_state, const Real a_dt) const noexcept
 Perform one leaping step using the midpoint method for ALL reactions over a time step a_dt.
 
void stepMidpoint (State &a_state, const ReactionList &a_reactions, const Real a_dt) const noexcept
 Perform one leaping step using the midpoint method for the input reactions over a time step a_dt.
 
void stepPRC (State &a_state, const Real a_dt) const noexcept
 Perform one leaping step using the PRC method for ALL reactions over a time step a_dt.
 
void stepPRC (State &a_state, const ReactionList &a_reactions, const Real a_dt) const noexcept
 Perform one leaping step using the PRC method for the input reactions over a time step a_dt.
 
void stepImplicitEuler (State &a_state, const Real a_dt) const noexcept
 Perform one implicit Euler tau-leaping step using ALL reactions.
 
void stepImplicitEuler (State &a_state, const ReactionList &a_reactions, const Real a_dt) const noexcept
 Perform one implicit Euler tau-leaping step over the input reactions using a time step a_dt.
 
void advanceTau (State &a_state, const Real &a_dt, const KMCLeapPropagator &a_leapPropagator=KMCLeapPropagator::ExplicitEuler) const noexcept
 Advance using a specified tau-leaping algorithm.
 
void advanceTau (State &a_state, const ReactionList &a_reactions, const Real &a_dt, const KMCLeapPropagator &a_leapPropagator=KMCLeapPropagator::ExplicitEuler) const noexcept
 Advance using a specified tau-leaping algorithm.
 
void advanceHybrid (State &a_state, const Real a_dt, const KMCLeapPropagator &a_leapPropagator=KMCLeapPropagator::ExplicitEuler) const noexcept
 Advance using Cao et. al. hybrid algorithm over the input time. This can end up using substepping.
 
void advanceHybrid (State &a_state, const ReactionList &a_reactions, const Real a_dt, const KMCLeapPropagator &a_leapPropagator=KMCLeapPropagator::ExplicitEuler) const noexcept
 Advance using Cao et. al. hybrid algorithm over the input time. This can end up using substepping.
 
void advanceHybrid (State &a_state, const ReactionList &a_reactions, const Real a_dt, const std::function< void(State &, const ReactionList &a_reactions, const Real a_dt)> &a_propagator) const noexcept
 Advance using Cao et. al. hybrid algorithm over the input time. This can end up using substepping.
 

Protected Attributes

ReactionList m_reactions
 List of reactions used when advancing states.
 
T m_Ncrit
 Definition of critical reactions.
 
T m_numSSA
 Maximum number of SSA steps to run when switching into SSA-based advancement for non-critical reactions.
 
T m_maxIter
 Maximum number of iterations for implicit solvers.
 
Real m_eps
 Maximum permitted change in propensities for non-critical reactions.
 
Real m_SSAlim
 Threshold for switching to SSA-based algorithm within the Cao algorithm.
 
Real m_exitTol
 Exit tolerance for implicit solvers.
 
std::vector< std::vector< int > > m_nu
 List of state changes for each reaction.
 

Detailed Description

template<typename R, typename State, typename T = long long>
class KMCSolver< R, State, T >

Class for running Kinetic Monte-Carlo simulations.

The template parameter State is the underlying state type that KMC operators on. There are rather simple required member functions on the State parameter, but the reaction type (template parameter R) MUST be able to operate on the state through the following functions:

  1. Real R::propensity(State) const -> Computes the reaction propensity for the input state.
  2. T R::computeCriticalNumberOfReactions(State) const -> Computes the minimum number of reactions that exhausts one of the reactants.
  3. void R::advanceState(State&, const T numReactions) const -> Advance state by numReactions
  4. std::<some_container> getReactants() const -> Get reactants involved in the reactions.
  5. static T R::population(const <some_type> reactant, const State& a_state) -> Get the population of the input reactant in the input state.

The template parameter T should agree across both both R, State, and KMCSolver. Note that this must be a signed integer type.

Constructor & Destructor Documentation

◆ KMCSolver()

template<typename R , typename State , typename T >
KMCSolver< R, State, T >::KMCSolver ( const ReactionList &  a_reactions)
inlinenoexcept

Full constructor.

Parameters
[in]a_reactionsList of reactions.

Member Function Documentation

◆ advanceHybrid() [1/3]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::advanceHybrid ( State &  a_state,
const ReactionList &  a_reactions,
const Real  a_dt,
const KMCLeapPropagator &  a_leapPropagator = KMCLeapPropagator::ExplicitEuler 
) const
inlinenoexcept

Advance using Cao et. al. hybrid algorithm over the input time. This can end up using substepping.

Parameters
[in,out]a_stateState vector to advance.
[in]a_reactionsReactions to advance.
[in]a_dtTime increment.
[in]a_leapPropagatorWhich leap propagator to use.
Note
Calls the other version with stepTau as the leap propagator.

◆ advanceHybrid() [2/3]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::advanceHybrid ( State &  a_state,
const ReactionList &  a_reactions,
const Real  a_dt,
const std::function< void(State &, const ReactionList &a_reactions, const Real a_dt)> &  a_propagator 
) const
inlinenoexcept

Advance using Cao et. al. hybrid algorithm over the input time. This can end up using substepping.

Parameters
[in,out]a_stateState vector to advance
[in]a_reactionsReactions to advance with
[in]a_dtTime increment
[in]a_leapPropagatorLeaping propagator

◆ advanceHybrid() [3/3]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::advanceHybrid ( State &  a_state,
const Real  a_dt,
const KMCLeapPropagator &  a_leapPropagator = KMCLeapPropagator::ExplicitEuler 
) const
inlinenoexcept

Advance using Cao et. al. hybrid algorithm over the input time. This can end up using substepping.

Parameters
[in,out]a_stateState vector to advance
[in]a_dtTime increment
[in]a_leapPropagatorWhich leap propagator to use.
Note
Calls the other version with m_reactions

◆ advanceSSA() [1/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::advanceSSA ( State &  a_state,
const ReactionList &  a_reactions,
const Real  a_dt 
) const
inlinenoexcept

Advance with the SSA over the input time. This can end up using substepping.

Parameters
[in,out]a_stateState vector to advance
[in]a_reactionsReactions to advance with
[in]a_dtTime increment

◆ advanceSSA() [2/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::advanceSSA ( State &  a_state,
const Real  a_dt 
) const
inlinenoexcept

Advance with the SSA over the input time. This can end up using substepping.

Parameters
[in,out]a_stateState vector to advance
[in]a_dtTime increment
Note
Calls the other version with m_reactions

◆ advanceTau() [1/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::advanceTau ( State &  a_state,
const ReactionList &  a_reactions,
const Real &  a_dt,
const KMCLeapPropagator &  a_leapPropagator = KMCLeapPropagator::ExplicitEuler 
) const
inlinenoexcept

Advance using a specified tau-leaping algorithm.

Parameters
[in,out]a_stateState vector to advance.
[in]a_reactionsReactions to advance.
[in]a_dtTime increment.
[in]a_leapPropagatorWhich leap propagator to use.
Note
Calls the other version with stepTau as the leap propagator.

◆ advanceTau() [2/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::advanceTau ( State &  a_state,
const Real &  a_dt,
const KMCLeapPropagator &  a_leapPropagator = KMCLeapPropagator::ExplicitEuler 
) const
inlinenoexcept

Advance using a specified tau-leaping algorithm.

Parameters
[in,out]a_stateState vector to advance
[in]a_dtTime increment
[in]a_leapPropagatorWhich leap propagator to use.
Note
Calls the other version with m_reactions

◆ computeDt()

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::computeDt ( const State &  a_state,
const ReactionList &  a_reactions,
const std::vector< Real > &  a_propensities,
const Real  a_epsilon 
) const
inlinenoexcept

Compute a time step, using the leap condition on the mean value with eps = 1.

This is just like getNonCriticalTimeStep EXCEPT that the leap condition is not applied for the deviation

Parameters
[in]a_stateState vector
[in]a_reactionsReactions
[in]a_propensitiesPropensities
[in]a_epsilonEpsilon value for maximum change in propensities

◆ define()

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::define ( const ReactionList &  a_reactions)
inlinenoexcept

Define function. Sets the reactions.

Parameters
[in]a_reactionsList of reactions.

◆ getCriticalTimeStep() [1/4]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::getCriticalTimeStep ( const Real &  a_totalPropensity) const
inlinenoexcept

Get the time to the next critical reaction.

Parameters
[in]a_totalPropensityTotal propensity

◆ getCriticalTimeStep() [2/4]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::getCriticalTimeStep ( const State &  a_state) const
inlinenoexcept

Get the time to the next critical reaction.

Parameters
[in]a_stateState vector
Note
Calls the other version with m_reactions

◆ getCriticalTimeStep() [3/4]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::getCriticalTimeStep ( const State &  a_state,
const ReactionList &  a_criticalReactions 
) const
inlinenoexcept

Get the time to the next critical reaction.

Parameters
[in]a_stateState vector
[in]a_criticalReactionsReaction list.
Note
Computes the total propensity and calls the other version.

◆ getCriticalTimeStep() [4/4]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::getCriticalTimeStep ( const std::vector< Real > &  a_propensities) const
inlinenoexcept

Get the time to the next critical reaction.

Parameters
[in]a_propensitiesReaction propensities

◆ getNonCriticalTimeStep() [1/3]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::getNonCriticalTimeStep ( const State &  a_state) const
inlinenoexcept

Get the non-critical time step.

Parameters
[in]a_stateState vector
Note
Calls the other version with m_reactions

◆ getNonCriticalTimeStep() [2/3]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::getNonCriticalTimeStep ( const State &  a_state,
const ReactionList &  a_nonCriticalReactions,
const std::vector< Real > &  a_nonCriticalPropensities 
) const
inlinenoexcept

Get the non-critical time step.

Parameters
[in]a_stateState vector
[in]a_nonCriticalReactionsNon-critical reactions
[in]a_nonCriticalPropensitiesNon-critical propensities

◆ getNonCriticalTimeStep() [3/3]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::getNonCriticalTimeStep ( const State &  a_state,
const ReactionList &  a_reactions 
) const
inlinenoexcept

Get the non-critical time step.

Parameters
[in]a_stateState vector
[in]a_reactionsReaction list
Note
Computes propensities and calls the other version.

◆ getNu()

template<typename R , typename State , typename T >
std::vector< std::vector< T > > KMCSolver< R, State, T >::getNu ( const State &  a_state,
const ReactionList &  a_reactions 
) const
inlinenoexcept

Compute the state vector changes for all reactions.

Parameters
[in]a_stateInput state.
[in]a_reactionsReaction list.
Note
The input parameter a_state can be a dummy state because the state change vector does not depend on it.

◆ partitionReactions() [1/2]

template<typename R , typename State , typename T >
std::pair< typename KMCSolver< R, State, T >::ReactionList, typename KMCSolver< R, State, T >::ReactionList > KMCSolver< R, State, T >::partitionReactions ( const State &  a_state) const
inlinenoexcept

Partition reactions into critical and non-critical reactions.

First member in the pair is the critical reactions, then the non-critical reactions.

Parameters
[in]a_stateState vector
Note
Calls the other version with m_reactions

◆ partitionReactions() [2/2]

template<typename R , typename State , typename T >
std::pair< typename KMCSolver< R, State, T >::ReactionList, typename KMCSolver< R, State, T >::ReactionList > KMCSolver< R, State, T >::partitionReactions ( const State &  a_state,
const ReactionList &  a_reactions 
) const
inlinenoexcept

Partition reactions into critical and non-critical reactions.

First member in the pair is the critical reactions, then the non-critical reactions

Parameters
[in]a_stateState vector
[in]a_reactionsReaction list to be partitioned

◆ propensities() [1/2]

template<typename R , typename State , typename T >
std::vector< Real > KMCSolver< R, State, T >::propensities ( const State &  a_state) const
inlinenoexcept

Compute propensities for ALL reactions.

Parameters
[in]a_stateState vector
Note
Calls the other version with m_reactions.

◆ propensities() [2/2]

template<typename R , typename State , typename T >
std::vector< Real > KMCSolver< R, State, T >::propensities ( const State &  a_state,
const ReactionList &  a_reactions 
) const
inlinenoexcept

Compute propensities for a subset of reactions.

Parameters
[in]a_stateState vector
[in]a_reactionsReaction list

◆ setSolverParameters()

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::setSolverParameters ( const T  a_numCrit,
const T  a_numSSA,
const T  a_maxIter,
const Real  a_eps,
const Real  a_SSAlim,
const Real  a_exitTol 
)
inlinenoexcept

Set solver parameters.

Parameters
[in]a_numCritDetermines critical reactions. This is the number of reactions that need to fire before depleting a reactant.
[in]a_numSSAMaximum number of SSA steps to run when switching from tau-leaping to SSA (hybrid algorithm only).
[in]a_maxIterMaximum permitted number of iterations for the implicit solver
[in]a_epsMaximum permitted change in propensities when performing tau-leaping for non-critical reactions.
[in]a_SSAlimThreshold for switching from tau-leaping of non-critical reactions to SSA for all reactions (hybrid algorithm only)
[in]a_exitTolExit tolerance for implicit solvres.

◆ stepExplicitEuler() [1/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepExplicitEuler ( State &  a_state,
const ReactionList &  a_reactions,
const Real  a_dt 
) const
inlinenoexcept

Perform one plain tau-leaping step over the input reactions using a time step a_dt.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_reactionsList of reactions to advance with
[in]a_dtTime increment

◆ stepExplicitEuler() [2/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepExplicitEuler ( State &  a_state,
const Real  a_dt 
) const
inlinenoexcept

Perform one plain tau-leaping step using ALL reactions.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_dtTime increment
Note
Calls the other version with m_reactions

◆ stepImplicitEuler() [1/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepImplicitEuler ( State &  a_state,
const ReactionList &  a_reactions,
const Real  a_dt 
) const
inlinenoexcept

Perform one implicit Euler tau-leaping step over the input reactions using a time step a_dt.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_reactionsList of reactions to advance with
[in]a_dtTime increment

◆ stepImplicitEuler() [2/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepImplicitEuler ( State &  a_state,
const Real  a_dt 
) const
inlinenoexcept

Perform one implicit Euler tau-leaping step using ALL reactions.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_dtTime increment
Note
Calls the other version with m_reactions

◆ stepMidpoint() [1/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepMidpoint ( State &  a_state,
const ReactionList &  a_reactions,
const Real  a_dt 
) const
inlinenoexcept

Perform one leaping step using the midpoint method for the input reactions over a time step a_dt.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_reactionsList of reactions to advance with
[in]a_dtTime increment

◆ stepMidpoint() [2/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepMidpoint ( State &  a_state,
const Real  a_dt 
) const
inlinenoexcept

Perform one leaping step using the midpoint method for ALL reactions over a time step a_dt.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_dtTime increment
Note
Calls the other version with m_reactions

◆ stepPRC() [1/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepPRC ( State &  a_state,
const ReactionList &  a_reactions,
const Real  a_dt 
) const
inlinenoexcept

Perform one leaping step using the PRC method for the input reactions over a time step a_dt.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_reactionsList of reactions to advance with
[in]a_dtTime increment

◆ stepPRC() [2/2]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepPRC ( State &  a_state,
const Real  a_dt 
) const
inlinenoexcept

Perform one leaping step using the PRC method for ALL reactions over a time step a_dt.

Parameters
[in,out]a_stateState vector to be advanced
[in]a_dtTime increment
Note
Calls the other version with m_reactions

◆ stepSSA() [1/3]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepSSA ( State &  a_state) const
inlinenoexcept

Perform a single SSA step.

Parameters
[in,out]a_stateState vector to advance
Note
Calls the other version with m_reactions

◆ stepSSA() [2/3]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepSSA ( State &  a_state,
const ReactionList &  a_reactions 
) const
inlinenoexcept

Perform a single SSA step.

Parameters
[in,out]a_stateState vector to advance
[in]a_reactionsReactions to advance with
Note
Computes propensities and calls the other version

◆ stepSSA() [3/3]

template<typename R , typename State , typename T >
void KMCSolver< R, State, T >::stepSSA ( State &  a_state,
const ReactionList &  a_reactions,
const std::vector< Real > &  a_propensities 
) const
inlinenoexcept

Perform a single SSA step. This version has pre-computed propensities (for optimization reasons)

Parameters
[in,out]a_stateState vector to advance
[in]a_reactionsReactions to advance with
[in]a_propensitiesPropensities for the reactions

◆ totalPropensity() [1/2]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::totalPropensity ( const State &  a_state) const
inlinenoexcept

Compute the total propensity for ALL reactions.

Parameters
[in]a_stateState vector
Note
Calls the other version with m_reactions.

◆ totalPropensity() [2/2]

template<typename R , typename State , typename T >
Real KMCSolver< R, State, T >::totalPropensity ( const State &  a_state,
const ReactionList &  a_reactions 
) const
inlinenoexcept

Compute the total propensity for a subset of reactions.

Parameters
[in]a_stateState vector
[in]a_reactionsReaction list

Member Data Documentation

◆ m_Ncrit

template<typename R , typename State , typename T = long long>
T KMCSolver< R, State, T >::m_Ncrit
protected

Definition of critical reactions.

A reaction is critical if it is m_Ncrit firings away from depleting a reactant.

◆ m_nu

template<typename R , typename State , typename T = long long>
std::vector<std::vector<int> > KMCSolver< R, State, T >::m_nu
protected

List of state changes for each reaction.

Note
The outer vector corresponds to the reactive index. The inner vector describes the change in the population of the state vector, with an indexing given by the State's linearIn/linearOut functions.
The documentation for this class was generated from the following files:
Generated by doxygen 1.9.8