Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Beta/MixMassFractionBeta.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 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). 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., MixMassFractionBetaCoeffConst::update().
19 : :
20 : : - Second, there two additional random variables computed, the same as also
21 : : computed by the mass-fraction beta equation, see also
22 : : DiffEq/MassFractionBeta.h.
23 : :
24 : : In a nutshell, the equation integrated governs a set of scalars,
25 : : \f$\color[HTML]{dcdcdc}0\!\le\!Y_\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}Y_\alpha(t) = \frac{b_\alpha}{2}\left(S_\alpha - Y_\alpha\right)
32 : : \mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha(1-Y_\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}Y_\alpha(t) = \frac{b_\alpha}{2}\left(S_\alpha - Y_\alpha\right)
40 : : \mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha(1-Y_\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{Y} (1-\irmean{Y}) ]\f$. The fluctuation about the mean,
52 : : \f$\color[HTML]{dcdcdc} \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
53 : : \irmean{Y} \f$, is defined as usual: \f$\color[HTML]{dcdcdc}
54 : : \newcommand{\irmean}[1]{{\langle{#1}\rangle}} x = Y - \irmean{Y} \f$, and
55 : : \f$\color[HTML]{dcdcdc}b'\f$ and \f$\color[HTML]{dcdcdc} \kappa'\f$ are
56 : : user-specified constants. Also,
57 : : \f$\color[HTML]{dcdcdc}\mathrm{d}W_\alpha(t)\f$ 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 : : In addition to integrating the above SDE, there are two additional functions
64 : : of \f$\color[HTML]{dcdcdc} Y_\alpha \f$ are computed as
65 : : \f[ \begin{aligned}
66 : : \rho(Y_\alpha) & = \frac{ \rho_{2\alpha} }{ 1 + r_\alpha Y_\alpha } \\
67 : : V(Y_\alpha) & = \frac{1}{ \rho(Y_\alpha) }
68 : : \end{aligned} \f]
69 : : These equations compute the instantaneous mixture density,
70 : : \f$\color[HTML]{dcdcdc} \rho \f$, and instantaneous specific volume,
71 : : \f$\color[HTML]{dcdcdc} V_\alpha \f$, for equation \f$\color[HTML]{dcdcdc}
72 : : \alpha \f$ in the system. These quantities are used in binary mixing of
73 : : variable-density turbulence between two fluids with constant densities,
74 : : \f$\color[HTML]{dcdcdc} \rho_1, \f$ and \f$\color[HTML]{dcdcdc} \rho_2 \f$. The
75 : : additional parameters, \f$\color[HTML]{dcdcdc} \rho_2 \f$ and
76 : : \f$\color[HTML]{dcdcdc} r \f$ are user input parameters and kept constant
77 : : during integration. Since we compute the above variables,
78 : : \f$\color[HTML]{dcdcdc}\rho,\f$ and \f$\color[HTML]{dcdcdc}V\f$, and call them
79 : : mixture density and specific volume, respectively, \f$\color[HTML]{dcdcdc}Y\f$,
80 : : governed by the beta SDE is a mass fraction.
81 : :
82 : : _All of this is unpublished, but will be linked in here once published_.
83 : : */
84 : : // *****************************************************************************
85 : : #ifndef MixMassFractionBeta_h
86 : : #define MixMassFractionBeta_h
87 : :
88 : : #include <vector>
89 : : #include <tuple>
90 : :
91 : : #include "InitPolicy.hpp"
92 : : #include "MixMassFractionBetaCoeffPolicy.hpp"
93 : : #include "RNG.hpp"
94 : : #include "Particles.hpp"
95 : : #include "Table.hpp"
96 : : #include "CoupledEq.hpp"
97 : : #include "HydroTimeScales.hpp"
98 : : #include "HydroProductions.hpp"
99 : : #include "Walker/Options/HydroTimeScales.hpp"
100 : : #include "Walker/Options/HydroProductions.hpp"
101 : :
102 : : namespace walker {
103 : :
104 : : extern ctr::InputDeck g_inputdeck;
105 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
106 : :
107 : : //! \brief MixMassFractionBeta SDE used polymorphically with DiffEq
108 : : //! \details The template arguments specify policies and are used to configure
109 : : //! the behavior of the class. The policies are:
110 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
111 : : //! - Coefficients - coefficients policy, see
112 : : //! DiffEq/MixMassFractionBetaCoeffPolicy.h
113 : : template< class Init, class Coefficients >
114 : : class MixMassFractionBeta {
115 : :
116 : : private:
117 : : using ncomp_t = tk::ctr::ncomp_t;
118 : : using eq = tag::mixmassfracbeta;
119 : :
120 : : public:
121 : : //! Number of derived variables computed by the SDE
122 : : //! \details Derived variables: density, specific volume, 1 - mass fraction
123 : : //! \warning If you change this, you must also change the constant in
124 : : //! Velocity::m_mixmassfracbeta_ncomp in DiffEq/Velocity/Velocity.hpp.
125 : : //! To see where, search for NUMDERIVED there.
126 : : static const std::size_t NUMDERIVED = 3;
127 : :
128 : : //! \brief Constructor
129 : : //! \param[in] c Index specifying which system of mix mass-fraction beta
130 : : //! SDEs to construct. There can be multiple mixmassfracbeta ... end blocks
131 : : //! in a control file. This index specifies which mix mass-fraction beta
132 : : //! SDE system to instantiate. The index corresponds to the order in which
133 : : //! the mixmassfracbeta ... end blocks are given the control file.
134 : 15 : explicit MixMassFractionBeta( ncomp_t c ) :
135 : : m_c( c ),
136 : 15 : m_depvar( g_inputdeck.get< tag::param, eq, tag::depvar >().at(c) ),
137 : : // divide by the number of derived variables computed, see derived()
138 : 15 : m_ncomp( g_inputdeck.get< tag::component >().get< eq >().at(c) /
139 : : (NUMDERIVED + 1) ),
140 : 15 : m_offset( g_inputdeck.get< tag::component >().offset< eq >(c) ),
141 [ + - ]: 15 : m_rng( g_rng.at( tk::ctr::raw(
142 : 15 : g_inputdeck.get< tag::param, eq, tag::rng >().at(c) ) ) ),
143 : 15 : m_solve( g_inputdeck.get< tag::param, eq, tag::solve >().at(c) ),
144 : 15 : m_velocity_coupled( coupled< eq, tag::velocity >( c ) ),
145 : 15 : m_velocity_depvar( depvar< eq, tag::velocity >( c ) ),
146 : 30 : m_velocity_offset( offset< eq, tag::velocity, tag::velocity_id >( c ) ),
147 : : m_velocity_solve(
148 : 15 : m_velocity_coupled ?
149 : 0 : g_inputdeck.get< tag::param, tag::velocity, tag::solve >().at(
150 : : system_id< eq, tag::velocity, tag::velocity_id >( c ) ) :
151 : : ctr::DepvarType::FULLVAR ),
152 : 15 : m_dissipation_coupled( coupled< eq, tag::dissipation >( c ) ),
153 : 15 : m_dissipation_depvar( depvar< eq, tag::dissipation >( c ) ),
154 : : m_dissipation_offset(
155 : 30 : offset< eq, tag::dissipation, tag::dissipation_id >( c ) ),
156 : 15 : m_dY( initScalarGradient() ),
157 : : m_bprime(),
158 : : m_S(),
159 : : m_kprime(),
160 : : m_rho2(),
161 : : m_r(),
162 : : m_b(),
163 : : m_k(),
164 : : coeff(
165 : 15 : m_ncomp,
166 [ + - ]: 15 : g_inputdeck.get< tag::param, eq, tag::bprime >().at(c),
167 [ + - ]: 15 : g_inputdeck.get< tag::param, eq, tag::S >().at(c),
168 [ + - ]: 15 : g_inputdeck.get< tag::param, eq, tag::kappaprime >().at(c),
169 [ + - ]: 15 : g_inputdeck.get< tag::param, eq, tag::rho2 >().at(c),
170 [ + - ]: 15 : g_inputdeck.get< tag::param, eq, tag::r >().at(c),
171 [ - + ][ + - ]: 120 : m_bprime, m_S, m_kprime, m_rho2, m_r, m_b, m_k )
172 : : {
173 : : // Zero prescribed scalar gradient if full variable is solved for
174 [ + - ][ + - ]: 15 : if (m_solve == ctr::DepvarType::FULLVAR) m_dY.fill( 0.0 );
175 : : // Populate inverse hydrodynamics time scales and hydrodyanmics
176 : : // production/dissipation extracted from DNS
177 [ + + ]: 15 : if ( Coefficients::type() == ctr::CoeffPolicyType::HYDROTIMESCALE ) {
178 : : // Configure inverse hydrodyanmics time scale from DNS
179 : : const auto& hts =
180 [ + - ]: 4 : g_inputdeck.get< tag::param, eq, tag::hydrotimescales >().at(c);
181 [ + - ]: 8 : ctr::HydroTimeScales ot;
182 : : // cppcheck-suppress useStlAlgorithm
183 [ + + ][ + - ]: 24 : for (auto t : hts) m_hts.push_back( ot.table(t) );
[ + - ]
184 [ - + ][ - - ]: 4 : Assert( m_hts.size() == m_ncomp, "Number of inverse hydro time scale "
[ - - ][ - - ]
185 : : "tables associated does not match the components integrated" );
186 : :
187 : : // Configure hydrodyanmics production/dissipation from DNS
188 : : const auto& hp =
189 [ + - ]: 4 : g_inputdeck.get< tag::param, eq, tag::hydroproductions >().at(c);
190 [ + - ]: 8 : ctr::HydroProductions op;
191 : : // cppcheck-suppress useStlAlgorithm
192 [ + + ][ + - ]: 24 : for (auto t : hp) m_hp.push_back( op.table(t) );
[ + - ]
193 [ - + ][ - - ]: 4 : Assert( m_hp.size() == m_ncomp, "Number of hydro "
[ - - ][ - - ]
194 : : "production/dissipation tables associated does not match the "
195 : : "components integrated" );
196 : : }
197 : 15 : }
198 : :
199 : : //! Initalize SDE, prepare for time integration
200 : : //! \param[in] stream Thread (or more precisely stream) ID
201 : : //! \param[in,out] particles Array of particle properties
202 : 51 : void initialize( int stream, tk::Particles& particles ) {
203 : : //! Set initial conditions using initialization policy
204 : : Init::template init< eq >
205 : 51 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
206 : : // Initialize values derived from primary prognostic variable
207 : 51 : const auto npar = particles.nunk();
208 [ + + ]: 304051 : for (auto p=decltype(npar){0}; p<npar; ++p)
209 [ + + ]: 1824000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
210 [ - + ][ - - ]: 1520000 : Assert( particles( p, i, m_offset ) > 0.0, "Beta IC out of bounds!" );
[ - - ][ - - ]
211 [ - + ][ - - ]: 1520000 : Assert( particles( p, i, m_offset ) < 1.0, "Beta IC out of bounds!" );
[ - - ][ - - ]
212 : 1520000 : derived( particles, p, i );
213 : : }
214 : 51 : }
215 : :
216 : : //! \brief Advance particles according to the system of mix mass-fraction
217 : : //! beta SDEs
218 : : //! \param[in,out] particles Array of particle properties
219 : : //! \param[in] stream Thread (or more precisely stream) ID
220 : : //! \param[in] dt Time step size
221 : : //! \param[in] t Physical time of the simulation
222 : : //! \param[in] moments Map of statistical moments
223 : 23536 : void advance( tk::Particles& particles,
224 : : int stream,
225 : : tk::real dt,
226 : : tk::real t,
227 : : const std::map< tk::ctr::Product, tk::real >& moments )
228 : : {
229 : : // Update SDE coefficients
230 : 23536 : coeff.update( m_depvar, m_dissipation_depvar, m_velocity_depvar,
231 : 23536 : m_velocity_solve, m_solve, m_ncomp, moments, m_bprime,
232 : 23536 : m_kprime, m_rho2, m_r, m_hts, m_hp, m_b, m_k, m_S, t );
233 : :
234 : 23536 : const auto eps = std::numeric_limits< tk::real >::epsilon();
235 : :
236 : : // Advance particles
237 : 23536 : const auto npar = particles.nunk();
238 [ + + ]: 4723536 : for (auto p=decltype(npar){0}; p<npar; ++p) {
239 : : // Generate Gaussian random numbers with zero mean and unit variance
240 [ + - ]: 9400000 : std::vector< tk::real > dW( m_ncomp );
241 [ + - ]: 4700000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
242 : :
243 : : // Access coupled particle velocity
244 : 4700000 : tk::real u = 0.0, v = 0.0, w = 0.0;
245 : : using std::abs;
246 [ + - ][ + - ]: 4700000 : if (abs(m_dY[0]) > eps || abs(m_dY[1]) > eps || abs(m_dY[2]) > eps) {
[ - + ][ - + ]
247 [ - - ]: 0 : u = particles( p, 0, m_velocity_offset );
248 [ - - ]: 0 : v = particles( p, 1, m_velocity_offset );
249 [ - - ]: 0 : w = particles( p, 2, m_velocity_offset );
250 : : }
251 : :
252 : : // Advance all m_ncomp scalars
253 [ + + ]: 28200000 : for (ncomp_t i=0; i<m_ncomp; ++i) {
254 [ + - ]: 23500000 : tk::real& Y = particles( p, i, m_offset );
255 : 23500000 : tk::real d = m_k[i] * Y * (1.0 - Y) * dt;
256 [ + + ]: 23500000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
257 : 23500000 : Y += 0.5*m_b[i]*(m_S[i] - Y)*dt + d*dW[i]
258 : 23500000 : - (m_dY[0]*u - m_dY[1]*v - m_dY[2]*w)*dt;
259 : : // Compute instantaneous values derived from updated Y
260 [ + - ]: 23500000 : derived( particles, p, i );
261 : : }
262 : : }
263 : 23536 : }
264 : :
265 : : private:
266 : : const ncomp_t m_c; //!< Equation system index
267 : : const char m_depvar; //!< Dependent variable
268 : : const ncomp_t m_ncomp; //!< Number of components
269 : : const ncomp_t m_offset; //!< Offset SDE operates from
270 : : const tk::RNG& m_rng; //!< Random number generator
271 : : const ctr::DepvarType m_solve; //!< Depndent variable to solve for
272 : :
273 : : const bool m_velocity_coupled; //!< True if coupled to velocity
274 : : const char m_velocity_depvar; //!< Coupled velocity dependent variable
275 : : const ncomp_t m_velocity_offset; //!< Offset of coupled velocity eq
276 : : //! Quantity the coupled velocity eq solves for
277 : : const ctr::DepvarType m_velocity_solve;
278 : :
279 : : const bool m_dissipation_coupled; //!< True if coupled to dissipation
280 : : const char m_dissipation_depvar; //!< Depvar of coupled dissipation eq
281 : : const ncomp_t m_dissipation_offset; //!< Offset of coupled dissipation eq
282 : :
283 : : std::array< tk::real, 3 > m_dY; //! Prescribed mean scalar gradient
284 : :
285 : : //! Coefficients
286 : : std::vector< kw::sde_bprime::info::expect::type > m_bprime;
287 : : std::vector< kw::sde_S::info::expect::type > m_S;
288 : : std::vector< kw::sde_kappaprime::info::expect::type > m_kprime;
289 : : std::vector< kw::sde_rho2::info::expect::type > m_rho2;
290 : : std::vector< kw::sde_r::info::expect::type > m_r;
291 : : std::vector< kw::sde_b::info::expect::type > m_b;
292 : : std::vector< kw::sde_kappa::info::expect::type > m_k;
293 : :
294 : : //! Coefficients policy
295 : : Coefficients coeff;
296 : :
297 : : //! Selected inverse hydrodynamics time scales (if used) for each component
298 : : //! \details This is only used if the coefficients policy is
299 : : //! MixMassFracBetaCoeffHydroTimeScale. See constructor.
300 : : std::vector< tk::Table<1> > m_hts;
301 : :
302 : : //! Selected hydrodynamics production/dissipation (if used) for each comp.
303 : : //! \details This is only used if the coefficients policy is
304 : : //! MixMassFracBetaCoeffHydroTimeScale. See constructor.
305 : : std::vector< tk::Table<1> > m_hp;
306 : :
307 : : //! \brief Return density for mass fraction
308 : : //! \details Functional wrapper around the dependent variable of the beta
309 : : //! SDE. This function returns the instantaneous density, rho,
310 : : //! based on the mass fraction, Y, and parameters rho2 and r'.
311 : : //! \param[in] Y Instantaneous value of the mass fraction, Y
312 : : //! \param[in] i Index specifying which (of multiple) parameters to use
313 : : //! \return Instantaneous value of the density, rho
314 : 25020000 : tk::real rho( tk::real Y, ncomp_t i ) const {
315 : 25020000 : return m_rho2[i] / ( 1.0 + m_r[i] * Y );
316 : : }
317 : :
318 : : //! \brief Return specific volume for mass fraction
319 : : //! \details Functional wrapper around the dependent variable of the beta
320 : : //! SDE. This function returns the instantaneous specific volume, V,
321 : : //! based on the mass fraction, Y, and parameters rho2 and r'.
322 : : //! \param[in] Y Instantaneous value of the mass fraction, Y
323 : : //! \param[in] i Index specifying which (of multiple) parameters to use
324 : : //! \return Instantaneous value of the specific volume, V
325 : 25020000 : tk::real vol( tk::real Y, ncomp_t i ) const {
326 : 25020000 : return ( 1.0 + m_r[i] * Y ) / m_rho2[i];
327 : : }
328 : :
329 : : //! Compute instantaneous values derived from updated Y
330 : : //! \param[in,out] particles Particle properties array
331 : : //! \param[in] p Particle index
332 : : //! \param[in] i Component index
333 : 25020000 : void derived( tk::Particles& particles, ncomp_t p, ncomp_t i ) const {
334 : 25020000 : tk::real& Y = particles( p, i, m_offset );
335 : 25020000 : particles( p, m_ncomp+i, m_offset ) = rho( Y, i );
336 : 25020000 : particles( p, m_ncomp*2+i, m_offset ) = vol( Y, i );
337 : 25020000 : particles( p, m_ncomp*3+i, m_offset ) = 1.0 - Y;
338 : 25020000 : }
339 : :
340 : : //! Initialize imposed mean scalar gradient from user input
341 : 15 : std::array< tk::real, 3 > initScalarGradient() {
342 : 15 : const auto& mg = g_inputdeck.get< tag::param, eq, tag::mean_gradient >();
343 : 15 : std::array< tk::real, 3 > dY{{ 0.0, 0.0, 0.0 }};
344 [ - + ]: 15 : if (mg.size() > m_c) {
345 : 0 : const auto& g = mg[ m_c ];
346 : 0 : dY = {{ g[0], g[1], g[2] }};
347 : : }
348 : : Assert( dY.size() == 3, "Mean scalar gradient vector size must be 3" );
349 : 15 : return dY;
350 : : }
351 : : };
352 : :
353 : : } // walker::
354 : :
355 : : #endif // MixMassFractionBeta_h
|
