Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Beta/MixMassFractionBetaCoeffPolicy.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 Mass-fraction beta SDE coefficients policies
9 : : \details This file defines coefficients policy classes for the mass-fraction
10 : : beta SDE, defined in DiffEq/Beta/MixMassFractionBeta.h. For general
11 : : requirements on mixture mass-fraction beta SDE coefficients policy
12 : : classes see the header file.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <iostream>
17 : :
18 : : #include "MixMassFractionBetaCoeffPolicy.hpp"
19 : : #include "TxtStatWriter.hpp"
20 : : #include "Walker/InputDeck/InputDeck.hpp"
21 : :
22 : : namespace walker {
23 : :
24 : : extern ctr::InputDeck g_inputdeck;
25 : :
26 : : } // ::walker
27 : :
28 [ - - ]: 0 : walker::MixMassFracBetaCoeffDecay::MixMassFracBetaCoeffDecay(
29 : : ncomp_t ncomp,
30 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime_,
31 : : const std::vector< kw::sde_S::info::expect::type >& S_,
32 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime_,
33 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2_,
34 : : const std::vector< kw::sde_r::info::expect::type >& r_,
35 : : std::vector< kw::sde_bprime::info::expect::type >& bprime,
36 : : std::vector< kw::sde_S::info::expect::type >& S,
37 : : std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
38 : : std::vector< kw::sde_rho2::info::expect::type >& rho2,
39 : : std::vector< kw::sde_r::info::expect::type >& r,
40 : : std::vector< kw::sde_b::info::expect::type >& b,
41 : : std::vector< kw::sde_kappa::info::expect::type >& k )
42 : : // *****************************************************************************
43 : : // Constructor: initialize coefficients
44 : : //! \param[in] ncomp Number of scalar components in this SDE system
45 : : //! \param[in] bprime_ Vector used to initialize coefficient vector bprime
46 : : //! \param[in] S_ Vector used to initialize coefficient vector S
47 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime
48 : : //! \param[in] rho2_ Vector used to initialize coefficient vector rho2
49 : : //! \param[in] r_ Vector used to initialize coefficient vector r
50 : : //! \param[in,out] bprime Coefficient vector to be initialized
51 : : //! \param[in,out] S Coefficient vector to be initialized
52 : : //! \param[in,out] kprime Coefficient vector to be initialized
53 : : //! \param[in,out] rho2 Coefficient vector to be initialized
54 : : //! \param[in,out] r Coefficient vector to be initialized
55 : : //! \param[in,out] b Coefficient vector to be initialized
56 : : //! \param[in,out] k Coefficient vector to be initialized
57 : : // *****************************************************************************
58 : : {
59 [ - - ][ - - ]: 0 : ErrChk( bprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
60 : : "Wrong number of mix mass-fraction beta SDE parameters 'b'");
61 [ - - ][ - - ]: 0 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
62 : : "Wrong number of mix mass-fraction beta SDE parameters 'S'");
63 [ - - ][ - - ]: 0 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
64 : : "Wrong number of mix mass-fraction beta SDE parameters 'kappa'");
65 [ - - ][ - - ]: 0 : ErrChk( rho2_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
66 : : "Wrong number of mix mass-fraction beta SDE parameters 'rho2'");
67 [ - - ][ - - ]: 0 : ErrChk( r_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
68 : : "Wrong number of mix mass-fraction beta SDE parameters 'r'");
69 : :
70 : 0 : bprime = bprime_;
71 : 0 : S = S_;
72 : 0 : kprime = kprime_;
73 : 0 : rho2 = rho2_;
74 : 0 : r = r_;
75 : :
76 : 0 : b.resize( bprime.size() );
77 : 0 : k.resize( kprime.size() );
78 : 0 : }
79 : :
80 : : void
81 : 0 : walker::MixMassFracBetaCoeffDecay::update(
82 : : char depvar,
83 : : char,
84 : : char,
85 : : ctr::DepvarType,
86 : : ctr::DepvarType /*scalar_solve*/,
87 : : ncomp_t ncomp,
88 : : const std::map< tk::ctr::Product, tk::real >& moments,
89 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime,
90 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
91 : : const std::vector< kw::sde_rho2::info::expect::type >&,
92 : : const std::vector< kw::sde_r::info::expect::type >&,
93 : : const std::vector< tk::Table<1> >&,
94 : : const std::vector< tk::Table<1> >&,
95 : : std::vector< kw::sde_b::info::expect::type >& b,
96 : : std::vector< kw::sde_kappa::info::expect::type >& k,
97 : : std::vector< kw::sde_S::info::expect::type >&,
98 : : tk::real ) const
99 : : // *****************************************************************************
100 : : // Update coefficients
101 : : //! \param[in] depvar Dependent variable
102 : : //! \param[in] ncomp Number of scalar components in this SDE system
103 : : //! \param[in] moments Map of statistical moments estimated
104 : : //! \param[in] bprime Coefficient vector b'
105 : : //! \param[in] kprime Coefficient vector kappa'
106 : : //! \param[in,out] b Coefficient vector to be updated
107 : : //! \param[in,out] k Coefficient vector to be updated
108 : : //! \details This where the mix mass-fraction beta SDE is made consistent
109 : : //! with the no-mix and fully mixed limits by specifying the SDE
110 : : //! coefficients, b and kappa as functions of b' and kappa'. We leave S
111 : : //! unchanged.
112 : : // *****************************************************************************
113 : : {
114 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
115 [ - - ]: 0 : tk::real m = tk::ctr::lookup( tk::ctr::mean(depvar,c), moments );
116 [ - - ]: 0 : tk::real v = tk::ctr::lookup( tk::ctr::variance(depvar,c), moments );
117 : :
118 [ - - ][ - - ]: 0 : if (m<1.0e-8 || m>1.0-1.0e-8) m = 0.5;
119 [ - - ][ - - ]: 0 : if (v<1.0e-8 || v>1.0-1.0e-8) v = 0.5;
120 : :
121 : 0 : b[c] = bprime[c] * (1.0 - v / m / ( 1.0 - m ));
122 : 0 : k[c] = kprime[c] * v;
123 : : }
124 : 0 : }
125 : :
126 [ - + ]: 11 : walker::MixMassFracBetaCoeffHomDecay::MixMassFracBetaCoeffHomDecay(
127 : : ncomp_t ncomp,
128 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime_,
129 : : const std::vector< kw::sde_S::info::expect::type >& S_,
130 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime_,
131 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2_,
132 : : const std::vector< kw::sde_r::info::expect::type >& r_,
133 : : std::vector< kw::sde_bprime::info::expect::type >& bprime,
134 : : std::vector< kw::sde_S::info::expect::type >& S,
135 : : std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
136 : : std::vector< kw::sde_rho2::info::expect::type >& rho2,
137 : : std::vector< kw::sde_r::info::expect::type >& r,
138 : : std::vector< kw::sde_b::info::expect::type >& b,
139 : : std::vector< kw::sde_kappa::info::expect::type >& k )
140 : : // *****************************************************************************
141 : : // Constructor: initialize coefficients
142 : : //! \param[in] ncomp Number of scalar components in this SDE system
143 : : //! \param[in] bprime_ Vector used to initialize coefficient vector bprime
144 : : //! \param[in] S_ Vector used to initialize coefficient vector S
145 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime
146 : : //! \param[in] rho2_ Vector used to initialize coefficient vector rho2
147 : : //! \param[in] r_ Vector used to initialize coefficient vector r
148 : : //! \param[in,out] bprime Coefficient vector to be initialized
149 : : //! \param[in,out] S Coefficient vector to be initialized
150 : : //! \param[in,out] kprime Coefficient vector to be initialized
151 : : //! \param[in,out] rho2 Coefficient vector to be initialized
152 : : //! \param[in,out] r Coefficient vector to be initialized
153 : : //! \param[in,out] b Coefficient vector to be initialized
154 : : //! \param[in,out] k Coefficient vector to be initialized
155 : : // *****************************************************************************
156 : : {
157 [ - + ][ - - ]: 11 : ErrChk( bprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
158 : : "Wrong number of mix mass-fraction beta SDE parameters 'b''");
159 [ - + ][ - - ]: 11 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
160 : : "Wrong number of mix mass-fraction beta SDE parameters 'S'");
161 [ - + ][ - - ]: 11 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
162 : : "Wrong number of mix mass-fraction beta SDE parameters 'kappa''");
163 [ - + ][ - - ]: 11 : ErrChk( rho2_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
164 : : "Wrong number of mix mass-fraction beta SDE parameters 'rho2'");
165 [ - + ][ - - ]: 11 : ErrChk( r_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
166 : : "Wrong number of mix mass-fraction beta SDE parameters 'r'");
167 : :
168 : 11 : bprime = bprime_;
169 : 11 : S = S_;
170 : 11 : kprime = kprime_;
171 : 11 : rho2 = rho2_;
172 : 11 : r = r_;
173 : :
174 : 11 : b.resize( bprime.size() );
175 : 11 : k.resize( kprime.size() );
176 : 11 : }
177 : :
178 : : void
179 : 23500 : walker::MixMassFracBetaCoeffHomDecay::update(
180 : : char depvar,
181 : : char,
182 : : char,
183 : : ctr::DepvarType,
184 : : ctr::DepvarType /*scalar_solve*/,
185 : : ncomp_t ncomp,
186 : : const std::map< tk::ctr::Product, tk::real >& moments,
187 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime,
188 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
189 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2,
190 : : const std::vector< kw::sde_r::info::expect::type >& r,
191 : : const std::vector< tk::Table<1> >&,
192 : : const std::vector< tk::Table<1> >&,
193 : : std::vector< kw::sde_b::info::expect::type >& b,
194 : : std::vector< kw::sde_kappa::info::expect::type >& k,
195 : : std::vector< kw::sde_S::info::expect::type >& S,
196 : : tk::real ) const
197 : : // *****************************************************************************
198 : : // Update coefficients
199 : : //! \param[in] depvar Dependent variable
200 : : //! \param[in] ncomp Number of scalar components in this SDE system
201 : : //! \param[in] moments Map of statistical moments estimated
202 : : //! \param[in] bprime Coefficient vector b'
203 : : //! \param[in] kprime Coefficient vector kappa'
204 : : //! \param[in] rho2 Coefficient vector rho2
205 : : //! \param[in] r Coefficient vector r
206 : : //! \param[in,out] b Coefficient vector to be updated
207 : : //! \param[in,out] k Coefficient vector to be updated
208 : : //! \param[in,out] S Coefficient vector to be updated
209 : : //! \details This where the mix mass-fraction beta SDE is made consistent
210 : : //! with the no-mix and fully mixed limits by specifying the SDE
211 : : //! coefficients, b and kappa as functions of b' and kappa'. We also
212 : : //! specify S to force d\<rho\>/dt = 0, where \<rho\> = rho_2/(1+rY).
213 : : // *****************************************************************************
214 : : {
215 : : using tk::ctr::lookup;
216 : : using tk::ctr::mean;
217 : : using tk::ctr::variance;
218 : : using tk::ctr::cen3;
219 : :
220 : : // statistics nomenclature:
221 : : // Y = instantaneous mass fraction,
222 : : // R = instantaneous density,
223 : : // y = Y - <Y>, mass fraction fluctuation about its mean,
224 : : // r = R - <R>, density fluctuation about its mean,
225 : : // <Y> = mean mass fraction,
226 : : // <R> = mean density,
227 : :
228 [ + + ]: 141000 : for (ncomp_t c=0; c<ncomp; ++c) {
229 [ + - ]: 117500 : tk::real m = lookup( mean(depvar,c), moments ); // <Y>
230 [ + - ]: 117500 : tk::real v = lookup( variance(depvar,c), moments ); // <y^2>
231 [ + - ]: 117500 : tk::real d = lookup( mean(depvar,c+ncomp), moments ); // <R>
232 [ + - ]: 117500 : tk::real d2 = lookup( variance(depvar,c+ncomp), moments ); // <r^2>
233 [ + - ]: 117500 : tk::real d3 = lookup( cen3(depvar,c+ncomp), moments ); // <r^3>
234 : :
235 [ + - ][ - + ]: 117500 : if (m<1.0e-8 || m>1.0-1.0e-8) m = 0.5;
236 [ + - ][ - + ]: 117500 : if (v<1.0e-8 || v>1.0-1.0e-8) v = 0.5;
237 [ - + ]: 117500 : b[c] = bprime[c] * (1.0 - v/m/(1.0-m));
238 : : //b[c] = bprime[c] * (1.0 - v/M[c]/(1.0-M[c]));
239 : 117500 : k[c] = kprime[c] * v;
240 : : //b[c] = 1.0;
241 : : //k[c] = 0.5*v/(m*(1.0-m));
242 : :
243 [ - + ]: 117500 : if (d < 1.0e-8) {
244 : 0 : std::cout << "d:" << d << " ";
245 : : d = 0.5;
246 : : }
247 [ + + ]: 117500 : tk::real R = 1.0 + d2/d/d;
248 : 117500 : tk::real B = -1.0/r[c]/r[c];
249 : 117500 : tk::real C = (2.0+r[c])/r[c]/r[c];
250 : 117500 : tk::real D = -(1.0+r[c])/r[c]/r[c];
251 : : tk::real diff =
252 : 117500 : B*d/rho2[c] +
253 : 117500 : C*d*d*R/rho2[c]/rho2[c] +
254 : 117500 : D*d*d*d*(1.0 + 3.0*d2/d/d + d3/d/d/d)/rho2[c]/rho2[c]/rho2[c];
255 [ + + ]: 117500 : S[c] = (rho2[c]/d/R +
256 : 117500 : 2.0*k[c]/b[c]*rho2[c]*rho2[c]/d/d*r[c]*r[c]/R*diff - 1.0) / r[c];
257 [ + + ][ - + ]: 117500 : if (S[c] < 0.0 || S[c] > 1.0) {
258 : 12 : std::cout << S[c] << " ";
259 : 12 : S[c] = 0.5;
260 : : }
261 : : }
262 : 23500 : }
263 : :
264 [ - - ]: 0 : walker::MixMassFracBetaCoeffMonteCarloHomDecay::
265 : : MixMassFracBetaCoeffMonteCarloHomDecay(
266 : : ncomp_t ncomp,
267 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime_,
268 : : const std::vector< kw::sde_S::info::expect::type >& S_,
269 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime_,
270 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2_,
271 : : const std::vector< kw::sde_r::info::expect::type >& r_,
272 : : std::vector< kw::sde_bprime::info::expect::type >& bprime,
273 : : std::vector< kw::sde_S::info::expect::type >& S,
274 : : std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
275 : : std::vector< kw::sde_rho2::info::expect::type >& rho2,
276 : : std::vector< kw::sde_r::info::expect::type >& r,
277 : : std::vector< kw::sde_b::info::expect::type >& b,
278 : : std::vector< kw::sde_kappa::info::expect::type >& k )
279 : : // *****************************************************************************
280 : : // Constructor: initialize coefficients
281 : : //! \param[in] ncomp Number of scalar components in this SDE system
282 : : //! \param[in] bprime_ Vector used to initialize coefficient vector bprime
283 : : //! \param[in] S_ Vector used to initialize coefficient vector S
284 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime
285 : : //! \param[in] rho2_ Vector used to initialize coefficient vector rho2
286 : : //! \param[in] r_ Vector used to initialize coefficient vector r
287 : : //! \param[in,out] bprime Coefficient vector to be initialized
288 : : //! \param[in,out] S Coefficient vector to be initialized
289 : : //! \param[in,out] kprime Coefficient vector to be initialized
290 : : //! \param[in,out] rho2 Coefficient vector to be initialized
291 : : //! \param[in,out] r Coefficient vector to be initialized
292 : : //! \param[in,out] b Coefficient vector to be initialized
293 : : //! \param[in,out] k Coefficient vector to be initialized
294 : : // *****************************************************************************
295 : : {
296 [ - - ][ - - ]: 0 : ErrChk( bprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
297 : : "Wrong number of mix mass-fraction beta SDE parameters 'b''");
298 [ - - ][ - - ]: 0 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
299 : : "Wrong number of mix mass-fraction beta SDE parameters 'S'");
300 [ - - ][ - - ]: 0 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
301 : : "Wrong number of mix mass-fraction beta SDE parameters 'kappa''");
302 [ - - ][ - - ]: 0 : ErrChk( rho2_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
303 : : "Wrong number of mix mass-fraction beta SDE parameters 'rho2'");
304 [ - - ][ - - ]: 0 : ErrChk( r_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
305 : : "Wrong number of mix mass-fraction beta SDE parameters 'r'");
306 : :
307 : 0 : bprime = bprime_;
308 : 0 : S = S_;
309 : 0 : kprime = kprime_;
310 : 0 : rho2 = rho2_;
311 : 0 : r = r_;
312 : :
313 : 0 : b.resize( bprime.size() );
314 : 0 : k.resize( kprime.size() );
315 : 0 : }
316 : :
317 : : void
318 : 0 : walker::MixMassFracBetaCoeffMonteCarloHomDecay::update(
319 : : char depvar,
320 : : char,
321 : : char,
322 : : ctr::DepvarType,
323 : : ctr::DepvarType /*scalar_solve*/,
324 : : ncomp_t ncomp,
325 : : const std::map< tk::ctr::Product, tk::real >& moments,
326 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime,
327 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
328 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2,
329 : : const std::vector< kw::sde_r::info::expect::type >& r,
330 : : const std::vector< tk::Table<1> >&,
331 : : const std::vector< tk::Table<1> >&,
332 : : std::vector< kw::sde_b::info::expect::type >& b,
333 : : std::vector< kw::sde_kappa::info::expect::type >& k,
334 : : std::vector< kw::sde_S::info::expect::type >& S,
335 : : tk::real ) 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] moments Map of statistical moments estimated
341 : : //! \param[in] bprime Coefficient vector b'
342 : : //! \param[in] kprime Coefficient vector kappa'
343 : : //! \param[in] rho2 Coefficient vector rho2
344 : : //! \param[in] r Coefficient vector r
345 : : //! \param[in,out] b Coefficient vector to be updated
346 : : //! \param[in,out] k Coefficient vector to be updated
347 : : //! \param[in,out] S Coefficient vector to be updated
348 : : //! \details This where the mix mass-fraction beta SDE is made consistent
349 : : //! with the no-mix and fully mixed limits by specifying the SDE
350 : : //! coefficients, b and kappa as functions of b' and kappa'. We also
351 : : //! specify S to force d\<rho\>/dt = 0, where \<rho\> = rho_2/(1+rY).
352 : : // *****************************************************************************
353 : : {
354 : : using tk::ctr::lookup;
355 : : using tk::ctr::mean;
356 : : using tk::ctr::variance;
357 : : using tk::ctr::ord2;
358 : : // statistics nomenclature:
359 : : // Y = instantaneous mass fraction,
360 : : // R = instantaneous density,
361 : : // y = Y - <Y>, mass fraction fluctuation about its mean,
362 : : // r = R - <R>, density fluctuation about its mean,
363 : : // <Y> = mean mass fraction,
364 : : // <R> = mean density,
365 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
366 [ - - ]: 0 : tk::real m = lookup( mean(depvar,c), moments ); // <Y>
367 [ - - ]: 0 : tk::real v = lookup( variance(depvar,c), moments ); // <y^2>
368 [ - - ]: 0 : tk::real r2 = lookup( ord2(depvar,c+ncomp), moments ); // <R^2>
369 : :
370 : 0 : const tk::ctr::Term Y( static_cast<char>(std::toupper(depvar)),
371 : : c,
372 : 0 : tk::ctr::Moment::ORDINARY );
373 : 0 : const tk::ctr::Term R( static_cast<char>(std::toupper(depvar)),
374 : : c+ncomp,
375 : 0 : tk::ctr::Moment::ORDINARY );
376 : 0 : const tk::ctr::Term OneMinusY( static_cast<char>(std::toupper(depvar)),
377 : 0 : c+3*ncomp,
378 : 0 : tk::ctr::Moment::ORDINARY );
379 : :
380 [ - - ]: 0 : const auto YR2 = tk::ctr::Product( { Y, R, R } );
381 [ - - ][ - - ]: 0 : const auto Y1MYR3 = tk::ctr::Product( { Y, OneMinusY, R, R, R } );
[ - - ]
382 : :
383 [ - - ]: 0 : tk::real yr2 = lookup( YR2, moments ); // <RY^2>
384 [ - - ]: 0 : tk::real y1myr3 = lookup( Y1MYR3, moments ); // <Y(1-Y)R^3>
385 : :
386 [ - - ][ - - ]: 0 : if (m<1.0e-8 || m>1.0-1.0e-8) m = 0.5;
387 [ - - ][ - - ]: 0 : if (v<1.0e-8 || v>1.0-1.0e-8) v = 0.5;
388 [ - - ]: 0 : b[c] = bprime[c] * (1.0 - v/m/(1.0-m));
389 : 0 : k[c] = kprime[c] * v;
390 : : //b[c] = 1.0;
391 : : //k[c] = 0.5*v/(m*(1.0-m));
392 : :
393 [ - - ]: 0 : if (r2 < 1.0e-8) {
394 [ - - ]: 0 : std::cout << "r2:" << r2 << " ";
395 : : r2 = 0.5;
396 : : }
397 [ - - ]: 0 : S[c] = (yr2 + 2.0*k[c]/b[c]*r[c]/rho2[c]*y1myr3) / r2;
398 [ - - ][ - - ]: 0 : if (S[c] < 0.0 || S[c] > 1.0) {
399 [ - - ][ - - ]: 0 : std::cout << "S:" << S[c] << " ";
400 : 0 : S[c] = 0.5;
401 : : }
402 : : }
403 : 0 : }
404 : :
405 : 4 : walker::MixMassFracBetaCoeffHydroTimeScale::
406 : : MixMassFracBetaCoeffHydroTimeScale(
407 : : ncomp_t ncomp,
408 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime_,
409 : : const std::vector< kw::sde_S::info::expect::type >& S_,
410 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime_,
411 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2_,
412 : : const std::vector< kw::sde_r::info::expect::type >& r_,
413 : : std::vector< kw::sde_bprime::info::expect::type >& bprime,
414 : : std::vector< kw::sde_S::info::expect::type >& S,
415 : : std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
416 : : std::vector< kw::sde_rho2::info::expect::type >& rho2,
417 : : std::vector< kw::sde_r::info::expect::type >& r,
418 : : std::vector< kw::sde_b::info::expect::type >& b,
419 [ - + ]: 4 : std::vector< kw::sde_kappa::info::expect::type >& k )
420 : : // *****************************************************************************
421 : : // Constructor: initialize coefficients
422 : : //! \param[in] ncomp Number of scalar components in this SDE system
423 : : //! \param[in] bprime_ Vector used to initialize coefficient vector bprime
424 : : //! \param[in] S_ Vector used to initialize coefficient vector S
425 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime
426 : : //! \param[in] rho2_ Vector used to initialize coefficient vector rho2
427 : : //! \param[in] r_ Vector used to initialize coefficient vector r
428 : : //! \param[in,out] bprime Coefficient vector to be initialized
429 : : //! \param[in,out] S Coefficient vector to be initialized
430 : : //! \param[in,out] kprime Coefficient vector to be initialized
431 : : //! \param[in,out] rho2 Coefficient vector to be initialized
432 : : //! \param[in,out] r Coefficient vector to be initialized
433 : : //! \param[in,out] b Coefficient vector to be initialized
434 : : //! \param[in,out] k Coefficient vector to be initialized
435 : : // *****************************************************************************
436 : : {
437 [ - + ][ - - ]: 4 : ErrChk( bprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
438 : : "Wrong number of mix mass-fraction beta SDE parameters 'b''");
439 [ - + ][ - - ]: 4 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
440 : : "Wrong number of mix mass-fraction beta SDE parameters 'S'");
441 [ - + ][ - - ]: 4 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
442 : : "Wrong number of mix mass-fraction beta SDE parameters 'kappa''");
443 [ - + ][ - - ]: 4 : ErrChk( rho2_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
444 : : "Wrong number of mix mass-fraction beta SDE parameters 'rho2'");
445 [ - + ][ - - ]: 4 : ErrChk( r_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
446 : : "Wrong number of mix mass-fraction beta SDE parameters 'r'");
447 : :
448 [ + - ]: 4 : bprime = bprime_;
449 [ + - ]: 4 : S = S_;
450 [ + - ]: 4 : kprime = kprime_;
451 [ + - ]: 4 : rho2 = rho2_;
452 [ + - ]: 4 : r = r_;
453 : :
454 [ + - ]: 4 : b.resize( bprime.size() );
455 [ + - ]: 4 : k.resize( kprime.size() );
456 : :
457 : : // Extra output besides normal statistics output
458 : : m_extra_out_filename = "coeff";
459 [ + - ]: 4 : tk::TxtStatWriter sw( m_extra_out_filename );
460 : 4 : std::vector< std::string > names;
461 [ + + ]: 24 : for (ncomp_t c=0; c<ncomp; ++c) {
462 [ + - ][ + - ]: 40 : names.push_back( "b" + std::to_string(c+1) );
[ - + ][ - - ]
463 [ + - ][ + - ]: 40 : names.push_back( "S" + std::to_string(c+1) );
[ - + ][ - - ]
464 [ + - ][ + - ]: 40 : names.push_back( "k" + std::to_string(c+1) );
[ - + ][ - - ]
465 : : }
466 [ + - ]: 4 : sw.header( names, {}, {} );
467 : :
468 [ + - ]: 4 : m_s.resize(ncomp);
469 : 4 : }
470 : :
471 : : void
472 : 36 : walker::MixMassFracBetaCoeffHydroTimeScale::update(
473 : : char depvar,
474 : : char,
475 : : char,
476 : : ctr::DepvarType,
477 : : ctr::DepvarType /*scalar_solve*/,
478 : : ncomp_t ncomp,
479 : : const std::map< tk::ctr::Product, tk::real >& moments,
480 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime,
481 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
482 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2,
483 : : const std::vector< kw::sde_r::info::expect::type >& r,
484 : : const std::vector< tk::Table<1> >& hts,
485 : : const std::vector< tk::Table<1> >& hp,
486 : : std::vector< kw::sde_b::info::expect::type >& b,
487 : : std::vector< kw::sde_kappa::info::expect::type >& k,
488 : : std::vector< kw::sde_S::info::expect::type >& S,
489 : : tk::real t ) const
490 : : // *****************************************************************************
491 : : // Update coefficients
492 : : //! \param[in] depvar Dependent variable
493 : : //! \param[in] ncomp Number of scalar components in this SDE system
494 : : //! \param[in] moments Map of statistical moments estimated
495 : : //! \param[in] bprime Coefficient vector b'
496 : : //! \param[in] kprime Coefficient vector kappa'
497 : : //! \param[in] rho2 Coefficient vector rho2
498 : : //! \param[in] r Coefficient vector r
499 : : //! \param[in] hts (Inverse) hydrodynamics time scale (as input from DNS)
500 : : //! \param[in] hp Production/dissipation (as input from DNS)
501 : : //! \param[in,out] b Coefficient vector to be updated
502 : : //! \param[in,out] k Coefficient vector to be updated
503 : : //! \param[in,out] S Coefficient vector to be updated
504 : : //! \param[in] t Physical time at which to update coefficients
505 : : //! \details This where the mix mass-fraction beta SDE is made consistent
506 : : //! with the no-mix and fully mixed limits by specifying the SDE
507 : : //! coefficients, b and kappa as functions of b' and kappa'. Additionally,
508 : : //! we pull in a hydrodynamic timescale from an external function. We also
509 : : //! specify S to force d\<rho\>/dt = 0, where \<rho\> = rho_2/(1+rY).
510 : : // *****************************************************************************
511 : : {
512 : : using tk::ctr::lookup;
513 : : using tk::ctr::mean;
514 : : using tk::ctr::variance;
515 : : using tk::ctr::cen3;
516 : : using tk::ctr::Product;
517 : :
518 [ + + ][ + + ]: 56 : if (m_it == 0) for (ncomp_t c=0; c<ncomp; ++c) m_s[c] = S[c];
519 : :
520 : : // Extra output besides normal statistics output
521 : 36 : tk::TxtStatWriter sw( m_extra_out_filename,
522 : : g_inputdeck.get< tag::flformat, tag::stat >(),
523 : : g_inputdeck.get< tag::prec, tag::stat >(),
524 : 36 : std::ios_base::app );
525 : :
526 [ + - ]: 36 : std::vector< tk::real > coeffs( ncomp * 3 );
527 : :
528 : : // statistics nomenclature:
529 : : // Y = instantaneous mass fraction,
530 : : // R = instantaneous density,
531 : : // y = Y - <Y>, mass fraction fluctuation about its mean,
532 : : // r = R - <R>, density fluctuation about its mean,
533 : : // <Y> = mean mass fraction,
534 : : // <R> = mean density,
535 [ + + ]: 216 : for (ncomp_t c=0; c<ncomp; ++c) {
536 : :
537 : 180 : const tk::ctr::Term Y( static_cast<char>(std::toupper(depvar)),
538 : : c,
539 : 180 : tk::ctr::Moment::ORDINARY );
540 : 180 : const tk::ctr::Term dens( static_cast<char>(std::toupper(depvar)),
541 : : c+ncomp,
542 : 180 : tk::ctr::Moment::ORDINARY );
543 : 180 : const tk::ctr::Term s1( static_cast<char>(std::tolower(depvar)),
544 : : c+ncomp,
545 : 180 : tk::ctr::Moment::CENTRAL );
546 : 180 : const tk::ctr::Term s2( static_cast<char>(std::tolower(depvar)),
547 : 180 : c+ncomp*2,
548 : 180 : tk::ctr::Moment::CENTRAL );
549 : :
550 [ + - ][ + - ]: 180 : const auto RY = tk::ctr::Product( { dens, Y } );
[ - - ]
551 [ + - ]: 180 : tk::real ry = lookup( RY, moments ); // <RY>
552 [ + - ][ + - ]: 180 : const auto dscorr = tk::ctr::Product( { s1, s2 } );
[ - - ]
553 [ + - ]: 180 : tk::real ds = -lookup( dscorr, moments ); // b = -<rv>
554 [ + - ][ + - ]: 180 : tk::real d = lookup( mean(depvar,c+ncomp), moments ); // <R>
555 [ + - ][ + - ]: 180 : tk::real d2 = lookup( variance(depvar,c+ncomp), moments ); // <r^2>
556 [ + - ][ + - ]: 180 : tk::real d3 = lookup( cen3(depvar,c+ncomp), moments ); // <r^3>
[ - - ]
557 : 180 : tk::real yt = ry/d;
558 : :
559 : : // Sample hydrodynamics timescale and prod/diss at time t
560 : : auto ts = hydrotimescale( t, hts[c] ); // eps/k
561 : : auto pe = hydroproduction( t, hp[c] ); // P/eps = (dk/dt+eps)/eps
562 : :
563 : 180 : tk::real a = r[c]/(1.0+r[c]*yt);
564 : 180 : tk::real bnm = a*a*yt*(1.0-yt);
565 : 180 : tk::real thetab = 1.0 - ds/bnm;
566 : : tk::real f2 =
567 [ + - ]: 180 : 1.0 / std::pow(1.0 + std::pow(pe-1.0,2.0)*std::pow(ds,0.25),0.5);
568 : 180 : tk::real b1 = m_s[0];
569 : 180 : tk::real b2 = m_s[1];
570 : 180 : tk::real b3 = m_s[2];
571 : 180 : tk::real eta = d2/d/d/ds;
572 : 180 : tk::real beta2 = b2*(1.0+eta*ds);
573 : 180 : tk::real Thetap = thetab*0.5*(1.0+eta/(1.0+eta*ds));
574 : 180 : tk::real beta3 = b3*(1.0+eta*ds);
575 : 180 : tk::real beta10 = b1 * (1.0+ds)/(1.0+eta*ds);
576 : 180 : tk::real beta1 = bprime[c] * 2.0/(1.0+eta+eta*ds) *
577 : 180 : (beta10 + beta2*Thetap*f2 + beta3*Thetap*(1.0-Thetap)*f2);
578 : 180 : b[c] = beta1 * ts;
579 : 180 : k[c] = kprime[c] * beta1 * ts * ds * ds;
580 : : //b[c] = bprime[c];
581 : : //k[c] = kprime[c];
582 : : //b[c] = bprime[c] + 0.25*std::sin(10.0*t);
583 : : //k[c] = 1.0 + 0.25*std::sin(10.0*t);
584 : : //k[c] = -(1.0 + std::sin(t)) * (S[c] - 1.0);
585 : :
586 : 180 : tk::real R = 1.0 + d2/d/d;
587 : 180 : tk::real B = -1.0/r[c]/r[c];
588 : 180 : tk::real C = (2.0+r[c])/r[c]/r[c];
589 : 180 : tk::real D = -(1.0+r[c])/r[c]/r[c];
590 : : tk::real diff =
591 : 180 : B*d/rho2[c] +
592 : 180 : C*d*d*R/rho2[c]/rho2[c] +
593 : 180 : D*d*d*d*(1.0 + 3.0*d2/d/d + d3/d/d/d)/rho2[c]/rho2[c]/rho2[c];
594 : 180 : S[c] = (rho2[c]/d/R +
595 : 180 : 2.0*k[c]/b[c]*rho2[c]*rho2[c]/d/d*r[c]*r[c]/R*diff - 1.0) / r[c];
596 : : //S[c] = 0.5 + 0.25*std::sin(10.0*t);
597 : :
598 : : // Implementation of a constraint for MixDirichlet for S_al to keep
599 : : // d<rho>/dt = 0 to verify its behavior for beta (binary case). As input
600 : : // file use mixmassfractbeta_mmS_A0.75.q.
601 : : // auto R2 = lookup( Product({dens,dens}), moments ); // <R^2>
602 : : // auto R2Yc = lookup( Product({dens,dens,Y}), moments ); // <R^2Yc>
603 : : // auto R3Yc = lookup( Product({dens,dens,dens,Y}), moments ); // <R^3Yc>
604 : : // auto R3Y2c = lookup( Product({dens,dens,dens,Y,Y}), moments ); // <R^3Y^2>
605 : : // tk::real drYc = -r[c]/rho2[c]*R2;
606 : : // tk::real drYcYc = -r[c]/rho2[c]*R2Yc;
607 : : // tk::real drYc2YcYN = 2.0*std::pow(r[c]/rho2[c],2.0)*(R3Yc-R3Y2c);
608 : : // S[c] = (drYcYc - k[c]/b[c]*drYc2YcYN) / drYc;
609 : :
610 : 180 : coeffs[ 3*c+0 ] = b[c];
611 : 180 : coeffs[ 3*c+1 ] = S[c];
612 [ + - ]: 180 : coeffs[ 3*c+2 ] = k[c];
613 : : }
614 : :
615 : : // Extra "stat" output of coefficients
616 [ + - ][ - + ]: 72 : sw.stat( 0, t, coeffs, {}, {} );
[ - - ][ - - ]
617 : :
618 [ + - ]: 36 : ++m_it;
619 : 36 : }
620 : :
621 : 0 : walker::MixMassFracBetaCoeffInstVel::MixMassFracBetaCoeffInstVel(
622 : : ncomp_t ncomp,
623 : : const std::vector< kw::sde_bprime::info::expect::type >& bprime_,
624 : : const std::vector< kw::sde_S::info::expect::type >& S_,
625 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime_,
626 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2_,
627 : : const std::vector< kw::sde_r::info::expect::type >& r_,
628 : : std::vector< kw::sde_bprime::info::expect::type >& bprime,
629 : : std::vector< kw::sde_S::info::expect::type >& S,
630 : : std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
631 : : std::vector< kw::sde_rho2::info::expect::type >& rho2,
632 : : std::vector< kw::sde_r::info::expect::type >& r,
633 : : std::vector< kw::sde_b::info::expect::type >& b,
634 [ - - ]: 0 : std::vector< kw::sde_kappa::info::expect::type >& k )
635 : : // *****************************************************************************
636 : : // Constructor: initialize coefficients
637 : : //! \param[in] ncomp Number of scalar components in this SDE system
638 : : //! \param[in] bprime_ Vector used to initialize coefficient vector bprime
639 : : //! \param[in] S_ Vector used to initialize coefficient vector S
640 : : //! \param[in] kprime_ Vector used to initialize coefficient vector kprime
641 : : //! \param[in] rho2_ Vector used to initialize coefficient vector rho2
642 : : //! \param[in] r_ Vector used to initialize coefficient vector r
643 : : //! \param[in,out] bprime Coefficient vector to be initialized
644 : : //! \param[in,out] S Coefficient vector to be initialized
645 : : //! \param[in,out] kprime Coefficient vector to be initialized
646 : : //! \param[in,out] rho2 Coefficient vector to be initialized
647 : : //! \param[in,out] r Coefficient vector to be initialized
648 : : //! \param[in,out] b Coefficient vector to be initialized
649 : : //! \param[in,out] k Coefficient vector to be initialized
650 : : // *****************************************************************************
651 : : {
652 [ - - ][ - - ]: 0 : ErrChk( bprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
653 : : "Wrong number of mix mass-fraction beta SDE parameters 'b''");
654 [ - - ][ - - ]: 0 : ErrChk( S_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
655 : : "Wrong number of mix mass-fraction beta SDE parameters 'S'");
656 [ - - ][ - - ]: 0 : ErrChk( kprime_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
657 : : "Wrong number of mix mass-fraction beta SDE parameters 'kappa''");
658 [ - - ][ - - ]: 0 : ErrChk( rho2_.size() == ncomp,
[ - - ][ - - ]
[ - - ]
659 : : "Wrong number of mix mass-fraction beta SDE parameters 'rho2'");
660 [ - - ][ - - ]: 0 : ErrChk( r_.size() == ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
661 : : "Wrong number of mix mass-fraction beta SDE parameters 'r'");
662 : :
663 [ - - ]: 0 : bprime = bprime_;
664 [ - - ]: 0 : S = S_;
665 [ - - ]: 0 : kprime = kprime_;
666 [ - - ]: 0 : rho2 = rho2_;
667 [ - - ]: 0 : r = r_;
668 : :
669 [ - - ]: 0 : b.resize( bprime.size() );
670 [ - - ]: 0 : k.resize( kprime.size() );
671 [ - - ]: 0 : m_s.resize(ncomp);
672 : 0 : }
673 : :
674 : : void
675 : 0 : walker::MixMassFracBetaCoeffInstVel::update(
676 : : char depvar,
677 : : char dissipation_depvar,
678 : : char /*velocity_depvar*/,
679 : : ctr::DepvarType /*velocity_solve*/,
680 : : ctr::DepvarType solve,
681 : : ncomp_t ncomp,
682 : : const std::map< tk::ctr::Product, tk::real >& moments,
683 : : const std::vector< kw::sde_bprime::info::expect::type >& /*bprime*/,
684 : : const std::vector< kw::sde_kappaprime::info::expect::type >& kprime,
685 : : const std::vector< kw::sde_rho2::info::expect::type >& rho2,
686 : : const std::vector< kw::sde_r::info::expect::type >& r,
687 : : const std::vector< tk::Table<1> >&,
688 : : const std::vector< tk::Table<1> >&,
689 : : std::vector< kw::sde_b::info::expect::type >& b,
690 : : std::vector< kw::sde_kappa::info::expect::type >& k,
691 : : std::vector< kw::sde_S::info::expect::type >& S,
692 : : tk::real ) const
693 : : // *****************************************************************************
694 : : // Update coefficients
695 : : //! \param[in] depvar Dependent variable
696 : : //! \param[in] dissipation_depvar Dependent variable of coupled dissipation eq
697 : : //! \param[in] solve Enum selecting whether the full variable or its
698 : : //! fluctuation is solved for
699 : : //! \param[in] ncomp Number of scalar components in this SDE system
700 : : //! \param[in] moments Map of statistical moments estimated
701 : : //! \param[in] kprime Coefficient vector kappa'
702 : : //! \param[in] rho2 Coefficient vector rho2
703 : : //! \param[in] r Coefficient vector r
704 : : //! \param[in,out] b Coefficient vector to be updated
705 : : //! \param[in,out] k Coefficient vector to be updated
706 : : //! \param[in,out] S Coefficient vector to be updated
707 : : //! \details This where the mix mass-fraction beta SDE is made consistent
708 : : //! with the no-mix and fully mixed limits by specifying the SDE
709 : : //! coefficients, b and kappa as functions of b' and kappa'. Additionally,
710 : : //! we specify the hydrodynamic timescale by coupling to anothe SDE, computing
711 : : //! the velocity field. We also specify S to force d\<rho\>/dt = 0, where
712 : : //! \<rho\> = rho_2/(1+rY).
713 : : // *****************************************************************************
714 : : {
715 : : using tk::ctr::lookup;
716 : : using tk::ctr::mean;
717 : : using tk::ctr::variance;
718 : : using tk::ctr::cen3;
719 : :
720 [ - - ][ - - ]: 0 : if (m_it == 0) for (ncomp_t c=0; c<ncomp; ++c) m_s[c] = S[c];
721 : :
722 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c) {
723 : :
724 : : tk::real y2 = 0.0;
725 : : tk::real ts = 1.0;
726 : :
727 [ - - ]: 0 : if (solve == ctr::DepvarType::FULLVAR) {
728 : :
729 [ - - ]: 0 : y2 = lookup( variance(depvar,c), moments ); // <y^2>
730 : :
731 : : // Access mean turbulence frequency from coupled dissipation model
732 : : // hydroptimescale: eps/k = <O>
733 [ - - ]: 0 : if (dissipation_depvar != '-') { // only if dissipation is coupled
734 [ - - ]: 0 : ts = lookup( mean(dissipation_depvar,0), moments );
735 : : }
736 : :
737 [ - - ]: 0 : } else if (solve == ctr::DepvarType::FLUCTUATION) {
738 : :
739 : : // Since we are solving for the fluctuating scalar, the "ordinary"
740 : : // moments are really central moments
741 : 0 : auto d = static_cast< char >( std::toupper( depvar ) );
742 : : tk::ctr::Term y( d, c, tk::ctr::Moment::ORDINARY );
743 [ - - ][ - - ]: 0 : y2 = lookup( tk::ctr::Product( {y,y} ), moments ); // <y^2>
744 : :
745 : : // Access mean turbulence frequency from coupled dissipation model
746 : : // hydroptimescale: eps/k = <O>
747 [ - - ]: 0 : if (dissipation_depvar != '-') { // only if dissipation is coupled
748 [ - - ]: 0 : ts = lookup( mean(dissipation_depvar,0), moments );
749 : : }
750 : :
751 [ - - ][ - - ]: 0 : } else Throw( "Depvar type not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
752 : :
753 : : // simple decay for now
754 : : tk::real beta1 = 2.0;
755 : 0 : b[c] = beta1 * ts;
756 : 0 : k[c] = kprime[c] * beta1 * ts * y2;
757 : :
758 [ - - ]: 0 : tk::real d = lookup( mean(depvar,c+ncomp), moments ); // <R>
759 [ - - ]: 0 : tk::real d2 = lookup( variance(depvar,c+ncomp), moments ); // <r^2>
760 [ - - ]: 0 : tk::real d3 = lookup( cen3(depvar,c+ncomp), moments ); // <r^3>
761 : :
762 : : // force d\<rho\>/dt = 0
763 : 0 : tk::real R = 1.0 + d2/d/d;
764 : 0 : tk::real B = -1.0/r[c]/r[c];
765 : 0 : tk::real C = (2.0+r[c])/r[c]/r[c];
766 : 0 : tk::real D = -(1.0+r[c])/r[c]/r[c];
767 : : tk::real diff =
768 : 0 : B*d/rho2[c] +
769 : 0 : C*d*d*R/rho2[c]/rho2[c] +
770 : 0 : D*d*d*d*(1.0 + 3.0*d2/d/d + d3/d/d/d)/rho2[c]/rho2[c]/rho2[c];
771 : 0 : S[c] = (rho2[c]/d/R +
772 : 0 : 2.0*k[c]/b[c]*rho2[c]*rho2[c]/d/d*r[c]*r[c]/R*diff - 1.0) / r[c];
773 : : }
774 : :
775 : 0 : ++m_it;
776 : 0 : }
|
