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