Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Dirichlet/Dirichlet.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 Dirichlet SDE
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs), whose invariant is the [Dirichlet
11 : : distribution](http://en.wikipedia.org/wiki/Dirichlet_distribution).
12 : :
13 : : In a nutshell, the equation integrated governs a set of scalars,
14 : : \f$\color[HTML]{dcdcdc}0\!\le\!Y_\alpha\f$,
15 : : \f$\color[HTML]{dcdcdc}\alpha\!=\!1,\dots,N-1\f$,
16 : : \f$\color[HTML]{dcdcdc}\sum_{\alpha=1}^{N-1}Y_\alpha\!\le\!1\f$, as
17 : :
18 : : @m_class{m-show-m}
19 : :
20 : : \f[
21 : : \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha Y_N -
22 : : (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha
23 : : Y_N}\mathrm{d}W_\alpha(t), \qquad \alpha=1,\dots,N-1
24 : : \f]
25 : :
26 : : @m_class{m-hide-m}
27 : :
28 : : \f[
29 : : \begin{split}
30 : : \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha Y_N -
31 : : (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha
32 : : Y_N}\mathrm{d}W_\alpha(t), \\ \alpha=1,\dots,N-1
33 : : \end{split}
34 : : \f]
35 : :
36 : : with parameter vectors \f$\color[HTML]{dcdcdc}b_\alpha > 0\f$,
37 : : \f$\color[HTML]{dcdcdc}\kappa_\alpha > 0\f$, and \f$\color[HTML]{dcdcdc}0 <
38 : : S_\alpha < 1\f$, and \f$\color[HTML]{dcdcdc}Y_N =
39 : : 1-\sum_{\alpha=1}^{N-1}Y_\alpha\f$. Here
40 : : \f$\color[HTML]{dcdcdc}\mathrm{d}W_\alpha(t)\f$ is an isotropic vector-valued
41 : : [Wiener
42 : : process](http://en.wikipedia.org/wiki/Wiener_process) with independent
43 : : increments. The invariant distribution is the Dirichlet distribution,
44 : : provided the parameters of the drift and diffusion terms satisfy
45 : : \f[
46 : : (1-S_1) b_1 / \kappa_1 = \dots = (1-S_{N-1}) b_{N-1} / \kappa_{N-1}.
47 : : \f]
48 : : To keep the invariant distribution Dirichlet, the above constraint on the
49 : : coefficients must be satisfied. For more details on the Dirichlet SDE, see
50 : : https://doi.org/10.1155/2013/842981.
51 : : */
52 : : // *****************************************************************************
53 : : #ifndef Dirichlet_h
54 : : #define Dirichlet_h
55 : :
56 : : #include <vector>
57 : : #include <cmath>
58 : :
59 : : #include "InitPolicy.hpp"
60 : : #include "DirichletCoeffPolicy.hpp"
61 : : #include "RNG.hpp"
62 : : #include "Particles.hpp"
63 : :
64 : : namespace walker {
65 : :
66 : : extern ctr::InputDeck g_inputdeck;
67 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
68 : :
69 : : //! \brief Dirichlet SDE used polymorphically with DiffEq
70 : : //! \details The template arguments specify policies and are used to configure
71 : : //! the behavior of the class. The policies are:
72 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
73 : : //! - Coefficients - coefficients policy, see DiffEq/DirCoeffPolicy.h
74 : : template< class Init, class Coefficients >
75 : : class Dirichlet {
76 : :
77 : : private:
78 : : using ncomp_t = tk::ctr::ncomp_t;
79 : :
80 : : public:
81 : : //! \brief Constructor
82 : : //! \param[in] c Index specifying which Dirichlet SDE to construct. There
83 : : //! can be multiple dirichlet ... end blocks in a control file. This index
84 : : //! specifies which Dirichlet SDE to instantiate. The index corresponds to
85 : : //! the order in which the dirichlet ... end blocks are given the control
86 : : //! file.
87 : 33 : explicit Dirichlet( ncomp_t c ) :
88 : : m_c( c ),
89 : : m_depvar(
90 : 33 : g_inputdeck.get< tag::param, tag::dirichlet, tag::depvar >().at(c) ),
91 : : m_ncomp(
92 : 33 : g_inputdeck.get< tag::component >().get< tag::dirichlet >().at(c) ),
93 : : m_offset(
94 : 33 : g_inputdeck.get< tag::component >().offset< tag::dirichlet >(c) ),
95 [ + - ]: 33 : m_rng( g_rng.at( tk::ctr::raw(
96 : 33 : g_inputdeck.get< tag::param, tag::dirichlet, tag::rng >().at(c) ) ) ),
97 : : m_b(),
98 : : m_S(),
99 : : m_k(),
100 : 33 : coeff( m_ncomp,
101 [ + - ]: 33 : g_inputdeck.get< tag::param, tag::dirichlet, tag::b >().at(c),
102 [ + - ]: 33 : g_inputdeck.get< tag::param, tag::dirichlet, tag::S >().at(c),
103 [ + - ]: 33 : g_inputdeck.get< tag::param, tag::dirichlet, tag::kappa >().at(c),
104 [ + - ]: 198 : m_b, m_S, m_k ) {}
105 : :
106 : : //! Initalize SDE, prepare for time integration
107 : : //! \param[in] stream Thread (or more precisely stream) ID
108 : : //! \param[in,out] particles Array of particle properties
109 : 141 : void initialize( int stream, tk::Particles& particles ) {
110 : : //! Set initial conditions using initialization policy
111 : : Init::template
112 : : init< tag::dirichlet >
113 : 141 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
114 : 141 : }
115 : :
116 : : //! \brief Advance particles according to the Dirichlet SDE
117 : : //! \param[in,out] particles Array of particle properties
118 : : //! \param[in] stream Thread (or more precisely stream) ID
119 : : //! \param[in] dt Time step size
120 : 394800 : void advance( tk::Particles& particles,
121 : : int stream,
122 : : tk::real dt,
123 : : tk::real,
124 : : const std::map< tk::ctr::Product, tk::real >& )
125 : : {
126 : 394800 : const auto npar = particles.nunk();
127 [ + + ]: 336394800 : for (auto p=decltype(npar){0}; p<npar; ++p) {
128 : : // Compute Nth scalar
129 [ + - ]: 336000000 : tk::real yn = 1.0 - particles(p, 0, m_offset);
130 [ + + ]: 672000000 : for (ncomp_t i=1; i<m_ncomp; ++i)
131 [ + - ]: 336000000 : yn -= particles( p, i, m_offset );
132 : :
133 : : // Generate Gaussian random numbers with zero mean and unit variance
134 [ + - ]: 672000000 : std::vector< tk::real > dW( m_ncomp );
135 [ + - ]: 336000000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
136 : :
137 : : // Advance first m_ncomp (K=N-1) scalars
138 [ + + ]: 1008000000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
139 [ + - ]: 672000000 : tk::real& par = particles( p, i, m_offset );
140 : 672000000 : tk::real d = m_k[i] * par * yn * dt;
141 [ + + ]: 672000000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
142 : 672000000 : par += 0.5*m_b[i]*( m_S[i]*yn - (1.0-m_S[i]) * par )*dt + d*dW[i];
143 : : }
144 : : }
145 : 394800 : }
146 : :
147 : : private:
148 : : const ncomp_t m_c; //!< Equation system index
149 : : const char m_depvar; //!< Dependent variable
150 : : const ncomp_t m_ncomp; //!< Number of components
151 : : const ncomp_t m_offset; //!< Offset SDE operates from
152 : : const tk::RNG& m_rng; //!< Random number generator
153 : :
154 : : //! Coefficients
155 : : std::vector< kw::sde_b::info::expect::type > m_b;
156 : : std::vector< kw::sde_S::info::expect::type > m_S;
157 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
158 : :
159 : : //! Coefficients policy
160 : : Coefficients coeff;
161 : : };
162 : :
163 : : } // walker::
164 : :
165 : : #endif // Dirichlet_h
|
