Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/SkewNormal/SkewNormal.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 skew-normal SDEs
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs), whose invariant is the joint [skew-normal
11 : : distribution](http://en.wikipedia.org/wiki/Skew_normal_distribution).
12 : :
13 : : In a nutshell, the equation integrated governs a set of scalars,
14 : : \f$\color[HTML]{dcdcdc}x_\alpha\f$,
15 : : \f$\color[HTML]{dcdcdc}\alpha\!=\!1,\dots,N\f$, as
16 : :
17 : : @m_class{m-show-m}
18 : :
19 : : \f[
20 : : \mathrm{d}x_\alpha(t) = -\frac{1}{T_\alpha}\left[x_\alpha -
21 : : \lambda_\alpha\sigma^2_\alpha\sqrt{\frac{2}{\pi}} \cdot
22 : : \frac{\exp{\left(-\lambda_\alpha^2x^2_\alpha/2\right)}}{1+\mathrm{erf}
23 : : \left( \lambda_\alpha x_\alpha/\sqrt{2}\right)} \right] \mathrm{d}t +
24 : : \sqrt{\frac{2\sigma^2_\alpha}{T_\alpha}}\mathrm{d}W_\alpha(t).
25 : : \f]
26 : :
27 : : @m_class{m-hide-m}
28 : :
29 : : \f[ \begin{split}
30 : : \mathrm{d}x_\alpha(t) = -\frac{1}{T_\alpha}\left[x_\alpha -
31 : : \lambda_\alpha\sigma^2_\alpha\sqrt{\frac{2}{\pi}} \cdot
32 : : \frac{\exp{\left(-\lambda_\alpha^2x^2_\alpha/2\right)}}{1+\mathrm{erf}
33 : : \left( \lambda_\alpha x_\alpha/\sqrt{2}\right)} \right] \mathrm{d}t \\ +
34 : : \sqrt{\frac{2\sigma^2_\alpha}{T_\alpha}}\mathrm{d}W_\alpha(t).
35 : : \end{split} \f]
36 : :
37 : : The invariant distribution is the joint skew-normal distribution
38 : :
39 : : @m_class{m-show-m}
40 : :
41 : : \f[
42 : : p(x_\alpha) = \frac{1}{\sigma_\alpha\sqrt{2\pi}} \exp\left(
43 : : -\frac{x^2_\alpha}{2\sigma^2_\alpha} \right) \left[1 + \mathrm{erf}\left(
44 : : \frac{\lambda_\alpha x_\alpha}{\sqrt{2}}\right) \right].
45 : : \f]
46 : :
47 : : @m_class{m-hide-m}
48 : :
49 : : \f[ \begin{split}
50 : : p(x_\alpha) = \frac{1}{\sigma_\alpha\sqrt{2\pi}} \exp\left(
51 : : -\frac{x^2_\alpha}{2\sigma^2_\alpha} \right) \\ \times
52 : : \left[1 + \mathrm{erf}\left(
53 : : \frac{\lambda_\alpha x_\alpha}{\sqrt{2}}\right) \right].
54 : : \end{split} \f]
55 : :
56 : : Here \f$\color[HTML]{dcdcdc}\mathrm{erf}(y) = 2/\sqrt{\pi} \int_0^y
57 : : \exp(-u^2) \mathrm{d}u\f$, \f$\color[HTML]{dcdcdc}T_\alpha\f$ are time scales,
58 : : \f$\color[HTML]{dcdcdc}\sigma_\alpha\f$ are variance parameters,
59 : : \f$\color[HTML]{dcdcdc}\mathrm{d}W_\alpha(t)\f$ is an isotropic [Wiener
60 : : process](http://en.wikipedia.org/wiki/Wiener_process) with independent
61 : : increments, and \f$\color[HTML]{dcdcdc}\lambda_\alpha\f$ are the parameters
62 : : that control the asymmetry and skewness for variable
63 : : \f$\color[HTML]{dcdcdc}x_\alpha\f$: \f$\color[HTML]{dcdcdc}\lambda_\alpha<0\f$
64 : : and \f$\color[HTML]{dcdcdc}\lambda_\alpha>0\f$ give negative and positive
65 : : skewness, respectively, while \f$\color[HTML]{dcdcdc}\lambda_\alpha=0\f$
66 : : reduces the system to a set of
67 : : independent Ornstein-Uhlenbeck processes. See also
68 : : DiffEq/DiagOrnsteinUhlenbeck.h. For more details on the skew-normal
69 : : distribution, see http://www.jstor.org/stable/2337278.
70 : : */
71 : : // *****************************************************************************
72 : : #ifndef SkewNormal_h
73 : : #define SkewNormal_h
74 : :
75 : : #include <vector>
76 : : #include <cmath>
77 : : #include <cfenv>
78 : :
79 : : #include "InitPolicy.hpp"
80 : : #include "SkewNormalCoeffPolicy.hpp"
81 : : #include "RNG.hpp"
82 : : #include "Particles.hpp"
83 : :
84 : : namespace walker {
85 : :
86 : : extern ctr::InputDeck g_inputdeck;
87 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
88 : :
89 : : //! \brief Skew-normal SDE used polymorphically with DiffEq
90 : : //! \details The template arguments specify policies and are used to configure
91 : : //! the behavior of the class. The policies are:
92 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
93 : : //! - Coefficients - coefficients policy, see DiffEq/SkewNormalCoeffPolicy.h
94 : : template< class Init, class Coefficients >
95 : : class SkewNormal {
96 : :
97 : : private:
98 : : using ncomp_t = tk::ctr::ncomp_t;
99 : :
100 : : public:
101 : : //! \brief Constructor
102 : : //! \param[in] c Index specifying which system of skew-normal SDEs to
103 : : //! construct. There can be multiple skew-normal ... end blocks in a
104 : : //! control file. This index specifies which skew-normal SDE system to
105 : : //! instantiate. The index corresponds to the order in which the
106 : : //! skew-normal ... end blocks are given the control file.
107 : 45 : explicit SkewNormal( ncomp_t c ) :
108 : : m_c( c ),
109 : : m_depvar(
110 : 45 : g_inputdeck.get< tag::param, tag::skewnormal, tag::depvar >().at(c) ),
111 : : m_ncomp(
112 : 45 : g_inputdeck.get< tag::component >().get< tag::skewnormal >().at(c) ),
113 : : m_offset(
114 : 45 : g_inputdeck.get< tag::component >().offset< tag::skewnormal >(c) ),
115 [ + - ]: 45 : m_rng( g_rng.at( tk::ctr::raw(
116 : 45 : g_inputdeck.get< tag::param, tag::skewnormal, tag::rng >().at(c) ) ) ),
117 : : m_T(),
118 : : m_sigmasq(),
119 : : m_lambda(),
120 : : coeff(
121 : 45 : m_ncomp,
122 [ + - ]: 45 : g_inputdeck.get< tag::param, tag::skewnormal, tag::timescale >().at(c),
123 [ + - ]: 45 : g_inputdeck.get< tag::param, tag::skewnormal, tag::sigmasq >().at(c),
124 [ + - ]: 45 : g_inputdeck.get< tag::param, tag::skewnormal, tag::lambda >().at(c),
125 [ + - ]: 270 : m_T, m_sigmasq, m_lambda ) {}
126 : :
127 : : //! Initalize SDE, prepare for time integration
128 : : //! \param[in] stream Thread (or more precisely stream) ID
129 : : //! \param[in,out] particles Array of particle properties
130 : 153 : void initialize( int stream, tk::Particles& particles ) {
131 : : //! Set initial conditions using initialization policy
132 : : Init::template
133 : : init< tag::skewnormal >
134 : 153 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
135 : 153 : }
136 : :
137 : : //! \brief Advance particles according to the system of skew-normal SDEs
138 : : //! \param[in,out] particles Array of particle properties
139 : : //! \param[in] stream Thread (or more precisely stream) ID
140 : : //! \param[in] dt Time step size
141 : 1530000 : void advance( tk::Particles& particles,
142 : : int stream,
143 : : tk::real dt,
144 : : tk::real,
145 : : const std::map< tk::ctr::Product, tk::real >& )
146 : : {
147 : : fenv_t fe;
148 : 1530000 : feholdexcept( &fe );
149 : :
150 : 1530000 : const auto npar = particles.nunk();
151 [ + + ]: 501530000 : for (auto p=decltype(npar){0}; p<npar; ++p) {
152 : : // Generate Gaussian random numbers with zero mean and unit variance
153 [ + - ]: 1000000000 : std::vector< tk::real > dW( m_ncomp );
154 [ + - ]: 500000000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
155 : :
156 : : // Advance all m_ncomp scalars
157 [ + + ]: 1500000000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
158 [ + - ]: 1000000000 : tk::real& x = particles( p, i, m_offset );
159 : 1000000000 : tk::real d = 2.0 * m_sigmasq[i] / m_T[i] * dt;
160 [ + - ]: 1000000000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
161 : 1000000000 : x += - ( x - m_lambda[i] * m_sigmasq[i]
162 : 1000000000 : * std::sqrt( 2.0 / M_PI )
163 : 1000000000 : * std::exp( - m_lambda[i] * m_lambda[i] * x * x / 2.0 )
164 : 1000000000 : / ( 1.0 + std::erf( m_lambda[i] * x / std::sqrt(2.0) ) )
165 : 1000000000 : ) / m_T[i] * dt
166 : 1000000000 : + d*dW[i];
167 : : }
168 : : }
169 : :
170 : 1530000 : feclearexcept( FE_UNDERFLOW );
171 : 1530000 : feupdateenv( &fe );
172 : 1530000 : }
173 : :
174 : : private:
175 : : const ncomp_t m_c; //!< Equation system index
176 : : const char m_depvar; //!< Dependent variable
177 : : const ncomp_t m_ncomp; //!< Number of components
178 : : const ncomp_t m_offset; //!< Offset SDE operates from
179 : : const tk::RNG& m_rng; //!< Random number generator
180 : :
181 : : //! Coefficients
182 : : std::vector< kw::sde_T::info::expect::type > m_T;
183 : : std::vector< kw::sde_sigmasq::info::expect::type > m_sigmasq;
184 : : std::vector< kw::sde_lambda::info::expect::type > m_lambda;
185 : :
186 : : //! Coefficients policy
187 : : Coefficients coeff;
188 : : };
189 : :
190 : : } // walker::
191 : :
192 : : #endif // SkewNormal_h
|
