TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Statistiques_dns_ijk_FT.cpp
1/****************************************************************************
2* Copyright (c) 2015 - 2016, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#include <Statistiques_dns_ijk_FT.h>
17#include <IJK_Field_vector.h>
18#include <Domaine_IJK.h>
19#include <TRUSTTab.h>
20#include <communications.h>
21#include <Param.h>
22#include <IJK_Navier_Stokes_tools.h>
23#include <Probleme_FTD_IJK_base.h>
24
25#define CREATE_STR(sum_id) #sum_id // turn the enum value into a string.
26
27// #define STAT_VERBOSE
28// #define STAT_EXIT
29#define PRECISION_DERIVEES (100.1) // (1e-10)
30#define PRECISION_DIVU (1e-10)
31#define PRECISION_DDIVU (1e-10)
32#define PRECISION_LAPLP (1e-10)
33
34Implemente_instanciable_sans_constructeur(Statistiques_dns_ijk_FT,"Statistiques_dns_ijk_FT",Statistiques_dns_ijk);
35#include <Statistiques_defines.h>
36#include <Statistiques_defines_thermique.h>
37// Pour generer la liste XXX_MOY, envoyer la liste des #define dans
38// sed -e '/\/\//d' -e 's/#define /"/g;s/_MOY.*/",/g' /tmp/titi >/tmp/resultat
39// et copier ci-dessous:
40// Ou inversement :
41// sed -e "s/\"//g" -e "s/,/_MOY/g" list | awk '{a=NR+476;print "#define " $0 " " a}' > list_define
42
44 check_stats_(0),
46 nvalt_(0)
47{
48// I : Indicatrice.
49// Iv : Indicatrice vapeur (1. - chi)
50 const char *noms_moyennes_prov[] =
51 {
52 "I",
53 "UI",
54 "VI",
55 "WI",
56 "PI",
57 "UIv",
58 "VIv",
59 "WIv",
60 "PIv",
61 "UUI",
62 "VVI",
63 "WWI",
64 "UVI",
65 "VWI",
66 "UWI",
67 "UUUI",
68 "UUVI",
69 "UUWI",
70 "UVVI",
71 "UVWI",
72 "UWWI",
73 "VVVI",
74 "VVWI",
75 "VWWI",
76 "WWWI",
77 "UPI",
78 "VPI",
79 "WPI",
80 "UdIdx",
81 "VdIdx",
82 "WdIdx",
83 "UdIdy",
84 "VdIdy",
85 "WdIdy",
86 "UdIdz",
87 "VdIdz",
88 "WdIdz",
89 "UPdIdx",
90 "VPdIdx",
91 "WPdIdx",
92 "UPdIdy",
93 "VPdIdy",
94 "WPdIdy",
95 "UPdIdz",
96 "VPdIdz",
97 "WPdIdz",
98 "UUdIdx",
99 "VVdIdx",
100 "WWdIdx",
101 "UUdIdy",
102 "VVdIdy",
103 "WWdIdy",
104 "UUdIdz",
105 "VVdIdz",
106 "WWdIdz",
107 "UVdIdx",
108 "UWdIdx",
109 "VWdIdx",
110 "UVdIdy",
111 "UWdIdy",
112 "VWdIdy",
113 "UVdIdz",
114 "UWdIdz",
115 "VWdIdz",
116 "PdIdx",
117 "PdIdy",
118 "PdIdz",
119 "IdIdx",
120 "IdIdy",
121 "IdIdz",
122 "IdUdx",
123 "IdVdx",
124 "IdWdx",
125 "IdUdy",
126 "IdVdy",
127 "IdWdy",
128 "IdUdz",
129 "IdVdz",
130 "IdWdz",
131 "IPdUdx",
132 "IPdVdx",
133 "IPdWdx",
134 "IPdUdy",
135 "IPdVdy",
136 "IPdWdy",
137 "IPdUdz",
138 "IPdVdz",
139 "IPdWdz",
140 "dUdxdIdx",
141 "dVdxdIdx",
142 "dWdxdIdx",
143 "dUdydIdy",
144 "dVdydIdy",
145 "dWdydIdy",
146 "dUdzdIdz",
147 "dVdzdIdz",
148 "dWdzdIdz",
149 "UdUdxdIdx",
150 "UdVdxdIdx",
151 "UdWdxdIdx",
152 "UdUdydIdy",
153 "UdVdydIdy",
154 "UdWdydIdy",
155 "UdUdzdIdz",
156 "UdVdzdIdz",
157 "UdWdzdIdz",
158 "VdUdxdIdx",
159 "VdVdxdIdx",
160 "VdWdxdIdx",
161 "VdUdydIdy",
162 "VdVdydIdy",
163 "VdWdydIdy",
164 "VdUdzdIdz",
165 "VdVdzdIdz",
166 "VdWdzdIdz",
167 "WdUdxdIdx",
168 "WdVdxdIdx",
169 "WdWdxdIdx",
170 "WdUdydIdy",
171 "WdVdydIdy",
172 "WdWdydIdy",
173 "WdUdzdIdz",
174 "WdVdzdIdz",
175 "WdWdzdIdz",
176 "IdUdxdUdx",
177 "IdUdxdUdy",
178 "IdUdxdUdz",
179 "IdUdxdVdx",
180 "IdUdxdVdy",
181 "IdUdxdVdz",
182 "IdUdxdWdx",
183 "IdUdxdWdy",
184 "IdUdxdWdz",
185 "IdUdydUdy",
186 "IdUdydUdz",
187 "IdUdydVdx",
188 "IdUdydVdy",
189 "IdUdydVdz",
190 "IdUdydWdx",
191 "IdUdydWdy",
192 "IdUdydWdz",
193 "IdUdzdUdz",
194 "IdUdzdVdx",
195 "IdUdzdVdy",
196 "IdUdzdVdz",
197 "IdUdzdWdx",
198 "IdUdzdWdy",
199 "IdUdzdWdz",
200 "IdVdxdVdx",
201 "IdVdxdVdy",
202 "IdVdxdVdz",
203 "IdVdxdWdx",
204 "IdVdxdWdy",
205 "IdVdxdWdz",
206 "IdVdydVdy",
207 "IdVdydVdz",
208 "IdVdydWdx",
209 "IdVdydWdy",
210 "IdVdydWdz",
211 "IdVdzdVdz",
212 "IdVdzdWdx",
213 "IdVdzdWdy",
214 "IdVdzdWdz",
215 "IdWdxdWdx",
216 "IdWdxdWdy",
217 "IdWdxdWdz",
218 "IdWdydWdy",
219 "IdWdydWdz",
220 "IdWdzdWdz",
221 "IdPdx",
222 "IdPdy",
223 "IdPdz",
224 "IUdPdx",
225 "IUdPdy",
226 "IUdPdz",
227 "IVdPdx",
228 "IVdPdy",
229 "IVdPdz",
230 "IWdPdx",
231 "IWdPdy",
232 "IWdPdz",
233 "IUdUdx",
234 "IUdVdx",
235 "IUdWdx",
236 "IUdUdy",
237 "IUdVdy",
238 "IUdWdy",
239 "IUdUdz",
240 "IUdVdz",
241 "IUdWdz",
242 "IVdUdx",
243 "IVdVdx",
244 "IVdWdx",
245 "IVdUdy",
246 "IVdVdy",
247 "IVdWdy",
248 "IVdUdz",
249 "IVdVdz",
250 "IVdWdz",
251 "IWdUdx",
252 "IWdVdx",
253 "IWdWdx",
254 "IWdUdy",
255 "IWdVdy",
256 "IWdWdy",
257 "IWdUdz",
258 "IWdVdz",
259 "IWdWdz",
260 "IddUdxx",
261 "IddVdxx",
262 "IddWdxx",
263 "IddUdyy",
264 "IddVdyy",
265 "IddWdyy",
266 "IddUdzz",
267 "IddVdzz",
268 "IddWdzz",
269 "Fx",
270 "Fy",
271 "Fz",
272 "Frx",
273 "Fry",
274 "Frz",
275 "Frax",
276 "Fray",
277 "Fraz",
278 "DISSIP",
279 "UUIv",
280 "VVIv",
281 "WWIv",
282 "UVIv",
283 "VWIv",
284 "UWIv",
285 "UUIbv",
286 "VVIbv",
287 "WWIbv",
288 "UUIb",
289 "VVIb",
290 "WWIb",
291 "UUUIb",
292 "UUVIb",
293 "UUWIb",
294 "UVVIb",
295 "UWWIb",
296 "VVVIb",
297 "VVWIb",
298 "VWWIb",
299 "WWWIb",
300 "AI",
301 "NK",
302 "IdUdxdUdxdUdx",
303 "IdUdydUdydUdx",
304 "IdUdzdUdzdUdx",
305 "IdUdxdVdxdUdy",
306 "IdUdydVdydUdy",
307 "IdUdzdVdzdUdy",
308 "IdUdxdWdxdUdz",
309 "IdUdydWdydUdz",
310 "IdUdzdWdzdUdz",
311 "IdVdxdUdxdVdx",
312 "IdVdydUdydVdx",
313 "IdVdzdUdzdVdx",
314 "IdVdxdVdxdVdy",
315 "IdVdydVdydVdy",
316 "IdVdzdVdzdVdy",
317 "IdVdxdWdxdVdz",
318 "IdVdydWdydVdz",
319 "IdVdzdWdzdVdz",
320 "IdWdxdUdxdWdx",
321 "IdWdydUdydWdx",
322 "IdWdzdUdzdWdx",
323 "IdWdxdVdxdWdy",
324 "IdWdydVdydWdy",
325 "IdWdzdVdzdWdy",
326 "IdWdxdWdxdWdz",
327 "IdWdydWdydWdz",
328 "IdWdzdWdzdWdz",
329 "IdUdxdUdxW",
330 "IdUdydUdyW",
331 "IdUdzdUdzW",
332 "IdVdxdVdxW",
333 "IdVdydVdyW",
334 "IdVdzdVdzW",
335 "IdWdxdWdxW",
336 "IdWdydWdyW",
337 "IdWdzdWdzW",
338 "IdUdxddPdxdx",
339 "IdUdyddPdxdy",
340 "IdUdzddPdxdz",
341 "IdVdxddPdydx",
342 "IdVdyddPdydy",
343 "IdVdzddPdydz",
344 "IdWdxddPdzdx",
345 "IdWdyddPdzdy",
346 "IdWdzddPdzdz",
347 "IddUdxdxddUdxdx",
348 "IddUdxdyddUdxdy",
349 "IddUdxdzddUdxdz",
350 "IddUdydxddUdydx",
351 "IddUdydyddUdydy",
352 "IddUdydzddUdydz",
353 "IddUdzdxddUdzdx",
354 "IddUdzdyddUdzdy",
355 "IddUdzdzddUdzdz",
356 "IddVdxdxddVdxdx",
357 "IddVdxdyddVdxdy",
358 "IddVdxdzddVdxdz",
359 "IddVdydxddVdydx",
360 "IddVdydyddVdydy",
361 "IddVdydzddVdydz",
362 "IddVdzdxddVdzdx",
363 "IddVdzdyddVdzdy",
364 "IddVdzdzddVdzdz",
365 "IddWdxdxddWdxdx",
366 "IddWdxdyddWdxdy",
367 "IddWdxdzddWdxdz",
368 "IddWdydxddWdydx",
369 "IddWdydyddWdydy",
370 "IddWdydzddWdydz",
371 "IddWdzdxddWdzdx",
372 "IddWdzdyddWdzdy",
373 "IddWdzdzddWdzdz",
374 "dIdxddUdxdxdUdx",
375 "dIdxddUdxdydUdy",
376 "dIdxddUdxdzdUdz",
377 "dIdyddUdydxdUdx",
378 "dIdyddUdydydUdy",
379 "dIdyddUdydzdUdz",
380 "dIdzddUdzdxdUdx",
381 "dIdzddUdzdydUdy",
382 "dIdzddUdzdzdUdz",
383 "dIdxddVdxdxdVdx",
384 "dIdxddVdxdydVdy",
385 "dIdxddVdxdzdVdz",
386 "dIdyddVdydxdVdx",
387 "dIdyddVdydydVdy",
388 "dIdyddVdydzdVdz",
389 "dIdzddVdzdxdVdx",
390 "dIdzddVdzdydVdy",
391 "dIdzddVdzdzdVdz",
392 "dIdxddWdxdxdWdx",
393 "dIdxddWdxdydWdy",
394 "dIdxddWdxdzdWdz",
395 "dIdyddWdydxdWdx",
396 "dIdyddWdydydWdy",
397 "dIdyddWdydzdWdz",
398 "dIdzddWdzdxdWdx",
399 "dIdzddWdzdydWdy",
400 "dIdzddWdzdzdWdz",
401 "dIdxddUdxdz",
402 "dIdyddUdydz",
403 "dIdzddUdzdz",
404 "dIdzdUdxdUdx",
405 "dIdzdUdydUdy",
406 "dIdzdUdzdUdz",
407 "dIdzdVdxdVdx",
408 "dIdzdVdydVdy",
409 "dIdzdVdzdVdz",
410 "dIdzdWdxdWdx",
411 "dIdzdWdydWdy",
412 "dIdzdWdzdWdz",
413 "Nx",
414 "Ny",
415 "Nz",
416 "UaiNx",
417 "VaiNx",
418 "WaiNx",
419 "UaiNy",
420 "VaiNy",
421 "WaiNy",
422 "UaiNz",
423 "VaiNz",
424 "WaiNz",
425 "UPaiNx",
426 "VPaiNx",
427 "WPaiNx",
428 "UPaiNy",
429 "VPaiNy",
430 "WPaiNy",
431 "UPaiNz",
432 "VPaiNz",
433 "WPaiNz",
434 "UUaiNx",
435 "VVaiNx",
436 "WWaiNx",
437 "UUaiNy",
438 "VVaiNy",
439 "WWaiNy",
440 "UUaiNz",
441 "VVaiNz",
442 "WWaiNz",
443 "UVaiNx",
444 "UWaiNx",
445 "VWaiNx",
446 "UVaiNy",
447 "UWaiNy",
448 "VWaiNy",
449 "UVaiNz",
450 "UWaiNz",
451 "VWaiNz",
452 "PaiNx",
453 "PaiNy",
454 "PaiNz",
455 "IaiNx",
456 "IaiNy",
457 "IaiNz",
458 "dUdxaiNx",
459 "dVdxaiNx",
460 "dWdxaiNx",
461 "dUdyaiNy",
462 "dVdyaiNy",
463 "dWdyaiNy",
464 "dUdzaiNz",
465 "dVdzaiNz",
466 "dWdzaiNz",
467 "UdUdxaiNx",
468 "UdVdxaiNx",
469 "UdWdxaiNx",
470 "UdUdyaiNy",
471 "UdVdyaiNy",
472 "UdWdyaiNy",
473 "UdUdzaiNz",
474 "UdVdzaiNz",
475 "UdWdzaiNz",
476 "VdUdxaiNx",
477 "VdVdxaiNx",
478 "VdWdxaiNx",
479 "VdUdyaiNy",
480 "VdVdyaiNy",
481 "VdWdyaiNy",
482 "VdUdzaiNz",
483 "VdVdzaiNz",
484 "VdWdzaiNz",
485 "WdUdxaiNx",
486 "WdVdxaiNx",
487 "WdWdxaiNx",
488 "WdUdyaiNy",
489 "WdVdyaiNy",
490 "WdWdyaiNy",
491 "WdUdzaiNz",
492 "WdVdzaiNz",
493 "WdWdzaiNz",
494 "aiNxddUdxdxdUdx",
495 "aiNxddUdxdydUdy",
496 "aiNxddUdxdzdUdz",
497 "aiNyddUdydydUdy",
498 "aiNyddUdydzdUdz",
499 "aiNzddUdzdzdUdz",
500 "aiNxddVdxdxdVdx",
501 "aiNxddVdxdydVdy",
502 "aiNxddVdxdzdVdz",
503 "aiNyddVdydydVdy",
504 "aiNyddVdydzdVdz",
505 "aiNzddVdzdzdVdz",
506 "aiNxddWdxdxdWdx",
507 "aiNxddWdxdydWdy",
508 "aiNxddWdxdzdWdz",
509 "aiNyddWdydydWdy",
510 "aiNyddWdydzdWdz",
511 "aiNzddWdzdzdWdz",
512 "aiNxddUdxdz",
513 "aiNyddUdydz",
514 "aiNzddUdzdz",
515 "aiNzdUdxdUdx",
516 "aiNzdUdydUdy",
517 "aiNzdUdzdUdz",
518 "aiNzdVdxdVdx",
519 "aiNzdVdydVdy",
520 "aiNzdVdzdVdz",
521 "aiNzdWdxdWdx",
522 "aiNzdWdydWdy",
523 "aiNzdWdzdWdz",
524 "aiNx",
525 "aiNy",
526 "aiNz",
527 "kaiNx",
528 "kaiNy",
529 "kaiNz",
530 "aiNxx",
531 "aiNxy",
532 "aiNxz",
533 "aiNyy",
534 "aiNyz",
535 "aiNzz",
536 "aiNNxxxx",
537 "aiNNxxxy",
538 "aiNNxxxz",
539 "aiNNxxyy",
540 "aiNNxxyz",
541 "aiNNxxzz",
542 "aiNNxyyy",
543 "aiNNxyyz",
544 "aiNNxyzz",
545 "aiNNxzyy",
546 "aiNNxzyz",
547 "aiNNxzzz",
548 "aiNNyyyy",
549 "aiNNyyyz",
550 "aiNNyyzz",
551 "aiNNyzzz",
552 "aiNNzzzz",
553 "kai",
554// };
555// const char *noms_moyennes_prov_fixe[] =
556// {
557 "UI_INT",
558 "VI_INT",
559 "WI_INT",
560 "UI_INTUI_INT",
561 "UI_INTVI_INT",
562 "UI_INTWI_INT",
563 "VI_INTVI_INT",
564 "VI_INTWI_INT",
565 "WI_INTWI_INT",
566 "UI_INTUI_INTUI_INT",
567 "UI_INTUI_INTVI_INT",
568 "UI_INTUI_INTWI_INT",
569 "UI_INTVI_INTUI_INT",
570 "UI_INTVI_INTVI_INT",
571 "UI_INTVI_INTWI_INT",
572 "UI_INTWI_INTUI_INT",
573 "UI_INTWI_INTVI_INT",
574 "UI_INTWI_INTWI_INT",
575 "VI_INTUI_INTUI_INT",
576 "VI_INTUI_INTVI_INT",
577 "VI_INTUI_INTWI_INT",
578 "VI_INTVI_INTUI_INT",
579 "VI_INTVI_INTVI_INT",
580 "VI_INTVI_INTWI_INT",
581 "VI_INTWI_INTUI_INT",
582 "VI_INTWI_INTVI_INT",
583 "VI_INTWI_INTWI_INT",
584 "WI_INTUI_INTUI_INT",
585 "WI_INTUI_INTVI_INT",
586 "WI_INTUI_INTWI_INT",
587 "WI_INTVI_INTUI_INT",
588 "WI_INTVI_INTVI_INT",
589 "WI_INTVI_INTWI_INT",
590 "WI_INTWI_INTUI_INT",
591 "WI_INTWI_INTVI_INT",
592 "WI_INTWI_INTWI_INT",
593 "UI_INTUIUI",
594 "UI_INTUIVI",
595 "UI_INTUIWI",
596 "UI_INTVIUI",
597 "UI_INTVIVI",
598 "UI_INTVIWI",
599 "UI_INTWIUI",
600 "UI_INTWIVI",
601 "UI_INTWIWI",
602 "VI_INTUIUI",
603 "VI_INTUIVI",
604 "VI_INTUIWI",
605 "VI_INTVIUI",
606 "VI_INTVIVI",
607 "VI_INTVIWI",
608 "VI_INTWIUI",
609 "VI_INTWIVI",
610 "VI_INTWIWI",
611 "WI_INTUIUI",
612 "WI_INTUIVI",
613 "WI_INTUIWI",
614 "WI_INTVIUI",
615 "WI_INTVIVI",
616 "WI_INTVIWI",
617 "WI_INTWIUI",
618 "WI_INTWIVI",
619 "WI_INTWIWI",
620 "VI_INTP_INT",
621 "WI_INTP_INT",
622 "dUINTdxdUINTdx",
623 "dUINTdxdUINTdy",
624 "dUINTdxdUINTdz",
625 "dUINTdydUINTdy",
626 "dUINTdydUINTdz",
627 "dUINTdzdUINTdz",
628 "dUINTdxdVINTdx",
629 "dUINTdxdVINTdy",
630 "dUINTdxdVINTdz",
631 "dUINTdydVINTdy",
632 "dUINTdydVINTdz",
633 "dUINTdzdVINTdz",
634 "dUINTdxdWINTdx",
635 "dUINTdxdWINTdy",
636 "dUINTdxdWINTdz",
637 "dUINTdydWINTdy",
638 "dUINTdydWINTdz",
639 "dUINTdzdWINTdz",
640 "dVINTdxdUINTdx",
641 "dVINTdxdUINTdy",
642 "dVINTdxdUINTdz",
643 "dVINTdydUINTdy",
644 "dVINTdydUINTdz",
645 "dVINTdzdUINTdz",
646 "dVINTdxdVINTdx",
647 "dVINTdxdVINTdy",
648 "dVINTdxdVINTdz",
649 "dVINTdydVINTdy",
650 "dVINTdydVINTdz",
651 "dVINTdzdVINTdz",
652 "dVINTdxdWINTdx",
653 "dVINTdxdWINTdy",
654 "dVINTdxdWINTdz",
655 "dVINTdydWINTdy",
656 "dVINTdydWINTdz",
657 "dVINTdzdWINTdz",
658 "dWINTdxdUINTdx",
659 "dWINTdxdUINTdy",
660 "dWINTdxdUINTdz",
661 "dWINTdydUINTdy",
662 "dWINTdydUINTdz",
663 "dWINTdzdUINTdz",
664 "dWINTdxdVINTdx",
665 "dWINTdxdVINTdy",
666 "dWINTdxdVINTdz",
667 "dWINTdydVINTdy",
668 "dWINTdydVINTdz",
669 "dWINTdzdVINTdz",
670 "dWINTdxdWINTdx",
671 "dWINTdxdWINTdy",
672 "dWINTdxdWINTdz",
673 "dWINTdydWINTdy",
674 "dWINTdydWINTdz",
675 "dWINTdzdWINTdz",
676 "P_INTdUINTdx",
677 "P_INTdUINTdy",
678 "P_INTdUINTdz",
679 "P_INTdVINTdx",
680 "P_INTdVINTdy",
681 "P_INTdVINTdz",
682 "P_INTdWINTdx",
683 "P_INTdWINTdy",
684 "P_INTdWINTdz",
685 "dissip_int",
686 "P_NOPERTURBE",
687 "U_NOPERTURBE",
688 "V_NOPERTURBE",
689 "W_NOPERTURBE",
690 "I_NP",
691 "UI_INTP_INT",
692 // pression
693 "PNP",
694 "PIUNP",
695 "PIVNP",
696 "PIWNP",
697 "IPdUdxNP",
698 "IPdVdxNP",
699 "IPdWdxNP",
700 "IPdUdyNP",
701 "IPdVdyNP",
702 "IPdWdyNP",
703 "IPdUdzNP",
704 "IPdVdzNP",
705 "IPdWdzNP",
706 "IdPdxNP",
707 "IdPdyNP",
708 "IdPdzNP",
709 "PIseuil",
710 "PIUseuil",
711 "PIVseuil",
712 "PIWseuil",
713 "IPdUdxseuil",
714 "IPdVdxseuil",
715 "IPdWdxseuil",
716 "IPdUdyseuil",
717 "IPdVdyseuil",
718 "IPdWdyseuil",
719 "IPdUdzseuil",
720 "IPdVdzseuil",
721 "IPdWdzseuil",
722 "IdPdxseuil",
723 "IdPdyseuil",
724 "IdPdzseuil",
725 "PNPseuil",
726 "PIUNPseuil",
727 "PIVNPseuil",
728 "PIWNPseuil",
729 "IPdUdxNPseuil",
730 "IPdVdxNPseuil",
731 "IPdWdxNPseuil",
732 "IPdUdyNPseuil",
733 "IPdVdyNPseuil",
734 "IPdWdyNPseuil",
735 "IPdUdzNPseuil",
736 "IPdVdzNPseuil",
737 "IPdWdzNPseuil",
738 "IdPdxNPseuil",
739 "IdPdyNPseuil",
740 "IdPdzNPseuil",
741 "Ivrappelx",
742 "Ivrappely",
743 "Ivrappelz",
744 "IP_INT",
745 "UIvrappelx",
746 "VIvrappelx",
747 "WIvrappelx",
748 "UIvrappely",
749 "VIvrappely",
750 "WIvrappely",
751 "UIvrappelz",
752 "VIvrappelz",
753 "WIvrappelz",
754
755 // termini aggiuntivi per il calcolo della parte deviatorica del tensore degli sforzi
756 "dUdyaiNx",
757 "dUdzaiNx",
758 "dVdxaiNy",
759 "dVdzaiNy",
760 "dWdxaiNz",
761 "dWdyaiNz",
762
763 // Viscous stress in the vapor phase
764 "dUdxIv",
765 "dUdyIv",
766 "dUdzIv",
767 "dVdxIv",
768 "dVdyIv",
769 "dVdzIv",
770 "dWdxIv",
771 "dWdyIv",
772 "dWdzIv",
773
774 // Attempt for the molecular diffusion
775 "IddUdxdxU",
776 "IddVdxdxU",
777 "IddWdxdxU",
778 "IddUdxdxV",
779 "IddVdxdxV",
780 "IddWdxdxV",
781 "IddUdxdxW",
782 "IddVdxdxW",
783 "IddWdxdxW",
784 "IddUdydyU",
785 "IddVdydyU",
786 "IddWdydyU",
787 "IddUdydyV",
788 "IddVdydyV",
789 "IddWdydyV",
790 "IddUdydyW",
791 "IddVdydyW",
792 "IddWdydyW",
793 "IddUdzdzU",
794 "IddVdzdzU",
795 "IddWdzdzU",
796 "IddUdzdzV",
797 "IddVdzdzV",
798 "IddWdzdzV",
799 "IddUdzdzW",
800 "IddVdzdzW",
801 "IddWdzdzW",
802
803 // To deal with the numerical pressure ..
804 "Ix",
805 "xaiNx",
806 "xaiNy",
807 "xaiNz",
808 "UIx",
809 "VIx",
810 "WIx",
811// Redistribution corrective term
812 "dUdxIx",
813 "dUdyIx",
814 "dUdzIx",
815 "dVdxIx",
816 "dVdyIx",
817 "dVdzIx",
818 "dWdxIx",
819 "dWdyIx",
820 "dWdzIx",
821// Interfacial terms
822 "xUaiNx",
823 "xVaiNx",
824 "xWaiNx",
825 "xUaiNy",
826 "xVaiNy",
827 "xWaiNy",
828 "xUaiNz",
829 "xVaiNz",
830 "xWaiNz",
831 //extended pressure
832
833 "P_VAP_Iv",
834 "P_VAP_aiNx",
835 "P_VAP_aiNy",
836 "P_VAP_aiNz",
837 "P_LIQ_I",
838 "P_LIQ_aiNx",
839 "P_LIQ_aiNy",
840 "P_LIQ_aiNz",
841 "UP_LIQ_aiNx",
842 "VP_LIQ_aiNx",
843 "WP_LIQ_aiNx",
844 "UP_LIQ_aiNy",
845 "VP_LIQ_aiNy",
846 "WP_LIQ_aiNy",
847 "UP_LIQ_aiNz",
848 "VP_LIQ_aiNz",
849 "WP_LIQ_aiNz",
850 "IP_LIQ_dUdx",
851 "IP_LIQ_dUdy",
852 "IP_LIQ_dUdz",
853 "IP_LIQ_dVdx",
854 "IP_LIQ_dVdy",
855 "IP_LIQ_dVdz",
856 "IP_LIQ_dWdx",
857 "IP_LIQ_dWdy",
858 "IP_LIQ_dWdz",
859 "UP_LIQ_I",
860 "VP_LIQ_I",
861 "WP_LIQ_I",
862
863 // GAB DISSIP
864 "TRUE_DISSIP_MOY",
865 "DISSIP_VAP_MOY",
866 "TRUE_DISSIP_VAP_MOY",
867 "POWER_INTERF"
868 //"IdP_LIQ_dx",
869 //"IdP_LIQ_dy",
870 //"IdP_LIQ_dz",
871 //"IUdP_LIQ_dx",
872 //"IUdP_LIQ_dy",
873 //"IUdP_LIQ_dz",
874 //"IVdP_LIQ_dx",
875 //"IVdP_LIQ_dy",
876 //"IVdP_LIQ_dz",
877 //"IWdP_LIQ_dx",
878 //"IWdP_LIQ_dy",
879 //"IWdP_LIQ_dz"
880
881 };
882
883 // int nb_val_1=502;
884 //int nb_val_2=196; // Bulles fixes.
885 //nval_=nb_val_1+nb_val_2;
886
887 // FONCTIONNE
888// nval_=502+196+27+6+9+25+29;//+6;
889 // A VOIR
890 nval_=502+196+27+6+9+25+29+3+1;
891
892 noms_moyennes_.dimensionner_force(nval_);
893 for (int i=0; i<nval_; i++)
894 noms_moyennes_[i]=noms_moyennes_prov[i];
895
896 const char *noms_moyennes_provt[] =
897 {
898 "TI",
899 "TTI",
900 "ITU",
901 "ITV",
902 "ITW",
903 "ITUU",
904 "ITUV",
905 "ITUW",
906 "ITVV",
907 "ITVW",
908 "ITWW",
909 "ITdPdx",
910 "ITdPdy",
911 "ITdPdz",
912 "ITddUdxdx",
913 "ITddUdydy",
914 "ITddUdzdz",
915 "ITddVdxdx",
916 "ITddVdydy",
917 "ITddVdzdz",
918 "ITddWdxdx",
919 "ITddWdydy",
920 "ITddWdzdz",
921 "IUddTdxdx",
922 "IUddTdydy",
923 "IUddTdzdz",
924 "IVddTdxdx",
925 "IVddTdydy",
926 "IVddTdzdz",
927 "IWddTdxdx",
928 "IWddTdydy",
929 "IWddTdzdz",
930 "ITbulles"
931 };
932 nvalt_ = 33;
933 noms_moyennes_temperature_.dimensionner_force(nvalt_);
934 for (int i=0; i<nvalt_; i++)
935 noms_moyennes_temperature_[i]=noms_moyennes_provt[i];
936
937}
938
939
941{
942 ref_ijk_ft_ = ijk_ft;
943}
944
945
946// Les champs doivent tous etre sur le domaine ns normalement...
947// y compris ceux-ci :
948// o field_ai
949// o field_kappa_ai
950// o normale_cell
951// On pourrait avoir des champs ft tant qu'on ne les multiplie pas avec des champs de splitting different.
952// Mais ce n'est pas le cas car on veut multiplier par P et U...
953
954double Calculer_valeur_seuil(double p_seuil, double pressionijk, double p_seuil_max, double p_seuil_min)
955{
956 if (pressionijk>p_seuil_max)
957 {
958 p_seuil = p_seuil_max;
959 }
960 else if (pressionijk<p_seuil_min)
961 {
962 p_seuil = p_seuil_min;
963 }
964 else
965 {
966 p_seuil = pressionijk;
967 }
968 return p_seuil;
969}
970
972{
973 Navier_Stokes_FTD_IJK& ns = cas.eq_ns();
974
975 IJK_Field_vector3_double& vitesse=ns.velocity_;
976 // GR262753 : travail des forces d'interfaces
977 IJK_Field_vector3_double& force_interfaces=ns.terme_source_interfaces_ft_;
978 const IJK_Field_vector3_double& repulsion_interfaces=ns.terme_repulsion_interfaces_ft_;
979 // ---
980 IJK_Field_double& extended_pressure_liq= (cas.get_post().extended_pressure_computed_) ? cas.get_post().extended_pl_: ns.pressure_;
981 IJK_Field_double& extended_pressure_vap= (cas.get_post().extended_pressure_computed_) ? cas.get_post().extended_pv_: ns.pressure_;
982 IJK_Field_double& pression=ns.pressure_;
983 const IJK_Field_double& indicatrice=cas.get_interface().I();
984
985 // Nombre total de mailles en K
986 const int nktot = pression.get_domaine().get_nb_items_global(Domaine_IJK::ELEM, DIRECTION_K);
987
988// For each slice in the domain, the integral is divided by the length of the domain
989// facteur 1./(ni*nj) car sommation de ni*nj valeurs sur des mailles de meme taille
990// OU facteur delta_z / taille_totale_en_z si mailles non uniformes en z
991 const int nijtot = pression.get_domaine().get_nb_elem_tot(DIRECTION_I)
992 * pression.get_domaine().get_nb_elem_tot(DIRECTION_J);
993 double facteur = 1./(double)(nijtot);
994
995 // Nombre local de mailles en K
996 const int imax = pression.ni();
997 const int jmax = pression.nj();
998 const int kmax = pression.nk();
999 const int offset = pression.get_domaine().get_offset_local(DIRECTION_K);
1000 const int periok = pression.get_domaine().get_periodic_flag(DIRECTION_K);
1001
1002 const bool dipha = (!Option_IJK::DISABLE_DIPHASIQUE);
1003 IJK_Field_local_double zero, zeros;
1004 zero.allocate(imax, jmax, kmax, 1 /* ghost pour le grad */, 1. /* additional_layers */, 1 /*nb compo*/);
1005 zeros.allocate(imax, jmax, kmax, 1 /* ghost pour le grad */, 1. /* additional_layers */, 3 /*nb compo*/);
1006 //allocate_velocity()
1007
1008 IJK_Field_vector3_double& gradI=cas.get_post().grad_I_ns_;
1009 IJK_Field_vector3_double& sourceI=ns.terme_source_interfaces_ns_;
1010 IJK_Field_vector3_double& repuls=ns.terme_repulsion_interfaces_ns_;
1011 IJK_Field_vector3_double& absrepuls=ns.terme_abs_repulsion_interfaces_ns_;
1012 IJK_Field_double& field_ai=cas.get_post().ai_ns_;
1013 IJK_Field_double& field_kappa_ai=cas.get_post().kappa_ai_ns_;
1014 IJK_Field_vector3_double& normale_cell=cas.get_post().normale_cell_ns_;
1015
1016 // Compute pressure gradient :
1017 IJK_Field_vector3_double& gradP=vect_post_fields_.at("dPd");
1018 for (int dir = 0; dir < 3; dir++)
1019 gradP[dir].data() = 0.;
1020
1021 add_gradient_times_constant(pression, 1. /*constant*/, gradP[0], gradP[1], gradP[2]);
1022 for (int dir = 0; dir < 3; dir++)
1023 gradP[dir].echange_espace_virtuel(1);
1024
1025 IJK_Field_vector3_double& gradU=vect_post_fields_.at("dUd");
1026 IJK_Field_vector3_double& gradV=vect_post_fields_.at("dVd");
1027 IJK_Field_vector3_double& gradW=vect_post_fields_.at("dWd");
1028 // TODO: Help for GB : check_stats_ ?
1029 IJK_Field_vector3_double& grad2Pi = vect_post_fields_.at("ddPdd");
1030 IJK_Field_vector3_double& grad2Pc = vect_post_fields_.at("ddPddc");
1031 IJK_Field_vector3_double& grad2Ui = vect_post_fields_.at("ddUdd");
1032 IJK_Field_vector3_double& grad2Uc = vect_post_fields_.at("ddUddc");
1033 IJK_Field_vector3_double& grad2Vi = vect_post_fields_.at("ddVdd");
1034 IJK_Field_vector3_double& grad2Vc = vect_post_fields_.at("ddVddc");
1035 IJK_Field_vector3_double& grad2Wi = vect_post_fields_.at("ddWdd");
1036 IJK_Field_vector3_double& grad2Wc = vect_post_fields_.at("ddWddc");
1037 const Motcles& liste_post_instantanes = ref_ijk_ft_->get_post().get_liste_post_instantanes();
1038
1039 double coef_immobilisation_=ns.coef_immobilisation_;
1040 IJK_Field_double& indicatrice_np=cas.get_post().indicatrice_non_perturbe_;
1041 IJK_Field_vector3_double& integral_vitesse=cas.get_post().integrated_velocity_;
1042 IJK_Field_double& integral_pression=cas.get_post().integrated_pressure_;
1043 IJK_Field_vector3_double force_rappel=ns.force_rappel_;
1044 double p_seuil_max = ns.p_seuil_max_;
1045 double p_seuil_min = ns.p_seuil_min_;
1046 const IJK_Field_double& integral_vitesse_i = integral_vitesse[0];
1047 const IJK_Field_double& integral_vitesse_j = integral_vitesse[1];
1048 const IJK_Field_double& integral_vitesse_k = integral_vitesse[2];
1049
1050 ArrOfDouble coord_i, coord_j, coord_k;
1051 build_local_coords(pression, coord_i, coord_j, coord_k);
1052
1053 if (elem_coord_.size_array() == 0)
1054 {
1055 Cerr << "Erreur dans Statistiques_dns_ijk::update_stat: non initialise" << finl;
1056 Process::exit();
1057 }
1058
1059 const IJK_Field_double& vitesse_i = vitesse[0];
1060 const IJK_Field_double& vitesse_j = vitesse[1];
1061 const IJK_Field_double& vitesse_k = vitesse[2];
1062
1063 // GR262753 : travail forces interfaces
1065 compute_vecA_minus_vecB_in_vecA(force_interfaces, repulsion_interfaces);
1066
1067
1068 // Calcul le gradient de U aux cellules a partir de la vitesse aux faces :
1069 compute_and_store_gradU_cell(vitesse_i, vitesse_j, vitesse_k);
1070
1071#ifdef STAT_VERBOSE
1072 // Pour verifications :
1073 double erreur, erreur_ddP, divU, laplP, d_divU_dx, d_divU_dy, d_divU_dz;
1074 double L1_erreur=0., L1_erreur_ddP=0., L1_divU=0., L1_laplP=0., L1_d_divU_dx=0., L1_d_divU_dy=0., L1_d_divU_dz=0.;
1075 double L2_erreur=0., L2_erreur_ddP=0., L2_divU=0., L2_laplP=0., L2_d_divU_dx=0., L2_d_divU_dy=0., L2_d_divU_dz=0.;
1076#endif
1077
1078 DoubleTab tmp(nktot, nval_);
1079
1080 //Invalid value for the interpolation
1081 const double p_invalid = 1.e19;
1082 //Number of invalid cells on a single slice
1083 IntTab occ(nktot,2); // 2 cols : 0:liq, 1:vap
1084 occ=0;
1085 DoubleTab tmpt(nktot, nvalt_,nb_thermal_fields_);
1086 //Here starts the loop in space
1087 for (int k = 0; k < kmax; k++)
1088 {
1089 const double dz = tab_dz_[k+offset];
1090 bool on_the_first_cell = false;
1091 bool on_the_last_cell = false;
1092 const int kglob=k + offset;
1093 if (!periok)
1094 {
1095 if (kglob == 0)
1096 on_the_first_cell = true;
1097 if (kglob == nktot - 1)
1098 on_the_last_cell = true;
1099 }
1100 // Calcul des moyennes spatiales sur le plan ij
1101
1102 // On y stocke la somme de toutes les valeurs sur le plan ij, on divisera apres
1103 ArrOfDouble moy(nval_);
1104 for (int i = 0; i < nval_; i++)
1105 {
1106 moy[i] = 0.;
1107 }
1108 DoubleTab moyv(nvalt_,nb_thermal_fields_);
1109 moyv = 0.; // Table initialization.
1110 for (int j = 0; j < jmax; j++)
1111 {
1112 for (int i = 0; i < imax; i++)
1113 {
1114 double x = coord_i[i]; //This function find the x-coordinate of the assigned point of the plane
1115 // Vitesses au centre de la maille i,j,k
1116 // On peut ici interpoler plus finement si on veut:
1117 double u = (vitesse_i(i,j,k) + vitesse_i(i+1, j, k)) * 0.5;
1118 double v = (vitesse_j(i,j,k) + vitesse_j(i, j+1, k)) * 0.5;
1119 double w = (vitesse_k(i,j,k) + vitesse_k(i, j, k+1)) * 0.5;
1120 double p = pression(i,j,k);
1121 double pl = extended_pressure_liq(i,j,k);
1122 if (extended_pressure_liq(i,j,k) > p_invalid)
1123 {
1124 pl = 0;
1125 occ(kglob,0) +=1;
1126 }
1127 double pv = extended_pressure_vap(i,j,k);
1128 if (extended_pressure_vap(i,j,k) > p_invalid)
1129 {
1130 pv = 0;
1131 occ(kglob,1) +=1;
1132 }
1133 double chi = indicatrice(i,j,k);
1134 double chiv = 1. - chi;
1135
1136 // Gradient de pression au centre de la maille i,j,k
1137 double dPdx = (gradP[0](i,j,k) + gradP[0](i+1, j, k)) * 0.5;
1138 double dPdy = (gradP[1](i,j,k) + gradP[1](i, j+1, k)) * 0.5;
1139 double dPdz = (gradP[2](i,j,k) + gradP[2](i, j, k+1)) * 0.5;
1140
1141 // Produit vitesse par gradient de pression au centre de la maille i,j,k
1142 double UdPdx = (vitesse_i(i,j,k)*gradP[0](i,j,k) + vitesse_i(i+1, j, k)*gradP[0](i+1, j, k)) * 0.5;
1143 double VdPdy = (vitesse_j(i,j,k)*gradP[1](i,j,k) + vitesse_j(i, j+1, k)*gradP[1](i, j+1, k)) * 0.5;
1144 double WdPdz = (vitesse_k(i,j,k)*gradP[2](i,j,k) + vitesse_k(i, j, k+1)*gradP[2](i, j, k+1)) * 0.5;
1145
1146 // Produits vitesse/vitesse au centre de la maille i,j,k
1147 // (pour la partie symetrique du tenseur uu, on regarde une autre interpolation) :
1148 // ((ui+uj)/2)**2 != (ui**2 + uj**2)/2
1149 double uu = (vitesse_i(i,j,k)*vitesse_i(i,j,k) + vitesse_i(i+1, j, k)*vitesse_i(i+1, j, k)) * 0.5;
1150 double vv = (vitesse_j(i,j,k)*vitesse_j(i,j,k) + vitesse_j(i, j+1, k)*vitesse_j(i, j+1, k)) * 0.5;
1151 double ww = (vitesse_k(i,j,k)*vitesse_k(i,j,k) + vitesse_k(i, j, k+1)*vitesse_k(i, j, k+1)) * 0.5;
1152
1153 double uuu = (vitesse_i(i,j,k)*vitesse_i(i,j,k)*vitesse_i(i,j,k) + vitesse_i(i+1, j, k)*vitesse_i(i+1, j, k)*vitesse_i(i+1, j, k)) * 0.5;
1154 double vvv = (vitesse_j(i,j,k)*vitesse_j(i,j,k)*vitesse_j(i,j,k) + vitesse_j(i, j+1, k)*vitesse_j(i, j+1, k)*vitesse_j(i, j+1, k)) * 0.5;
1155 double www = (vitesse_k(i,j,k)*vitesse_k(i,j,k)*vitesse_k(i,j,k) + vitesse_k(i, j, k+1)*vitesse_k(i, j, k+1)*vitesse_k(i, j, k+1)) * 0.5;
1156
1157 // Gradient de vitesse au centre de la maille i,j,k
1158 double dUdx=0., dVdx=0., dWdx=0., dUdy=0., dVdy=0., dWdy=0., dUdz=0., dVdz=0., dWdz=0.;
1159 double ddUdxdx=0., ddVdxdx=0., ddWdxdx=0., ddUdydy=0., ddVdydy=0., ddWdydy=0., ddUdzdz=0., ddVdzdz=0., ddWdzdz=0.;
1160 double pseudo_dissip = calculer_gradients_vitesse(vitesse_i, vitesse_j, vitesse_k,
1161 i,j,k,
1162 dz,
1163 dUdx, dVdx, dWdx,
1164 dUdy, dVdy, dWdy,
1165 dUdz, dVdz, dWdz,
1166 ddUdxdx, ddVdxdx, ddWdxdx,
1167 ddUdydy, ddVdydy, ddWdydy,
1168 ddUdzdz, ddVdzdz, ddWdzdz,
1169 on_the_first_cell, on_the_last_cell);
1170 double true_dissip = calculer_vraie_dissipation(pseudo_dissip,
1171 dUdx, dUdy, dUdz,
1172 dVdx, dVdy, dVdz,
1173 dWdx, dWdy, dWdz);
1174 double travail_Force_interfaces_vitesse = 0.;
1176 {
1177 travail_Force_interfaces_vitesse = calculer_produit_scalaire_faces_to_center (
1178 vitesse_i, vitesse_j, vitesse_k,
1179 force_interfaces[0], force_interfaces[1], force_interfaces[2],
1180 i,j,k
1181 );
1182 }
1183
1184
1185 // Pour verifier un peu les stats ou pour le post-traitement, on stocke (systematiquement) les grad de vitesse :
1186 //if ((check_stats_) || (postraitement des gradients demande?))
1187 //{
1188 gradU[0](i,j,k) = dUdx;
1189 gradU[1](i,j,k) = dUdy;
1190 gradU[2](i,j,k) = dUdz;
1191 gradV[0](i,j,k) = dVdx;
1192 gradV[1](i,j,k) = dVdy;
1193 gradV[2](i,j,k) = dVdz;
1194 gradW[0](i,j,k) = dWdx;
1195 gradW[1](i,j,k) = dWdy;
1196 gradW[2](i,j,k) = dWdz;
1197 //}
1198#ifdef STAT_VERBOSE
1199 // Debug pour verifier que je m'emmele pas les pinceaux dans les variables :
1200 erreur = std::fabs(gradU[0](i,j,k) - dUdx)
1201 + std::fabs(gradV[1](i,j,k) - dVdy)
1202 + std::fabs(gradW[0](i,j,k) - dWdx)
1203 + std::fabs(gradU[2](i,j,k) - dUdz)
1204 + std::fabs(gradV[2](i,j,k) - dVdz)
1205 + std::fabs(gradW[2](i,j,k) - dWdz);
1206
1207 if (erreur > PRECISION_DERIVEES)
1208 {
1209 Cerr << "Statistiques_dns_ijk_FT::update_stat -- Erreur (gradU) superieure au seuil : "
1210 << i << " " << j << " " << k << " "
1211 << erreur
1212 << finl;
1213 Cerr << "dudx: " << field_dudx(i,j,k) << " " << dUdx << " " << field_dudx(i,j,k) - dUdx << finl;
1214 Cerr << "dvdy: " << field_dvdy(i,j,k) << " " << dVdy << " " << field_dvdy(i,j,k) - dVdy << finl;
1215 Cerr << "dwdx: " << field_dwdx(i,j,k) << " " << dWdx << " " << field_dwdx(i,j,k) - dWdx << finl;
1216 Cerr << "dudz: " << field_dudz(i,j,k) << " " << dUdz << " " << field_dudz(i,j,k) - dUdz << finl;
1217 Cerr << "dvdz: " << field_dvdz(i,j,k) << " " << dVdz << " " << field_dvdz(i,j,k) - dVdz << finl;
1218 Cerr << "dwdz: " << field_dwdz(i,j,k) << " " << dWdz << " " << field_dwdz(i,j,k) - dWdz << finl;
1219 Cerr << "On proc " << Process::me() << " Cell: " << i << " " << j << " " << k << finl;
1220#ifdef STAT_EXIT
1221 Process::exit();
1222#endif
1223 }
1224
1225 divU = dUdx + dVdy + dWdz;
1226 if (std::fabs(divU) > PRECISION_DIVU)
1227 {
1228 Cerr << "Statistiques_dns_ijk_FT::update_stat -- divU : " << divU << finl;
1229 Cerr << "dUdx,dVdy,dWdz : " << dUdx << " " << dVdy << " " << dWdz << finl;
1230 Cerr << "On proc " << Process::me() << " Cell: " << i << " " << j << " " << k << finl;
1231#ifdef STAT_EXIT
1232 Process::exit();
1233#endif
1234 }
1235
1236 // Compute L1 and L2 norms :
1237 L1_erreur += erreur;
1238 L2_erreur += erreur*erreur;
1239 L1_divU += std::fabs(divU);
1240 L2_divU += divU*divU;
1241#endif
1242
1243 // Second gradients croises :
1244 double ddUdxdy = 0.;
1245 double ddUdxdz = 0.;
1246 //double ddUdydx = ddUdxdy;
1247
1248 double ddUdydz = 0.;
1249 //double ddUdzdx = ddUdxdz;
1250 //double ddUdzdy = ddUdydz;
1251
1252 double ddVdxdy = 0.;
1253 double ddVdxdz = 0.;
1254 //double ddVdydx = ddVdxdy;
1255
1256 double ddVdydz = 0.;
1257 //double ddVdzdx = ddVdxdz;
1258 //double ddVdzdy = ddVdydz;
1259
1260 double ddWdxdy = 0.;
1261 double ddWdxdz = 0.;
1262 //double ddWdydx = ddWdxdy;
1263
1264 double ddWdydz = 0.;
1265 //double ddWdzdx = ddWdxdz;
1266 //double ddWdzdy = ddWdydz;
1268 gradU[0], gradV[1], gradW[0],
1269 gradU[2], gradV[2], gradW[2], /* field_dudz, field_dvdz, field_dwdz, */
1270 /* Et les outputs en ref aussi!! */
1271 ddUdxdy, ddUdxdz, ddUdydz,
1272 ddVdxdy, ddVdxdz, ddVdydz,
1273 ddWdxdy, ddWdxdz, ddWdydz);
1274
1275 // Pour verifier un peu les stats, on peut stocker le double grad :
1276 if (check_stats_)
1277 {
1278 grad2Ui[0](i,j,k) = ddUdxdx;
1279 grad2Ui[1](i,j,k) = ddUdydy;
1280 grad2Ui[2](i,j,k) = ddUdzdz;
1281 // Et les croisees :
1282 grad2Uc[0](i,j,k) = ddUdxdy;
1283 grad2Uc[1](i,j,k) = ddUdxdz;
1284 grad2Uc[2](i,j,k) = ddUdydz;
1285 //
1286 grad2Vi[0](i,j,k) = ddVdxdx;
1287 grad2Vi[1](i,j,k) = ddVdydy;
1288 grad2Vi[2](i,j,k) = ddVdzdz;
1289 // Et les croisees :
1290 grad2Vc[0](i,j,k) = ddVdxdy;
1291 grad2Vc[1](i,j,k) = ddVdxdz;
1292 grad2Vc[2](i,j,k) = ddVdydz;
1293 //
1294 grad2Wi[0](i,j,k) = ddWdxdx;
1295 grad2Wi[1](i,j,k) = ddWdydy;
1296 grad2Wi[2](i,j,k) = ddWdzdz;
1297 // Et les croisees :
1298 grad2Wc[0](i,j,k) = ddWdxdy;
1299 grad2Wc[1](i,j,k) = ddWdxdz;
1300 grad2Wc[2](i,j,k) = ddWdydz;
1301 //
1302 }
1303
1304 // Les croises non calcules sont en fait symetriques :
1305 double ddUdydx = ddUdxdy;
1306 double ddUdzdx = ddUdxdz;
1307 double ddUdzdy = ddUdydz;
1308 double ddVdydx = ddVdxdy;
1309 double ddVdzdx = ddVdxdz;
1310 double ddVdzdy = ddVdydz;
1311 double ddWdydx = ddWdxdy;
1312 double ddWdzdx = ddWdxdz;
1313 double ddWdzdy = ddWdydz;
1314#ifdef STAT_VERBOSE
1315 d_divU_dx = ddUdxdx + ddVdydx + ddWdzdx;
1316 d_divU_dy = ddUdxdy + ddVdydy + ddWdzdy;
1317 d_divU_dz = ddUdxdz + ddVdydz + ddWdzdz;
1318 if ((std::fabs(d_divU_dx) > PRECISION_DDIVU)
1319 || (std::fabs(d_divU_dy) > PRECISION_DDIVU)
1320 || (std::fabs(d_divU_dz) > PRECISION_DDIVU) )
1321 {
1322 Cerr << "Statistiques_dns_ijk_FT::update_stat -- grad(divU) : "
1323 << d_divU_dx << " "
1324 << d_divU_dy << " "
1325 << d_divU_dz << " "
1326 << finl;
1327 Cerr << "ddUdxdz,ddVdydz,ddWdzdz : " << ddUdxdz << " " << ddVdydz << " " << ddWdzdz << finl;
1328 Cerr << "On proc " << Process::me() << " Cell: " << i << " " << j << " " << k << finl;
1329#ifdef STAT_EXIT
1330 Process::exit();
1331#endif
1332 }
1333
1334 // Compute L1 and L2 norms :
1335 L1_d_divU_dx += std::fabs(d_divU_dx);
1336 L2_d_divU_dx += d_divU_dx*d_divU_dx;
1337 L1_d_divU_dy += std::fabs(d_divU_dy);
1338 L2_d_divU_dy += d_divU_dy*d_divU_dy;
1339 L1_d_divU_dz += std::fabs(d_divU_dz);
1340 L2_d_divU_dz += d_divU_dz*d_divU_dz;
1341#endif
1342
1343 // Derivee seconde de la pression :
1344 double ddPdxdx = 0.;
1345 double ddPdxdy = 0.;
1346 double ddPdxdz = 0.;
1347 double ddPdydx = 0.;
1348 double ddPdydy = 0.;
1349 double ddPdydz = 0.;
1350 double ddPdzdx = 0.;
1351 double ddPdzdy = 0.;
1352 double ddPdzdz = 0.;
1353 face_to_cell_gradient(gradP[0],gradP[1],gradP[2],
1354 i,j,k,
1355 dz,
1356 ddPdxdx,ddPdydx,ddPdzdx,
1357 ddPdxdy,ddPdydy,ddPdzdy,
1358 ddPdxdz,ddPdydz,ddPdzdz,
1359 on_the_first_cell, on_the_last_cell,
1360 1 /* bc_type for gradP at the wall :
1361 assumed equal to faces values */);
1362
1363 // Pour verifier un peu les stats, on peut stocker le double grad :
1364 if (check_stats_)
1365 {
1366 grad2Pi[0](i,j,k) = ddPdxdx;
1367 grad2Pi[1](i,j,k) = ddPdydy;
1368 grad2Pi[2](i,j,k) = ddPdzdz;
1369 // Et les croisees :
1370 grad2Pc[0](i,j,k) = ddPdxdy;
1371 grad2Pc[1](i,j,k) = ddPdxdz;// ou ddPdzdx (sont equivalents d'ordre 2)
1372 grad2Pc[2](i,j,k) = ddPdzdy;// ou ddPdzdy (sont equivalents d'ordre 2)
1373 }
1374
1375#ifdef STAT_VERBOSE
1376 // Debug pour verifier que je m'emmele pas les pinceaux dans les variables :
1377 erreur_ddP = std::fabs(ddPdydx - ddPdxdy)
1378 + std::fabs(ddPdzdx - ddPdxdz)
1379 + std::fabs(ddPdzdy - ddPdydz);
1380
1381 if (erreur_ddP > PRECISION_DERIVEES)
1382 {
1383 Cerr << "Statistiques_dns_ijk_FT::update_stat -- Erreur_ddP superieure au seuil : "
1384 << i << " " << j << " " << k << " "
1385 << erreur_ddP
1386 << finl;
1387 Cerr << "We should have symmetric behaviour for crossed derivatives : \n"
1388 << "ddPdydx: " << ddPdydx << " " << ddPdxdy << " diff: " << ddPdydx-ddPdxdy << "\n";
1389 Cerr << "ddPdzdx: " << ddPdzdx << " " << ddPdxdz << " diff: " << ddPdzdx-ddPdxdz << "\n"
1390 << "ddPdydz: " << ddPdydz << " " << ddPdzdy << " diff: " << ddPdydz-ddPdzdy << "\n"
1391 << finl;
1392 Cerr << "On proc " << Process::me() << " Cell: " << i << " " << j << " " << k << finl;
1393 Process::exit();
1394 }
1395
1396 laplP = ddPdxdx + ddPdydy + ddPdzdz;
1397 if (std::fabs(dIdx)+std::fabs(dIdy)+std::fabs(dIdz) < PRECISION_LAPLP )
1398 {
1399 // On a gradI nul sur toutes les faces donc la cellule est "franchement" monophasique.
1400 // On ne devrait donc pas avoir de surprise, rho aux faces devrait etre constant
1401 // et par consequent, on devrait avoir : laplP = 0
1402 if (std::fabs(laplP) > PRECISION_LAPLP)
1403 {
1404 Cerr << "Statistiques_dns_ijk_FT::update_stat -- Laplacien P : " << laplP << finl;
1405 Cerr << "On proc " << Process::me() << " Cell: " << i << " " << j << " " << k << finl;
1406#ifdef STAT_EXIT
1407 Process::exit();
1408#endif
1409 }
1410 }
1411
1412 // Compute L1 and L2 norms :
1413 L1_erreur_ddP += erreur_ddP;
1414 L2_erreur_ddP += erreur_ddP*erreur_ddP;
1415 L1_laplP += std::fabs(laplP);
1416 L2_laplP += laplP*laplP;
1417#endif
1418
1419 double ai=0., kai=0., Nx=0.,Ny=0.,Nz=0., aiNx=0.,aiNy=0.,aiNz=0., Nxx=0.,Nxy=0.,Nxz=0.;
1420 double Frx=0.,Fry=0.,Frz=0., Frax=0.,Fray=0.,Fraz=0., Fx=0.,Fy=0.,Fz=0., Nyy=0.,Nyz=0.,Nzz=0.;
1421 double dIdx=0., dIdy=0., dIdz=0.;
1422 if (dipha)
1423 {
1424 ai = field_ai(i,j,k);
1425 kai = field_kappa_ai(i,j,k);
1426
1427 // Normale a l'interface au centre de l'elem :
1428 Nx = normale_cell[0](i,j,k);
1429 Ny = normale_cell[1](i,j,k);
1430 Nz = normale_cell[2](i,j,k);
1431
1432 aiNx = ai*Nx;
1433 aiNy = ai*Ny;
1434 aiNz = ai*Nz;
1435
1436 Nxx = Nx*Nx;
1437 Nxy = Nx*Ny;
1438 Nxz = Nx*Nz;
1439 Nyy = Ny*Ny;
1440 Nyz = Ny*Nz;
1441 Nzz = Nz*Nz;
1442
1443 // Force de repulsion :
1444 Fx = (sourceI[0](i,j,k) + sourceI[0](i+1, j, k)) * 0.5;
1445 Fy = (sourceI[1](i,j,k) + sourceI[1](i, j+1, k)) * 0.5;
1446 Fz = (sourceI[2](i,j,k) + sourceI[2](i, j, k+1)) * 0.5;
1447
1448 // Force de repulsion :
1449 Frx = (repuls[0](i,j,k) + repuls[0](i+1, j, k)) * 0.5;
1450 Fry = (repuls[1](i,j,k) + repuls[1](i, j+1, k)) * 0.5;
1451 Frz = (repuls[2](i,j,k) + repuls[2](i, j, k+1)) * 0.5;
1452
1453 // Force de repulsion :
1454 Frax = (absrepuls[0](i,j,k) + absrepuls[0](i+1, j, k)) * 0.5;
1455 Fray = (absrepuls[1](i,j,k) + absrepuls[1](i, j+1, k)) * 0.5;
1456 Fraz = (absrepuls[2](i,j,k) + absrepuls[2](i, j, k+1)) * 0.5;
1457
1458 // Gradient de l'indicatrice au centre de la maille i,j,k
1459 dIdx = (gradI[0](i,j,k) + gradI[0](i+1, j, k)) * 0.5;
1460 dIdy = (gradI[1](i,j,k) + gradI[1](i, j+1, k)) * 0.5;
1461 dIdz = (gradI[2](i,j,k) + gradI[2](i, j, k+1)) * 0.5;
1462 }
1463
1464 if (ref_ijk_ft_->has_thermals())
1465 {
1466 int idx =0;
1467 for (auto& itr : ref_ijk_ft_->get_ijk_thermals().get_liste_eqs())
1468 {
1469 const IJK_Field_double& temperature = *itr->get_temperature();
1470 const double T = temperature(i,j,k);
1471
1472 double T_adim_bulles = 0.;
1473 if (liste_post_instantanes.contient_("TEMPERATURE_ADIM_BULLES"))
1474 {
1475 const IJK_Field_double& temperature_adim_bulles = itr->get_temperature_adim_bulles();
1476 T_adim_bulles = temperature_adim_bulles(i,j,k);
1477 }
1478
1479
1480 // Derivee seconde de la temperature :
1481 const IJK_Field_vector3_double& gradT = itr->get_gradient_temperature();
1482 double ddTdxdx = 0.;
1483 double ddTdxdy = 0.;
1484 double ddTdxdz = 0.;
1485 double ddTdydx = 0.;
1486 double ddTdydy = 0.;
1487 double ddTdydz = 0.;
1488 double ddTdzdx = 0.;
1489 double ddTdzdy = 0.;
1490 double ddTdzdz = 0.;
1491 face_to_cell_gradient(gradT[0],gradT[1],gradT[2],
1492 i,j,k,
1493 dz,
1494 ddTdxdx,ddTdydx,ddTdzdx,
1495 ddTdxdy,ddTdydy,ddTdzdy,
1496 ddTdxdz,ddTdydz,ddTdzdz,
1497 on_the_first_cell, on_the_last_cell,
1498 1 /* bc_type for gradT at the wall :
1499 assumed equal to faces values */);
1500#define AJOUTT(somme,val) moyv(somme,idx) += val
1501 AJOUTT(TI_MOY,chi*T);
1502 AJOUTT(TTI_MOY, chi*T*T);
1503 AJOUTT(ITU_MOY, chi*T*u);
1504 AJOUTT(ITV_MOY, chi*T*v);
1505 AJOUTT(ITW_MOY, chi*T*w);
1506 AJOUTT(ITUU_MOY, chi*T*u*u);
1507 AJOUTT(ITUV_MOY, chi*T*u*v);
1508 AJOUTT(ITUW_MOY, chi*T*u*w);
1509 AJOUTT(ITVV_MOY, chi*T*v*v);
1510 AJOUTT(ITVW_MOY, chi*T*v*w);
1511 AJOUTT(ITWW_MOY, chi*T*w*w);
1512 AJOUTT(ITdPdx_MOY, chi*T*dPdx);
1513 AJOUTT(ITdPdy_MOY, chi*T*dPdy);
1514 AJOUTT(ITdPdz_MOY, chi*T*dPdz);
1515 AJOUTT(ITddUdxdx_MOY, chi*T*ddUdxdx);
1516 AJOUTT(ITddUdydy_MOY, chi*T*ddUdydy);
1517 AJOUTT(ITddUdzdz_MOY, chi*T*ddUdzdz);
1518 AJOUTT(ITddVdxdx_MOY, chi*T*ddVdxdx);
1519 AJOUTT(ITddVdydy_MOY, chi*T*ddVdydy);
1520 AJOUTT(ITddVdzdz_MOY, chi*T*ddVdzdz);
1521 AJOUTT(ITddWdxdx_MOY, chi*T*ddWdxdx);
1522 AJOUTT(ITddWdydy_MOY, chi*T*ddWdydy);
1523 AJOUTT(ITddWdzdz_MOY, chi*T*ddWdzdz);
1524 AJOUTT(IUddTdxdx_MOY, chi*u*ddTdxdx);
1525 AJOUTT(IUddTdydy_MOY, chi*u*ddTdydy);
1526 AJOUTT(IUddTdzdz_MOY, chi*u*ddTdzdz);
1527 AJOUTT(IVddTdxdx_MOY, chi*v*ddTdxdx);
1528 AJOUTT(IVddTdydy_MOY, chi*v*ddTdydy);
1529 AJOUTT(IVddTdzdz_MOY, chi*v*ddTdzdz);
1530 AJOUTT(IWddTdxdx_MOY, chi*w*ddTdxdx);
1531 AJOUTT(IWddTdydy_MOY, chi*w*ddTdydy);
1532 AJOUTT(IWddTdzdz_MOY, chi*w*ddTdzdz);
1533 AJOUTT(ITbulles_MOY, chi*T_adim_bulles);
1534#undef AJOUTT
1535 idx++;
1536 }
1537 }
1538
1539#define AJOUT(somme,val) moy[somme] += val
1540 // moyennes
1541 AJOUT(I_MOY,chi);
1542 AJOUT(UI_MOY,u*chi);
1543 AJOUT(VI_MOY,v*chi);
1544 AJOUT(WI_MOY,w*chi);
1545 AJOUT(PI_MOY,p*chi);
1546 AJOUT(UIv_MOY,u*chiv);
1547 AJOUT(VIv_MOY,v*chiv);
1548 AJOUT(WIv_MOY,w*chiv);
1549 AJOUT(PIv_MOY,p*chiv);
1550 // correlations 2 eme ordre vitesse :
1551 AJOUT(UUI_MOY,u*u*chi);
1552 AJOUT(VVI_MOY,v*v*chi);
1553 AJOUT(WWI_MOY,w*w*chi);
1554 AJOUT(UVI_MOY,u*v*chi);
1555 AJOUT(VWI_MOY,v*w*chi);
1556 AJOUT(UWI_MOY,u*w*chi);
1557 // 3 eme ordre (10 termes)
1558 AJOUT(UUUI_MOY,u*u*u*chi);
1559 AJOUT(UUVI_MOY,u*u*v*chi);
1560 AJOUT(UUWI_MOY,u*u*w*chi);
1561 AJOUT(UVVI_MOY,u*v*v*chi);
1562 AJOUT(UVWI_MOY,u*v*w*chi);
1563 AJOUT(UWWI_MOY,u*w*w*chi);
1564 AJOUT(VVVI_MOY,v*v*v*chi);
1565 AJOUT(VVWI_MOY,v*v*w*chi);
1566 AJOUT(VWWI_MOY,v*w*w*chi);
1567 AJOUT(WWWI_MOY,w*w*w*chi);
1568 // correlation vitesse/pression :
1569 AJOUT(UPI_MOY,u*p*chi);
1570 AJOUT(VPI_MOY,v*p*chi);
1571 AJOUT(WPI_MOY,w*p*chi);
1572 // Maintenant, ce qui fait intervenir les derivees :
1573 // 1. de l'indicatrice :
1574 //^^^^^^^^^^^^^^^^^^^^^^^
1575 // ui*dI/dxb : (9)
1576 AJOUT(UdIdx_MOY,u*dIdx);
1577 AJOUT(VdIdx_MOY,v*dIdx);
1578 AJOUT(WdIdx_MOY,w*dIdx);
1579 AJOUT(UdIdy_MOY,u*dIdy);
1580 AJOUT(VdIdy_MOY,v*dIdy);
1581 AJOUT(WdIdy_MOY,w*dIdy);
1582 AJOUT(UdIdz_MOY,u*dIdz);
1583 AJOUT(VdIdz_MOY,v*dIdz);
1584 AJOUT(WdIdz_MOY,w*dIdz);
1585 //^^^^^^^^^^^^^^^^^^^^^^^
1586 // ui*p*dI/dxb : (9)
1587 AJOUT(UPdIdx_MOY,u*p*dIdx);
1588 AJOUT(VPdIdx_MOY,v*p*dIdx);
1589 AJOUT(WPdIdx_MOY,w*p*dIdx);
1590 AJOUT(UPdIdy_MOY,u*p*dIdy);
1591 AJOUT(VPdIdy_MOY,v*p*dIdy);
1592 AJOUT(WPdIdy_MOY,w*p*dIdy);
1593 AJOUT(UPdIdz_MOY,u*p*dIdz);
1594 AJOUT(VPdIdz_MOY,v*p*dIdz);
1595 AJOUT(WPdIdz_MOY,w*p*dIdz);
1596 //^^^^^^^^^^^^^^^^^^^^^^^
1597 // ui*uj*dI/dxb : (18)
1598 AJOUT(UUdIdx_MOY,u*u*dIdx);
1599 AJOUT(VVdIdx_MOY,v*v*dIdx);
1600 AJOUT(WWdIdx_MOY,w*w*dIdx);
1601 AJOUT(UUdIdy_MOY,u*u*dIdy);
1602 AJOUT(VVdIdy_MOY,v*v*dIdy);
1603 AJOUT(WWdIdy_MOY,w*w*dIdy);
1604 AJOUT(UUdIdz_MOY,u*u*dIdz);
1605 AJOUT(VVdIdz_MOY,v*v*dIdz);
1606 AJOUT(WWdIdz_MOY,w*w*dIdz);
1607 AJOUT(UVdIdx_MOY,u*v*dIdx);
1608 AJOUT(UWdIdx_MOY,u*w*dIdx);
1609 AJOUT(VWdIdx_MOY,v*w*dIdx);
1610 AJOUT(UVdIdy_MOY,u*v*dIdy);
1611 AJOUT(UWdIdy_MOY,u*w*dIdy);
1612 AJOUT(VWdIdy_MOY,v*w*dIdy);
1613 AJOUT(UVdIdz_MOY,u*v*dIdz);
1614 AJOUT(UWdIdz_MOY,u*w*dIdz);
1615 AJOUT(VWdIdz_MOY,v*w*dIdz);
1616 //^^^^^^^^^^^^^^^^^^^^^^
1617 // p*dI/dxb : (3)
1618 AJOUT(PdIdx_MOY,p*dIdx);
1619 AJOUT(PdIdy_MOY,p*dIdy);
1620 AJOUT(PdIdz_MOY,p*dIdz);
1621 //^^^^^^^^^^^^^^^^^^^^^^
1622 // I*dI/dxb : (3)
1623 AJOUT(IdIdx_MOY,chi*dIdx);
1624 AJOUT(IdIdy_MOY,chi*dIdy);
1625 AJOUT(IdIdz_MOY,chi*dIdz);
1626 // 2. de la vitesse :
1627 //^^^^^^^^^^^^^^^^^^^^^^^
1628 // I*dUidxb : (9)
1629 AJOUT(IdUdx_MOY,chi*dUdx);
1630 AJOUT(IdUdy_MOY,chi*dUdy);
1631 AJOUT(IdUdz_MOY,chi*dUdz);
1632 AJOUT(IdVdx_MOY,chi*dVdx);
1633 AJOUT(IdVdy_MOY,chi*dVdy);
1634 AJOUT(IdVdz_MOY,chi*dVdz);
1635 AJOUT(IdWdx_MOY,chi*dWdx);
1636 AJOUT(IdWdy_MOY,chi*dWdy);
1637 AJOUT(IdWdz_MOY,chi*dWdz);
1638 // I*P*dUidxb : (9)
1639 AJOUT(IPdUdx_MOY,chi*p*dUdx);
1640 AJOUT(IPdUdy_MOY,chi*p*dUdy);
1641 AJOUT(IPdUdz_MOY,chi*p*dUdz);
1642 AJOUT(IPdVdx_MOY,chi*p*dVdx);
1643 AJOUT(IPdVdy_MOY,chi*p*dVdy);
1644 AJOUT(IPdVdz_MOY,chi*p*dVdz);
1645 AJOUT(IPdWdx_MOY,chi*p*dWdx);
1646 AJOUT(IPdWdy_MOY,chi*p*dWdy);
1647 AJOUT(IPdWdz_MOY,chi*p*dWdz);
1648 // 3. de l'indicatrice et de la vitesse :
1649 //^^^^^^^^^^^^^^^^^^^^^^^
1650 // dUidxb*dIdxb : (9)
1651 AJOUT(dUdxdIdx_MOY,dUdx*dIdx);
1652 AJOUT(dUdydIdy_MOY,dUdy*dIdy);
1653 AJOUT(dUdzdIdz_MOY,dUdz*dIdz);
1654 AJOUT(dVdxdIdx_MOY,dVdx*dIdx);
1655 AJOUT(dVdydIdy_MOY,dVdy*dIdy);
1656 AJOUT(dVdzdIdz_MOY,dVdz*dIdz);
1657 AJOUT(dWdxdIdx_MOY,dWdx*dIdx);
1658 AJOUT(dWdydIdy_MOY,dWdy*dIdy);
1659 AJOUT(dWdzdIdz_MOY,dWdz*dIdz);
1660 // Ui*dUjdxb*dIdxb : (27)
1661 AJOUT(UdUdxdIdx_MOY,u*dUdx*dIdx);
1662 AJOUT(UdVdxdIdx_MOY,u*dVdx*dIdx);
1663 AJOUT(UdWdxdIdx_MOY,u*dWdx*dIdx);
1664 AJOUT(UdUdydIdy_MOY,u*dUdy*dIdy);
1665 AJOUT(UdVdydIdy_MOY,u*dVdy*dIdy);
1666 AJOUT(UdWdydIdy_MOY,u*dWdy*dIdy);
1667 AJOUT(UdUdzdIdz_MOY,u*dUdz*dIdz);
1668 AJOUT(UdVdzdIdz_MOY,u*dVdz*dIdz);
1669 AJOUT(UdWdzdIdz_MOY,u*dWdz*dIdz);
1670 AJOUT(VdUdxdIdx_MOY,v*dUdx*dIdx);
1671 AJOUT(VdVdxdIdx_MOY,v*dVdx*dIdx);
1672 AJOUT(VdWdxdIdx_MOY,v*dWdx*dIdx);
1673 AJOUT(VdUdydIdy_MOY,v*dUdy*dIdy);
1674 AJOUT(VdVdydIdy_MOY,v*dVdy*dIdy);
1675 AJOUT(VdWdydIdy_MOY,v*dWdy*dIdy);
1676 AJOUT(VdUdzdIdz_MOY,v*dUdz*dIdz);
1677 AJOUT(VdVdzdIdz_MOY,v*dVdz*dIdz);
1678 AJOUT(VdWdzdIdz_MOY,v*dWdz*dIdz);
1679 AJOUT(WdUdxdIdx_MOY,w*dUdx*dIdx);
1680 AJOUT(WdVdxdIdx_MOY,w*dVdx*dIdx);
1681 AJOUT(WdWdxdIdx_MOY,w*dWdx*dIdx);
1682 AJOUT(WdUdydIdy_MOY,w*dUdy*dIdy);
1683 AJOUT(WdVdydIdy_MOY,w*dVdy*dIdy);
1684 AJOUT(WdWdydIdy_MOY,w*dWdy*dIdy);
1685 AJOUT(WdUdzdIdz_MOY,w*dUdz*dIdz);
1686 AJOUT(WdVdzdIdz_MOY,w*dVdz*dIdz);
1687 AJOUT(WdWdzdIdz_MOY,w*dWdz*dIdz);
1688 // 4. de la vitesse et de la vitesse (fois l'indic) :
1689 //^^^^^^^^^^^^^^^^^^^^^^^
1690 AJOUT(IdUdxdUdx_MOY,chi*dUdx*dUdx);
1691 AJOUT(IdUdxdUdy_MOY,chi*dUdx*dUdy);
1692 AJOUT(IdUdxdUdz_MOY,chi*dUdx*dUdz);
1693 AJOUT(IdUdxdVdx_MOY,chi*dUdx*dVdx);
1694 AJOUT(IdUdxdVdy_MOY,chi*dUdx*dVdy);
1695 AJOUT(IdUdxdVdz_MOY,chi*dUdx*dVdz);
1696 AJOUT(IdUdxdWdx_MOY,chi*dUdx*dWdx);
1697 AJOUT(IdUdxdWdy_MOY,chi*dUdx*dWdy);
1698 AJOUT(IdUdxdWdz_MOY,chi*dUdx*dWdz);
1699 AJOUT(IdUdydUdy_MOY,chi*dUdy*dUdy);
1700 AJOUT(IdUdydUdz_MOY,chi*dUdy*dUdz);
1701 AJOUT(IdUdydVdx_MOY,chi*dUdy*dVdx);
1702 AJOUT(IdUdydVdy_MOY,chi*dUdy*dVdy);
1703 AJOUT(IdUdydVdz_MOY,chi*dUdy*dVdz);
1704 AJOUT(IdUdydWdx_MOY,chi*dUdy*dWdx);
1705 AJOUT(IdUdydWdy_MOY,chi*dUdy*dWdy);
1706 AJOUT(IdUdydWdz_MOY,chi*dUdy*dWdz);
1707 AJOUT(IdUdzdUdz_MOY,chi*dUdz*dUdz);
1708 AJOUT(IdUdzdVdx_MOY,chi*dUdz*dVdx);
1709 AJOUT(IdUdzdVdy_MOY,chi*dUdz*dVdy);
1710 AJOUT(IdUdzdVdz_MOY,chi*dUdz*dVdz);
1711 AJOUT(IdUdzdWdx_MOY,chi*dUdz*dWdx);
1712 AJOUT(IdUdzdWdy_MOY,chi*dUdz*dWdy);
1713 AJOUT(IdUdzdWdz_MOY,chi*dUdz*dWdz);
1714 AJOUT(IdVdxdVdx_MOY,chi*dVdx*dVdx);
1715 AJOUT(IdVdxdVdy_MOY,chi*dVdx*dVdy);
1716 AJOUT(IdVdxdVdz_MOY,chi*dVdx*dVdz);
1717 AJOUT(IdVdxdWdx_MOY,chi*dVdx*dWdx);
1718 AJOUT(IdVdxdWdy_MOY,chi*dVdx*dWdy);
1719 AJOUT(IdVdxdWdz_MOY,chi*dVdx*dWdz);
1720 AJOUT(IdVdydVdy_MOY,chi*dVdy*dVdy);
1721 AJOUT(IdVdydVdz_MOY,chi*dVdy*dVdz);
1722 AJOUT(IdVdydWdx_MOY,chi*dVdy*dWdx);
1723 AJOUT(IdVdydWdy_MOY,chi*dVdy*dWdy);
1724 AJOUT(IdVdydWdz_MOY,chi*dVdy*dWdz);
1725 AJOUT(IdVdzdVdz_MOY,chi*dVdz*dVdz);
1726 AJOUT(IdVdzdWdx_MOY,chi*dVdz*dWdx);
1727 AJOUT(IdVdzdWdy_MOY,chi*dVdz*dWdy);
1728 AJOUT(IdVdzdWdz_MOY,chi*dVdz*dWdz);
1729 AJOUT(IdWdxdWdx_MOY,chi*dWdx*dWdx);
1730 AJOUT(IdWdxdWdy_MOY,chi*dWdx*dWdy);
1731 AJOUT(IdWdxdWdz_MOY,chi*dWdx*dWdz);
1732 AJOUT(IdWdydWdy_MOY,chi*dWdy*dWdy);
1733 AJOUT(IdWdydWdz_MOY,chi*dWdy*dWdz);
1734 AJOUT(IdWdzdWdz_MOY,chi*dWdz*dWdz);
1735 //^^^^^^^^^^^^^^^^^^^^^^^
1736 // Other correlations to avoid using gradI :
1737 AJOUT(IdPdx_MOY, chi*dPdx);
1738 AJOUT(IdPdy_MOY, chi*dPdy);
1739 AJOUT(IdPdz_MOY, chi*dPdz);
1740 AJOUT(IUdPdx_MOY, chi*UdPdx);
1741 AJOUT(IUdPdy_MOY, chi*u*dPdy);
1742 AJOUT(IUdPdz_MOY, chi*u*dPdz);
1743 AJOUT(IVdPdx_MOY, chi*v*dPdx);
1744 AJOUT(IVdPdy_MOY, chi*VdPdy);
1745 AJOUT(IVdPdz_MOY, chi*v*dPdz);
1746 AJOUT(IWdPdx_MOY, chi*w*dPdx);
1747 AJOUT(IWdPdy_MOY, chi*w*dPdy);
1748 AJOUT(IWdPdz_MOY, chi*WdPdz);
1749 //
1750 AJOUT(IddUdxx_MOY, chi*ddUdxdx);
1751 AJOUT(IddVdxx_MOY, chi*ddVdxdx);
1752 AJOUT(IddWdxx_MOY, chi*ddWdxdx);
1753 AJOUT(IddUdyy_MOY, chi*ddUdydy);
1754 AJOUT(IddVdyy_MOY, chi*ddVdydy);
1755 AJOUT(IddWdyy_MOY, chi*ddWdydy);
1756 AJOUT(IddUdzz_MOY, chi*ddUdzdz);
1757 AJOUT(IddVdzz_MOY, chi*ddVdzdz);
1758 AJOUT(IddWdzz_MOY, chi*ddWdzdz);
1759 //
1760 AJOUT(IUdUdx_MOY, chi*u*dUdx);
1761 AJOUT(IUdVdx_MOY, chi*u*dVdx);
1762 AJOUT(IUdWdx_MOY, chi*u*dWdx);
1763 AJOUT(IUdUdy_MOY, chi*u*dUdy);
1764 AJOUT(IUdVdy_MOY, chi*u*dVdy);
1765 AJOUT(IUdWdy_MOY, chi*u*dWdy);
1766 AJOUT(IUdUdz_MOY, chi*u*dUdz);
1767 AJOUT(IUdVdz_MOY, chi*u*dVdz);
1768 AJOUT(IUdWdz_MOY, chi*u*dWdz);
1769 //
1770 AJOUT(IVdUdx_MOY, chi*v*dUdx);
1771 AJOUT(IVdVdx_MOY, chi*v*dVdx);
1772 AJOUT(IVdWdx_MOY, chi*v*dWdx);
1773 AJOUT(IVdUdy_MOY, chi*v*dUdy);
1774 AJOUT(IVdVdy_MOY, chi*v*dVdy);
1775 AJOUT(IVdWdy_MOY, chi*v*dWdy);
1776 AJOUT(IVdUdz_MOY, chi*v*dUdz);
1777 AJOUT(IVdVdz_MOY, chi*v*dVdz);
1778 AJOUT(IVdWdz_MOY, chi*v*dWdz);
1779 //
1780 AJOUT(IWdUdx_MOY, chi*w*dUdx);
1781 AJOUT(IWdVdx_MOY, chi*w*dVdx);
1782 AJOUT(IWdWdx_MOY, chi*w*dWdx);
1783 AJOUT(IWdUdy_MOY, chi*w*dUdy);
1784 AJOUT(IWdVdy_MOY, chi*w*dVdy);
1785 AJOUT(IWdWdy_MOY, chi*w*dWdy);
1786 AJOUT(IWdUdz_MOY, chi*w*dUdz);
1787 AJOUT(IWdVdz_MOY, chi*w*dVdz);
1788 AJOUT(IWdWdz_MOY, chi*w*dWdz);
1789 //
1790 //^^^^^^^^^^^^^^^^^^^^^^^
1791 // Other quantities :
1792 //
1793 // Force : pot * dIdxb
1794 AJOUT(Fx_MOY,Fx);
1795 AJOUT(Fy_MOY,Fy);
1796 AJOUT(Fz_MOY,Fz);
1797 //
1798 // Repulsion : pot_rep * dIdxb
1799 AJOUT(Frx_MOY,Frx);
1800 AJOUT(Fry_MOY,Fry);
1801 AJOUT(Frz_MOY,Frz);
1802 //
1803 // Repulsion : pot_rep * dIdxb
1804 AJOUT(Frax_MOY,Frax);
1805 AJOUT(Fray_MOY,Fray);
1806 AJOUT(Fraz_MOY,Fraz);
1807 //
1808 // La dissipation :
1809 AJOUT(DISSIP_MOY,chi*pseudo_dissip);
1810 AJOUT(DISSIP_VAP_MOY,chiv*pseudo_dissip);
1811 AJOUT(TRUE_DISSIP_MOY,chi*true_dissip);
1812 AJOUT(TRUE_DISSIP_VAP_MOY,chiv*true_dissip);
1813 AJOUT(POWER_INTERF,travail_Force_interfaces_vitesse);
1814 //
1815 // correlations 2 eme ordre vitesse (cote vapeur) :
1816 AJOUT(UUIv_MOY,u*u*chiv);
1817 AJOUT(VVIv_MOY,v*v*chiv);
1818 AJOUT(WWIv_MOY,w*w*chiv);
1819 AJOUT(UVIv_MOY,u*v*chiv);
1820 AJOUT(VWIv_MOY,v*w*chiv);
1821 AJOUT(UWIv_MOY,u*w*chiv);
1822 // b : bis (autre interpolation)
1823 AJOUT(UUIbv_MOY,uu*chiv);
1824 AJOUT(VVIbv_MOY,vv*chiv);
1825 AJOUT(WWIbv_MOY,ww*chiv);
1826 // correlations 2 eme ordre vitesse (liquide, autre interpolation) :
1827 AJOUT(UUIb_MOY,uu*chi);
1828 AJOUT(VVIb_MOY,vv*chi);
1829 AJOUT(WWIb_MOY,ww*chi);
1830 // Il manque une seule correl triple qui aurait la meme interp :
1831 // u*v*w car il faut interpoler chaque compo avant le produit...
1832 AJOUT(UUUIb_MOY,uuu*chi);
1833 AJOUT(UUVIb_MOY,uu*v*chi);
1834 AJOUT(UUWIb_MOY,uu*w*chi);
1835 AJOUT(UVVIb_MOY,u*vv*chi);
1836 AJOUT(UWWIb_MOY,u*ww*chi);
1837 AJOUT(VVVIb_MOY,vvv*chi);
1838 AJOUT(VVWIb_MOY,vv*w*chi);
1839 AJOUT(VWWIb_MOY,v*ww*chi);
1840 AJOUT(WWWIb_MOY,www*chi);
1841 //
1842 // ai (et nk?):
1843 AJOUT(AI_MOY, ai);
1844 AJOUT(Nx_MOY, Nx);
1845 AJOUT(Ny_MOY, Ny);
1846 AJOUT(Nz_MOY, Nz);
1847 //
1848 // AJOUT(NK_MOY, ???);
1849 //
1850 // Les stats manquantes pour l'equation d'epsilon :
1851 // Attention, elles ne sont pas toutes implementees!!!
1852 AJOUT(IdUdxdUdxdUdx_MOY,chi * dUdx * dUdx * dUdx);
1853 AJOUT(IdUdydUdydUdx_MOY,chi * dUdy * dUdy * dUdx);
1854 AJOUT(IdUdzdUdzdUdx_MOY,chi * dUdz * dUdz * dUdx);
1855 AJOUT(IdUdxdVdxdUdy_MOY,chi * dUdx * dVdx * dUdy);
1856 AJOUT(IdUdydVdydUdy_MOY,chi * dUdy * dVdy * dUdy);
1857 AJOUT(IdUdzdVdzdUdy_MOY,chi * dUdz * dVdz * dUdy);
1858 AJOUT(IdUdxdWdxdUdz_MOY,chi * dUdx * dWdx * dUdz);
1859 AJOUT(IdUdydWdydUdz_MOY,chi * dUdy * dWdy * dUdz);
1860 AJOUT(IdUdzdWdzdUdz_MOY,chi * dUdz * dWdz * dUdz);
1861 AJOUT(IdVdxdUdxdVdx_MOY,chi * dVdx * dUdx * dVdx);
1862 AJOUT(IdVdydUdydVdx_MOY,chi * dVdy * dUdy * dVdx);
1863 AJOUT(IdVdzdUdzdVdx_MOY,chi * dVdz * dUdz * dVdx);
1864 AJOUT(IdVdxdVdxdVdy_MOY,chi * dVdx * dVdx * dVdy);
1865 AJOUT(IdVdydVdydVdy_MOY,chi * dVdy * dVdy * dVdy);
1866 AJOUT(IdVdzdVdzdVdy_MOY,chi * dVdz * dVdz * dVdy);
1867 AJOUT(IdVdxdWdxdVdz_MOY,chi * dVdx * dWdx * dVdz);
1868 AJOUT(IdVdydWdydVdz_MOY,chi * dVdy * dWdy * dVdz);
1869 AJOUT(IdVdzdWdzdVdz_MOY,chi * dVdz * dWdz * dVdz);
1870 AJOUT(IdWdxdUdxdWdx_MOY,chi * dWdx * dUdx * dWdx);
1871 AJOUT(IdWdydUdydWdx_MOY,chi * dWdy * dUdy * dWdx);
1872 AJOUT(IdWdzdUdzdWdx_MOY,chi * dWdz * dUdz * dWdx);
1873 AJOUT(IdWdxdVdxdWdy_MOY,chi * dWdx * dVdx * dWdy);
1874 AJOUT(IdWdydVdydWdy_MOY,chi * dWdy * dVdy * dWdy);
1875 AJOUT(IdWdzdVdzdWdy_MOY,chi * dWdz * dVdz * dWdy);
1876 AJOUT(IdWdxdWdxdWdz_MOY,chi * dWdx * dWdx * dWdz);
1877 AJOUT(IdWdydWdydWdz_MOY,chi * dWdy * dWdy * dWdz);
1878 AJOUT(IdWdzdWdzdWdz_MOY,chi * dWdz * dWdz * dWdz);
1879 AJOUT(IdUdxdUdxW_MOY,chi * dUdx * dUdx * w);
1880 AJOUT(IdUdydUdyW_MOY,chi * dUdy * dUdy * w);
1881 AJOUT(IdUdzdUdzW_MOY,chi * dUdz * dUdz * w);
1882 AJOUT(IdVdxdVdxW_MOY,chi * dVdx * dVdx * w);
1883 AJOUT(IdVdydVdyW_MOY,chi * dVdy * dVdy * w);
1884 AJOUT(IdVdzdVdzW_MOY,chi * dVdz * dVdz * w);
1885 AJOUT(IdWdxdWdxW_MOY,chi * dWdx * dWdx * w);
1886 AJOUT(IdWdydWdyW_MOY,chi * dWdy * dWdy * w);
1887 AJOUT(IdWdzdWdzW_MOY,chi * dWdz * dWdz * w);
1888 AJOUT(IdUdxddPdxdx_MOY,chi * dUdx * ddPdxdx);
1889 AJOUT(IdUdyddPdxdy_MOY,chi * dUdy * ddPdxdy);
1890 AJOUT(IdUdzddPdxdz_MOY,chi * dUdz * ddPdxdz);
1891 AJOUT(IdVdxddPdydx_MOY,chi * dVdx * ddPdydx);
1892 AJOUT(IdVdyddPdydy_MOY,chi * dVdy * ddPdydy);
1893 AJOUT(IdVdzddPdydz_MOY,chi * dVdz * ddPdydz);
1894 AJOUT(IdWdxddPdzdx_MOY,chi * dWdx * ddPdzdx);
1895 AJOUT(IdWdyddPdzdy_MOY,chi * dWdy * ddPdzdy);
1896 AJOUT(IdWdzddPdzdz_MOY,chi * dWdz * ddPdzdz);
1897 AJOUT(IddUdxdxddUdxdx_MOY,chi * ddUdxdx * ddUdxdx);
1898 AJOUT(IddUdxdyddUdxdy_MOY,chi * ddUdxdy * ddUdxdy);
1899 AJOUT(IddUdxdzddUdxdz_MOY,chi * ddUdxdz * ddUdxdz);
1900 AJOUT(IddUdydxddUdydx_MOY,chi * ddUdydx * ddUdydx);
1901 AJOUT(IddUdydyddUdydy_MOY,chi * ddUdydy * ddUdydy);
1902 AJOUT(IddUdydzddUdydz_MOY,chi * ddUdydz * ddUdydz);
1903 AJOUT(IddUdzdxddUdzdx_MOY,chi * ddUdzdx * ddUdzdx);
1904 AJOUT(IddUdzdyddUdzdy_MOY,chi * ddUdzdy * ddUdzdy);
1905 AJOUT(IddUdzdzddUdzdz_MOY,chi * ddUdzdz * ddUdzdz);
1906 AJOUT(IddVdxdxddVdxdx_MOY,chi * ddVdxdx * ddVdxdx);
1907 AJOUT(IddVdxdyddVdxdy_MOY,chi * ddVdxdy * ddVdxdy);
1908 AJOUT(IddVdxdzddVdxdz_MOY,chi * ddVdxdz * ddVdxdz);
1909 AJOUT(IddVdydxddVdydx_MOY,chi * ddVdydx * ddVdydx);
1910 AJOUT(IddVdydyddVdydy_MOY,chi * ddVdydy * ddVdydy);
1911 AJOUT(IddVdydzddVdydz_MOY,chi * ddVdydz * ddVdydz);
1912 AJOUT(IddVdzdxddVdzdx_MOY,chi * ddVdzdx * ddVdzdx);
1913 AJOUT(IddVdzdyddVdzdy_MOY,chi * ddVdzdy * ddVdzdy);
1914 AJOUT(IddVdzdzddVdzdz_MOY,chi * ddVdzdz * ddVdzdz);
1915 AJOUT(IddWdxdxddWdxdx_MOY,chi * ddWdxdx * ddWdxdx);
1916 AJOUT(IddWdxdyddWdxdy_MOY,chi * ddWdxdy * ddWdxdy);
1917 AJOUT(IddWdxdzddWdxdz_MOY,chi * ddWdxdz * ddWdxdz);
1918 AJOUT(IddWdydxddWdydx_MOY,chi * ddWdydx * ddWdydx);
1919 AJOUT(IddWdydyddWdydy_MOY,chi * ddWdydy * ddWdydy);
1920 AJOUT(IddWdydzddWdydz_MOY,chi * ddWdydz * ddWdydz);
1921 AJOUT(IddWdzdxddWdzdx_MOY,chi * ddWdzdx * ddWdzdx);
1922 AJOUT(IddWdzdyddWdzdy_MOY,chi * ddWdzdy * ddWdzdy);
1923 AJOUT(IddWdzdzddWdzdz_MOY,chi * ddWdzdz * ddWdzdz);
1924 AJOUT(dIdxddUdxdxdUdx_MOY,dIdx * ddUdxdx * dUdx);
1925 AJOUT(dIdxddUdxdydUdy_MOY,dIdx * ddUdxdy * dUdy);
1926 AJOUT(dIdxddUdxdzdUdz_MOY,dIdx * ddUdxdz * dUdz);
1927 AJOUT(dIdyddUdydxdUdx_MOY,dIdy * ddUdydx * dUdx);
1928 AJOUT(dIdyddUdydydUdy_MOY,dIdy * ddUdydy * dUdy);
1929 AJOUT(dIdyddUdydzdUdz_MOY,dIdy * ddUdydz * dUdz);
1930 AJOUT(dIdzddUdzdxdUdx_MOY,dIdz * ddUdzdx * dUdx);
1931 AJOUT(dIdzddUdzdydUdy_MOY,dIdz * ddUdzdy * dUdy);
1932 AJOUT(dIdzddUdzdzdUdz_MOY,dIdz * ddUdzdz * dUdz);
1933 AJOUT(dIdxddVdxdxdVdx_MOY,dIdx * ddVdxdx * dVdx);
1934 AJOUT(dIdxddVdxdydVdy_MOY,dIdx * ddVdxdy * dVdy);
1935 AJOUT(dIdxddVdxdzdVdz_MOY,dIdx * ddVdxdz * dVdz);
1936 AJOUT(dIdyddVdydxdVdx_MOY,dIdy * ddVdydx * dVdx);
1937 AJOUT(dIdyddVdydydVdy_MOY,dIdy * ddVdydy * dVdy);
1938 AJOUT(dIdyddVdydzdVdz_MOY,dIdy * ddVdydz * dVdz);
1939 AJOUT(dIdzddVdzdxdVdx_MOY,dIdz * ddVdzdx * dVdx);
1940 AJOUT(dIdzddVdzdydVdy_MOY,dIdz * ddVdzdy * dVdy);
1941 AJOUT(dIdzddVdzdzdVdz_MOY,dIdz * ddVdzdz * dVdz);
1942 AJOUT(dIdxddWdxdxdWdx_MOY,dIdx * ddWdxdx * dWdx);
1943 AJOUT(dIdxddWdxdydWdy_MOY,dIdx * ddWdxdy * dWdy);
1944 AJOUT(dIdxddWdxdzdWdz_MOY,dIdx * ddWdxdz * dWdz);
1945 AJOUT(dIdyddWdydxdWdx_MOY,dIdy * ddWdydx * dWdx);
1946 AJOUT(dIdyddWdydydWdy_MOY,dIdy * ddWdydy * dWdy);
1947 AJOUT(dIdyddWdydzdWdz_MOY,dIdy * ddWdydz * dWdz);
1948 AJOUT(dIdzddWdzdxdWdx_MOY,dIdz * ddWdzdx * dWdx);
1949 AJOUT(dIdzddWdzdydWdy_MOY,dIdz * ddWdzdy * dWdy);
1950 AJOUT(dIdzddWdzdzdWdz_MOY,dIdz * ddWdzdz * dWdz);
1951 AJOUT(dIdxddUdxdz_MOY,dIdx * ddUdxdz);
1952 AJOUT(dIdyddUdydz_MOY,dIdy * ddUdydz);
1953 AJOUT(dIdzddUdzdz_MOY,dIdz * ddUdzdz);
1954 AJOUT(dIdzdUdxdUdx_MOY,dIdz * dUdx * dUdx);
1955 AJOUT(dIdzdUdydUdy_MOY,dIdz * dUdy * dUdy);
1956 AJOUT(dIdzdUdzdUdz_MOY,dIdz * dUdz * dUdz);
1957 AJOUT(dIdzdVdxdVdx_MOY,dIdz * dVdx * dVdx);
1958 AJOUT(dIdzdVdydVdy_MOY,dIdz * dVdy * dVdy);
1959 AJOUT(dIdzdVdzdVdz_MOY,dIdz * dVdz * dVdz);
1960 AJOUT(dIdzdWdxdWdx_MOY,dIdz * dWdx * dWdx);
1961 AJOUT(dIdzdWdydWdy_MOY,dIdz * dWdy * dWdy);
1962 AJOUT(dIdzdWdzdWdz_MOY,dIdz * dWdz * dWdz);
1963 // Ajout de tout ce qui avait du dIdx qui est mainetenant aiNx...
1964 AJOUT(UaiNx_MOY,u*aiNx);
1965 AJOUT(VaiNx_MOY,v*aiNx);
1966 AJOUT(WaiNx_MOY,w*aiNx);
1967 AJOUT(UaiNy_MOY,u*aiNy);
1968 AJOUT(VaiNy_MOY,v*aiNy);
1969 AJOUT(WaiNy_MOY,w*aiNy);
1970 AJOUT(UaiNz_MOY,u*aiNz);
1971 AJOUT(VaiNz_MOY,v*aiNz);
1972 AJOUT(WaiNz_MOY,w*aiNz);
1973 AJOUT(UPaiNx_MOY,u*p*aiNx);
1974 AJOUT(VPaiNx_MOY,v*p*aiNx);
1975 AJOUT(WPaiNx_MOY,w*p*aiNx);
1976 AJOUT(UPaiNy_MOY,u*p*aiNy);
1977 AJOUT(VPaiNy_MOY,v*p*aiNy);
1978 AJOUT(WPaiNy_MOY,w*p*aiNy);
1979 AJOUT(UPaiNz_MOY,u*p*aiNz);
1980 AJOUT(VPaiNz_MOY,v*p*aiNz);
1981 AJOUT(WPaiNz_MOY,w*p*aiNz);
1982 AJOUT(UUaiNx_MOY,uu*aiNx);
1983 AJOUT(VVaiNx_MOY,vv*aiNx);
1984 AJOUT(WWaiNx_MOY,ww*aiNx);
1985 AJOUT(UUaiNy_MOY,uu*aiNy);
1986 AJOUT(VVaiNy_MOY,vv*aiNy);
1987 AJOUT(WWaiNy_MOY,ww*aiNy);
1988 AJOUT(UUaiNz_MOY,uu*aiNz);
1989 AJOUT(VVaiNz_MOY,vv*aiNz);
1990 AJOUT(WWaiNz_MOY,ww*aiNz);
1991 AJOUT(UVaiNx_MOY,u*v*aiNx);
1992 AJOUT(UWaiNx_MOY,u*w*aiNx);
1993 AJOUT(VWaiNx_MOY,v*w*aiNx);
1994 AJOUT(UVaiNy_MOY,u*v*aiNy);
1995 AJOUT(UWaiNy_MOY,u*w*aiNy);
1996 AJOUT(VWaiNy_MOY,v*w*aiNy);
1997 AJOUT(UVaiNz_MOY,u*v*aiNz);
1998 AJOUT(UWaiNz_MOY,u*w*aiNz);
1999 AJOUT(VWaiNz_MOY,v*w*aiNz);
2000 AJOUT(PaiNx_MOY,p*aiNx);
2001 AJOUT(PaiNy_MOY,p*aiNy);
2002 AJOUT(PaiNz_MOY,p*aiNz);
2003 AJOUT(IaiNx_MOY,chi*aiNx);
2004 AJOUT(IaiNy_MOY,chi*aiNy);
2005 AJOUT(IaiNz_MOY,chi*aiNz);
2006 AJOUT(dUdxaiNx_MOY,dUdx*aiNx);
2007 AJOUT(dVdxaiNx_MOY,dVdx*aiNx);
2008 AJOUT(dWdxaiNx_MOY,dWdx*aiNx);
2009 AJOUT(dUdyaiNy_MOY,dUdy*aiNy);
2010 AJOUT(dVdyaiNy_MOY,dVdy*aiNy);
2011 AJOUT(dWdyaiNy_MOY,dWdy*aiNy);
2012 AJOUT(dUdzaiNz_MOY,dUdz*aiNz);
2013 AJOUT(dVdzaiNz_MOY,dVdz*aiNz);
2014 AJOUT(dWdzaiNz_MOY,dWdz*aiNz);
2015 AJOUT(UdUdxaiNx_MOY,u*dUdx*aiNx);
2016 AJOUT(UdVdxaiNx_MOY,u*dVdx*aiNx);
2017 AJOUT(UdWdxaiNx_MOY,u*dWdx*aiNx);
2018 AJOUT(UdUdyaiNy_MOY,u*dUdy*aiNy);
2019 AJOUT(UdVdyaiNy_MOY,u*dVdy*aiNy);
2020 AJOUT(UdWdyaiNy_MOY,u*dWdy*aiNy);
2021 AJOUT(UdUdzaiNz_MOY,u*dUdz*aiNz);
2022 AJOUT(UdVdzaiNz_MOY,u*dVdz*aiNz);
2023 AJOUT(UdWdzaiNz_MOY,u*dWdz*aiNz);
2024 AJOUT(VdUdxaiNx_MOY,v*dUdx*aiNx);
2025 AJOUT(VdVdxaiNx_MOY,v*dVdx*aiNx);
2026 AJOUT(VdWdxaiNx_MOY,v*dWdx*aiNx);
2027 AJOUT(VdUdyaiNy_MOY,v*dUdy*aiNy);
2028 AJOUT(VdVdyaiNy_MOY,v*dVdy*aiNy);
2029 AJOUT(VdWdyaiNy_MOY,v*dWdy*aiNy);
2030 AJOUT(VdUdzaiNz_MOY,v*dUdz*aiNz);
2031 AJOUT(VdVdzaiNz_MOY,v*dVdz*aiNz);
2032 AJOUT(VdWdzaiNz_MOY,v*dWdz*aiNz);
2033 AJOUT(WdUdxaiNx_MOY,w*dUdx*aiNx);
2034 AJOUT(WdVdxaiNx_MOY,w*dVdx*aiNx);
2035 AJOUT(WdWdxaiNx_MOY,w*dWdx*aiNx);
2036 AJOUT(WdUdyaiNy_MOY,w*dUdy*aiNy);
2037 AJOUT(WdVdyaiNy_MOY,w*dVdy*aiNy);
2038 AJOUT(WdWdyaiNy_MOY,w*dWdy*aiNy);
2039 AJOUT(WdUdzaiNz_MOY,w*dUdz*aiNz);
2040 AJOUT(WdVdzaiNz_MOY,w*dVdz*aiNz);
2041 AJOUT(WdWdzaiNz_MOY,w*dWdz*aiNz);
2042 AJOUT(aiNxddUdxdxdUdx_MOY,aiNx * ddUdxdx * dUdx);
2043 AJOUT(aiNxddUdxdydUdy_MOY,aiNx * ddUdxdy * dUdy);
2044 AJOUT(aiNxddUdxdzdUdz_MOY,aiNx * ddUdxdz * dUdz);
2045 AJOUT(aiNyddUdydydUdy_MOY,aiNy * ddUdydy * dUdy);
2046 AJOUT(aiNyddUdydzdUdz_MOY,aiNy * ddUdydz * dUdz);
2047 AJOUT(aiNzddUdzdzdUdz_MOY,aiNz * ddUdzdz * dUdz);
2048 AJOUT(aiNxddVdxdxdVdx_MOY,aiNx * ddVdxdx * dVdx);
2049 AJOUT(aiNxddVdxdydVdy_MOY,aiNx * ddVdxdy * dVdy);
2050 AJOUT(aiNxddVdxdzdVdz_MOY,aiNx * ddVdxdz * dVdz);
2051 AJOUT(aiNyddVdydydVdy_MOY,aiNy * ddVdydy * dVdy);
2052 AJOUT(aiNyddVdydzdVdz_MOY,aiNy * ddVdydz * dVdz);
2053 AJOUT(aiNzddVdzdzdVdz_MOY,aiNz * ddVdzdz * dVdz);
2054 AJOUT(aiNxddWdxdxdWdx_MOY,aiNx * ddWdxdx * dWdx);
2055 AJOUT(aiNxddWdxdydWdy_MOY,aiNx * ddWdxdy * dWdy);
2056 AJOUT(aiNxddWdxdzdWdz_MOY,aiNx * ddWdxdz * dWdz);
2057 AJOUT(aiNyddWdydydWdy_MOY,aiNy * ddWdydy * dWdy);
2058 AJOUT(aiNyddWdydzdWdz_MOY,aiNy * ddWdydz * dWdz);
2059 AJOUT(aiNzddWdzdzdWdz_MOY,aiNz * ddWdzdz * dWdz);
2060 AJOUT(aiNxddUdxdz_MOY,aiNx * ddUdxdz);
2061 AJOUT(aiNyddUdydz_MOY,aiNy * ddUdydz);
2062 AJOUT(aiNzddUdzdz_MOY,aiNz * ddUdzdz);
2063 AJOUT(aiNzdUdxdUdx_MOY,aiNz * dUdx * dUdx);
2064 AJOUT(aiNzdUdydUdy_MOY,aiNz * dUdy * dUdy);
2065 AJOUT(aiNzdUdzdUdz_MOY,aiNz * dUdz * dUdz);
2066 AJOUT(aiNzdVdxdVdx_MOY,aiNz * dVdx * dVdx);
2067 AJOUT(aiNzdVdydVdy_MOY,aiNz * dVdy * dVdy);
2068 AJOUT(aiNzdVdzdVdz_MOY,aiNz * dVdz * dVdz);
2069 AJOUT(aiNzdWdxdWdx_MOY,aiNz * dWdx * dWdx);
2070 AJOUT(aiNzdWdydWdy_MOY,aiNz * dWdy * dWdy);
2071 AJOUT(aiNzdWdzdWdz_MOY,aiNz * dWdz * dWdz);
2072
2073 AJOUT(aiNx_MOY,aiNx);
2074 AJOUT(aiNy_MOY,aiNy);
2075 AJOUT(aiNz_MOY,aiNz);
2076 AJOUT(kaiNx_MOY,kai*Nx);
2077 AJOUT(kaiNy_MOY,kai*Ny);
2078 AJOUT(kaiNz_MOY,kai*Nz);
2079
2080 AJOUT(aiNxx_MOY,ai*Nxx);
2081 AJOUT(aiNxy_MOY,ai*Nxy);
2082 AJOUT(aiNxz_MOY,ai*Nxz);
2083 AJOUT(aiNyy_MOY,ai*Nyy);
2084 AJOUT(aiNyz_MOY,ai*Nyz);
2085 AJOUT(aiNzz_MOY,ai*Nzz);
2086
2087 AJOUT(aiNNxxxx_MOY,ai*Nxx*Nxx);
2088 AJOUT(aiNNxxxy_MOY,ai*Nxx*Nxy);
2089 AJOUT(aiNNxxxz_MOY,ai*Nxx*Nxz);
2090 AJOUT(aiNNxxyy_MOY,ai*Nxx*Nyy);
2091 AJOUT(aiNNxxyz_MOY,ai*Nxx*Nyz);
2092 AJOUT(aiNNxxzz_MOY,ai*Nxx*Nzz);
2093
2094 AJOUT(aiNNxyyy_MOY,ai*Nxy*Nyy);
2095 AJOUT(aiNNxyyz_MOY,ai*Nxy*Nyz);
2096 AJOUT(aiNNxyzz_MOY,ai*Nxy*Nzz);
2097 AJOUT(aiNNxzyy_MOY,ai*Nxz*Nyy);
2098 AJOUT(aiNNxzyz_MOY,ai*Nxz*Nyz);
2099 AJOUT(aiNNxzzz_MOY,ai*Nxz*Nzz);
2100
2101 AJOUT(aiNNyyyy_MOY,ai*Nyy*Nyy);
2102 AJOUT(aiNNyyyz_MOY,ai*Nyy*Nyz);
2103 AJOUT(aiNNyyzz_MOY,ai*Nyy*Nzz);
2104 AJOUT(aiNNyzzz_MOY,ai*Nyz*Nzz);
2105 AJOUT(aiNNzzzz_MOY,ai*Nzz*Nzz);
2106
2107 AJOUT(kai_MOY,kai);
2108
2109
2110
2111
2112
2113 if(coef_immobilisation_>1e-16)
2114 {
2115 double intu = integral_vitesse_i(i,j,k);
2116 double intv = integral_vitesse_j(i,j,k);
2117 double intw = integral_vitesse_k(i,j,k);
2118 double intp = integral_pression(i,j,k);
2119 double rappelx=force_rappel[0](i,j,k);
2120 double rappely=force_rappel[1](i,j,k);
2121 double rappelz=force_rappel[2](i,j,k);
2122 double p_seuil =0., p_seuil_isup=0., p_seuil_iinf=0., p_seuil_jsup=0., p_seuil_jinf=0. , p_seuil_ksup=0. , p_seuil_kinf=0.;
2123 p_seuil = Calculer_valeur_seuil(p_seuil, pression(i,j,k), p_seuil_max, p_seuil_min);
2124 p_seuil_isup = Calculer_valeur_seuil(p_seuil_isup, pression(i+1,j,k), p_seuil_max, p_seuil_min);
2125 p_seuil_jsup = Calculer_valeur_seuil(p_seuil_jsup, pression(i,j+1,k), p_seuil_max, p_seuil_min);
2126 p_seuil_ksup = Calculer_valeur_seuil(p_seuil_ksup, pression(i,j,k+1), p_seuil_max, p_seuil_min);
2127 p_seuil_iinf = Calculer_valeur_seuil(p_seuil_iinf, pression(i-1,j,k), p_seuil_max, p_seuil_min);
2128 p_seuil_jinf = Calculer_valeur_seuil(p_seuil_jinf, pression(i,j-1,k), p_seuil_max, p_seuil_min);
2129 p_seuil_kinf = Calculer_valeur_seuil(p_seuil_kinf, pression(i,j,k-1), p_seuil_max, p_seuil_min);
2130 double dPdx_seuil;
2131 double dPdy_seuil;
2132 double dPdz_seuil;
2133 dPdx_seuil=((p_seuil_isup-p_seuil)*(1./dx_)+(p_seuil-p_seuil_iinf)*(1./dx_))*0.5;
2134 dPdy_seuil=((p_seuil_jsup-p_seuil)*(1./dy_)+(p_seuil-p_seuil_jinf)*(1./dy_))*0.5;
2135 dPdz_seuil=((p_seuil_ksup-p_seuil)*(1./dz)+(p_seuil-p_seuil_kinf)*(1./dz))*0.5;
2136 chi = floor(indicatrice(i,j,k));
2137 chiv = floor(1.-indicatrice(i,j,k));
2138 double chinp = 1.0-indicatrice_np(i,j,k);
2139
2140 double dUintdx=0., dVintdx=0., dWintdx=0., dUintdy=0., dVintdy=0., dWintdy=0., dUintdz=0., dVintdz=0., dWintdz=0.;
2141 double ddUintdxdx=0., ddVintdxdx=0., ddWintdxdx=0., ddUintdydy=0., ddVintdydy=0., ddWintdydy=0., ddUintdzdz=0., ddVintdzdz=0., ddWintdzdz=0.;
2142 double pseudo_dissipint = calculer_gradients_vitesse(integral_vitesse_i, integral_vitesse_j, integral_vitesse_k,
2143 i,j,k,
2144 dz,
2145 dUintdx, dVintdx, dWintdx,
2146 dUintdy, dVintdy, dWintdy,
2147 dUintdz, dVintdz, dWintdz,
2148 ddUintdxdx, ddVintdxdx, ddWintdxdx,
2149 ddUintdydy, ddVintdydy, ddWintdydy,
2150 ddUintdzdz, ddVintdzdz, ddWintdzdz,
2151 on_the_first_cell, on_the_last_cell);
2152// double true_dissipint = calculer_vraie_dissipation(pseudo_dissipint,
2153// dUintdx, dUintdy, dUintdz,
2154// dVintdx, dVintdy, dVintdz,
2155// dWintdx, dWintdy, dWintdz);
2156
2157 // ajout des stats pour bulle fixe et moyenne tmp
2158 AJOUT(UI_INT,intu*chi);
2159 AJOUT(VI_INT,intv*chi);
2160 AJOUT(WI_INT,intw*chi);
2161 AJOUT(UI_INTUI_INT, intu*intu*chi);
2162 AJOUT(UI_INTVI_INT, intu*intv*chi);
2163 AJOUT(UI_INTWI_INT, intu*intw*chi);
2164 AJOUT(VI_INTVI_INT, intv*intv*chi );
2165 AJOUT(VI_INTWI_INT, intv*intw*chi);
2166 AJOUT(WI_INTWI_INT, intw*intw*chi);
2167 AJOUT(UI_INTUI_INTUI_INT, intu*intu*intu*chi);
2168 AJOUT(UI_INTUI_INTVI_INT, intu*intu*intv*chi );
2169 AJOUT(UI_INTUI_INTWI_INT, intu*intu*intw*chi );
2170 AJOUT(UI_INTVI_INTUI_INT, intu*intv*intu*chi );
2171 AJOUT(UI_INTVI_INTVI_INT, intu*intv*intv*chi );
2172 AJOUT(UI_INTVI_INTWI_INT, intu*intv*intw*chi );
2173 AJOUT(UI_INTWI_INTUI_INT, intu*intw*intu*chi );
2174 AJOUT(UI_INTWI_INTVI_INT, intu*intw*intv*chi );
2175 AJOUT(UI_INTWI_INTWI_INT, intu*intw*intw*chi );
2176 AJOUT(VI_INTUI_INTUI_INT, intv*intu*intu*chi);
2177 AJOUT(VI_INTUI_INTVI_INT, intv*intu*intv*chi );
2178 AJOUT(VI_INTUI_INTWI_INT, intv*intu*intw*chi );
2179 AJOUT(VI_INTVI_INTUI_INT, intv*intv*intu*chi );
2180 AJOUT(VI_INTVI_INTVI_INT, intv*intv*intv*chi );
2181 AJOUT(VI_INTVI_INTWI_INT, intv*intv*intw*chi );
2182 AJOUT(VI_INTWI_INTUI_INT, intv*intw*intu*chi );
2183 AJOUT(VI_INTWI_INTVI_INT, intv*intw*intv*chi );
2184 AJOUT(VI_INTWI_INTWI_INT, intv*intw*intw*chi );
2185 AJOUT(WI_INTUI_INTUI_INT, intw*intu*intu*chi);
2186 AJOUT(WI_INTUI_INTVI_INT, intw*intu*intv*chi );
2187 AJOUT(WI_INTUI_INTWI_INT, intw*intu*intw*chi );
2188 AJOUT(WI_INTVI_INTUI_INT, intw*intv*intu*chi );
2189 AJOUT(WI_INTVI_INTVI_INT, intw*intv*intv*chi );
2190 AJOUT(WI_INTVI_INTWI_INT, intw*intv*intw*chi );
2191 AJOUT(WI_INTWI_INTUI_INT, intw*intw*intu*chi );
2192 AJOUT(WI_INTWI_INTVI_INT, intw*intw*intv*chi );
2193 AJOUT(WI_INTWI_INTWI_INT, intw*intw*intw*chi );
2194 AJOUT( UI_INTUIUI, intu*u*u*chi );
2195 AJOUT( UI_INTUIVI, intu*u*v*chi );
2196 AJOUT( UI_INTUIWI, intu*u*w*chi );
2197 AJOUT( UI_INTVIUI, intu*v*u*chi );
2198 AJOUT( UI_INTVIVI, intu*v*v*chi );
2199 AJOUT( UI_INTVIWI, intu*v*w*chi );
2200 AJOUT( UI_INTWIUI, intu*w*u*chi );
2201 AJOUT( UI_INTWIVI, intu*w*v*chi );
2202 AJOUT( UI_INTWIWI, intu*w*w*chi );
2203 AJOUT( VI_INTUIUI, intv*u*u*chi );
2204 AJOUT( VI_INTUIVI, intv*u*v*chi );
2205 AJOUT( VI_INTUIWI, intv*u*w*chi );
2206 AJOUT( VI_INTVIUI, intv*v*u*chi );
2207 AJOUT( VI_INTVIVI, intv*v*v*chi );
2208 AJOUT( VI_INTVIWI, intv*v*w*chi );
2209 AJOUT( VI_INTWIUI, intv*w*u*chi );
2210 AJOUT( VI_INTWIVI, intv*w*v*chi );
2211 AJOUT( VI_INTWIWI, intv*w*w*chi );
2212 AJOUT( WI_INTUIUI, intw*u*u*chi );
2213 AJOUT( WI_INTUIVI, intw*u*v*chi );
2214 AJOUT( WI_INTUIWI, intw*u*w*chi );
2215 AJOUT( WI_INTVIUI, intw*v*u*chi );
2216 AJOUT( WI_INTVIVI, intw*v*v*chi );
2217 AJOUT( WI_INTVIWI, intw*v*w*chi );
2218 AJOUT( WI_INTWIUI, intw*w*u*chi );
2219 AJOUT( WI_INTWIVI, intw*w*v*chi );
2220 AJOUT( WI_INTWIWI, intw*w*w*chi );
2221 AJOUT(UI_INTP_INT, intu*intp*chi );
2222 AJOUT(VI_INTP_INT, intv*intp*chi );
2223 AJOUT(WI_INTP_INT, intw*intp*chi );
2224
2225 AJOUT(dUINTdxdUINTdx, dUintdx*dUintdx*chi );
2226 AJOUT(dUINTdxdUINTdy, dUintdx*dUintdy*chi );
2227 AJOUT(dUINTdxdUINTdz, dUintdx*dUintdz*chi );
2228 AJOUT(dUINTdydUINTdy, dUintdy*dUintdy*chi );
2229 AJOUT(dUINTdydUINTdz, dUintdy*dUintdz*chi );
2230 AJOUT(dUINTdzdUINTdz, dUintdz*dUintdz*chi );
2231 AJOUT(dUINTdxdVINTdx, dUintdx*dVintdx*chi );
2232 AJOUT(dUINTdxdVINTdy, dUintdx*dVintdy*chi );
2233 AJOUT(dUINTdxdVINTdz, dUintdx*dVintdz*chi );
2234 AJOUT(dUINTdydVINTdy, dUintdy*dVintdy*chi );
2235 AJOUT(dUINTdydVINTdz, dUintdy*dVintdz*chi );
2236 AJOUT(dUINTdzdVINTdz, dUintdz*dVintdz*chi );
2237 AJOUT(dUINTdxdWINTdx, dUintdx*dWintdx*chi );
2238 AJOUT(dUINTdxdWINTdy, dUintdx*dWintdy*chi );
2239 AJOUT(dUINTdxdWINTdz, dUintdx*dWintdz*chi );
2240 AJOUT(dUINTdydWINTdy, dUintdy*dWintdy*chi );
2241 AJOUT(dUINTdydWINTdz, dUintdy*dWintdz*chi );
2242 AJOUT(dUINTdzdWINTdz, dUintdz*dWintdz*chi );
2243
2244
2245 AJOUT(dVINTdxdUINTdx, dVintdx*dUintdx*chi );
2246 AJOUT(dVINTdxdUINTdy, dVintdx*dUintdy*chi );
2247 AJOUT(dVINTdxdUINTdz, dVintdx*dUintdz*chi );
2248 AJOUT(dVINTdydUINTdy, dVintdy*dUintdy*chi );
2249 AJOUT(dVINTdydUINTdz, dVintdy*dUintdz*chi );
2250 AJOUT(dVINTdzdUINTdz, dVintdz*dUintdz*chi );
2251 AJOUT(dVINTdxdVINTdx, dVintdx*dVintdx*chi );
2252 AJOUT(dVINTdxdVINTdy, dVintdx*dVintdy*chi );
2253 AJOUT(dVINTdxdVINTdz, dVintdx*dVintdz*chi );
2254 AJOUT(dVINTdydVINTdy, dVintdy*dVintdy*chi );
2255 AJOUT(dVINTdydVINTdz, dVintdy*dVintdz*chi );
2256 AJOUT(dVINTdzdVINTdz, dVintdz*dVintdz*chi );
2257 AJOUT(dVINTdxdWINTdx, dVintdx*dWintdx*chi );
2258 AJOUT(dVINTdxdWINTdy, dVintdx*dWintdy*chi );
2259 AJOUT(dVINTdxdWINTdz, dVintdx*dWintdz*chi );
2260 AJOUT(dVINTdydWINTdy, dVintdy*dWintdy*chi );
2261 AJOUT(dVINTdydWINTdz, dVintdy*dWintdz*chi );
2262 AJOUT(dVINTdzdWINTdz, dVintdz*dWintdz*chi );
2263
2264 AJOUT(dWINTdxdUINTdx, dWintdx*dUintdx*chi );
2265 AJOUT(dWINTdxdUINTdy, dWintdx*dUintdy*chi );
2266 AJOUT(dWINTdxdUINTdz, dWintdx*dUintdz*chi );
2267 AJOUT(dWINTdydUINTdy, dWintdy*dUintdy*chi );
2268 AJOUT(dWINTdydUINTdz, dWintdy*dUintdz*chi );
2269 AJOUT(dWINTdzdUINTdz, dWintdz*dUintdz*chi );
2270 AJOUT(dWINTdxdVINTdx, dWintdx*dVintdx*chi );
2271 AJOUT(dWINTdxdVINTdy, dWintdx*dVintdy*chi );
2272 AJOUT(dWINTdxdVINTdz, dWintdx*dVintdz*chi );
2273 AJOUT(dWINTdydVINTdy, dWintdy*dVintdy*chi );
2274 AJOUT(dWINTdydVINTdz, dWintdy*dVintdz*chi );
2275 AJOUT(dWINTdzdVINTdz, dWintdz*dVintdz*chi );
2276 AJOUT(dWINTdxdWINTdx, dWintdx*dWintdx*chi );
2277 AJOUT(dWINTdxdWINTdy, dWintdx*dWintdy*chi );
2278 AJOUT(dWINTdxdWINTdz, dWintdx*dWintdz*chi );
2279 AJOUT(dWINTdydWINTdy, dWintdy*dWintdy*chi );
2280 AJOUT(dWINTdydWINTdz, dWintdy*dWintdz*chi );
2281 AJOUT(dWINTdzdWINTdz, dWintdz*dWintdz*chi );
2282 AJOUT(P_INTdUINTdx, intp*dUintdx*chi);
2283 AJOUT(P_INTdUINTdy, intp*dUintdy*chi);
2284 AJOUT(P_INTdUINTdz, intp*dUintdz*chi);
2285 AJOUT(P_INTdVINTdx, intp*dVintdx*chi);
2286 AJOUT(P_INTdVINTdy, intp*dVintdy*chi);
2287 AJOUT(P_INTdVINTdz, intp*dVintdz*chi);
2288 AJOUT(P_INTdWINTdx, intp*dWintdx*chi);
2289 AJOUT(P_INTdWINTdy, intp*dWintdy*chi);
2290 AJOUT(P_INTdWINTdz, intp*dWintdz*chi);
2291 AJOUT(DISSIP_INT,chi*pseudo_dissipint);
2292// AJOUT(DISSIP_VAP_INT,chiv*pseudo_dissipint);
2293// AJOUT(TRUE_DISSIP_INT,chi*true_dissipint);
2294// AJOUT(TRUE_DISSIP_VAP_INT,chiv*true_dissipint);
2295 AJOUT(U_NOPERTURBE,u*chinp);
2296 AJOUT(V_NOPERTURBE,v*chinp);
2297 AJOUT(W_NOPERTURBE,w*chinp);
2298 AJOUT(P_NOPERTURBE,p*chinp);
2299 AJOUT(I_NP,chinp);
2300 AJOUT(P_NP, p*chinp);
2301 AJOUT(PIU_NP, u*p*chinp);
2302 AJOUT(PIV_NP, v*p*chinp);
2303 AJOUT(PIW_NP, w*p*chinp);
2304 AJOUT(IPdUdx_NP,chinp*dUdx*p );
2305 AJOUT(IPdVdx_NP,chinp*dVdx*p );
2306 AJOUT(IPdWdx_NP,chinp*dWdx*p );
2307 AJOUT(IPdUdy_NP,chinp*dUdy*p );
2308 AJOUT(IPdVdy_NP,chinp*dVdy*p );
2309 AJOUT(IPdWdy_NP,chinp*dWdy*p );
2310 AJOUT(IPdUdz_NP,chinp*dUdz*p );
2311 AJOUT(IPdVdz_NP,chinp*dVdz*p );
2312 AJOUT(IPdWdz_NP,chinp*dWdz*p );
2313 AJOUT(IdPdx_NP, chinp*dPdx );
2314 AJOUT(IdPdy_NP, chinp*dPdy);
2315 AJOUT(IdPdz_NP, chinp*dPdz);
2316 AJOUT(PI_seuil, p_seuil*chi);
2317 AJOUT(PIU_seuil,u*p_seuil*chi );
2318 AJOUT(PIV_seuil,v*p_seuil*chi );
2319 AJOUT(PIW_seuil,w*p_seuil*chi );
2320 AJOUT(IPdUdx_seuil, chi*dUdx*p_seuil);
2321 AJOUT(IPdVdx_seuil, chi*dVdx*p_seuil);
2322 AJOUT(IPdWdx_seuil, chi*dWdx*p_seuil);
2323 AJOUT(IPdUdy_seuil, chi*dUdy*p_seuil);
2324 AJOUT(IPdVdy_seuil, chi*dVdy*p_seuil);
2325 AJOUT(IPdWdy_seuil, chi*dWdy*p_seuil);
2326 AJOUT(IPdUdz_seuil, chi*dUdz*p_seuil);
2327 AJOUT(IPdVdz_seuil, chi*dVdz*p_seuil);
2328 AJOUT(IPdWdz_seuil, chi*dWdz*p_seuil);
2329 AJOUT(IdPdx_seuil, dPdx_seuil*chi );
2330 AJOUT(IdPdy_seuil, dPdy_seuil*chi );
2331 AJOUT(IdPdz_seuil, dPdz_seuil*chi );
2332 AJOUT(P_NP_seuil, p_seuil*chinp);
2333 AJOUT(PIU_NP_seuil,u*p_seuil*chinp );
2334 AJOUT(PIV_NP_seuil,v*p_seuil*chinp );
2335 AJOUT(PIW_NP_seuil,w*p_seuil*chinp );
2336 AJOUT(IPdUdx_NP_seuil, chinp*dUdx*p_seuil);
2337 AJOUT(IPdVdx_NP_seuil, chinp*dVdx*p_seuil);
2338 AJOUT(IPdWdx_NP_seuil, chinp*dWdx*p_seuil);
2339 AJOUT(IPdUdy_NP_seuil, chinp*dUdy*p_seuil);
2340 AJOUT(IPdVdy_NP_seuil, chinp*dVdy*p_seuil);
2341 AJOUT(IPdWdy_NP_seuil, chinp*dWdy*p_seuil);
2342 AJOUT(IPdUdz_NP_seuil, chinp*dUdz*p_seuil);
2343 AJOUT(IPdVdz_NP_seuil, chinp*dVdz*p_seuil);
2344 AJOUT(IPdWdz_NP_seuil, chinp*dWdz*p_seuil);
2345 AJOUT(IdPdx_NP_seuil, dPdx_seuil*chinp );
2346 AJOUT(IdPdy_NP_seuil, dPdy_seuil*chinp );
2347 AJOUT(IdPdz_NP_seuil, dPdz_seuil*chinp );
2348 AJOUT(Ivrappelx, chiv*rappelx);
2349 AJOUT(Ivrappely, chiv*rappely);
2350 AJOUT(Ivrappelz, chiv*rappelz);
2351 AJOUT(IP_INT, chi*intp );
2352 AJOUT(UIvrappelx, chiv*rappelx*u);
2353 AJOUT(VIvrappelx, chiv*rappelx*v);
2354 AJOUT(WIvrappelx, chiv*rappelx*w);
2355 AJOUT(UIvrappely, chiv*rappely*u);
2356 AJOUT(VIvrappely, chiv*rappely*v);
2357 AJOUT(WIvrappely, chiv*rappely*w);
2358 AJOUT(UIvrappelz, chiv*rappelz*u);
2359 AJOUT(VIvrappelz, chiv*rappelz*v);
2360 AJOUT(WIvrappelz, chiv*rappelz*w);
2361 }
2362
2363 //NEW TERMS FOR THE DEVIATORIC PART OF THE VISCOUS STRESS TENSOR
2364 AJOUT(dUdyaiNx_MOY, dUdy*aiNx);
2365 AJOUT(dUdzaiNx_MOY, dUdz*aiNx);
2366 AJOUT(dVdxaiNy_MOY, dVdx*aiNy);
2367 AJOUT(dVdzaiNy_MOY, dVdz*aiNy);
2368 AJOUT(dWdxaiNz_MOY, dWdx*aiNz);
2369 AJOUT(dWdyaiNz_MOY, dWdy*aiNz);
2370
2371 // Iv*dUidxb : (9)
2372 AJOUT(dUdxIv_MOY,chiv*dUdx);
2373 AJOUT(dUdyIv_MOY,chiv*dUdy);
2374 AJOUT(dUdzIv_MOY,chiv*dUdz);
2375 AJOUT(dVdxIv_MOY,chiv*dVdx);
2376 AJOUT(dVdyIv_MOY,chiv*dVdy);
2377 AJOUT(dVdzIv_MOY,chiv*dVdz);
2378 AJOUT(dWdxIv_MOY,chiv*dWdx);
2379 AJOUT(dWdyIv_MOY,chiv*dWdy);
2380 AJOUT(dWdzIv_MOY,chiv*dWdz);
2381 // Attempt for viscous diffusion
2382 AJOUT(IddUdxdxU_MOY, chi*ddUdxdx*u);
2383 AJOUT(IddVdxdxU_MOY, chi*ddVdxdx*u);
2384 AJOUT(IddWdxdxU_MOY, chi*ddWdxdx*u);
2385 AJOUT(IddUdxdxV_MOY, chi*ddUdxdx*v);
2386 AJOUT(IddVdxdxV_MOY, chi*ddVdxdx*v);
2387 AJOUT(IddWdxdxV_MOY, chi*ddWdxdx*v);
2388 AJOUT(IddUdxdxW_MOY, chi*ddUdxdx*w);
2389 AJOUT(IddVdxdxW_MOY, chi*ddVdxdx*w);
2390 AJOUT(IddWdxdxW_MOY, chi*ddWdxdx*w);
2391
2392 AJOUT(IddUdydyU_MOY, chi*ddUdydy*u);
2393 AJOUT(IddVdydyU_MOY, chi*ddVdydy*u);
2394 AJOUT(IddWdydyU_MOY, chi*ddWdydy*u);
2395 AJOUT(IddUdydyV_MOY, chi*ddUdydy*v);
2396 AJOUT(IddVdydyV_MOY, chi*ddVdydy*v);
2397 AJOUT(IddWdydyV_MOY, chi*ddWdydy*v);
2398 AJOUT(IddUdydyW_MOY, chi*ddUdydy*w);
2399 AJOUT(IddVdydyW_MOY, chi*ddVdydy*w);
2400 AJOUT(IddWdydyW_MOY, chi*ddWdydy*w);
2401
2402 AJOUT(IddUdzdzU_MOY, chi*ddUdzdz*u);
2403 AJOUT(IddVdzdzU_MOY, chi*ddVdzdz*u);
2404 AJOUT(IddWdzdzU_MOY, chi*ddWdzdz*u);
2405 AJOUT(IddUdzdzV_MOY, chi*ddUdzdz*v);
2406 AJOUT(IddVdzdzV_MOY, chi*ddVdzdz*v);
2407 AJOUT(IddWdzdzV_MOY, chi*ddWdzdz*v);
2408 AJOUT(IddUdzdzW_MOY, chi*ddUdzdz*w);
2409 AJOUT(IddVdzdzW_MOY, chi*ddVdzdz*w);
2410 AJOUT(IddWdzdzW_MOY, chi*ddWdzdz*w);
2411
2412 AJOUT(Ix_MOY , x*chi);
2413 AJOUT(xaiNx_MOY, x*aiNx);
2414 AJOUT(xaiNy_MOY, x*aiNy);
2415 AJOUT(xaiNz_MOY, x*aiNz);
2416
2417 AJOUT(dUdxIx_MOY, x*chi*dUdx);
2418 AJOUT(dUdyIx_MOY, x*chi*dUdy);
2419 AJOUT(dUdzIx_MOY, x*chi*dUdz);
2420 AJOUT(dVdxIx_MOY, x*chi*dVdx);
2421 AJOUT(dVdyIx_MOY, x*chi*dVdy);
2422 AJOUT(dVdzIx_MOY, x*chi*dVdz);
2423 AJOUT(dWdxIx_MOY, x*chi*dWdx);
2424 AJOUT(dWdyIx_MOY, x*chi*dWdy);
2425 AJOUT(dWdzIx_MOY, x*chi*dWdz);
2426
2427 AJOUT(xUaiNx_MOY, x*u*aiNx);
2428 AJOUT(xVaiNx_MOY, x*v*aiNx);
2429 AJOUT(xWaiNx_MOY, x*w*aiNx);
2430 AJOUT(xUaiNy_MOY, x*u*aiNy);
2431 AJOUT(xVaiNy_MOY, x*v*aiNy);
2432 AJOUT(xWaiNy_MOY, x*w*aiNy);
2433 AJOUT(xUaiNz_MOY, x*u*aiNz);
2434 AJOUT(xVaiNz_MOY, x*v*aiNz);
2435 AJOUT(xWaiNz_MOY, x*w*aiNz);
2436
2437 //NEW TERMS FOR THE EXTENDED PRESSURE FIELDS ((Need to avoid the invalid cells in this esteem))
2438
2439 //if (chi != 0 && chi != 1){ Inserire tuti quelli definiti solo sull'interfaccia`
2440 //nel caso di pl vanno comunque escluse dal calcolo le celle invalide
2441 //}
2442 AJOUT(P_VAP_Iv_MOY,pv*chiv);
2443 //Momentum equation for the vapor phase
2444 AJOUT(P_VAP_aiNx_MOY,pv*aiNx); //definito solo sull'interfaccia
2445 AJOUT(P_VAP_aiNy_MOY,pv*aiNy); //definito solo sull'interfaccia
2446 AJOUT(P_VAP_aiNz_MOY,pv*aiNz); //definito solo sull'interfaccia
2447
2448 AJOUT(P_LIQ_I_MOY,pl*chi);
2449 //Momentum equation for the liquid phase
2450 AJOUT(P_LIQ_aiNx_MOY,pl*aiNx); //definito solo sull'interfaccia
2451 AJOUT(P_LIQ_aiNy_MOY,pl*aiNy); //definito solo sull'interfaccia
2452 AJOUT(P_LIQ_aiNz_MOY,pl*aiNz); //definito solo sull'interfaccia
2453 // Reynolds stress equation (only liquid phase)
2454 AJOUT(UP_LIQ_aiNx_MOY, u*pl*aiNx); //definito solo sull'interfaccia
2455 AJOUT(VP_LIQ_aiNx_MOY, v*pl*aiNx); //definito solo sull'interfaccia
2456 AJOUT(WP_LIQ_aiNx_MOY, w*pl*aiNx); //definito solo sull'interfaccia
2457 AJOUT(UP_LIQ_aiNy_MOY, u*pl*aiNy); //definito solo sull'interfaccia
2458 AJOUT(VP_LIQ_aiNy_MOY, v*pl*aiNy); //definito solo sull'interfaccia
2459 AJOUT(WP_LIQ_aiNy_MOY, w*pl*aiNy); //definito solo sull'interfaccia
2460 AJOUT(UP_LIQ_aiNz_MOY, u*pl*aiNz); //definito solo sull'interfaccia
2461 AJOUT(VP_LIQ_aiNz_MOY, v*pl*aiNz); //definito solo sull'interfaccia
2462 AJOUT(WP_LIQ_aiNz_MOY, w*pl*aiNz); //definito solo sull'interfaccia
2463 //Puo essere utile definirli su tutto il dominio, come correzione ai valori precedenti??
2464 AJOUT(IP_LIQ_dUdx_MOY, chi*pl*dUdx);
2465 AJOUT(IP_LIQ_dUdy_MOY, chi*pl*dUdy);
2466 AJOUT(IP_LIQ_dUdz_MOY, chi*pl*dUdz);
2467 AJOUT(IP_LIQ_dVdx_MOY, chi*pl*dVdx);
2468 AJOUT(IP_LIQ_dVdy_MOY, chi*pl*dVdy);
2469 AJOUT(IP_LIQ_dVdz_MOY, chi*pl*dVdz);
2470 AJOUT(IP_LIQ_dWdx_MOY, chi*pl*dWdx);
2471 AJOUT(IP_LIQ_dWdy_MOY, chi*pl*dWdy);
2472 AJOUT(IP_LIQ_dWdz_MOY, chi*pl*dWdz);
2473
2474 AJOUT(UP_LIQ_I_MOY, chi*u*pl);
2475 AJOUT(VP_LIQ_I_MOY, chi*v*pl);
2476 AJOUT(WP_LIQ_I_MOY, chi*w*pl);
2477
2478 //New stats to compare all pressure terms
2479 //AJOUT(IdP_LIQ_dx_MOY, chi*pl*);
2480
2481#undef AJOUT
2482 }
2483 }
2484 //rang_begin_liq_ = BEG_LIQ;
2485 for (int i = 0; i < nval_; i++)
2486 tmp(kglob, i) = moy[i] * facteur;
2487
2488 // For thermal variables :
2489 for (int i = 0; i < nvalt_; i++)
2490 for (int j = 0; j < nb_thermal_fields_; j++)
2491 tmpt(kglob, i,j) = moyv(i,j) * facteur;
2492 } // end loop k
2493
2494// Somme sur tous les processeurs:
2495 mp_sum_for_each_item(occ); // contains the number of invalid cells within the full slice
2496 // The second index is 0 for liq, 1 for vap.
2499
2500// Sur processeur 0, ajout de la contribution a l'integrale temporelle:
2502 {
2503 Nom str_occ = Nom(" ");
2504 Nom str_occ_vap = Nom(" ");
2505 for (int k = 0; k < nktot; k++)
2506 {
2507 double oc_liq = occ(k,0);
2508 double oc_vap = occ(k,1);
2509 str_occ += Nom(oc_liq)+Nom(" ");
2510 str_occ_vap += Nom(oc_vap)+Nom(" ");
2511 double facteur_liq = nijtot/(double)((nijtot)-oc_liq);
2512 double facteur_vap = nijtot/(double)((nijtot)-oc_vap);
2513 // Correct the value in tmp to account for invalid extended pressure cells :
2514 // BEG and END are included
2515 for (int i = BEG_LIQ; i <= END_LIQ; i++)
2516 tmp(k, i) *= facteur_liq;
2517
2518 for (int i = BEG_VAP; i <= END_VAP; i++)
2519 tmp(k, i) *= facteur_vap;
2520
2521 for (int i = 0; i < nval_; i++)
2522 {
2523 integrale_temporelle_[i][k] += tmp(k, i) * dt; //Numero di righe e di colonne e invertito riapetto a tmp
2524 moyenne_spatiale_instantanee_[i][k] = tmp(k, i);
2525 }
2526 for (int i = 0; i < nvalt_; i++)
2527 for (int j = 0; j < nb_thermal_fields_; j++)
2528 {
2529 integrale_temporelle_temperature_(i,k,j) += tmpt(k, i,j) * dt;
2530 moyenne_spatiale_instantanee_temperature_(i,k,j) = tmpt(k, i,j);
2531 }
2532 }
2533 Cout<< "Number of invalid cells for liquid: " << str_occ << finl;
2534 Cout<< "Number of invalid cells for vapor: " << str_occ_vap << finl;
2535 t_integration_ += dt;
2536 }
2537
2538#ifdef STAT_VERBOSE
2539// Pour verifications : impressions des normes :
2540 Journal() << "VERBOSE update_state: " << finl;
2541 Journal() << " erreur, erreur_ddP, divU, laplP, d_divU_dx, d_divU_dy, d_divU_dz" << finl;
2542 Journal() << " L1 "
2543 << L1_erreur << " "
2544 << L1_erreur_ddP << " "
2545 << L1_divU << " ";
2546 Journal() << L1_laplP << " "
2547 << L1_d_divU_dx << " "
2548 << L1_d_divU_dy << " "
2549 << L1_d_divU_dz << finl;
2550 Journal() << " L2 "
2551 << L2_erreur << " "
2552 << L2_erreur_ddP << " "
2553 << L2_divU << " ";
2554 Journal() << L2_laplP << " "
2555 << L2_d_divU_dx << " "
2556 << L2_d_divU_dy << " "
2557 << L2_d_divU_dz << finl;
2558#endif
2559 Cout << "Statistiques_dns_ijk_FT::update_stat" << finl;
2560}
2561
2562// This method is called only on master proc.
2563void Statistiques_dns_ijk_FT::postraiter_thermique(const double current_time) const
2564{
2566 return;
2567
2568 const int nz = elem_coord_.size_array(); // nombre de points en K
2569 for (int ith = 0; ith < nb_thermal_fields_; ith++)
2570 {
2571 Nom n("");
2572 for (int flag_valeur_instantanee=0; flag_valeur_instantanee<2; flag_valeur_instantanee++)
2573 {
2575 n="monophasique_";
2576 else
2577 n="diphasique_";
2578
2579 if (flag_valeur_instantanee == 0)
2580 n+="statistiques_thermique_";
2581 else
2582 n+="moyenne_spatiale_thermique_";
2583
2584 n += Nom(ith)+Nom("_")+Nom(current_time)+Nom(".txt");
2585 SFichier os(n);
2586 os.setf(ios::scientific);
2587 os.precision(15);
2588 if (flag_valeur_instantanee == 0)
2589 {
2590 os << "# temps_integration " << t_integration_ << finl;
2591 os << "# Impression des moyennes temporelles" << finl;
2592 }
2593 else
2594 os << "# Impression des moyennes spatiales instantanee" << finl;
2595
2596 os << "# colonne 1 : coordonnee_K" << finl;
2597 for (int i = 0; i < nvalt_; i++)
2598 os << "# colonne " << i+2 << " : " << noms_moyennes_temperature_[i] << finl;
2599
2600 for (int j = 0; j < nz; j++)
2601 {
2602 char s[100];
2603 snprintf(s, 100, "%16.16e ", elem_coord_[j]);
2604 os << s;
2605 for (int i = 0; i < nvalt_; i++)
2606 {
2607 double x;
2608 if (flag_valeur_instantanee == 0)
2610 else
2612 snprintf(s, 100, "%16.16e ", x);
2613 os << s;
2614 }
2615 os << finl;
2616 }
2617 }
2618 }
2619}
2620
2621void Statistiques_dns_ijk_FT::postraiter(Sortie& os, int flag_valeur_instantanee) const
2622{
2623 // Mother-class method
2624 Statistiques_dns_ijk::postraiter(os,flag_valeur_instantanee);
2625 // Completer par la thermique
2626 //if (nb_thermal_fields_)
2627 // posttraiter_thermique a besoin du temps mais sinon on aurrait pu l'appeler la.
2628}
2629
2630// Impression des donnees pour reprise
2632{
2634 return os;
2635}
2636
2638{
2640 {
2641 const int nz = elem_coord_.size_array(); // nombre de points en K
2642 ArrOfDouble tmp_integrale_temporelle(nz);
2643 for (int ith = 0; ith < nb_thermal_fields_; ith++)
2644 {
2645 os << " thermal_fields " << ith<< " { " << "\n";
2646 for (int i = 0; i < nvalt_; i++)
2647 {
2648 for (int k=0; k< nz; k++)
2649 tmp_integrale_temporelle[k] = integrale_temporelle_temperature_(i,k,ith);
2650
2651 os << noms_moyennes_temperature_[i] << " " << tmp_integrale_temporelle << finl;
2652 }
2653 os << " }" << finl;
2654 }
2655 Cout << "Statistical thermal fields written for restart: t_integration=" << t_integration_ << finl;
2656 }
2657 return os;
2658}
2659
2660// Reprise des donnees stat dans un fichier reprise
2661// Attention, cette methode peut etre appelee avant initialize() !
2663{
2665 return is;
2666}
2667
2669{
2670 Cerr << "[Statistiques_dns_ijk_FT] Motcle : " << mot << finl;
2671 Param param(que_suis_je()+Nom("_fields"));
2672 int ret_val = 1;
2673 if (mot=="thermal_fields")
2674 {
2675 VECT(ArrOfDouble) tmp;
2676 tmp.dimensionner(nval_);
2677 int ith = -1;
2678 is >> ith;
2679
2680 for (int i = 0; i < nvalt_; i++)
2681 {
2682 param.ajouter(noms_moyennes_temperature_[i], &tmp[i]); // que mettre.. Ne rien faire ??? il faudrait un VECT(VECT(A
2683 }
2684 param.lire_avec_accolades(is);
2685
2686 // Copy the tmp into the good cell of the target :
2687 const int nz = tmp[0].size_array();
2688 // At last, we know the sizes and can dimension the tables correctly :
2689 if (ith ==0)
2690 {
2691 // no need to repeat it several times, and I'm afraid the resize
2692 // might erase the first field if no option like NO_INIT is used...
2695 }
2696 for (int i = 0; i < nvalt_; i++)
2697 for (int k = 0; k< nz; k++)
2698 integrale_temporelle_temperature_(i,k,ith) = tmp[i][k] ; // (nvalt_,nb_elem_k_tot,nb_thermal_fields_)
2699
2700 Cerr << "End of reading 'thermal_fields' in Statistiques_dns_ijk_FT" << finl;
2701 }
2702 else
2703 {
2704 Cerr << "Keyword " << mot << " does not belong to Statistiques_dns_ijk_FT, "
2705 << "resorting to the mother-class for resolution" << finl;
2707 }
2708 return ret_val;
2709}
2710
2712{
2713 // At this stage, nz is 0, we haven't read the initialize.
2714 // Thus, we cannot resize integrale_temporelle_temperature_ or moyenne_spatiale_instantanee_temperature_ correctly.
2715
2716 param.ajouter_non_std("thermal_fields",(this));
2717 // ou pour lire quand meme les anciens directement sans accolade :
2718 // Ce format n'a normalement jamais existe pour la temperature.
2719 /* for (int i = 0; i < nvalt_; i++)
2720 {
2721 param.ajouter(noms_moyennes_temperature_[i], &integrale_temporelle_temperature_[i]); // que mettre.. Ne rien faire ??? il faudrait un VECT(VECT(A
2722 } */
2723 return;
2724}
2725
2726// Attention, cette methode est appelee apres readOn(),
2727// il ne faut pas casser les donnees lues
2729{
2730 ref_ijk_ft_=ijk_ft;
2731 dx_=geom.get_constant_delta(0); //modif AT 20/06/2013
2732 dy_=geom.get_constant_delta(1); //modif AT 20/06/2013
2733 tab_dz_=geom.get_delta(2); //modif AT 20/06/2013
2734 const int n = geom.get_nb_elem_tot(DIRECTION_K);
2735 elem_coord_.resize_array(n);
2736 const ArrOfDouble& coord_z = geom.get_node_coordinates(DIRECTION_K);
2737 for (int i = 0; i < n; i++)
2738 elem_coord_[i] = (coord_z[i] + coord_z[i+1]) * 0.5;
2739
2740 if (ref_ijk_ft_->has_thermals())
2741 {
2742 nb_thermal_fields_ = ref_ijk_ft_->get_ijk_thermals().size();
2743 Cerr << que_suis_je() << " initialised with nb_thermal_fields=" << nb_thermal_fields_ << finl;
2744 }
2746 {
2747 moyenne_spatiale_instantanee_.dimensionner(nval_);
2748 for (int i = 0; i < nval_; i++)
2749 {
2750 moyenne_spatiale_instantanee_[i].resize_array(n);
2751 }
2752 // moyenne_spatiale_instantanee_[WFACE_MOY].resize_array(n+1);//En effet, il y a une face de plus que d'elem et ce tableau est au face
2755 }
2756
2757 int flag = 0;
2758 if (integrale_temporelle_.size() == 0)
2759 {
2760 integrale_temporelle_.dimensionner(nval_);
2761 t_integration_ = 0.;
2762 }
2763 else
2764 flag = 1;
2765 for (int i = 0; i < nval_; i++)
2766 {
2767 if (flag)
2768 {
2769 if (integrale_temporelle_[i].size_array() != n)
2770 {
2771 Cerr << "Erreur dans Statistiques_dns_ijk_FT::initialize: reprise avec mauvais nombre de mailles en z" << finl;
2772 Cerr << "integrale_temporelle_[" << i << "].size_array() = " << integrale_temporelle_[i].size_array() << finl;
2773 Process::exit();
2774 }
2775 }
2776 else
2777 {
2778 integrale_temporelle_[i].resize_array(n);
2779
2780 }
2781 }
2782//modif AT 20/06/2013
2783 if (check_converge_)
2784 {
2785 //on recupere les vitesses moyennes aux faces pour calculer les fluctuations de vitesse.
2786 vit_moy_.dimensionner(3);
2787 for (int i = 0; i < 3; i++)
2788 vit_moy_[i].resize_array(n);
2789
2790 vit_moy_[0]=integrale_temporelle_[0];
2791 vit_moy_[1]=integrale_temporelle_[1];
2792 vit_moy_[2]=integrale_temporelle_[40];//car on veut w aux faces !
2793 }
2794 check_stats_ = 0;
2795}
2796
2798{
2799 // TODO (teo boutin) only allocate these fields when needed, but right now I don't have the time to determine when that is
2800 // Pour verification des stats :
2801 {
2802 const int ghost = 1; // 1 ghost cell necessary to permit face-to-cell interpolation in post-pro...
2803 // Face-vectors (std vectors)
2804 const std::vector<std::string> noms_vecto_faces = {"dPd"};
2805 for (const auto& nam : noms_vecto_faces)
2806 {
2807 {
2808 vect_post_fields_[nam] = IJK_Field_vector3_double();
2809 allocate_velocity(vect_post_fields_.at(nam), domaine_ijk_, ghost, nam);
2810 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(nam));
2811 }
2812 const Nom& fld_nam = Nom("ANA_")+Nom(nam);
2813 const std::string& ana_nam = fld_nam.getString();
2814 {
2815 vect_post_fields_[ana_nam] = IJK_Field_vector3_double();
2816 allocate_velocity(vect_post_fields_.at(ana_nam), domaine_ijk_, 0, ana_nam);
2817 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(ana_nam));
2818 }
2819 }
2820
2821 // Cell-vectors
2822 const std::vector<std::string> noms_vectoriels = {"dUd", "dVd", "dWd"};
2823 for (const auto& nam : noms_vectoriels)
2824 {
2825 {
2826 vect_post_fields_[nam] = IJK_Field_vector3_double();
2827 allocate_cell_vector(vect_post_fields_.at(nam), domaine_ijk_, ghost, nam);
2828 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(nam));
2829 }
2830 const Nom fld_nam = Nom("ANA_")+Nom(nam);
2831 const std::string ana_nam = fld_nam.getString();
2832 {
2833 vect_post_fields_[ana_nam] = IJK_Field_vector3_double();
2834 allocate_cell_vector(vect_post_fields_.at(ana_nam), domaine_ijk_, 0, ana_nam);
2835 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(ana_nam));
2836 }
2837 }
2838
2839 // TODO: Deal with tensors for crossed compo "ddPd", "ddUd"
2840 /* Temorary fix : allocate a second cell-vector for cross-components and register it. For example, compos are in this order:
2841 * ddPdxdy, ddPdxdz, ddPdydz
2842 * But their names are :
2843 * ddPddc_x, ddPddc_y, ddPddc_z -> c: cross-component
2844 */
2845 // TODO: rename components ?
2846 const std::vector<std::string> noms_tensoriels = {"ddPdd", "ddUdd", "ddVdd", "ddWdd"};
2847 for (const auto& nam : noms_tensoriels)
2848 {
2849 {
2850 vect_post_fields_[nam] = IJK_Field_vector3_double();
2851 allocate_cell_vector(vect_post_fields_.at(nam), domaine_ijk_, 0, nam);
2852 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(nam));
2853 }
2854 const std::string nam2 = Nom(nam+"c").getString(); // For off-diagonal components
2855 {
2856 vect_post_fields_[nam2] = IJK_Field_vector3_double();
2857 allocate_cell_vector(vect_post_fields_.at(nam2), domaine_ijk_, 0, nam2);
2858 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(nam2));
2859 }
2860 const Nom& fld_nam = Nom("ANA_")+Nom(nam);
2861 const std::string ana_nam = fld_nam.getString();
2862 {
2863 vect_post_fields_[ana_nam] = IJK_Field_vector3_double();
2864 allocate_cell_vector(vect_post_fields_.at(ana_nam), domaine_ijk_, 0, ana_nam);
2865 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(ana_nam));
2866 }
2867 const std::string ana_nam2 = Nom(fld_nam+"c").getString();
2868 {
2869 vect_post_fields_[ana_nam2] = IJK_Field_vector3_double();
2870 allocate_cell_vector(vect_post_fields_.at(ana_nam2), domaine_ijk_, 0, ana_nam2);
2871 champs_compris_.ajoute_champ_vectoriel(vect_post_fields_.at(ana_nam2));
2872 }
2873 }
2874 }
2875}
2876
2878 const int check_stats)
2879{
2880 initialize(ijk_ft,geom);
2881 // veut-on verifier les derivees :
2882 check_stats_ = check_stats;
2883 alloc_fields();
2884}
2885
2886/** Retrieve requested field for postprocessing, potentially updating it.
2887 */
2888const IJK_Field_double& Statistiques_dns_ijk_FT::get_IJK_field(const Motcle& nom)
2889{
2890 if (!has_champ(nom))
2891 {
2892 Cerr << "ERROR in Statistiques_dns_ijk::get_IJK_field : " << finl;
2893 Cerr << "Requested field '" << nom << "' is not recognized by Statistiques_dns_ijk::get_IJK_field()." << finl;
2894 throw;
2895 }
2896
2897 return champs_compris_.get_champ(nom);
2898}
2899
2900const IJK_Field_vector3_double& Statistiques_dns_ijk_FT::get_IJK_field_vector(const Motcle& nom)
2901{
2902 if (!has_champ_vectoriel(nom))
2903 {
2904 Cerr << "ERROR in Statistiques_dns_ijk::get_IJK_field_vector : " << finl;
2905 Cerr << "Requested field '" << nom << "' is not recognized by Statistiques_dns_ijk::get_IJK_field_vector()." << finl;
2906 throw;
2907 }
2908
2909 double current_time = ref_ijk_ft_->schema_temps_ijk().get_current_time();
2910 if (nom == "ANA_dPd")
2911 {
2912 for (int i = 0; i < 3; i++)
2913 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_gradP_analytique_[i], current_time);
2914 }
2915 if (nom == "ANA_dUd")
2916 {
2917 for (int i = 0; i < 3; i++)
2918 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_gradU_analytique_[i], current_time);
2919 }
2920 if (nom == "ANA_dVd")
2921 {
2922 for (int i = 0; i < 3; i++)
2923 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_gradV_analytique_[i], current_time);
2924 }
2925 if (nom == "ANA_dWd")
2926 {
2927 for (int i = 0; i < 3; i++)
2928 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_gradW_analytique_[i], current_time);
2929 }
2930 // Pour les deriv secondes :
2931 if (nom == "ANA_ddPdd")
2932 {
2933 for (int i = 0; i < 3; i++)
2934 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2P_analytique_[i], current_time);
2935 }
2936 if (nom == "ANA_ddUdd")
2937 {
2938 for (int i = 0; i < 3; i++)
2939 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2U_analytique_[i], current_time);
2940 }
2941 if (nom == "ANA_ddVdd")
2942 {
2943 for (int i = 0; i < 3; i++)
2944 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2V_analytique_[i], current_time);
2945 }
2946 if (nom == "ANA_ddWdd")
2947 {
2948 for (int i = 0; i < 3; i++)
2949 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2W_analytique_[i], current_time);
2950 }
2951 // Et pour les croise des deriv secondes : (Attention au +3 dans les compo du parser
2952 if (nom == "ANA_ddPddc")
2953 {
2954 for (int i = 0; i < 3; i++)
2955 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2P_analytique_[i+3], current_time);
2956 }
2957 if (nom == "ANA_ddUddc")
2958 {
2959 for (int i = 0; i < 3; i++)
2960 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2U_analytique_[i+3], current_time);
2961 }
2962 if (nom == "ANA_ddVddc")
2963 {
2964 for (int i = 0; i < 3; i++)
2965 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2V_analytique_[i+3], current_time);
2966 }
2967 if (nom == "ANA_ddWddc")
2968 {
2969 for (int i = 0; i < 3; i++)
2970 set_field_data(vect_post_fields_.at(nom)[i], ref_ijk_ft_->get_post().expression_grad2W_analytique_[i+3], current_time);
2971 }
2972 return champs_compris_.get_champ_vectoriel(nom);
2973}
2974
2976 const double portee_wall_repulsion) const
2977{
2979 {
2980 const double zmin = geom_NS.get_origin(DIRECTION_K);
2981 const double zmax = zmin+geom_NS.get_domain_length(DIRECTION_K);
2982 const int n = elem_coord_.size_array(); // nombre de points en K
2983 const int i =0 ; // LA CASE DU TAUX DE VIDE...
2984 double almoy_g = 0., almoy_d = 0.;
2985 int ng = 0, nd = 0;
2986 for (int j = 0; j < n; j++)
2987 {
2988 const double z = elem_coord_[j];
2989 if (z - zmin < portee_wall_repulsion)
2990 {
2991 // Dans le domaine de repulsion gauche
2992 almoy_g +=moyenne_spatiale_instantanee_[i][j];
2993 ng++;
2994 }
2995 else if (zmax - z < portee_wall_repulsion )
2996 {
2997 // Dans le domaine de repulsion droite
2998 almoy_d +=moyenne_spatiale_instantanee_[i][j];
2999 nd++;
3000 }
3001 }
3002 if (ng>0)
3003 almoy_g /= ng;
3004 if (nd>0)
3005 almoy_d /= nd;
3006 return (almoy_g - almoy_d);
3007
3008 }
3009 else
3010 {
3011 Cerr << " Statistiques_dns_ijk_ft::compute_desequil_alpha() should be called on master" << finl;
3012 Process::exit();
3013 return -1;
3014 }
3015}
This class encapsulates all the information related to the eulerian mesh for TrioIJK.
Definition Domaine_IJK.h:47
double get_domain_length(int direction) const
Returns the length of the entire domain in requested direction.
double get_constant_delta(int direction) const
Returns the size of cells in a direction.
const ArrOfDouble & get_delta(int direction) const
Returns the array of mesh cell sizes in requested direction.
double get_origin(int direction) const
Returns the coordinate of the first node (global) of the mesh in the requested direction.
int get_nb_elem_tot(int direction) const
Returns the total (global) number of mesh cells in requested direction.
const ArrOfDouble & get_node_coordinates(int direction) const
Returns an array with the coordinates of all nodes in the mesh in requested direction.
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
void allocate(int ni, int nj, int nk, int ghosts, int additional_k_layers=0, int nb_compo=1, bool external_storage=false)
const IJK_Field_double & I() const
Une chaine de caractere (Nom) en majuscules.
Definition Motcle.h:26
Un tableau d'objets de la classe Motcle.
Definition Motcle.h:63
int contient_(const char *const ch) const
Definition Motcle.cpp:333
IJK_Field_vector3_double terme_repulsion_interfaces_ft_
IJK_Field_vector3_double velocity_
IJK_Field_vector3_double terme_repulsion_interfaces_ns_
IJK_Field_vector3_double terme_source_interfaces_ns_
IJK_Field_vector3_double force_rappel_
IJK_Field_vector3_double terme_abs_repulsion_interfaces_ns_
IJK_Field_vector3_double terme_source_interfaces_ft_
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
const std::string & getString() const
Definition Nom.h:92
friend class Entree
Definition Objet_U.h:76
friend class Sortie
Definition Objet_U.h:75
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
static bool DISABLE_DIPHASIQUE
Definition Option_IJK.h:26
Helper class to factorize the readOn method of Objet_U classes.
Definition Param.h:112
void ajouter(const char *keyword, const int *value, Param::Nature nat=Param::OPTIONAL)
Register an integer parameter.
Definition Param.cpp:364
void ajouter_non_std(const char *keyword, const Objet_U *value, Param::Nature nat=Param::OPTIONAL)
Register a keyword handled by Objet_U::lire_motcle_non_standard.
Definition Param.cpp:489
int lire_avec_accolades(Entree &is)
Alias of lire_avec_accolades_depuis.
Definition Param.h:577
IJK_Field_double ai_ns_
IJK_Field_double kappa_ai_ns_
IJK_Field_double integrated_pressure_
IJK_Field_double extended_pv_
IJK_Field_double indicatrice_non_perturbe_
IJK_Field_vector3_double integrated_velocity_
IJK_Field_vector3_double grad_I_ns_
IJK_Field_vector3_double normale_cell_ns_
IJK_Field_double extended_pl_
const IJK_Interfaces & get_interface() const
const Postprocessing_IJK & get_post() const
const Navier_Stokes_FTD_IJK & eq_ns() const
static void mp_sum_for_each_item(TRUSTArray< _TYPE_ > &x, int n=-1)
Definition Process.cpp:193
static Sortie & Journal(int message_level=0)
Renvoie un objet statique de type Sortie qui sert de journal d'evenements.
Definition Process.cpp:588
static int me()
renvoie mon rang dans le groupe de communication courant.
Definition Process.cpp:125
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Definition Process.cpp:86
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
Definition SFichier.h:27
void precision(int pre) override
void setf(IOS_FORMAT code) override
Classe de base des flux de sortie.
Definition Sortie.h:52
Sortie & completer_print(Sortie &os) const override
const IJK_Field_double & get_IJK_field(const Motcle &nom) override
int lire_motcle_non_standard(const Motcle &mot, Entree &is) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
const IJK_Field_vector3_double & get_IJK_field_vector(const Motcle &nom) override
void postraiter(Sortie &, int flag_valeur_instantanee=0) const
void completer_read(Param &param) override
VECT(Nom) noms_moyennes_temperature_
void associer_probleme(const Probleme_FTD_IJK_base &)
void initialize(const Probleme_FTD_IJK_base &ijk_ft, const Domaine_IJK &)
void postraiter_thermique(const double t) const
void update_stat(Probleme_FTD_IJK_base &cas, const double dt)
double compute_desequil_alpha(const Domaine_IJK &geom_NS, const double portee_wall_repulsion) const
double calculer_vraie_dissipation(const double &pseudo_dissip, const double &duidx, const double &duidy, const double &duidz, const double &dujdx, const double &dujdy, const double &dujdz, const double &dukdx, const double &dukdy, const double &dukdz) const
bool has_champ_vectoriel(const Motcle &nom) const override
bool has_champ(const Motcle &nom) const override
void compute_vecA_minus_vecB_in_vecA(IJK_Field_vector3_double &vecA, const IJK_Field_vector3_double &vecB)
void postraiter(Sortie &, int flag_valeur_instantanee=0) const
double face_to_cell_gradient(const IJK_Field_double &vitesse_i, const IJK_Field_double &vitesse_j, const IJK_Field_double &vitesse_k, const int i, const int j, const int k, const double dz, double &duidx, double &dujdx, double &dukdx, double &duidy, double &dujdy, double &dukdy, double &duidz, double &dujdz, double &dukdz, const bool on_the_first_cell, const bool on_the_last_cell, const int bc_type) const
double calculer_produit_scalaire_faces_to_center(const IJK_Field_double &ui, const IJK_Field_double &uj, const IJK_Field_double &uk, const IJK_Field_double &vi, const IJK_Field_double &vj, const IJK_Field_double &vk, const int i, const int j, const int k)
std::map< Motcle, IJK_Field_vector3_double > vect_post_fields_
double calculer_gradients_vitesse(const IJK_Field_double &vitesse_i, const IJK_Field_double &vitesse_j, const IJK_Field_double &vitesse_k, const int i, const int j, const int k, const double dz, double &duidx, double &dujdx, double &dukdx, double &duidy, double &dujdy, double &dukdy, double &duidz, double &dujdz, double &dukdz, double &dduidxx, double &ddujdxx, double &ddukdxx, double &dduidyy, double &ddujdyy, double &ddukdyy, double &dduidzz, double &ddujdzz, double &ddukdzz, const bool on_the_first_cell, const bool on_the_last_cell) const
int lire_motcle_non_standard(const Motcle &mot, Entree &is) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
Champs_compris_IJK champs_compris_
the actual fields registered and managed by the post-processing part (=all the extra fields,...
void compute_and_store_gradU_cell(const IJK_Field_double &vitesse_i, const IJK_Field_double &vitesse_j, const IJK_Field_double &vitesse_k)
void cell_to_cell_gradient(const int i, const int j, const int k, const IJK_Field_double &dudx, const IJK_Field_double &dvdy, const IJK_Field_double &dwdx, const IJK_Field_double &dudz, const IJK_Field_double &dvdz, const IJK_Field_double &dwdz, double &ddudxy, double &ddudxz, double &ddudyz, double &ddvdxy, double &ddvdxz, double &ddvdyz, double &ddwdxy, double &ddwdxz, double &ddwdyz) const