Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Beta/NumberFractionBeta.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 System of number-fraction beta SDEs
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs) with linear drift and quadratic diagonal
11 : : diffusion, whose invariant is the joint [beta
12 : : distribution](http://en.wikipedia.org/wiki/Beta_distribution). The main
13 : : difference compared to the plain beta SDE (see DiffEq/Beta.h), is that in
14 : : the number-fraction beta SDE the dependent variable, there are two
15 : : additional stochastic variables computed from the beta variables.
16 : :
17 : : In a nutshell, the equation integrated governs a set of scalars,
18 : : \f$0\!\le\!X_\alpha\f$, \f$\alpha\!=\!1,\dots,N\f$, as
19 : : \f[
20 : : \mathrm{d}X_\alpha(t) = \frac{b_\alpha}{2}\left(S_\alpha - X_\alpha\right)
21 : : \mathrm{d}t + \sqrt{\kappa_\alpha X_\alpha(1-X_\alpha)}
22 : : \mathrm{d}W_\alpha(t), \qquad \alpha=1,\dots,N
23 : : \f]
24 : : with parameter vectors \f$b_\alpha > 0\f$, \f$\kappa_\alpha > 0\f$, and \f$0
25 : : < S_\alpha < 1\f$. This is the same as in DiffEq/Beta.h. Here
26 : : \f$\mathrm{d}W_\alpha(t)\f$ is an isotropic vector-valued [Wiener
27 : : process](http://en.wikipedia.org/wiki/Wiener_process) with independent
28 : : increments. The invariant distribution is the joint beta distribution. This
29 : : system of SDEs consists of N independent equations. For
30 : : more on the beta SDE, see https://doi.org/10.1080/14685248.2010.510843.
31 : :
32 : : In addition to integrating the above SDE, there are two additional functions
33 : : of \f$ X_\alpha \f$ are computed as
34 : : \f[ \begin{aligned}
35 : : \rho(X_\alpha) & = \rho_{2\alpha} ( 1 - r'_\alpha X_\alpha ) \\
36 : : V(X_\alpha) & = \frac{1}{ \rho(X\alpha) }
37 : : \end{aligned} \f]
38 : : These equations compute the instantaneous mixture density, \f$ \rho \f$, and
39 : : instantaneous specific volume, \f$ V_\alpha \f$, for equation \f$ \alpha \f$
40 : : in the system. These quantities are used in binary mixing of
41 : : variable-density turbulence between two fluids with constant densities, \f$
42 : : \rho_1, \f$ and \f$ \rho_2 \f$. The additional parameters, \f$ \rho_2 \f$
43 : : and \f$ r' \f$ are user input parameters and kept constant during
44 : : integration. Since we compute the above variables, \f$\rho,\f$ and \f$V\f$,
45 : : and call them mixture density and specific volume, respectively, \f$X\f$,
46 : : governed by the beta SDE is a number (or mole) fraction, hence the name
47 : : number-fraction beta.
48 : :
49 : : _All of this is unpublished, but will be linked in here once published_.
50 : : */
51 : : // *****************************************************************************
52 : : #ifndef NumberFractionBeta_h
53 : : #define NumberFractionBeta_h
54 : :
55 : : #include <vector>
56 : : #include <cmath>
57 : :
58 : : #include "InitPolicy.hpp"
59 : : #include "NumberFractionBetaCoeffPolicy.hpp"
60 : : #include "RNG.hpp"
61 : : #include "Particles.hpp"
62 : :
63 : : namespace walker {
64 : :
65 : : extern ctr::InputDeck g_inputdeck;
66 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
67 : :
68 : : //! \brief NumberFractionBeta SDE used polymorphically with DiffEq
69 : : //! \details The template arguments specify policies and are used to configure
70 : : //! the behavior of the class. The policies are:
71 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
72 : : //! - Coefficients - coefficients policy, see
73 : : //! DiffEq/NumberFractionBetaCoeffPolicy.h
74 : : template< class Init, class Coefficients >
75 : : class NumberFractionBeta {
76 : :
77 : : private:
78 : : using ncomp_t = tk::ctr::ncomp_t;
79 : :
80 : : public:
81 : : //! \brief Constructor
82 : : //! \param[in] c Index specifying which system of number-fraction beta SDEs
83 : : //! to construct. There can be multiple numfracbeta ... end blocks in a
84 : : //! control file. This index specifies which number-fraction beta SDE
85 : : //! system to instantiate. The index corresponds to the order in which the
86 : : //! numfracbeta ... end blocks are given the control file.
87 : 11 : explicit NumberFractionBeta( ncomp_t c ) :
88 : : m_c( c ),
89 : : m_depvar(
90 : 11 : g_inputdeck.get< tag::param, tag::numfracbeta, tag::depvar >().at(c) ),
91 : : m_ncomp(
92 : 11 : g_inputdeck.get< tag::component >().get< tag::numfracbeta >().at(c) / 3 ),
93 : : m_offset(
94 : 11 : g_inputdeck.get< tag::component >().offset< tag::numfracbeta >(c) ),
95 [ + - ]: 11 : m_rng( g_rng.at( tk::ctr::raw(
96 : 11 : g_inputdeck.get< tag::param, tag::numfracbeta, tag::rng >().at(c) ) ) ),
97 : : m_b(),
98 : : m_S(),
99 : : m_k(),
100 : : m_rho2(),
101 : : m_rcomma(),
102 : 11 : coeff( m_ncomp,
103 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::numfracbeta, tag::b >().at(c),
104 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::numfracbeta, tag::S >().at(c),
105 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::numfracbeta, tag::kappa >().at(c),
106 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::numfracbeta, tag::rho2 >().at(c),
107 [ + - ]: 11 : g_inputdeck.get< tag::param, tag::numfracbeta, tag::rcomma >().at(c),
108 [ + - ]: 66 : m_b, m_S, m_k, m_rho2, m_rcomma ) {}
109 : :
110 : : //! Initalize SDE, prepare for time integration
111 : : //! \param[in] stream Thread (or more precisely stream) ID
112 : : //! \param[in,out] particles Array of particle properties
113 : 40 : void initialize( int stream, tk::Particles& particles ) {
114 : : //! Set initial conditions using initialization policy
115 : : Init::template
116 : : init< tag::numfracbeta >
117 : 40 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
118 : 40 : }
119 : :
120 : : //! \brief Advance particles according to the system of number-fraction beta
121 : : //! SDEs
122 : : //! \param[in,out] particles Array of particle properties
123 : : //! \param[in] stream Thread (or more precisely stream) ID
124 : : //! \param[in] dt Time step size
125 : 500000 : void advance( tk::Particles& particles,
126 : : int stream,
127 : : tk::real dt,
128 : : tk::real,
129 : : const std::map< tk::ctr::Product, tk::real >& )
130 : : {
131 : : // Advance particles
132 : 500000 : const auto npar = particles.nunk();
133 [ + + ]: 5487500 : for (auto p=decltype(npar){0}; p<npar; ++p) {
134 : : // Generate Gaussian random numbers with zero mean and unit variance
135 [ + - ]: 9975000 : std::vector< tk::real > dW( m_ncomp );
136 [ + - ]: 4987500 : m_rng.gaussian( stream, m_ncomp, dW.data() );
137 : : // Advance all m_ncomp scalars
138 [ + + ]: 29925000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
139 [ + - ]: 24937500 : tk::real& X = particles( p, i, m_offset );
140 : 24937500 : tk::real d = m_k[i] * X * (1.0 - X) * dt;
141 [ + + ]: 24937500 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
142 : 24937500 : X += 0.5*m_b[i]*(m_S[i] - X)*dt + d*dW[i];
143 : : // Compute instantaneous values derived from updated X
144 [ + - ]: 24937500 : particles( p, m_ncomp+i, m_offset ) = rho( X, i );
145 [ + - ]: 24937500 : particles( p, m_ncomp*2+i, m_offset ) = vol( X, i );
146 : : }
147 : : }
148 : 500000 : }
149 : :
150 : : private:
151 : : const ncomp_t m_c; //!< Equation system index
152 : : const char m_depvar; //!< Dependent variable
153 : : const ncomp_t m_ncomp; //!< Number of components
154 : : const ncomp_t m_offset; //!< Offset SDE operates from
155 : : const tk::RNG& m_rng; //!< Random number generator
156 : :
157 : : //! Coefficients
158 : : std::vector< kw::sde_b::info::expect::type > m_b;
159 : : std::vector< kw::sde_S::info::expect::type > m_S;
160 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
161 : : std::vector< kw::sde_rho2::info::expect::type > m_rho2;
162 : : std::vector< kw::sde_rcomma::info::expect::type > m_rcomma;
163 : :
164 : : //! Coefficients policy
165 : : Coefficients coeff;
166 : :
167 : : //! \brief Return density for mole fraction
168 : : //! \details Functional wrapper around the dependent variable of the beta
169 : : //! SDE. This function returns the instantaneous density, rho,
170 : : //! based on the number fraction, X, and parameters rho2 and r'.
171 : : //! \param[in] X Instantaneous value of the mole fraction, X
172 : : //! \param[in] i Index specifying which (of multiple) parameters to use
173 : : //! \return Instantaneous value of the density, rho
174 : 49875000 : tk::real rho( tk::real X, ncomp_t i ) const {
175 : 49875000 : return m_rho2[i] * ( 1.0 - m_rcomma[i] * X );
176 : : }
177 : :
178 : : //! \brief Return specific volume for mole fraction
179 : : //! \details Functional wrapper around the dependent variable of the beta
180 : : //! SDE. This function returns the instantaneous specific volume, V,
181 : : //! based on the number fraction, X, and parameters rho2 and r'.
182 : : //! \param[in] X Instantaneous value of the mole fraction, X
183 : : //! \param[in] i Index specifying which (of multiple) parameters to use
184 : : //! \return Instantaneous value of the specific volume, V
185 : 24937500 : tk::real vol( tk::real X, ncomp_t i ) const {
186 : 24937500 : return 1.0 / rho( X, i );
187 : : }
188 : : };
189 : :
190 : : } // walker::
191 : :
192 : : #endif // NumberFractionBeta_h
|
