chombo-discharge
Loading...
Searching...
No Matches
CD_KMCSolver.H
Go to the documentation of this file.
1/*
2 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
3 *
4 * SPDX-License-Identifier: GPL-3.0-or-later
5 */
6
13#ifndef CD_KMCSOLVER_H
14#define CD_KMCSOLVER_H
15
16// Std includes
17#include <functional>
18#include <memory>
19#include <vector>
20
21// Chombo includes
22#include <REAL.H>
23
24// Our includes
25#include <CD_NamespaceHeader.H>
26
35{
37 Midpoint,
38 PRC,
40};
41
59template <typename R, typename State, typename T = long long>
61{
62public:
66 using ReactionList = std::vector<std::shared_ptr<const R>>;
67
71 KMCSolver() noexcept;
72
76 KMCSolver(const KMCSolver&) = default;
77
81 KMCSolver(const KMCSolver&&) = delete;
82
87 inline KMCSolver(const ReactionList& a_reactions) noexcept;
88
92 virtual ~KMCSolver() noexcept;
93
100 operator=(const KMCSolver& a_other) = default;
101
105 KMCSolver&
106 operator=(const KMCSolver&&) = delete;
107
112 inline void
113 define(const ReactionList& a_reactions) noexcept;
114
128 inline void
129 setSolverParameters(T a_numCrit, T a_numSSA, T a_maxIter, Real a_eps, Real a_SSAlim, Real a_exitTol) noexcept;
130
138 inline std::vector<std::vector<T>>
139 getNu(const State& a_state, const ReactionList& a_reactions) const noexcept;
140
147 inline std::vector<Real>
148 propensities(const State& a_state) const noexcept;
149
156 inline std::vector<Real>
157 propensities(const State& a_state, const ReactionList& a_reactions) const noexcept;
158
165 inline Real
166 totalPropensity(const State& a_state) const noexcept;
167
174 inline Real
175 totalPropensity(const State& a_state, const ReactionList& a_reactions) const noexcept;
176
184 inline std::pair<ReactionList, ReactionList>
185 partitionReactions(const State& a_state) const noexcept;
186
194 inline std::pair<ReactionList, ReactionList>
195 partitionReactions(const State& a_state, const ReactionList& a_reactions) const noexcept;
196
203 inline Real
204 getCriticalTimeStep(const State& a_state) const noexcept;
205
213 inline Real
214 getCriticalTimeStep(const State& a_state, const ReactionList& a_criticalReactions) const noexcept;
215
221 inline Real
222 getCriticalTimeStep(const std::vector<Real>& a_propensities) const noexcept;
223
229 inline Real
230 getCriticalTimeStep(const Real& a_totalPropensity) const noexcept;
231
238 inline Real
239 getNonCriticalTimeStep(const State& a_state) const noexcept;
240
248 inline Real
249 getNonCriticalTimeStep(const State& a_state, const ReactionList& a_reactions) const noexcept;
250
273 inline Real
274 getNonCriticalTimeStep(const State& a_state,
275 const ReactionList& a_reactions,
276 const ReactionList& a_nonCriticalReactions,
277 const std::vector<Real>& a_nonCriticalPropensities) const noexcept;
278
288 inline Real
289 computeDt(const State& a_state,
290 const ReactionList& a_reactions,
291 const std::vector<Real>& a_propensities,
292 Real a_epsilon) const noexcept;
293
299 inline void
300 stepSSA(State& a_state) const noexcept;
301
308 inline void
309 stepSSA(State& a_state, const ReactionList& a_reactions) const noexcept;
310
317 inline void
318 stepSSA(State& a_state, const ReactionList& a_reactions, const std::vector<Real>& a_propensities) const noexcept;
319
326 inline void
327 advanceSSA(State& a_state, Real a_dt) const noexcept;
328
335 inline void
336 advanceSSA(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
337
344 inline void
345 stepExplicitEuler(State& a_state, Real a_dt) const noexcept;
346
353 inline void
354 stepExplicitEuler(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
355
362 inline void
363 stepMidpoint(State& a_state, Real a_dt) const noexcept;
364
371 inline void
372 stepMidpoint(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
373
380 inline void
381 stepPRC(State& a_state, Real a_dt) const noexcept;
382
389 inline void
390 stepPRC(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
391
398 inline void
399 stepImplicitEuler(State& a_state, Real a_dt) const noexcept;
400
407 inline void
408 stepImplicitEuler(State& a_state, const ReactionList& a_reactions, Real a_dt) const noexcept;
409
417 inline void
418 advanceTau(State& a_state,
419 const Real& a_dt,
420 const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
421
433 inline void
434 advanceTau(State& a_state,
435 const ReactionList& a_reactions,
436 const Real& a_dt,
437 const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
438
446 inline void
447 advanceHybrid(State& a_state,
448 Real a_dt,
449 const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
450
459 inline void
460 advanceHybrid(State& a_state,
461 const ReactionList& a_reactions,
462 Real a_dt,
463 const KMCLeapPropagator& a_leapPropagator = KMCLeapPropagator::ExplicitEuler) const noexcept;
464
480 inline void
482 State& a_state,
483 const ReactionList& a_reactions,
484 Real a_dt,
485 const std::function<void(State&, const ReactionList& a_reactions, const Real a_dt)>& a_propagator) const noexcept;
486
487protected:
492
498
503
508
512 Real m_eps;
513
518
523
529 std::vector<std::vector<int>> m_nu;
530
541 mutable std::vector<Real> m_muScratch;
542
548 mutable std::vector<Real> m_sigmaScratch;
549
556 mutable std::vector<Real> m_lossScratch;
557
563 mutable std::vector<size_t> m_reactantScratch;
564
570 mutable std::vector<char> m_seenScratch;
571
580 inline void
581 gatherDistinctReactants(const ReactionList& a_reactions) const noexcept;
582};
583
584#include <CD_NamespaceFooter.H>
585
586#include <CD_KMCSolverImplem.H>
587
588#endif
Implementation of CD_KMCSolver.H.
KMCLeapPropagator
Supported propagators for hybrid tau leaping.
Definition CD_KMCSolver.H:35
@ ImplicitEuler
Implicit Euler tau leaping.
@ Midpoint
Gillespie's midpoint method.
@ ExplicitEuler
Regular tau leaping.
@ PRC
Hu and Li's Poisson random correction method.
Class for running Kinetic Monte-Carlo simulations.
Definition CD_KMCSolver.H:61
std::vector< std::shared_ptr< const R > > ReactionList
Alias for the list of reactions.
Definition CD_KMCSolver.H:66
std::vector< size_t > m_reactantScratch
Scratch holding the distinct reactants of the reaction list currently being processed.
Definition CD_KMCSolver.H:563
Real m_eps
Maximum permitted change in propensities for non-critical reactions.
Definition CD_KMCSolver.H:512
void setSolverParameters(T a_numCrit, T a_numSSA, T a_maxIter, Real a_eps, Real a_SSAlim, Real a_exitTol) noexcept
Set solver parameters.
Definition CD_KMCSolverImplem.H:59
std::vector< Real > m_sigmaScratch
Scratch accumulator for the per-reactant propensity variance, indexed by species.
Definition CD_KMCSolver.H:548
Real computeDt(const State &a_state, const ReactionList &a_reactions, const std::vector< Real > &a_propensities, Real a_epsilon) const noexcept
Compute a time step using the leap condition on the mean value.
Definition CD_KMCSolverImplem.H:399
std::vector< std::vector< int > > m_nu
List of state changes for each reaction.
Definition CD_KMCSolver.H:529
void define(const ReactionList &a_reactions) noexcept
Define function. Sets the reactions.
Definition CD_KMCSolverImplem.H:49
KMCSolver() noexcept
Default constructor. Must subsequently call define.
Definition CD_KMCSolverImplem.H:32
std::vector< std::vector< T > > getNu(const State &a_state, const ReactionList &a_reactions) const noexcept
Compute the state vector changes for all reactions.
Definition CD_KMCSolverImplem.H:76
ReactionList m_reactions
List of reactions used when advancing states.
Definition CD_KMCSolver.H:491
Real m_exitTol
Exit tolerance for implicit solvers.
Definition CD_KMCSolver.H:522
std::vector< Real > m_muScratch
Scratch accumulator for the per-reactant expected state change, indexed by species.
Definition CD_KMCSolver.H:541
Real getCriticalTimeStep(const State &a_state) const noexcept
Get the time to the next critical reaction.
Definition CD_KMCSolverImplem.H:190
T m_maxIter
Maximum number of iterations for implicit solvers.
Definition CD_KMCSolver.H:507
void advanceSSA(State &a_state, Real a_dt) const noexcept
Advance with the SSA over the input time. This can end up using substepping.
Definition CD_KMCSolverImplem.H:526
void gatherDistinctReactants(const ReactionList &a_reactions) const noexcept
Fill m_reactantScratch with the distinct reactants of the input reactions.
Definition CD_KMCSolverImplem.H:260
T m_numSSA
Maximum number of SSA steps to run when switching into SSA-based advancement for non-critical reactio...
Definition CD_KMCSolver.H:502
T m_Ncrit
Definition of critical reactions.
Definition CD_KMCSolver.H:497
void stepImplicitEuler(State &a_state, Real a_dt) const noexcept
Perform one implicit Euler tau-leaping step using ALL reactions.
Definition CD_KMCSolverImplem.H:664
void advanceHybrid(State &a_state, 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.
Definition CD_KMCSolverImplem.H:943
void stepPRC(State &a_state, Real a_dt) const noexcept
Perform one leaping step using the PRC method for ALL reactions.
Definition CD_KMCSolverImplem.H:625
void stepMidpoint(State &a_state, Real a_dt) const noexcept
Perform one leaping step using the midpoint method for ALL reactions.
Definition CD_KMCSolverImplem.H:590
std::vector< char > m_seenScratch
Membership stamps used to build m_reactantScratch without sorting, indexed by species.
Definition CD_KMCSolver.H:570
std::pair< ReactionList, ReactionList > partitionReactions(const State &a_state) const noexcept
Partition reactions into critical and non-critical reactions.
Definition CD_KMCSolverImplem.H:153
void stepSSA(State &a_state) const noexcept
Perform a single SSA step.
Definition CD_KMCSolverImplem.H:465
void advanceTau(State &a_state, const Real &a_dt, const KMCLeapPropagator &a_leapPropagator=KMCLeapPropagator::ExplicitEuler) const noexcept
Advance using a specified tau-leaping algorithm.
Definition CD_KMCSolverImplem.H:863
Real getNonCriticalTimeStep(const State &a_state) const noexcept
Get the non-critical time step.
Definition CD_KMCSolverImplem.H:240
std::vector< Real > m_lossScratch
Scratch accumulator for the per-reactant expected consumption, indexed by species.
Definition CD_KMCSolver.H:556
Real totalPropensity(const State &a_state) const noexcept
Compute the total propensity for ALL reactions.
Definition CD_KMCSolverImplem.H:131
Real m_SSAlim
Threshold for switching to SSA-based algorithm within the Cao algorithm.
Definition CD_KMCSolver.H:517
void stepExplicitEuler(State &a_state, Real a_dt) const noexcept
Perform one plain tau-leaping step using ALL reactions.
Definition CD_KMCSolverImplem.H:561
std::vector< Real > propensities(const State &a_state) const noexcept
Compute propensities for ALL reactions.
Definition CD_KMCSolverImplem.H:109