Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/OrnsteinUhlenbeck/OrnsteinUhlenbeck.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 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 diffusion,
11 : : whose invariant is the [joint normal
12 : : distribution](http://en.wikipedia.org/wiki/Multivariate_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 + \sum_{\beta=1}^N \sigma_{\alpha\beta}\mathrm{d}W_\beta(t),
23 : : \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) = \theta_\alpha\left(\mu_\alpha - Y_\alpha\right)
30 : : \mathrm{d}t + \sum_{\beta=1}^N \sigma_{\alpha\beta}\mathrm{d}W_\beta(t),
31 : : \\ \alpha=1,\dots,N
32 : : \end{split} \f]
33 : :
34 : : with parameter vectors \f$\color[HTML]{dcdcdc}\theta_\alpha > 0\f$,
35 : : \f$\color[HTML]{dcdcdc}\mu_\alpha\f$, and symmetric positive semi-definite
36 : : diffusion matrix \f$\color[HTML]{dcdcdc}\sigma_{\alpha\beta}\f$. Here
37 : : \f$\color[HTML]{dcdcdc}\mathrm{d}W_\beta(t)\f$ is an isotropic
38 : : vector-valued [Wiener process](http://en.wikipedia.org/wiki/Wiener_process)
39 : : with independent increments. The invariant distribution is the joint normal
40 : : distribution. This system of SDEs consists of N coupled equations, each
41 : : being a single-variate [Ornstein-Uhlenbeck
42 : : process](http://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process).
43 : :
44 : : From the Fokker-Planck equation, equivalent to the SDE above, the equations
45 : : governing the means, \f$\color[HTML]{dcdcdc} \langle Y_\alpha \rangle\f$, are
46 : : \f[
47 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
48 : : \frac{\partial\irmean{Y_\alpha}}{\partial t} =
49 : : \theta_\alpha\left(\mu_\alpha - \irmean{Y_\alpha}\right)
50 : : \f]
51 : : while the equation governing the covariance matrix, \f$\color[HTML]{dcdcdc}
52 : : \langle y_\alpha y_\beta \rangle \equiv \left\langle (Y_\alpha - \langle
53 : : Y_\alpha \rangle) (Y_\beta - \langle Y_\beta\rangle) \right\rangle \f$, is
54 : :
55 : : @m_class{m-show-m}
56 : :
57 : : \f[
58 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
59 : : \newcommand{\irv}[1]{\langle{#1^2}\rangle}
60 : : \frac{\partial\irmean{y_\alpha y_\beta}}{\partial t} =
61 : : -\left(\theta_\alpha+\theta_\beta\right)\irmean{y_\alpha y_\beta}
62 : : +\sum_{\gamma=1}^N \sigma_{\alpha\gamma} \sigma_{\gamma\beta}.
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 : : +\sum_{\gamma=1}^N \sigma_{\alpha\gamma} \sigma_{\gamma\beta}.
73 : : \end{split} \f]
74 : : */
75 : : // *****************************************************************************
76 : : #ifndef OrnsteinUhlenbeck_h
77 : : #define OrnsteinUhlenbeck_h
78 : :
79 : : #include <vector>
80 : : #include <cmath>
81 : :
82 : : #include "WalkerBuildConfig.hpp"
83 : :
84 : : #ifdef HAS_MKL
85 : : #include <mkl_lapacke.h>
86 : : #else
87 : : #include <lapacke.h>
88 : : #endif
89 : :
90 : : #include "InitPolicy.hpp"
91 : : #include "OrnsteinUhlenbeckCoeffPolicy.hpp"
92 : : #include "RNG.hpp"
93 : : #include "Particles.hpp"
94 : :
95 : : namespace walker {
96 : :
97 : : extern ctr::InputDeck g_inputdeck;
98 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
99 : :
100 : : //! \brief Ornstein-Uhlenbeck SDE used polymorphically with DiffEq
101 : : //! \details The template arguments specify policies and are used to configure
102 : : //! the behavior of the class. The policies are:
103 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
104 : : //! - Coefficients - coefficients policy, see
105 : : //! DiffEq/OrnsteinUhlenbeckCoeffPolicy.h
106 : : template< class Init, class Coefficients >
107 : : class OrnsteinUhlenbeck {
108 : :
109 : : private:
110 : : using ncomp_t = tk::ctr::ncomp_t;
111 : :
112 : : public:
113 : : //! \brief Constructor
114 : : //! \param[in] c Index specifying which system of Ornstein-Uhlenbeck SDEs to
115 : : //! construct. There can be multiple ornstein-uhlenbeck ... end blocks in
116 : : //! a control file. This index specifies which Ornstein-Uhlenbeck SDE
117 : : //! system to instantiate. The index corresponds to the order in which the
118 : : //! ornstein-uhlenbeck ... end blocks are given the control file.
119 : 34 : explicit OrnsteinUhlenbeck( ncomp_t c ) :
120 : : m_c( c ),
121 : 34 : m_depvar( g_inputdeck.get< tag::param, tag::ou, tag::depvar >().at(c) ),
122 : 34 : m_ncomp( g_inputdeck.get< tag::component >().get< tag::ou >().at(c) ),
123 : 34 : m_offset( g_inputdeck.get< tag::component >().offset< tag::ou >(c) ),
124 [ + - ]: 34 : m_rng( g_rng.at( tk::ctr::raw(
125 : 34 : g_inputdeck.get< tag::param, tag::ou, tag::rng >().at(c) ) ) ),
126 : : m_sigma(),
127 : : m_theta(),
128 : : m_mu(),
129 : 34 : coeff( m_ncomp,
130 [ + - ]: 34 : g_inputdeck.get< tag::param, tag::ou, tag::sigmasq >().at(c),
131 [ + - ]: 34 : g_inputdeck.get< tag::param, tag::ou, tag::theta >().at(c),
132 [ + - ]: 34 : g_inputdeck.get< tag::param, tag::ou, tag::mu >().at(c),
133 [ + - ]: 204 : m_sigma, m_theta, m_mu )
134 : : {
135 : : // Compute diffusion matrix using Cholesky-decomposition
136 : 34 : lapack_int n = static_cast< lapack_int >( m_ncomp );
137 : : #ifndef NDEBUG
138 : : lapack_int info =
139 : : #endif
140 [ + - ]: 34 : LAPACKE_dpotrf( LAPACK_ROW_MAJOR, 'U', n, m_sigma.data(), n );
141 [ - + ][ - - ]: 34 : Assert( info == 0, "Error in Cholesky-decomposition" );
[ - - ][ - - ]
142 : 34 : }
143 : :
144 : : //! Initalize SDE, prepare for time integration
145 : : //! \param[in] stream Thread (or more precisely stream) ID
146 : : //! \param[in,out] particles Array of particle properties
147 : 106 : void initialize( int stream, tk::Particles& particles ) {
148 : : //! Set initial conditions using initialization policy
149 : : Init::template
150 : : init< tag::ou >
151 : 106 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
152 : 106 : }
153 : :
154 : : //! \brief Advance particles according to the system of Orsntein-Uhlenbeck
155 : : //! SDEs
156 : : //! \param[in,out] particles Array of particle properties
157 : : //! \param[in] stream Thread (or more precisely stream) ID
158 : : //! \param[in] dt Time step size
159 : 53000 : void advance( tk::Particles& particles,
160 : : int stream,
161 : : tk::real dt,
162 : : tk::real,
163 : : const std::map< tk::ctr::Product, tk::real >& )
164 : : {
165 : 53000 : const auto npar = particles.nunk();
166 [ + + ]: 105053000 : for (auto p=decltype(npar){0}; p<npar; ++p) {
167 : : // Generate Gaussian random numbers with zero mean and unit variance
168 [ + - ]: 210000000 : std::vector< tk::real > dW( m_ncomp );
169 [ + - ]: 105000000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
170 : :
171 : : // Advance all m_ncomp scalars
172 [ + + ]: 420000000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
173 [ + - ]: 315000000 : tk::real& par = particles( p, i, m_offset );
174 : 315000000 : par += m_theta[i]*(m_mu[i] - par)*dt;
175 [ + + ]: 1260000000 : for (ncomp_t j=0; j<m_ncomp; ++j) {
176 : 945000000 : tk::real d = m_sigma[ j*m_ncomp+i ] * sqrt(dt); // use transpose
177 : 945000000 : par += d*dW[j];
178 : : }
179 : : }
180 : : }
181 : 53000 : }
182 : :
183 : : private:
184 : : const ncomp_t m_c; //!< Equation system index
185 : : const char m_depvar; //!< Dependent variable
186 : : const ncomp_t m_ncomp; //!< Number of components
187 : : const ncomp_t m_offset; //!< Offset SDE operates from
188 : : const tk::RNG& m_rng; //!< Random number generator
189 : :
190 : : //! Coefficients
191 : : std::vector< kw::sde_sigmasq::info::expect::type > m_sigma;
192 : : std::vector< kw::sde_theta::info::expect::type > m_theta;
193 : : std::vector< kw::sde_mu::info::expect::type > m_mu;
194 : :
195 : : //! Coefficients policy
196 : : Coefficients coeff;
197 : : };
198 : :
199 : : } // walker::
200 : :
201 : : #endif // OrnsteinUhlenbeck_h
|
