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