TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Terme_Source_Canal_RANS_LES_VDF_Face.cpp
1/****************************************************************************
2* Copyright (c) 2024, 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 <Terme_Source_Canal_RANS_LES_VDF_Face.h>
17#include <math.h>
18#include <Champ_Uniforme.h>
19#include <Domaine_VDF.h>
20#include <Domaine_Cl_VDF.h>
21#include <Pb_Hydraulique.h>
22#include <Pb_Thermohydraulique.h>
23#include <EFichier.h>
24#include <Interprete.h>
25#include <SFichier.h>
26
27Implemente_instanciable_sans_destructeur(Terme_Source_Canal_RANS_LES_VDF_Face,"Canal_RANS_LES_VDF_Face",Source_base);
28
29// Terme_Source_Canal_RANS_LES_VDF_Face::Terme_Source_Canal_RANS_LES_VDF_Face()
30// {
31// dernier_tps_calc = new double();
32// }
33
35{
36 //Le destructeur est appele a l'initialisation alors
37 //la sauvegarde du champ se fait hors initialisation
38
39 if(umoy.size()!=0)
40 {
41 SFichier vit_sauv ("utemp_sum.dat");
42 vit_sauv << utemp_sum;
43 SFichier vit_sauv2 ("umoy.dat");
44 vit_sauv2 << umoy;
45 }
46}
47
48//// printOn
49//
50
52{
53 return s << que_suis_je() ;
54}
55
56
57//// readOn
58//
59
61{
62 int compteur = 0;
63 Motcle mot_lu;
64 Motcle acc_ouverte("{");
65 Motcle acc_fermee("}");
66
67 nom_pb_rans="non_couple";
68 // 0 => moyenne spatiale
69 // 1 => moyenne temporelle glissante (moyenne en alpha)
70 // 2 => moyenne temporelle commence a t=f_start-t_av
71 // 3 => moyenne temporelle commence a t=0 (mettre valeur positive a t_av dans .data)
72
73 Motcles les_mots(6);
74 {
75 les_mots[0]="alpha_tau"; //coeff de relaxation du terme source
76 les_mots[1]="Ly"; //Hauteur du canal plan (utile pour la moyenne spatiale)
77 les_mots[2]="f_start"; //temps a partir duquel le terme source est active
78 les_mots[3]="t_av"; //temps de prise de moyenne (moyenne temporelle glissante)
79 les_mots[4]="type_moyenne"; //temps de prise de moyenne (moyenne temporelle glissante)
80 les_mots[5]="nom_pb_rans";
81 }
82 is >> mot_lu;
83 if(mot_lu != acc_ouverte)
84 {
85 Cerr << "On attendait { a la place de " << mot_lu
86 << " lors de la lecture des parametres de la loi de paroi " << finl;
87 }
88 is >> mot_lu;
89 while(mot_lu != acc_fermee)
90 {
91 int rang=les_mots.search(mot_lu);
92 switch(rang)
93 {
94 case 0 :
95 is >> alpha_tau;
96 Cerr << "alpha_tau = " << alpha_tau << finl;
97 compteur++;
98 break;
99 case 1 :
100 is >> Ly;
101 Cerr << "Ly = "<< Ly << finl;
102 compteur++;
103 break;
104 case 2 :
105 is >> f_start;
106 Cerr << "f_start = "<< f_start << finl;
107 compteur++;
108 break;
109 case 3 :
110 is >> t_av;
111 Cerr << "t_av = " << t_av << finl;
112 Cerr << "type_moyenne : " << moyenne << finl;
113 if((moyenne==1)&&(t_av<=0))
114 {
115 Cerr << "La periode en temps de la moyenne temporelle glissante est non precisee !"
116 << finl;
117 exit();
118 }
119 compteur++;
120 break;
121 case 4 :
122 is >> moyenne;
123 Cerr << "type_moyenne : " << moyenne << finl;
124 if((moyenne==1)&&(t_av<=0))
125 {
126 Cerr << "La periode en temps de la moyenne temporelle glissante est non precisee !"
127 << finl;
128 exit();
129 }
130 compteur++;
131 break;
132 case 5 :
133 is >> nom_pb_rans;
134
135 compteur++;
136 break;
137 default :
138 {
139 Cerr << mot_lu << " n'est pas un mot compris" << finl;
140 Cerr << "Les mots compris sont : " << les_mots << finl;
141 exit();
142 }
143 }
144 is >> mot_lu;
145 }
146 Cerr << "nom_pb_rans = " << nom_pb_rans << finl;
147 Cerr << compteur << " motcles ont ete lus dans le readOn" << finl;
148
149 init();
150
151 return is;
152
153}
154
156 const Domaine_Cl_dis_base& domaine_Cl_dis)
157{
158 le_dom_VDF = ref_cast(Domaine_VDF, domaine_dis);
159 le_dom_Cl_VDF = ref_cast(Domaine_Cl_VDF, domaine_Cl_dis);
160}
161
166
168{
169 const Domaine_dis_base& zdisbase=mon_equation->inconnue().domaine_dis_base();
170 const Domaine_VDF& domaine_VDF=ref_cast(Domaine_VDF, zdisbase);
171 int nb_faces = domaine_VDF.nb_faces();
172 const double tps = mon_equation->schema_temps().temps_courant();
173
174
175 int nb_elems = domaine_VDF.nb_elem();
176 int nb_faces_y = domaine_VDF.nb_faces_Y();
177 int Ny = nb_elems/(nb_faces_y-nb_elems);
178
179 //SFichier fic_moyx("moyx.dat");
180 //SFichier fic_moyy("moyy.dat");
181 //SFichier fic_moyz("moyz.dat");
182 //SFichier fic_ransx("ransx.dat");
183 //SFichier fic_ransy("ransy.dat");
184 //SFichier fic_ransz("ransz.dat");
185
186 tau.resize(nb_faces);
187 U_RANS.resize(nb_faces);
188 umoy_spat.resize(nb_faces);
189 // champ_spat.resize(nb_faces);
190 umoy_spat.resize(nb_faces);
191 umoy_spat.resize(nb_faces);
192
193 if (nom_pb_rans == "non_couple")
194 {
195 SFichier fic_verif("verif.RANS");
196 EFichier fic_vit("vitesse_RANS.dat");
197
198 for(int num_face=0 ; num_face<nb_faces ; num_face++)
199 {
200 int face;
201 fic_vit >> face ;
202 fic_vit >> U_RANS(face) ;
203 fic_verif << face << " " << U_RANS(face) << finl;
204 }
205 }
206 for(int num_face=0 ; num_face<nb_faces ; num_face++)
207 {
208 tau(num_face) = alpha_tau;
209 }
210
211
212 umoy_x.resize(Ny);
213 umoy_y.resize(Ny+1);
214 umoy_z.resize(Ny);
215
216 umoy_x = 0.;
217 umoy_y = 0.;
218 umoy_z = 0.;
219
221
222 utemp_gliss.resize(nb_faces);
223 utemp.resize(nb_faces);
224
225 utemp_gliss = 0.;
226
227 utemp_gliss.resize(nb_faces);
228 utemp.resize(nb_faces);
229 utemp_sum.resize(nb_faces);
230
231 utemp_gliss = 0.;
232 utemp = 0.;
233 utemp_sum = 0.;
234
235 umoy.resize(nb_faces);
236
237 //Reprise
238 if (tps > f_start)
239 {
240 EFichier vit_umoy ("utemp_sum.dat");
241 EFichier vit_umoy2 ("umoy.dat");
242 SFichier vit_reprise ("LES.reprise");
243 SFichier vit_reprise2 ("LES2.reprise");
244
245 int trash;
246 vit_umoy >> trash;
247 Cerr << "trash = " << trash << finl;
248 vit_umoy2 >> trash;
249 Cerr << "trash = " << trash << finl;
250
251 for(int num_face = 0 ; num_face<nb_faces ; num_face++)
252 {
253 vit_umoy >> utemp_sum(num_face);
254 vit_umoy2 >> umoy(num_face);
255
256 vit_reprise << num_face << " " << utemp_sum(num_face) << finl;
257 vit_reprise2 << num_face << " " << umoy(num_face) << finl;
258 }
259 }
260}//fin init
261
263{
264 Cerr << "Debut initialisation des tableaux de moyenne spatiale..." << finl;
265
266 // On va initialiser les differents parametres membres de la classe
267 // utiles au calcul des differentes moyennes
268 // Initialisation de : Yu,Yv,Yw + compt_x,compt_y,compt_z
269 // + corresp_u,corresp_v,corresp_w
270 const Domaine_dis_base& zdisbase=mon_equation->inconnue().domaine_dis_base();
271 const Domaine_VDF& domaine_VDF=ref_cast(Domaine_VDF, zdisbase);
272
273 const DoubleTab& xv = domaine_VDF.xv();
274
275 int nb_faces = domaine_VDF.nb_faces();
276 int nb_elems = domaine_VDF.domaine().nb_elem();
277 int nb_faces_y = domaine_VDF.nb_faces_Y();
278
279 int Ny; //Nombre de maille en y + 1 (=nombre de faces)
280 Ny = nb_elems/(nb_faces_y-nb_elems);
281 Nxy = Ny;
282 Nyy = Ny+1;
283 Nzy = Ny;
284
285 const IntVect& orientation = domaine_VDF.orientation();
286
287 int num_face,j,ori,indic,indicv,indicw,trouve;
288 double y;
289 j=0;
290 indic = 0;
291 indicv = 0;
292 indicw = 0;
293
294 // dimensionnement aux valeurs rentrees dans le jeu de donnees
295 Yu.resize(Nxy);
296 Yv.resize(Nyy);
297 Yw.resize(Nzy);
298
299 compt_x.resize(Nxy);
300 compt_y.resize(Nyy);
301 compt_z.resize(Nzy);
302
303 // champ_spat.resize(nb_faces);
304
305 // champ_spat_x.resize(Nxy);
306 // champ_spat_y.resize(Nyy);
307 // champ_spat_z.resize(Nzy);
308
309 // champ_spat_x_2.resize(Nxy);
310 // champ_spat_y_2.resize(Nyy);
311 // champ_spat_z_2.resize(Nzy);
312
313 corresp_u.resize(nb_faces);
314 corresp_v.resize(nb_faces);
315 corresp_w.resize(nb_faces);
316
317 // initialisation
318 Yu = -100.;
319 Yv = -100.;
320 Yw = -100.;
321
322 compt_x = 0;
323 compt_y = 0;
324 compt_z = 0;
325
326 corresp_u = -1;
327 corresp_v = -1;
328 corresp_w = -1;
329
330 // remplissage des tableaux ci-dessus
331
332 // Pour le calcul de u, v, w sur les plans d hmogeneite
333 for (num_face=0; num_face<nb_faces; num_face++)
334 {
335 ori = orientation(num_face);
336 switch(ori)
337 {
338 case 0 :
339 {
340 y = xv(num_face,1);
341 trouve = 0;
342 for (j=0; j<indic+1; j++)
343 {
344 if(y==Yu[j])
345 {
346 corresp_u[num_face] = j;
347 compt_x[j] ++;
348 j=indic+1;
349 trouve = 1;
350 break;
351 }
352 }
353 if (trouve==0)
354 {
355 corresp_u[num_face] = indic;
356 Yu[indic]=y;
357 compt_x[indic] ++;
358 indic++;
359 }
360 break;
361 }
362 case 1 :
363 {
364 y = xv(num_face,1);
365 trouve = 0;
366 for (j=0; j<indicv+1; j++)
367 {
368 if(y==Yv[j])
369 {
370 corresp_v[num_face] = j;
371 compt_y[j] ++;
372 j = indicv+1;
373 trouve = 1;
374 break;
375 }
376 }
377 if (trouve==0)
378 {
379 corresp_v[num_face] = indicv;
380 Yv[indicv]=y;
381 compt_y[indicv] ++;
382 indicv++;
383 }
384 break;
385 }
386 case 2 :
387 {
388 y = xv(num_face,1);
389 trouve = 0;
390 for (j=0; j<indicw+1; j++)
391 {
392 if(y==Yw[j])
393 {
394 corresp_w[num_face] = j;
395 compt_z[j] ++;
396 j = indicw+1;
397 trouve = 1;
398 break;
399 }
400 }
401 if (trouve==0)
402 {
403 Yw[indicw]=y;
404 corresp_w[num_face] = indicw;
405 compt_z[indicw] ++;
406 indicw++;
407 }
408 break;
409 }
410 default :
411 {
412 Cerr << "Cas de figure impossible!!" << finl;
413 exit();
414 break;
415 }
416 }
417 }
418 Nxy = indic; // nombre de y pour umoy
419 Nyy = indicv;
420 Nzy = indicw;
421
422 Cerr << "Initialisation des tableaux de moyenne spatiale : OK" << finl;
423
424}
425
427{
428 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
429 int nb_elems = domaine_VDF.nb_elem();
430 int nb_faces = domaine_VDF.nb_faces();
431 const IntTab& face_voisins = domaine_VDF.face_voisins();
432 const IntTab& elem_faces = domaine_VDF.elem_faces();
433 DoubleVect norme_vit_elem(nb_elems);
434 norme_vit_elem.resize(nb_faces);
435 double vit1, vit2, vit3;
436 DoubleTab norme(nb_faces);
437
438 for(int num_elem = 0 ; num_elem<nb_elems ; num_elem++)
439 {
440 //WARNING 3D ONLY !!!!!!
441 int face1 = elem_faces(num_elem,0);
442 int face2 = elem_faces(num_elem,3);
443 int face3 = elem_faces(num_elem,1);
444 int face4 = elem_faces(num_elem,4);
445 int face5 = elem_faces(num_elem,2);
446 int face6 = elem_faces(num_elem,5);
447
448 // Cerr << "face1 = "<<face1<<finl;
449 // Cerr << "face2 = "<<face2<<finl;
450 // Cerr << "U_RANS(face1) = "<< U_RANS(face1) <<finl;
451 // Cerr << "U_RANS(face2) = "<< U_RANS(face2) <<finl;
452
453
454 vit1 = 0.5*(U_RANS(face1)+U_RANS(face2));
455 vit2 = 0.5*(U_RANS(face3)+U_RANS(face4));
456 vit3 = 0.5*(U_RANS(face5)+U_RANS(face6));
457
458 norme_vit_elem(num_elem) = sqrt(vit1*vit1 + vit2*vit2 + vit3*vit3);
459 }
460
461 //Redistribution de la norme sur les faces
462
463 for(int num_face = 0 ; num_face<nb_faces ; num_face++)
464 {
465 int elem0 = face_voisins(num_face,0);
466 int elem1 = face_voisins(num_face,1);
467
468 if(elem0==-1)
469 norme(num_face) = norme_vit_elem(elem1);
470 else if(elem1==-1)
471 norme(num_face) = norme_vit_elem(elem0);
472 else
473 norme(num_face) = 0.5*(norme_vit_elem(elem0)+norme_vit_elem(elem1));
474 }
475
476 return norme;
477}
478
479void Terme_Source_Canal_RANS_LES_VDF_Face::moy_spat(DoubleVect& champ, DoubleVect& champ_spat_x,
480 DoubleVect& champ_spat_y, DoubleVect& champ_spat_z,
481 DoubleVect& champ_spat_x_2, DoubleVect& champ_spat_y_2,
482 DoubleVect& champ_spat_z_2)
483{
484 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
485 const IntVect& orientation = domaine_VDF.orientation();
486 int nb_faces = domaine_VDF.nb_faces();
487
488 int j=-1;
489
490 compt_x=0;
491 compt_y=0;
492 compt_z=0;
493
494 for (int num_face=0; num_face<nb_faces; num_face++)
495 {
496 int ori = orientation[num_face];
497
498 if(ori==0)
499 {
500 j=corresp_u[num_face];
501 champ_spat_x[j] += champ(num_face);
502 champ_spat_x_2[j] += champ(num_face)*champ(num_face);
503 compt_x[j]++;
504 }
505
506 if(ori==1)
507 {
508 j=corresp_v[num_face];
509 champ_spat_y[j] += champ(num_face);
510 champ_spat_y_2[j] += champ(num_face)*champ(num_face);
511 compt_y[j]++;
512 }
513
514 if(ori==2)
515 {
516 j=corresp_w[num_face];
517 champ_spat_z[j] += champ(num_face);
518 champ_spat_z_2[j] += champ(num_face)*champ(num_face);
519 compt_z[j]++;
520 }
521 }
522
523 for(j=0; j<Nxy; j++)
524 {
525 champ_spat_x[j] = champ_spat_x[j]/compt_x[j];
526 champ_spat_x_2[j] = champ_spat_x_2[j]/compt_x[j];
527 }
528
529 for(j=0; j<Nyy; j++)
530 {
531 champ_spat_y[j] = champ_spat_y[j]/compt_y[j];
532 champ_spat_y_2[j] = champ_spat_y_2[j]/compt_y[j];
533 }
534
535 for(j=0; j<Nzy; j++)
536 {
537 champ_spat_z[j] = champ_spat_z[j]/compt_z[j];
538 champ_spat_z_2[j] = champ_spat_z_2[j]/compt_z[j];
539 }
540}
541
542void Terme_Source_Canal_RANS_LES_VDF_Face::calculer_integrale_temporelle(DoubleVect& moy_temp, const DoubleVect& moy)
543{
544 double dt_v = mon_equation->schema_temps().pas_de_temps();
545 if(je_suis_maitre())
546 moy_temp.ajoute(dt_v,moy);
547}
548
549void Terme_Source_Canal_RANS_LES_VDF_Face::ecriture_moy_spat(DoubleVect& champ_spat_x, DoubleVect& champ_spat_y,
550 DoubleVect& champ_spat_z, DoubleVect& champ_spat_x_2,
551 DoubleVect& champ_spat_y_2, DoubleVect& champ_spat_z_2)
552{
553 const double tps = mon_equation->schema_temps().temps_courant();
554 Nom temps = Nom(tps);
555
556 Nom nom_ficx ="moy_spat_force_x_";
557 nom_ficx+= temps;
558 nom_ficx+= ".dat";
559 SFichier ficx(nom_ficx);
560
561 Nom nom_ficy ="moy_spat_force_y_";
562 nom_ficy+= temps;
563 nom_ficy+= ".dat";
564 SFichier ficy(nom_ficy);
565
566 Nom nom_ficz ="moy_spat_force_z_";
567 nom_ficz+= temps;
568 nom_ficz+= ".dat";
569 SFichier ficz(nom_ficz);
570
571 for(int i=0 ; i<Nxy ; i++)
572 ficx << Yu[i] << " " << champ_spat_x[i] << " " << sqrt(champ_spat_x_2[i]-champ_spat_x[i]*champ_spat_x[i]) << finl;
573 for(int i=0 ; i<Nyy ; i++)
574 ficy << Yv[i] << " " << champ_spat_y[i] << " " << sqrt(champ_spat_y_2[i]-champ_spat_y[i]*champ_spat_y[i]) << finl;
575 for(int i=0 ; i<Nzy ; i++)
576 ficz << Yw[i] << " " << champ_spat_z[i] << " " << sqrt(champ_spat_z_2[i]-champ_spat_z[i]*champ_spat_z[i]) << finl;
577}
578
579void Terme_Source_Canal_RANS_LES_VDF_Face::ecriture_moy_temp(DoubleVect& champ_temp_x, DoubleVect& champ_temp_y,
580 DoubleVect& champ_temp_z, DoubleVect& champ_temp_x_2,
581 DoubleVect& champ_temp_y_2, DoubleVect& champ_temp_z_2,
582 const double dt)
583{
584 const double tps = mon_equation->schema_temps().temps_courant();
585 Nom temps = Nom(tps);
586
587 Nom nom_ficx ="moy_temp_force_x_";
588 nom_ficx+= temps;
589 nom_ficx+= ".dat";
590 SFichier ficx(nom_ficx);
591
592 Nom nom_ficy ="moy_temp_force_y_";
593 nom_ficy+= temps;
594 nom_ficy+= ".dat";
595 SFichier ficy(nom_ficy);
596
597 Nom nom_ficz ="moy_temp_force_z_";
598 nom_ficz+= temps;
599 nom_ficz+= ".dat";
600 SFichier ficz(nom_ficz);
601
602 for(int i=0 ; i<Nxy ; i++)
603 ficx << Yu[i] << " " << champ_temp_x[i]/dt << " " << sqrt(champ_temp_x_2[i]-champ_temp_x[i]*champ_temp_x[i])/dt << finl;
604 for(int i=0 ; i<Nyy ; i++)
605 ficy << Yv[i] << " " << champ_temp_y[i]/dt << " " << sqrt(champ_temp_y_2[i]-champ_temp_y[i]*champ_temp_y[i])/dt << finl;
606 for(int i=0 ; i<Nzy ; i++)
607 ficz << Yw[i] << " " << champ_temp_z[i]/dt << " " << sqrt(champ_temp_z_2[i]-champ_temp_z[i]*champ_temp_z[i])/dt << finl;
608}
609
611{
612 //Cerr << "Je suis dans le mettre_a_jour" << finl;
613 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
614
615 const DoubleTab& vitesse = mon_equation->inconnue().valeurs();
616 const double dt = mon_equation->schema_temps().pas_de_temps();
617 const double tps = mon_equation->schema_temps().temps_courant();
618 const double dt_min = mon_equation->schema_temps().pas_temps_min();
619
620 const IntVect& orientation = domaine_VDF.orientation();
621 int nb_elems = domaine_VDF.nb_elem();
622 int nb_faces = domaine_VDF.nb_faces();
623
624 int nb_faces_x = domaine_VDF.nb_faces_X();
625 int nb_faces_y = domaine_VDF.nb_faces_Y();
626 int nb_faces_z = domaine_VDF.nb_faces_Z();
627
628 //static double duree = 0;
629
630 double vit;
631 SFichier fic_moyx("moyx.dat");
632 SFichier fic_moyy("moyy.dat");
633 SFichier fic_moyz("moyz.dat");
634 SFichier fic_f("f.dat", ios::app);
635 SFichier fic_fface("force_face.dat", ios::app);
636 const DoubleTab& xv = domaine_VDF.xv();
637
638
639
640 //****************************************************
641 //******* MaJ de la vitesse cible (RANS) ***********
642 //**************************************************
643 if(nom_pb_rans != "non_couple")
644 {
645 OBS_PTR(Probleme_base) pb_rans;
646 Objet_U& obj=Interprete::objet(nom_pb_rans);
647
648 if( sub_type(Probleme_base, obj) )
649 {
650 pb_rans = ref_cast(Probleme_base, obj);
651 }
652
653 U_RANS = pb_rans->equation(0).inconnue().valeurs();
654
655 }
656 //***************************************************
657
658 if(moyenne==0)
659 {
660 if(tps>(f_start))
661 {
662 //******************************************************
663 //*************** MOYENNE SPATIALE *********************
664 //******************************************************
665
666 int num_face,Nx, Ny, Nz, j;
667 umoy_x = 0;
668 umoy_y = 0;
669 umoy_z = 0;
670
671 Nx = nb_elems/(nb_faces_x-nb_elems); //nbre de maille
672 Ny = nb_elems/(nb_faces_y-nb_elems);
673 Nz = nb_elems/(nb_faces_z-nb_elems);
674
675 for (num_face=0; num_face<nb_faces; num_face++)
676 {
677 int ori = orientation[num_face];
678 vit = vitesse[num_face];
679 if (ori==0)
680 {
681 j=corresp_u[num_face];
682 umoy_x[j] += vit/((Nx+1)*Nz);
683 }
684 else if (ori==1)
685 {
686 j=corresp_v[num_face];
687 umoy_y[j] += vit/(Nx*Nz);
688 }
689 else if (ori==2)
690 {
691 j=corresp_w[num_face];
692 umoy_z[j] += vit/(Nx*(Nz+1));
693 }
694 }
695
696 // // Ecriture des profils moyens
697 int i;
698 for (i=0 ; i<Ny ; i++)
699 {
700 fic_moyx << i << " " << umoy_x(i) << finl;
701 }
702 for (i=0 ; i<Ny+1 ; i++)
703 {
704 fic_moyy << i << " " << umoy_y(i) << finl;
705 }
706 for (i=0 ; i<Ny ; i++)
707 {
708 fic_moyz << i << " " << umoy_z(i) << finl;
709 }
710
711 // Redistribution du profil moyen spatial sur tout le champ de vitesse
712
713 for(num_face = 0 ; num_face<nb_faces ; num_face++)
714 {
715 int ori = orientation(num_face);
716
717 if (ori==0)
718 {
719 j = corresp_u(num_face);
720 umoy(num_face) = umoy_x(j);
721 }
722 else if (ori==1)
723 {
724 j = corresp_v(num_face);
725 umoy(num_face) = umoy_y(j);
726 }
727 else if (ori==2)
728 {
729 j = corresp_w(num_face);
730 umoy(num_face) = umoy_z(j);
731 }
732 else
733 Cerr << "Erreur orientation dans la redistribution de la moy. spat. de SouCanRANS_LESVDFFa" << finl;
734
735 }
736
737 //******************************************************
738 //************** FIN MOYENNE SPATIALE ******************
739 //******************************************************
740 }
741 }//fin if moyenne = 0
742
743 else if(moyenne==1)
744 {
745 //****************************************************************
746 //*************** MOYENNE TEMPORELLE GLISSANTE *******************
747 //****************************************************************
748
749 //Initialisation de la moyenne temp glissante
750 //avec une moyenne temporelle classique
751
752
753 // if(/*(tps>(f_start-t_av))&&*/(tps<f_start))
754 // {
755 // if (cpt==0)
756 // {
757 // for (int num_face=0;num_face<nb_faces;num_face++)
758 // {
759 // utemp_gliss(num_face) = vitesse(num_face)*(tps-(f_start-t_av));
760 // umoy(num_face) = vitesse(num_face);
761 // }
762 // cpt = 1;
763 // }
764 // else
765 // {
766 // for (int num_face=0;num_face<nb_faces;num_face++)
767 // {
768 // utemp_gliss(num_face) += vitesse(num_face)*dt;
769 // umoy(num_face) = utemp_gliss(num_face)/(tps-(f_start-t_av));
770 // }
771 // }
772
773 // }
774
775
776 //Calcul d'une premiere moyenne temporelle significative
777
778 // if((tps>(f_start-t_av))&&(tps<f_start))
779 // {
780 // if (cpt==0)
781 // {
782 // for (int num_face=0;num_face<nb_faces;num_face++)
783 // {
784 // utemp_sum(num_face) = vitesse(num_face)*(tps-(f_start-t_av)+dt);
785 // umoy(num_face) = vitesse(num_face);
786 // }
787 // cpt = 3;
788 // }
789 // else
790 // {
791 // for (int num_face=0;num_face<nb_faces;num_face++)
792 // {
793 // utemp_sum(num_face) += vitesse(num_face)*dt;
794 // umoy(num_face) = utemp_sum(num_face)/(tps-(f_start-t_av)+dt);
795 // }
796 // }
797 // //fic_f << tps << " " << utemp_sum(123) << " " << umoy(123) << finl;
798
799 // }
800 //Calcul de la moyenne temporelle glissante
801 //t_av est le temps d'integration de la moyenne temp glissante
802
803 // if(tps>=f_start)
804 // {
805 // for (int num_face=0 ; num_face<nb_faces ; num_face++)
806 // {
807 // if(cpt==1)
808 // {
809 // umoy(num_face) = vitesse(num_face);
810 // }
811 // else if(tps <= t_av)
812 // {
813 // umoy(num_face) = (dt/tps)*vitesse(num_face)
814 // +(1-(dt/tps))*umoy(num_face);
815 // cpt=2;
816 // }
817 // else
818 // {
819 // umoy(num_face) = (dt/t_av)*vitesse(num_face)
820 // +(1-(dt/t_av))*umoy(num_face);
821 // }
822 // }
823 // }
824
825
826 if(tps>dt_min)
827 {
828 if (cpt==0)
829 {
830 for (int num_face=0; num_face<nb_faces; num_face++)
831 {
832 utemp_sum(num_face) = vitesse(num_face)*dt;
833 umoy(num_face) = vitesse(num_face);
834 }
835 cpt++;
836 }
837 else if (tps <= t_av)
838 {
839 for (int num_face=0; num_face<nb_faces; num_face++)
840 {
841 utemp_sum(num_face) += vitesse(num_face)*dt;
842 umoy(num_face) = utemp_sum(num_face)/(tps-dt_min);
843 }
844 }
845 else
846 {
847 for (int num_face=0; num_face<nb_faces; num_face++)
848 {
849 umoy(num_face) = (dt/t_av)*vitesse(num_face)
850 +(1-(dt/t_av))*umoy(num_face);
851 }
852 }
853
854 fic_f << tps << " " << dt << " " << tps-dt_min << " "<< utemp_sum(123);
855 fic_f << " " << umoy(123) << " " << vitesse(123) << " " << U_RANS(123) << finl;
856
857 }
858
859 //****************************************************************
860 //*************** FIN MOYENNE TEMPORELLE GLISSANTE ***************
861 //****************************************************************
862
863 }// fin moyenne = 1
864 else if(moyenne==2)
865 {
866 //******************************************************
867 //*************** MOYENNE TEMPORELLE *******************
868 //******************************************************
869
870 //Calcul d'une premiere moyenne temporelle significative
871
872 // if((tps>=f_start)/*&&(tps<f_start)*/)
873 // {
874 // if (cpt==0)
875 // {
876 // for (int num_face=0;num_face<nb_faces;num_face++)
877 // {
878 // //utemp_sum(num_face) = vitesse(num_face)*(tps-f_start+dt);
879 // umoy(num_face) = vitesse(num_face);
880 // }
881 // cpt = 3;
882 // }
883 // else
884 // {
885 // for (int num_face=0;num_face<nb_faces;num_face++)
886 // {
887 // utemp_sum(num_face) += vitesse(num_face)*dt;
888 // umoy(num_face) = utemp_sum(num_face)/(tps-f_start+dt);
889 // }
890 // }
891 //fic_f << tps << " " << utemp_sum(123) << " " << umoy(123) << finl;
892
893 // }
894
895 if(tps>dt_min)
896 {
897 if (cpt==0)
898 {
899 for (int num_face=0; num_face<nb_faces; num_face++)
900 {
901 utemp_sum(num_face) = vitesse(num_face)*dt;
902 umoy(num_face) = vitesse(num_face);
903 }
904 cpt++;
905 }
906 else
907 {
908 for (int num_face=0; num_face<nb_faces; num_face++)
909 {
910 utemp_sum(num_face) += vitesse(num_face)*dt;
911 umoy(num_face) = utemp_sum(num_face)/(tps-dt_min);
912 }
913 }
914 fic_f << tps << " " << dt << " " << tps-dt_min << " "<< utemp_sum(123);
915 fic_f << " " << umoy(123) << " " << vitesse(123) << " " << U_RANS(123) << finl;
916
917 }
918
919
920
921
922 // if(tps>=f_start)
923 // {
924 // for (int num_face=0;num_face<nb_faces;num_face++)
925 // {
926 // utemp_sum(num_face) += dt*vitesse(num_face);
927 // umoy(num_face)=utemp_sum(num_face)/(tps-f_start+dt);
928 // }
929 // cpt=4;
930 // }
931
932
933 //***********************************************************
934 //*************** FIN MOYENNE TEMPORELLE ***************
935 //*******************************************************
936
937
938 }
939 else if(moyenne==3)
940 {
941 if (tps>0)
942 {
943 if (cpt==0)
944 {
945 for (int num_face=0; num_face<nb_faces; num_face++)
946 {
947 // duree=100*dt;
948 // utemp_sum(num_face) = vitesse(num_face)*duree;
949 // duree-=dt;
950 // umoy(num_face) = utemp_sum(num_face)/(tps-dt_min+duree);
951 // Cerr << "nb_faces = "<< nb_faces << finl;
952 // Cerr << "num_face = "<< num_face << finl;
953 utemp_sum(num_face) = vitesse(num_face)*dt;
954 umoy(num_face) = vitesse(num_face);
955 }
956 cpt ++;
957 }
958 else
959 {
960 for (int num_face=0; num_face<nb_faces; num_face++)
961 {
962 utemp_sum(num_face) += vitesse(num_face)*dt;
963 umoy(num_face) = utemp_sum(num_face)/(tps-dt_min+dt);
964 }
965 }
966 }
967 }
968 else
969 {
970 Cerr << "moyenne = " << moyenne << finl;
971 Cerr << "Probleme pour le choix du type de moyenne" << finl;
972 }
973
974 compteur_reprise++;
975
976 const DoubleVect& porosite_surf = equation().milieu().porosite_face();
977 const DoubleVect& volumes_entrelaces = domaine_VDF.volumes_entrelaces();
978
979 DoubleVect force(nb_faces);
980
981 DoubleVect force_x(nb_faces);
982 DoubleVect force_y(nb_faces);
983 DoubleVect force_z(nb_faces);
984
985 DoubleVect force_x_2(nb_faces);
986 DoubleVect force_y_2(nb_faces);
987 DoubleVect force_z_2(nb_faces);
988
989 DoubleVect force_x_temp(nb_faces);
990 DoubleVect force_y_temp(nb_faces);
991 DoubleVect force_z_temp(nb_faces);
992
993 DoubleVect force_x_temp_2(nb_faces);
994 DoubleVect force_y_temp_2(nb_faces);
995 DoubleVect force_z_temp_2(nb_faces);
996
997 double vol=0;
998 DoubleTab norm=norme_vit();
999
1000 for(int num_face = 0 ; num_face<nb_faces ; num_face++)
1001 {
1002 vol = volumes_entrelaces(num_face)*porosite_surf(num_face);
1003
1004 force(num_face) = ((U_RANS(num_face)-umoy(num_face))/(dt))
1005 *(U_RANS(num_face)/norm(num_face))*vol;
1006
1007 }
1008
1009 moy_spat(force, force_x, force_y, force_z, force_x_2, force_y_2, force_z_2);
1010
1011 calculer_integrale_temporelle(force_x_temp, force_x);
1012 calculer_integrale_temporelle(force_y_temp, force_y);
1013 calculer_integrale_temporelle(force_z_temp, force_z);
1014 calculer_integrale_temporelle(force_x_temp_2, force_x_2);
1015 calculer_integrale_temporelle(force_y_temp_2, force_y_2);
1016 calculer_integrale_temporelle(force_z_temp_2, force_z_2);
1017
1018 // Voir Schema_Temps_base::limpr pour information sur epsilon et modf
1019 double i, j, epsilon = 1.e-8;
1020 double dt_post_inst=2.;
1021 double tps_stats=10.;
1022
1023 modf((tps+dt)/dt_post_inst + epsilon, &i);
1024 modf(tps/dt_post_inst + epsilon, &j);
1025 if(i>j)
1026 {
1027 ecriture_moy_spat(force_x, force_y, force_z, force_x_2, force_y_2, force_z_2);
1028 if(tps>tps_stats)
1029 ecriture_moy_temp(force_x_temp, force_y_temp, force_z_temp,
1030 force_x_temp_2, force_y_temp_2, force_z_temp_2, tps-tps_stats);
1031 }
1032
1033 //////////////////////////////////////////
1034 //developpement a nettoyer apres la these
1035 //////////////////////////////////////////
1036
1037 int num_face = 123;
1038 static int cpt_sonde=0;
1039 if(cpt_sonde==0)
1040 {
1041 for(int face = 0 ; face<nb_faces ; face++)
1042 {
1043 if((xv(face,1) > 1.1)&&(xv(face,1) < 0.9)&&(orientation(face)==0))
1044 num_face=face;
1045 }
1046 fic_fface << "# " << num_face << " "<< orientation(num_face)
1047 << " " << xv(num_face,0) << " " << xv(num_face,1) << " " << xv(num_face,2) << finl;
1048 fic_f << "# " << num_face << " "<< orientation(num_face)
1049 << " " << xv(num_face,0) << " " << xv(num_face,1) << " " << xv(num_face,2) << finl;
1050 cpt_sonde++;
1051 }
1052
1053 fic_fface << tps << " " << vitesse(num_face) << " " << force(num_face) << finl;
1054 //////////////////////////////////////////////////////////
1055
1056}//fin mettre_a_jour
1057
1058
1059
1060void Terme_Source_Canal_RANS_LES_VDF_Face::ajouter_blocs(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl) const
1061{
1062 const Domaine_VDF& domaine_VDF = le_dom_VDF.valeur();
1063 int nb_faces = domaine_VDF.nb_faces();
1064 const DoubleVect& porosite_surf = equation().milieu().porosite_face();
1065 const DoubleVect& volumes_entrelaces = domaine_VDF.volumes_entrelaces();
1066 const double tps = mon_equation->schema_temps().temps_courant();
1067 const double dt = mon_equation->schema_temps().pas_de_temps();
1068 const double dt_min = mon_equation->schema_temps().pas_temps_min();
1069
1070 // const DoubleTab& vitesse = mon_equation->inconnue().valeurs();
1071
1072 double vol=0;
1073 //SFichier fic_f("f.dat", ios::app);
1074 //SFichier fic_finst("finst.dat", ios::app);
1075
1076 //double mbf = 0.; //terme de forcage instantanne ponctuel
1077 // double mbf2 = 0.; //terme de forcage moyenne
1078
1079 // int cpt2=0;
1080 // static double beta=1000;
1081
1082 // Calcul de la norme des vitesses au centre des elements
1083
1084 DoubleTab norm=norme_vit();
1085
1086 if((tps>dt_min)&&(tps>f_start))
1087 {
1088
1089 for(int num_face = 0 ; num_face<nb_faces ; num_face++)
1090 {
1091 vol = volumes_entrelaces(num_face)*porosite_surf(num_face);
1092
1093 secmem(num_face) += ((U_RANS(num_face)-umoy(num_face))/(tau(num_face)*dt))
1094 *(U_RANS(num_face)/norm(num_face))*vol;
1095
1096 // mbf2 += (((U_RANS(num_face)-umoy(num_face))
1097 // *(U_RANS(num_face)/norm(num_face)))
1098 // /(tau(num_face)*dt))*vol;
1099 // cpt2++;
1100 }
1101
1102 //mbf2 = mbf2/cpt2;
1103
1104 }//fin if f_start
1105
1106}
1107
1108DoubleTab& Terme_Source_Canal_RANS_LES_VDF_Face::calculer(DoubleTab& resu) const
1109{
1110 resu = 0;
1111 return ajouter(resu);
1112}
1113
int_t nb_elem() const
Definition Domaine.h:131
class Domaine_Cl_VDF
classe Domaine_Cl_dis_base Les objets Domaine_Cl_dis_base representent les conditions aux limites
class Domaine_VDF
Definition Domaine_VDF.h:64
int nb_faces_Z() const
int nb_faces_Y() const
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
int nb_faces_X() const
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
DoubleVect & volumes_entrelaces()
Definition Domaine_VF.h:99
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
int elem_faces(int i, int j) const
renvoie le numero de le ieme face de la maille num_elem la facon dont ces faces sont numerotees est
Definition Domaine_VF.h:543
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
const Domaine & domaine() const
Fichier en lecture Cette classe est a la classe C++ ifstream ce que la classe Entree est a la.
Definition EFichier.h:29
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
virtual const Milieu_base & milieu() const =0
static Objet_U & objet(const Nom &)
Voir Interprete_bloc::objet_global() BM: la classe Interprete n'est pas le meilleur endroit pour cett...
DoubleVect & porosite_face()
Definition Milieu_base.h:62
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
class Nom Une chaine de caractere pour nommer les objets de TRUST
Definition Nom.h:31
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
Objet_U()
Constructeur par defaut : attribue un numero d'identifiant unique a l'objet (object_id_),...
Definition Objet_U.cpp:55
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
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
Classe de base des flux de sortie.
Definition Sortie.h:52
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
Definition Source_base.h:42
virtual DoubleTab & ajouter(DoubleTab &) const
void ajoute(_SCALAR_TYPE_ alpha, const TRUSTVect &y, Mp_vect_options opt=VECT_ALL_ITEMS)
Definition TRUSTVect.tpp:52
void resize(_SIZE_, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTVect.tpp:91
class Terme_Source_Canal_RANS_LES_VDF_Face Cette classe concerne un terme source calcule en partie gr...
void moy_spat(DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &)
void calculer_integrale_temporelle(DoubleVect &, const DoubleVect &)
void ajouter_blocs(matrices_t matrices, DoubleTab &secmem, const tabs_t &semi_impl) const override
void associer_domaines(const Domaine_dis_base &, const Domaine_Cl_dis_base &) override
void mettre_a_jour(double) override
DOES NOTHING - to override in derived classes.
OBS_PTR(Domaine_VDF) le_dom_VDF
void ecriture_moy_temp(DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, const double)
void ecriture_moy_spat(DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &, DoubleVect &)