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