Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Beta/MixNumberFractionBeta.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 mix 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). There are two
13 : : differences compared to the plain beta SDE (see DiffEq/Beta.h):
14 : :
15 : : - First, the parameters, b, and kappa are specified via functions that
16 : : constrain the beta SDE to be consistent with the turbulent mixing process.
17 : : In particular, the SDE is made consistent with the no-mix and fully mixed
18 : : limits. See, e.g., MixNumberFractionBetaCoeffConst::update().
19 : :
20 : : - Second, there two additional random variables computed, the same as also
21 : : computed by the number-fraction beta equation, see also
22 : : DiffEq/NumberFractionBeta.h.
23 : :
24 : : In a nutshell, the equation integrated governs a set of scalars,
25 : : \f$\color[HTML]{dcdcdc}0\!\le\!X_\alpha\f$,
26 : : \f$\color[HTML]{dcdcdc}\alpha\!=\!1,\dots,N\f$, as
27 : :
28 : : @m_class{m-show-m}
29 : :
30 : : \f[
31 : : \mathrm{d}X_\alpha(t) = \frac{b_\alpha}{2}\left(S_\alpha - X_\alpha\right)
32 : : \mathrm{d}t + \sqrt{\kappa_\alpha X_\alpha(1-X_\alpha)}
33 : : \mathrm{d}W_\alpha(t), \qquad \alpha=1,\dots,N
34 : : \f]
35 : :
36 : : @m_class{m-hide-m}
37 : :
38 : : \f[ \begin{split}
39 : : \mathrm{d}X_\alpha(t) = \frac{b_\alpha}{2}\left(S_\alpha - X_\alpha\right)
40 : : \mathrm{d}t + \sqrt{\kappa_\alpha X_\alpha(1-X_\alpha)}
41 : : \mathrm{d}W_\alpha(t), \\ \alpha=1,\dots,N
42 : : \end{split} \f]
43 : :
44 : : with parameter vectors \f$\color[HTML]{dcdcdc}b_\alpha = \Theta b'_\alpha >
45 : : 0\f$, \f$\color[HTML]{dcdcdc} \newcommand{\irv}[1]{\langle{#1^2}\rangle}
46 : : \kappa_\alpha = \kappa' \irv{x} > 0\f$, and \f$\color[HTML]{dcdcdc}0 < S_\alpha
47 : : < 1\f$. This is similar to DiffEq/Beta.h, but the parameters,
48 : : \f$\color[HTML]{dcdcdc}b\f$ and \f$\color[HTML]{dcdcdc}\kappa\f$ constrained.
49 : : Here \f$\color[HTML]{dcdcdc} \newcommand{\irv}[1]{\langle{#1^2}\rangle}
50 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}} \Theta = 1 - \irv{x} / [
51 : : \irmean{X} (1-\irmean{X}) ]\f$. The fluctuation about the mean,
52 : : \f$\color[HTML]{dcdcdc} \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
53 : : \irmean{X} \f$, is defined as usual: \f$\color[HTML]{dcdcdc}
54 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}} x = X - \irmean{X} \f$, and
55 : : \f$\color[HTML]{dcdcdc}b'\f$ and \f$\color[HTML]{dcdcdc} \kappa'\f$ are
56 : : user-specified constants. Also, \f$\color[HTML]{dcdcdc}\mathrm{d}W_\alpha(t)\f$
57 : : is an isotropic vector-valued
58 : : [Wiener process](http://en.wikipedia.org/wiki/Wiener_process) with
59 : : independent increments. The invariant distribution is the joint beta
60 : : distribution. This system of SDEs consists of N independent equations. For
61 : : more on the beta SDE, see https://doi.org/10.1080/14685248.2010.510843.
62 : :
63 : : Similar to the number-fraction beta SDE (DiffEq/NumberFractionBeta.h), in
64 : : addition to integrating the above SDE, there are two additional functions
65 : : of \f$\color[HTML]{dcdcdc} X_\alpha \f$ are computed as
66 : : \f[ \begin{aligned}
67 : : \rho(X_\alpha) & = \rho_{2\alpha} ( 1 - r'_\alpha X_\alpha ) \\
68 : : V(X_\alpha) & = \frac{1}{ \rho_{2\alpha} ( 1 - r'_\alpha X_\alpha ) }
69 : : \end{aligned} \f]
70 : : These equations compute the instantaneous mixture density,
71 : : \f$\color[HTML]{dcdcdc} \rho \f$, and instantaneous specific volume,
72 : : \f$\color[HTML]{dcdcdc} V_\alpha \f$, for equation \f$\color[HTML]{dcdcdc}
73 : : \alpha \f$ in the system. These quantities are used in binary mixing of
74 : : variable-density turbulence between two fluids with constant densities,
75 : : \f$\color[HTML]{dcdcdc} \rho_1, \f$ and \f$\color[HTML]{dcdcdc} \rho_2 \f$. The
76 : : additional parameters, \f$\color[HTML]{dcdcdc} \rho_2 \f$ and
77 : : \f$\color[HTML]{dcdcdc} r' \f$ are user input parameters and kept constant
78 : : during integration. Since we compute the above variables,
79 : : \f$\color[HTML]{dcdcdc}\rho,\f$ and \f$\color[HTML]{dcdcdc}V\f$, and call them
80 : : mixture density and specific volume, respectively, \f$\color[HTML]{dcdcdc}X\f$,
81 : : governed by the beta SDE is a number (or mole) fraction.
82 : :
83 : : _All of this is unpublished, but will be linked in here once published_.
84 : : */
85 : : // *****************************************************************************
86 : : #ifndef MixNumberFractionBeta_h
87 : : #define MixNumberFractionBeta_h
88 : :
89 : : #include <vector>
90 : : #include <cmath>
91 : :
92 : : #include "InitPolicy.hpp"
93 : : #include "MixNumberFractionBetaCoeffPolicy.hpp"
94 : : #include "RNG.hpp"
95 : : #include "Particles.hpp"
96 : :
97 : : namespace walker {
98 : :
99 : : extern ctr::InputDeck g_inputdeck;
100 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
101 : :
102 : : //! \brief MixNumberFractionBeta SDE used polymorphically with DiffEq
103 : : //! \details The template arguments specify policies and are used to configure
104 : : //! the behavior of the class. The policies are:
105 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
106 : : //! - Coefficients - coefficients policy, see
107 : : //! DiffEq/MixNumberFractionBetaCoeffPolicy.h
108 : : template< class Init, class Coefficients >
109 : : class MixNumberFractionBeta {
110 : :
111 : : private:
112 : : using ncomp_t = tk::ctr::ncomp_t;
113 : :
114 : : public:
115 : : //! \brief Constructor
116 : : //! \param[in] c Index specifying which system of mix number-fraction beta
117 : : //! SDEs to construct. There can be multiple mixnumfracbeta ... end blocks
118 : : //! in a control file. This index specifies which mix number-fraction beta
119 : : //! SDE system to instantiate. The index corresponds to the order in which
120 : : //! the mixnumfracbeta ... end blocks are given the control file.
121 : 0 : explicit MixNumberFractionBeta( ncomp_t c ) :
122 : : m_c( c ),
123 : : m_depvar(
124 : 0 : g_inputdeck.get< tag::param, tag::mixnumfracbeta, tag::depvar >().at(c)
125 : : ),
126 : : m_ncomp(
127 : 0 : g_inputdeck.get< tag::component >().get< tag::mixnumfracbeta >().at(c)/3
128 : : ),
129 : : m_offset(
130 : 0 : g_inputdeck.get< tag::component >().offset< tag::mixnumfracbeta >(c) ),
131 [ - - ]: 0 : m_rng( g_rng.at( tk::ctr::raw(
132 : 0 : g_inputdeck.get< tag::param, tag::mixnumfracbeta, tag::rng >().at(c) ) )
133 : : ),
134 : : m_bprime(),
135 : : m_S(),
136 : : m_kprime(),
137 : : m_rho2(),
138 : : m_rcomma(),
139 : : m_b(),
140 : : m_k(),
141 : : coeff(
142 : 0 : m_ncomp,
143 [ - - ]: 0 : g_inputdeck.get< tag::param, tag::mixnumfracbeta, tag::bprime >().at(c),
144 [ - - ]: 0 : g_inputdeck.get< tag::param, tag::mixnumfracbeta, tag::S >().at(c),
145 [ - - ]: 0 : g_inputdeck.get< tag::param, tag::mixnumfracbeta, tag::kappaprime >().at(c),
146 [ - - ]: 0 : g_inputdeck.get< tag::param, tag::mixnumfracbeta, tag::rho2 >().at(c),
147 [ - - ]: 0 : g_inputdeck.get< tag::param, tag::mixnumfracbeta, tag::rcomma >().at(c),
148 [ - - ]: 0 : m_bprime, m_S, m_kprime, m_rho2, m_rcomma, m_b, m_k ) {}
149 : :
150 : : //! Initalize SDE, prepare for time integration
151 : : //! \param[in] stream Thread (or more precisely stream) ID
152 : : //! \param[in,out] particles Array of particle properties
153 : 0 : void initialize( int stream, tk::Particles& particles ) {
154 : : //! Set initial conditions using initialization policy
155 : : Init::template
156 : : init< tag::mixnumfracbeta >
157 : 0 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
158 : 0 : }
159 : :
160 : : //! \brief Advance particles according to the system of mix number-fraction
161 : : //! beta SDEs
162 : : //! \param[in,out] particles Array of particle properties
163 : : //! \param[in] stream Thread (or more precisely stream) ID
164 : : //! \param[in] dt Time step size
165 : : //! \param[in] moments Map of statistical moments
166 : 0 : void advance( tk::Particles& particles,
167 : : int stream,
168 : : tk::real dt,
169 : : tk::real,
170 : : const std::map< tk::ctr::Product, tk::real >& moments )
171 : : {
172 : : // Update SDE coefficients
173 : 0 : coeff.update( m_depvar, m_ncomp, moments, m_bprime, m_kprime, m_b, m_k );
174 : : // Advance particles
175 : 0 : const auto npar = particles.nunk();
176 [ - - ]: 0 : for (auto p=decltype(npar){0}; p<npar; ++p) {
177 : : // Generate Gaussian random numbers with zero mean and unit variance
178 [ - - ]: 0 : std::vector< tk::real > dW( m_ncomp );
179 [ - - ]: 0 : m_rng.gaussian( stream, m_ncomp, dW.data() );
180 : : // Advance all m_ncomp scalars
181 [ - - ]: 0 : for (ncomp_t i=0; i<m_ncomp; ++i) {
182 [ - - ]: 0 : tk::real& X = particles( p, i, m_offset );
183 : 0 : tk::real d = m_k[i] * X * (1.0 - X) * dt;
184 [ - - ]: 0 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
185 : 0 : X += 0.5*m_b[i]*(m_S[i] - X)*dt + d*dW[i];
186 : : // Compute instantaneous values derived from updated X
187 [ - - ]: 0 : particles( p, m_ncomp+i, m_offset ) = rho( X, i );
188 [ - - ]: 0 : particles( p, m_ncomp*2+i, m_offset ) = vol( X, i );
189 : : }
190 : : }
191 : 0 : }
192 : :
193 : : private:
194 : : const ncomp_t m_c; //!< Equation system index
195 : : const char m_depvar; //!< Dependent variable
196 : : const ncomp_t m_ncomp; //!< Number of components
197 : : const ncomp_t m_offset; //!< Offset SDE operates from
198 : : const tk::RNG& m_rng; //!< Random number generator
199 : :
200 : : //! Coefficients
201 : : std::vector< kw::sde_bprime::info::expect::type > m_bprime;
202 : : std::vector< kw::sde_S::info::expect::type > m_S;
203 : : std::vector< kw::sde_kappaprime::info::expect::type > m_kprime;
204 : : std::vector< kw::sde_rho2::info::expect::type > m_rho2;
205 : : std::vector< kw::sde_rcomma::info::expect::type > m_rcomma;
206 : : std::vector< kw::sde_b::info::expect::type > m_b;
207 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
208 : :
209 : : //! Coefficients policy
210 : : Coefficients coeff;
211 : :
212 : : //! \brief Return density for mole fraction
213 : : //! \details Functional wrapper around the dependent variable of the beta
214 : : //! SDE. This function returns the instantaneous density, rho,
215 : : //! based on the number fraction, X, and parameters rho2 and r'.
216 : : //! \param[in] X Instantaneous value of the mole fraction, X
217 : : //! \param[in] i Index specifying which (of multiple) parameters to use
218 : : //! \return Instantaneous value of the density, rho
219 : 0 : tk::real rho( tk::real X, ncomp_t i ) const {
220 : 0 : return m_rho2[i] * ( 1.0 - m_rcomma[i] * X );
221 : : }
222 : :
223 : : //! \brief Return specific volume for mole fraction
224 : : //! \details Functional wrapper around the dependent variable of the beta
225 : : //! SDE. This function returns the instantaneous specific volume, V,
226 : : //! based on the number fraction, X, and parameters rho2 and r'.
227 : : //! \param[in] X Instantaneous value of the mole fraction, X
228 : : //! \param[in] i Index specifying which (of multiple) parameters to use
229 : : //! \return Instantaneous value of the specific volume, V
230 : 0 : tk::real vol( tk::real X, ncomp_t i ) const {
231 : 0 : return 1.0 / rho( X, i );
232 : : }
233 : : };
234 : :
235 : : } // walker::
236 : :
237 : : #endif // MixNumberFractionBeta_h
|
