Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Gamma/Gamma.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief System of gamma SDEs
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs), with linear drift and linear diagonal
11 : : diffusion, whose invariant is the joint [gamma
12 : : distribution](http://en.wikipedia.org/wiki/Gamma_distribution).
13 : :
14 : : In a nutshell, the equation integrated governs a set of scalars,
15 : : \f$\color[HTML]{dcdcdc}0\!\le\!Y_\alpha\f$,
16 : : \f$\color[HTML]{dcdcdc}\alpha\!=\!1,\dots,N\f$, as
17 : :
18 : : @m_class{m-show-m}
19 : :
20 : : \f[
21 : : \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha -
22 : : (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha}
23 : : \mathrm{d}W_\alpha(t), \qquad \alpha=1,\dots,N
24 : : \f]
25 : :
26 : : @m_class{m-hide-m}
27 : :
28 : : \f[ \begin{split}
29 : : \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha -
30 : : (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha}
31 : : \mathrm{d}W_\alpha(t), \\ \alpha=1,\dots,N
32 : : \end{split} \f]
33 : :
34 : : with parameter vectors \f$\color[HTML]{dcdcdc}b_\alpha > 0\f$,
35 : : \f$\color[HTML]{dcdcdc}\kappa_\alpha > 0\f$, and \f$\color[HTML]{dcdcdc}0 <
36 : : S_\alpha < 1\f$. Here \f$\color[HTML]{dcdcdc}\mathrm{d}W_\alpha(t)\f$ is an
37 : : isotropic
38 : : vector-valued [Wiener process](http://en.wikipedia.org/wiki/Wiener_process)
39 : : with independent increments. The invariant distribution is the joint gamma
40 : : distribution. This system of SDEs consists of N independent equations.
41 : : */
42 : : // *****************************************************************************
43 : : #ifndef Gamma_h
44 : : #define Gamma_h
45 : :
46 : : #include <vector>
47 : : #include <cmath>
48 : :
49 : : #include "InitPolicy.hpp"
50 : : #include "GammaCoeffPolicy.hpp"
51 : : #include "RNG.hpp"
52 : : #include "Particles.hpp"
53 : :
54 : : namespace walker {
55 : :
56 : : extern ctr::InputDeck g_inputdeck;
57 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
58 : :
59 : : //! \brief Gamma SDE used polymorphically with DiffEq
60 : : //! \details The template arguments specify policies and are used to configure
61 : : //! the behavior of the class. The policies are:
62 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
63 : : //! - Coefficients - coefficients policy, see DiffEq/GammaCoeffPolicy.h
64 : : template< class Init, class Coefficients >
65 : : class Gamma {
66 : :
67 : : private:
68 : : using ncomp_t = tk::ctr::ncomp_t;
69 : :
70 : : public:
71 : : //! \brief Constructor
72 : : //! \param[in] c Index specifying which system of gamma SDEs to construct.
73 : : //! There can be multiple gamma ... end blocks in a control file. This
74 : : //! index specifies which gamma SDE system to instantiate. The index
75 : : //! corresponds to the order in which the gamma ... end blocks are given
76 : : //! the control file.
77 : 11 : explicit Gamma( ncomp_t c ) :
78 : : m_c( c ),
79 : 11 : m_depvar( g_inputdeck.get< tag::param, tag::gamma, tag::depvar >().at(c) ),
80 : 11 : m_ncomp( g_inputdeck.get< tag::component >().get< tag::gamma >().at(c) ),
81 : 11 : m_offset( g_inputdeck.get< tag::component >().offset< tag::gamma >(c) ),
82 [ + - ]: 11 : m_rng( g_rng.at( tk::ctr::raw(
83 : 11 : g_inputdeck.get< tag::param, tag::gamma, tag::rng >().at(c) ) ) ),
84 : : m_b(),
85 : : m_S(),
86 : : m_k(),
87 : 11 : coeff( m_ncomp,
88 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::gamma, tag::b >().at(c),
89 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::gamma, tag::S >().at(c),
90 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::gamma, tag::kappa >().at(c),
91 [ + - ]: 66 : m_b, m_S, m_k ) {}
92 : :
93 : : //! Initalize SDE, prepare for time integration
94 : : //! \param[in] stream Thread (or more precisely stream) ID
95 : : //! \param[in,out] particles Array of particle properties
96 : 47 : void initialize( int stream, tk::Particles& particles ) {
97 : : //! Set initial conditions using initialization policy
98 : : Init::template
99 : : init< tag::gamma >
100 : 47 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
101 : 47 : }
102 : :
103 : : //! \brief Advance particles according to the system of gamma SDEs
104 : : //! \param[in,out] particles Array of particle properties
105 : : //! \param[in] stream Thread (or more precisely stream) ID
106 : : //! \param[in] dt Time step size
107 : 164453 : void advance( tk::Particles& particles,
108 : : int stream,
109 : : tk::real dt,
110 : : tk::real,
111 : : const std::map< tk::ctr::Product, tk::real >& )
112 : : {
113 : 164453 : const auto npar = particles.nunk();
114 [ + + ]: 70144453 : for (auto p=decltype(npar){0}; p<npar; ++p) {
115 : : // Generate Gaussian random numbers with zero mean and unit variance
116 [ + - ]: 139960000 : std::vector< tk::real > dW( m_ncomp );
117 [ + - ]: 69980000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
118 : :
119 : : // Advance all m_ncomp scalars
120 [ + + ]: 209940000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
121 [ + - ]: 139960000 : tk::real& par = particles( p, i, m_offset );
122 : 139960000 : tk::real d = m_k[i] * par * dt;
123 [ + + ]: 139960000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
124 : 139960000 : par += 0.5*m_b[i]*(m_S[i] - (1.0 - m_S[i])*par)*dt + d*dW[i];
125 : : }
126 : : }
127 : 164453 : }
128 : :
129 : : private:
130 : : const ncomp_t m_c; //!< Equation system index
131 : : const char m_depvar; //!< Dependent variable
132 : : const ncomp_t m_ncomp; //!< Number of components
133 : : const ncomp_t m_offset; //!< Offset SDE operates from
134 : : const tk::RNG& m_rng; //!< Random number generator
135 : :
136 : : //! Coefficients
137 : : std::vector< kw::sde_b::info::expect::type > m_b;
138 : : std::vector< kw::sde_S::info::expect::type > m_S;
139 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
140 : :
141 : : //! Coefficients policy
142 : : Coefficients coeff;
143 : : };
144 : :
145 : : } // walker::
146 : :
147 : : #endif // Gamma_h
|
