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