Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Dirichlet/MixDirichletCoeffPolicy.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 Mixture Dirichlet coefficients policies
9 : : \details This file defines coefficients policy classes for the mixture
10 : : Dirichlet SDE, defined in DiffEq/Dirichlet/MixDirichlet.h.
11 : : For general requirements on the mixture Dirichlet SDE coefficients
12 : : policy classes see the header file.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <iostream>
17 : :
18 : : #include "MixDirichletCoeffPolicy.hpp"
19 : :
20 : : std::vector< kw::sde_r::info::expect::type >
21 : 40 : walker::MixDir_r( const std::vector< kw::sde_rho::info::expect::type >& rho,
22 : : ctr::NormalizationType norm )
23 : : // *****************************************************************************
24 : : // Compute parameter vector r based on r_i = rho_N/rho_i - 1
25 : : //! \param[in] rho Parameter vector rho to MixDirichlet
26 : : //! \param[in] norm Normalization type (N=heavy or N=light)
27 : : //! \return Parameter vector r, determined by parameter vector rho
28 : : // *****************************************************************************
29 : : {
30 [ - + ][ - - ]: 40 : Assert( rho.size() > 1, "Parameter vector rho must not be empty" );
[ - - ][ - - ]
31 : :
32 [ + - ]: 40 : std::vector< kw::sde_r::info::expect::type > r( rho.size()-1 );
33 : :
34 [ + + ]: 120 : for (std::size_t i=0; i<rho.size()-1; ++i) {
35 [ + + ]: 80 : if (norm == ctr::NormalizationType::LIGHT) // rhoN = rhoL
36 : 50 : r[i] = rho.back()/rho[i] + 1.0;
37 : : else // rhoN = rhoH
38 : 30 : r[i] = rho.back()/rho[i] - 1.0;
39 : : }
40 : :
41 : : //if (norm == ctr::NormalizationType::LIGHT)
42 : : // r.push_back( 2.0 );
43 : :
44 : 40 : return r;
45 : : }
46 : :
47 : 8 : walker::MixDirichletCoeffConst::MixDirichletCoeffConst(
48 : : ncomp_t ncomp,
49 : : ctr::NormalizationType norm,
50 : : const std::vector< kw::sde_b::info::expect::type >& b_,
51 : : const std::vector< kw::sde_S::info::expect::type >& S_,
52 : : const std::vector< kw::sde_kappa::info::expect::type >& kprime_,
53 : : const std::vector< kw::sde_rho::info::expect::type >& rho_,
54 : : std::vector< kw::sde_b::info::expect::type >& b,
55 : : std::vector< kw::sde_S::info::expect::type >& S,
56 : : std::vector< kw::sde_kappa::info::expect::type >& kprime,
57 : : std::vector< kw::sde_rho::info::expect::type >& rho,
58 : : std::vector< kw::sde_r::info::expect::type >& r,
59 : : std::vector< kw::sde_kappa::info::expect::type >& k )
60 : : // *****************************************************************************
61 : : // Constructor: initialize coefficients
62 : : //! \param[in] ncomp Number of scalar components in this SDE system
63 : : //! \param[in] norm Normalization type (N=heavy or N=light)
64 : : //! \param[in] b_ Vector used to initialize coefficient vector b
65 : : //! \param[in] S_ Vector used to initialize coefficient vector S
66 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime and k
67 : : //! \param[in] rho_ Vector used to initialize coefficient vector rho and r
68 : : //! \param[in,out] b Coefficient vector to be initialized
69 : : //! \param[in,out] S Coefficient vector to be initialized
70 : : //! \param[in,out] kprime Coefficient vector to be initialized
71 : : //! \param[in,out] rho Coefficient vector to be initialized
72 : : //! \param[in,out] r Coefficient vector to be initialized
73 : : //! \param[in,out] k Coefficient vector to be initialized
74 : : // *****************************************************************************
75 : : {
76 [ - + ][ - - ]: 8 : ErrChk( b_.size() == ncomp,
[ - - ][ - - ]
77 : : "Wrong number of MixDirichlet SDE parameters 'b'");
78 [ - + ][ - - ]: 8 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
79 : : "Wrong number of MixDirichlet SDE parameters 'S'");
80 [ - + ][ - - ]: 8 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
81 : : "Wrong number of MixDirichlet SDE parameters 'kappaprime'");
82 [ - + ][ - - ]: 8 : ErrChk( rho_.size() == ncomp+1,
[ - - ][ - - ]
83 : : "Wrong number of MixDirichlet SDE parameters 'rho'");
84 : :
85 : 8 : b = b_;
86 : 8 : S = S_;
87 : 8 : kprime = kprime_;
88 : 8 : rho = rho_;
89 [ + - ]: 8 : k.resize( kprime.size(), 0.0 );
90 : :
91 : : // Compute parameter vector r based on r_i = rho_N/rho_i - 1
92 [ - + ][ - - ]: 8 : Assert( r.empty(), "Parameter vector r must be empty" );
[ - - ][ - - ]
93 : 8 : r = MixDir_r( rho, norm );
94 : 8 : }
95 : :
96 : : void
97 : 2392 : walker::MixDirichletCoeffConst::update(
98 : : char /*depvar*/,
99 : : ncomp_t ncomp,
100 : : ctr::NormalizationType /*norm*/,
101 : : std::size_t /* density_offset */,
102 : : std::size_t /* volume_offset */,
103 : : const std::map< tk::ctr::Product, tk::real >& /*moments*/,
104 : : const std::vector< kw::sde_rho::info::expect::type >& /*rho*/,
105 : : const std::vector< kw::sde_r::info::expect::type >& /*r*/,
106 : : const std::vector< kw::sde_kappa::info::expect::type >& kprime,
107 : : const std::vector< kw::sde_b::info::expect::type >& /*b*/,
108 : : std::vector< kw::sde_kappa::info::expect::type >& k,
109 : : std::vector< kw::sde_kappa::info::expect::type >& S ) const
110 : : // *****************************************************************************
111 : : // Update coefficients
112 : : //! \param[in] depvar Dependent variable
113 : : //! \param[in] ncomp Number of scalar components in this SDE system
114 : : //! \param[in] norm Normalization type (N=heavy or N=light)
115 : : //! \param[in] density_offset Offset of particle density in solution array
116 : : //! relative to YN
117 : : //! \param[in] volume_offset Offset of particle specific volume in solution
118 : : //! array relative to YN
119 : : //! \param[in] moments Map of statistical moments estimated
120 : : //! \param[in] rho Coefficient vector
121 : : //! \param[in] r Coefficient Vector
122 : : //! \param[in] kprime Coefficient vector
123 : : //! \param[in] b Coefficient vector
124 : : //! \param[in,out] k Coefficient vector to be updated
125 : : //! \param[in,out] S Coefficient vector to be updated
126 : : // *****************************************************************************
127 : : {
128 : : using tk::ctr::lookup;
129 : : using tk::ctr::mean;
130 : : using tk::ctr::variance;
131 : : using tk::ctr::Term;
132 : : using tk::ctr::Moment;
133 : : using tk::ctr::Product;
134 : :
135 [ + + ]: 7176 : for (ncomp_t c=0; c<ncomp; ++c) {
136 : 4784 : k[c] = kprime[c];
137 : : }
138 : :
139 [ + + ]: 7176 : for (ncomp_t c=0; c<ncomp; ++c) {
140 [ + - ][ - + ]: 4784 : if (S[c] < 0.0 || S[c] > 1.0) {
[ - + ]
141 : 0 : std::cout << "S[" << c << "] bounds violated: " << S[c] << '\n';
142 : : }
143 : : }
144 : 2392 : }
145 : :
146 : 24 : walker::MixDirichletHomogeneous::MixDirichletHomogeneous(
147 : : ncomp_t ncomp,
148 : : ctr::NormalizationType norm,
149 : : const std::vector< kw::sde_b::info::expect::type >& b_,
150 : : const std::vector< kw::sde_S::info::expect::type >& S_,
151 : : const std::vector< kw::sde_kappa::info::expect::type >& kprime_,
152 : : const std::vector< kw::sde_rho::info::expect::type >& rho_,
153 : : std::vector< kw::sde_b::info::expect::type >& b,
154 : : std::vector< kw::sde_S::info::expect::type >& S,
155 : : std::vector< kw::sde_kappa::info::expect::type >& kprime,
156 : : std::vector< kw::sde_rho::info::expect::type >& rho,
157 : : std::vector< kw::sde_r::info::expect::type >& r,
158 : : std::vector< kw::sde_kappa::info::expect::type >& k )
159 : : // *****************************************************************************
160 : : // Constructor: initialize coefficients
161 : : //! \param[in] ncomp Number of scalar components in this SDE system
162 : : //! \param[in] norm Normalization type (N=heavy or N=light)
163 : : //! \param[in] b_ Vector used to initialize coefficient vector b
164 : : //! \param[in] S_ Vector used to initialize coefficient vector S
165 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime and k
166 : : //! \param[in] rho_ Vector used to initialize coefficient vector rho and r
167 : : //! \param[in,out] b Coefficient vector to be initialized
168 : : //! \param[in,out] S Coefficient vector to be initialized
169 : : //! \param[in,out] kprime Coefficient vector to be initialized
170 : : //! \param[in,out] rho Coefficient vector to be initialized
171 : : //! \param[in,out] r Coefficient vector to be initialized
172 : : //! \param[in,out] k Coefficient vector to be initialized
173 : : // *****************************************************************************
174 : : {
175 [ - + ][ - - ]: 24 : ErrChk( b_.size() == ncomp,
[ - - ][ - - ]
176 : : "Wrong number of MixDirichlet SDE parameters 'b'");
177 [ - + ][ - - ]: 24 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
178 : : "Wrong number of MixDirichlet SDE parameters 'S'");
179 [ - + ][ - - ]: 24 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
180 : : "Wrong number of MixDirichlet SDE parameters 'kappaprime'");
181 [ - + ][ - - ]: 24 : ErrChk( rho_.size() == ncomp+1,
[ - - ][ - - ]
182 : : "Wrong number of MixDirichlet SDE parameters 'rho'");
183 : :
184 : 24 : b = b_;
185 : 24 : S = S_;
186 : 24 : kprime = kprime_;
187 : 24 : rho = rho_;
188 [ + - ]: 24 : k.resize( kprime.size(), 0.0 );
189 : :
190 : : // Compute parameter vector r based on r_i = rho_N/rho_i - 1
191 [ - + ][ - - ]: 24 : Assert( r.empty(), "Parameter vector r must be empty" );
[ - - ][ - - ]
192 : 24 : r = MixDir_r( rho, norm );
193 : 24 : }
194 : :
195 : : void
196 : 7176 : walker::MixDirichletHomogeneous::update(
197 : : char depvar,
198 : : ncomp_t ncomp,
199 : : ctr::NormalizationType norm,
200 : : std::size_t density_offset,
201 : : std::size_t /*volume_offset*/,
202 : : const std::map< tk::ctr::Product, tk::real >& moments,
203 : : const std::vector< kw::sde_rho::info::expect::type >& rho,
204 : : const std::vector< kw::sde_r::info::expect::type >& r,
205 : : const std::vector< kw::sde_kappa::info::expect::type >& kprime,
206 : : const std::vector< kw::sde_b::info::expect::type >& b,
207 : : std::vector< kw::sde_kappa::info::expect::type >& k,
208 : : std::vector< kw::sde_kappa::info::expect::type >& S ) const
209 : : // *****************************************************************************
210 : : // Update coefficients
211 : : //! \param[in] depvar Dependent variable
212 : : //! \param[in] ncomp Number of scalar components in this SDE system
213 : : //! \param[in] norm Normalization type (N=heavy or N=light)
214 : : //! \param[in] density_offset Offset of particle density in solution array
215 : : //! relative to YN
216 : : //! \param[in] volume_offset Offset of particle specific volume in solution
217 : : //! array relative to YN
218 : : //! \param[in] moments Map of statistical moments estimated
219 : : //! \param[in] rho Coefficient vector
220 : : //! \param[in] r Coefficient Vector
221 : : //! \param[in] kprime Coefficient vector
222 : : //! \param[in] b Coefficient vector
223 : : //! \param[in,out] k Coefficient vector to be updated
224 : : //! \param[in,out] S Coefficient vector to be updated
225 : : // *****************************************************************************
226 : : {
227 : : using tk::ctr::lookup;
228 : : using tk::ctr::mean;
229 : : using tk::ctr::Term;
230 : : using tk::ctr::Moment;
231 : : using tk::ctr::Product;
232 : :
233 : : // Shorthands for dependent variable, Y, used to construct statistics
234 : 7176 : auto Yv = static_cast< char >( std::toupper(depvar) );
235 : :
236 : 7176 : Term tR( Yv, ncomp+density_offset, Moment::ORDINARY );
237 : 7176 : Term tYN( Yv, ncomp, Moment::ORDINARY );
238 : :
239 [ + - ][ + - ]: 7176 : auto R2YN = lookup( Product({tR,tR,tYN}), moments );
240 : :
241 [ + - ]: 14352 : std::vector< tk::real > R2Y( ncomp, 0.0 );
242 [ + - ]: 14352 : std::vector< tk::real > R3YNY( ncomp, 0.0 );
243 [ + + ]: 21528 : for (ncomp_t c=0; c<ncomp; ++c) {
244 : 14352 : Term tYc( Yv, c, Moment::ORDINARY );
245 [ + - ][ + - ]: 14352 : R2Y[c] = lookup( Product({tR,tR,tYc}), moments ); // <R^2Yc>
246 [ + - ][ + - ]: 14352 : R3YNY[c] = lookup( Product({tR,tR,tR,tYc,tYN}), moments ); // <R^3YNYc>
247 : : }
248 : :
249 : : // Assume heavy-fluid normalization by default: rhoN = rhoH
250 : 7176 : tk::real rhoL = rho[0], rhoH = rho[ncomp];
251 : : // Overwrite if light-fluid normalization is configured
252 [ + + ]: 7176 : if (norm == ctr::NormalizationType::LIGHT) { // rhoN = rhoL
253 : 3588 : rhoL = rho[ncomp];
254 : 3588 : rhoH = rho[0];
255 : : }
256 : :
257 [ + + ]: 21528 : for (ncomp_t c=0; c<ncomp; ++c) {
258 : 14352 : k[c] = kprime[c];
259 : 14352 : auto rcp = rhoL/rho[c] + 1.0;
260 : 14352 : auto rc = r[c]; // heavy-fluid normalization by default
261 : : // Overwrite if light-fluid normalization is configured
262 [ + + ]: 14352 : if (norm == ctr::NormalizationType::LIGHT) rc = (rcp - 2.0) * rhoH/rhoL;
263 : 14352 : S[c] = (R2Y[c] + 2.0*k[c]/b[c]*rc/rhoH*R3YNY[c]) / (R2Y[c] + R2YN);
264 : : }
265 : :
266 [ + + ]: 21528 : for (ncomp_t c=0; c<ncomp; ++c) {
267 [ + - ][ - + ]: 14352 : if (S[c] < 0.0 || S[c] > 1.0) {
[ - + ]
268 [ - - ][ - - ]: 0 : std::cout << "S[" << c << "] bounds violated: " << S[c] << '\n';
[ - - ][ - - ]
[ - - ]
269 : : }
270 : : }
271 : 7176 : }
272 : :
273 : 0 : walker::MixDirichletHydroTimeScale::MixDirichletHydroTimeScale(
274 : : tk::ctr::ncomp_t ncomp,
275 : : ctr::NormalizationType norm,
276 : : const std::vector< kw::sde_b::info::expect::type >& b_,
277 : : const std::vector< kw::sde_S::info::expect::type >& S_,
278 : : const std::vector< kw::sde_kappa::info::expect::type >& kprime_,
279 : : const std::vector< kw::sde_rho::info::expect::type >& rho_,
280 : : std::vector< kw::sde_b::info::expect::type >& b,
281 : : std::vector< kw::sde_S::info::expect::type >& S,
282 : : std::vector< kw::sde_kappa::info::expect::type >& kprime,
283 : : std::vector< kw::sde_rho::info::expect::type >& rho,
284 : : std::vector< kw::sde_r::info::expect::type >& r,
285 : : std::vector< kw::sde_kappa::info::expect::type >& k )
286 : : // *****************************************************************************
287 : : // Constructor: initialize coefficients
288 : : //! \param[in] ncomp Number of scalar components in this SDE system
289 : : //! \param[in] norm Normalization type (N=heavy or N=light)
290 : : //! \param[in] b_ Vector used to initialize coefficient vector b
291 : : //! \param[in] S_ Vector used to initialize coefficient vector S
292 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime and k
293 : : //! \param[in] rho_ Vector used to initialize coefficient vector rho and r
294 : : //! \param[in,out] b Coefficient vector to be initialized
295 : : //! \param[in,out] S Coefficient vector to be initialized
296 : : //! \param[in,out] kprime Coefficient vector to be initialized
297 : : //! \param[in,out] rho Coefficient vector to be initialized
298 : : //! \param[in,out] r Coefficient vector to be initialized
299 : : //! \param[in,out] k Coefficient vector to be initialized
300 : : // *****************************************************************************
301 : : {
302 [ - - ][ - - ]: 0 : ErrChk( b_.size() == ncomp,
[ - - ][ - - ]
303 : : "Wrong number of MixDirichlet SDE parameters 'b'");
304 [ - - ][ - - ]: 0 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
305 : : "Wrong number of MixDirichlet SDE parameters 'S'");
306 [ - - ][ - - ]: 0 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
307 : : "Wrong number of MixDirichlet SDE parameters 'kappaprime'");
308 [ - - ][ - - ]: 0 : ErrChk( rho_.size() == ncomp+1,
[ - - ][ - - ]
309 : : "Wrong number of MixDirichlet SDE parameters 'rho'");
310 : :
311 : 0 : b = b_;
312 : 0 : S = S_;
313 : 0 : kprime = kprime_;
314 : 0 : rho = rho_;
315 [ - - ]: 0 : k.resize( kprime.size(), 0.0 );
316 : :
317 : : // Compute parameter vector r based on r_i = rho_N/rho_i - 1
318 [ - - ][ - - ]: 0 : Assert( r.empty(), "Parameter vector r must be empty" );
[ - - ][ - - ]
319 : 0 : r = MixDir_r( rho, norm );
320 : 0 : }
321 : :
322 : : void
323 : 0 : walker::MixDirichletHydroTimeScale::update(
324 : : char depvar,
325 : : ncomp_t ncomp,
326 : : ctr::NormalizationType norm,
327 : : std::size_t density_offset,
328 : : std::size_t volume_offset,
329 : : const std::map< tk::ctr::Product, tk::real >& moments,
330 : : const std::vector< kw::sde_rho::info::expect::type >& rho,
331 : : const std::vector< kw::sde_r::info::expect::type >& r,
332 : : const std::vector< kw::sde_kappa::info::expect::type >& kprime,
333 : : const std::vector< kw::sde_b::info::expect::type >& b,
334 : : std::vector< kw::sde_kappa::info::expect::type >& k,
335 : : std::vector< kw::sde_kappa::info::expect::type >& S ) const
336 : : // *****************************************************************************
337 : : // Update coefficients
338 : : //! \param[in] depvar Dependent variable
339 : : //! \param[in] ncomp Number of scalar components in this SDE system
340 : : //! \param[in] norm Normalization type (N=heavy or N=light)
341 : : //! \param[in] density_offset Offset of particle density in solution array
342 : : //! relative to YN
343 : : //! \param[in] volume_offset Offset of particle specific volume in solution
344 : : //! array relative to YN
345 : : //! \param[in] moments Map of statistical moments estimated
346 : : //! \param[in] rho Coefficient vector
347 : : //! \param[in] r Coefficient Vector
348 : : //! \param[in] kprime Coefficient vector
349 : : //! \param[in] b Coefficient vector
350 : : //! \param[in,out] k Coefficient vector to be updated
351 : : //! \param[in,out] S Coefficient vector to be updated
352 : : // *****************************************************************************
353 : : {
354 : : using tk::ctr::lookup;
355 : : using tk::ctr::mean;
356 : : using tk::ctr::variance;
357 : : using tk::ctr::Term;
358 : : using tk::ctr::Moment;
359 : : using tk::ctr::Product;
360 : :
361 : : // Shorthands for dependent variables, Y and y = Y - <Y>, used to construct
362 : : // statistics
363 : 0 : auto Yv = static_cast< char >( std::toupper(depvar) );
364 : 0 : auto yv = static_cast< char >( std::tolower(depvar) );
365 : :
366 : : // <R>
367 [ - - ][ - - ]: 0 : tk::real R = lookup( mean(depvar,ncomp+density_offset), moments );
368 : : //if (R < 1.0e-8) R = 1.0;
369 : :
370 : : // b = -<rv>, density-specific-volume covariance
371 : 0 : Term rhoprime( yv, ncomp+density_offset, Moment::CENTRAL );
372 : 0 : Term vprime( yv, ncomp+volume_offset, Moment::CENTRAL );
373 : : //auto ds = -lookup( Product({rhoprime,vprime}), moments );
374 : :
375 : : // b. = -<ry.>/<R>
376 [ - - ]: 0 : std::vector< tk::real > bc( ncomp, 0.0 );
377 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
378 : 0 : Term tr( yv, ncomp+density_offset, Moment::CENTRAL );
379 : 0 : Term ty( yv, c, Moment::CENTRAL );
380 [ - - ][ - - ]: 0 : bc[c] = -lookup( Product({tr,ty}), moments ) / R; // -<ryc>/<R>
381 : : }
382 : :
383 : 0 : Term tR( Yv, ncomp+density_offset, Moment::ORDINARY );
384 : 0 : Term tYN( Yv, ncomp, Moment::ORDINARY );
385 : :
386 [ - - ][ - - ]: 0 : auto R2YN = lookup( Product({tR,tR,tYN}), moments );
387 : :
388 [ - - ]: 0 : std::vector< tk::real > y2( ncomp, 0.0 );
389 [ - - ]: 0 : std::vector< tk::real > RY( ncomp, 0.0 );
390 [ - - ]: 0 : std::vector< tk::real > R2Y( ncomp, 0.0 );
391 [ - - ]: 0 : std::vector< tk::real > R3YNY( ncomp, 0.0 );
392 [ - - ]: 0 : std::vector< tk::real > R3Y2( ncomp*ncomp, 0.0 );
393 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
394 [ - - ][ - - ]: 0 : y2[c] = lookup( variance(yv,c), moments );
395 : 0 : Term tYc( Yv, c, Moment::ORDINARY );
396 [ - - ][ - - ]: 0 : RY[c] = lookup( Product({tR,tYc}), moments ); // <RYc>
397 [ - - ][ - - ]: 0 : R2Y[c] = lookup( Product({tR,tR,tYc}), moments ); // <R^2Yc>
398 [ - - ][ - - ]: 0 : R3YNY[c] = lookup( Product({tR,tR,tR,tYc,tYN}), moments ); // <R^3YNYc>
399 [ - - ]: 0 : for (ncomp_t d=0; d<ncomp; ++d) {
400 : 0 : Term tYd( Yv, d, Moment::ORDINARY );
401 : : // <R^3YcYd>
402 [ - - ][ - - ]: 0 : R3Y2[c*ncomp+d] = lookup( Product({tR,tR,tR,tYc,tYd}), moments );
403 : : }
404 : : //std::cout << "R2Y: " << R2Y[c] << ' ';
405 : : }
406 : : //std::cout << std::endl;
407 : :
408 : : // Reynolds means
409 : :
410 : : // Reynolds means, Yc
411 : : //std::vector< tk::real > Y( ncomp, 0.0 );
412 : : //for (ncomp_t c=0; c<ncomp; ++c) {
413 : : // Y[c] = lookup( mean(depvar,c), moments );
414 : : // //std::cout << "Y: " << Y[c] << ' ';
415 : : //}
416 : : //std::cout << std::endl;
417 : :
418 : : // sum of Yc
419 : : //tk::real sumY = 0.0;
420 : : //for (ncomp_t c=0; c<ncomp; ++c) sumY += Y[c];
421 : :
422 : : //// Y|Kc
423 : : //std::vector< tk::real > YK( ncomp, 0.0 );
424 : : //for (ncomp_t c=0; c<ncomp; ++c) {
425 : : // YK[c] = sumY - lookup( mean(depvar,c), moments );
426 : : // //std::cout << "YK: " << YK[c] << ' ';
427 : : //}
428 : : //std::cout << std::endl;
429 : :
430 : : // Favre means
431 : :
432 : : // Ytc
433 : : //std::vector< tk::real > Yt( ncomp, 0.0 );
434 : : //for (ncomp_t c=0; c<ncomp; ++c) {
435 : : // Yt[c] = RY[c] / R;
436 : : // //std::cout << "Yt: " << Yt[c] << ' ';
437 : : //}
438 : : //std::cout << std::endl;
439 : :
440 : : // sum of Ytc
441 : : //tk::real sumYt = 0.0;
442 : : //for (ncomp_t c=0; c<ncomp; ++c) sumYt += Yt[c];
443 : : //std::cout << "sumYt: " << sumYt << '\n';
444 : :
445 : : // Yt|Kc
446 : : //std::vector< tk::real > YtK( ncomp, 0.0 );
447 : : //for (ncomp_t c=0; c<ncomp; ++c) {
448 : : // YtK[c] = sumYt - Yt[c];
449 : : // //std::cout << "YtK: " << YtK[c] << ' ';
450 : : //}
451 : : //std::cout << std::endl;
452 : :
453 : : // sum of <R^2Yc>
454 : : //tk::real sumR2Y = 0.0;
455 : : //std::vector< tk::real > sumR3Y2( ncomp, 0.0 );
456 : : //for (ncomp_t c=0; c<ncomp; ++c) {
457 : : // sumR2Y += R2Y[c];
458 : : // for (ncomp_t d=0; d<ncomp; ++d) sumR3Y2[c] += R3Y2[c*ncomp+d];
459 : : //}
460 : : //std::cout << "sumR2Y: " << sumR2Y << std::endl;
461 : :
462 : : // <r^2>, density variance
463 : : //auto rhovar = lookup( variance(yv,ncomp), moments );
464 : : //std::cout << "<r^2>: " << rhovar << std::endl;
465 : :
466 : : // <R^2>
467 : : //auto R2 = lookup( Product({tR,tR}), moments );
468 : :
469 : : // Assume heavy-fluid normalization by default: rhoN = rhoH
470 : 0 : tk::real rhoL = rho[0], rhoH = rho[ncomp];
471 : : // Overwrite if light-fluid normalization is configured
472 [ - - ]: 0 : if (norm == ctr::NormalizationType::LIGHT) { // rhoN = rhoL
473 : 0 : rhoL = rho[ncomp];
474 : 0 : rhoH = rho[0];
475 : : }
476 : :
477 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
478 : :
479 : : //k[c] = kprime[c] * bc[c];
480 : : //k[c] = kprime[c] * ds;
481 : 0 : k[c] = kprime[c];
482 : : //k[c] = kprime[c] * y2[c];
483 : : //if (k[c] < 0.0)
484 : : // std::cout << "Positivity of k[" << c << "] violated: "
485 : : // << k[c] << '\n';
486 : :
487 : : //S[c] = 1.0/(1.0-YK[c]) - (1.0-Yt[c])/(1.0-YtK[c]);
488 : : //S[c] = YK[c]/(1.0-YK[c]) - (1.0-Yt[c])*YtK[c]/(1.0-YtK[c]) + Yt[c];
489 : : //S[c] = Yt[c] / ( 1.0 - sumYt + Yt[c] );
490 : : //S[c] = ( -2.0*(r[c]/rho[ncomp]*R2Y[c])*(1.0-sumYt) +
491 : : // (r[c]/rho[ncomp]*(rhovar-sumR2Y))*Yt[c] ) /
492 : : // ( -2.0*(r[c]/rho[ncomp]*R2Y[c])*(1.0-sumYt) -
493 : : // (1.0-sumYt-Yt[c])*(r[c]/rho[ncomp]*(rhovar-sumR2Y)) );
494 : :
495 : : // correlation of density gradient wrt Y_alpha and Y_alpha (Y_alpha = Yc)
496 : : //tk::real drYcYc = -r[c]/rho[ncomp]*R2Y[c];
497 : : //tk::real drYcYN = -r[c]/rho[ncomp]*(R2-sumR2Y);
498 : : //tk::real drYc2YcYN = 2.0*std::pow(r[c]/rho[ncomp],2.0)*(R3Y[c]-sumR3Y2[c]);
499 : : //S[c] = (drYcYc - k[c]/b[c]*drYc2YcYN) / (drYcYN + drYcYc);
500 : : //S[c] = Yt[c] / (1.0 - sumYt + Yt[c]) - k[c]/b[c]*drYc2YcYN / (drYcYN + drYcYc);
501 : :
502 : : //S[c] = Yt[c] / (1.0 - sumYt + Yt[c]); // S_infty
503 : :
504 : 0 : auto rcp = rhoL/rho[c] + 1.0;
505 : 0 : auto rc = r[c]; // heavy-fluid normalization by default
506 : : // Overwrite if light-fluid normalization is configured
507 [ - - ]: 0 : if (norm == ctr::NormalizationType::LIGHT) rc = (rcp - 2.0) * rhoH/rhoL;
508 : 0 : S[c] = (R2Y[c] + 2.0*k[c]/b[c]*rc/rhoH*R3YNY[c]) / (R2Y[c] + R2YN);
509 : :
510 : : //std::cout << "S[" << c << "] = " << S[c] << ", b/k(1-S) = "
511 : : // << b[c]/k[c]*(1.0-S[c]) << '\n';
512 : :
513 : : }
514 : : //std::cout << std::endl;
515 : :
516 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
517 [ - - ][ - - ]: 0 : if (S[c] < 0.0 || S[c] > 1.0) {
[ - - ]
518 [ - - ][ - - ]: 0 : std::cout << "S[" << c << "] bounds violated: " << S[c] << '\n';
[ - - ][ - - ]
[ - - ]
519 : : //S[c] = 0.5;
520 : : }
521 : : }
522 : : //std::cout << std::endl;
523 : 0 : }
|
