Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/OrnsteinUhlenbeck/DiagOrnsteinUhlenbeck.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 diagonal Ornstein-Uhlenbeck SDEs
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs), with linear drift and constant diagonal
11 : : diffusion, whose invariant is the joint [normal
12 : : distribution](http://en.wikipedia.org/wiki/Normal_distribution).
13 : :
14 : : In a nutshell, the equation integrated governs a set of scalars,
15 : : \f$\color[HTML]{dcdcdc}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) = \theta_\alpha\left(\mu_\alpha - Y_\alpha\right)
22 : : \mathrm{d}t + \sigma_\alpha\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) = \theta_\alpha\left(\mu_\alpha - Y_\alpha\right)
29 : : \mathrm{d}t + \sigma_\alpha\mathrm{d}W_\alpha(t), \\ \alpha=1,\dots,N
30 : : \end{split} \f]
31 : :
32 : : with parameter vectors \f$\color[HTML]{dcdcdc}\theta_\alpha > 0\f$,
33 : : \f$\color[HTML]{dcdcdc}\mu_\alpha\f$, and \f$\color[HTML]{dcdcdc}\sigma_\alpha
34 : : > 0\f$. Here \f$\color[HTML]{dcdcdc}\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 normal
37 : : distribution. This system of SDEs consists of N independent equations, each
38 : : being a single-variate [Ornstein-Uhlenbeck
39 : : process](http://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process).
40 : :
41 : : From the Fokker-Planck equation, equivalent to the SDE above, the equations
42 : : governing the means, \f$\color[HTML]{dcdcdc} \langle Y_\alpha \rangle\f$, are
43 : : \f[
44 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
45 : : \frac{\partial\irmean{Y_\alpha}}{\partial t} =
46 : : \theta_\alpha\left(\mu_\alpha - \irmean{Y_\alpha}\right)
47 : : \f]
48 : : while the equation governing the covariance matrix, \f$\color[HTML]{dcdcdc}
49 : : \langle y_\alpha y_\beta \rangle \equiv \left\langle (Y_\alpha - \langle
50 : : Y_\alpha \rangle) (Y_\beta - \langle Y_\beta\rangle) \right\rangle \f$, is
51 : :
52 : : @m_class{m-show-m}
53 : :
54 : : \f[
55 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
56 : : \newcommand{\irv}[1]{\langle{#1^2}\rangle}
57 : : \frac{\partial\irmean{y_\alpha y_\beta}}{\partial t} =
58 : : -\left(\theta_\alpha+\theta_\beta\right)\irmean{y_\alpha y_\beta} +
59 : : \left\{ \begin{array}{lr}
60 : : \sigma_\alpha^2 & \mathrm{for} \quad \alpha = \beta,\\
61 : : 0 & \mathrm{for} \quad \alpha \ne \beta
62 : : \end{array} \right..
63 : : \f]
64 : :
65 : : @m_class{m-hide-m}
66 : :
67 : : \f[ \begin{split}
68 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
69 : : \newcommand{\irv}[1]{\langle{#1^2}\rangle}
70 : : \frac{\partial\irmean{y_\alpha y_\beta}}{\partial t} =
71 : : -\left(\theta_\alpha+\theta_\beta\right)\irmean{y_\alpha y_\beta} \\ +
72 : : \left\{ \begin{array}{lr}
73 : : \sigma_\alpha^2 & \mathrm{for} \quad \alpha = \beta,\\
74 : : 0 & \mathrm{for} \quad \alpha \ne \beta
75 : : \end{array} \right..
76 : : \end{split} \f]
77 : : */
78 : : // *****************************************************************************
79 : : #ifndef DiagOrnsteinUhlenbeck_h
80 : : #define DiagOrnsteinUhlenbeck_h
81 : :
82 : : #include <vector>
83 : : #include <cmath>
84 : :
85 : : #include "InitPolicy.hpp"
86 : : #include "DiagOrnsteinUhlenbeckCoeffPolicy.hpp"
87 : : #include "RNG.hpp"
88 : : #include "Particles.hpp"
89 : :
90 : : namespace walker {
91 : :
92 : : extern ctr::InputDeck g_inputdeck;
93 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
94 : :
95 : : //! \brief Diagonal Ornstein-Uhlenbeck SDE used polymorphically with DiffEq
96 : : //! \details The template arguments specify policies and are used to configure
97 : : //! the behavior of the class. The policies are:
98 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
99 : : //! - Coefficients - coefficients policy, see
100 : : //! DiffEq/DiagOrnsteinUhlenbeckCoeffPolicy.h
101 : : template< class Init, class Coefficients >
102 : : class DiagOrnsteinUhlenbeck {
103 : :
104 : : private:
105 : : using ncomp_t = tk::ctr::ncomp_t;
106 : :
107 : : public:
108 : : //! \brief Constructor
109 : : //! \param[in] c Index specifying which system of diagonal
110 : : //! Ornstein-Uhlenbeck SDEs to construct. There can be multiple diag_ou
111 : : //! ... end blocks in a control file. This index specifies which diagonal
112 : : //! Ornstein-Uhlenbeck SDE system to instantiate. The index corresponds to
113 : : //! the order in which the diag_ou ... end blocks are given the control
114 : : //! file.
115 : 45 : explicit DiagOrnsteinUhlenbeck( ncomp_t c ) :
116 : : m_c( c ),
117 : 45 : m_depvar( g_inputdeck.get< tag::param, tag::diagou, tag::depvar >().at(c) ),
118 : 45 : m_ncomp( g_inputdeck.get< tag::component >().get< tag::diagou >().at(c) ),
119 : 45 : m_offset( g_inputdeck.get< tag::component >().offset< tag::diagou >(c) ),
120 [ + - ]: 45 : m_rng( g_rng.at( tk::ctr::raw(
121 : 45 : g_inputdeck.get< tag::param, tag::diagou, tag::rng >().at(c) ) ) ),
122 : : m_sigmasq(),
123 : : m_theta(),
124 : : m_mu(),
125 : 45 : coeff( m_ncomp,
126 [ + - ]: 45 : g_inputdeck.get< tag::param, tag::diagou, tag::sigmasq >().at(c),
127 [ + - ]: 45 : g_inputdeck.get< tag::param, tag::diagou, tag::theta >().at(c),
128 [ + - ]: 45 : g_inputdeck.get< tag::param, tag::diagou, tag::mu >().at(c),
129 [ + - ]: 270 : m_sigmasq, m_theta, m_mu ) {}
130 : :
131 : : //! Initalize SDE, prepare for time integration
132 : : //! \param[in] stream Thread (or more precisely stream) ID
133 : : //! \param[in,out] particles Array of particle properties
134 : 153 : void initialize( int stream, tk::Particles& particles ) {
135 : : //! Set initial conditions using initialization policy
136 : : Init::template
137 : : init< tag::diagou >
138 : 153 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
139 : 153 : }
140 : :
141 : : //! \brief Advance particles according to the system of diagonal
142 : : //! Orsntein-Uhlenbeck SDEs
143 : : //! \param[in,out] particles Array of particle properties
144 : : //! \param[in] stream Thread (or more precisely stream) ID
145 : : //! \param[in] dt Time step size
146 : 1530000 : void advance( tk::Particles& particles,
147 : : int stream,
148 : : tk::real dt,
149 : : tk::real,
150 : : const std::map< tk::ctr::Product, tk::real >& )
151 : : {
152 : 1530000 : const auto npar = particles.nunk();
153 [ + + ]: 521530000 : for (auto p=decltype(npar){0}; p<npar; ++p) {
154 : : // Generate Gaussian random numbers with zero mean and unit variance
155 [ + - ]: 1040000000 : std::vector< tk::real > dW( m_ncomp );
156 [ + - ]: 520000000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
157 : :
158 : : // Advance all m_ncomp scalars
159 [ + + ]: 1560000000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
160 [ + - ]: 1040000000 : tk::real& par = particles( p, i, m_offset );
161 : 1040000000 : tk::real d = m_sigmasq[i] * dt;
162 [ + - ]: 1040000000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
163 : 1040000000 : par += m_theta[i]*(m_mu[i] - par)*dt + d*dW[i];
164 : : }
165 : : }
166 : 1530000 : }
167 : :
168 : : private:
169 : : const ncomp_t m_c; //!< Equation system index
170 : : const char m_depvar; //!< Dependent variable
171 : : const ncomp_t m_ncomp; //!< Number of components
172 : : const ncomp_t m_offset; //!< Offset SDE operates from
173 : : const tk::RNG& m_rng; //!< Random number generator
174 : :
175 : : //! Coefficients
176 : : std::vector< kw::sde_sigmasq::info::expect::type > m_sigmasq;
177 : : std::vector< kw::sde_theta::info::expect::type > m_theta;
178 : : std::vector< kw::sde_mu::info::expect::type > m_mu;
179 : :
180 : : //! Coefficients policy
181 : : Coefficients coeff;
182 : : };
183 : :
184 : : } // walker::
185 : :
186 : : #endif // DiagOrnsteinUhlenbeck_h
|
