Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Dirichlet/MixDirichlet.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 Mixture 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), with
12 : : various constraints to model multi-material mixing process in turbulent
13 : : flows.
14 : :
15 : : In a nutshell, the equation integrated governs a set of scalars,
16 : : \f$\color[HTML]{dcdcdc}0\!\le\!Y_\alpha\f$,
17 : : \f$\color[HTML]{dcdcdc}\alpha\!=\!1,\dots,N-1\f$,
18 : : \f$\color[HTML]{dcdcdc}\sum_{\alpha=1}^{N-1}Y_\alpha\!\le\!1\f$, as
19 : :
20 : : @m_class{m-show-m}
21 : :
22 : : \f[
23 : : \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha Y_N -
24 : : (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha
25 : : Y_N}\mathrm{d}W_\alpha(t), \qquad \alpha=1,\dots,N-1
26 : : \f]
27 : :
28 : : @m_class{m-hide-m}
29 : :
30 : : \f[
31 : : \begin{split}
32 : : \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha Y_N -
33 : : (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha
34 : : Y_N}\mathrm{d}W_\alpha(t), \\ \alpha=1,\dots,N-1
35 : : \end{split}
36 : : \f]
37 : :
38 : : with parameter vectors \f$\color[HTML]{dcdcdc}b_\alpha > 0\f$,
39 : : \f$\color[HTML]{dcdcdc}\kappa_\alpha > 0\f$, and \f$\color[HTML]{dcdcdc}0 <
40 : : S_\alpha < 1\f$, and \f$\color[HTML]{dcdcdc}Y_N =
41 : : 1-\sum_{\alpha=1}^{N-1}Y_\alpha\f$. Here
42 : : \f$\color[HTML]{dcdcdc}\mathrm{d}W_\alpha(t)\f$ is an isotropic vector-valued
43 : : [Wiener
44 : : process](http://en.wikipedia.org/wiki/Wiener_process) with independent
45 : : increments. The invariant distribution is the Dirichlet distribution,
46 : : provided the parameters of the drift and diffusion terms satisfy
47 : : \f[
48 : : (1-S_1) b_1 / \kappa_1 = \dots = (1-S_{N-1}) b_{N-1} / \kappa_{N-1}.
49 : : \f]
50 : : To keep the invariant distribution Dirichlet, the above constraint on the
51 : : coefficients must be satisfied. For more details on the Dirichlet SDE, see
52 : : https://doi.org/10.1155/2013/842981.
53 : : */
54 : : // *****************************************************************************
55 : : #ifndef MixDirichlet_h
56 : : #define MixDirichlet_h
57 : :
58 : : #include <vector>
59 : : #include <cmath>
60 : : #include <cfenv>
61 : :
62 : : #include "InitPolicy.hpp"
63 : : #include "MixDirichletCoeffPolicy.hpp"
64 : : #include "RNG.hpp"
65 : : #include "Particles.hpp"
66 : :
67 : : namespace walker {
68 : :
69 : : extern ctr::InputDeck g_inputdeck;
70 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
71 : :
72 : : //! \brief MixDirichlet SDE used polymorphically with DiffEq
73 : : //! \details The template arguments specify policies and are used to configure
74 : : //! the behavior of the class. The policies are:
75 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
76 : : //! - Coefficients - coefficients policy, see DiffEq/DirCoeffPolicy.h
77 : : template< class Init, class Coefficients >
78 : : class MixDirichlet {
79 : :
80 : : private:
81 : : using ncomp_t = tk::ctr::ncomp_t;
82 : : using eq = tag::mixdirichlet;
83 : :
84 : : public:
85 : : //! Number of derived variables computed by the MixDirichlet SDE
86 : : //! \details Derived variables: Nth scalar, density, specific volume.
87 : : static const std::size_t NUMDERIVED = 3;
88 : :
89 : : //! \brief Constructor
90 : : //! \param[in] c Index specifying which MixDirichlet SDE to construct. There
91 : : //! can be multiple dirichlet ... end blocks in a control file. This index
92 : : //! specifies which MixDirichlet SDE to instantiate. The index corresponds
93 : : //! to the order in which the dirichlet ... end blocks are given the
94 : : //! control file.
95 : 32 : explicit MixDirichlet( ncomp_t c ) :
96 : : m_c( c ),
97 : 32 : m_depvar( g_inputdeck.get< tag::param, eq, tag::depvar >().at(c) ),
98 : : // subtract the number of derived variables computed, see advance()
99 : 32 : m_ncomp( g_inputdeck.get< tag::component >().get< eq >().at(c) -
100 : : NUMDERIVED ),
101 : : m_offset(
102 : 32 : g_inputdeck.get< tag::component >().offset< eq >(c) ),
103 [ + - ]: 32 : m_rng( g_rng.at( tk::ctr::raw(
104 : 32 : g_inputdeck.get< tag::param, eq, tag::rng >().at(c) ) ) ),
105 : 32 : m_norm( g_inputdeck.get< tag::param, eq, tag::normalization >().at(c) ),
106 : : m_b(),
107 : : m_S(),
108 : : m_kprime(),
109 : : m_k(),
110 : : m_rho(),
111 : : m_r(),
112 : 32 : coeff( m_ncomp,
113 : 32 : m_norm,
114 [ + - ]: 32 : g_inputdeck.get< tag::param, eq, tag::b >().at(c),
115 [ + - ]: 32 : g_inputdeck.get< tag::param, eq, tag::S >().at(c),
116 [ + - ]: 32 : g_inputdeck.get< tag::param, eq, tag::kappaprime >().at(c),
117 [ + - ]: 32 : g_inputdeck.get< tag::param, eq, tag::rho >().at(c),
118 [ + - ]: 224 : m_b, m_S, m_kprime, m_rho, m_r, m_k ) {}
119 : :
120 : : //! Initalize SDE, prepare for time integration
121 : : //! \param[in] stream Thread (or more precisely stream) ID
122 : : //! \param[in,out] particles Array of particle properties
123 : 32 : void initialize( int stream, tk::Particles& particles ) {
124 : :
125 : : // Ensure the size of the parameter vector holding the pure-fluid
126 : : // densities is N = K+1 = m_ncomp+1
127 [ - + ][ - - ]: 32 : Assert( m_rho.size() == m_ncomp+1, "Size mismatch" );
[ - - ][ - - ]
128 : :
129 : : //! Set initial conditions using initialization policy
130 : 32 : Init::template init< eq >( g_inputdeck, m_rng, stream, particles, m_c,
131 [ + - ]: 32 : m_ncomp, m_offset );
132 : :
133 : : fenv_t fe;
134 : 32 : feholdexcept( &fe );
135 : :
136 : 32 : const auto npar = particles.nunk();
137 [ + + ]: 400032 : for (auto p=decltype(npar){0}; p<npar; ++p) {
138 : : // Violating boundedness here is a hard error as indicates a problem
139 : : // with the initial conditions
140 [ + + ]: 1600000 : for (ncomp_t i=0; i<m_ncomp+1; ++i) {
141 [ + - ]: 1200000 : auto y = particles( p, i, m_offset );
142 [ + - ][ - + ]: 1200000 : if (y < 0.0 || y > 1.0) Throw( "IC scalar out of bounds" );
[ - - ][ - - ]
[ - - ]
143 : : }
144 : : // Initialize derived instantaneous variables
145 [ + - ]: 400000 : derived( particles, p );
146 : : }
147 : :
148 : 32 : feclearexcept( FE_UNDERFLOW );
149 : 32 : feupdateenv( &fe );
150 : 32 : }
151 : :
152 : : //! \brief Advance particles according to the MixDirichlet SDE
153 : : //! \param[in,out] particles Array of particle properties
154 : : //! \param[in] stream Thread (or more precisely stream) ID
155 : : //! \param[in] dt Time step size
156 : : //! \param[in] moments Map of statistical moments
157 : 9568 : void advance( tk::Particles& particles,
158 : : int stream,
159 : : tk::real dt,
160 : : tk::real,
161 : : const std::map< tk::ctr::Product, tk::real >& moments )
162 : : {
163 : : // Update SDE coefficients
164 : 9568 : coeff.update( m_depvar, m_ncomp, m_norm, DENSITY_OFFSET, VOLUME_OFFSET,
165 [ + - ]: 9568 : moments, m_rho, m_r, m_kprime, m_b, m_k, m_S );
166 : :
167 : : fenv_t fe;
168 : 9568 : feholdexcept( &fe );
169 : :
170 : : // Advance particles
171 : 9568 : const auto npar = particles.nunk();
172 [ + + ]: 119609568 : for (auto p=decltype(npar){0}; p<npar; ++p) {
173 : : // Generate Gaussian random numbers with zero mean and unit variance
174 [ + - ]: 239200000 : std::vector< tk::real > dW( m_ncomp );
175 [ + - ]: 119600000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
176 : :
177 : : // Advance all m_ncomp (=N=K+1) scalars
178 [ + - ]: 119600000 : auto& yn = particles( p, m_ncomp, m_offset );
179 [ + + ]: 358800000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
180 [ + - ]: 239200000 : auto& y = particles( p, i, m_offset );
181 : 239200000 : tk::real d = m_k[i] * y * yn * dt;
182 [ + + ]: 239200000 : if (d < 0.0) d = 0.0;
183 : 239200000 : d = std::sqrt( d );
184 : 239200000 : auto dy = 0.5*m_b[i]*( m_S[i]*yn - (1.0-m_S[i])*y )*dt + d*dW[i];
185 : 239200000 : y += dy;
186 : 239200000 : yn -= dy;
187 : : }
188 : : // Compute derived instantaneous variables
189 [ + - ]: 119600000 : derived( particles, p );
190 : : }
191 : :
192 : 9568 : feclearexcept( FE_UNDERFLOW );
193 : 9568 : feupdateenv( &fe );
194 : 9568 : }
195 : :
196 : : private:
197 : : //! Offset of particle density in solution array relative to YN
198 : : static const std::size_t DENSITY_OFFSET = 1;
199 : : //! Offset of particle specific volume in solution array relative to YN
200 : : static const std::size_t VOLUME_OFFSET = 2;
201 : :
202 : : const ncomp_t m_c; //!< Equation system index
203 : : const char m_depvar; //!< Dependent variable
204 : : const ncomp_t m_ncomp; //!< Number of components, K = N-1
205 : : const ncomp_t m_offset; //!< Offset SDE operates from
206 : : const tk::RNG& m_rng; //!< Random number generator
207 : : const ctr::NormalizationType m_norm;//!< Normalization type
208 : :
209 : : //! Coefficients
210 : : std::vector< kw::sde_b::info::expect::type > m_b;
211 : : std::vector< kw::sde_S::info::expect::type > m_S;
212 : : std::vector< kw::sde_kappa::info::expect::type > m_kprime;
213 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
214 : : std::vector< kw::sde_rho::info::expect::type > m_rho;
215 : : std::vector< kw::sde_r::info::expect::type > m_r;
216 : :
217 : : //! Coefficients policy
218 : : Coefficients coeff;
219 : :
220 : : //! \brief Return density for mass fractions
221 : : //! \details This function returns the instantaneous density, rho,
222 : : //! based on the multiple mass fractions, Y_c.
223 : : //! \param[in] particles Array of particle properties
224 : : //! \param[in] p Particle index
225 : : //! \return Instantaneous value of the density, rho
226 : : //! \details This is computed based 1/rho = sum_{i=1}^N Y_i/R_i, where R_i
227 : : //! are the constant pure-fluid densities and the Y_i are the mass
228 : : //! fractions, of the N materials.
229 : 120000000 : tk::real rho( const tk::Particles& particles, ncomp_t p ) const {
230 : : // start computing density
231 : 120000000 : tk::real d = 0.0;
232 [ + + ]: 480000000 : for (ncomp_t i=0; i<m_ncomp+1; ++i)
233 : 360000000 : d += particles( p, i, m_offset ) / m_rho[i];
234 : : // return particle density
235 : 120000000 : return 1.0/d;
236 : : }
237 : :
238 : : //! Compute instantaneous values derived from particle mass fractions
239 : : //! \param[in,out] particles Particle properties array
240 : : //! \param[in] p Particle index
241 : 120000000 : void derived( tk::Particles& particles, ncomp_t p ) const {
242 : : // compute instantaneous fluid-density based on particle mass fractions
243 : 120000000 : auto density = rho( particles, p );
244 : : //// Compute and store instantaneous density
245 : 120000000 : particles( p, m_ncomp+DENSITY_OFFSET, m_offset ) = density;
246 : : // Store instantaneous specific volume
247 : 120000000 : particles( p, m_ncomp+VOLUME_OFFSET, m_offset ) = 1.0/density;
248 : 120000000 : }
249 : :
250 : : };
251 : :
252 : : } // walker::
253 : :
254 : : #endif // MixDirichlet_h
|
