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