17#include <ParoiVDF_TBLE.h>
18#include <Domaine_Cl_VDF.h>
19#include <Dirichlet_paroi_fixe.h>
20#include <Dirichlet_paroi_defilante.h>
21#include <Champ_Uniforme.h>
22#include <Fluide_base.h>
25#include <Schema_Temps_base.h>
26#include <Operateur_Grad.h>
27#include <Probleme_base.h>
28#include <Modele_turbulence_scal_base.h>
29#include <Paroi_TBLE_scal_VDF.h>
30#include <Terme_Boussinesq_VDF_Face.h>
32#include <EcrFicPartage.h>
33#include <Navier_Stokes_std.h>
34#include <Modele_turbulence_hyd_base.h>
49 Cerr <<
"Reading of data for a " <<
que_suis_je() <<
" wall law" << finl;
52 param.lire_avec_accolades_depuis(is);
90 Cerr <<
"debut init_lois_paroi" << finl;
93 const IntVect& orientation = domaine_VDF.
orientation();
95 const IntTab& elem_faces = domaine_VDF.
elem_faces();
96 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
99 int compteur_faces_paroi = 0;
108 corresp.resize(le_dom_dis_->nb_faces_bord());
110 compteur_faces_paroi = 0;
112 for (
int n_bord = 0; n_bord < domaine_VDF.
nb_front_Cl(); n_bord++)
114 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
120 const int nfin = ndeb + le_bord.
nb_faces();
122 for (
int num_face = ndeb; num_face < nfin; num_face++)
124 corresp[num_face]=compteur_faces_paroi;
125 compteur_faces_paroi++;
130 compteur_faces_paroi = 0;
132 for (
int n_bord = 0; n_bord < domaine_VDF.
nb_front_Cl(); n_bord++)
138 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
144 const int nfin = ndeb + le_bord.
nb_faces();
149 for (
int num_face = ndeb; num_face < nfin; num_face++)
152 int ori = orientation(num_face);
154 if ((elem = face_voisins(num_face, 0)) == -1)
155 elem = face_voisins(num_face, 1);
166 eq_vit[corresp[num_face]].set_y0(0.);
167 eq_vit[corresp[num_face]].set_yn(dist);
170 eq_vit[corresp[num_face]].set_u_y0(0,0.);
172 vmoy = 0.5*(vit(elem_faces(elem, (ori + 1) % 4)) + vit(elem_faces(elem, (ori + 3) % 4)));
174 eq_vit[corresp[num_face]].set_u_yn(0, vmoy);
178 eq_vit[corresp[num_face]].set_Uinit_lin(0, 0., vmoy);
180 for (
int itble = 0; itble <
nb_pts + 1; itble++)
181 eq_vit[corresp[num_face]].set_Uinit(0, itble,
valeurs_reprises(corresp[num_face], 0, itble));
183 compteur_faces_paroi++;
189 for (
int num_face = ndeb; num_face < nfin; num_face++)
192 const int ori = orientation(num_face);
194 if ((elem = face_voisins(num_face, 0)) == -1)
195 elem = face_voisins(num_face, 1);
203 eq_vit[corresp[num_face]].set_y0(0.);
204 eq_vit[corresp[num_face]].set_yn(dist);
208 eq_vit[corresp[num_face]].set_u_y0(0, 0.);
209 eq_vit[corresp[num_face]].set_u_y0(1, 0.);
211 vmoy = 0.5*(vit(elem_faces(elem, (ori + 1) % 6)) + vit(elem_faces(elem, (ori + 4) % 6)));
213 eq_vit[corresp[num_face]].set_u_yn(0, vmoy);
217 eq_vit[corresp[num_face]].set_Uinit_lin(0, 0., vmoy);
219 for (
int itble = 0; itble <
nb_pts + 1; itble++)
220 eq_vit[corresp[num_face]].set_Uinit(0, itble,
valeurs_reprises(corresp[num_face], 0, itble));
222 vmoy = 0.5*(vit(elem_faces(elem, (ori + 2) % 6)) + vit(elem_faces(elem, (ori + 5) % 6)));
223 eq_vit[corresp[num_face]].set_u_yn(1, vmoy);
225 eq_vit[corresp[num_face]].set_Uinit_lin(1, 0., vmoy);
227 for (
int itble = 0; itble <
nb_pts + 1; itble++)
228 eq_vit[corresp[num_face]].set_Uinit(1, itble,
valeurs_reprises(corresp[num_face], 1, itble));
230 compteur_faces_paroi++;
243 int nb_sources=les_sources.size();
245 for (
int j=0; j<nb_sources; j++)
254 const DoubleTab& tab_champ_beta_t = ch_beta_t.
valeurs();
257 beta_t = std::max(tab_champ_beta_t(0,0),DMINFLOAT);
261 Cerr <<
" On ne sait pas traiter un beta_t non uniforme dans TBLE !!!"<< finl;
272 Cerr <<
"fin init_lois_paroi" << finl;
300 const IntVect& orientation = domaine_VDF.
orientation();
301 const IntTab& face_voisins = domaine_VDF.
face_voisins();
302 const IntTab& elem_faces = domaine_VDF.
elem_faces();
304 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
308 const DoubleTab& tab_visco = ch_visco_cin.
valeurs();
316 Cerr <<
" On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
325 visco = std::max(tab_visco(0,0),DMINFLOAT);
335 Cerr <<
"In ParoiVDF_TBLE::calculer_hyd_BiK : visco = " << tab_visco.
local_min_vect() <<
" <= 0 ? " << finl;
347 double gradient_de_pression0 = 0., gradient_de_pression1 = 0., vmoy = 0., ts0 =0., ts1=0.;
352 DoubleTab grad_p(vit);
355 gradient.calculer(p, grad_p);
356 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
359 DoubleTab termes_sources;
364 int nb_sources=les_sources.size();
365 for (
int j=0; j<nb_sources; j++)
379 int compteur_faces_paroi = 0;
384 compteur_faces_paroi = 0;
386 for (
int n_bord=0; n_bord<domaine_VDF.
nb_front_Cl(); n_bord++)
393 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
413 for (
int num_face=ndeb; num_face<nfin; num_face++)
415 ori = orientation(num_face);
420 if ((elem = face_voisins(num_face,0)) == -1)
422 elem = face_voisins(num_face,1);
429 eq_vit[corresp[num_face]].set_dt(1.e12);
431 eq_vit[corresp[num_face]].set_dt(dt);
436 int face1 = elem_faces(elem,(ori+1));
437 int face2 = elem_faces(elem,(ori+3)%4);
444 DoubleVect ts_boussi(
nb_pts+1);
449 const DoubleVect& gravite = ch_gravite.
valeurs();
463 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
470 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
473 double g_t=gt_vect(1-ori);
475 for(
int i=0; i<
nb_pts+1; i++)
486 vmoy = 0.5*(vit(face1) + vit(face2));
490 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
491 + termes_sources(face2)/volumes_entrelaces(face2));
493 eq_vit[corresp[num_face]].set_u_y0(0,0.);
495 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
501 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
503 for(
int i=0 ; i<
nb_pts+1; i++)
504 eq_vit[compteur_faces_paroi].set_F(0, i, -gradient_de_pression0+ts0+ts_boussi(i));
511 eq_vit[compteur_faces_paroi].aller_au_temps(tps);
516 if(itmax<eq_vit[corresp[num_face]].get_it())
518 itmax = eq_vit[compteur_faces_paroi].get_it();
520 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
555 d_visco = tab_visco[elem];
557 tab_u_star_(num_face) = pow(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0),0.25);
568 else if (y_plus>8 && y_plus<30)
572 double deriv = Fdypar_direct(lm_plus);
573 double x = lm_plus*u_star*deriv;
574 tab_k(elem) = x*x/sqrt(
Cmu_) ;
575 tab_eps(elem) = (tab_k(elem)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
580 tab_k(elem) =us2/sqrt(
Cmu_);
585 compteur_faces_paroi++;
601 for (
int num_face=ndeb; num_face<nfin; num_face++)
604 ori = orientation(num_face);
611 if ((elem = face_voisins(num_face,0)) == -1)
613 elem = face_voisins(num_face,1);
620 eq_vit[corresp[num_face]].set_dt(1.e12);
622 eq_vit[corresp[num_face]].set_dt(dt);
626 int face1 = elem_faces(elem,(ori+1));
627 int face2 = elem_faces(elem,(ori+4)%6);
629 int face3 = elem_faces(elem,(ori+2));
630 int face4 = elem_faces(elem,(ori+5)%6);
637 vmoy = 0.5*(vit(face1) + vit(face2));
641 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
642 + termes_sources(face2)/volumes_entrelaces(face2));
644 eq_vit[corresp[num_face]].set_u_y0(0,0.);
646 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
652 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
658 vmoy = 0.5*(vit(face3) + vit(face4));
660 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
661 + termes_sources(face4)/volumes_entrelaces(face4));
663 eq_vit[corresp[num_face]].set_u_y0(1,0.);
665 eq_vit[corresp[num_face]].set_u_yn(1,vmoy);
673 gradient_de_pression1 = 0.5*(grad_p(face3)/volumes_entrelaces(face3)+grad_p(face4)/volumes_entrelaces(face4));
676 DoubleTab ts_boussi(
nb_pts+1,2);
681 const DoubleVect& gravite = ch_gravite.
valeurs();
693 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
700 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
703 for(
int i=0; i<
nb_pts+1; i++)
714 calculer_convection(num_face,face1, face2, face3, face4, elem, ndeb, nfin, ori, gradient_de_pression0, ts0, gradient_de_pression1, ts1);
717 for(
int i=0 ; i<
nb_pts+1 ; i++)
719 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 + ts_boussi(i,0));
720 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 + ts_boussi(i,1));
736 eq_vit[corresp[num_face]].aller_au_temps(tps);
740 if(itmax<eq_vit[corresp[num_face]].get_it())
742 itmax = eq_vit[corresp[num_face]].get_it();
744 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
783 d_visco = tab_visco[elem];
785 tab_u_star_(num_face) = sqrt(sqrt(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0)
786 +eq_vit[corresp[num_face]].get_cis(1)*eq_vit[corresp[num_face]].get_cis(1)));
798 else if (y_plus>8 && y_plus<30)
802 double deriv = Fdypar_direct(lm_plus);
803 double x = lm_plus*u_star*deriv;
804 tab_k(elem) = x*x/sqrt(
Cmu_);
805 tab_eps(elem) = (tab_k(elem)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
810 tab_k(elem)=us2/sqrt(
Cmu_);
814 compteur_faces_paroi++;
824 fic << tps <<
" " << itmax << finl;
848 const IntVect& orientation = domaine_VDF.
orientation();
849 const IntTab& face_voisins = domaine_VDF.
face_voisins();
850 const IntTab& elem_faces = domaine_VDF.
elem_faces();
852 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
856 const DoubleTab& tab_visco = ch_visco_cin.
valeurs();
864 Cerr <<
" On ne sait pas traiter la gravite non uniforme dans TBLE !!!"<< finl;
873 visco = std::max(tab_visco(0,0),DMINFLOAT);
883 Cerr <<
"In ParoiVDF_TBLE::calculer_hyd : visco = " << tab_visco.
local_min_vect() <<
" <= 0 ? " << finl;
895 double gradient_de_pression0 = 0., gradient_de_pression1 = 0., vmoy = 0., ts0 =0., ts1=0.;
900 DoubleTab grad_p(vit);
903 gradient.calculer(p, grad_p);
904 const DoubleTab& visco_turb = mon_modele_turb_hyd->viscosite_turbulente().valeurs();
907 DoubleTab termes_sources;
912 int nb_sources=les_sources.size();
913 for (
int j=0; j<nb_sources; j++)
927 int compteur_faces_paroi = 0;
932 compteur_faces_paroi = 0;
934 for (
int n_bord=0; n_bord<domaine_VDF.
nb_front_Cl(); n_bord++)
941 const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
961 for (
int num_face=ndeb; num_face<nfin; num_face++)
963 ori = orientation(num_face);
968 if ((elem = face_voisins(num_face,0)) == -1)
970 elem = face_voisins(num_face,1);
977 eq_vit[corresp[num_face]].set_dt(1.e12);
979 eq_vit[corresp[num_face]].set_dt(dt);
984 int face1 = elem_faces(elem,(ori+1));
985 int face2 = elem_faces(elem,(ori+3)%4);
992 DoubleVect ts_boussi(
nb_pts+1);
997 const DoubleVect& gravite = ch_gravite.
valeurs();
1009 for (
int idim=0; idim<
dimension; idim++)
1011 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1016 for (
int idim=0; idim<
dimension; idim++)
1018 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1021 double g_t=gt_vect(1-ori);
1023 for(
int i=0; i<
nb_pts+1; i++)
1037 vmoy = 0.5*(vit(face1) + vit(face2));
1041 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1042 + termes_sources(face2)/volumes_entrelaces(face2));
1044 eq_vit[corresp[num_face]].set_u_y0(0,0.);
1046 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
1052 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
1054 for(
int i=0 ; i<
nb_pts+1; i++)
1055 eq_vit[compteur_faces_paroi].set_F(0, i, -gradient_de_pression0+ts0+ts_boussi(i));
1062 eq_vit[compteur_faces_paroi].aller_au_temps(tps);
1067 if(itmax<eq_vit[corresp[num_face]].get_it())
1069 itmax = eq_vit[compteur_faces_paroi].get_it();
1071 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1106 d_visco = tab_visco[elem];
1108 tab_u_star_(num_face) = pow(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0),0.25);
1118 tab1(elem)=eq_vit[corresp[num_face]].get_nu_t_yn();
1129 else if (y_plus>8 && y_plus<30)
1133 double deriv = Fdypar_direct(lm_plus);
1134 double x = lm_plus*u_star*deriv;
1135 tab1(elem,0) = x*x/sqrt(
Cmu_) ;
1136 tab1(elem,1) = (tab1(elem,0)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
1141 tab1(elem,0)=us2/sqrt(
Cmu_);
1146 compteur_faces_paroi++;
1162 for (
int num_face=ndeb; num_face<nfin; num_face++)
1165 ori = orientation(num_face);
1172 if ((elem = face_voisins(num_face,0)) == -1)
1174 elem = face_voisins(num_face,1);
1181 eq_vit[corresp[num_face]].set_dt(1.e12);
1183 eq_vit[corresp[num_face]].set_dt(dt);
1187 int face1 = elem_faces(elem,(ori+1));
1188 int face2 = elem_faces(elem,(ori+4)%6);
1190 int face3 = elem_faces(elem,(ori+2));
1191 int face4 = elem_faces(elem,(ori+5)%6);
1198 vmoy = 0.5*(vit(face1) + vit(face2));
1202 ts0 = 0.5*(termes_sources(face1)/volumes_entrelaces(face1)
1203 + termes_sources(face2)/volumes_entrelaces(face2));
1205 eq_vit[corresp[num_face]].set_u_y0(0,0.);
1207 eq_vit[corresp[num_face]].set_u_yn(0,vmoy);
1213 gradient_de_pression0 = 0.5*(grad_p(face1)/volumes_entrelaces(face1)+grad_p(face2)/volumes_entrelaces(face2));
1219 vmoy = 0.5*(vit(face3) + vit(face4));
1221 ts1 = 0.5*(termes_sources(face3)/volumes_entrelaces(face3)
1222 + termes_sources(face4)/volumes_entrelaces(face4));
1224 eq_vit[corresp[num_face]].set_u_y0(1,0.);
1226 eq_vit[corresp[num_face]].set_u_yn(1,vmoy);
1234 gradient_de_pression1 = 0.5*(grad_p(face3)/volumes_entrelaces(face3)+grad_p(face4)/volumes_entrelaces(face4));
1237 DoubleTab ts_boussi(
nb_pts+1,2);
1242 const DoubleVect& gravite = ch_gravite.
valeurs();
1252 for (
int idim=0; idim<
dimension; idim++)
1254 gn+=gravite(idim)*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1259 for (
int idim=0; idim<
dimension; idim++)
1261 gt_vect(idim) = gravite(idim)-gn*domaine_VDF.
face_normales(num_face,idim)/norm_n;
1264 for(
int i=0; i<
nb_pts+1; i++)
1281 calculer_convection(num_face,face1, face2, face3, face4, elem, ndeb, nfin, ori, gradient_de_pression0, ts0, gradient_de_pression1, ts1);
1284 for(
int i=0 ; i<
nb_pts+1 ; i++)
1286 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 + ts_boussi(i,0));
1287 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 + ts_boussi(i,1));
1303 eq_vit[corresp[num_face]].aller_au_temps(tps);
1307 if(itmax<eq_vit[corresp[num_face]].get_it())
1309 itmax = eq_vit[corresp[num_face]].get_it();
1311 Cerr <<
"WARNING : TOO MANY ITERATIONS ARE NEEDED IN THE TBLE MESH"
1350 d_visco = tab_visco[elem];
1352 tab_u_star_(num_face) = sqrt(sqrt(eq_vit[corresp[num_face]].get_cis(0)*eq_vit[corresp[num_face]].get_cis(0)
1353 +eq_vit[corresp[num_face]].get_cis(1)*eq_vit[corresp[num_face]].get_cis(1)));
1361 tab1(elem)=eq_vit[corresp[num_face]].get_nu_t_yn();
1375 else if (y_plus>8 && y_plus<30)
1379 double deriv = Fdypar_direct(lm_plus);
1380 double x = lm_plus*u_star*deriv;
1381 tab1(elem,0) = x*x/sqrt(
Cmu_);
1382 tab1(elem,1) = (tab1(elem,0)*u_star*u_star*deriv)*sqrt(
Cmu_)/d_visco;
1387 tab1(elem,0)=us2/sqrt(
Cmu_);
1392 compteur_faces_paroi++;
1401 SFichier fic(
"iter.dat", ios::app);
1402 fic << tps <<
" " << itmax << finl;
1417int ParoiVDF_TBLE::calculer_stats()
1420 const IntVect& orientation = domaine_VDF.
orientation();
1422 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1436 for(
int i=0 ; i<
nb_pts+1 ; i++)
1451 for(
int i=0 ; i<
nb_pts+1 ; i++)
1463 for(
int i=0 ; i<
nb_pts+1 ; i++)
1475 for(
int i=0 ; i<
nb_pts+1 ; i++)
1494 const Equation_base& eqn_hydr = mon_modele_turb_hyd->equation();
1498 const IntVect& orientation = domaine_VDF.
orientation();
1513 for(
int i=0; i<
nb_pts+1; i++)
1519 for(
int i=0; i<
nb_pts+1; i++)
1523 for(
int i=0; i<
nb_pts+1; i++)
1527 for(
int i=0; i<
nb_pts+1; i++)
1542 double tps = mon_modele_turb_hyd->equation().inconnue().temps();
1552 double tps_reprise = mon_modele_turb_hyd->equation().schema_temps().temps_courant();
1564void ParoiVDF_TBLE::calculer_convection(
int num_face,
int face1,
int face2,
int face3,
int face4,
int elem,
int ndeb,
int nfin,
int ori,
double gradient_de_pression0,
double ts0,
double gradient_de_pression1,
double ts1)
1567 const IntTab& face_voisins = domaine_VDF.
face_voisins();
1568 const IntTab& elem_faces = domaine_VDF.
elem_faces();
1571 double d00d0 = 0., d0vdy= 0., d01d1= 0., d10d0= 0., d1vdy= 0., d11d1=0.;
1572 int elem_gauche0, elem_gauche1, elem_droite0, elem_droite1;
1573 int face_gauche0, face_gauche1, face_droite0, face_droite1;
1582 elem_gauche0 = face_voisins(face1,0);
1583 elem_droite0 = face_voisins(face2,1);
1584 if (elem_gauche0 == elem)
1586 elem_droite0 = face_voisins(face1,1);
1587 elem_gauche0 = face_voisins(face2,0);
1590 elem_gauche1 = face_voisins(face3,0);
1591 elem_droite1 = face_voisins(face4,1);
1592 if (elem_gauche1 == elem)
1594 elem_droite1 = face_voisins(face3,1);
1595 elem_gauche1 = face_voisins(face4,0);
1599 if((elem_droite1!=-1)&&(elem_droite0!=-1)&&(elem_gauche1!=-1)&&(elem_gauche0!=-1))
1602 face_gauche0 = elem_faces(elem_gauche0, ori);
1603 if ((face_voisins(face_gauche0,0) != -1)&&(face_voisins(face_gauche0,1) != -1))
1605 face_gauche0 = elem_faces(elem_gauche0, ori+3);
1607 face_droite0 = elem_faces(elem_droite0, ori);
1608 if ((face_voisins(face_droite0,0) != -1)&&(face_voisins(face_droite0,1) != -1))
1610 face_droite0 = elem_faces(elem_droite0, ori+3);
1613 face_gauche1 = elem_faces(elem_gauche1, ori);
1614 if ((face_voisins(face_gauche1,0) != -1)&&(face_voisins(face_gauche1,1) != -1))
1616 face_gauche1 = elem_faces(elem_gauche1, ori+3);
1618 face_droite1 = elem_faces(elem_droite1, ori);
1619 if ((face_voisins(face_droite1,0) != -1)&&(face_voisins(face_droite1,1) != -1))
1621 face_droite1 = elem_faces(elem_droite1, ori+3);
1625 if((face_droite1<nfin)&&(face_droite0<nfin)&&(face_gauche1<nfin)&&(face_gauche0<nfin)&&
1626 (ndeb<=face_droite1)&&(ndeb<=face_droite0)&&(ndeb<=face_gauche1)&&(ndeb<=face_gauche0))
1635 if(eq_vit[corresp[num_face]].get_Unp1(0,1)>0)
1636 d0d0_1 = (eq_vit[corresp[face_droite0]].get_Unp1(0,1)-eq_vit[corresp[num_face]].get_Unp1(0,1))
1639 d0d0_1 = (eq_vit[corresp[num_face]].get_Unp1(0,1)-eq_vit[corresp[face_gauche0]].get_Unp1(0,1))
1642 if(eq_vit[corresp[num_face]].get_Unp1(0,1)>0)
1643 d1d1_1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,1)-eq_vit[corresp[num_face]].get_Unp1(1,1))
1646 d1d1_1 = (eq_vit[corresp[num_face]].get_Unp1(1,1)-eq_vit[corresp[face_gauche1]].get_Unp1(1,1))
1649 v += -(d0d0_1+d1d1_1)
1650 *0.5*(eq_vit[corresp[num_face]].get_yc(1)-eq_vit[corresp[num_face]].get_yc(0));
1652 eq_vit[corresp[num_face]].set_v(i,v);
1658 for(i=2 ; i<
nb_pts ; i++)
1661 (eq_vit[corresp[face_droite0]].get_Unp1(0,i)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1664 (eq_vit[corresp[face_droite1]].get_Unp1(1,i)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1667 *0.5*(eq_vit[corresp[num_face]].get_yc(i)-eq_vit[corresp[num_face]].get_yc(i-1))
1669 (eq_vit[corresp[face_droite0]].get_Unp1(0,i-1)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i-1))
1672 (eq_vit[corresp[face_droite1]].get_Unp1(1,i-1)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i-1))
1675 *0.5*(eq_vit[corresp[num_face]].get_yc(i-1)-eq_vit[corresp[num_face]].get_yc(i-2));
1677 eq_vit[corresp[num_face]].set_v(i,v);
1688 (eq_vit[corresp[face_droite0]].get_Unp1(0,i)-eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1691 (eq_vit[corresp[face_droite1]].get_Unp1(1,i)-eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1694 *0.5*(eq_vit[corresp[num_face]].get_yc(
nb_pts-1)-eq_vit[corresp[num_face]].get_yc(
nb_pts-2));
1696 eq_vit[corresp[num_face]].set_v(
nb_pts,v);
1707 eq_vit[corresp[num_face]].set_F(0, 0, - gradient_de_pression0 + ts0);
1708 eq_vit[corresp[num_face]].set_F(1, 0, - gradient_de_pression1 + ts1);
1710 for(i=1 ; i<
nb_pts ; i++)
1712 d00d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(0,i)
1713 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1717 d01d1 = (eq_vit[corresp[face_droite1]].get_Unp1(0,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1718 -eq_vit[corresp[face_gauche1]].get_Unp1(0,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1722 d10d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(1,i)
1723 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(1,i))
1727 d11d1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1728 -eq_vit[corresp[face_gauche1]].get_Unp1(1,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1732 d0vdy = (eq_vit[corresp[num_face]].get_Unp1(0,i+1)*eq_vit[corresp[num_face]].get_v(i+1)
1733 -eq_vit[corresp[num_face]].get_Unp1(0,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1734 /(eq_vit[corresp[num_face]].get_y(i+1)-eq_vit[corresp[num_face]].get_y(i-1));
1736 d1vdy = (eq_vit[corresp[num_face]].get_Unp1(1,i+1)*eq_vit[corresp[num_face]].get_v(i+1)
1737 -eq_vit[corresp[num_face]].get_Unp1(1,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1738 /(eq_vit[corresp[num_face]].get_y(i+1)-eq_vit[corresp[num_face]].get_y(i-1));
1740 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 - (d00d0+d01d1+d0vdy));
1741 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 - (d11d1+d10d0+d1vdy));
1747 d00d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(0,i)
1748 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(0,i))
1752 d01d1 = (eq_vit[corresp[face_droite1]].get_Unp1(0,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1753 -eq_vit[corresp[face_gauche1]].get_Unp1(0,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1757 d10d0 = (eq_vit[corresp[face_droite0]].get_Unp1(0,i)*eq_vit[corresp[face_droite0]].get_Unp1(1,i)
1758 -eq_vit[corresp[face_gauche0]].get_Unp1(0,i)*eq_vit[corresp[face_gauche0]].get_Unp1(1,i))
1762 d11d1 = (eq_vit[corresp[face_droite1]].get_Unp1(1,i)*eq_vit[corresp[face_droite1]].get_Unp1(1,i)
1763 -eq_vit[corresp[face_gauche1]].get_Unp1(1,i)*eq_vit[corresp[face_gauche1]].get_Unp1(1,i))
1767 d0vdy = (eq_vit[corresp[num_face]].get_Unp1(0,i)*eq_vit[corresp[num_face]].get_v(i)
1768 -eq_vit[corresp[num_face]].get_Unp1(0,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1769 /(eq_vit[corresp[num_face]].get_y(i)-eq_vit[corresp[num_face]].get_y(i-1));
1771 d1vdy = (eq_vit[corresp[num_face]].get_Unp1(1,i)*eq_vit[corresp[num_face]].get_v(i)
1772 -eq_vit[corresp[num_face]].get_Unp1(1,i-1)*eq_vit[corresp[num_face]].get_v(i-1))
1773 /(eq_vit[corresp[num_face]].get_y(i)-eq_vit[corresp[num_face]].get_y(i-1));
1779 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0 - (d00d0+d01d1+d0vdy));
1780 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1 - (d11d1+d10d0+d1vdy));
1788 for(
int i=0 ; i<
nb_pts+1 ; i++)
1790 eq_vit[corresp[num_face]].set_F(0, i, - gradient_de_pression0 + ts0);
1791 eq_vit[corresp[num_face]].set_F(1, i, - gradient_de_pression1 + ts1);
classe Champ_Don_base classe de base des Champs donnes (non calcules)
DoubleTab & valeurs() override
Surcharge Champ_base::valeurs() Renvoie le tableau des valeurs.
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
double temps() const
Renvoie le temps du champ.
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Classe Diffu_totale_hyd_base Classe abstraite calculant la diffusivite totale (somme diffusivite.
classe Dirichlet_paroi_fixe Represente une paroi immobile dans une equation de type Navier_Stokes.
double dist_norm_bord_axi(int num_face) const
int orientation(int) const override
inline DoubleVect& Domaine_VDF::porosite_face() {
double dist_elem_period(int, int, int) const
double face_normales(int, int) const override
double dist_norm_bord(int num_face) const override
virtual const DoubleVect & face_surfaces() const
int nb_faces() const
renvoie le nombre global de faces.
DoubleVect & volumes_entrelaces()
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
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Sortie & syncfile() override
Provoque l'ecriture sur disque des donnees accumulees sur les differents processeurs depuis le dernie...
Class defining operators and methods for all reading operation in an input flow (file,...
classe Equation_base Le role d'une equation est le calcul d'un ou plusieurs champs....
virtual const Milieu_base & milieu() const =0
virtual const RefObjU & get_modele(Type_modele type) const
Sources & sources()
Renvoie les termes sources asssocies a l'equation.
virtual const Champ_Inc_base & inconnue() const =0
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
classe Fluide_base Cette classe represente un d'un fluide incompressible ainsi que
const Champ_Don_base & viscosite_cinematique() const
int num_premiere_face() const
virtual const Champ_Don_base & beta_t() const
Renvoie beta_t du milieu.
virtual const Champ_Don_base & gravite() const
Renvoie la gravite du milieu si elle a ete associe provoque une erreur sinon.
Classe Modele_turbulence_scal_base Cette classe represente un modele de turbulence pour une equation ...
const Turbulence_paroi_scal_base & loi_paroi() const
Renvoie la loi de turbulence sur la paroi (version const).
Une chaine de caractere (Nom) en majuscules.
classe Navier_Stokes_std Cette classe porte les termes de l'equation de la dynamique
Operateur_Grad & operateur_gradient()
Renvoie l'operateur de calcul du gradient associe a l'equation.
Champ_Inc_base & pression()
class Nom Une chaine de caractere pour nommer les objets de TRUST
virtual void set_param(Param &) const
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
virtual const Nom & le_nom() const
Donne le nom de l'Objet_U Methode a surcharger : renvoie "neant" dans cette implementation.
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Classe Operateur_Grad Classe generique de la hierarchie des operateurs calculant le gradient.
Helper class to factorize the readOn method of Objet_U classes.
void ajouter_flag(const char *keyword, const bool *value)
Register a boolean flag whose mere presence switches it to true.
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.
void set_param(Param ¶m) const override
int reprendre(Entree &) override
Reprise d'un Objet_U sur un flot d'entree Methode a surcharger.
int calculer_hyd_BiK(DoubleTab &, DoubleTab &) override
void imprimer_ustar(Sortie &) const override
int lire_motcle_non_standard(const Motcle &, Entree &) override
Lecture des parametres de type non simple d'un objet_U a partir d'un flot d'entree.
int init_lois_paroi() override
int sauvegarder(Sortie &) const override
Sauvegarde d'un Objet_U sur un flot de sortie Methode a surcharger.
int calculer_hyd(DoubleTab &) override
const Probleme_base & getPbBase() const override
int est_initialise() const
Eq_couch_lim & get_eq_couche_lim_T(int i)
double tps_start_stat_nu_t_dyn
void calculer_stat(int indice_post, int indice_maillage, double Fx, double Fy, double Fz, double u, double v, double w, double dt)
void imprimer_stat(Sortie &, double) const
void set_param(Param ¶m) const
int reprendre(Entree &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps)
IntTab num_global_faces_post
int sauvegarder(Sortie &, const Domaine_dis_base &, const Domaine_Cl_dis_base &, double tps) const
DoubleTab valeurs_reprises
int lire_motcle_non_standard(const Motcle &, Entree &)
int init_lois_paroi(const Domaine_VF &, const Domaine_Cl_dis_base &)
CLASS: Paroi_ODVM_scal_VDF.
void set_param(Param ¶m) const
double calcul_lm_plus(double d_plus)
classe Probleme_base C'est un Probleme_U qui n'est pas un couplage.
virtual const Equation_base & equation(int) const =0
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
static int je_suis_maitre()
renvoie 1 si on est sur le processeur maitre du groupe courant (c'est a dire me() == 0),...
Cette classe est a la classe C++ ofstream ce que la classe Sortie est a la classe C++ ostream Elle re...
double temps_courant() const
Renvoie le temps courant.
double pas_temps_min() const
Renvoie le pas de temps minimum.
double pas_de_temps() const
Renvoie le pas de temps (delta_t) courant.
Classe de base des flux de sortie.
classe Source_base Un objet Source_base est un terme apparaissant au second membre d'une
virtual DoubleTab & ajouter(DoubleTab &) const
class Sources Sources represente une liste de Source.
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
_TYPE_ local_min_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
const Objet_U & valeur() const
class Terme_Boussinesq_scalaire_VDF_Face
Classe Terme_Boussinesq_base Cette classe represente le terme de gravite qui figure dans l'equation.
double Scalaire0(int i) const
const DoubleVect & tab_d_plus() const
const DoubleVect & tab_u_star() const
DoubleTab Cisaillement_paroi_
Classe Turbulence_paroi_scal_base Classe de base pour la hierarchie des classes representant les mode...