Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Dissipation/Dissipation.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 A dissipation model for Lagrangian particles
9 : : \details This file implements the time integration of a system of
10 : : stochastic differential equations to model fluctuating
11 : : particle time scales (from which the turbulent kinetic energy dissipation
12 : : rate can be computed).
13 : : */
14 : : // *****************************************************************************
15 : : #ifndef Dissipation_h
16 : : #define Dissipation_h
17 : :
18 : : #include <array>
19 : : #include <vector>
20 : : #include <cmath>
21 : :
22 : : #include "InitPolicy.hpp"
23 : : #include "DissipationCoeffPolicy.hpp"
24 : : #include "RNG.hpp"
25 : : #include "Particles.hpp"
26 : : #include "CoupledEq.hpp"
27 : :
28 : : namespace walker {
29 : :
30 : : extern ctr::InputDeck g_inputdeck;
31 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
32 : :
33 : : //! \brief Dissipation equation used polymorphically with DiffEq
34 : : //! \details The template arguments specify policies and are used to configure
35 : : //! the behavior of the class. The policies are:
36 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
37 : : //! - Coefficients - coefficients policy, see DiffEq/DissipationCoeffPolicy.h
38 : : template< class Init, class Coefficients >
39 : : class Dissipation {
40 : :
41 : : private:
42 : : using ncomp_t = tk::ctr::ncomp_t;
43 : : using eq = tag::dissipation;
44 : :
45 : : public:
46 : : //! \brief Constructor
47 : : //! \param[in] c Index specifying which system of dissipation SDEs to
48 : : //! construct. There can be multiple dissipation ... end blocks in a
49 : : //! control file. This index specifies which dissipation SDE system to
50 : : //! instantiate. The index corresponds to the order in which the
51 : : //! dissipation ... end blocks are given the control file.
52 : 18 : explicit Dissipation( ncomp_t c ) :
53 : : m_c( c ),
54 : 18 : m_depvar( g_inputdeck.get< tag::param, eq, tag::depvar >().at(c) ),
55 : 18 : m_ncomp( g_inputdeck.get< tag::component >().get< eq >().at(c) ),
56 : 18 : m_offset( g_inputdeck.get< tag::component >().offset< eq >(c) ),
57 [ + - ]: 18 : m_rng( g_rng.at( tk::ctr::raw(
58 : 18 : g_inputdeck.get< tag::param, eq, tag::rng >().at(c) ) ) ),
59 : 18 : m_velocity_coupled( coupled< eq, tag::velocity >( c ) ),
60 : 18 : m_velocity_depvar( depvar< eq, tag::velocity >( c ) ),
61 : 36 : m_velocity_offset( offset< eq, tag::velocity, tag::velocity_id >( c ) ),
62 : : m_coeff(
63 : 18 : g_inputdeck.get< tag::param, eq, tag::c3 >().at(c),
64 : 18 : g_inputdeck.get< tag::param, eq, tag::c4 >().at(c),
65 : 18 : g_inputdeck.get< tag::param, eq, tag::com1 >().at(c),
66 : 18 : g_inputdeck.get< tag::param, eq, tag::com2 >().at(c),
67 : 18 : m_c3, m_c4, m_com1, m_com2 ),
68 : 18 : m_O( tk::ctr::mean( m_depvar, 0 ) ),
69 : 18 : m_R( {{ tk::ctr::variance( m_velocity_depvar, 0 ),
70 : 18 : tk::ctr::variance( m_velocity_depvar, 1 ),
71 : 18 : tk::ctr::variance( m_velocity_depvar, 2 ),
72 : 18 : tk::ctr::covariance( m_velocity_depvar, 0, m_velocity_depvar, 1 )
73 [ + - ][ + - ]: 162 : }} )
[ + - ][ + - ]
74 : : {
75 [ - + ][ - - ]: 18 : Assert( m_ncomp == 1, "Dissipation eq number of components must be 1" );
[ - - ][ - - ]
76 : 18 : }
77 : :
78 : : //! Initalize SDE, prepare for time integration
79 : : //! \param[in] stream Thread (or more precisely stream) ID
80 : : //! \param[in,out] particles Array of particle properties
81 : 90 : void initialize( int stream, tk::Particles& particles ) {
82 : : //! Set initial conditions using initialization policy
83 : : Init::template init< eq >
84 : 90 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
85 : 90 : }
86 : :
87 : : //! \brief Advance particles according to the dissipation SDE
88 : : //! \param[in,out] particles Array of particle properties
89 : : //! \param[in] stream Thread (or more precisely stream) ID
90 : : //! \param[in] dt Time step size
91 : : //! \param[in] moments Map of statistical moments
92 : 6750 : void advance( tk::Particles& particles,
93 : : int stream,
94 : : tk::real dt,
95 : : tk::real,
96 : : const std::map< tk::ctr::Product, tk::real >& moments )
97 : : {
98 : : using tk::ctr::lookup;
99 : :
100 : : // Access mean turbulence frequency
101 [ + - ]: 6750 : tk::real O = lookup( m_O, moments );
102 : :
103 : 6750 : tk::ctr::Term u( static_cast<char>(std::toupper(m_velocity_depvar)), 0,
104 : : tk::ctr::Moment::ORDINARY );
105 : 6750 : tk::ctr::Term v( static_cast<char>(std::toupper(m_velocity_depvar)), 1,
106 : : tk::ctr::Moment::ORDINARY );
107 : 6750 : tk::ctr::Term w( static_cast<char>(std::toupper(m_velocity_depvar)), 2,
108 : : tk::ctr::Moment::ORDINARY );
109 [ + - ]: 13500 : auto r11 = tk::ctr::Product( { u, u } );
110 [ + - ]: 13500 : auto r22 = tk::ctr::Product( { v, v } );
111 [ + - ]: 13500 : auto r33 = tk::ctr::Product( { w, w } );
112 [ + - ]: 13500 : auto r12 = tk::ctr::Product( { u, v } );
113 : :
114 : : // Compute turbulent kinetic energy
115 [ + - ]: 6750 : tk::real k = ( lookup(r11,moments) +
116 [ + - ]: 6750 : lookup(r22,moments) +
117 [ + - ]: 6750 : lookup(r33,moments) ) / 2.0;
118 : :
119 : : // Production of turbulent kinetic energy
120 : 6750 : tk::real S = 1.0; // prescribed shear: hard-coded in a single direction
121 [ + - ]: 6750 : tk::real P = -lookup(r12,moments)*S;
122 : :
123 : : // Source for turbulent frequency
124 : 6750 : tk::real Som = m_com2 - m_com1*P/(O*k);
125 : :
126 : : // Update source based on coefficients policy
127 : 6750 : Coefficients::src( Som );
128 : :
129 : 6750 : const auto npar = particles.nunk();
130 [ + + ]: 31506750 : for (auto p=decltype(npar){0}; p<npar; ++p) {
131 : : // Generate a Gaussian random number with zero mean and unit variance
132 : : tk::real dW;
133 [ + - ]: 31500000 : m_rng.gaussian( stream, m_ncomp, &dW );
134 : : // Advance particle frequency
135 [ + - ]: 31500000 : tk::real& Op = particles( p, 0, m_offset );
136 : 31500000 : tk::real d = 2.0*m_c3*m_c4*O*O*Op*dt;
137 [ + + ]: 31500000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
138 : 31500000 : Op += (-m_c3*(Op-O) - Som*Op)*O*dt + d*dW;
139 : : }
140 : 6750 : }
141 : :
142 : : private:
143 : : const ncomp_t m_c; //!< Equation system index
144 : : const char m_depvar; //!< Dependent variable
145 : : const ncomp_t m_ncomp; //!< Number of components
146 : : const ncomp_t m_offset; //!< Offset SDE operates from
147 : : const tk::RNG& m_rng; //!< Random number generator
148 : :
149 : : const bool m_velocity_coupled; //!< True if coupled to velocity
150 : : const char m_velocity_depvar; //!< Coupled velocity dependent variable
151 : : const ncomp_t m_velocity_offset; //!< Offset of coupled velocity eq
152 : :
153 : :
154 : : //! Coefficients policy
155 : : Coefficients m_coeff;
156 : :
157 : : // Model coefficients
158 : : tk::real m_c3;
159 : : tk::real m_c4;
160 : : tk::real m_com1;
161 : : tk::real m_com2;
162 : :
163 : : //! tk::Product used to access the mean of turbulence frequency
164 : : tk::ctr::Product m_O;
165 : : //! Array of tk::Product used to access the the Reynolds stress
166 : : std::array< tk::ctr::Product, 4 > m_R;
167 : : };
168 : :
169 : : } // walker::
170 : :
171 : : #endif // Dissipation_h
|
