TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Traitement_particulier_NS_Brech_VEF.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 <Traitement_particulier_NS_Brech_VEF.h>
17#include <Navier_Stokes_std.h>
18#include <Fluide_Incompressible.h>
19#include <Probleme_base.h>
20#include <TRUSTTrav.h>
21#include <Periodique.h>
22#include <Champ_P1NC.h>
23
24Implemente_instanciable_sans_constructeur(Traitement_particulier_NS_Brech_VEF,"Traitement_particulier_NS_Brech_VEF",Traitement_particulier_NS_VEF);
25
29
30Sortie& Traitement_particulier_NS_Brech_VEF::printOn(Sortie& is) const
31{
32 return is;
33}
34
35Entree& Traitement_particulier_NS_Brech_VEF::readOn(Entree& is)
36{
37 return is;
38}
39
40/*
41 inline int Traitement_particulier_NS_Brech_VEF::comprend_champ(const Motcle& mot) const
42 {
43 if ( (mot == "Richardson") || (mot == "Pression_porosite") ) return 1 ;
44 return 0 ;
45 }
46
47 int Traitement_particulier_NS_Brech_VEF::a_pour_Champ_Fonc(const Motcle& mot,
48 OBS_PTR(Champ_base)& ch_ref) const
49 {
50 if (mot == "Richardson")
51 {
52 ch_ref = ch_ri;
53 return 1 ;
54 } else if (mot == "Pression_porosite")
55 {
56 ch_ref = ch_p;
57 return 1 ;
58 }
59 return 0 ;
60 }
61*/
62Entree& Traitement_particulier_NS_Brech_VEF::lire(Entree& is)
63{
64 Motcle accouverte = "{" , accfermee = "}" ;
65 Motcle motbidon, motlu;
66
67 is >> motbidon ;
68
69 Motcles les_mots(3);
70 {
71 les_mots[0] = "calcul_flux";
72 les_mots[1] = "Richardson";
73 les_mots[2] = "Pression_porosite";
74 }
75
76
77 if (motbidon == accouverte)
78 {
79 // is >> motlu;
80
81 while(motlu != accfermee)
82 {
83 is >> motlu;
84 int rang=les_mots.search(motlu);
85 switch(rang)
86 {
87 case 0 :
88 {
89 is >> motbidon ;
90 if (motbidon != accouverte)
91 {
92 Cerr << " expected { instead of " << motbidon << finl ;
93 }
94
95 is >> nb ;
96
97 C_trans.resize(nb,3) ;
98 R_loc.resize(nb,3) ;
99 r_int.resize(nb) ;
100 r_out.resize(nb) ;
101 delta_r.resize(nb) ;
102 delta_teta.resize(nb);
103
104 for (int i=0; i<nb; i++)
105 {
106
107 Cerr << " Reading cross-section plane in x " << finl;
108
109 is >> R_loc(i,0) ;
110
111 is >> R_loc(i,1) ;
112
113 is >> R_loc(i,2) ;
114
115 is >> C_trans(i,0) ;
116
117 is >> C_trans(i,1) ;
118
119 is >> C_trans(i,2) ;
120
121 is >> r_int(i);
122
123 is >> r_out(i);
124
125 is >> delta_r(i);
126
127 is >> delta_teta(i) ;
128 }
129
130 c_flux = 1 ;
131
132 is >> motbidon ;
133 if (motbidon != accfermee)
134 {
135 Cerr << " expected } instead of " << motbidon << finl ;
136 }
137 break ;
138
139 }
140 case 1 :
141 {
142 Cerr << " Reading Richardson " << finl;
143 const Domaine_dis_base& zdis=mon_equation->inconnue().domaine_dis_base();
144 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, zdis);
145 // const Probleme_base& pb = mon_equation->probleme();
146 const int nb_faces = domaine_VEF.nb_faces() ;
147
148 ch_ri.associer_domaine_dis_base(zdis);
149 ch_ri.nommer("Richardson");
150 ch_ri.fixer_nb_comp(1);
151 ch_ri.fixer_nb_valeurs_nodales(nb_faces);
152 ch_ri.fixer_unite("-");
153 c_Richardson = 1 ;
154 champs_compris_.ajoute_champ(ch_ri);
155 break ;
156
157 }
158 case 2 :
159 {
160 Cerr << " Reading Pression_porosite " << finl;
161 const Domaine_dis_base& zdis=mon_equation->inconnue().domaine_dis_base();
162 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, zdis);
163 // const Probleme_base& pb = mon_equation->probleme();
164 const int nb_elem = domaine_VEF.nb_elem() ;
165
166 ch_p.associer_domaine_dis_base(zdis);
167 ch_p.nommer("Pression_porosite");
168 ch_p.fixer_nb_comp(1);
169 ch_p.fixer_nb_valeurs_nodales(nb_elem);
170 ch_p.fixer_unite("Pa.m3/kg");
171 c_pression = 1 ;
172 champs_compris_.ajoute_champ(ch_p);
173 break ;
174 }
175 default :
176 {
177 if (motlu == accfermee)
178 {
179 break ;
180 }
181 else
182 {
183 Cerr << "Error in reading Traitement_particulier_Brech_VEF" << finl;
184 Cerr << "Possible keywords are: calcul_flux or Richardson " << finl;
185 Cerr << "Read keyword: " << motlu << finl;
186 exit();
187 break;
188 }
189 }
190 }
191 }
192 is >> motlu;
193 if (motlu != accfermee)
194 {
195 Cerr << "Error reading Traitement_particulier_NS_Brech_VEF 1 ";
196 Cerr << "Expected a } but got " << motlu << finl;
197 exit();
198 }
199 // is >> motlu;
200 // if (motlu != accfermee)
201 // {
202 // Cerr << "Error reading Traitement_particulier_NS_Brech_VEF 2 ";
203 // Cerr << "Expected a } but got " << motlu << finl;
204 // exit();
205 // }
206 // is >> motlu;
207 // if (motlu != accfermee)
208 // {
209 // Cerr << "Error reading Traitement_particulier_NS_Brech_VEF 3 ";
210 // Cerr << "Expected a } but got " << motlu << finl;
211 // exit();
212 // }
213 }
214
215 return is;
216}
217
218
219void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier()
220{
221
222 if (c_flux == 1 ) post_traitement_particulier_calcul_flux() ;
223 if (c_Richardson == 1 ) post_traitement_particulier_Richardson() ;
224 if (c_pression == 1 ) post_traitement_particulier_calcul_pression() ;
225}
226
227void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier_calcul_flux()
228{
229
230 for (int ii =0; ii<nb; ii++)
231 {
232 int taille = 1 ;
233 DoubleTab coord_trace (taille,3) ;
234 DoubleTab Surf_trace(taille,3) ;
235 DoubleTab valeurs_(taille,3) ;
236
237 // VB modifications to account for rho and compute enthalpy flux + Tmoy
238 OBS_PTR(Champ_base) rch1 ;
239 OBS_PTR(Champ_Inc_base) l_inco ;
240 const Probleme_base& pb = mon_equation->probleme();
241 /*
242 pb.a_pour_Champ_Inc("Temperature", rch1 ) ;
243 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
244 const Champ_Inc_base& temp = l_inco.valeur();
245 */
246 rch1 = pb.get_champ("Temperature");
247 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
248 const Champ_Inc_base& temp = l_inco.valeur();
249
250 DoubleTab temper_(taille,1);
251 DoubleTab rho_(taille,1) ;
252 DoubleTab cp_(taille,1) ;
253 // End VB modifications
254 //
255 double rad, teta ;
256 double dr = ( r_out(ii) - r_int(ii) ) /delta_r(ii) ;
257 double dteta = 2.*3.14159/delta_teta(ii) ;
258
259 taille = (int)(delta_r(ii)) * (int)(delta_teta(ii)) ;
260 coord_trace.resize(taille,3) ;
261 Surf_trace.resize(taille,3) ;
262 valeurs_.resize(taille,3);
263
264 // VB modifications to account for rho and compute enthalpy flux + Tmoy
265 rho_.resize(taille,1) ;
266 temper_.resize(taille,1);
267 cp_.resize(taille,1);
268 // End VB modifications
269
270 double a,b,c ;
271
272 // VB modifications to compute the flux on an inclined face
273 // anglex is the angle between the face normal and the x axis (acute angle)
274 double anglex = acos(C_trans(ii,0)/sqrt( C_trans(ii,0) * C_trans(ii,0)
275 + C_trans(ii,1) * C_trans(ii,1)
276 + C_trans(ii,2) * C_trans(ii,2) )) ;
277 // End VB modifications
278
279 int i =-1 ;
280
281 for ( rad = (r_int(ii)+(dr/2.)) ; rad < r_out(ii) ; rad += dr )
282 for (teta = dteta/2. ; teta < 2.*3.14159 ; teta += dteta )
283 {
284 i++ ;
285 // VB modifications to compute the flux on an inclined face
286 // a = 0. ;
287 // b = rad * cos(teta) ;
288 a = - rad * cos(teta) * sin(anglex) ;
289 b = rad * cos(teta) * cos(anglex) ;
290 // End VB modifications
291 c = rad * sin(teta) ;
292 coord_trace(i,0) = a + R_loc(ii,0) ;
293 coord_trace(i,1) = b + R_loc(ii,1) ;
294 coord_trace(i,2) = c + R_loc(ii,2) ;
295 double ra = (dr/2.+rad) * (dr/2.+rad) ;
296 double ri = (rad - (dr/2.))* (rad - (dr/2.)) ;
297 Surf_trace(i,0) = 3.14159 * (ra -ri) / delta_teta(ii) ;
298 }
299
300 double flux_pos = 0. ;
301 double flux_neg = 0. ;
302 double flux ;
303 // VB modifications to account for rho and compute enthalpy flux + Tmoy
304 double tempmoy = 0. ;
305 double massvol = 0. ;
306 double fluxE_pos = 0. ;
307 double fluxE_neg = 0. ;
308 double fluxE ;
309 // end VB modifications
310
311 const Domaine_dis_base& zdis=mon_equation->inconnue().domaine_dis_base();
312 const Domaine& domaine=zdis.domaine();
313 IntVect les_polys(coord_trace.dimension(0));
314 domaine.chercher_elements(coord_trace, les_polys);
315 mon_equation->inconnue().valeur_aux_elems(coord_trace, les_polys, valeurs_);
316
317 // VB modifications to account for rho and compute enthalpy flux + Tmoy
318 mon_equation->fluide().masse_volumique().valeur_aux_elems(coord_trace, les_polys, rho_);
319 mon_equation->fluide().capacite_calorifique().valeur_aux_elems(coord_trace, les_polys, cp_);
320 temp.valeur_aux_elems(coord_trace, les_polys, temper_);
321 // End VB modifications
322
323 // Zero out the virtual part
324 for (i=0 ; i<taille; i++)
325 if (les_polys(i)>=domaine.nb_elem())
326 {
327 for (int k=0; k<dimension; k++)
328 valeurs_(i,k)=0;
329 // VB modifications to account for rho and compute enthalpy flux + Tmoy
330 rho_(i)=0.;
331 temper_(i)=0.;
332 cp_(i)=0.;
333 }
334 // End VB modifications
335
336 for (i=0 ; i<taille; i++ )
337 {
338 // VB modifications to account for rho and compute enthalpy flux + Tmoy
339 // flux = valeurs_(i,0)*Surf_trace(i,0) ;
340 massvol += rho_(i) ;
341 tempmoy += temper_(i) ;
342 flux = ( valeurs_(i,0)*cos(anglex) + valeurs_(i,1)*sin(anglex) )
343 * Surf_trace(i,0)*rho_(i) ;
344 fluxE = ( valeurs_(i,0)*cos(anglex) + valeurs_(i,1)*sin(anglex) )
345 * Surf_trace(i,0)*rho_(i)*cp_(i)*temper_(i) ;
346 flux=Process::mp_sum(flux);
347 fluxE=Process::mp_sum(fluxE);
348 if ( flux > 0. ) flux_pos += flux ;
349 if ( flux < 0. ) flux_neg += flux ;
350 if ( fluxE > 0. ) fluxE_pos += fluxE ;
351 if ( fluxE < 0. ) fluxE_neg += fluxE ;
352 }
353
354 tempmoy=Process::mp_sum(tempmoy);
355 massvol=Process::mp_sum(massvol);
356 tempmoy /= taille ;
357 massvol /= taille ;
358
359 if (je_suis_maitre())
360 {
361 Cout << finl ;
362 Cout << " ----------------------- " << finl ;
363 Cout << " Traitement_particulier flux in x (Plane Oyz) at: "<< R_loc(ii,0) << " "
364 << R_loc(ii,1) << " "
365 << R_loc(ii,2) << finl ;
366 Cout << " angle of the plane normal with respect to x = " << anglex/3.1415*180. << " (deg) " << finl ;
367 Cout << " flux_pos = " << flux_pos << " flux_neg = " << flux_neg << " (kg/s) " << finl ;
368 Cout << " fluxEnthalpique_pos = " << fluxE_pos << " fluxEnthalpique_neg = " << fluxE_neg << " (Watts) " << finl ;
369 Cout << " average mass density on the cross-section plane: " << massvol << finl ;
370 Cout << " average temperature on the cross-section plane: " << tempmoy << finl ;
371 }
372 // End VB modifications
373 }
374}
375
376void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier_Richardson()
377{
378 const Domaine_dis_base& zdis=mon_equation->domaine_dis();
379 const Domaine_VEF& domaine_VEF=ref_cast(Domaine_VEF, zdis);
380 const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,mon_equation->domaine_Cl_dis() );
381
382 OBS_PTR(Champ_base) rch1 ;
383 OBS_PTR(Champ_Inc_base) l_inco ;
384
385 const Probleme_base& pb = mon_equation->probleme();
386 /*
387 pb.a_pour_Champ_Inc("Temperature", rch1 ) ;
388
389 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
390 */
391 rch1 = pb.get_champ("Temperature");
392 l_inco = ref_cast(Champ_Inc_base, rch1.valeur()) ;
393
394 const Champ_Inc_base& temp = l_inco.valeur();
395 const DoubleTab& temperature = temp.valeurs();
396
397 const DoubleTab& vitesse = mon_equation->inconnue().valeurs();
398 const DoubleVect& gravite = mon_equation->fluide().gravite().valeurs();
399 const DoubleVect& beta = mon_equation->fluide().beta_t().valeurs();
400
401 const int nb_faces = domaine_VEF.nb_faces();
402
403 DoubleVect P(nb_faces) ;
404 DoubleVect G(nb_faces) ;
405
406 DoubleVect beta_i(nb_faces) ;
407 beta_i = beta(0) ;
408
409 calculer_terme_production_K( domaine_VEF, domaine_Cl_VEF,
410 P, vitesse);
411
412 calculer_terme_destruction_K( domaine_VEF, domaine_Cl_VEF,
413 G, temperature,
414 beta_i, gravite) ;
415
416 DoubleVect& richard_loc = ch_ri.valeurs();
417 richard_loc = 0. ;
418
419 for ( int i=0; i<nb_faces; i++ )
420 {
421 if (std::fabs (P[i]) > 1.e-9 )richard_loc [i] = -G[i]/P[i] ;
422 else richard_loc [i] = 0. ;
423 }
424
425 ch_ri.changer_temps(mon_equation->inconnue().temps());
426
427}
428
429void Traitement_particulier_NS_Brech_VEF::post_traitement_particulier_calcul_pression()
430{
431 const Domaine_VEF& zvef=ref_cast(Domaine_VEF, mon_equation->domaine_dis());
432 const DoubleVect& porosite_face = mon_equation->milieu().porosite_face();
433 int i,comp;
434 int nb_face = zvef.nb_faces();
435 Operateur_Div divergence = mon_equation->operateur_divergence();
436 Operateur_Grad gradient = mon_equation->operateur_gradient();
437 SolveurSys solveur_pression_ = mon_equation->solveur_pression();
438
439 DoubleTab& pression=mon_equation->pression().valeurs();
440 DoubleTab& vitesse=mon_equation->vitesse().valeurs();
441 //DoubleTrav gradP(la_vitesse.valeurs());
442 //DoubleTrav gradP(vitesse);
443 //DoubleTrav inc_pre(pression);
444 DoubleTab gradP(vitesse);
445 DoubleTab inc_pre(pression);
446
447 // DoubleVect& inc_pre = ch_p.valeurs();
448 inc_pre = 0. ;
449 //DoubleTrav secmem(pression);
450 DoubleTab secmem(pression);
451
452 gradient.calculer(mon_equation->pression().valeurs(),gradP);
453 //gradient.calculer(la_pression,gradP);
454 // Cerr << "gradP " << gradP << finl;
455 // we want BM-1Bt(psi*Pression)
456 mon_equation->solv_masse().appliquer(gradP);
457
458 DoubleTab grad_temp(vitesse);
459 for(i=0; i<nb_face; i++)
460 {
461 for (comp=0; comp<dimension; comp++)
462 grad_temp(i,comp) = gradP(i,comp)/porosite_face(i);
463 }
464
465 divergence.calculer(grad_temp, secmem);
466 secmem *= -1; // because div = -B
467 solveur_pression_.resoudre_systeme(mon_equation->matrice_pression().valeur(),secmem, inc_pre);
468
469 DoubleVect& la_pression_porosite = ch_p.valeurs();
470 la_pression_porosite = inc_pre ;
471 Cerr << "pressure " << mon_equation->pression().valeurs() << finl;
472 Cerr << "ch_p.valeurs() " << ch_p.valeurs() << finl;
473}
474
475
476////////////////////////////////////////////////////////////////////////////
477//
478// Implementation of member functions of class Calcul_Production_K_VEF
479//
480/////////////////////////////////////////////////////////////////////////////
481
482
483void Traitement_particulier_NS_Brech_VEF::
484calculer_terme_production_K(const Domaine_VEF& domaine_VEF,const Domaine_Cl_VEF& zcl_VEF,
485 DoubleVect& P,
486 const DoubleTab& vit) const
487{
488 // P is discretised like K and Eps, i.e. at element centres
489 //
490 // P(elem) = -(2/3)*k(i)*div_U(i) + nu_t(i) * F(u,v,w)
491 //
492 // 2 2 2
493 // with F(u,v,w) = 2[(du/dx) + (dv/dy) + (dw/dz) ] +
494 //
495 // 2 2 2
496 // (du/dy+dv/dx) + (du/dz+dw/dx) + (dw/dy+dv/dz)
497 //
498 // Note: we work in the incompressible framework so we neglect
499 // the term (2/3)*k(i)*div_U(i)
500
501 P= 0;
502
503 // Compute F(u,v,w):
504
505 // const Domaine& domaine = domaine_VEF.domaine();
506 int nb_elem = domaine_VEF.nb_elem();
507 // const IntTab& elem_faces = domaine_VEF.elem_faces();
508 const IntTab& face_voisins = domaine_VEF.face_voisins();
509 // const int nb_faces = domaine_VEF.nb_faces();
510 const DoubleVect& volumes = domaine_VEF.volumes();
511 // int premiere_face_entier = domaine_VEF.premiere_face_int();
512 // const IntTab& les_Polys = domaine.les_elems();
513
514 int fac=0, poly1, poly2;
515 int nb_faces_ = domaine_VEF.nb_faces();
516
517 // (du/dx du/dy dv/dx dv/dy ...) for a polyhedron
518
519 ///////////////////////////////////////////////////////////////////////////////////////////////
520 // <
521 // compute gradients; < [ Ujp*np/vol(j) ]
522 // j
523 ////////////////////////////////////////////////////////////////////////////////////////////////
524 int n_bord;
525
526 DoubleTab gradient_elem(nb_elem,dimension,dimension);
527 Champ_P1NC::calcul_gradient(vit,gradient_elem,zcl_VEF);
528
529 // We have gradient_elem per element
530 ////////////////////////////////////////////////////////////////////////////////////
531
532 double du_dx;
533 double du_dy;
534 double du_dz;
535 double dv_dx;
536 double dv_dy;
537 double dv_dz;
538 double dw_dx;
539 double dw_dy;
540 double dw_dz;
541
542 /////////////////////////////////////////////////////////////////////////////////////////////////////////
543 // Compute du/dx, dv/dy and the cross-derivatives on the faces of each element in the 2D case
544 /////////////////////////////////////////////////////////////////////////////////////////////////////////
545
546 // Loop over boundaries to process boundary faces,
547 // distinguishing the periodic case
548
549 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
550 {
551 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
552 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
553 int ndeb = le_bord.num_premiere_face();
554 int nfin = ndeb + le_bord.nb_faces();
555
556 if (sub_type(Periodique,la_cl.valeur()))
557 {
558 for (fac=ndeb; fac<nfin; fac++)
559 {
560 poly1 = face_voisins(fac,0);
561 poly2 = face_voisins(fac,1);
562 double a=volumes(poly1)/(volumes(poly1)+volumes(poly2));
563 double b=volumes(poly2)/(volumes(poly1)+volumes(poly2));
564
565 du_dx =a*gradient_elem(poly1,0,0) + b*gradient_elem(poly2,0,0);
566 du_dy =a*gradient_elem(poly1,0,1) + b*gradient_elem(poly2,0,1);
567 dv_dx =a*gradient_elem(poly1,1,0) + b*gradient_elem(poly2,1,0);
568 dv_dy =a*gradient_elem(poly1,1,1) + b*gradient_elem(poly2,1,1);
569
570 // Compute the production term
571
572 P(fac) = ( 2*(du_dx *du_dx + dv_dy *dv_dy) + ((du_dy+dv_dx)*(du_dy+dv_dx)));
573
574 if (dimension == 3)
575 {
576 du_dz=a*gradient_elem(poly1,0,2) + b*gradient_elem(poly2,0,2);
577 dv_dz=a*gradient_elem(poly1,1,2) + b*gradient_elem(poly2,1,2);
578 dw_dx=a*gradient_elem(poly1,2,0) + b*gradient_elem(poly2,2,0);
579 dw_dy=a*gradient_elem(poly1,2,1) + b*gradient_elem(poly2,2,1);
580 dw_dz=a*gradient_elem(poly1,2,2) + b*gradient_elem(poly2,2,2);
581
582 // Compute the production term
583
584 P(fac) = (2*( du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz )
585 + ( (du_dy+dv_dx)*(du_dy+dv_dx)
586 + (du_dz+dw_dx)*(du_dz+dw_dx)
587 + (dw_dy+dv_dz)*(dw_dy+dv_dz)));
588 }
589 }
590 }
591 else
592 {
593 for (fac=ndeb; fac<nfin; fac++)
594 {
595 poly1 = face_voisins(fac,0);
596
597 du_dx=gradient_elem(poly1,0,0);
598 du_dy=gradient_elem(poly1,0,1);
599 dv_dx=gradient_elem(poly1,1,0);
600 dv_dy=gradient_elem(poly1,1,1);
601 // Compute the production term
602 P(fac) = ( 2*(du_dx*du_dx + dv_dy*dv_dy) + ((du_dy+dv_dx)*(du_dy+dv_dx)));
603
604
605 if (dimension == 3)
606 {
607 du_dz=gradient_elem(poly1,0,2);
608 dv_dz=gradient_elem(poly1,1,2);
609 dw_dx=gradient_elem(poly1,2,0);
610 dw_dy=gradient_elem(poly1,2,1);
611 dw_dz=gradient_elem(poly1,2,2);
612
613 // Compute the production term
614
615 P(fac) = (2*( du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz )
616 + ( (du_dy+dv_dx)*(du_dy+dv_dx)
617 + (du_dz+dw_dx)*(du_dz+dw_dx)
618 + (dw_dy+dv_dz)*(dw_dy+dv_dz)));
619 }
620 }
621 }
622 }
623
624
625 // Process internal faces
626
627 for (; fac<nb_faces_; fac++)
628 {
629 poly1 = face_voisins(fac,0);
630 poly2 = face_voisins(fac,1);
631 double a=volumes(poly1)/(volumes(poly1)+volumes(poly2));
632 double b=volumes(poly2)/(volumes(poly1)+volumes(poly2));
633
634 du_dx=a*gradient_elem(poly1,0,0) + b*gradient_elem(poly2,0,0);
635 du_dy=a*gradient_elem(poly1,0,1) + b*gradient_elem(poly2,0,1);
636 dv_dx=a*gradient_elem(poly1,1,0) + b*gradient_elem(poly2,1,0);
637 dv_dy=a*gradient_elem(poly1,1,1) + b*gradient_elem(poly2,1,1);
638
639 // Determination of the production term
640
641 P(fac) = ( 2*(du_dx*du_dx + dv_dy*dv_dy) + ((du_dy+dv_dx)*(du_dy+dv_dx)));
642
643 if (dimension == 3)
644 {
645 du_dz=a*gradient_elem(poly1,0,2) + b*gradient_elem(poly2,0,2);
646 dv_dz=a*gradient_elem(poly1,1,2) + b*gradient_elem(poly2,1,2);
647 dw_dx=a*gradient_elem(poly1,2,0) + b*gradient_elem(poly2,2,0);
648 dw_dy=a*gradient_elem(poly1,2,1) + b*gradient_elem(poly2,2,1);
649 dw_dz=a*gradient_elem(poly1,2,2) + b*gradient_elem(poly2,2,2);
650
651 // Determination of the production term
652
653 P(fac) = (2*( du_dx*du_dx + dv_dy*dv_dy + dw_dz*dw_dz )
654 + ( (du_dy+dv_dx)*(du_dy+dv_dx)
655 + (du_dz+dw_dx)*(du_dz+dw_dx)
656 + (dw_dy+dv_dz)*(dw_dy+dv_dz)));
657 }
658 }
659}
660
661
662void Traitement_particulier_NS_Brech_VEF::
663calculer_terme_destruction_K(const Domaine_VEF& domaine_VEF,
664 const Domaine_Cl_VEF& zcl_VEF,
665 DoubleVect& G,const DoubleTab& temper,
666 const DoubleVect& beta,const DoubleVect& gravite) const
667{
668 //
669 // G is discretised like K and Eps, i.e. at face centres
670 //
671 // --> ----->
672 // G(face) = beta alpha_t(face) G . gradT(face)
673 //
674 int nb_elem = domaine_VEF.nb_elem();
675 int nb_faces= domaine_VEF.nb_faces();
676 DoubleTrav u_teta(nb_faces,dimension);
677 DoubleTrav gradient_elem(nb_elem,dimension);
678 u_teta=0;
679 gradient_elem=0;
680
681 // const Domaine& le_dom=domaine_VEF.domaine();
682 // int nb_faces_elem = le_dom.nb_faces_elem();
683 // const IntTab& elem_faces = domaine_VEF.elem_faces();
684 const IntTab& face_voisins = domaine_VEF.face_voisins();
685 const DoubleTab& face_normales = domaine_VEF.face_normales();
686 int premiere_face_entier = domaine_VEF.premiere_face_int();
687 const DoubleVect& volumes = domaine_VEF.volumes();
688 int nb_faces_ = domaine_VEF.nb_faces();
689 int elem1,elem2,fac;
690
691 //DoubleVect coef(Objet_U::dimension);
692 // const IntTab& les_elem_faces = domaine_VEF.elem_faces();
693
694 // Compute the temperature gradient:
695
696 // Process boundary faces
697 int n_bord;
698 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
699 {
700 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
701 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
702 int ndeb = le_bord.num_premiere_face();
703 int nfin = ndeb + le_bord.nb_faces();
704
705 if (sub_type(Periodique,la_cl.valeur()))
706 {
707 for (fac=ndeb; fac<nfin; fac++)
708 {
709 elem1=face_voisins(fac,0);
710 elem2=face_voisins(fac,1);
711 for (int i=0; i<dimension; i++)
712 {
713 gradient_elem(elem1,i) += 0.5*face_normales(fac,i)*temper(fac);
714 gradient_elem(elem2,i) -= 0.5*face_normales(fac,i)*temper(fac);
715 }
716 }
717 }
718 else
719 {
720 for (int fac2=ndeb; fac2<nfin; fac2++)
721 {
722 elem1=face_voisins(fac2,0);
723 for (int i=0; i<dimension; i++)
724 gradient_elem(elem1,i) += face_normales(fac2,i)*temper(fac2);
725 }
726 }
727 }
728
729 // Process internal faces
730 for (fac = premiere_face_entier ; fac<nb_faces_; fac++)
731 {
732 elem1=face_voisins(fac,0);
733 elem2=face_voisins(fac,1);
734 for (int i=0; i<dimension; i++)
735 {
736 gradient_elem(elem1,i) += face_normales(fac,i)*temper(fac);
737 gradient_elem(elem2,i) -= face_normales(fac,i)*temper(fac);
738 }
739 }
740
741 for (int elem=0; elem<nb_elem; elem++)
742 for (int i=0; i<dimension; i++)
743 gradient_elem(elem,i) /= volumes(elem);
744
745 // Compute u_teta:
746
747 // Process boundary faces
748
749 for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
750 {
751 const Cond_lim& la_cl = zcl_VEF.les_conditions_limites(n_bord);
752 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
753 int ndeb = le_bord.num_premiere_face();
754 int nfin = ndeb + le_bord.nb_faces();
755
756 if (sub_type(Periodique,la_cl.valeur()))
757 {
758 for (fac=ndeb; fac<nfin; fac++)
759 {
760 elem1 = face_voisins(fac,0);
761 elem2 = face_voisins(fac,1);
762 double a=volumes(elem1)/(volumes(elem1)+volumes(elem2));
763 double b=volumes(elem2)/(volumes(elem1)+volumes(elem2));
764 for (int i=0; i<dimension; i++)
765 u_teta(fac,i)=a*beta(elem1)*gradient_elem(elem1,i)
766 +b*beta(elem2)*gradient_elem(elem2,i);
767 }
768 }
769 else
770 {
771 for (fac=ndeb; fac< nfin; fac++)
772 {
773 elem1 = face_voisins(fac,0);
774 for (int i=0; i<dimension; i++)
775 u_teta(fac,i)=beta(elem1)*gradient_elem(elem1,i);
776 }
777 }
778 }
779
780 for (; fac<nb_faces_; fac++)
781 {
782 elem1 = face_voisins(fac,0);
783 elem2 = face_voisins(fac,1);
784 double a=volumes(elem1)/(volumes(elem1)+volumes(elem2));
785 double b=volumes(elem2)/(volumes(elem1)+volumes(elem2));
786 for (int i=0; i<dimension; i++)
787 u_teta(fac,i)=a*beta(elem1)*gradient_elem(elem1,i)+b*beta(elem2)*gradient_elem(elem2,i);
788 }
789 // -------> ------>
790 // Compute gravity . u_teta
791 //
792 G = 0;
793 for (fac=0; fac< nb_faces_; fac++)
794 {
795 if (dimension == 2)
796 G[fac] = gravite(0)*u_teta(fac,0) + gravite(1)*u_teta(fac,1);
797 else if (dimension == 3)
798 G[fac] = gravite(0)*u_teta(fac,0) + gravite(1)*u_teta(fac,1) + gravite(2)*u_teta(fac,2);
799 }
800}
801
Champ_Inc_base
Class Champ_Inc_base.
Definition Champ_Inc_base.h:53
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Returns the array of field values at the current time.
Definition Champ_Inc_base.h:77
Champ_P1NC::calcul_gradient
static DoubleTab & calcul_gradient(const DoubleTab &, DoubleTab &, const Domaine_Cl_VEF &)
Definition Champ_P1NC.cpp:912
Champ_base
class Champ_base This class is the base of the fields hierarchy.
Definition Champ_base.h:43
Champ_base::valeur_aux_elems
virtual DoubleTab & valeur_aux_elems(const DoubleTab &positions, const IntVect &les_polys, DoubleTab &valeurs) const
Causes an error! Must be overridden by derived classes.
Definition Champ_base.cpp:245
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:38
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Returns the i-th boundary condition.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:53
Domaine_VF::nb_faces
int nb_faces() const
Returns the total number of faces.
Definition Domaine_VF.h:471
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::volumes
double volumes(int i) const
Definition Domaine_VF.h:113
Domaine_VF::premiere_face_int
int premiere_face_int() const
A face is internal if and only if it separates two elements.
Definition Domaine_VF.h:463
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
Returns the neighbouring element of num_face in direction i.
Definition Domaine_VF.h:418
Domaine_dis_base
class Domaine_dis_base This class is the base of the hierarchy of discretized domains.
Definition Domaine_dis_base.h:39
Domaine_dis_base::nb_elem
int nb_elem() const
Definition Domaine_dis_base.h:49
Domaine_dis_base::nb_front_Cl
int nb_front_Cl() const
Definition Domaine_dis_base.h:53
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Front_VF::nb_faces
int nb_faces() const
Definition Front_VF.h:53
Front_VF::num_premiere_face
int num_premiere_face() const
Definition Front_VF.h:63
Motcle
A character string (Nom) in uppercase.
Definition Motcle.h:26
Motcles
An array of Motcle objects.
Definition Motcle.h:63
Motcles::search
int search(const Motcle &t) const
Definition Motcle.cpp:319
Objet_U::Entree
friend class Entree
Definition Objet_U.h:71
Objet_U::dimension
static int dimension
Definition Objet_U.h:94
Objet_U::readOn
virtual Entree & readOn(Entree &)
Reads an Objet_U from an input stream. Virtual method to override.
Definition Objet_U.cpp:289
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Writes the object to an output stream. Virtual method to override.
Definition Objet_U.cpp:278
Probleme_base
class Probleme_base It is a Probleme_U that is not a coupling.
Definition Probleme_base.h:55
Probleme_base::get_champ
const Champ_base & get_champ(const Motcle &nom) const override
Definition Probleme_base.cpp:836
Process::mp_sum
static double mp_sum(double)
Computes the sum of x over all processors in the current group.
Definition Process.cpp:145
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Process::je_suis_maitre
static int je_suis_maitre()
Returns 1 if on the master processor of the current group (i.e. me() == 0), 0 otherwise.
Definition Process.cpp:82
SolveurSys::resoudre_systeme
int resoudre_systeme(const Matrice_Base &matrice, const DoubleVect &secmem, DoubleVect &solution)
Definition SolveurSys.cpp:48
Sortie
Base class for output streams.
Definition Sortie.h:52
Traitement_particulier_NS_Brech_VEF
class Traitement_particulier_Brech_VEF This class performs special post-processing
Definition Traitement_particulier_NS_Brech_VEF.h:35
Traitement_particulier_NS_Brech_VEF::lire
Entree & lire(Entree &) override
Definition Traitement_particulier_NS_Brech_VEF.cpp:62
Traitement_particulier_NS_Brech_VEF::Traitement_particulier_NS_Brech_VEF
Traitement_particulier_NS_Brech_VEF()
Definition Traitement_particulier_NS_Brech_VEF.cpp:26
Traitement_particulier_NS_Brech_VEF::ch_p
Champ_Fonc_P1NC ch_p
Definition Traitement_particulier_NS_Brech_VEF.h:46
Traitement_particulier_NS_Brech_VEF::post_traitement_particulier
void post_traitement_particulier() override
Definition Traitement_particulier_NS_Brech_VEF.cpp:219
Traitement_particulier_NS_Brech_VEF::ch_ri
Champ_Fonc_P1NC ch_ri
Definition Traitement_particulier_NS_Brech_VEF.h:45
Traitement_particulier_NS_VEF
class Traitement_particulier_VEF This class performs no special post-processing
Definition Traitement_particulier_NS_VEF.h:31
Traitement_particulier_NS_base::OBS_PTR
OBS_PTR(Navier_Stokes_std) mon_equation
Traitement_particulier_NS_base::champs_compris_
Champs_compris champs_compris_
Definition Traitement_particulier_NS_base.h:68