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