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