Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Beta/Beta.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 beta SDEs
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs) with linear drift and quadratic diagonal
11 : : diffusion, whose invariant is the joint [beta
12 : : distribution](http://en.wikipedia.org/wiki/Beta_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}\left(S_\alpha - Y_\alpha\right)
22 : : \mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha(1-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}\left(S_\alpha - Y_\alpha\right)
30 : : \mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha(1-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 vector-valued [Wiener
38 : : process](http://en.wikipedia.org/wiki/Wiener_process) with independent
39 : : increments. The invariant distribution is the joint beta distribution. This
40 : : system of SDEs consists of N independent equations. For more on the beta SDE,
41 : : see https://doi.org/10.1080/14685248.2010.510843.
42 : : */
43 : : // *****************************************************************************
44 : : #ifndef Beta_h
45 : : #define Beta_h
46 : :
47 : : #include <vector>
48 : : #include <cmath>
49 : :
50 : : #include "InitPolicy.hpp"
51 : : #include "BetaCoeffPolicy.hpp"
52 : : #include "RNG.hpp"
53 : : #include "Particles.hpp"
54 : :
55 : : namespace walker {
56 : :
57 : : extern ctr::InputDeck g_inputdeck;
58 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
59 : :
60 : : //! \brief Beta SDE used polymorphically with DiffEq
61 : : //! \details The template arguments specify policies and are used to configure
62 : : //! the behavior of the class. The policies are:
63 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
64 : : //! - Coefficients - coefficients policy, see DiffEq/BetaCoeffPolicy.h
65 : : template< class Init, class Coefficients >
66 : : class Beta {
67 : :
68 : : private:
69 : : using ncomp_t = tk::ctr::ncomp_t;
70 : :
71 : : public:
72 : : //! \brief Constructor
73 : : //! \param[in] c Index specifying which system of beta SDEs to construct.
74 : : //! There can be multiple beta ... end blocks in a control file. This
75 : : //! index specifies which beta SDE system to instantiate. The index
76 : : //! corresponds to the order in which the beta ... end blocks are given
77 : : //! the control file.
78 : 11 : explicit Beta( ncomp_t c ) :
79 : : m_c( c ),
80 : 11 : m_depvar( g_inputdeck.get< tag::param, tag::beta, tag::depvar >().at(c) ),
81 : 11 : m_ncomp( g_inputdeck.get< tag::component >().get< tag::beta >().at(c) ),
82 : 11 : m_offset( g_inputdeck.get< tag::component >().offset< tag::beta >(c) ),
83 [ + - ]: 11 : m_rng( g_rng.at( tk::ctr::raw(
84 : 11 : g_inputdeck.get< tag::param, tag::beta, tag::rng >().at(c) ) ) ),
85 : : m_b(),
86 : : m_S(),
87 : : m_k(),
88 : 11 : coeff( m_ncomp,
89 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::beta, tag::b >().at(c),
90 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::beta, tag::S >().at(c),
91 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::beta, tag::kappa >().at(c),
92 [ + - ]: 66 : m_b, m_S, m_k ) {}
93 : :
94 : : //! Initalize SDE, prepare for time integration
95 : : //! \param[in] stream Thread (or more precisely stream) ID
96 : : //! \param[in,out] particles Array of particle properties
97 : 47 : void initialize( int stream, tk::Particles& particles ) {
98 : : //! Set initial conditions using initialization policy
99 : : Init::template
100 : : init< tag::beta >
101 : 47 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
102 : 47 : }
103 : :
104 : : //! \brief Advance particles according to the system of beta SDEs
105 : : //! \param[in,out] particles Array of particle properties
106 : : //! \param[in] stream Thread (or more precisely stream) ID
107 : : //! \param[in] dt Time step size
108 : 469953 : void advance( tk::Particles& particles,
109 : : int stream,
110 : : tk::real dt,
111 : : tk::real,
112 : : const std::map< tk::ctr::Product, tk::real >& )
113 : : {
114 : 469953 : const auto npar = particles.nunk();
115 [ + + ]: 40465953 : for (auto p=decltype(npar){0}; p<npar; ++p) {
116 : : // Generate Gaussian random numbers with zero mean and unit variance
117 [ + - ]: 79992000 : std::vector< tk::real > dW( m_ncomp );
118 [ + - ]: 39996000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
119 : :
120 : : // Advance all m_ncomp scalars
121 [ + + ]: 239976000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
122 [ + - ]: 199980000 : tk::real& par = particles( p, i, m_offset );
123 : 199980000 : tk::real d = m_k[i] * par * (1.0 - par) * dt;
124 [ + + ]: 199980000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
125 : 199980000 : par += 0.5*m_b[i]*(m_S[i] - par)*dt + d*dW[i];
126 : : }
127 : : }
128 : 469953 : }
129 : :
130 : : private:
131 : : const ncomp_t m_c; //!< Equation system index
132 : : const char m_depvar; //!< Dependent variable
133 : : const ncomp_t m_ncomp; //!< Number of components
134 : : const ncomp_t m_offset; //!< Offset SDE operates from
135 : : const tk::RNG& m_rng; //!< Random number generator
136 : :
137 : : //! Coefficients
138 : : std::vector< kw::sde_b::info::expect::type > m_b;
139 : : std::vector< kw::sde_S::info::expect::type > m_S;
140 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
141 : :
142 : : //! Coefficients policy
143 : : Coefficients coeff;
144 : : };
145 : :
146 : : } // walker::
147 : :
148 : : #endif // Beta_h
|
