Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Dirichlet/GeneralizedDirichlet.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 Lochner's generalized Dirichlet SDE
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs) whose invariant is Lochner's [generalized
11 : : Dirichlet distribution](http://en.wikipedia.org/wiki/Generalized_Dirichlet_distribution).
12 : :
13 : : In a nutshell, the equation integrated governs a set of scalars,
14 : : \f$\color{default}0 \le Y_i\f$, \f$i=1,\dots,K\f$, \f$\sum_{i=1}^KY_i\le1\f$,
15 : : as
16 : :
17 : : @m_class{m-show-m}
18 : :
19 : : \f[ \begin{split}
20 : : \mathrm{d}Y_i(t) = \frac{\mathcal{U}_i}{2}\left\{ b_i\Big[S_i
21 : : \mathcal{Y}_K - (1-S_i)Y_i\Big] + Y_i\mathcal{Y}_K
22 : : \sum_{j=i}^{K-1}\frac{c_{ij}}{\mathcal{Y}_j}\right\}\mathrm{d}t +
23 : : \sqrt{\kappa_i Y_i \mathcal{Y}_K \mathcal{U}_i}\mathrm{d}W_i(t), \\
24 : : \qquad i=1,\dots,K, \end{split}
25 : : \f]
26 : :
27 : : @m_class{m-hide-m}
28 : :
29 : : \f[ \begin{split}
30 : : \mathrm{d}Y_i(t) = \frac{\mathcal{U}_i}{2}\left\{ b_i\Big[S_i
31 : : \mathcal{Y}_K - (1-S_i)Y_i\Big] + Y_i\mathcal{Y}_K
32 : : \sum_{j=i}^{K-1}\frac{c_{ij}}{\mathcal{Y}_j}\right\}\mathrm{d}t \\ +
33 : : \sqrt{\kappa_i Y_i \mathcal{Y}_K \mathcal{U}_i}\mathrm{d}W_i(t),
34 : : \qquad i=1,\dots,K, \end{split}
35 : : \f]
36 : :
37 : : where \f$\mathrm{d}W_i(t)\f$ is an isotropic vector-valued [Wiener
38 : : process](http://en.wikipedia.org/wiki/Wiener_process) with independent
39 : : increments. The statistically stationary solution of the above coupled
40 : : system of nonlinear stochastic differential equations is the generalized
41 : : Dirichlet distribution,
42 : :
43 : : @m_class{m-show-m}
44 : :
45 : : \f[
46 : : \newcommand{\bv}[1]{{\mbox{$\mathbf{#1}$}}}
47 : : G(\bv{Y},\bv{\alpha},\bv{\beta}) =
48 : : \prod_{i=1}^K\frac{\Gamma(\alpha_i+\beta_i)}{\Gamma(\alpha_i)
49 : : \Gamma(\beta_i)}Y_i^{\alpha_i-1} \mathcal{Y}_i^{\gamma_i} \qquad
50 : : \mathrm{with} \qquad \mathcal{Y}_i = 1-\sum_{k=1}^i Y_k,
51 : : \f]
52 : :
53 : : @m_class{m-hide-m}
54 : :
55 : : \f[ \begin{split}
56 : : \newcommand{\bv}[1]{{\mbox{$\mathbf{#1}$}}}
57 : : G(\bv{Y},\bv{\alpha},\bv{\beta}) =
58 : : \prod_{i=1}^K\frac{\Gamma(\alpha_i+\beta_i)}{\Gamma(\alpha_i)
59 : : \Gamma(\beta_i)}Y_i^{\alpha_i-1} \mathcal{Y}_i^{\gamma_i} \\
60 : : \mathrm{with} \qquad \mathcal{Y}_i = 1-\sum_{k=1}^i Y_k,
61 : : \end{split} \f]
62 : :
63 : : provided the coefficients, \f$b_i\!>\!0\f$, \f$\kappa_i\!>\!0\f$,
64 : : \f$0\!<\!S_i\!<\!1\f$, and \f$c_{ij}\f$, with \f$c_{ij}\!=\!0\f$ for
65 : : \f$i\!>\!j\f$, \f$i,j\!=\!1,\dots,K\!-\!1\f$, satisfy
66 : : \f[ \begin{split}
67 : : \alpha_i & = \frac{b_i}{\kappa_i}S_i, \qquad i=1,\dots,K,\\
68 : : 1-\gamma_i & = \frac{c_{1i}}{\kappa_1} = \dots = \frac{c_{ii}}{\kappa_i},
69 : : \qquad i=1,\dots,K-1,\\
70 : : 1+\gamma_K & = \frac{b_1}{\kappa_1}(1-S_1) = \dots =
71 : : \frac{b_K}{\kappa_K}(1-S_K). \end{split}
72 : : \f]
73 : :
74 : : Here \f$\mathcal{U}_i = \prod_{j=1}^{K-i}\mathcal{Y}_{K-j}^{-1}\f$,
75 : : \f$\alpha_i>0\f$, and \f$\beta_i>0\f$ are parameters, while
76 : : \f$\gamma_i=\beta_i-\alpha_{i+1}-\beta_{i+1}\f$ for \f$i=1,\dots,K-1\f$, and
77 : : \f$\gamma_K=\beta_K-1\f$. \f$\Gamma(\cdot)\f$ denotes the [gamma
78 : : function](http://en.wikipedia.org/wiki/Gamma_function). To keep the
79 : : invariant distribution generalized Dirichlet, the above set of constraints
80 : : on the coefficients must be satisfied. For more details on the generalized
81 : : Dirichlet SDE, see https://doi.org/10.1063/1.4822416.
82 : : */
83 : : // *****************************************************************************
84 : : #ifndef GeneralizedDirichlet_h
85 : : #define GeneralizedDirichlet_h
86 : :
87 : : #include <vector>
88 : : #include <cmath>
89 : :
90 : : #include "InitPolicy.hpp"
91 : : #include "GeneralizedDirichletCoeffPolicy.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 Lochner's generalized Dirichlet 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/GeneralizedDirichletCoeffPolicy.h
106 : : template< class Init, class Coefficients >
107 : : class GeneralizedDirichlet {
108 : :
109 : : private:
110 : : using ncomp_t = tk::ctr::ncomp_t;
111 : :
112 : : public:
113 : : //! \brief Constructor
114 : : //! \param[in] c Index specifying which generalized Dirichlet system of SDEs
115 : : //! to construct. There can be multiple gendir ... end blocks in a control
116 : : //! file. This index specifies which generalized Dirichlet SDE system to
117 : : //! instantiate. The index corresponds to the order in which the gendir
118 : 11 : //! ... end blocks are given the control file.
119 : : explicit GeneralizedDirichlet( ncomp_t c ) :
120 : : m_c( c ),
121 : : m_depvar( g_inputdeck.get< tag::param, tag::gendir, tag::depvar >().at(c) ),
122 : 11 : m_ncomp( g_inputdeck.get< tag::component >().get< tag::gendir >().at(c) ),
123 : 11 : m_offset( g_inputdeck.get< tag::component >().offset< tag::gendir >(c) ),
124 : : m_rng( g_rng.at( tk::ctr::raw(
125 : : g_inputdeck.get< tag::param, tag::gendir, tag::rng >().at(c) ) ) ),
126 : : m_b(),
127 : : m_S(),
128 : : m_k(),
129 : 11 : m_cij(),
130 : : coeff( m_ncomp,
131 : : g_inputdeck.get< tag::param, tag::gendir, tag::b >().at(c),
132 : : g_inputdeck.get< tag::param, tag::gendir, tag::S >().at(c),
133 : : g_inputdeck.get< tag::param, tag::gendir, tag::kappa >().at(c),
134 [ - + ][ - + ]: 33 : g_inputdeck.get< tag::param, tag::gendir, tag::c >().at(c),
[ - + ][ - + ]
[ + - ]
135 : : m_b, m_S, m_k, m_cij ) {}
136 : :
137 : : //! Initalize SDE, prepare for time integration
138 : : //! \param[in] stream Thread (or more precisely stream) ID
139 : : //! \param[in,out] particles Array of particle properties
140 : : void initialize( int stream, tk::Particles& particles ) {
141 : : //! Set initial conditions using initialization policy
142 : : Init::template
143 : 0 : init< tag::gendir >
144 : : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
145 : : }
146 : :
147 : : //! \brief Advance particles according to the generalized Dirichlet SDE
148 : : //! \param[in,out] particles Array of particle properties
149 : : //! \param[in] stream Thread (or more precisely stream) ID
150 : 263153 : //! \param[in] dt Time step size
151 : : void advance( tk::Particles& particles,
152 : : int stream,
153 : : tk::real dt,
154 : : tk::real,
155 : : const std::map< tk::ctr::Product, tk::real >& )
156 : : {
157 [ + + ]: 22659153 : const auto npar = particles.nunk();
158 : : for (auto p=decltype(npar){0}; p<npar; ++p) {
159 : 22396000 : // Y_i = 1 - sum_{k=1}^{i} y_k
160 : 22396000 : std::vector< tk::real > Y( m_ncomp );
161 [ + + ]: 44792000 : Y[0] = 1.0 - particles( p, 0, m_offset );
162 : 22396000 : for (ncomp_t i=1; i<m_ncomp; ++i)
163 : : Y[i] = Y[i-1] - particles( p, i, m_offset );
164 : :
165 [ + - ]: 22396000 : // U_i = prod_{j=1}^{K-i} 1/Y_{K-j}
166 : 22396000 : std::vector< tk::real > U( m_ncomp );
167 [ + + ]: 44792000 : U[m_ncomp-1] = 1.0;
168 : 22396000 : for (long i=static_cast<long>(m_ncomp)-2; i>=0; --i) {
169 : 22396000 : auto j = static_cast< std::size_t >( i );
170 : : U[j] = U[j+1]/Y[j];
171 : : }
172 : :
173 [ + - ]: 22396000 : // Generate Gaussian random numbers with zero mean and unit variance
174 [ + - ]: 22396000 : std::vector< tk::real > dW( m_ncomp );
175 : : m_rng.gaussian( stream, m_ncomp, dW.data() );
176 : :
177 : : // Advance first m_ncomp (K=N-1) scalars
178 [ + + ]: 67188000 : ncomp_t k=0;
179 [ + + ]: 44792000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
180 : 44792000 : tk::real& par = particles( p, i, m_offset );
181 [ + + ]: 44792000 : tk::real d = m_k[i] * par * Y[m_ncomp-1] * U[i] * dt;
182 : : d = (d > 0.0 ? std::sqrt(d) : 0.0);
183 [ + + ]: 67188000 : tk::real a=0.0;
184 : 44792000 : for (ncomp_t j=i; j<m_ncomp-1; ++j) a += m_cij[k++]/Y[j];
185 : 44792000 : par += U[i]/2.0*( m_b[i]*( m_S[i]*Y[m_ncomp-1] - (1.0-m_S[i])*par ) +
186 : : par*Y[m_ncomp-1]*a )*dt + d*dW[i];
187 : : }
188 : 263153 : }
189 : : }
190 : :
191 : : private:
192 : : const ncomp_t m_c; //!< Equation system index
193 : : const char m_depvar; //!< Dependent variable
194 : : const ncomp_t m_ncomp; //!< Number of components
195 : : const ncomp_t m_offset; //!< Offset SDE operates from
196 : : const tk::RNG& m_rng; //!< Random number generator
197 : :
198 : : //! Coefficients
199 : : std::vector< kw::sde_b::info::expect::type > m_b;
200 : : std::vector< kw::sde_S::info::expect::type > m_S;
201 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
202 : : std::vector< kw::sde_c::info::expect::type > m_cij;
203 : :
204 : : //! Coefficients policy
205 : : Coefficients coeff;
206 : : };
207 : :
208 : : } // walker::
209 : :
210 : : #endif // GeneralizedDirichlet_h
|
