TrioCFD 1.9.8
TrioCFD documentation
Loading...
Searching...
No Matches
Op_Conv_EF_VEF_P1NC_Stab.cpp
1/****************************************************************************
2* Copyright (c) 2026, 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 <Op_Conv_EF_VEF_P1NC_Stab.h>
17#include <Champ_P1NC.h>
18#include <BilanQdmVEF.h>
19#include <Sous_Domaine.h>
20#include <Sous_domaine_dis_base.h>
21#include <Schema_Temps_base.h>
22#include <Debog.h>
23#include <Porosites_champ.h>
24#include <Sous_domaine_VF.h>
25#include <Probleme_base.h>
26#include <ArrOfBit.h>
27#include <TRUSTVects.h>
28#include <SFichier.h>
29#include <TRUSTList.h>
30#include <TRUSTTabs.h>
31#include <Neumann_sortie_libre.h>
32#include <Dirichlet_homogene.h>
33#include <Neumann_homogene.h>
34#include <Periodique.h>
35#include <Symetrie.h>
36#include <Echange_impose_base.h>
37#include <Array_tools.h>
38
39Implemente_instanciable(Op_Conv_EF_VEF_P1NC_Stab,"Op_Conv_EF_Stab_VEF_P1NC",Op_Conv_VEF_Face);
40
41// XD sous_zone_valeur objet_lecture nul NO_BRACE Two words.
42// XD attr sous_zone ref_sous_zone sous_zone REQ sous zone
43// XD attr valeur floattant valeur REQ value
44
45// XD listsous_zone_valeur listobj nul NO_BRACE sous_zone_valeur NO_COMMA List of groups of two words.
46
47// XD convection_ef_stab convection_deriv ef_stab BRACE Keyword for a VEF convective scheme.
48// XD attr alpha floattant alpha OPT To weight the scheme centering with the factor double (between 0 (full centered)
49// XD_CONT and 1 (mix between upwind and centered), by default 1). For scalar equation, it is adviced to use alpha=1 and
50// XD_CONT for the momentum equation, alpha=0.2 is adviced.
51// XD attr test entier test OPT Developer option to compare old and new version of EF_stab
52// XD attr tdivu rien tdivu OPT To have the convective operator calculated as div(TU)-TdivU(=UgradT).
53// XD attr old rien old OPT To use old version of EF_stab scheme (default no).
54// XD attr volumes_etendus rien volumes_etendus OPT Option for the scheme to use the extended volumes (default, yes).
55// XD attr volumes_non_etendus rien volumes_non_etendus OPT Option for the scheme to not use the extended volumes
56// XD_CONT (default, no).
57// XD attr amont_sous_zone ref_sous_zone amont_sous_zone OPT Option to degenerate EF_stab scheme into Amont (upwind)
58// XD_CONT scheme in the sub zone of name sz_name. The sub zone may be located arbitrarily in the domain but the more
59// XD_CONT often this option will be activated in a zone where EF_stab scheme generates instabilities as for free outlet
60// XD_CONT for example.
61// XD attr alpha_sous_zone listsous_zone_valeur alpha_sous_zone OPT Option to change locally the alpha value on N
62// XD_CONT sub-zones named sub_zone_name_I. Generally, it is used to prevent from a local divergence by increasing
63// XD_CONT locally the alpha parameter.
64
65Sortie& Op_Conv_EF_VEF_P1NC_Stab::printOn(Sortie& s ) const
66{
67 return s << que_suis_je();
68}
69
70Entree& Op_Conv_EF_VEF_P1NC_Stab::readOn(Entree& s )
71{
72 //Les mots a reconnaitre
73 Motcle motlu, accouverte = "{" , accfermee = "}" ;
74 Motcles les_mots(11);
75 {
76 les_mots[0] = "alpha";
77 les_mots[1] = "test";
78 les_mots[2] = "TdivU";
79 les_mots[3] = "old";
80 les_mots[4] = "volumes_etendus";
81 les_mots[5] = "volumes_non_etendus";
82 les_mots[6] = "amont_sous_domaine";
83 les_mots[7] = "amont_sous_zone";
84 les_mots[8] = "nouvelle_matrice_implicite";
85 les_mots[9] = "alpha_sous_domaine";
86 les_mots[10] = "alpha_sous_zone";
87 }
88
89 s >> motlu;
90 if (motlu!=accouverte)
91 {
92 Cerr << "Erreur Op_Conv_EF_VEF_P1NC_Stab::readOn()" << finl;
93 Cerr << "Depuis la 1.5.3, la syntaxe du mot cle EF_stab a change." << finl;
94 Cerr << "Il faut commencer par une accolade ouvrante {" << finl;
95 Cerr << "et les options eventuelles sont entre les accolades:" << finl;
96 Cerr << "Convection { EF_stab } -> Convection { EF_stab { } }" << finl;
97 Process::exit();
98 }
99 s >> motlu;
100
101 while(motlu!=accfermee)
102 {
103 int rang=les_mots.search(motlu);
104
105 switch(rang)
106 {
107 case 0 :
108
109 s >> alpha_;
110 break;
111
112 case 1 :
113
114 test_=1;
115 break;
116
117 case 2 :
118
119 is_compressible_=1;
120 break;
121
122 case 3 :
123
124 old_=1;
125 break;
126
127 case 4 :
128
129 volumes_etendus_=1;
130 break;
131
132 case 5 :
133
134 volumes_etendus_=0;
135 break;
136
137 case 6 :
138 case 7 :
139 sous_domaine=true;
140 s >> nom_sous_domaine;
141 break;
142 case 8 :
143 s >> new_jacobienne_;
144 break;
145 case 9 :
146 case 10 :
147 ssz_alpha=true;
148 s >> nb_ssz_alpha;
149 noms_ssz_alpha.dimensionner(nb_ssz_alpha);
150 alpha_ssz.resize(nb_ssz_alpha);
151 for (int i=0; i<nb_ssz_alpha; i++)
152 {
153 s>>noms_ssz_alpha[i];
154 s>>alpha_ssz(i);
155 }
156 break;
157 default :
158
159 Cerr << "Erreur Op_Conv_EF_VEF_P1NC_Stab::readOn()" << finl;
160 Cerr << "Mot clef " << motlu << " non connu." << finl;
161 Cerr << "Sortie du programme." << finl;
162 Process::exit();
163
164 }//fin du switch
165
166 //Suite de la lecture
167 s >> motlu;
168
169 }//fin du while
170
171 return s ;
172}
173
174
175static KOKKOS_INLINE_FUNCTION double maximum(const double x,
176 const double y)
177{
178 return x<y ? y : x;
179}
180
181static KOKKOS_INLINE_FUNCTION double maximum(const double x,
182 const double y,
183 const double z)
184{
185 return maximum(maximum(x,y),z);
186}
187
188static KOKKOS_INLINE_FUNCTION double minimum(const double x,
189 const double y)
190{
191 return x>y ? y : x;
192}
193
194static KOKKOS_INLINE_FUNCTION double Dij(int elem,
195 int face_loc_i,
196 int face_loc_j,
197 CDoubleTabView3 Kij)
198{
199 const double kij=Kij(elem,face_loc_i,face_loc_j);
200 const double kji=Kij(elem,face_loc_j,face_loc_i);
201 return maximum(-kij,-kji,0);
202}
203
204static inline double Dij(int elem,
205 int face_loc_i,
206 int face_loc_j,
207 const DoubleTab& Kij)
208{
209 const double kij=Kij(elem,face_loc_i,face_loc_j);
210 const double kji=Kij(elem,face_loc_j,face_loc_i);
211 return maximum(-kij,-kji,0);
212}
213
214static KOKKOS_INLINE_FUNCTION double limiteur(double r)
215{
216 return r<=0 ? 0 : Kokkos::fmax(Kokkos::fmin(2,r),Kokkos::fmin(1,2*r));//SuperBee
217}
218
219KOKKOS_INLINE_FUNCTION double formule_Id_2D(int n)
220{
221 return 1./3;
222}
223
224KOKKOS_INLINE_FUNCTION double formule_Id_3D(int n)
225{
226 return 0.25;
227}
228
229KOKKOS_INLINE_FUNCTION double formule_2D(int n)
230{
231 switch(n)
232 {
233 case 0:
234 return 1./3;
235 case 1:
236 return 0.5;
237 case 2:
238 return 1.;
239 default:
240 Process::Kokkos_exit("Erreur Op_Conv_EF_VEF_P1NC_Stab::formule_2D() ");
241 }
242 return 0.;
243}
244
245KOKKOS_INLINE_FUNCTION double formule_3D(int n)
246{
247 switch(n)
248 {
249 case 0:
250 return 0.25;
251 case 1:
252 return 1./3;
253 case 2:
254 return 0.5;
255 case 3:
256 return 1.;
257 default:
258 Process::Kokkos_exit("Erreur Op_Conv_EF_VEF_P1NC_Stab::formule_3D() ");
259 }
260 return 0.;
261}
262
263
264
265////////////////////////////////////////////////////////////////////
266//
267// Implementation des fonctions
268//
269// de la classe Op_Conv_EF_VEF_P1NC_Stab
270//
271////////////////////////////////////////////////////////////////////
272
273void Op_Conv_EF_VEF_P1NC_Stab::reinit_conv_pour_Cl(const DoubleTab& transporte,const IntList& faces, const DoubleTabs& valeurs_faces, const DoubleTab& tab_vitesse, DoubleTab& resu) const
274{
275 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
276 const Domaine_Cl_VEF& domaine_Cl_VEF=la_zcl_vef.valeur();
277 const DoubleTab& face_normales=domaine_VEF.face_normales();
278 const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
279 const int nb_comp=transporte.line_size();
280 int n_bord=0, num1=0, num2=0, ind_face=0,facei=0, dim=0;
281 double psc=0.;
282
283 const DoubleVect& transporteV= transporte;
284
285 //Boucle pour dimensionner la liste "faces"
286 for (n_bord=0; n_bord<nb_bord; n_bord++)
287 {
288 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
289 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
290 num1 = 0;
291 num2 = le_bord.nb_faces_tot();
292
293 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
294 {
295 const Neumann_sortie_libre& la_sortie_libre
296 = ref_cast(Neumann_sortie_libre, la_cl.valeur());
297 ToDo_Kokkos("critical");
298 for (ind_face=num1; ind_face<num2; ind_face++)
299 {
300 facei = le_bord.num_face(ind_face);
301
302 psc=0.;
303 for (dim=0; dim<dimension; dim++)
304 psc-=tab_vitesse(facei,dim)*face_normales(facei,dim);
305
306 if (psc>0)
307 for (dim=0; dim<nb_comp; dim++)
308 resu(facei,dim)+=psc*(la_sortie_libre.val_ext(facei-num1,dim)-transporteV[facei*nb_comp+dim]);
309
310 }
311 }//fin du if sur "Neumann_sortie_libre"
312 }
313}
314
315void Op_Conv_EF_VEF_P1NC_Stab::calculer_coefficients_operateur_centre(DoubleTab& tab_Kij,const int nb_comp, const DoubleTab& tab_vitesse) const
316{
317 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
318 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
319 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
320 const int nb_faces_elem = domaine_VEF.elem_faces().line_size();
321 int dim = Objet_U::dimension;
322
323 assert(tab_Kij.nb_dim()==3);
324 assert(tab_Kij.dimension(0)==nb_elem_tot);
325 assert(tab_Kij.dimension(1)==tab_Kij.dimension(2));
326 assert(tab_Kij.dimension(1)==nb_faces_elem);
327
328 //
329 //Calcul des coefficients de l'operateur
330 //
331 CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
332 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
333 CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
334 CDoubleTabView vitesse = tab_vitesse.view_ro();
335 DoubleTabView3 Kij = tab_Kij.view_rw<3>();
336 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
337 range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
338 const int elem, const int face_loci)
339 {
340 int face_i=elem_faces(elem,face_loci);
341
342 double signei=1.0;
343 if(face_voisins(face_i,0)!=elem) signei=-1.0;
344
345 double psci=0;
346 for(int comp=0; comp<dim; comp++)
347 psci+=vitesse(face_i,comp)*face_normales(face_i,comp);
348 psci*=signei;
349
350 //Kij(elem,face_loci,face_loci)=0.;
351 for(int face_locj=face_loci+1; face_locj<nb_faces_elem; face_locj++)
352 {
353 int face_j=elem_faces(elem,face_locj);
354 double signej=1.0;
355 if(face_voisins(face_j,0)!=elem)
356 signej=-1.0;
357
358 double pscj=0;
359 //psci=0;
360 for(int comp=0; comp<dim; comp++)
361 pscj+=vitesse(face_j,comp)*face_normales(face_j,comp);
362 pscj*=signej;
363
364 Kokkos::atomic_add(&Kij(elem,face_loci,face_locj), -1./nb_faces_elem*pscj);
365 Kokkos::atomic_add(&Kij(elem,face_loci,face_loci), +1./nb_faces_elem*pscj);
366 Kokkos::atomic_add(&Kij(elem,face_locj,face_loci), -1./nb_faces_elem*psci);
367 Kokkos::atomic_add(&Kij(elem,face_locj,face_locj), +1./nb_faces_elem*psci);
368 }//fin du for sur "face_locj"
369 });//fin du for sur "elem"
370 end_gpu_timer(__KERNEL_NAME__);
371 //
372 // Correction des Kij pour Dirichlet !
373 //
374 {
375 int nb_bord=domaine_Cl_VEF.nb_cond_lim();
376 for (int n_bord=0; n_bord<nb_bord; n_bord++)
377 {
378 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
379
380 if ( (sub_type(Dirichlet,la_cl.valeur()))
381 || (sub_type(Dirichlet_homogene,la_cl.valeur()))
382 )
383 {
384 // GF on retire le test sur nb_comp car sinon en scalaire on fait tjs du volume etendu
385 if ( volumes_etendus_ )
386 {
387 CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
388 CIntTabView elem_faces_dirichlet_v = elem_faces_dirichlet_.view_ro();
389 CIntArrView elem_nb_faces_dirichlet_v = elem_nb_faces_dirichlet_.view_ro();
390 CIntArrView elem_faces_frontiere_v = elem_faces_frontiere[n_bord].view_ro();
391 const int elem_faces_frontiere_size = elem_faces_frontiere[n_bord].size_array();
392 //
393 //Modification des coefficients de la matrice
394 //
395 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
396 Kokkos::RangePolicy<>(0,elem_faces_frontiere_size), KOKKOS_LAMBDA(
397 const int elem_ind)
398 {
399 const int elem=elem_faces_frontiere_v(elem_ind);
400 assert(elem!=-1);
401 const int nb_faces_bord = elem_nb_faces_dirichlet_v(elem);
402
403 //
404 //Calcul du coefficient ponderateur
405 //
406
407 const double coeff = (dim==2) ?
408 nb_faces_bord/2. : nb_faces_bord*nb_faces_bord/6.-nb_faces_bord/3.+1./2;
409
410 //
411 //Fin du calcul du coefficient ponderateur
412 //
413
414 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
415 {
416 const int face_i = elem_faces(elem,face_loc_i);
417
418 /* On determine si "face_i" est une face de Dirichlet */
419 bool is_not_on_boundary = true;
420 for (int f_loc=0; f_loc<nb_faces_bord; f_loc++)
421 is_not_on_boundary&=(face_i!=elem_faces_dirichlet_v(elem,f_loc));
422
423 if (is_not_on_boundary) /* si "face_i" n'est pas de Dirichlet */
424 {
425 for (int face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
426 {
427 /* Modification des coefficients de la matrice */
428 /* selon les fonctions de formes etendues */
429 for (int f_loc=0; f_loc<nb_faces_bord; f_loc++)
430 {
431 const int face_bord = elem_faces_dirichlet_v(elem,f_loc);
432 const int face_loc_k = num_fac_loc(face_bord,0);
433 assert(face_loc_k>=0);
434 assert(face_loc_k<nb_faces_elem);
435
436 const double kkj = Kij(elem,face_loc_k,face_loc_j);
437 Kij(elem,face_loc_i,face_loc_j) += coeff*kkj;
438 }//fin du for sur "f_loc"
439
440 }//fin du for sur "face_loc_k"
441
442 }//fin du if
443
444 }//fin du for sur "face_loc_i"
445
446
447 //
448 //Remise a zero des coefficients associes aux noeuds
449 //qui se situent sur la frontiere de Dirichlet
450 //
451 {
452 for (int f_loc=0; f_loc<nb_faces_bord; f_loc++)
453 {
454 const int face_bord = elem_faces_dirichlet_v(elem,f_loc);
455 const int face_loc_j = num_fac_loc(face_bord,0);
456 assert(face_loc_j>=0);
457 assert(face_loc_j<nb_faces_elem);
458
459 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
460 Kij(elem,face_loc_j,face_loc_i)=0;
461 }//fin du for sur "f_loc"
462 }
463 //
464 //Fin de la remise a zero
465 //
466
467 //
468 //Pour retrouver le schema LED
469 //
470 {
471 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
472 {
473 double sum=0.;
474 for (int face_loc_k=0; face_loc_k<nb_faces_elem; face_loc_k++)
475 {
476 sum+=Kij(elem,face_loc_i,face_loc_k);
477 }
478 //Cerr << "somme apres : " << sum << finl;
479 Kij(elem,face_loc_i,face_loc_i) -= sum;//car div(u)=0!
480 }
481 }
482 //
483 //Fin du schema LED
484 //
485
486 });// for face
487
488 end_gpu_timer(__KERNEL_NAME__);
489
490 //
491 //Fin de la modification des coefficients de la matrice
492 //
493
494 }//fin du if sur "volumes_etendus_"
495
496 }// sub_type Dirichlet
497
498
499 else if (sub_type(Neumann,la_cl.valeur()) || sub_type(Neumann_homogene,la_cl.valeur()) || sub_type(Neumann_val_ext,la_cl.valeur()))
500 {
501 //On ne fait rien
502 }//fin du if sur Neumann
503
504 else if (sub_type(Symetrie,la_cl.valeur()))
505 {
506 //On ne fait rien
507 }//fin du if sur Symetrie
508
509 else if (sub_type(Periodique,la_cl.valeur()))
510 {
511 //On ne fait rien
512 }//fin du if sur Periodique
513
514 else if (sub_type(Echange_impose_base,la_cl.valeur()))
515 {
516 //On ne fait rien
517 }//fin du if sur Echange_impose_base
518
519 else
520 {
521 Cerr << "Erreur Op_Conv_EF_VEF_P1NC_Stab::calculer_coefficients_operateur_centre()" << finl;
522 Cerr << "Condition aux limites " << la_cl.que_suis_je() << " non codee." << finl;
523 Cerr << "Sortie du programme." << finl;
524 Process::exit();
525 }//fin du else sur les autres conditions aux limites
526
527 }//fin des conditions aux limites
528 }
529
530 //
531 // Fin de la correction des Kij
532 //
533}
534void Op_Conv_EF_VEF_P1NC_Stab::remplir_fluent() const
535{
536 // calcul de la CFL.
537 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
538 const int nb_faces = domaine_VEF.nb_faces();
539 int nb_comp = Objet_U::dimension;
540 CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
541 CDoubleTabView tab_vitesse = vitesse_->valeurs().view_ro();
542 DoubleArrView fluent = fluent_.view_rw();
543 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
544 Kokkos::RangePolicy<>(0, nb_faces), KOKKOS_LAMBDA(
545 const int num_face)
546 {
547 double psc=0.;
548 for (int i=0; i<nb_comp; i++)
549 psc+=tab_vitesse(num_face,i)*face_normales(num_face,i);
550 fluent(num_face)=std::fabs(psc);
551 });
552 end_gpu_timer(__KERNEL_NAME__);
553}
554
555
556DoubleTab& Op_Conv_EF_VEF_P1NC_Stab::ajouter(const DoubleTab& transporte_2,
557 DoubleTab& resu) const
558{
559 DoubleTab sauv(resu);
560 resu=0;
561
562 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
563 const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
564 const DoubleTab& vitesse_2=la_vitesse.valeurs();
565 const DoubleVect& porosite_face = equation().milieu().porosite_face();
566
567 const IntTab& elem_faces = domaine_VEF.elem_faces();
568
569 const int marq=phi_u_transportant(equation());
570 const int nb_faces_elem=elem_faces.dimension(1);
571 assert(nb_faces_elem==(dimension+1));
572 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
573 int nb_comp=resu.line_size();
574
575 DoubleTab transporte_;
576 DoubleTab vitesse_face_;
577 // soit on a transporte=phi*transporte_ et vitesse=vitesse_
578 // soit transporte=transporte_ et vitesse=phi*vitesse_
579 // cela depend si on transporte avec phi u ou avec u.
580 const DoubleTab& tab_vitesse=modif_par_porosite_si_flag(vitesse_2,vitesse_face_,marq,porosite_face);
581
582 // soit on a transporte=phi*transporte_ et vitesse=vitesse_
583 // soit transporte=transporte_ et vitesse=phi*vitesse_
584 // cela depend si on transporte avec phi u ou avec u.
585 const DoubleTab& transporte=modif_par_porosite_si_flag(transporte_2,transporte_,!marq,porosite_face);
586
587 DoubleTrav Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
588 //Initialisation du tableau flux_bords_ pour le calcul des pertes de charge
589 flux_bords_.resize(domaine_VEF.nb_faces_bord(),nb_comp);
590 calculer_flux_bords(Kij,tab_vitesse,transporte);
591
592 if (!old_)
593 {
594 calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
595 if (is_compressible_) ajouter_partie_compressible(transporte,resu,tab_vitesse);
596
597 ajouter_operateur_centre(Kij,transporte,resu);
598
599 ajouter_diffusion(Kij,transporte,resu);
600 ajouter_antidiffusion(Kij,transporte,resu);
601 mettre_a_jour_pour_periodicite(resu);
602
603 }
604 else
605 {
606 assert(old_==1);
607 ajouter_old(transporte,resu,tab_vitesse);
608 }
609
610 //Pour tenir compte des conditions de Neumann sortie libre
611 IntList NeumannFaces;
612 DoubleTabs ValeursNeumannFaces;
613 reinit_conv_pour_Cl(transporte,NeumannFaces,ValeursNeumannFaces,tab_vitesse,resu);
614
615 resu+=sauv;
616
617 if (test_) test(transporte,resu,tab_vitesse);
618 modifier_flux(*this);
619 return resu;
620}
621
622//Correction pour le poreux : on rajoute la partie en T div(u)
623//Variable transportee : T
624//Variable transportante : u
625//REMARQUE : il ne FAUT SURTOUT PAS utiliser le tableau Kij car par
626//construction celui-ci est telle que sum_{j} Kij =0 ce qui revient a
627//imposer une vitesse a divergence nulle par element. Ce qui est
628//problematique quand on est en compressible
629DoubleTab& Op_Conv_EF_VEF_P1NC_Stab::ajouter_partie_compressible(const DoubleTab& transporte,
630 DoubleTab& resu, const DoubleTab& vitesse_2) const
631{
632 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
633 const IntTab& elem_faces=domaine_VEF.elem_faces();
634 const IntTab& face_voisins = domaine_VEF.face_voisins();
635 const DoubleTab& face_normales=domaine_VEF.face_normales();
636 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
637 const int nb_faces_elem=elem_faces.line_size();
638
639 //Pour tenir compte de la porosite
640 const int marq = phi_u_transportant(equation());
641 const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
642 const DoubleVect& porosite_face = equation().milieu().porosite_face();
643
644 DoubleTab tab_vitesse(vitesse_->valeurs());
645 DoubleTabView tab_vitesse_v = tab_vitesse.view_rw();
646 CDoubleArrView porosite_face_v = porosite_face.view_ro();
647
648 const int vit_size0 = tab_vitesse.dimension(0);
649 const int vit_size1 = tab_vitesse.dimension(1);
650 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
651 Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {vit_size0, vit_size1}), KOKKOS_LAMBDA(
652 const int i, const int j)
653 {
654 tab_vitesse_v(i,j)*=porosite_face_v(i);
655 });
656 end_gpu_timer(__KERNEL_NAME__);
657 const int nb_comp=transporte.line_size();
658 int nb_dim = dimension;
659 int vol_etendus = volumes_etendus_;
660
661 // read only
662 CDoubleTabView face_normales_v = face_normales.view_ro();
663 CIntTabView face_voisins_v = face_voisins.view_ro();
664 CIntTabView elem_faces_v = elem_faces.view_ro();
665 CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
666 CIntArrView elem_nb_faces_dirichlet_v = elem_nb_faces_dirichlet_.view_ro();
667 CDoubleArrView porosite_elem_v = porosite_elem.view_ro();
668 // Write only
669 DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
670
671 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
672 Kokkos::RangePolicy<>(0, nb_elem_tot), KOKKOS_LAMBDA(
673 const int elem)
674 {
675 //Type de l'element : le nombre de faces de Dirichlet
676 //qu'il contient
677 int type_elem=elem_nb_faces_dirichlet_v(elem);
678 double coeff;
679 if (!vol_etendus)
680 coeff = (nb_dim==2) ? formule_Id_2D(type_elem) : formule_Id_3D(type_elem);
681 else
682 coeff = (nb_dim==2) ? formule_2D(type_elem) : formule_3D(type_elem);
683
684 int facei, facei_loc, dim;
685
686 //Calcul de la divergence par element
687 double div=0.;
688 for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
689 {
690 facei=elem_faces_v(elem,facei_loc);
691 double signe=(face_voisins_v(facei,0)==elem)? 1.:-1.;
692
693 for (dim=0; dim<nb_dim; dim++)
694 div+=signe*face_normales_v(facei,dim)*tab_vitesse_v(facei,dim);
695 }
696 div*=coeff;
697 if (!marq) div/=porosite_elem_v(elem);
698
699 //Calcul de la partie compressible
700 for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
701 {
702 facei=elem_faces_v(elem,facei_loc);
703 for (dim=0; dim<nb_comp; dim++)
704 {
705 int ligne=facei*nb_comp+dim;
706 double delta = div*transporteV[ligne];
707 Kokkos::atomic_sub(&resuV[ligne], delta);
708 }
709 }
710 });
711 end_gpu_timer(__KERNEL_NAME__);
712
713 return resu;
714}
715
716void Op_Conv_EF_VEF_P1NC_Stab::calculer_flux_bords(const DoubleTab& Kij, const DoubleTab& tab_vitesse, const DoubleTab& transporte) const
717{
718 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
719 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
720 const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
721 const int nb_comp=transporte.line_size();
722
723 for (int n_bord=0; n_bord<nb_bord; n_bord++)
724 {
725 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
726 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
727 int num1 = 0;
728 int num2 = le_bord.nb_faces();//il ne faut boucler que sur les faces reelles ici
729
730 if ( sub_type(Dirichlet_homogene,la_cl.valeur()) )
731 {
732 //On ne calcule pas le flux aux bords Dirichlet_homogene
733 }
734 else if ( sub_type(Neumann,la_cl.valeur())
735 || sub_type(Neumann_val_ext,la_cl.valeur())
736 || sub_type(Neumann_homogene,la_cl.valeur())
737 || sub_type(Symetrie,la_cl.valeur())
738 || sub_type(Echange_impose_base,la_cl.valeur())
739 || sub_type(Dirichlet,la_cl.valeur())
740 || sub_type(Periodique,la_cl.valeur()) )
741 {
742 CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
743 CIntArrView num_face = le_bord.num_face().view_ro();
744 CDoubleTabView vitesse = tab_vitesse.view_ro();
745 CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
746 DoubleTabView flux_bords = flux_bords_.view_wo();
747 const int nb_dim=Objet_U::dimension;
748 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
749 Kokkos::RangePolicy<>(num1, num2), KOKKOS_LAMBDA(
750 const int ind_face)
751 {
752 int facei = num_face(ind_face);
753 double psc=0.;
754 for (int dim=0; dim<nb_dim; dim++) //tempo fix to compile
755 psc-=vitesse(facei,dim)*face_normales(facei,dim);
756
757 for (int dim=0; dim<nb_comp; dim++)
758 flux_bords(facei,dim)=psc*transporteV[facei*nb_comp+dim];
759 });
760 end_gpu_timer(__KERNEL_NAME__);
761
762 }//fin du if sur "Neumann", "Neumann_homogene", "Symetrie", "Echange_impose_base"
763 else
764 {
765 Cerr << "Erreur Op_Conv_EF_VEF_P1NC_Stab::calculer_flux_bords()" << finl;
766 Cerr << "Condition aux limites " << la_cl.que_suis_je() << " non codee." << finl;
767 Cerr << "Sortie du programme." << finl;
768 Process::exit();
769 }//fin du else sur les autres conditions aux limites
770 }
771}
772
773DoubleTab&
774Op_Conv_EF_VEF_P1NC_Stab::ajouter_operateur_centre(const DoubleTab& tab_Kij, const DoubleTab& transporte, DoubleTab& resu) const
775{
776 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
777 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
778 const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
779 const int nb_comp=transporte.line_size();
780
781 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
782 CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
783 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
784 DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
785 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
786 range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
787 const int elem, const int facei_loc)
788 {
789 int facei = elem_faces(elem, facei_loc);
790 for (int facej_loc = facei_loc + 1; facej_loc < nb_faces_elem; facej_loc++)
791 {
792 int facej = elem_faces(elem, facej_loc);
793 double kij = Kij(elem, facei_loc, facej_loc);
794 double kji = Kij(elem, facej_loc, facei_loc);
795
796 for (int dim = 0; dim < nb_comp; dim++)
797 {
798 int ligne = facei * nb_comp + dim;
799 int colonne = facej * nb_comp + dim;
800 double delta = transporteV[colonne] - transporteV[ligne];
801 Kokkos::atomic_add(&resuV[ligne], kij * delta);
802 Kokkos::atomic_sub(&resuV[colonne], kji * delta);
803 }
804 }
805 });
806 end_gpu_timer(__KERNEL_NAME__);
807 return resu;
808}
809
810DoubleTab&
811Op_Conv_EF_VEF_P1NC_Stab::ajouter_diffusion(const DoubleTab& tab_Kij, const DoubleTab& transporte, DoubleTab& resu) const
812{
813 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
814 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
815 const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
816 const int nb_comp=transporte.line_size();
817
818 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
819 CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
820 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
821 CDoubleArrView alpha_tab = alpha_tab_.view_ro();
822 DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
823 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
824 range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
825 const int elem, const int facei_loc)
826 {
827 int facei=elem_faces(elem,facei_loc);
828 for (int facej_loc=facei_loc+1; facej_loc<nb_faces_elem; facej_loc++)
829 {
830 int facej=elem_faces(elem,facej_loc);
831 double dij=Dij(elem,facei_loc,facej_loc,Kij);
832 double coeffij=alpha_tab[facei]*dij;
833 double coeffji=alpha_tab[facej]*dij;
834
835 for (int dim=0; dim<nb_comp; dim++)
836 {
837 int ligne=facei*nb_comp+dim;
838 int colonne=facej*nb_comp+dim;
839 double delta=transporteV[colonne]-transporteV[ligne];
840
841 //ATTENTION AU SIGNE : ici on code +div(uT)
842 //REMARQUE : on utilise la symetrie de l'operateur
843 Kokkos::atomic_add(&resuV[ligne], coeffij*delta);
844 Kokkos::atomic_sub(&resuV[colonne], coeffji*delta);
845 }
846 }
847 });
848 end_gpu_timer(__KERNEL_NAME__);
849 return resu;
850}
851
852DoubleTab&
853Op_Conv_EF_VEF_P1NC_Stab::ajouter_antidiffusion(const DoubleTab& tab_Kij, const DoubleTab& transporte, DoubleTab& resu) const
854{
855 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
856 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
857 const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
858 const int nb_comp=transporte.line_size();
859 if (nb_comp>3) Process::exit("EF_stab is not coded for more than 3 components for array transporte.");
860
861 //Pour le limiteur
862
863 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
864 CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
865 CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
866 CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
867 CDoubleArrView alpha_tab = alpha_tab_.view_ro();
868 CDoubleArrView beta = beta_.view_ro();
869 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
870 DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
871 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
872 range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
873 const int elem, const int facei_loc)
874 {
875 double P_plus[3],P_moins[3],Q_plus[3],Q_moins[3];
876 int facei = elem_faces(elem, facei_loc);
877 calculer_senseur(Kij, transporteV, nb_comp, facei, elem_faces, face_voisins, num_fac_loc, P_plus, P_moins,
878 Q_plus, Q_moins);
879
880 for (int facej_loc = 0; facej_loc < nb_faces_elem; facej_loc++)
881 if (facej_loc != facei_loc)
882 {
883 int facej = elem_faces(elem, facej_loc);
884
885 double kij = Kij(elem, facei_loc, facej_loc);
886 double kji = Kij(elem, facej_loc, facei_loc);
887 double dij = Dij(elem, facei_loc, facej_loc, Kij);
888 double lij = kij + dij;
889 double lji = kji + dij;
890 assert(lij >= 0);
891 assert(lji >= 0);
892
893 if (lij <= lji) //facei est amont
894 {
895 int face_amont = facei;
896 int face_aval = facej;
897
898 //Si lij==lji, on passe deux foix dans la boucle
899 //d'ou la presence du coefficient 1/2
900 double coeff = 1. * (lij < lji) + 0.5 * (lij == lji);
901 assert(coeff == 1. || coeff == 0.5);
902
903 // Registers for performance:
904 double alpha_beta_amont = alpha_tab[face_amont] * beta[face_amont];
905 double alpha_beta_aval = alpha_tab[face_aval] * beta[face_aval];
906
907 for (int dim = 0; dim < nb_comp; dim++)
908 {
909 int ligne = face_aval * nb_comp + dim;
910 int colonne = face_amont * nb_comp + dim;
911
912 double delta = transporteV[colonne] - transporteV[ligne];
913
914 //Limiteur de pente
915 double R;
916 if (delta >= 0.) R = (Kokkos::fabs(P_plus[dim]) < DMINFLOAT) ? 0. : Q_plus[dim] / P_plus[dim];
917 else R = (Kokkos::fabs(P_moins[dim]) < DMINFLOAT) ? 0. : Q_moins[dim] / P_moins[dim];
918
919 double limit = limiteur(R);
920 //double daij = minimum(limit * dij, lji);
921 double daij = Kokkos::fmin(limit * dij, lji);
922 assert(daij >= 0);
923 assert(daij <= lji);
924
925 double coeffij = alpha_beta_amont * daij * coeff * delta;
926 double coeffji = alpha_beta_aval * daij * coeff * delta;
927
928 //Calcul de resu
929 Kokkos::atomic_add(&resuV[colonne], + coeffij);
930 Kokkos::atomic_add(&resuV[ligne], - coeffji);
931 }
932 }
933 }
934 });
935 end_gpu_timer(__KERNEL_NAME__);
936 return resu;
937}
938
939//ATTENTION : suppose les parametres P_plus, P_moins, Q_plus, Q_moins nuls en entree
940KOKKOS_INLINE_FUNCTION void
941Op_Conv_EF_VEF_P1NC_Stab::calculer_senseur(CDoubleTabView3 Kij, CDoubleArrView transporteV,
942 const int nb_comp, const int face_i,
943 CIntTabView elem_faces, CIntTabView face_voisins, CIntTabView num_fac_loc,
944 double* P_plus, double* P_moins,
945 double* Q_plus, double* Q_moins) const
946{
947 for (int i = 0; i < nb_comp; i++)
948 {
949 P_plus[i] = 0., P_moins[i] = 0.;
950 Q_plus[i] = 0., Q_moins[i] = 0.;
951 }
952 const int nb_faces_elem=(int)elem_faces.extent(1);
953 for (int elem_voisin=0; elem_voisin<2; elem_voisin++)
954 {
955 int elem = face_voisins(face_i,elem_voisin);
956 if (elem!=-1)
957 {
958 int face_i_loc = num_fac_loc(face_i,elem_voisin);
959 //On travaille sur les faces de "elem"
960 for (int face_k_loc=0; face_k_loc<nb_faces_elem; face_k_loc++)
961 {
962 int face_k=elem_faces(elem,face_k_loc);
963 double kik=Kij(elem,face_i_loc,face_k_loc);
964 //
965 //Calcul des variables intermediaires
966 //
967 for (int dim=0; dim<nb_comp; dim++)
968 {
969 double deltaki = transporteV[face_k*nb_comp+dim]-transporteV[face_i*nb_comp+dim];
970 //Calcul de P_plus et P_moins
971 //Calcul de Q_plus et Q_moins
972 /* Codage initial:
973 P_plus(dim)+=minimum(0.,kik)*minimum(0.,deltaki);
974 P_moins(dim)+=minimum(0.,kik)*maximum(0.,deltaki);
975 Q_plus(dim)+=maximum(0.,kik)*maximum(0.,deltaki);
976 Q_moins(dim)+=maximum(0.,kik)*minimum(0.,deltaki);
977 */
978 // Codage optimise:
979 double tmp = kik*deltaki;
980 if (kik>0)
981 {
982 if (tmp>0) Q_plus[dim]+=tmp;
983 else Q_moins[dim]+=tmp;
984 }
985 else
986 {
987 if (tmp>0) P_plus[dim]+=tmp;
988 else P_moins[dim]+=tmp;
989 }
990 }//fin du for sur "dim"
991 //
992 //Fin du calcul des variables intermediaires
993 //
994 }//fin du for sur "face_k_loc"
995 }//fin du if sur "elem!=-1"
996 }//fin du for sur "elem_voisin"
997}
998
999//ATTENTION : suppose les parametres P_plus, P_moins, Q_plus, Q_moins nuls en entree
1000inline void
1001Op_Conv_EF_VEF_P1NC_Stab::calculer_senseur(const DoubleTab& Kij, const DoubleVect& transporteV,
1002 const int nb_comp, const int face_i,
1003 const IntTab& elem_faces, const IntTab& face_voisins, const IntTab& num_fac_loc,
1004 ArrOfDouble& P_plus, ArrOfDouble& P_moins,
1005 ArrOfDouble& Q_plus, ArrOfDouble& Q_moins) const
1006{
1007 assert(P_plus.size_array()==nb_comp);
1008 assert(Q_plus.size_array()==nb_comp);
1009 assert(P_moins.size_array()==nb_comp);
1010 assert(Q_moins.size_array()==nb_comp);
1011 const int nb_faces_elem=elem_faces.dimension(1);
1012 for (int elem_voisin=0; elem_voisin<2; elem_voisin++)
1013 {
1014 int elem = face_voisins(face_i,elem_voisin);
1015 if (elem!=-1)
1016 {
1017 int face_i_loc = num_fac_loc(face_i,elem_voisin);
1018 assert(face_i_loc>=0);
1019 assert(face_i_loc<nb_faces_elem);
1020 //On travaille sur les faces de "elem"
1021 for (int face_k_loc=0; face_k_loc<nb_faces_elem; face_k_loc++)
1022 {
1023 int face_k=elem_faces(elem,face_k_loc);
1024 double kik=Kij(elem,face_i_loc,face_k_loc);
1025 //
1026 //Calcul des variables intermediaires
1027 //
1028 for (int dim=0; dim<nb_comp; dim++)
1029 {
1030 double deltaki = transporteV[face_k*nb_comp+dim]-transporteV[face_i*nb_comp+dim];
1031 //Calcul de P_plus et P_moins
1032 //Calcul de Q_plus et Q_moins
1033 /* Codage initial:
1034 P_plus(dim)+=minimum(0.,kik)*minimum(0.,deltaki);
1035 P_moins(dim)+=minimum(0.,kik)*maximum(0.,deltaki);
1036 Q_plus(dim)+=maximum(0.,kik)*maximum(0.,deltaki);
1037 Q_moins(dim)+=maximum(0.,kik)*minimum(0.,deltaki);
1038 */
1039 // Codage optimise:
1040 double tmp = kik*deltaki;
1041 if (kik>0)
1042 {
1043 if (tmp>0) Q_plus[dim]+=tmp;
1044 else Q_moins[dim]+=tmp;
1045 }
1046 else
1047 {
1048 if (tmp>0) P_plus[dim]+=tmp;
1049 else P_moins[dim]+=tmp;
1050 }
1051 assert(P_plus[dim]>=0);
1052 assert(Q_plus[dim]>=0);
1053 assert(P_moins[dim]<=0);
1054 assert(Q_moins[dim]<=0);
1055 }//fin du for sur "dim"
1056 //
1057 //Fin du calcul des variables intermediaires
1058 //
1059 }//fin du for sur "face_k_loc"
1060 }//fin du if sur "elem!=-1"
1061 }//fin du for sur "elem_voisin"
1062}
1063
1064void Op_Conv_EF_VEF_P1NC_Stab::test(const DoubleTab& transporte, const DoubleTab& resu, const DoubleTab& tab_vitesse) const
1065{
1066 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1067 const IntTab& elem_faces = domaine_VEF.elem_faces();
1068 const int nb_faces_elem=elem_faces.dimension(1);
1069 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
1070
1071 DoubleTab Kij2(nb_elem_tot,nb_faces_elem, nb_faces_elem);
1072 Kij2=0.;
1073
1074 DoubleTab Kij_ancien(nb_elem_tot,nb_faces_elem, nb_faces_elem);
1075 Kij_ancien=0.;
1076
1077 test_difference_Kij(transporte,Kij2,Kij_ancien,tab_vitesse);
1078 test_difference_resu(Kij2,Kij_ancien,transporte,resu,tab_vitesse);
1079}
1080
1081void Op_Conv_EF_VEF_P1NC_Stab::test_difference_Kij(const DoubleTab& transporte, DoubleTab& Kij, DoubleTab& Kij_ancien, const DoubleTab& tab_vitesse ) const
1082{
1083 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1084 const IntTab& elem_faces = domaine_VEF.elem_faces();
1085 const IntTab& face_voisins = domaine_VEF.face_voisins();
1086 const DoubleTab& face_normales=domaine_VEF.face_normales();
1087 const int nb_faces_elem=elem_faces.dimension(1);
1088 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
1089 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1090
1091 const int nb_comp = (transporte.nb_dim()!=1) ? transporte.dimension(1) : 1;
1092
1093 calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
1094
1095 //
1096 // Calcul des Kij_ancien :
1097 //
1098 for(int elem0=0; elem0<nb_elem_tot; elem0++)
1099 {
1100 int face_locj=0;
1101 for(int face_loci=0; face_loci<nb_faces_elem; face_loci++)
1102 {
1103 int face_i0=elem_faces(elem0,face_loci);
1104 double signei=1.0;
1105 if(face_voisins(face_i0,0)!=elem0)
1106 signei=-1.0;
1107 double psci=0;
1108 for(int comp=0; comp<dimension; comp++)
1109 psci+=tab_vitesse(face_i0,comp)*face_normales(face_i0,comp);
1110 psci*=signei;
1111 //Kij_ancien(elem,face_loci,face_loci)=0.;
1112 for(face_locj=face_loci+1; face_locj<nb_faces_elem; face_locj++)
1113 {
1114 int face_j0=elem_faces(elem0,face_locj);
1115 double signej=1.0;
1116 if(face_voisins(face_j0,0)!=elem0)
1117 signej=-1.0;
1118
1119 double pscj=0;
1120 //psci=0;
1121 for(int comp=0; comp<dimension; comp++)
1122 pscj+=tab_vitesse(face_j0,comp)*face_normales(face_j0,comp);
1123 pscj*=signej;
1124 Kij_ancien(elem0,face_loci,face_locj)=-1./nb_faces_elem*pscj;
1125 Kij_ancien(elem0,face_loci,face_loci)+=1./nb_faces_elem*pscj;
1126 Kij_ancien(elem0,face_locj,face_loci)=-1./nb_faces_elem*psci;
1127 Kij_ancien(elem0,face_locj,face_locj)+=1./nb_faces_elem*psci;
1128 }
1129 }
1130 }
1131 // Correction des Kij_ancien pour Dirichlet !
1132 {
1133 int nb_bord=domaine_Cl_VEF.nb_cond_lim();
1134 double coeff=1./dimension;
1135 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1136 {
1137 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1138 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1139 int nb_faces_tot=le_bord.nb_faces_tot();
1140
1141 if ( (sub_type(Dirichlet,la_cl.valeur()))
1142 || (sub_type(Dirichlet_homogene,la_cl.valeur()))
1143 )
1144 {
1145 if ( volumes_etendus_ )
1146 {
1147 for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
1148 {
1149 int face = le_bord.num_face(ind_face);
1150 int elem=face_voisins(face,0);
1151 assert(elem!=-1);
1152 int face_loc_j;
1153 int face_j=-1;
1154 for (face_loc_j=0; (face_loc_j<nb_faces_elem && face_j!=face); face_loc_j++)
1155 {
1156 face_j=elem_faces(elem,face_loc_j);
1157 }
1158 face_loc_j--;
1159 assert(face_loc_j>=0);
1160 assert(face_loc_j<nb_faces_elem);
1161 assert(elem_faces(elem,face_loc_j)==face);
1162 const double kjj=Kij_ancien(elem,face_loc_j,face_loc_j);
1163 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
1164 {
1165 int face_i=elem_faces(elem,face_loc_i);
1166 if(face_i!=face)
1167 {
1168 double& kii=Kij_ancien(elem,face_loc_i,face_loc_i);
1169 const double kji=Kij_ancien(elem,face_loc_j,face_loc_i);
1170 kii+=coeff*kji;
1171 double& kij=Kij_ancien(elem,face_loc_i,face_loc_j);
1172 kij+=coeff*kjj;
1173 for (int face_loc_k=(face_loc_i+1); face_loc_k<nb_faces_elem; face_loc_k++)
1174 {
1175 int face_k=elem_faces(elem,face_loc_k);
1176 if(face_k!=face)
1177 {
1178 double& kik=Kij_ancien(elem,face_loc_i,face_loc_k);
1179 const double kjk=Kij_ancien(elem,face_loc_j,face_loc_k);
1180 double& kki=Kij_ancien(elem,face_loc_k,face_loc_i);
1181 kik+=coeff*kjk;
1182 kki+=coeff*kji;
1183 }
1184 }
1185 }
1186 }
1187 {
1188 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
1189 Kij_ancien(elem,face_loc_j,face_loc_i)=0;
1190 }
1191 {
1192 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
1193 {
1194 double sum=0.;
1195 for (int face_loc_k=0; face_loc_k<nb_faces_elem; face_loc_k++)
1196 {
1197 sum+=Kij_ancien(elem,face_loc_i,face_loc_k);
1198 }
1199 //Cerr << "somme apres : " << sum << finl;
1200 Kij_ancien(elem,face_loc_i,face_loc_i)-=sum;//car div(u)=0!
1201 }
1202 }
1203 }// for face
1204
1205 }//fin du if sur "volumes_etendus_"
1206
1207 }// sub_type Dirichlet
1208 }
1209 }
1210
1211 Kij-=Kij_ancien;
1212 const double max_kij = local_max_abs_vect(Kij);
1213 if (max_kij > 1.e-15)
1214 {
1215 Cerr << "Erreur dans le calcul des Kij : " << max_kij << finl;
1216 Cerr << "Sortie du programme" << finl;
1217 Process::exit();
1218 }
1219
1220 Kij+=Kij_ancien;
1221}
1222
1223void Op_Conv_EF_VEF_P1NC_Stab::test_difference_resu(const DoubleTab& Kij, const DoubleTab& Kij_ancien,
1224 const DoubleTab& transporte,const DoubleTab& resu, const DoubleTab& tab_vitesse) const
1225{
1226 DoubleTab resu1(resu);
1227 resu1=0;
1228
1229 if (is_compressible_) ajouter_partie_compressible(transporte,resu1,tab_vitesse);
1230 ajouter_operateur_centre(Kij,transporte,resu1);
1231 ajouter_diffusion(Kij,transporte,resu1);
1232 ajouter_antidiffusion(Kij,transporte,resu1);
1233 mettre_a_jour_pour_periodicite(resu1);
1234
1235 DoubleTab resu2(resu);
1236 resu2=0;
1237
1238 //
1239 // Calcul de resu2
1240 //
1241
1242 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1243
1244 const IntTab& elem_faces = domaine_VEF.elem_faces();
1245 const IntTab& face_voisins = domaine_VEF.face_voisins();
1246 const int nb_faces_elem=elem_faces.line_size();
1247 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
1248 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1249 const int nb_comp=resu.line_size();
1250 const int nb_faces0 = transporte.dimension(0);
1251
1252 ArrOfDouble pplusi(nb_comp);
1253 ArrOfDouble qplusi(nb_comp);
1254 ArrOfDouble pmoinsi(nb_comp);
1255 ArrOfDouble qmoinsi(nb_comp);
1256
1257 for(int elem=0; elem<nb_elem_tot; elem++)
1258 {
1259 int face_locj=0;
1260
1261 for(int face_loci=0; face_loci<nb_faces_elem; face_loci++)
1262 {
1263 int face_i0=elem_faces(elem,face_loci);
1264 for(int comp0=0; comp0<nb_comp; comp0++)
1265 resu2(face_i0,comp0)+=Kij_ancien(elem,face_loci,face_loci)*transporte(face_i0,comp0);
1266
1267 pplusi=0.;
1268 qplusi=0.;
1269 pmoinsi=0.;
1270 qmoinsi=0.;
1271
1272 int elem1=face_voisins(face_i0,0), elem2=face_voisins(face_i0,1);
1273 int face_loc_i=0, face_loc_j=0;
1274 double dT,min_dT,max_dT, K,min_K,max_K;
1275
1276 //
1277 // dans elem1 :
1278 //
1279 while((face_loc_i<nb_faces_elem)&&(elem_faces(elem1,face_loc_i)!=face_i0))
1280 face_loc_i++;
1281 if(face_loc_i==nb_faces_elem)
1282 {
1283 //Periodique!!
1284 assert(elem2!=-1);
1285 face_loc_i=0;
1286 while( (face_loc_i<nb_faces_elem) &&
1287 (face_voisins(elem_faces(elem1,face_loc_i),0)!=elem2) &&
1288 (face_voisins(elem_faces(elem1,face_loc_i),1)!=elem2) )
1289 face_loc_i++;
1290 }
1291
1292 assert(face_loc_i<nb_faces_elem);
1293 for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
1294 {
1295 int face_j=elem_faces(elem1,face_loc_j);
1296 if(face_j!=face_i0)
1297 {
1298 K=Kij_ancien(elem1,face_loc_i,face_loc_j);
1299
1300 if(K>0.)
1301 {
1302 max_K=K ;
1303 min_K=0.;
1304 }
1305 else
1306 {
1307 max_K=0.;
1308 min_K=K ;
1309 }
1310
1311 for(int comp0=0; comp0<nb_comp; comp0++)
1312 {
1313 dT =transporte(face_j,comp0);
1314 dT-=transporte(face_i0,comp0);
1315
1316 if(dT>0.)
1317 {
1318 max_dT=dT ;
1319 min_dT=0 ;
1320 }
1321 else
1322 {
1323 max_dT=0. ;
1324 min_dT=dT ;
1325 }
1326
1327 pplusi[comp0] +=min_K*min_dT;
1328 pmoinsi[comp0]+=min_K*max_dT;
1329 qplusi[comp0] +=max_K*max_dT;
1330 qmoinsi[comp0]+=max_K*min_dT;
1331 }
1332 }
1333 }
1334
1335 //
1336 // dans elem2 :
1337 //
1338
1339 if(elem2!=-1)
1340 {
1341 face_loc_i=0;
1342 while((face_loc_i<nb_faces_elem)&&(elem_faces(elem2,face_loc_i)!=face_i0))
1343 face_loc_i++;
1344 if(face_loc_i==nb_faces_elem)
1345 {
1346 //Periodique!!
1347 face_loc_i=0;
1348 while( (face_loc_i<nb_faces_elem) &&
1349 (face_voisins(elem_faces(elem2,face_loc_i),0)!=elem1) &&
1350 (face_voisins(elem_faces(elem2,face_loc_i),1)!=elem1) )
1351 face_loc_i++;
1352 }
1353 assert(face_loc_i<nb_faces_elem);
1354 for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
1355 {
1356 int face_j=elem_faces(elem2,face_loc_j);
1357 if(face_j!=face_i0)
1358 {
1359 K=Kij_ancien(elem2,face_loc_i,face_loc_j);
1360
1361 if(K>0.)
1362 {
1363 max_K=K ;
1364 min_K=0.;
1365 }
1366 else
1367 {
1368 max_K=0.;
1369 min_K=K ;
1370 }
1371
1372 for(int comp0=0; comp0<nb_comp; comp0++)
1373 {
1374 dT =transporte(face_j,comp0);
1375 dT-=transporte(face_i0,comp0);
1376
1377 if(dT>0.)
1378 {
1379 max_dT=dT ;
1380 min_dT=0 ;
1381 }
1382 else
1383 {
1384 max_dT=0. ;
1385 min_dT=dT ;
1386 }
1387
1388 pplusi[comp0] +=min_K*min_dT;
1389 pmoinsi[comp0]+=min_K*max_dT;
1390 qplusi[comp0] +=max_K*max_dT;
1391 qmoinsi[comp0]+=max_K*min_dT;
1392 }
1393 }
1394 }
1395 }
1396
1397 for(face_locj=0; face_locj<nb_faces_elem; face_locj++)
1398 if(face_locj!=face_loci)
1399 {
1400 int face_j0=elem_faces(elem,face_locj);
1401 const double kij=Kij_ancien(elem,face_loci,face_locj);
1402 const double kji=Kij_ancien(elem,face_locj,face_loci);
1403 double dij=Dij(elem,face_loci,face_locj,Kij_ancien);
1404 double lji=kji+dij;
1405 double lij=kij+dij;
1406 assert(lij>=0);
1407 assert(lji>=0);
1408
1409 for(int comp0=0; comp0<nb_comp; comp0++)
1410 {
1411 const double Ti=transporte(face_i0,comp0);
1412 const double Tj=transporte(face_j0,comp0);
1413 double deltaij=Ti-Tj;
1414 double Fij=0;
1415 if(lij<=lji)
1416 {
1417 double coef=1;
1418 if (lij==lji) coef=.5;
1419 if(deltaij)
1420 {
1421 if(Ti >= Tj)
1422 {
1423 if(pplusi[comp0])
1424 {
1425 double R=qplusi[comp0]/pplusi[comp0];
1426 Fij=minimum(limiteur(R)*dij,lji);
1427 }
1428 }
1429 else if(pmoinsi[comp0])
1430 {
1431 double R=qmoinsi[comp0]/pmoinsi[comp0];
1432 Fij=minimum(limiteur(R)*dij,lji);
1433 }
1434
1435 assert(Fij*dij>=0);
1436 Fij-=dij;
1437 Fij*=deltaij;
1438 }
1439 resu2(face_i0,comp0)+=coef*(kij*Tj+Fij);
1440 resu2(face_j0,comp0)+=coef*(kji*Ti-Fij);
1441 }
1442 }
1443 }
1444 }
1445 }
1446
1447 // Pour periodique
1448
1449 int nb_bord=domaine_Cl_VEF.nb_cond_lim();
1450 int face;
1451 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1452 {
1453 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1454 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1455 int num1 = le_bord.num_premiere_face();
1456 int nb_faces_b=le_bord.nb_faces();
1457 int num2 = num1 + nb_faces_b;
1458 if (sub_type(Periodique,la_cl.valeur()))
1459 {
1460 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
1461 int face_associee;
1462 IntVect fait(nb_faces_b);
1463 fait = 0;
1464
1465 for (face=num1; face<num2; face++)
1466 {
1467 if (fait(face-num1) == 0)
1468 {
1469 fait(face-num1) = 1;
1470 face_associee=la_cl_perio.face_associee(face-num1);
1471 fait(face_associee) = 1;
1472 for (int comp=0; comp<nb_comp; comp++)
1473 resu2(face_associee+num1, comp)=(resu2(face,comp)+=resu2(face_associee+num1,comp));
1474 }// if fait
1475 }// for face
1476 }// sub_type Perio
1477 }
1478
1479 resu1-=resu2;
1480
1481 //Dans le cas des faces de Dirichlet, notre algorithme ne calcule aucune
1482 //valeur sur les faces de Dirichlet car de toute facon, elles sont ecrasees
1483 //par la matrice de masse. Mais ce n'est pas le cas de l'ancien algorithme
1484 //d'ou la modification apportee ici pour eviter de ne voir que les erreurs
1485 //commises sur les bords de Dirichlet qui n'ont aucune importance
1486 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1487 {
1488 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1489 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1490 int num1 = le_bord.num_premiere_face();
1491 int nb_faces=le_bord.nb_faces();
1492 int num2 = num1 + nb_faces;
1493
1494 if (sub_type(Dirichlet,la_cl.valeur()) || sub_type(Dirichlet_homogene,la_cl.valeur()))
1495 {
1496 for (face=num1; face<num2; face++)
1497 for (int dim=0; dim<nb_comp; dim++)
1498 resu1(face,dim)=0;
1499 }//fin du if sur Dirichlet
1500 }//fin du for sur "n_bord"
1501
1502 const double max_abs_resu1 = local_max_abs_vect(resu1);
1503 Journal() << "local_max_abs_vect(resu1) = " << max_abs_resu1
1504 << " " << equation().schema_temps().temps_courant() << finl;
1505
1506 if (max_abs_resu1 > 1.e-15)
1507 {
1508 Cerr << "Erreur dans le calcul des resu : " << max_abs_resu1 << finl;
1509 Cerr << "Affichage des faces de bord" << finl;
1510
1511 /**************************************************/
1512 Cerr << "Affichage des faces de bord." << finl;
1513
1514 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1515 {
1516 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1517 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1518 int num1 = le_bord.num_premiere_face();
1519 int nb_faces=le_bord.nb_faces();
1520 int num2 = num1 + nb_faces;
1521
1522 if (sub_type(Periodique,la_cl.valeur()))
1523 {
1524 Cerr << "Bord periodique : ";
1525 for (face=num1; face<num2; face++)
1526 {
1527 Cerr << face << ",";
1528 }// for face
1529 Cerr << finl;
1530 }// sub_type Perio
1531
1532 else if (sub_type(Dirichlet,la_cl.valeur()) || (sub_type(Dirichlet_homogene,la_cl.valeur())) )
1533 {
1534 Cerr << "Bord de Dirichlet : ";
1535 for (face=num1; face<num2; face++)
1536 {
1537 Cerr << face << ",";
1538 }// for face
1539 Cerr << finl;
1540 }
1541
1542 }//fin du for sur "nbord"
1543 /**************************************************/
1544
1545 /**************************************************/
1546 Cerr << "Affichage des faces qui posent probleme : " << finl;
1547 if (nb_comp==1)
1548 {
1549 for (int face_i=0; face_i<nb_faces0; face_i++)
1550 {
1551 if (resu1(face_i)>1.e-15)
1552 Cerr << face_i << "(" << face_voisins(face_i,0) << ","
1553 << face_voisins(face_i,1) << ") ; ";
1554 }//fin du for sur "face_i"
1555 }
1556 else
1557 {
1558 for (int face_i=0; face_i<nb_faces0; face_i++)
1559 {
1560 Cerr << face_i << "(" << face_voisins(face_i,0) << ","
1561 << face_voisins(face_i,1) << ") ";
1562
1563 for (int dim=0; dim<nb_comp; dim++)
1564 {
1565 Cerr << ", resu1("
1566 << face_i << "," << dim << ")= "
1567 << resu1(face_i,dim);
1568 }
1569 Cerr << finl;
1570 }
1571 }
1572 Cerr << finl;
1573 /**************************************************/
1574
1575 /**************************************************/
1576 resu1+=resu2;
1577
1578 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1579 {
1580 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1581 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1582 int num1 = le_bord.num_premiere_face();
1583 int nb_faces=le_bord.nb_faces();
1584 int num2 = num1 + nb_faces;
1585
1586 if (sub_type(Periodique,la_cl.valeur()))
1587 {
1588 Cerr << "Affichage des valeurs des resu au bord perio" << finl;
1589 for (face=num1; face<num2; face++)
1590 {
1591 Cerr << "resu1(" << face << ") : " << resu1(face) << finl;
1592 Cerr << "resu2(" << face << ") : " << resu2(face) << finl;
1593 }// for face
1594 Cerr << finl;
1595 }// sub_type Perio
1596
1597 }//fin du for sur "nbord"
1598
1599 /**************************************************/
1600
1601 static int count = 0;
1602 count++;
1603 if (count==2)
1604 {
1605 Cerr << "Sortie du programme" << finl;
1606 Process::exit();
1607 }
1608 }
1609}
1610
1611void Op_Conv_EF_VEF_P1NC_Stab::mettre_a_jour_pour_periodicite(DoubleTab& resu) const
1612{
1613 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1614 const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
1615 const int nb_comp = (resu.nb_dim()==1) ? 1 : resu.dimension(1);
1616
1617 //Faces de bord
1618 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1619 {
1620 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1621 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1622 int num1=0;
1623 int num2=le_bord.nb_faces();
1624
1625 if (sub_type(Periodique,la_cl.valeur()))
1626 {
1627 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
1628 CIntArrView le_bord_num_face = static_cast<const ArrOfInt&>(le_bord.num_face()).view_ro();
1629 CIntArrView face_associee = static_cast<const ArrOfInt&>(la_cl_perio.face_associee()).view_ro();
1630 DoubleArrView resuV = static_cast<ArrOfDouble&>(resu).view_rw();
1631 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(num1, num2), KOKKOS_LAMBDA(const int ind_face)
1632 {
1633 int facei = le_bord_num_face(ind_face);
1634 int ind_face_associee = face_associee(ind_face);
1635 int faceiAss = le_bord_num_face(ind_face_associee);
1636
1637 if (facei<faceiAss)
1638 for (int dim=0; dim<nb_comp; dim++)
1639 {
1640 int ligne=facei*nb_comp+dim;
1641 int ligneAss=faceiAss*nb_comp+dim;
1642
1643 Kokkos::atomic_add(&resuV[ligneAss],resuV[ligne]);
1644 Kokkos::atomic_store(&resuV[ligne],resuV[ligneAss]);
1645 }
1646
1647 });//fin du for sur "face_i"
1648 end_gpu_timer(__KERNEL_NAME__);
1649 }//fin du if sur "Periodique"
1650
1651 }//fin du for sur "n_bord"
1652}
1653
1654void Op_Conv_EF_VEF_P1NC_Stab::ajouter_old(const DoubleTab& transporte, DoubleTab& resu, const DoubleTab& tab_vitesse ) const
1655{
1656 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1657 // const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
1658 // const DoubleTab& vitesse=la_vitesse.valeurs();
1659
1660 const IntTab& elem_faces = domaine_VEF.elem_faces();
1661 const IntTab& face_voisins = domaine_VEF.face_voisins();
1662 const DoubleTab& face_normales=domaine_VEF.face_normales();
1663 const int nb_faces_elem=elem_faces.dimension(1);
1664 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
1665 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1666
1667 assert(nb_faces_elem==(dimension+1));
1668
1669 int nb_comp=resu.line_size();
1670
1671 int face_i0, face_j0, elem0, comp0;
1672 DoubleTab Kij(nb_elem_tot,nb_faces_elem, nb_faces_elem);
1673 //
1674 // Calcul des Kij :
1675 //
1676 for(elem0=0; elem0<nb_elem_tot; elem0++)
1677 {
1678 int face_loci=0;
1679 int face_locj=0;
1680 for(; face_loci<nb_faces_elem; face_loci++)
1681 {
1682 face_i0=elem_faces(elem0,face_loci);
1683 double signei=1.0;
1684 if(face_voisins(face_i0,0)!=elem0)
1685 signei=-1.0;
1686 double psci=0;
1687 for(comp0=0; comp0<dimension; comp0++)
1688 psci+=tab_vitesse(face_i0,comp0)*face_normales(face_i0,comp0);
1689 psci*=signei;
1690 //Kij(elem,face_loci,face_loci)=0.;
1691 for(face_locj=face_loci+1; face_locj<nb_faces_elem; face_locj++)
1692 {
1693 face_j0=elem_faces(elem0,face_locj);
1694 double signej=1.0;
1695 if(face_voisins(face_j0,0)!=elem0)
1696 signej=-1.0;
1697
1698 double pscj=0;
1699 //psci=0;
1700 for(comp0=0; comp0<dimension; comp0++)
1701 pscj+=tab_vitesse(face_j0,comp0)*face_normales(face_j0,comp0);
1702 pscj*=signej;
1703 Kij(elem0,face_loci,face_locj)=-1./nb_faces_elem*pscj;
1704 Kij(elem0,face_loci,face_loci)+=1./nb_faces_elem*pscj;
1705 Kij(elem0,face_locj,face_loci)=-1./nb_faces_elem*psci;
1706 Kij(elem0,face_locj,face_locj)+=1./nb_faces_elem*psci;
1707 }
1708 }
1709 }
1710 //
1711 // Correction des Kij pour Dirichlet !
1712 //
1713 {
1714 int nb_bord=domaine_Cl_VEF.nb_cond_lim();
1715 int face;
1716 double coeff=1./dimension;
1717 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1718 {
1719 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1720 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1721 int nb_faces_tot=le_bord.nb_faces_tot();
1722 if (( (sub_type(Dirichlet,la_cl.valeur())) || (sub_type(Dirichlet_homogene,la_cl.valeur())) )
1723 && ( volumes_etendus_))
1724 {
1725 for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
1726 {
1727 face=le_bord.num_face(ind_face);
1728 int elem=face_voisins(face,0);
1729 assert(elem!=-1);
1730 int face_loc_j;
1731 int face_j=-1;
1732 for (face_loc_j=0; (face_loc_j<nb_faces_elem && face_j!=face); face_loc_j++)
1733 {
1734 face_j=elem_faces(elem,face_loc_j);
1735 }
1736 face_loc_j--;
1737 assert(face_loc_j>=0);
1738 assert(face_loc_j<nb_faces_elem);
1739 assert(elem_faces(elem,face_loc_j)==face);
1740 const double kjj=Kij(elem,face_loc_j,face_loc_j);
1741 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
1742 {
1743 int face_i=elem_faces(elem,face_loc_i);
1744 if(face_i!=face)
1745 {
1746 double& kii=Kij(elem,face_loc_i,face_loc_i);
1747 const double kji=Kij(elem,face_loc_j,face_loc_i);
1748 kii+=coeff*kji;
1749 double& kij=Kij(elem,face_loc_i,face_loc_j);
1750 kij+=coeff*kjj;
1751 for (int face_loc_k=(face_loc_i+1); face_loc_k<nb_faces_elem; face_loc_k++)
1752 {
1753 int face_k=elem_faces(elem,face_loc_k);
1754 if(face_k!=face)
1755 {
1756 double& kik=Kij(elem,face_loc_i,face_loc_k);
1757 const double kjk=Kij(elem,face_loc_j,face_loc_k);
1758 double& kki=Kij(elem,face_loc_k,face_loc_i);
1759 kik+=coeff*kjk;
1760 kki+=coeff*kji;
1761 }
1762 }
1763 }
1764 }
1765 {
1766 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
1767 Kij(elem,face_loc_j,face_loc_i)=0;
1768 }
1769 {
1770 for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
1771 {
1772 double sum=0.;
1773 for (int face_loc_k=0; face_loc_k<nb_faces_elem; face_loc_k++)
1774 {
1775 sum+=Kij(elem,face_loc_i,face_loc_k);
1776 }
1777 //Cerr << "somme apres : " << sum << finl;
1778 Kij(elem,face_loc_i,face_loc_i)-=sum;//car div(u)=0!
1779 }
1780 }
1781 }// for face
1782 }// sub_type Dirichlet
1783 }
1784 }
1785
1786 //
1787 // Calcul de resu
1788 //
1789
1790 ArrOfDouble pplusi(nb_comp);
1791 ArrOfDouble qplusi(nb_comp);
1792 ArrOfDouble pmoinsi(nb_comp);
1793 ArrOfDouble qmoinsi(nb_comp);
1794
1795 for(elem0=0; elem0<nb_elem_tot; elem0++)
1796 {
1797 int face_loci=0;
1798 int face_locj=0;
1799
1800 for(; face_loci<nb_faces_elem; face_loci++)
1801 {
1802 face_i0=elem_faces(elem0,face_loci);
1803
1804 for(comp0=0; comp0<nb_comp; comp0++)
1805 resu(face_i0,comp0)+=Kij(elem0,face_loci,face_loci)*transporte(face_i0,comp0);
1806
1807 pplusi=0., qplusi=0., pmoinsi=0., qmoinsi=0.;
1808
1809 int elem1=face_voisins(face_i0,0);
1810 int elem2=face_voisins(face_i0,1);
1811 int face_loc_i=0;
1812 int face_loc_j=0;
1813 double dT,min_dT,max_dT;
1814 double K,min_K,max_K;
1815
1816 //
1817 // dans elem1 :
1818 //
1819 while((face_loc_i<nb_faces_elem)&&(elem_faces(elem1,face_loc_i)!=face_i0))
1820 face_loc_i++;
1821 if(face_loc_i==nb_faces_elem)
1822 {
1823 //Periodique!!
1824 assert(elem2!=-1);
1825 face_loc_i=0;
1826 while( (face_loc_i<nb_faces_elem) &&
1827 (face_voisins(elem_faces(elem1,face_loc_i),0)!=elem2) &&
1828 (face_voisins(elem_faces(elem1,face_loc_i),1)!=elem2) )
1829 face_loc_i++;
1830 }
1831 assert(face_loc_i<nb_faces_elem);
1832 for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
1833 {
1834 int face_j=elem_faces(elem1,face_loc_j);
1835 if(face_j!=face_i0)
1836 {
1837 K=Kij(elem1,face_loc_i,face_loc_j);
1838
1839 if(K>0.)
1840 {
1841 max_K=K ;
1842 min_K=0.;
1843 }
1844 else
1845 {
1846 max_K=0.;
1847 min_K=K ;
1848 }
1849
1850 for(comp0=0; comp0<nb_comp; comp0++)
1851 {
1852 dT =transporte(face_j,comp0);
1853 dT-=transporte(face_i0,comp0);
1854
1855 if(dT>0.)
1856 {
1857 max_dT=dT ;
1858 min_dT=0 ;
1859 }
1860 else
1861 {
1862 max_dT=0. ;
1863 min_dT=dT ;
1864 }
1865
1866 pplusi[comp0] +=min_K*min_dT;
1867 pmoinsi[comp0]+=min_K*max_dT;
1868 qplusi[comp0] +=max_K*max_dT;
1869 qmoinsi[comp0]+=max_K*min_dT;
1870 }
1871 }
1872 }
1873
1874 //
1875 // dans elem2 :
1876 //
1877
1878 if(elem2!=-1)
1879 {
1880 face_loc_i=0;
1881 while((face_loc_i<nb_faces_elem)&&(elem_faces(elem2,face_loc_i)!=face_i0))
1882 face_loc_i++;
1883 if(face_loc_i==nb_faces_elem)
1884 {
1885 //Periodique!!
1886 face_loc_i=0;
1887 while( (face_loc_i<nb_faces_elem) &&
1888 (face_voisins(elem_faces(elem2,face_loc_i),0)!=elem1) &&
1889 (face_voisins(elem_faces(elem2,face_loc_i),1)!=elem1) )
1890 face_loc_i++;
1891 }
1892 assert(face_loc_i<nb_faces_elem);
1893 for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
1894 {
1895 int face_j=elem_faces(elem2,face_loc_j);
1896 if(face_j!=face_i0)
1897 {
1898 K=Kij(elem2,face_loc_i,face_loc_j);
1899
1900 if(K>0.)
1901 {
1902 max_K=K ;
1903 min_K=0.;
1904 }
1905 else
1906 {
1907 max_K=0.;
1908 min_K=K ;
1909 }
1910
1911 for(comp0=0; comp0<nb_comp; comp0++)
1912 {
1913 dT =transporte(face_j,comp0);
1914 dT-=transporte(face_i0,comp0);
1915
1916 if(dT>0.)
1917 {
1918 max_dT=dT ;
1919 min_dT=0 ;
1920 }
1921 else
1922 {
1923 max_dT=0. ;
1924 min_dT=dT ;
1925 }
1926
1927 pplusi[comp0] +=min_K*min_dT;
1928 pmoinsi[comp0]+=min_K*max_dT;
1929 qplusi[comp0] +=max_K*max_dT;
1930 qmoinsi[comp0]+=max_K*min_dT;
1931 }
1932 }
1933 }
1934 }
1935
1936
1937
1938 for(face_locj=0; face_locj<nb_faces_elem; face_locj++)
1939 if(face_locj!=face_loci)
1940 {
1941 face_j0=elem_faces(elem0,face_locj);
1942 const double kij=Kij(elem0,face_loci,face_locj);
1943 const double kji=Kij(elem0,face_locj,face_loci);
1944 double dij=Dij(elem0,face_loci,face_locj,Kij);
1945 double lji=kji+dij;
1946 double lij=kij+dij;
1947 assert(lij>=0);
1948 assert(lji>=0);
1949
1950 for(comp0=0; comp0<nb_comp; comp0++)
1951 {
1952 const double Ti=transporte(face_i0,comp0);
1953 const double Tj=transporte(face_j0,comp0);
1954 double deltaij=Ti-Tj;
1955 double Fij=0;
1956 if(lij<=lji)
1957 {
1958 double coef=1;
1959 if (lij==lji) coef=.5;
1960 if(deltaij)
1961 {
1962 if(Ti >= Tj)
1963 {
1964 if(pplusi[comp0])
1965 {
1966 double R=qplusi[comp0]/pplusi[comp0];
1967 Fij=minimum(limiteur(R)*dij,lji);
1968 }
1969 }
1970 else if(pmoinsi[comp0])
1971 {
1972 double R=qmoinsi[comp0]/pmoinsi[comp0];
1973 Fij=minimum(limiteur(R)*dij,lji);
1974 }
1975 assert(Fij*dij>=0);
1976 Fij-=dij;
1977 Fij*=deltaij;
1978 }
1979 resu(face_i0,comp0)+=coef*(kij*Tj+Fij);
1980 resu(face_j0,comp0)+=coef*(kji*Ti-Fij);
1981 }
1982 }
1983 }
1984 }
1985 }
1986 //On reprend ici le traitement qui est fait pour la CL de Neumann_sortie_libre dans Op_Conv_VEF_Face
1987 int nb_bord=domaine_Cl_VEF.nb_cond_lim();
1988 int face;
1989 const int ncomp_ch_transporte = transporte.line_size();
1990
1991 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1992 {
1993 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1994 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1995 int nb_faces_tot=le_bord.nb_faces_tot();
1996
1997 double psc;
1998 int num_face,i;
1999
2000 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
2001 {
2002 const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
2003 //const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2004 int num1 = le_bord.num_premiere_face();
2005 int num2 = num1 + le_bord.nb_faces();
2006
2007 for (num_face=num1; num_face<num2; num_face++)
2008 {
2009 psc =0;
2010 for (i=0; i<dimension; i++)
2011 psc += tab_vitesse(num_face,i)*face_normales(num_face,i);
2012 if (psc>0)
2013 {
2014 for (i=0; i<ncomp_ch_transporte; i++)
2015 resu(num_face,i) -= psc*transporte(num_face,i);
2016 }
2017 else
2018 {
2019 for (i=0; i<ncomp_ch_transporte; i++)
2020 resu(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
2021 fluent_(num_face) -= psc;
2022 }
2023 }
2024 }
2025 // Pour periodique
2026 else if (sub_type(Periodique,la_cl.valeur()))
2027 {
2028 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
2029 int face_associee,ind_face_associee;
2030 IntVect fait(nb_faces_tot);
2031 fait = 0;
2032
2033 for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
2034 {
2035 face=le_bord.num_face(ind_face);
2036 if (fait(ind_face) == 0)
2037 {
2038 fait(ind_face) = 1;
2039 ind_face_associee=la_cl_perio.face_associee(ind_face);
2040 fait(ind_face_associee) = 1;
2041 face_associee=le_bord.num_face(ind_face_associee);
2042 for (int comp=0; comp<nb_comp; comp++)
2043 resu(face_associee, comp)=(resu(face,comp)+=resu(face_associee,comp));
2044 }// if fait
2045 }// for face
2046 }// sub_type Perio
2047 }
2048 /*
2049 ArrOfDouble bilan(nb_comp);
2050 BilanQdmVEF::bilan_qdm(resu, domaine_Cl_VEF, bilan);
2051 if(nb_comp==1)
2052 Cout << "Scalaire Bilan Convectif : " << bilan[0] << finl;
2053 else
2054 for (int comp=0; comp<nb_comp; comp++)
2055 Cout << "Vecteur Bilan Convectif " << comp << " : " << bilan[comp] << finl;
2056 bilan=0;
2057 BilanQdmVEF::bilan_energie(resu, transporte, domaine_Cl_VEF, bilan);
2058 if(nb_comp==1)
2059 Cout << "Scalaire Bilan Convectif Energie : " << bilan[0] << finl;
2060 else
2061 for (int comp=0; comp<nb_comp; comp++)
2062 Cout << "Vecteur Bilan Convectif Energie " << comp << " : " << bilan[comp] << finl;
2063 if(nb_comp==1)
2064 {
2065 Cout << "min = " << min(transporte);
2066 Cout << " max = " << std::max(transporte) << finl;
2067 }
2068 // else
2069 // for (int comp=0; comp<nb_comp; comp++)
2070 // {
2071 // Cout << "min("<<comp<<") = " << transporte.min(comp);
2072 // Cout << " std::max("<<comp<<") = " << transporte.max(comp) << finl;
2073 // }
2074 Cout << " Ratio Antidiffusion/Diffusion = " << sigma_fija/sigma_fijd << finl;
2075 */
2076}
2077
2078//Fonction qui initialise les attributs "elem_nb_faces_dirichlet_"
2079//et "elem_faces_dirichlet_"
2080//REMARQUE : "elem_nb_faces_dirichlet_" contient le nombre de faces de Dirichlet
2081//pour chaque element du maillage
2082//REMARQUE : "elem_faces_dirichlet_" le numero global des faces de Dirichlet
2083//contenu dans un element quelconque du maillage
2084void Op_Conv_EF_VEF_P1NC_Stab::calculer_data_pour_dirichlet()
2085{
2086 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
2087 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
2088
2089 const IntTab& face_voisins = domaine_VEF.face_voisins();
2090 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
2091 const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
2092
2093 //Dimensionnement et initialisation des attributs
2094 //REMARQUE : un element ne peut pas avoir plus de
2095 //(dimension) faces de Dirichlet
2096 elem_nb_faces_dirichlet_.resize(nb_elem_tot);
2097 elem_faces_dirichlet_.resize(nb_elem_tot,Objet_U::dimension);
2098 elem_nb_faces_dirichlet_=0;
2099 elem_faces_dirichlet_=-1;
2100 elem_faces_frontiere.dimensionner(nb_bord);
2101
2102 for (int n_bord=0; n_bord<nb_bord; n_bord++)
2103 {
2104 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
2105 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2106 int face, nb_faces_tot=le_bord.nb_faces_tot();
2107
2108 if ( (sub_type(Dirichlet,la_cl.valeur()))
2109 || (sub_type(Dirichlet_homogene,la_cl.valeur()))
2110 )
2111 {
2112 //
2113 //Remplissage des tableaux
2114 //
2115 for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
2116 {
2117 face = le_bord.num_face(ind_face);
2118 const int elem=face_voisins(face,0);
2119 assert(elem!=-1);
2120 elem_faces_frontiere[n_bord].append_array(elem);
2121
2122 elem_nb_faces_dirichlet_(elem)+=1;
2123 assert(elem_nb_faces_dirichlet_(elem)<=Objet_U::dimension);
2124
2125 if (elem_faces_dirichlet_(elem,0)==-1) elem_faces_dirichlet_(elem,0)=face;
2126 else if (elem_faces_dirichlet_(elem,1)==-1) elem_faces_dirichlet_(elem,1)=face;
2127 else if (Objet_U::dimension==3 && elem_faces_dirichlet_(elem,2)==-1) elem_faces_dirichlet_(elem,2)=face;
2128 else
2129 {
2130 Cerr << "Erreur Op_Conv_EF_VEF_P1NC_Stab::calculer_data_pour_dirichlet()" << finl;
2131 Cerr << "L'element numero " << elem << " contient plus de "
2132 << Objet_U::dimension << " faces de Dirichlet" << finl;
2133 Cerr << "Sortie du programme." << finl;
2134 Process::exit();
2135 }
2136 }//fin du for sur "face"
2137 //
2138 //Fin du remplissage des tableaux
2139 //
2140
2141 }//fin du if sur "Dirichlet"
2142 array_trier_retirer_doublons(elem_faces_frontiere[n_bord]);
2143 }//fin du for sur "n_bord"
2144}
2145
2146void Op_Conv_EF_VEF_P1NC_Stab::completer()
2147{
2148 Op_Conv_VEF_Face::completer();
2149 calculer_data_pour_dirichlet();
2150
2151 // int nb_comp=1;
2152 // if (equation().inconnue().valeurs().nb_dim()>1)
2153 // nb_comp=equation().inconnue().valeurs().dimension(1);
2154
2155 // limiteurs_.resize(le_dom_vef->nb_faces_tot(),nb_comp);
2156 // limiteurs_=0.;
2157
2158 alpha_tab_.resize_array(le_dom_vef->nb_faces_tot());
2159 alpha_tab_ = alpha_;
2160 beta_.resize_array(le_dom_vef->nb_faces_tot());
2161 beta_=1.;
2162
2163 if (ssz_alpha)
2164 {
2165 for (int i=0; i<nb_ssz_alpha; i++)
2166 {
2167 OBS_PTR(Sous_domaine_VF) la_ssz;
2168 const Sous_Domaine& le_sous_domaine=equation().probleme().domaine().ss_domaine(noms_ssz_alpha[i]);
2169 const Domaine_dis_base& le_domaine_dis=le_dom_vef.valeur();
2170 bool trouve=false;
2171 for (int ssz=0; ssz<le_domaine_dis.nombre_de_sous_domaines_dis(); ssz++)
2172 {
2173 if (le_domaine_dis.sous_domaine_dis(ssz).sous_domaine().est_egal_a(le_sous_domaine))
2174 {
2175 trouve=true;
2176 la_ssz=ref_cast(Sous_domaine_VF,le_domaine_dis.sous_domaine_dis(ssz));
2177 }
2178 }
2179
2180 if(!trouve)
2181 {
2182 Cerr << "On ne trouve pas le sous_domaine discretise associe a " << noms_ssz_alpha[i] << finl;
2183 Process::exit();
2184 }
2185 const Sous_domaine_VF& ssz=la_ssz.valeur();
2186 int nb_faces = ssz.les_faces().size();
2187
2188 for (int face=0; face<nb_faces; face++)
2189 {
2190 int la_face=ssz.les_faces()[face];
2191 beta_[la_face] = 1.;
2192 alpha_tab_[la_face] = alpha_ssz(i);
2193 }
2194 }
2195 }
2196
2197
2198
2199
2200 if (sous_domaine)
2201 {
2202 sous_domaine=false;
2203 const Sous_Domaine& le_sous_domaine=equation().probleme().domaine().ss_domaine(nom_sous_domaine);
2204 const Domaine_dis_base& le_domaine_dis=le_dom_vef.valeur();
2205 for (int ssz=0; ssz<le_domaine_dis.nombre_de_sous_domaines_dis(); ssz++)
2206 {
2207 if (le_domaine_dis.sous_domaine_dis(ssz).sous_domaine().est_egal_a(le_sous_domaine))
2208 {
2209 sous_domaine=true;
2210 le_sous_domaine_dis=ref_cast(Sous_domaine_VF,le_domaine_dis.sous_domaine_dis(ssz));
2211 }
2212 }
2213
2214 if(!sous_domaine)
2215 {
2216 Cerr << "On ne trouve pas le sous_domaine discretise associe a " << nom_sous_domaine << finl;
2217 Process::exit();
2218 }
2219
2220 const Sous_domaine_VF& ssz=le_sous_domaine_dis.valeur();
2221 int nb_faces = ssz.les_faces().size();
2222
2223 for (int face=0; face<nb_faces; face++)
2224 {
2225 int la_face=ssz.les_faces()[face];
2226 beta_[la_face] = 0.;
2227 alpha_tab_[la_face] = 1.;
2228 }
2229 }
2230
2231}
2232
2233void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution(const DoubleTab& transporte_2, Matrice_Morse& matrice) const
2234{
2235 if (new_jacobienne_==0)
2236 {
2237 Op_Conv_VEF_Face::ajouter_contribution(transporte_2, matrice) ;
2238 return;
2239 }
2240 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
2241
2242 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
2243 const int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
2244 const int nb_comp=transporte_2.line_size();
2245
2246 DoubleTrav Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
2247 Kij=0.;
2248
2249 //
2250 //Pour tenir compte de la porosite
2251 //
2252 const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
2253 const DoubleTab& vitesse_2=la_vitesse.valeurs();
2254 const DoubleVect& porosite_face = equation().milieu().porosite_face();
2255 DoubleTrav transporte_;
2256 DoubleTrav vitesse_face_;
2257
2258 // soit on a transporte=phi*transporte_ et vitesse=vitesse_
2259 // soit transporte=transporte_ et vitesse=phi*vitesse_
2260 // cela depend si on transporte avec phi u ou avec u.
2261 const int marq=phi_u_transportant(equation());
2262 const DoubleTab& transporte=modif_par_porosite_si_flag(transporte_2,transporte_,!marq,porosite_face);
2263 const DoubleTab& tab_vitesse=modif_par_porosite_si_flag(vitesse_2,vitesse_face_,marq,porosite_face);
2264
2265 calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
2266 if (is_compressible_) ajouter_contribution_partie_compressible(transporte,tab_vitesse,matrice);
2267 ajouter_contribution_operateur_centre(Kij,transporte,matrice);
2268 ajouter_contribution_diffusion(Kij,transporte,matrice);
2269
2270 if (test_) test_implicite();
2271}
2272
2273void Op_Conv_EF_VEF_P1NC_Stab::modifier_pour_Cl (Matrice_Morse& matrice, DoubleTab& secmem) const
2274{
2275 Op_Conv_VEF_Face::modifier_pour_Cl(matrice,secmem);
2276}
2277
2278void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_operateur_centre(const DoubleTab& tab_Kij, const DoubleTab& transporte, Matrice_Morse& matrice_morse) const
2279{
2280 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
2281 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
2282
2283
2284 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
2285 const int nb_faces_elem=domaine_VEF.elem_faces().dimension(1);
2286 const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
2287
2288 const int nb_comp=transporte.line_size();
2289
2290 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
2291 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
2292 Matrice_Morse_View matrice;
2293 matrice.set(matrice_morse);
2294 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
2295 range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
2296 const int elem, const int facei_loc)
2297 {
2298 int facei = elem_faces(elem, facei_loc);
2299
2300 for (int facej_loc = facei_loc+1; facej_loc < nb_faces_elem; facej_loc++)
2301 {
2302 int facej = elem_faces(elem, facej_loc);
2303
2304 double kij = Kij(elem, facei_loc, facej_loc);
2305 double kji = Kij(elem, facej_loc, facei_loc);
2306
2307 for (int dim = 0; dim < nb_comp; dim++)
2308 {
2309 int ligne = facei*nb_comp + dim;
2310 int colonne = facej*nb_comp + dim;
2311
2312 //ATTENTION AU SIGNE : ici on code +div(uT)
2313 matrice.atomic_add(ligne, ligne, kij);
2314 matrice.atomic_add(ligne, colonne, -kij);
2315 matrice.atomic_add(colonne, colonne, kji);
2316 matrice.atomic_add(colonne, ligne, -kji);
2317 }
2318 }
2319 });
2320 end_gpu_timer(__KERNEL_NAME__);
2321
2322 //
2323 //Pour la periodicite
2324 //
2325 const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
2326 for (int n_bord=0; n_bord<nb_bord; n_bord++)
2327 {
2328 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
2329 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2330 int num1 = 0;
2331 int num2=le_bord.nb_faces_tot();//et surtout pas nb_faces sinon on oublie certains coefficients
2332
2333 if (sub_type(Periodique,la_cl.valeur()))
2334 {
2335 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
2336 int faceiAss=0,ind_faceiAss=0;
2337 const IntTab& face_voisins = domaine_VEF.face_voisins();
2338 for (int ind_face=num1; ind_face<num2; ind_face++)
2339 {
2340 ind_faceiAss=la_cl_perio.face_associee(ind_face);
2341
2342 int facei=le_bord.num_face(ind_face);
2343 faceiAss=le_bord.num_face(ind_faceiAss);
2344
2345 //Pour ne parcourir qu'une seule fois les faces perio
2346 if (facei<faceiAss)
2347 for (int elem_loc=0; elem_loc<2; elem_loc++)
2348 {
2349 int elem=face_voisins(facei,elem_loc);
2350 assert(elem!=-1);
2351
2352 //Calcul du numero local de la face dans "elem"
2353 int facei_loc=num_fac_loc(facei,elem_loc);
2354 int faceToComplete;
2355 if (facei_loc!=-1)
2356 faceToComplete=faceiAss;
2357 else
2358 {
2359 faceToComplete=facei;
2360 facei_loc=num_fac_loc(faceiAss,elem_loc);
2361 assert(facei_loc!=-1);
2362 }
2363
2364 //Calcul des coefficients de la matrice dus a "elem"
2365 for (int facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
2366 {
2367 int facej=elem_faces(elem,facej_loc);
2368
2369 if (facej_loc!=facei_loc)
2370 {
2371 double kij=Kij(elem,facei_loc,facej_loc);
2372 //double kji=Kij(elem,facej_loc,facei_loc);
2373
2374 for (int dim=0; dim<nb_comp; dim++)
2375 {
2376 int ligne=faceToComplete*nb_comp+dim;
2377 int colonne=facej*nb_comp+dim;
2378
2379 //ATTENTION AU SIGNE : ici on code +div(uT)
2380 matrice_morse(ligne,ligne)+=kij;
2381 matrice_morse(ligne,colonne)-=kij;
2382 }
2383 }
2384 }
2385 }
2386 }
2387 }
2388 }
2389}
2390
2391void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_diffusion(const DoubleTab& tab_Kij, const DoubleTab& transporte, Matrice_Morse& matrice_morse) const
2392{
2393 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
2394 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
2395 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
2396 const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
2397 const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
2398 const int nb_comp=transporte.line_size();
2399
2400 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
2401 CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
2402 CDoubleArrView alpha_tab = static_cast<const ArrOfDouble&>(alpha_tab_).view_ro();
2403 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
2404 Matrice_Morse_View matrice;
2405 matrice.set(matrice_morse);
2406 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
2407 range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
2408 const int elem, const int facei_loc)
2409 {
2410 int facei = elem_faces(elem, facei_loc);
2411
2412 for (int facej_loc = facei_loc+1; facej_loc < nb_faces_elem; facej_loc++)
2413 {
2414 int facej = elem_faces(elem, facej_loc);
2415
2416 double dij = Dij(elem, facei_loc, facej_loc, Kij);
2417
2418 double coeffij = alpha_tab(facei)*dij;
2419 double coeffji = alpha_tab(facej)*dij;
2420
2421 for (int dim = 0; dim < nb_comp; dim++)
2422 {
2423 int ligne = facei*nb_comp + dim;
2424 int colonne = facej*nb_comp + dim;
2425
2426 //ATTENTION AU SIGNE : ici on code +div(uT)
2427 //REMARQUE : on utilise la symetrie de l'operateur
2428 matrice.atomic_add(ligne, ligne, coeffij);
2429 matrice.atomic_add(ligne, colonne, -coeffij);
2430 matrice.atomic_add(colonne, colonne, coeffji);
2431 matrice.atomic_add(colonne, ligne, -coeffji);
2432 }
2433 }
2434 });
2435 end_gpu_timer(__KERNEL_NAME__);
2436
2437 //
2438 //Pour la periodicite
2439 //
2440 const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
2441 for (int n_bord=0; n_bord<nb_bord; n_bord++)
2442 {
2443 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
2444 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2445 int num1 = 0;
2446 int num2=le_bord.nb_faces_tot();//et surtout pas nb_faces() sinon on oublie certains coefficiens
2447
2448 if (sub_type(Periodique,la_cl.valeur()))
2449 {
2450 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
2451 int faceiAss=0,ind_faceiAss=0;
2452
2453 for (int ind_face=num1; ind_face<num2; ind_face++)
2454 {
2455 int facei=le_bord.num_face(ind_face);
2456 ind_faceiAss=la_cl_perio.face_associee(ind_face);
2457 faceiAss=le_bord.num_face(ind_faceiAss);
2458
2459 //Pour ne parcourir qu'une seule fois les faces perio
2460 if (facei<faceiAss)
2461 for (int elem_loc=0; elem_loc<2; elem_loc++)
2462 {
2463 int elem=face_voisins(facei,elem_loc);
2464 assert(elem!=-1);
2465
2466 //Calcul du numero local de la face dans "elem"
2467 int facei_loc=num_fac_loc(facei,elem_loc);
2468 int faceToComplete;
2469 if (facei_loc!=-1)
2470 faceToComplete=faceiAss;
2471 else
2472 {
2473 faceToComplete=facei;
2474 facei_loc=num_fac_loc(faceiAss,elem_loc);
2475 assert(facei_loc!=-1);
2476 }
2477
2478 //Calcul des coefficients de la matrice dus a "elem"
2479 for (int facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
2480 {
2481 int facej=elem_faces(elem,facej_loc);
2482
2483 if (facej_loc!=facei_loc)
2484 {
2485 double dij=Dij(elem,facei_loc,facej_loc,tab_Kij);
2486 assert(dij>=0);
2487
2488 double coeffij=alpha_tab_[faceToComplete]*dij;
2489 //double coeffji=alpha_tab_[facej]*dij;
2490
2491 for (int dim=0; dim<nb_comp; dim++)
2492 {
2493 int ligne=faceToComplete*nb_comp+dim;
2494 int colonne=facej*nb_comp+dim;
2495
2496 //ATTENTION AU SIGNE : ici on code +div(uT)
2497 matrice_morse(ligne,ligne)+=coeffij;
2498 matrice_morse(ligne,colonne)-=coeffij;
2499 }
2500 }
2501 }
2502 }
2503 }
2504 }
2505 }
2506}
2507
2508//Correction pour le poreux : on rajoute la partie en T div(u)
2509//Variable transportee : T
2510//Variable transportante : u
2511//REMARQUE : il ne FAUT SURTOUT PAS utiliser le tableau Kij car par
2512//construction celui-ci est telle que sum_{j} Kij =0 ce qui revient a
2513//imposer une vitesse a divergence nulle par element. Ce qui est
2514//problematique quand on est en compressible
2515void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_partie_compressible(const DoubleTab& transporte, const DoubleTab& vitesse_2,
2516 Matrice_Morse& matrice) const
2517{
2518 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
2519 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
2520 const IntTab& elem_faces=domaine_VEF.elem_faces();
2521 const IntTab& face_voisins = domaine_VEF.face_voisins();
2522 const DoubleTab& face_normales=domaine_VEF.face_normales();
2523 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
2524 const int nb_faces_elem=elem_faces.line_size();
2525 const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
2526
2527 //Pour tenir compte de la porosite
2528 const int marq = phi_u_transportant(equation());
2529 const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
2530 const DoubleVect& porosite_face = equation().milieu().porosite_face();
2531
2532 DoubleTrav tab_vitesse(vitesse_->valeurs());
2533 for (int i=0; i<tab_vitesse.dimension(0); i++)
2534 for (int j=0; j<tab_vitesse.line_size(); j++)
2535 tab_vitesse(i,j)*=porosite_face(i);
2536
2537 const int nb_comp=transporte.line_size();
2538
2539 double (*formule)(int);
2540
2541 if (!volumes_etendus_)
2542 formule= (dimension==2) ? &formule_Id_2D : &formule_Id_3D;
2543 else
2544 formule= (dimension==2) ? &formule_2D : &formule_3D;
2545
2546 ToDo_Kokkos("critical");
2547 for (int elem=0; elem<nb_elem_tot; elem++)
2548 {
2549 //Type de l'element : le nombre de faces de Dirichlet
2550 //qu'il contient
2551 int type_elem=elem_nb_faces_dirichlet_(elem);
2552 double coeff=formule(type_elem);
2553
2554 //Calcul de la divergence par element
2555 double div=0.;
2556 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
2557 {
2558 int facei=elem_faces(elem,facei_loc);
2559 int signe=(face_voisins(facei,0)==elem)? 1.:-1.;
2560
2561 for (int dim=0; dim<dimension; dim++)
2562 div+=signe*face_normales(facei,dim)*tab_vitesse(facei,dim);
2563 }
2564 div*=coeff;
2565 if (!marq) div/=porosite_elem(elem);
2566
2567 //Calcul de la partie compressible
2568 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
2569 {
2570 int facei=elem_faces(elem,facei_loc);
2571
2572 for (int dim=0; dim<nb_comp; dim++)
2573 {
2574 int ligne=facei*nb_comp+dim;
2575 matrice(ligne,ligne)+=div;
2576 }
2577 }
2578 }
2579
2580 //
2581 //Pour la periodicite
2582 //
2583 const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
2584 for (int n_bord=0; n_bord<nb_bord; n_bord++)
2585 {
2586 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
2587 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2588 int num1 = 0;
2589 int num2 = le_bord.nb_faces();//pour ne parcourir que les faces reelles
2590
2591 if (sub_type(Periodique,la_cl.valeur()))
2592 {
2593 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
2594
2595 for (int ind_face=num1; ind_face<num2; ind_face++)
2596 {
2597 int facei = le_bord.num_face(ind_face);
2598 int ind_faceiAss = la_cl_perio.face_associee(ind_face);
2599 int faceiAss = le_bord.num_face(ind_faceiAss);
2600
2601 //Pour ne parcourir qu'une seule fois les faces perio
2602 if (facei<faceiAss)
2603 for (int elem_loc=0; elem_loc<2; elem_loc++)
2604 {
2605 int elem = face_voisins(facei,elem_loc);
2606 assert(elem!=-1);
2607
2608 //Calcul du numero local de la face dans "elem"
2609 int facei_loc=num_fac_loc(facei,elem_loc);
2610 int faceToComplete;
2611 if (facei_loc!=-1)
2612 faceToComplete=faceiAss;
2613 else
2614 {
2615 faceToComplete=facei;
2616 facei_loc=num_fac_loc(faceiAss,elem_loc);
2617 assert(facei_loc!=-1);
2618 }
2619
2620 //Type de l'element : le nombre de faces de Dirichlet
2621 //qu'il contient
2622 int type_elem=elem_nb_faces_dirichlet_(elem);
2623 double coeff=formule(type_elem);
2624
2625 //Calcul de la divergence par element
2626 double div=0.;
2627 for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
2628 {
2629 facei=elem_faces(elem,facei_loc);
2630 int signe=(face_voisins(facei,0)==elem)? 1.:-1.;
2631
2632 for (int dim=0; dim<dimension; dim++)
2633 div+=signe*face_normales(facei,dim)*tab_vitesse(facei,dim);
2634 }
2635 div*=coeff;
2636 if (!marq) div/=porosite_elem(elem);
2637
2638 //Calcul de la partie compressible
2639 for (int dim=0; dim<nb_comp; dim++)
2640 {
2641 int ligne=faceToComplete*nb_comp+dim;
2642 matrice(ligne,ligne)+=div;
2643 }
2644 }
2645 }
2646 }
2647 }
2648}
2649
2650void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_antidiffusion(const DoubleTab& Kij, const DoubleTab& transporte,
2651 Matrice_Morse& matrice) const
2652{
2653 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
2654 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
2655 const IntTab& elem_faces=domaine_VEF.elem_faces();
2656 const IntTab& face_voisins = domaine_VEF.face_voisins();
2657 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
2658 const int nb_faces_elem=elem_faces.line_size();
2659 const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
2660 const int nb_comp=transporte.line_size();
2661
2662 int elem=0, elem_loc=0, facei=0,facei_loc=0, faceiAss=0, ind_face=0,ind_faceiAss=0;
2663 int facej=0,facej_loc=0, ligne=0,colonne=0, dim=0, face_amont=0,face_aval=0;
2664 int faceToComplete=0, num1=0,num2=0, n_bord=0;
2665 double kij=0.,kji=0.,dij=0., lij=0.,lji=0., daij=0.;
2666 double delta=0., coeffij=0.,coeffji=0., coeff=0., R=0.;
2667
2668 //Pour le limiteur
2669 ArrOfDouble P_plus(nb_comp),P_moins(nb_comp);
2670 ArrOfDouble Q_plus(nb_comp),Q_moins(nb_comp);
2671 P_plus=0., P_moins=0., Q_plus=0., Q_moins=0.;
2672
2673 const DoubleVect& transporteV = transporte;
2674 const ArrOfDouble& alpha_tab = alpha_tab_;
2675 const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
2676 ToDo_Kokkos("critical");
2677 for (elem=0; elem<nb_elem_tot; elem++)
2678 for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
2679 {
2680 facei=elem_faces(elem,facei_loc);
2681 P_plus=0., P_moins=0., Q_plus=0., Q_moins=0.;
2682 calculer_senseur(Kij,transporteV,nb_comp,facei,elem_faces,face_voisins,num_fac_loc,P_plus,P_moins,Q_plus,Q_moins);
2683 for (facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
2684 if (facej_loc!=facei_loc)
2685 {
2686 facej=elem_faces(elem,facej_loc);
2687
2688 kij = Kij(elem,facei_loc,facej_loc);
2689 kji = Kij(elem,facej_loc,facei_loc);
2690 dij = Dij(elem,facei_loc,facej_loc,Kij);
2691 lij = kij+dij;
2692 lji = kji+dij;
2693 assert(lij>=0);
2694 assert(lji>=0);
2695
2696 if (lij<=lji) //facei est amont
2697 {
2698 face_amont = facei;
2699 face_aval = facej;
2700
2701 //Si lij==lji, on passe deux foix dans la boucle
2702 //d'ou la presence du coefficient 1/2
2703 coeff = 1.*(lij<lji)+0.5*(lij==lji);
2704 assert(coeff==1. || coeff==0.5);
2705
2706 for (dim=0; dim<nb_comp; dim++)
2707 {
2708 ligne=face_amont*nb_comp+dim;
2709 colonne=face_aval*nb_comp+dim;
2710
2711 delta=transporteV[ligne]-transporteV[colonne];
2712
2713 //Limiteur de pente
2714 // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/P_plus(dim);
2715 // else R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/P_moins(dim);
2716
2717 // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/(P_plus(dim)+DMINFLOAT);
2718 // else R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/(P_moins(dim)+DMINFLOAT);
2719
2720 if (delta>=0.) R=(std::fabs(P_plus[dim])<DMINFLOAT) ? 0. : Q_plus[dim]/P_plus[dim];
2721 else R=(std::fabs(P_moins[dim])<DMINFLOAT) ? 0. : Q_moins[dim]/P_moins[dim];
2722
2723
2724 daij=minimum(limiteur(R)*dij,lji);
2725 assert(daij>=0);
2726 assert(daij<=lji);
2727 coeffij=alpha_tab_[face_amont]*beta_[face_amont]*daij;
2728 coeffji=alpha_tab_[face_aval]*beta_[face_aval]*daij;
2729
2730 //Calcul de la matrice
2731 matrice(ligne,ligne)-=coeffij*coeff;
2732 matrice(ligne,colonne)+=coeffij*coeff;
2733 matrice(colonne,colonne)-=coeffji*coeff;
2734 matrice(colonne,ligne)+=coeffji*coeff;
2735 }
2736 }
2737 }
2738 }
2739
2740 //
2741 //Pour la periodicite
2742 //
2743 for (n_bord=0; n_bord<nb_bord; n_bord++)
2744 {
2745 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
2746 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2747 num1 = 0;
2748 num2=le_bord.nb_faces();//pour ne parcourir que les faces reelles
2749
2750 if (sub_type(Periodique,la_cl.valeur()))
2751 {
2752 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
2753
2754 //Pour le limiteur
2755 ArrOfDouble Pj_plus(nb_comp),Pj_moins(nb_comp);
2756 ArrOfDouble Qj_plus(nb_comp),Qj_moins(nb_comp);
2757 Pj_plus=0., Pj_moins=0.;
2758 Qj_plus=0., Qj_moins=0.;
2759
2760 for (ind_face=num1; ind_face<num2; ind_face++)
2761 {
2762 facei=le_bord.num_face(ind_face);
2763 ind_faceiAss=la_cl_perio.face_associee(ind_face);
2764 faceiAss=le_bord.num_face(ind_faceiAss);
2765
2766 //Pour ne parcourir qu'une seule fois les faces perio
2767 if (facei<faceiAss)
2768 for (elem_loc=0; elem_loc<2; elem_loc++)
2769 {
2770 elem=face_voisins(facei,elem_loc);
2771 assert(elem!=-1);
2772
2773 //Calcul du numero local de la face dans "elem"
2774 facei_loc=num_fac_loc(facei,elem_loc);
2775 if (facei_loc!=-1)
2776 faceToComplete=faceiAss;
2777 else
2778 {
2779 faceToComplete=facei;
2780 facei_loc=num_fac_loc(faceiAss,elem_loc);
2781 assert(facei_loc!=-1);
2782 }
2783
2784 //Calcul du coefficient a rajouter dans la matrice
2785 P_plus=0., P_moins=0.;
2786 Q_plus=0., Q_moins=0.;
2787 calculer_senseur(Kij,transporteV,nb_comp,faceToComplete,elem_faces,face_voisins,num_fac_loc,P_plus,P_moins,Q_plus,Q_moins);
2788
2789 for (facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
2790 if (facej_loc!=facei_loc)
2791 {
2792 facej=elem_faces(elem,facej_loc);
2793
2794 kij = Kij(elem,facei_loc,facej_loc);
2795 kji = Kij(elem,facej_loc,facei_loc);
2796 dij = Dij(elem,facei_loc,facej_loc,Kij);
2797 lij = kij+dij;
2798 lji = kji+dij;
2799 assert(lij>=0);
2800 assert(lji>=0);
2801
2802 if (lij<=lji) //faceToComplete est amont
2803 {
2804 face_amont=faceToComplete;
2805 face_aval=facej;
2806
2807 //Si lij==lji, on passe deux foix dans la boucle
2808 //d'ou la presence du coefficient 1/2
2809 coeff = 1.*(lij<lji)+0.5*(lij==lji);
2810 assert(coeff==1. || coeff==0.5);
2811
2812 for (dim=0; dim<nb_comp; dim++)
2813 {
2814 ligne=face_amont*nb_comp+dim;
2815 colonne=face_aval*nb_comp+dim;
2816 delta=transporteV[ligne]-transporteV[colonne];
2817
2818 //Limiteur de pente
2819 // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/P_plus(dim);
2820 // else R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/P_moins(dim);
2821
2822 // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/(P_plus(dim)+DMINFLOAT);
2823 // else R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/(P_moins(dim)+DMINFLOAT);
2824
2825 if (delta>=0.) R=(std::fabs(P_plus[dim])<DMINFLOAT) ? 0. : Q_plus[dim]/P_plus[dim];
2826 else R=(std::fabs(P_moins[dim])<DMINFLOAT) ? 0. : Q_moins[dim]/P_moins[dim];
2827
2828 daij=minimum(limiteur(R)*dij,lji);
2829 assert(daij>=0);
2830 assert(daij<=lji);
2831 coeffij=alpha_tab[face_amont]*beta_[face_amont]*daij;
2832
2833 //Calcul de la matrice
2834 matrice(ligne,ligne)-=coeffij*coeff;
2835 matrice(ligne,colonne)+=coeffij*coeff;
2836 }
2837 }
2838 else //faceToComplete est aval
2839 {
2840 face_aval=faceToComplete;
2841 face_amont=facej;
2842 coeff=1.;
2843 Pj_plus=0., Pj_moins=0., Qj_plus=0., Qj_moins=0.;
2844 calculer_senseur(Kij,transporteV,nb_comp,facej,elem_faces,face_voisins,num_fac_loc,Pj_plus,Pj_moins,Qj_plus,Qj_moins);
2845
2846 for (dim=0; dim<nb_comp; dim++)
2847 {
2848 ligne=face_amont*nb_comp+dim;
2849 colonne=face_aval*nb_comp+dim;
2850
2851 delta=transporteV[ligne]-transporteV[colonne];
2852
2853 //Limiteur de pente
2854 // if (delta>=0.) R=(Pj_plus(dim)==0.) ? 0. : Qj_plus(dim)/Pj_plus(dim);
2855 // else R=(Pj_moins(dim)==0.) ? 0. : Qj_moins(dim)/Pj_moins(dim);
2856
2857 // if (delta>=0.) R=(Pj_plus(dim)==0.) ? 0. : Qj_plus(dim)/(Pj_plus(dim)+DMINFLOAT);
2858 // else R=(Pj_moins(dim)==0.) ? 0. : Qj_moins(dim)/(Pj_moins(dim)+DMINFLOAT);
2859
2860 if (delta>=0.) R=(std::fabs(Pj_plus[dim])<DMINFLOAT) ? 0. : Qj_plus[dim]/Pj_plus[dim];
2861 else R=(std::fabs(Pj_moins[dim])<DMINFLOAT) ? 0. : Qj_moins[dim]/Pj_moins[dim];
2862
2863 daij=minimum(limiteur(R)*dij,lij);
2864 assert(daij>=0);
2865 assert(daij<=lij);
2866 coeffij=alpha_tab[face_aval]*beta_[face_aval]*daij;
2867
2868 //Calcul de la matrice
2869 matrice(colonne,colonne)-=coeffij*coeff;
2870 matrice(colonne,ligne)+=coeffij*coeff;
2871 }
2872 }
2873 }
2874 }
2875 }
2876 }
2877 }
2878}
2879
2880void Op_Conv_EF_VEF_P1NC_Stab::test_implicite() const
2881{
2882 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
2883 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
2884
2885 const DoubleTab& unknown=equation().inconnue().valeurs();
2886 const DoubleTab& tab_vitesse=vitesse_->valeurs();
2887
2888 DoubleTab tab_test(unknown);
2889 DoubleVect& test2 = tab_test;
2890 test2 = 0.;
2891
2892 DoubleTab resuExp(unknown);
2893 DoubleVect& resu2Exp = resuExp;
2894 resu2Exp = 0.;
2895
2896 DoubleTab resuImp(unknown);
2897 DoubleVect& resu2Imp = resuImp;
2898 resu2Imp = 0.;
2899
2900 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
2901 const int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
2902 const int nb_faces_tot=domaine_VEF.nb_faces_tot();
2903 const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
2904
2905 DoubleTab Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
2906 Kij=0.;
2907
2908 const int nb_comp=unknown.line_size();
2909 int size=unknown.dimension(0);
2910 int face=0,face2=0, faceAss=0, ind_face=0,ind_faceAss=0, n_bord=0, num1=0,num2=0;
2911
2912 SFichier testResu("test.txt");
2913 SFichier testMat("matrice.txt");
2914
2915 //
2916 //Pour la periodicite
2917 //
2918 IntTab faces_associees(nb_faces_tot);
2919 for (face=0; face<nb_faces_tot; face++)
2920 faces_associees(face)=face;
2921
2922 for (n_bord=0; n_bord<nb_bord; n_bord++)
2923 {
2924 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
2925 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
2926 num1 = 0;
2927 num2=le_bord.nb_faces_tot();//pour ne pas en oublier
2928
2929 if (sub_type(Periodique,la_cl.valeur()))
2930 {
2931 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
2932
2933 for (ind_face=num1; ind_face<num2; ind_face++)
2934 {
2935 face=le_bord.num_face(ind_face);
2936 ind_faceAss=la_cl_perio.face_associee(ind_face);
2937 faceAss=le_bord.num_face(ind_faceAss);
2938
2939 //Pour ne parcourir que la moitie des faces periodiques
2940 if (face<faceAss)
2941 {
2942 faces_associees(face)=faceAss;
2943 faces_associees(faceAss)=face;
2944 }
2945 }
2946 }
2947 }
2948 //
2949 //Fin du traitement pour la periodicite
2950 //
2951
2952 calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
2953
2954 //
2955 //Construction de la matrice
2956 //
2957 Matrice_Morse matrice;
2958 dimensionner(matrice);
2959 if (is_compressible_)
2960 ajouter_contribution_partie_compressible(unknown,tab_vitesse,matrice);
2961 ajouter_contribution_operateur_centre(Kij,unknown,matrice);
2962 ajouter_contribution_diffusion(Kij,unknown,matrice);
2963 ajouter_contribution_antidiffusion(Kij,unknown,matrice);
2964 matrice.imprimer_formatte(testMat);
2965 //
2966 //Fin de la construction de la matrice
2967 //
2968
2969 //
2970 //Calcul de l'operateur explicite et comparaison
2971 //
2972 for (face=0; face<size; face++)
2973 {
2974 test2[face]=1.;
2975 test2[faces_associees(face)]=1.;
2976
2977 /* Calcul de l'operateur explicite */
2978 resuExp=0.;
2979 if (is_compressible_)
2980 ajouter_partie_compressible(tab_test,resuExp,tab_vitesse);
2981 ajouter_operateur_centre(Kij,tab_test,resuExp);
2982 ajouter_diffusion(Kij,tab_test,resuExp);
2983 ajouter_antidiffusion(Kij,tab_test,resuExp);
2984 mettre_a_jour_pour_periodicite(resuExp);
2985
2986 /* Calcul de l'operateur implicite */
2987 resuImp=0.;
2988 matrice.ajouter_multvect_(tab_test,resuImp);
2989
2990 /* Calcul de la difference */
2991 resuExp+=resuImp;
2992
2993 /* Affichage de la difference */
2994 testResu<<"*************************"<<finl;
2995 testResu<<"Face test : "<<face<<finl;
2996 for (face2=0; face2<size; face2++)
2997 if (resu2Exp[face2]<=1.e-13)
2998 testResu<<face2<<" OK"<<finl;
2999 else
3000 testResu<<face2<<" residu : "<<resu2Exp[face2]<<finl;
3001 testResu<<"*************************"<<finl;
3002
3003 test2[face]=0.;
3004 test2[faces_associees(face)]=0.;
3005 }
3006 //
3007 //Fin du calcul de l'operateur explicite et comparaison
3008 //
3009
3010 Process::exit();
3011}
3012
Champ_Inc_base::valeurs
DoubleTab & valeurs() override
Renvoie le tableau des valeurs du champ au temps courant.
Definition Champ_Inc_base.h:77
Champ_P1NC
Definition Champ_P1NC.h:35
Cond_lim
classe Cond_lim Classe generique servant a representer n'importe quelle classe
Definition Cond_lim.h:31
Dirichlet_homogene
Classe Dirichlet_homogene Cette classe est la classe de base de la hierarchie des conditions aux limi...
Definition Dirichlet_homogene.h:31
Dirichlet
classe Dirichlet Cette classe est la classe de base de la hierarchie des conditions aux limites de ty...
Definition Dirichlet.h:31
Domaine_32_64::ss_domaine
const Sous_Domaine_t & ss_domaine(int i) const
Definition Domaine.h:290
Domaine_32_64::nb_faces_elem
int nb_faces_elem(int=0) const
Renvoie le nombre de face de type i des elements geometriques constituants le domaine.
Definition Domaine.h:484
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_dis_base::nb_cond_lim
int nb_cond_lim() const
Renvoie le nombre de conditions aux limites.
Definition Domaine_Cl_dis_base.cpp:425
Domaine_Cl_dis_base::les_conditions_limites
const Cond_lim & les_conditions_limites(int) const
Renvoie la i-ieme condition aux limites.
Definition Domaine_Cl_dis_base.cpp:386
Domaine_VEF
class Domaine_VEF
Definition Domaine_VEF.h:54
Domaine_VF::nb_faces
int nb_faces() const
renvoie le nombre global de faces.
Definition Domaine_VF.h:471
Domaine_VF::nb_faces_tot
int nb_faces_tot() const
renvoie le nombre total de faces.
Definition Domaine_VF.h:481
Domaine_VF::face_normales
virtual double face_normales(int face, int comp) const
Definition Domaine_VF.h:47
Domaine_VF::get_num_fac_loc
const IntTab & get_num_fac_loc() const
Definition Domaine_VF.h:140
Domaine_VF::elem_faces
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
Domaine_VF::face_voisins
int face_voisins(int num_face, int i) const
renvoie l'element voisin de numface dans la direction i.
Definition Domaine_VF.h:418
Domaine_VF::nb_faces_bord
int nb_faces_bord() const
renvoie le nombre de faces sur lesquelles sont appliquees les conditions limites :
Definition Domaine_VF.h:513
Domaine_dis_base
classe Domaine_dis_base Cette classe est la base de la hierarchie des domaines discretisees.
Definition Domaine_dis_base.h:39
Domaine_dis_base::nombre_de_sous_domaines_dis
int nombre_de_sous_domaines_dis() const
Definition Domaine_dis_base.cpp:36
Domaine_dis_base::nb_elem_tot
int nb_elem_tot() const
Definition Domaine_dis_base.h:50
Domaine_dis_base::sous_domaine_dis
const Sous_domaine_dis_base & sous_domaine_dis(int i) const
Definition Domaine_dis_base.cpp:41
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Echange_impose_base
classe Echange_impose_base: Cette condition limite sert uniquement pour l'equation d'energie.
Definition Echange_impose_base.h:41
Entree
Class defining operators and methods for all reading operation in an input flow (file,...
Definition Entree.h:42
Equation_base::milieu
virtual const Milieu_base & milieu() const =0
Equation_base::inconnue
virtual const Champ_Inc_base & inconnue() const =0
Equation_base::probleme
Probleme_base & probleme()
Renvoie le probleme associe a l'equation.
Definition Equation_base.cpp:747
Equation_base::schema_temps
Schema_Temps_base & schema_temps()
Renvoie le schema en temps associe a l'equation.
Definition Equation_base.cpp:865
Front_VF
class Front_VF
Definition Front_VF.h:36
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
Front_VF::nb_faces_tot
int nb_faces_tot() const
Definition Front_VF.h:58
Front_VF::num_face
int num_face(const int) const
Definition Front_VF.h:68
Matrice_Morse
Classe Matrice_Morse Represente une matrice M (creuse), non necessairement carree.
Definition Matrice_Morse.h:50
Matrice_Morse::imprimer_formatte
Sortie & imprimer_formatte(Sortie &s) const override
Definition Matrice_Morse.cpp:78
Matrice_Morse::ajouter_multvect_
DoubleVect & ajouter_multvect_(const DoubleVect &, DoubleVect &) const override
Operation de multiplication-accumulation (saxpy) matrice vecteur.
Definition Matrice_Morse.cpp:877
Milieu_base::porosite_elem
DoubleVect & porosite_elem()
Definition Milieu_base.h:58
Milieu_base::porosite_face
DoubleVect & porosite_face()
Definition Milieu_base.h:62
MorEqn::equation
const Equation_base & equation() const
Renvoie la reference sur l'equation pointe par MorEqn::mon_equation.
Definition MorEqn.h:62
Neumann_homogene
Classe Neumann_homogene Cette classe est la classe de base de la hierarchie des conditions aux limite...
Definition Neumann_homogene.h:31
Neumann_sortie_libre::val_ext
double val_ext(int i) const override
Renvoie la valeur de la i-eme composante du champ impose a l'exterieur de la frontiere.
Definition Neumann_sortie_libre.cpp:120
Neumann_val_ext
Classe Neumann_val_ext Cette classe est la classe de base de la hierarchie des conditions.
Definition Neumann_val_ext.h:32
Neumann
Classe Neumann Cette classe est la classe de base de la hierarchie des conditions aux limites de type...
Definition Neumann.h:31
Objet_U::dimension
static int dimension
Definition Objet_U.h:99
Objet_U::que_suis_je
const Nom & que_suis_je() const
renvoie la chaine identifiant la classe.
Definition Objet_U.cpp:104
Objet_U::readOn
virtual Entree & readOn(Entree &)
Lecture d'un Objet_U sur un flot d'entree Methode a surcharger.
Definition Objet_U.cpp:293
Objet_U::est_egal_a
virtual int est_egal_a(const Objet_U &) const
Renvoie 1 si l'objet x et *this sont une seule et meme instance (meme adresse en memoire).
Definition Objet_U.cpp:301
Objet_U::printOn
virtual Sortie & printOn(Sortie &) const
Ecriture de l'objet sur un flot de sortie Methode a surcharger.
Definition Objet_U.cpp:282
Op_Conv_EF_VEF_P1NC_Stab
class Op_Conv_EF_VEF_P1NC_Stab
Definition Op_Conv_EF_VEF_P1NC_Stab.h:44
Op_Conv_EF_VEF_P1NC_Stab::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:556
Op_Conv_EF_VEF_P1NC_Stab::ajouter_diffusion
DoubleTab & ajouter_diffusion(const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:811
Op_Conv_EF_VEF_P1NC_Stab::remplir_fluent
void remplir_fluent() const override
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:534
Op_Conv_EF_VEF_P1NC_Stab::completer
void completer() override
Associe l'operateur au domaine_dis, le domaine_Cl_dis, et a l'inconnue de son equation.
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:2146
Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_diffusion
void ajouter_contribution_diffusion(const DoubleTab &, const DoubleTab &, Matrice_Morse &) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:2391
Op_Conv_EF_VEF_P1NC_Stab::calculer_coefficients_operateur_centre
void calculer_coefficients_operateur_centre(DoubleTab &, const int, const DoubleTab &vitesse) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:315
Op_Conv_EF_VEF_P1NC_Stab::calculer_flux_bords
public_for_cuda void calculer_flux_bords(const DoubleTab &, const DoubleTab &, const DoubleTab &) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:716
Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_operateur_centre
void ajouter_contribution_operateur_centre(const DoubleTab &, const DoubleTab &, Matrice_Morse &) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:2278
Op_Conv_EF_VEF_P1NC_Stab::ajouter_operateur_centre
DoubleTab & ajouter_operateur_centre(const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:774
Op_Conv_EF_VEF_P1NC_Stab::ajouter_partie_compressible
DoubleTab & ajouter_partie_compressible(const DoubleTab &, DoubleTab &, const DoubleTab &vitesse) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:629
Op_Conv_EF_VEF_P1NC_Stab::ajouter_antidiffusion
DoubleTab & ajouter_antidiffusion(const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:853
Op_Conv_EF_VEF_P1NC_Stab::modifier_pour_Cl
void modifier_pour_Cl(Matrice_Morse &, DoubleTab &) const override
DOES NOTHING - to override in derived classes.
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:2273
Op_Conv_EF_VEF_P1NC_Stab::mettre_a_jour_pour_periodicite
void mettre_a_jour_pour_periodicite(DoubleTab &) const
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:1611
Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution
void ajouter_contribution(const DoubleTab &, Matrice_Morse &) const override
Definition Op_Conv_EF_VEF_P1NC_Stab.cpp:2233
Op_Conv_VEF_Face
class Op_Conv_VEF_Face
Definition Op_Conv_VEF_Face.h:35
Op_Conv_VEF_Face::dimensionner
void dimensionner(Matrice_Morse &) const override
on dimensionne notre matrice au moyen de la methode dimensionner de la classe Op_VEF_Face.
Definition Op_Conv_VEF_Face.h:81
Op_Conv_VEF_Face::ajouter_contribution
virtual void ajouter_contribution(const DoubleTab &, Matrice_Morse &) const
Definition Op_Conv_VEF_Face.cpp:1356
Op_Conv_VEF_Face::modifier_pour_Cl
void modifier_pour_Cl(Matrice_Morse &, DoubleTab &) const override
On modifie le second membre et la matrice dans le cas des conditions de dirichlet.
Definition Op_Conv_VEF_Face.h:89
Op_Conv_VEF_Face::completer
void completer() override
Associe l'operateur au domaine_dis, le domaine_Cl_dis, et a l'inconnue de son equation.
Definition Op_Conv_VEF_Face.cpp:146
Op_Conv_VEF_base::fluent_
DoubleVect fluent_
Definition Op_Conv_VEF_base.h:61
Op_Conv_VEF_base::phi_u_transportant
int phi_u_transportant(const Equation_base &eq) const
definit si l'on convecte psi avec phi*u ou avec u
Definition Op_Conv_VEF_base.cpp:44
Op_Conv_VEF_base::vitesse
const Champ_Inc_base & vitesse() const
Definition Op_Conv_VEF_base.h:41
Op_VEF_Face::modifier_flux
void modifier_flux(const Operateur_base &) const
Definition Op_VEF_Face.cpp:338
Operateur_base::flux_bords_
DoubleTab flux_bords_
Definition Operateur_base.h:140
Operateur_base::flux_bords
DoubleTab & flux_bords()
Definition Operateur_base.h:107
Periodique
classe Periodique Cette classe represente une condition aux limites periodique.
Definition Periodique.h:31
Periodique::face_associee
int face_associee(int i) const
Definition Periodique.h:35
Probleme_base::domaine
const Domaine & domaine() const
Renvoie le domaine associe au probleme.
Definition Probleme_base.cpp:613
Process::Kokkos_exit
static KOKKOS_INLINE_FUNCTION void Kokkos_exit(const char *)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.h:173
Process::Journal
static Sortie & Journal(int message_level=0)
Renvoie un objet statique de type Sortie qui sert de journal d'evenements.
Definition Process.cpp:588
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Schema_Temps_base::temps_courant
double temps_courant() const
Renvoie le temps courant.
Definition Schema_Temps_base.h:445
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
Sous_domaine_VF
Cette classe abstraite contient les informations geometrique de sous-domaine communes aux methodes de...
Definition Sous_domaine_VF.h:51
Sous_domaine_VF::les_faces
const IntTab & les_faces() const
Definition Sous_domaine_VF.h:58
Sous_domaine_dis_base::sous_domaine
const Sous_Domaine & sous_domaine() const
Definition Sous_domaine_dis_base.cpp:34
Symetrie
classe Symetrie Sur les faces de symetrie on a les proprietes suivantes:
Definition Symetrie.h:37
TRUSTArray::size_array
_SIZE_ size_array() const
Definition TRUSTArray.tpp:187
TRUSTTab::nb_dim
int nb_dim() const
Definition TRUSTTab.h:199
TRUSTTab::view_ro
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, ConstView< _TYPE_, _SHAPE_ > > view_ro() const
Definition TRUSTTab.h:261
TRUSTTab::view_rw
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, View< _TYPE_, _SHAPE_ > > view_rw()
Definition TRUSTTab.h:291
TRUSTTab::dimension
_SIZE_ dimension(int d) const
Definition TRUSTTab.tpp:133
TRUSTVect::size
_SIZE_ size() const
Definition TRUSTVect.tpp:45
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67