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