Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Velocity/Velocity.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 model for velocity in variable-density turbulence
9 : : \details This file implements the time integration of a system of stochastic
10 : : differential equations (SDEs) to model the fluctuating velocity components
11 : : in homogeneous variable-density turbulence. This model is an extension of
12 : : the generalized Langevin (GLM) model for constant-density flows by Haworth &
13 : : Pope (https://doi.org/10.1063/1.865723). The extension is roughly along the
14 : : lines of https://doi.org/10.1080/14685248.2011.554419.
15 : : */
16 : : // *****************************************************************************
17 : : #ifndef Velocity_h
18 : : #define Velocity_h
19 : :
20 : : #include <array>
21 : : #include <vector>
22 : : #include <cmath>
23 : :
24 : : #include "InitPolicy.hpp"
25 : : #include "VelocityCoeffPolicy.hpp"
26 : : #include "RNG.hpp"
27 : : #include "Particles.hpp"
28 : : #include "CoupledEq.hpp"
29 : :
30 : : namespace walker {
31 : :
32 : : extern ctr::InputDeck g_inputdeck;
33 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
34 : :
35 : : //! \brief Velocity SDE used polymorphically with DiffEq
36 : : //! \details The template arguments specify policies and are used to configure
37 : : //! the behavior of the class. The policies are:
38 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
39 : : //! - Coefficients - coefficients policy, see DiffEq/VelocityCoeffPolicy.h
40 : : template< class Init, class Coefficients >
41 : : class Velocity {
42 : :
43 : : private:
44 : : using ncomp_t = tk::ctr::ncomp_t;
45 : : using eq = tag::velocity;
46 : :
47 : : public:
48 : : //! \brief Constructor
49 : : //! \param[in] c Index specifying which system of velocity SDEs to construct
50 : : //! There can be multiple velocity ... end blocks in a control file. This
51 : : //! index specifies which velocity SDE system to instantiate. The index
52 : : //! corresponds to the order in which the velocity ... end blocks are
53 : : //! given the control file.
54 : 22 : explicit Velocity( ncomp_t c ) :
55 : : m_c( c ),
56 : : m_depvar( g_inputdeck.get< tag::param, eq, tag::depvar >().at(c) ),
57 : : m_solve( g_inputdeck.get< tag::param, eq, tag::solve >().at(c) ),
58 : : m_mixmassfracbeta_coupled( coupled< eq, tag::mixmassfracbeta >( c ) ),
59 : : m_mixmassfracbeta_depvar( depvar< eq, tag::mixmassfracbeta >( c ) ),
60 : : m_mixmassfracbeta_offset(
61 : 22 : offset< eq, tag::mixmassfracbeta, tag::mixmassfracbeta_id >( c ) ),
62 : : m_mixmassfracbeta_ncomp(
63 : : // The magic number, 4, below is MixMassFractionBeta::NUMDERIVED + 1,
64 : : // but cannot be given as such, because that would lead to circular
65 : : // dependencies of Velocity depending on MixMassfractionBeta, and vice
66 : : // versa.
67 [ - + ]: 22 : m_mixmassfracbeta_coupled ?
68 : 0 : ncomp< eq, tag::mixmassfracbeta, tag::mixmassfracbeta_id >( c ) / 4 :
69 : : 0 ),
70 : : m_numderived( numderived() ),
71 : 22 : m_ncomp( g_inputdeck.get< tag::component, eq >().at(c) - m_numderived ),
72 : : m_offset(
73 : 22 : g_inputdeck.get< tag::component >().offset< eq >(c) ),
74 : 22 : m_rng( g_rng.at( tk::ctr::raw(
75 : : g_inputdeck.get< tag::param, eq, tag::rng >().at(c) ) ) ),
76 : : m_position_coupled( coupled< eq, tag::position >( c ) ),
77 : : m_position_depvar( depvar< eq, tag::position >( c ) ),
78 [ + + ]: 22 : m_position_offset( offset< eq, tag::position, tag::position_id >( c ) ),
79 : : m_dissipation_coupled( coupled< eq, tag::dissipation >( c ) ),
80 : : m_dissipation_depvar( depvar< eq, tag::dissipation >( c ) ),
81 : : m_dissipation_offset(
82 : 22 : offset< eq, tag::dissipation, tag::dissipation_id >( c ) ),
83 : 22 : m_U( {{ tk::ctr::mean( m_depvar, 0 ),
84 : 22 : tk::ctr::mean( m_depvar, 1 ),
85 : 22 : tk::ctr::mean( m_depvar, 2 ) }} ),
86 : : m_variant( g_inputdeck.get< tag::param, eq, tag::variant >().at(c) ),
87 : : m_c0(),
88 : : m_G(),
89 [ - + ]: 22 : m_coeff( g_inputdeck.get< tag::param, eq, tag::c0 >().at(c), m_c0, m_dU ),
90 [ - + ][ - + ]: 118 : m_gravity( { 0.0, 0.0, 0.0 } )
[ + - ][ - + ]
[ - + ][ - + ]
[ + + ][ + + ]
[ - + ][ + - ]
91 : : {
92 : : Assert( m_ncomp == 3, "Velocity eq number of components must be 3" );
93 : : // Zero prescribed mean velocity gradient if full variable is solved for
94 [ + + ]: 22 : if (m_solve == ctr::DepvarType::FULLVAR) m_dU.fill( 0.0 );
95 : : // Populate inverse hydrodynamics time scales extracted from DNS
96 : : if ( Coefficients::type() == ctr::CoeffPolicyType::HYDROTIMESCALE ) {
97 : : // Configure inverse hydrodyanmics time scale from DNS
98 : : const auto& hts =
99 : : g_inputdeck.get< tag::param, eq, tag::hydrotimescales >().at(c);
100 : : Assert( hts.size() == 1,
101 : : "Velocity eq Hydrotimescales vector size must be 1" );
102 [ - - ][ - - ]: 0 : m_hts = ctr::HydroTimeScales().table( hts[0] );
[ - - ][ - - ]
103 : : }
104 : : // Initialize gravity body force if configured
105 : : const auto& gravity =
106 : : g_inputdeck.get< tag::param, eq, tag::gravity >().at(c);
107 [ + - ]: 22 : if (!gravity.empty()) {
108 : 22 : m_gravity[0] = gravity[0];
109 : 22 : m_gravity[1] = gravity[1];
110 : 22 : m_gravity[2] = gravity[2];
111 : : }
112 : 22 : }
113 : :
114 : : //! Compute number of derived variables
115 : : //! \return Number of derived variables computed
116 : : std::size_t numderived() const {
117 [ - - ][ - - ]: 22 : if (m_solve == ctr::DepvarType::PRODUCT ||
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - + ][ - + ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
118 : : m_solve == ctr::DepvarType::FLUCTUATING_MOMENTUM)
119 : : { // derived: 3 velocity components for each coupled mass fraction
120 : 0 : return m_mixmassfracbeta_ncomp * 3;
121 : : } else {
122 : : return 0;
123 : : }
124 : : }
125 : :
126 : : //! Initalize SDE, prepare for time integration
127 : : //! \param[in] stream Thread (or more precisely stream) ID
128 : : //! \param[in,out] particles Array of particle properties
129 : : void initialize( int stream, tk::Particles& particles ) {
130 : : // Set initial conditions using initialization policy
131 : : Init::template init< eq >
132 : 94 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
133 : : }
134 : :
135 : : //! \brief Advance particles according to the system of velocity SDEs
136 : : //! \param[in,out] particles Array of particle properties
137 : : //! \param[in] stream Thread (or more precisely stream) ID
138 : : //! \param[in] dt Time step size
139 : : //! \param[in] t Physical time of the simulation
140 : : //! \param[in] moments Map of statistical moments
141 : 7146 : void advance( tk::Particles& particles,
142 : : int stream,
143 : : tk::real dt,
144 : : tk::real t,
145 : : const std::map< tk::ctr::Product, tk::real >& moments )
146 : : {
147 : : using ctr::DepvarType;
148 : : const auto epsilon = std::numeric_limits< tk::real >::epsilon();
149 : :
150 : : // Update coefficients
151 : 7146 : tk::real eps = 0.0;
152 : 7146 : m_coeff.update( m_depvar, m_dissipation_depvar, moments, m_hts, m_solve,
153 : 7146 : m_variant, m_c0, t, eps, m_G );
154 : :
155 : : // Access mean velocity (if needed)
156 : : std::array< tk::real, 3 > U{{ 0.0, 0.0, 0.0 }};
157 [ + + ]: 7146 : if (m_solve == DepvarType::FULLVAR || m_solve == DepvarType::PRODUCT) {
158 : : using tk::ctr::lookup;
159 : 396 : U[0] = lookup( m_U[0], moments );
160 : 396 : U[1] = lookup( m_U[1], moments );
161 : 396 : U[2] = lookup( m_U[2], moments );
162 : : }
163 : :
164 : : // Modify G with the mean velocity gradient
165 [ + + ]: 71460 : for (std::size_t i=0; i<9; ++i) m_G[i] -= m_dU[i];
166 : :
167 : : // Access mean specific volume (if needed)
168 : : using tk::ctr::mean;
169 : : using tk::ctr::Term;
170 : : using tk::ctr::Product;
171 [ - + ]: 7146 : auto mixncomp = m_mixmassfracbeta_ncomp;
172 : : std::vector< tk::real > R;
173 : : std::vector< tk::real > RU;
174 : 7146 : auto Uc = static_cast< char >( std::toupper(m_depvar) );
175 [ - + ]: 7146 : if (m_solve == ctr::DepvarType::PRODUCT ||
176 : : m_solve == DepvarType::FLUCTUATING_MOMENTUM)
177 : : {
178 [ - - ]: 0 : R.resize( mixncomp, 0.0 );
179 [ - - ]: 0 : RU.resize( mixncomp*3, 0.0 );
180 [ - - ]: 0 : for (std::size_t c=0; c<mixncomp; ++c) {
181 [ - - ][ - - ]: 0 : R[c] = lookup(mean(m_mixmassfracbeta_depvar, c+mixncomp), moments);
[ - - ][ - - ]
182 : 0 : Term Rs( static_cast<char>(std::toupper(m_mixmassfracbeta_depvar)),
183 : : mixncomp + c,
184 : : tk::ctr::Moment::ORDINARY );
185 : 0 : std::array< Term, 3 > Us{
186 : 0 : Term( Uc, m_ncomp+(c*3)+0, tk::ctr::Moment::ORDINARY ),
187 : : Term( Uc, m_ncomp+(c*3)+1, tk::ctr::Moment::ORDINARY ),
188 : : Term( Uc, m_ncomp+(c*3)+2, tk::ctr::Moment::ORDINARY ) };
189 [ - - ][ - - ]: 0 : std::array< Product, 3 > RsUs{ Product( { Us[0], Rs } ),
[ - - ][ - - ]
[ - - ]
190 : : Product( { Us[1], Rs } ),
191 : : Product( { Us[2], Rs } ) };
192 [ - - ]: 0 : RU[ c*3+0 ] = lookup( RsUs[0], moments );
193 [ - - ][ - - ]: 0 : RU[ c*3+1 ] = lookup( RsUs[1], moments );
194 [ - - ]: 0 : RU[ c*3+2 ] = lookup( RsUs[2], moments );
195 : : }
196 : : }
197 : :
198 : : const auto npar = particles.nunk();
199 [ + + ]: 36457146 : for (auto p=decltype(npar){0}; p<npar; ++p) {
200 : : // Generate Gaussian random numbers with zero mean and unit variance
201 [ + - ]: 36450000 : std::vector< tk::real > dW( m_ncomp );
202 [ + - ]: 36450000 : m_rng.gaussian( stream, m_ncomp, dW.data() );
203 : : // Access particle velocity
204 [ + - ]: 36450000 : tk::real& Up = particles( p, 0, m_offset );
205 : : tk::real& Vp = particles( p, 1, m_offset );
206 : : tk::real& Wp = particles( p, 2, m_offset );
207 : : // Compute diffusion
208 : 36450000 : tk::real d = m_c0 * eps * dt;
209 [ + - ]: 36450000 : d = (d > 0.0 ? std::sqrt(d) : 0.0);
210 : : // Compute velocity fluctuation
211 : 36450000 : tk::real u = Up - U[0];
212 : 36450000 : tk::real v = Vp - U[1];
213 : 36450000 : tk::real w = Wp - U[2];
214 : : // Update particle velocity based on Langevin model
215 [ - + ]: 36450000 : Up += (m_G[0]*u + m_G[1]*v + m_G[2]*w)*dt + d*dW[0];
216 : 36450000 : Vp += (m_G[3]*u + m_G[4]*v + m_G[5]*w)*dt + d*dW[1];
217 : 36450000 : Wp += (m_G[6]*u + m_G[7]*v + m_G[8]*w)*dt + d*dW[2];
218 : : // Add gravity
219 [ - + ]: 36450000 : if (m_solve == ctr::DepvarType::PRODUCT ||
220 : : m_solve == ctr::DepvarType::FLUCTUATING_MOMENTUM)
221 : : {
222 [ - - ]: 0 : for (ncomp_t c=0; c<mixncomp; ++c) {
223 [ - - ]: 0 : auto rhoi = particles( p, mixncomp + c, m_mixmassfracbeta_offset );
224 [ - - ]: 0 : if (std::abs(rhoi) > epsilon) {
225 : : // add gravity force to particle momentum
226 : 0 : Up += (rhoi - R[c]) * m_gravity[0] * dt;
227 : 0 : Vp += (rhoi - R[c]) * m_gravity[1] * dt;
228 : 0 : Wp += (rhoi - R[c]) * m_gravity[2] * dt;
229 : : // compute derived particle velocity
230 : 0 : particles( p, m_ncomp+(c*3)+0, m_offset ) = (Up + RU[c*3+0])/rhoi;
231 : 0 : particles( p, m_ncomp+(c*3)+1, m_offset ) = (Vp + RU[c*3+1])/rhoi;
232 : 0 : particles( p, m_ncomp+(c*3)+2, m_offset ) = (Wp + RU[c*3+2])/rhoi;
233 : : }
234 : : }
235 : : } else {
236 : 36450000 : Up += m_gravity[0] * dt;
237 : 36450000 : Vp += m_gravity[1] * dt;
238 : 36450000 : Wp += m_gravity[2] * dt;
239 : : }
240 : : }
241 : 7146 : }
242 : :
243 : : private:
244 : : const ncomp_t m_c; //!< Equation system index
245 : : const char m_depvar; //!< Dependent variable
246 : : const ctr::DepvarType m_solve; //!< Dependent variable to solve for
247 : :
248 : : //! True if coupled to mixmassfracbeta
249 : : const bool m_mixmassfracbeta_coupled;
250 : : //! Depvar of coupled mixmassfracbeta eq
251 : : const char m_mixmassfracbeta_depvar;
252 : : //! Offset of coupled mixmassfracbeta eq
253 : : const ncomp_t m_mixmassfracbeta_offset;
254 : : //! Number of scalar components in coupled mixmassfracbeta eq
255 : : const ncomp_t m_mixmassfracbeta_ncomp;
256 : :
257 : : const ncomp_t m_numderived; //!< Number of derived variables
258 : : const ncomp_t m_ncomp; //!< Number of components
259 : : const ncomp_t m_offset; //!< Offset SDE operates from
260 : : const tk::RNG& m_rng; //!< Random number generator
261 : :
262 : : const bool m_position_coupled; //!< True if coupled to position
263 : : const char m_position_depvar; //!< Coupled position dependent variable
264 : : const ncomp_t m_position_offset; //!< Offset of coupled position eq
265 : :
266 : : const bool m_dissipation_coupled; //!< True if coupled to dissipation
267 : : const char m_dissipation_depvar; //!< Coupled dissipation dependent var
268 : : const ncomp_t m_dissipation_offset; //!< Offset of coupled dissipation eq
269 : :
270 : : //! Array of tk::ctr::Product used to access the mean velocity
271 : : const std::array< tk::ctr::Product, 3 > m_U;
272 : : //! Velocity model variant
273 : : const ctr::VelocityVariantType m_variant;
274 : :
275 : : //! Selected inverse hydrodynamics time scale (if used)
276 : : //! \details This is only used if the coefficients policy is
277 : : //! VelocityCoeffHydroTimeScale. See constructor.
278 : : tk::Table<1> m_hts;
279 : :
280 : : //! Coefficients
281 : : kw::sde_c0::info::expect::type m_c0;
282 : : std::array< tk::real, 9 > m_G;
283 : :
284 : : //! Coefficients policy
285 : : Coefficients m_coeff;
286 : :
287 : : //! (Optionally) prescribed mean velocity gradient
288 : : std::array< tk::real, 9 > m_dU;
289 : :
290 : : //! Optional gravity body force
291 : : std::array< tk::real, 3 > m_gravity;
292 : : };
293 : :
294 : : } // walker::
295 : :
296 : : #endif // Velocity_h
|
