Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Velocity/Langevin.cpp
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 Functionality implementing Langevin models for the velocity
9 : : \details Functionality implementing Langevin models for the velocity.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include <cctype>
14 : : #include <array>
15 : :
16 : : #include "StatCtr.hpp"
17 : : #include "Langevin.hpp"
18 : :
19 : : std::array< tk::real, 9 >
20 : 3375 : walker::slm( tk::real hts, tk::real C0 )
21 : : // *****************************************************************************
22 : : // Calculate the 2nd order tensor Gij based on the simplified Langevin model
23 : : //! \param[in] hts Inverse hydrodynamics time scale, e.g., eps/k
24 : : //! \param[in] C0 Coefficient C0 in SLM
25 : : //! \return Tensor Gij computed based on the simplified Langevin model
26 : : // *****************************************************************************
27 : : {
28 : : std::array< tk::real, 9 > G;
29 : : G.fill( 0.0 );
30 : 3375 : G[0] = G[4] = G[8] = -(0.5+0.75*C0) * hts;
31 : :
32 : 3375 : return G;
33 : : }
34 : :
35 : : std::array< tk::real, 9 >
36 : 3375 : walker::glm( tk::real hts,
37 : : tk::real C0,
38 : : const std::array< tk::real, 6 >& rs,
39 : : const std::array< tk::real, 9 >& dU )
40 : : // *****************************************************************************
41 : : // Calculate the 2nd order tensor Gij based on the generalized Langevin model
42 : : //! \param[in] hts Inverse hydrodynamics time scale, e.g., eps/k
43 : : //! \param[in] C0 Coefficient C0 in SLM
44 : : //! \param[in] rs Reynolds stress
45 : : //! \param[in] dU Mean velocity gradient
46 : : //! \return Tensor Gij computed based on the simplified Langevin model
47 : : // *****************************************************************************
48 : : {
49 : : // Generalized Langevion model coefficients
50 : 3375 : std::array< tk::real, 2 > ALPHA{{ -(0.5 + 0.75*C0), 3.7 }};
51 : : std::array< tk::real, 3 > BETA{{ -0.2, 0.8, -0.2 }};
52 : : std::array< tk::real, 6 > GAMMA{{ -1.28, 3.01, -2.18, 0.0, 4.29, -3.09 }};
53 : :
54 : : // Compute Reynolds stress anisotropy
55 : 3375 : tk::real tr = rs[0] + rs[1] + rs[2];
56 : 3375 : std::array< tk::real, 9 > b{{ rs[0]/tr-1.0/3.0,
57 : 3375 : rs[3]/tr,
58 : 3375 : rs[4]/tr,
59 : : rs[3]/tr,
60 : 3375 : rs[1]/tr-1.0/3.0,
61 : 3375 : rs[5]/tr,
62 : : rs[4]/tr,
63 : : rs[5]/tr,
64 : 3375 : rs[2]/tr-1.0/3.0 }};
65 : : // Compute Gij
66 : : std::array< tk::real, 9 > G;
67 : : G.fill( 0.0 );
68 [ + + ]: 13500 : for (std::size_t i=0; i<3; ++i ) {
69 : : // to main diagonal: hts * ALPHA1 * delta_ij + BETA1 * deltaij * d<Ul>/dxl
70 : 10125 : G[i*3+i] += hts*ALPHA[0] + BETA[0]*(dU[0] + dU[4] + dU[8]);
71 : : // to main diagonal: GAMMA1 * deltaij * bml * d<Um>/dxl
72 : : tk::real dtmp = 0.0;
73 [ + + ]: 40500 : for (std::size_t m=0; m<3; ++m)
74 [ + + ]: 121500 : for (std::size_t l=0; l<3; ++l )
75 : 91125 : dtmp += b[m*3+l]*dU[m*3+l];
76 : 10125 : G[i*3+i] += GAMMA[0]*dtmp;
77 : : // to main and off-diagonal
78 [ + + ]: 40500 : for (std::size_t j=0; j<3; ++j) {
79 : 30375 : G[i*3+j] += hts*ALPHA[1]*b[i*3+j] + // eps/k * ALPHA2 * bij
80 : 30375 : BETA[1]*dU[i*3+j] + // BETA2 * d<Ui>/dj
81 : 30375 : BETA[2]*dU[j*3+i] + // BETA3 * d<Uj>/di
82 : : // GAMMA4 * bij * d<Ul>/dxl
83 : 30375 : GAMMA[3]*b[i*3+j]*(dU[0] + dU[4] + dU[8]);
84 [ + + ]: 121500 : for (std::size_t l=0; l<3; ++l)
85 : 91125 : G[i*3+j] += GAMMA[1]*b[j*3+l]*dU[i*3+l] + // GAMMA2 * bjl * d<Ui>/dxl
86 : 91125 : GAMMA[2]*b[j*3+l]*dU[l*3+i] + // GAMMA3 * bjl * d<Ul>/dxi
87 : 91125 : GAMMA[4]*b[i*3+l]*dU[l*3+j] + // GAMMA5 * bil * d<Ul>/dxj
88 : 91125 : GAMMA[5]*b[i*3+l]*dU[j*3+l]; // GAMMA6 * bil * d<Uj>/dxl
89 : : }
90 : : }
91 : :
92 : 3375 : return G;
93 : : }
94 : :
95 : : std::array< tk::real, 6 >
96 [ - + ]: 6750 : walker::reynoldsStress( char depvar,
97 : : ctr::DepvarType solve,
98 : : const std::map< tk::ctr::Product, tk::real >& moments )
99 : : // *****************************************************************************
100 : : // Compute the Reynolds stress tensor
101 : : //! \param[in] depvar Dependent variable labeling a velocity eq
102 : : //! \param[in] solve Enum selecting what the velocity eq solved for
103 : : //! (fluctuating velocity of full/instantaneous veelocity)
104 : : //! \param[in] moments Map of statistical moments
105 : : //! \return Symmetric part of the Reynolds stress tensor
106 : : // *****************************************************************************
107 : : {
108 : : using tk::ctr::Product;
109 : :
110 : : // Extract diagonal of the Reynolds stress
111 : : Product r11, r22, r33, r12, r13, r23;
112 [ - + ]: 6750 : if (solve == ctr::DepvarType::FULLVAR) {
113 : :
114 : : using tk::ctr::variance;
115 : : using tk::ctr::covariance;
116 [ - - ]: 0 : r11 = variance( depvar, 0 );
117 [ - - ]: 0 : r22 = variance( depvar, 1 );
118 [ - - ]: 0 : r33 = variance( depvar, 2 );
119 [ - - ]: 0 : r12 = covariance( depvar, 0, depvar, 1 );
120 [ - - ]: 0 : r13 = covariance( depvar, 0, depvar, 2 );
121 [ - - ]: 0 : r23 = covariance( depvar, 1, depvar, 2 );
122 : :
123 [ + - ]: 6750 : } else if (solve == ctr::DepvarType::FLUCTUATION) {
124 : :
125 : : // Since we are solving for the fluctuating velocity, the "ordinary"
126 : : // moments, e.g., <U1U1>, are really central moments, i.e., <u1u1>.
127 : : using tk::ctr::Term;
128 : : using tk::ctr::Moment;
129 : 6750 : auto d = static_cast< char >( std::toupper( depvar ) );
130 : : Term u( d, 0, Moment::ORDINARY );
131 : : Term v( d, 1, Moment::ORDINARY );
132 : : Term w( d, 2, Moment::ORDINARY );
133 [ + - ]: 6750 : r11 = Product( { u, u } );
134 [ + - ]: 6750 : r22 = Product( { v, v } );
135 [ + - ]: 6750 : r33 = Product( { w, w } );
136 [ + - ]: 6750 : r12 = Product( { u, v } );
137 [ + - ]: 6750 : r13 = Product( { u, w } );
138 [ + - ][ - - ]: 13500 : r23 = Product( { v, w } );
139 : :
140 [ - - ][ - - ]: 0 : } else Throw( "Depvar type not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
141 : :
142 : : // Compute nonzero components of the Reynolds stress
143 [ + - ]: 6750 : return {{ lookup(r11,moments), lookup(r22,moments), lookup(r33,moments),
144 [ + - ][ + - ]: 13500 : lookup(r12,moments), lookup(r13,moments), lookup(r23,moments) }};
[ + - ][ + - ]
[ + - ][ + - ]
145 : : }
146 : :
147 : : tk::real
148 [ - - ]: 0 : walker::tke( char depvar,
149 : : ctr::DepvarType solve,
150 : : const std::map< tk::ctr::Product, tk::real >& moments )
151 : : // *****************************************************************************
152 : : // Compute the turbulent kinetic energy
153 : : //! \param[in] depvar Dependent variable labeling a velocity eq
154 : : //! \param[in] solve Enum selecting what the velocity eq solved for
155 : : //! (fluctuating velocity of full/instantaneous veelocity)
156 : : //! \param[in] moments Map of statistical moments
157 : : //! \return Turbulent kinetic energy
158 : : // *****************************************************************************
159 : : {
160 : : using tk::ctr::Product;
161 : :
162 : : // Extract diagonal of the Reynolds stress
163 : : Product r11, r22, r33;
164 [ - - ]: 0 : if (solve == ctr::DepvarType::FULLVAR) {
165 : :
166 : : using tk::ctr::variance;
167 : : using tk::ctr::covariance;
168 [ - - ]: 0 : r11 = variance( depvar, 0 );
169 [ - - ]: 0 : r22 = variance( depvar, 1 );
170 [ - - ]: 0 : r33 = variance( depvar, 2 );
171 : :
172 [ - - ]: 0 : } else if (solve == ctr::DepvarType::FLUCTUATION) {
173 : :
174 : : // Since we are solving for the fluctuating velocity, the "ordinary"
175 : : // moments, e.g., <U1U1>, are really central moments, i.e., <u1u1>.
176 : : using tk::ctr::Term;
177 : : using tk::ctr::Moment;
178 : 0 : auto d = static_cast< char >( std::toupper( depvar ) );
179 : : Term u( d, 0, Moment::ORDINARY );
180 : : Term v( d, 1, Moment::ORDINARY );
181 : : Term w( d, 2, Moment::ORDINARY );
182 [ - - ]: 0 : r11 = Product( { u, u } );
183 [ - - ]: 0 : r22 = Product( { v, v } );
184 [ - - ][ - - ]: 0 : r33 = Product( { w, w } );
185 : :
186 [ - - ][ - - ]: 0 : } else Throw( "Depvar type not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
187 : :
188 : : // Compute nonzero components of the Reynolds stress
189 [ - - ][ - - ]: 0 : return (lookup(r11,moments) + lookup(r22,moments) + lookup(r33,moments))/2.0;
[ - - ][ - - ]
190 : : }
|
