TRUST 1.9.8
HPC thermohydraulic platform
Loading...
Searching...
No Matches
Op_Conv_Muscl_New_VEF_Face.cpp
1/****************************************************************************
2* Copyright (c) 2025, 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_Muscl_New_VEF_Face.h>
17#include <Champ_P1NC.h>
18#include <Schema_Temps_base.h>
19#include <Porosites_champ.h>
20#include <Sous_domaine_VF.h>
21#include <Probleme_base.h>
22#include <ArrOfBit.h>
23#include <TRUSTTabs.h>
24#include <Neumann.h>
25#include <Symetrie.h>
26#include <Dirichlet_homogene.h>
27#include <Neumann_homogene.h>
28#include <Echange_impose_base.h>
29#include <Periodique.h>
30#include <Neumann_sortie_libre.h>
31#include <Tetra_VEF.h>
32#include <Tri_VEF.h>
33
34static KOKKOS_INLINE_FUNCTION double minmod(double r)
35{
36 if (r<=0.) return 0.;
37 if (r<=1.) return r;
38 else return 1.;
39}
40static KOKKOS_INLINE_FUNCTION double optimum(double a, double b)
41{
42 if (a>=0.)
43 if (b>=0.)
44 return (a<=b)?a:b;
45 else
46 return 0.;
47 else if (b>=0.)
48 return 0.;
49 else
50 return (a>=b)?a:b;
51}
52
53Implemente_instanciable(Op_Conv_Muscl_New_VEF_Face,"Op_Conv_Muscl_New_VEF_P1NC",Op_Conv_VEF_Face);
54// XD convection_muscl_new convection_deriv muscl_new NO_BRACE Only for VEF discretization.
55
56Sortie& Op_Conv_Muscl_New_VEF_Face::printOn(Sortie& s ) const
57{
58 return s << que_suis_je();
59}
60
61Entree& Op_Conv_Muscl_New_VEF_Face::readOn(Entree& s )
62{
63 //Les mots a reconnaitre
64 /*
65 Motcles les_mots(8);
66 {
67 les_mots[0] = "centre";
68 les_mots[1] = "amont";
69 les_mots[2] = "stab";
70 les_mots[3] = "alpha";
71 les_mots[4] = "centre_old";
72 les_mots[5] = "version";
73 les_mots[6] = "limiteur";
74 les_mots[7] = "facsec_auto";
75 }
76 Motcles les_mots2(3);
77 {
78 les_mots2[0] = "vanleer";
79 les_mots2[1] = "minmod";
80 les_mots2[2] = "superbee";
81 } */
82
83 //Les variables a instancier
84 alpha_=1.;
85 //limiteur_=&minmod;
86
87 return s;
88}
89
90////////////////////////////////////////////////////////////////////
91//
92// Implementation des fonctions
93//
94// de la classe Op_Conv_Muscl_New_VEF_Face
95//
96////////////////////////////////////////////////////////////////////
97
98void Op_Conv_Muscl_New_VEF_Face::calculer_coefficients_operateur_centre(DoubleTab& tab_Kij, DoubleTab& tab_Cij, DoubleTab& tab_Sij, DoubleTab& tab_Sij2, const int nb_comp, const DoubleTab& tab_velocity) const
99{
100 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
101 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
102 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
103 const int nb_faces_elem = domaine_VEF.elem_faces().dimension(1);
104 const int nfa7 = domaine_VEF.type_elem().nb_facette();
105 const int marq = phi_u_transportant(equation());
106 int nb_dim = Objet_U::dimension;
107
108 assert(tab_Kij.nb_dim()==3);
109 assert(tab_Kij.dimension(0)==nb_elem_tot);
110 assert(tab_Kij.dimension(1)==tab_Kij.dimension(2));
111 assert(tab_Kij.dimension(1)==nb_faces_elem);
112
113 assert(tab_Cij.nb_dim()==2);
114 assert(tab_Cij.dimension(0)==nb_elem_tot);
115 assert(tab_Cij.dimension(1)==nfa7);
116
117 assert(tab_Sij.nb_dim()==2);
118 assert(tab_Sij.dimension(0)==nb_elem_tot);
119 assert(tab_Sij.dimension(1)==nfa7);
120
121 if (nb_dim==3)
122 {
123 assert(tab_Sij2.nb_dim()==2);
124 assert(tab_Sij2.dimension(0)==nb_elem_tot);
125 assert(tab_Sij2.dimension(1)==nfa7);
126 }
127
128 //
129 //Calcul des coefficients de l'operateur
130 //
131 CDoubleTabView3 facette_normales = domaine_VEF.facette_normales().view_ro<3>();
132 CDoubleTabView3 facette_normales_Cl = domaine_Cl_VEF.normales_facettes_Cl().view_ro<3>();
133 CIntArrView rang_elem_non_std = domaine_VEF.rang_elem_non_std().view_ro();
134 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
135 CDoubleTabView velocity = tab_velocity.view_ro();
136 CDoubleTabView vitesse = vitesse_->valeurs().view_ro();
137 CDoubleArrView porosite_face = equation().milieu().porosite_face().view_ro();
138 CDoubleArrView porosite_elem = equation().milieu().porosite_elem().view_ro();
139 CIntArrView type_elem_Cl = domaine_Cl_VEF.type_elem_Cl().view_ro();
140 CIntTabView KEL = domaine_VEF.type_elem().KEL().view_ro();
141 DoubleTabView Cij = tab_Cij.view_wo();
142 DoubleTabView Sij = tab_Sij.view_wo();
143 DoubleTabView Sij2 = tab_Sij2.view_wo();
144 DoubleTabView3 Kij = tab_Kij.view_rw<3>();
145 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem_tot, KOKKOS_LAMBDA(
146 const int elem)
147 {
148 double vc[3];
149 double vs[3];
150 double vsom[12];
151 int face[4] = {-1, -1, -1, -1};
152 int rang=rang_elem_non_std(elem);
153 int itypcl=(rang==-1)?0:type_elem_Cl(rang);
154
155 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
156 face[facei_loc]=elem_faces(elem,facei_loc);
157
158 for (int dim=0; dim<nb_dim; dim++)
159 {
160 vc[dim] = 0;
161 vs[dim] = 0;
162 for (int facei_loc = 0; facei_loc < nb_faces_elem; facei_loc++)
163 vs[dim] += velocity(face[facei_loc], dim);
164 }
165 for (int dim=0; dim<nb_dim; dim++)
166 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
167 vsom[facei_loc*nb_dim+dim]=vs[dim]-nb_dim*velocity(face[facei_loc],dim);
168
169 //domaine_VEF.type_elem()->calcul_vc(face,vc,vs,vsom,vitesse,itypcl,porosite_face);
170 if (nb_dim==3)
171 calcul_vc_tetra_views(face,vc,vs,vsom,vitesse,itypcl,porosite_face);
172 else
173 calcul_vc_tri_views(face,vc,vs,vsom,vitesse,itypcl,porosite_face);
174
175 if (marq==0)
176 {
177 double porosite=1./porosite_elem(elem);
178 for (int dim=0; dim<nb_dim; dim++)
179 {
180 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
181 vsom[facei_loc*nb_dim+dim] *= porosite;
182 vs[dim] *= porosite;
183 vc[dim] *= porosite;
184 }
185 }
186 for (int fa7=0; fa7<nfa7; fa7++)
187 {
188 int facei_loc=KEL(0,fa7);
189 int facej_loc=KEL(1,fa7);
190
191 Cij(elem,fa7)=0.;
192 Sij(elem,fa7)=0.;
193 Sij2(elem,fa7)=0.;
194 for (int dim=0; dim<nb_dim; dim++)
195 {
196 double facette_normale = (rang==-1 ? facette_normales(elem, fa7, dim) : facette_normales_Cl(rang,fa7,dim));
197 Cij(elem, fa7) += vc[dim] * facette_normale;
198 Sij(elem, fa7) += vsom[KEL(2, fa7) * nb_dim + dim] * facette_normale;
199 if (nb_dim == 3) Sij2(elem, fa7) += vsom[KEL(3, fa7) * nb_dim + dim] * facette_normale;
200 }
201 double psc=Cij(elem,fa7)+Sij(elem,fa7)+Sij2(elem,fa7);
202 psc/=nb_dim;
203
204 Kij(elem,facei_loc,facej_loc)=psc;
205 Kij(elem,facei_loc,facei_loc)-=psc;//pour l'aspect LED
206 Kij(elem,facej_loc,facei_loc)=-psc;
207 Kij(elem,facej_loc,facej_loc)+=psc;//pour l'aspect LED
208 }
209 });
210 end_gpu_timer(__KERNEL_NAME__);
211
212 //
213 // Correction des Kij pour Dirichlet !
214 //
215 {
216 int nb_bord=domaine_Cl_VEF.nb_cond_lim();
217
218 for (int n_bord=0; n_bord<nb_bord; n_bord++)
219 {
220 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
221
222 if ( (sub_type(Dirichlet,la_cl.valeur()))
223 || (sub_type(Dirichlet_homogene,la_cl.valeur()))
224 )
225 {
226 //On ne fait rien
227 }// sub_type Dirichlet
228
229
230 else if (sub_type(Neumann,la_cl.valeur()) || sub_type(Neumann_homogene,la_cl.valeur()))
231 {
232 //On ne fait rien
233 }//fin du if sur Neumann
234
235 else if (sub_type(Symetrie,la_cl.valeur()))
236 {
237 //On ne fait rien
238 }//fin du if sur Symetrie
239
240 else if (sub_type(Periodique,la_cl.valeur()))
241 {
242 //On ne fait rien
243 }//fin du if sur Periodique
244
245 else if (sub_type(Echange_impose_base,la_cl.valeur()))
246 {
247 //On ne fait rien
248 }//fin du if sur Echange_impose_base
249
250 else
251 {
252 Cerr << "Erreur Op_Conv_Muscl_New_VEF_Face::calculer_coefficients_operateur_centre()" << finl;
253 Cerr << "Condition aux limites " << la_cl.que_suis_je() << " non codee." << finl;
254 Cerr << "Sortie du programme." << finl;
255 Process::exit();
256 }//fin du else sur les autres conditions aux limites
257
258 }//fin des conditions aux limites
259 }
260
261 //
262 // Fin de la correction des Kij
263 //
264}
265
266
267void Op_Conv_Muscl_New_VEF_Face::
268calculer_flux_operateur_centre(DoubleTab& tab_Fij,const DoubleTab& tab_Kij,const DoubleTab& tab_Cij, const DoubleTab& tab_Sij, const DoubleTab& tab_Sij2, const int nb_comp, const DoubleTab& tab_velocity, const DoubleTab& tab_transporte) const
269{
270 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
271 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
272 const Domaine& domaine = domaine_VEF.domaine();
273
274 const DoubleTab& tab_vecteur_face_facette = ref_cast_non_const(Domaine_VEF,domaine_VEF).vecteur_face_facette();
275 const DoubleTab& tab_vecteur_face_facette_Cl = domaine_Cl_VEF.vecteur_face_facette_Cl();
276 const DoubleTab& tab_coord_sommets = domaine.coord_sommets();
277 const DoubleTab& tab_xv = domaine_VEF.xv();
278
279 const IntTab& tab_elem_faces = domaine_VEF.elem_faces();
280 // const IntTab& face_voisins = domaine_VEF.face_voisins();
281 const IntVect& tab_rang_elem_non_std = domaine_VEF.rang_elem_non_std();
282 const IntTab& tab_KEL=domaine_VEF.type_elem().KEL();
283
284 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
285 const int nb_faces_elem=tab_elem_faces.dimension(1);
286 const int nfa7 = domaine_VEF.type_elem().nb_facette();
287 const int nb_dim = Objet_U::dimension;
288
289 //DoubleTab gradient_elem(nb_elem_tot,nb_comp,nb_dim); //!< (du/dx du/dy dv/dx dv/dy) pour un poly gradient_elem=0.;
290 if (gradient_elem_.size_array() == 0) gradient_elem_.resize(nb_elem_tot, nb_comp, nb_dim); // (du/dx du/dy dv/dx dv/dy) pour un poly
291
292 assert(tab_Fij.nb_dim()==4);
293 assert(tab_Fij.dimension(0)==nb_elem_tot);
294 assert(tab_Fij.dimension(1)==tab_Fij.dimension(2));
295 assert(tab_Fij.dimension(1)==nb_faces_elem);
296 assert(tab_Fij.dimension(3)==nb_comp);
297
298 //
299 //Calcul des flux de l'operateur
300 //
301 Champ_P1NC::calcul_gradient(tab_transporte,gradient_elem_,domaine_Cl_VEF);
302
303 //Declaration des vues sur les tableaux TRUSTS
304 CIntArrView rang_elem_non_std = tab_rang_elem_non_std.view_ro();
305 CIntTabView elem_faces = tab_elem_faces.view_ro();
306 CIntTabView KEL = tab_KEL.view_ro();
307 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
308 CDoubleTabView Cij = tab_Cij.view_ro();
309 CDoubleTabView Sij = tab_Sij.view_ro();
310 CDoubleTabView Sij2 = tab_Sij2.view_ro();
311 CDoubleTabView transporteVect = tab_transporte.view_ro();
312 CDoubleTabView4 vecteur_face_facette = tab_vecteur_face_facette.view_ro<4>();
313 CDoubleTabView4 vecteur_face_facette_Cl = tab_vecteur_face_facette_Cl.view_ro<4>();
314 CDoubleTabView coord_sommets = tab_coord_sommets.view_ro();
315 CDoubleTabView xv = tab_xv.view_ro();
316 DoubleTabView4 Fij = tab_Fij.view_wo<4>();
317 CDoubleTabView3 gradient_elem = gradient_elem_.view_ro<3>(); //Gradient elem declare plus haut, je le renome pas _tab dans tout le fichier, donc convection un peu froissee ici
318 CIntTabView sommet_elem = domaine.les_elems().view_ro();//On n'utilise plus la fonction sommet_elem(.,.) de domaine.h
319
320 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
321 Kokkos::RangePolicy<>(0, nb_elem_tot), KOKKOS_LAMBDA(const int elem)
322 {
323 int rang = rang_elem_non_std(elem);
324 int face[4]; // Declare face inside the lambda, with compile time size ! (sneaky bug otherwise)
325
326 for (int facei_loc = 0; facei_loc < nb_faces_elem; facei_loc++)
327 face[facei_loc] = elem_faces(elem, facei_loc);
328
329 for (int fa7 = 0; fa7 < nfa7; fa7++)
330 {
331 int facei_loc = KEL(0, fa7);
332 int facej_loc = KEL(1, fa7);
333
334 int facei = face[facei_loc];
335 int facej = face[facej_loc];
336
337 //REMARQUE : ON N'OUBLIE PAS QUE LES NORMALES FA7 SONT ORIENTeES DE FACE_LOCI VERS FACE_LOCJ PAR DEFAUT
338 //LES PRODUITS SCALAIRES PSC_M, PSC_C, PSC_S SONT DONC CALCULeS POUR OBTENIR LE FLUX REcU PAR LA FACE FACE_LOCI
339 //ET EMIS PAR LA FACE_LOCJ
340 //SI L'ON VEUT OBTENIR LE FLUX REcU PAR LA FACE FACE_LOCJ ET EMIS PAR LA FACE_LOCI, IL FAUT MULTIPLIER LES
341 //PRODUITS SCALAIRES PSC_M, PSC_C, PSC_S PAR -1 POUR IMPLICITEMENT ReORIENTER LES NORMALES FA7 DE FACE_LOCJ VERS
342 //FACE_LOCI
343 //MEME REMARQUE POUR LE SCALAIRE "FLUX"
344
345 // Compute psc_m, psc_c, psc_s, psc_s2
346 const double psc_m = Kij(elem, facei_loc, facej_loc);
347 const double psc_c = Cij(elem, fa7);
348
349 // Arrays to handle both dimensions
350 double psc_s_array[2];
351 int s_array[2];
352
353 if (nb_dim == 2)
354 {
355 psc_s_array[0] = Sij(elem, fa7);
356 psc_s_array[1] = 0;
357 s_array[0] = sommet_elem(elem, KEL(2, fa7));
358 }
359 else // nb_dim == 3
360 {
361 psc_s_array[0] = Sij(elem, fa7);
362 psc_s_array[1] = Sij2(elem, fa7);
363 s_array[0] = sommet_elem(elem, KEL(2, fa7));
364 s_array[1] = sommet_elem(elem, KEL(3, fa7));
365 }
366
367 int face_amont, dir;
368 if (psc_m >= 0.)
369 face_amont = facei, dir = 0;
370 else
371 face_amont = facej, dir = 1;
372
373 for (int comp = 0; comp < nb_comp; comp++)
374 {
375 //Calcul de la valeur de l'inconnue aux points d'integration de la formule de Simspon
376 double inco_m = transporteVect(face_amont, comp);
377 for (int dim = 0; dim < nb_dim; dim++)
378 {
379 if (rang == -1)
380 inco_m += gradient_elem(elem, comp, dim) * vecteur_face_facette(elem, fa7, dim, dir);
381 else
382 inco_m += gradient_elem(elem, comp, dim) * vecteur_face_facette_Cl(rang, fa7, dim, dir);
383 }
384
385 // Compute inco_s_array
386 double inco_s_array[2];
387 int num_integration_points = (nb_dim == 2) ? 1 : 2;
388 for (int ip = 0; ip < num_integration_points; ip++)
389 {
390 double inco_s = transporteVect(face_amont, comp);
391 for (int dim = 0; dim < nb_dim; dim++)
392 {
393 inco_s += gradient_elem(elem, comp, dim) *
394 (coord_sommets(s_array[ip], dim) - xv(face_amont, dim));
395 }
396 inco_s_array[ip] = inco_s;
397 }
398
399 // Compute inco_c
400 double inco_c = (nb_dim == 2 ? 2.0 : 3.0) * inco_m;
401 for (int ip = 0; ip < num_integration_points; ip++)
402 {
403 inco_c -= inco_s_array[ip];
404 }
405
406 //Calcul du flux final : formule d'integration de Simpson
407 double flux;
408 if (nb_dim == 2)
409 {
410 flux = inco_c * psc_c;
411 flux += 4 * inco_m * psc_m;
412 flux += inco_s_array[0] * psc_s_array[0];
413 flux /= 6.0;
414 }
415 else // nb_dim == 3
416 {
417 double sum1 = (inco_s_array[0] + inco_s_array[1]) * (psc_s_array[0] + psc_s_array[1]);
418 double sum2 = (inco_s_array[0] + inco_c) * (psc_s_array[0] + psc_c);
419 double sum3 = (inco_s_array[1] + inco_c) * (psc_s_array[1] + psc_c);
420 flux = (sum1 + sum2 + sum3) / 12.0;
421 }
422 //Fij(elem,facei_loc,facej_loc,comp) est le flux recu par la facei_loc dans elem
423 //envoye par la face facej_loc dans elem pour la composante comp de l'inconnue
424 Fij(elem, facei_loc, facej_loc, comp) = flux - psc_m * transporteVect(facei, comp);
425 Fij(elem, facej_loc, facei_loc, comp) = psc_m * transporteVect(facej, comp) - flux;
426 }
427 }
428 });
429 end_gpu_timer(__KERNEL_NAME__);
430}
431
432void Op_Conv_Muscl_New_VEF_Face::
433modifier_flux_operateur_centre(DoubleTab& Fij,const DoubleTab& Kij,const DoubleTab& Cij, const DoubleTab& Sij, const DoubleTab& Sij2, const int nb_comp, const DoubleTab& velocity, const DoubleTab& transporte) const
434{
435 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
436 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
437
438 const DoubleTab& vecteur_face_facette = ref_cast_non_const(Domaine_VEF,domaine_VEF).vecteur_face_facette();
439 const DoubleTab& vecteur_face_facette_Cl = domaine_Cl_VEF.vecteur_face_facette_Cl();
440
441 const DoubleVect& transporteVect = transporte;
442
443 const IntTab& elem_faces = domaine_VEF.elem_faces();
444 const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
445 const IntTab& KEL=domaine_VEF.type_elem().KEL();
446
447 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
448 const int nb_faces_elem=elem_faces.dimension(1);
449 const int nfa7 = domaine_VEF.type_elem().nb_facette();
450 int nb_dim = Objet_U::dimension;
451
452 //DoubleTab gradient_elem(nb_elem_tot,nb_comp,nb_dim); //!< (du/dx du/dy dv/dx dv/dy) pour un poly gradient_elem=0.;
453 if (gradient_elem_.size_array() == 0) gradient_elem_.resize(nb_elem_tot, nb_comp, dimension); // (du/dx du/dy dv/dx dv/dy) pour un poly
454 IntTab face(nb_dim+1);
455
456 assert(Fij.nb_dim()==4);
457 assert(Fij.dimension(0)==nb_elem_tot);
458 assert(Fij.dimension(1)==Fij.dimension(2));
459 assert(Fij.dimension(1)==nb_faces_elem);
460 assert(Fij.dimension(3)==nb_comp);
461
462 //
463 //Calcul des flux de l'operateur
464 //
465 Champ_P1NC::calcul_gradient(transporte,gradient_elem_,domaine_Cl_VEF);
466 ToDo_Kokkos("critical");
467 for(int elem=0; elem<nb_elem_tot; elem++)
468 {
469 int rang=rang_elem_non_std(elem);
470
471 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
472 face(facei_loc)=elem_faces(elem,facei_loc);
473
474 //RAZ des flux des elements a modifier
475 if (is_element_for_upwinding_(elem))
476 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
477 for (int facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
478 for (int comp=0; comp<nb_comp; comp++)
479 Fij(elem,facei_loc,facej_loc,comp)=0.;
480
481 if (is_element_for_upwinding_(elem))
482 {
483 if (rang==-1)
484 for (int fa7=0; fa7<nfa7; fa7++)
485 {
486 int facei_loc=KEL(0,fa7);
487 int facej_loc=KEL(1,fa7);
488
489 int facei=face(facei_loc);
490 int facej=face(facej_loc);
491
492 //REMARQUE : ON N'OUBLIE PAS QUE LES NORMALES FA7 SONT ORIENTeES DE FACE_LOCI VERS FACE_LOCJ PAR DEFAUT
493 //LE PRODUIT SCALAIRE PSC_M EST DONC CALCULe POUR OBTENIR LE FLUX REcU PAR LA FACE FACE_LOCI
494 //ET EMIS PAR LA FACE_LOCJ
495 //SI L'ON VEUT OBTENIR LE FLUX REcU PAR LA FACE FACE_LOCJ ET EMIS PAR LA FACE_LOCI, IL FAUT MULTIPLIER LE
496 //PRODUIT SCALAIRE PSC_M PAR -1 POUR IMPLICITEMENT ReORIENTER LES NORMALES FA7 DE FACE_LOCJ VERS
497 //FACE_LOCI
498 //MEME REMARQUE POUR LE SCALAIRE "FLUX"
499 const double psc_m = Kij(elem,facei_loc,facej_loc);
500
501 int face_amont, dir;
502 if (psc_m>=0.)
503 face_amont=facei,dir=0;
504 else
505 face_amont=facej,dir=1;
506
507 for (int comp=0; comp<nb_comp; comp++)
508 {
509 //Calcul de la valeur de l'inconnue aux points d'integration de la formule du point milieu
510 double inco_m=transporteVect[face_amont*nb_comp+comp];
511 for(int dim=0; dim<nb_dim; dim++)
512 inco_m+=gradient_elem_(elem,comp,dim)*vecteur_face_facette(elem,fa7,dim,dir);
513
514 //Calcul du flux final : formule du point milieu
515 double flux=inco_m*psc_m;
516
517 //Fij(elem,facei_loc,facej_loc,comp) est le flux recu par la facei_loc dans elem
518 //envoye par la face facej_loc dans elem pour la composante comp de l'inconnue
519 Fij(elem,facei_loc,facej_loc,comp)=flux-psc_m*transporteVect[facei*nb_comp+comp];
520 Fij(elem,facej_loc,facei_loc,comp)=psc_m*transporteVect[facej*nb_comp+comp]-flux;
521 }
522 }
523 }
524 else //rang!=-1
525 {
526 for (int fa7=0; fa7<nfa7; fa7++)
527 {
528 int facei_loc=KEL(0,fa7);
529 int facej_loc=KEL(1,fa7);
530
531 int facei=face(facei_loc);
532 int facej=face(facej_loc);
533
534 //REMARQUE : ON N'OUBLIE PAS QUE LES NORMALES FA7 SONT ORIENTeES DE FACE_LOCI VERS FACE_LOCJ PAR DEFAUT
535 //LE PRODUIT SCALAIRE PSC_M EST DONC CALCULe POUR OBTENIR LE FLUX REcU PAR LA FACE FACE_LOCI
536 //ET EMIS PAR LA FACE_LOCJ
537 //SI L'ON VEUT OBTENIR LE FLUX REcU PAR LA FACE FACE_LOCJ ET EMIS PAR LA FACE_LOCI, IL FAUT MULTIPLIER LE
538 //PRODUIT SCALAIRE PSC_M PAR -1 POUR IMPLICITEMENT ReORIENTER LES NORMALES FA7 DE FACE_LOCJ VERS
539 //FACE_LOCI
540 //MEME REMARQUE POUR LE SCALAIRE "FLUX"
541 const double psc_m = Kij(elem,facei_loc,facej_loc);
542
543 int face_amont, dir;
544 if (psc_m>=0.)
545 face_amont=facei,dir=0;
546 else
547 face_amont=facej,dir=1;
548
549 for (int comp=0; comp<nb_comp; comp++)
550 {
551 //Calcul de la valeur de l'inconnue aux points d'integration de la formule du point milieu
552 double inco_m=transporteVect[face_amont*nb_comp+comp];
553 for(int dim=0; dim<nb_dim; dim++)
554 inco_m+=gradient_elem_(elem,comp,dim)*vecteur_face_facette_Cl(rang,fa7,dim,dir);
555
556 //Calcul du flux final : formule d'integration du point milieu
557 double flux=inco_m*psc_m;
558
559 //Fij(elem,facei_loc,facej_loc,comp) est le flux recu par la facei_loc dans elem
560 //envoye par la face facej_loc dans elem pour la composante comp de l'inconnue
561 Fij(elem,facei_loc,facej_loc,comp)=flux-psc_m*transporteVect[facei*nb_comp+comp];
562 Fij(elem,facej_loc,facei_loc,comp)=psc_m*transporteVect[facej*nb_comp+comp]-flux;
563 }
564 }
565 }
566 }//for sur elem
567}
568
569void Op_Conv_Muscl_New_VEF_Face::remplir_fluent() const
570{
571 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
572 const int nb_faces = domaine_VEF.nb_faces();
573 const int dim = dimension;
574 // calcul de la CFL.
575 // On remet a zero le tableau qui sert pour
576 // le calcul du pas de temps de stabilite
577 CDoubleTabView face_normales=domaine_VEF.face_normales().view_ro();
578 CDoubleTabView velocity=vitesse_->valeurs().view_ro();
579 DoubleArrView fluent = fluent_.view_wo();
580 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
581 Kokkos::RangePolicy<>(0, nb_faces), KOKKOS_LAMBDA(
582 const int num_face)
583 {
584 double psc=0.;
585 for (int i=0; i<dim; i++)
586 psc+=velocity(num_face,i)*face_normales(num_face,i);
587 fluent(num_face)=std::fabs(psc);
588 });
589 end_gpu_timer(__KERNEL_NAME__);
590}
591
592
593DoubleTab& Op_Conv_Muscl_New_VEF_Face::ajouter(const DoubleTab& transporte_2,
594 DoubleTab& resu) const
595{
596 DoubleTrav sauv(resu);
597 sauv = resu;
598 resu = 0;
599
600 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
601 const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
602 const DoubleTab& vitesse_2=la_vitesse.valeurs();
603 const DoubleVect& porosite_face = equation().milieu().porosite_face();
604
605 const IntTab& elem_faces = domaine_VEF.elem_faces();
606
607 const int marq=phi_u_transportant(equation());
608 const int nb_faces_elem=elem_faces.dimension(1);
609 assert(nb_faces_elem==(dimension+1));
610 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
611 const int nfa7 = domaine_VEF.type_elem().nb_facette();
612
613 int nb_comp=1;
614 if(resu.nb_dim()!=1)
615 nb_comp=resu.dimension(1);
616
617 DoubleTrav transporte_;
618 DoubleTrav vitesse_face_;
619 // soit on a transporte=phi*transporte_ et vitesse=vitesse_
620 // soit transporte=transporte_ et vitesse=phi*vitesse_
621 // cela depend si on transporte avec phi u ou avec u.
622 const DoubleTab& velocity=modif_par_porosite_si_flag(vitesse_2,vitesse_face_,marq,porosite_face);
623
624 DoubleTrav Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
625 DoubleTrav Fij(nb_elem_tot,nb_faces_elem,nb_faces_elem,nb_comp);
626 DoubleTrav Cij(nb_elem_tot,nfa7);
627 DoubleTrav Sij(nb_elem_tot,nfa7);
628 DoubleTrav Sij2(nb_elem_tot,nfa7);
629
630 //Pour tenir compte des conditions de Neumann sortie libre
631
632 // soit on a transporte=phi*transporte_ et vitesse=vitesse_
633 // soit transporte=transporte_ et vitesse=phi*vitesse_
634 // cela depend si on transporte avec phi u ou avec u.
635 const DoubleTab& transporte=modif_par_porosite_si_flag(transporte_2,transporte_,!marq,porosite_face);
636
637 //Initialisation du tableau flux_bords_ pour le calcul des pertes de charge
638 flux_bords_.resize(domaine_VEF.nb_faces_bord(),nb_comp);
639 calculer_flux_bords(Kij,velocity,transporte);
640
641 calculer_coefficients_operateur_centre(Kij,Cij,Sij,Sij2,nb_comp,velocity);
642 calculer_flux_operateur_centre(Fij,Kij,Cij,Sij,Sij2,nb_comp,velocity,transporte);
643 if (old_centered_) modifier_flux_operateur_centre(Fij,Kij,Cij,Sij,Sij2,nb_comp,velocity,transporte); // default false
644 if (centered_) ajouter_operateur_centre(Kij,Fij,transporte,resu); // default true
645 if (upwind_) ajouter_diffusion(Kij,Fij,transporte,resu); // default true
646 if (stabilized_) ajouter_antidiffusion(Kij,Fij,transporte,resu); // default true
647
648 mettre_a_jour_pour_periodicite(Kij,transporte,resu);
649
650 resu*=-1.;//car l'operateur de convection est construit en tant que terme source
651 resu+=sauv;
652
653 modifier_flux(*this);
654
655 return resu;
656}
657
658//ALGO TRES GROSSIER MAIS FONCTIONNE FACILEMENT
659double Op_Conv_Muscl_New_VEF_Face::calculer_dt_stab() const
660{
661 if (!facsec_auto_)
662 return Op_Conv_VEF_Face::calculer_dt_stab(); // Default true
663 else
664 {
665 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
666 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
667 const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
668 const DoubleTab& vitesse_2=la_vitesse.valeurs();
669 const DoubleVect& porosite_face = equation().milieu().porosite_face();
670
671 const IntTab& elem_faces = domaine_VEF.elem_faces();
672
673 const int marq=phi_u_transportant(equation());
674 const int nb_faces_elem=elem_faces.dimension(1);
675 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
676 const int nfa7 = domaine_VEF.type_elem().nb_facette();
677 const int nb_faces_tot = domaine_VEF.nb_faces_tot();
678 const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
679
680 int nb_comp=1;
681 if(equation().inconnue().valeurs().nb_dim()!=1)
682 nb_comp=equation().inconnue().valeurs().dimension(1);
683
684 double dt_stab = Op_Conv_VEF_Face::calculer_dt_stab();
685 double security_coeff=0.95;
686 double dt_corrector=-1.;
687
688 // soit on a transporte=phi*transporte_ et vitesse=vitesse_
689 // soit transporte=transporte_ et vitesse=phi*vitesse_
690 // cela depend si on transporte avec phi u ou avec u.
691 ToDo_Kokkos("critical DoubleTab -> DoubleTrav");
692 DoubleTab vitesse_face_;
693 const DoubleTab& velocity=modif_par_porosite_si_flag(vitesse_2,vitesse_face_,marq,porosite_face);
694
695 DoubleTab Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
696 DoubleTab Cij(nb_elem_tot,nfa7);
697 DoubleTab Sij(nb_elem_tot,nfa7);
698 DoubleTab Sij2;
699 if (Objet_U::dimension==3)
700 {
701 Sij2.resize(nb_elem_tot,nfa7);
702 Sij2=0.;
703 }
704 calculer_coefficients_operateur_centre(Kij,Cij,Sij,Sij2,nb_comp,velocity);
705
706 //Debut du calcul
707 IntTab plus_tab(nb_faces_tot);
708 ToDo_Kokkos("critical");
709 for (int elem=0; elem<nb_elem_tot; elem++)
710 for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
711 {
712 int facei=elem_faces(elem,facei_loc);
713 for (int facej_loc=facei_loc+1; facej_loc<nb_faces_elem; facej_loc++)
714 {
715 double kij=Kij(elem,facei_loc,facej_loc);
716 int facej=elem_faces(elem,facej_loc);
717
718 if (kij>=0.)
719 plus_tab(facei)+=1;
720 else
721 plus_tab(facej)+=1;
722 }
723 }
724 for (int n_bord=0; n_bord<nb_bord; n_bord++)
725 {
726 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
727 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
728 int num2=le_bord.nb_faces_tot();
729
730 if (sub_type(Periodique,la_cl.valeur()))
731 {
732 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
733
734 for (int ind_face=0; ind_face<num2; ind_face++)
735 {
736 int facei=le_bord.num_face(ind_face);
737 int ind_face_associee=la_cl_perio.face_associee(ind_face);
738 int faceiAss=le_bord.num_face(ind_face_associee);
739
740 if (facei<faceiAss)
741 {
742 plus_tab[faceiAss]+=plus_tab[facei];
743 plus_tab[facei]=plus_tab[faceiAss];
744 }
745
746 }//fin du for sur "ind_face"
747
748 }//fin du if sur "Periodique"
749
750 }//fin du for sur "n_bord"
751 int max_int=plus_tab.mp_max_vect();
752 max_int=(int)Process::mp_max(max_int);
753
754 dt_corrector=1./(1+max_limiteur_*max_int);
755 dt_corrector*=security_coeff;
756 dt_stab*=dt_corrector;
757
758 Op_Conv_VEF_base& op = ref_cast_non_const(Op_Conv_VEF_base,*this);
759 op.fixer_dt_stab_conv(dt_stab);
760
761 return dt_stab;
762 }
763}
764
765void Op_Conv_Muscl_New_VEF_Face::calculer_flux_bords(const DoubleTab& Kij, const DoubleTab& tab_velocity, const DoubleTab& tab_transporte) const
766{
767 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
768 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
769 const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
770
771 int nb_comp=1;
772 if (tab_transporte.nb_dim()!=1) nb_comp=tab_transporte.dimension(1);
773 int nb_dim = Objet_U::dimension;
774 flux_bords_ = 0;
775 // Kokkos: example of kernels launched concurently on each boundary. The Kokkos::fence() is called in end_gpu_timer
776 CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
777 CDoubleArrView transporte = static_cast<const DoubleVect&>(tab_transporte).view_ro();
778 CDoubleTabView velocity = tab_velocity.view_ro();
779 DoubleTabView flux_bords = flux_bords_.view_wo();
780 for (int n_bord=0; n_bord<nb_bord; n_bord++)
781 {
782 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
783 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
784 int num2 = le_bord.nb_faces();//il ne faut boucler que sur les faces reelles ici
785 CIntArrView le_bord_num_face = le_bord.num_face().view_ro();
786 if ( sub_type(Dirichlet_homogene,la_cl.valeur()) )
787 {
788 //On ne calcule pas le flux aux bords de Dirichlet_homogene
789 }//fin du if sur "Dirichlet"
790 else if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
791 {
792 const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
793 CDoubleTabView val_ext = la_sortie_libre.val_ext().view_ro();
794 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), num2, KOKKOS_LAMBDA(const int ind_face)
795 {
796 int facei = le_bord_num_face(ind_face);
797
798 double psc=0.;
799 for (int dim=0; dim<nb_dim; dim++)
800 psc-=velocity(facei,dim)*face_normales(facei,dim);
801
802 for (int dim=0; dim<nb_comp; dim++)
803 {
804 double val = psc < 0. ? transporte[facei*nb_comp+dim] : (nb_comp==1 ? val_ext(ind_face,0) : val_ext(ind_face,dim));
805 flux_bords(facei,dim) = psc*val;
806 }
807 });
808 end_gpu_timer(__KERNEL_NAME__);
809 }
810 else if ( sub_type(Dirichlet,la_cl.valeur())
811 || sub_type(Neumann,la_cl.valeur())
812 || sub_type(Neumann_homogene,la_cl.valeur())
813 || sub_type(Symetrie,la_cl.valeur())
814 || sub_type(Echange_impose_base,la_cl.valeur())
815 )
816 {
817 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), num2, KOKKOS_LAMBDA(const int ind_face)
818 {
819 int facei = le_bord_num_face(ind_face);
820
821 double psc=0.;
822 for (int dim=0; dim<nb_dim; dim++)
823 psc-=velocity(facei,dim)*face_normales(facei,dim);
824
825 for (int dim=0; dim<nb_comp; dim++)
826 flux_bords(facei,dim)=psc*transporte[facei*nb_comp+dim];
827 });
828 end_gpu_timer(__KERNEL_NAME__);
829 }//fin du if sur "Neumann", "Neumann_homogene", "Symetrie", "Echange_impose_base"
830 else if (sub_type(Periodique,la_cl.valeur()))
831 {
832 const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
833 CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
834 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), num2, KOKKOS_LAMBDA(const int ind_face)
835 {
836 int facei = le_bord_num_face(ind_face);
837 int ind_face_voisine = face_associee(ind_face);
838 int facei_voisine = le_bord_num_face(ind_face_voisine);
839
840 double psc=0.;
841 for (int dim=0; dim<nb_dim; dim++)
842 psc-=velocity(facei,dim)*face_normales(facei,dim);
843
844 for (int dim=0; dim<nb_comp; dim++)
845 {
846 double flux = psc*transporte[facei*nb_comp+dim];
847 Kokkos::atomic_add(&flux_bords(facei,dim), 0.5 * flux);
848 Kokkos::atomic_add(&flux_bords(facei_voisine,dim), -0.5 * flux);
849 }
850 });
851 end_gpu_timer(__KERNEL_NAME__);
852 }
853 else
854 {
855 Cerr << "Erreur Op_Conv_Muscl_New_VEF_Face::calculer_flux_bords()" << finl;
856 Cerr << "Condition aux limites " << la_cl.que_suis_je() << " non codee." << finl;
857 Cerr << "Sortie du programme." << finl;
858 Process::exit();
859 }//fin du else sur les autres conditions aux limites
860 }
861}
862
863DoubleTab&
864Op_Conv_Muscl_New_VEF_Face::ajouter_operateur_centre(const DoubleTab& tab_Kij, const DoubleTab& tab_Fij, const DoubleTab& tab_transporte, DoubleTab& tab_resu) const
865{
866 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
867 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
868 const int nfa7 = domaine_VEF.type_elem().nb_facette();
869 const int premiere_face_int = domaine_VEF.premiere_face_int();
870 const int nb_faces = domaine_VEF.nb_faces();
871
872 int nb_comp=1;
873 if (tab_transporte.nb_dim()!=1) nb_comp=tab_transporte.dimension(1);
874
875 //Faces internes
876 {
877 CIntTabView KEL = domaine_VEF.type_elem().KEL().view_ro();
878 CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
879 CDoubleTabView4 Fij = tab_Fij.view_ro<4>();
880 DoubleArrView resuV = static_cast<DoubleVect&>(tab_resu).view_rw();
881 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
882 Kokkos::RangePolicy<>(0, nb_elem_tot), KOKKOS_LAMBDA(
883 const int elem)
884 {
885 for (int fa7 = 0; fa7 < nfa7; fa7++)
886 {
887 int facei_loc = KEL(0, fa7);
888 int facej_loc = KEL(1, fa7);
889
890 int facei = elem_faces(elem, facei_loc);
891 int facej = elem_faces(elem, facej_loc);
892
893 for (int dim = 0; dim < nb_comp; dim++)
894 {
895 double fij = Fij(elem, facei_loc, facej_loc, dim);
896 double fji = Fij(elem, facej_loc, facei_loc, dim);
897
898 int ligne = facei * nb_comp + dim;
899 int colonne = facej * nb_comp + dim;
900
901 Kokkos::atomic_add(&resuV[ligne], fij);
902 Kokkos::atomic_add(&resuV[colonne], fji);
903 }
904 }
905 });
906 end_gpu_timer(__KERNEL_NAME__);
907 }
908
909 if (old_centered_)
910 {
911 CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
912 CDoubleArrView transporteV = static_cast<const DoubleVect&>(tab_transporte).view_ro();
913 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
914 CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
915 DoubleArrView resuV = static_cast<DoubleVect&>(tab_resu).view_rw();
916 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
917 Kokkos::RangePolicy<>(premiere_face_int, nb_faces), KOKKOS_LAMBDA(
918 const int facei)
919 {
920 //Premier voisin
921 int elem = face_voisins(facei, 0);
922 int facei_loc = num_fac_loc(facei, 0);
923 double fij = Kij(elem, facei_loc, facei_loc);
924
925 //Deuxieme voisin
926 elem = face_voisins(facei, 1);
927 facei_loc = num_fac_loc(facei, 1);
928 fij += Kij(elem, facei_loc, facei_loc);
929
930 for (int dim = 0; dim < nb_comp; dim++)
931 {
932 int ligne = facei * nb_comp + dim;
933 resuV[ligne] -= fij * transporteV[ligne];
934 }
935 });
936 end_gpu_timer(__KERNEL_NAME__);
937 }
938
939 //Faces de bord
940 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
941
942 const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
943
944 for (int n_bord=0; n_bord<nb_bord; n_bord++)
945 {
946 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
947 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
948 int num2=le_bord.nb_faces();
949
950 if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
951 {
952 const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
953 CIntArrView num_face = le_bord.num_face().view_ro();
954 CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
955 CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
956 CDoubleArrView transporteV = static_cast<const DoubleVect&>(tab_transporte).view_ro();
957 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
958 CDoubleTabView val_ext = la_sortie_libre.val_ext().view_ro();
959 DoubleArrView resuV = static_cast<DoubleVect&>(tab_resu).view_rw();
960 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
961 Kokkos::RangePolicy<>(0, num2), KOKKOS_LAMBDA(
962 const int ind_face)
963 {
964 int facei=num_face(ind_face);
965 int facei_loc=num_fac_loc(facei,0);
966 int elem=face_voisins(facei,0);
967 double psc=Kij(elem,facei_loc,facei_loc);//ATTENTION : SUPPOSE QU'ON EST a DIVERGENCE NULLE
968
969 if (psc>=0.)
970 {
971 //On ne fait rien car on en a deja tenu compte dans l'operateur centre
972 }
973 else
974 for (int dim=0; dim<nb_comp; dim++)
975 {
976 //On modifie car on a tenu compte dans l'operateur centre de psc*tansporteV[ligne]
977 int ligne=facei*nb_comp+dim;
978 //PL: Wow, fix a bug, ind_face and not facei !
979 //All use of val_ext() in the code has the same mistake !
980 //But as the val_ext field is uniform in general and psc>=0, we have never seen the bug...
981 //ToDo: fix all the val_ext(face - num1, k) -> val_ext(ind_face, k);
982 //resuV[ligne]+=psc*(val_ext(facei,dim)-transporteV[ligne]);
983 resuV[ligne]+=psc*(val_ext(ind_face,dim)-transporteV[ligne]);
984 }
985
986 });//fin du for sur "face_i"
987 end_gpu_timer(__KERNEL_NAME__);
988
989 }//fin du if sur "Neumann"
990
991 }//fin du for sur "n_bord"
992
993
994 //Retour du resultat
995 return tab_resu;
996}
997
998DoubleTab&
999Op_Conv_Muscl_New_VEF_Face::ajouter_diffusion(const DoubleTab& tab_Kij,const DoubleTab& tab_Fij, const DoubleTab& tab_transporte, DoubleTab& tab_resu) const
1000{
1001 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
1002 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
1003 const int nfa7 = domaine_VEF.type_elem().nb_facette();
1004
1005 int nb_comp=1;
1006 if (tab_transporte.nb_dim()!=1) nb_comp=tab_transporte.dimension(1);
1007 double alpha = alpha_;
1008 //Pour les faces internes
1009 CIntTabView KEL=domaine_VEF.type_elem().KEL().view_ro();
1010 CIntTabView elem_faces=domaine_VEF.elem_faces().view_ro();
1011 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
1012 CDoubleTabView4 Fij = tab_Fij.view_ro<4>();
1013 DoubleArrView resuV = static_cast<DoubleVect&>(tab_resu).view_rw();
1014 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
1015 Kokkos::RangePolicy<>(0, nb_elem_tot), KOKKOS_LAMBDA(
1016 const int elem)
1017 {
1018 for (int fa7 = 0; fa7 < nfa7; fa7++)
1019 {
1020 int facei_loc = KEL(0, fa7);
1021 int facej_loc = KEL(1, fa7);
1022
1023 int facei = elem_faces(elem, facei_loc);
1024 int facej = elem_faces(elem, facej_loc);
1025
1026 double psc = Kij(elem, facei_loc, facej_loc);
1027
1028 for (int dim = 0; dim < nb_comp; dim++)
1029 {
1030 double fij = alpha * Fij(elem, facei_loc, facej_loc, dim);
1031 double fji = alpha * Fij(elem, facej_loc, facei_loc, dim);
1032
1033 int ligne = facei * nb_comp + dim;
1034 int colonne = facej * nb_comp + dim;
1035
1036 if (psc >= 0.)
1037 {
1038 Kokkos::atomic_sub(&resuV[ligne], fij);
1039 Kokkos::atomic_add(&resuV[colonne], fij);
1040 }
1041 else
1042 {
1043 Kokkos::atomic_add(&resuV[ligne], fji);
1044 Kokkos::atomic_sub(&resuV[colonne], fji);
1045 }
1046 }
1047 }
1048 });
1049 end_gpu_timer(__KERNEL_NAME__);
1050
1051 //Pour les faces de bord :
1052 //ON N'A RIEN a FAIRE
1053
1054 //Retour du resultat
1055 return tab_resu;
1056}
1057
1058DoubleTab&
1059Op_Conv_Muscl_New_VEF_Face::ajouter_antidiffusion(const DoubleTab& Kij, const DoubleTab& Fij,const DoubleTab& transporte, DoubleTab& resu) const
1060{
1061 switch(version_)
1062 {
1063 case 1 :
1064
1065 return ajouter_antidiffusion_v1(Kij,Fij,transporte,resu);
1066
1067 case 2 :
1068
1069 return ajouter_antidiffusion_v2(Kij,Fij,transporte,resu);
1070
1071 default :
1072 Cerr<<"Error in Op_Conv_Muscl_New_VEF_Face::ajouter_antidiffusion()"<<finl;
1073 Cerr<<"Version number "<<version_<<" of antidiffusive operator does not exist"<<finl;
1074 Process::exit();
1075 break;
1076 }
1077
1078 return resu;
1079}
1080
1081DoubleTab&
1082Op_Conv_Muscl_New_VEF_Face::ajouter_antidiffusion_v2(const DoubleTab& tab_Kij, const DoubleTab& tab_Fij,const DoubleTab& tab_transporte, DoubleTab& tab_resu) const
1083{
1084 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
1085 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
1086 const int nfa7 = domaine_VEF.type_elem().nb_facette();
1087
1088 int nb_comp=1;
1089 if (tab_transporte.nb_dim()!=1) nb_comp=tab_transporte.dimension(1);
1090
1091 double alpha = alpha_;
1092 //Pour les faces internes
1093 CIntTabView KEL=domaine_VEF.type_elem().KEL().view_ro();
1094 CIntTabView elem_faces=domaine_VEF.elem_faces().view_ro();
1095 CDoubleArrView transporteV = static_cast<const DoubleVect&>(tab_transporte).view_ro();
1096 CIntTabView face_voisins=domaine_VEF.face_voisins().view_ro();
1097 CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
1098 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
1099 CDoubleTabView4 Fij = tab_Fij.view_ro<4>();
1100 CIntArrView is_dirichlet_faces = is_dirichlet_faces_.view_ro();
1101 DoubleArrView resuV = static_cast<DoubleVect&>(tab_resu).view_rw();
1102 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
1103 Kokkos::RangePolicy<>(0, nb_elem_tot), KOKKOS_LAMBDA(
1104 const int elem)
1105 {
1106 for (int dim=0; dim<nb_comp; dim++)
1107 for (int fa7 = 0; fa7 < nfa7; fa7++)
1108 {
1109 int facei_loc = KEL(0, fa7);
1110 int facej_loc = KEL(1, fa7);
1111
1112 int facei = elem_faces(elem, facei_loc);
1113 int facej = elem_faces(elem, facej_loc);
1114
1115 double kij = Kij(elem, facei_loc, facej_loc);
1116
1117 int face = kij >= 0. ? facei : facej;
1118 double P_plus = 0;
1119 double P_moins = 0;
1120 double Q_plus = 0;
1121 double Q_moins = 0;
1122 calculer_senseur(Kij, Fij, transporteV, dim, nb_comp, face, elem_faces, face_voisins, num_fac_loc, P_plus,
1123 P_moins, Q_plus, Q_moins);
1124
1125 double fij = alpha * Fij(elem, facei_loc, facej_loc, dim);
1126 double fji = alpha * Fij(elem, facej_loc, facei_loc, dim);
1127
1128 int ligne = facei * nb_comp + dim;
1129 int colonne = facej * nb_comp + dim;
1130
1131 double fij_low = transporteV[colonne] - transporteV[ligne];
1132 fij_low *= kij;
1133 double fji_low = fij_low;
1134
1135 double R;
1136 if (kij >= 0.) //facei amont
1137 {
1138 if (fij >= 0.) R = (std::fabs(P_plus) < DMINFLOAT) ? 0. : Q_plus / P_plus;
1139 else R = (std::fabs(P_moins) < DMINFLOAT) ? 0. : Q_moins / P_moins;
1140
1141
1142 R = minmod(R);
1143 R *= fij;
1144
1145 double tmp;
1146 if (!is_dirichlet_faces(facej))
1147 tmp = optimum(R, fji_low);
1148 else
1149 tmp = R;
1150
1151 Kokkos::atomic_add(&resuV[ligne], tmp);
1152 Kokkos::atomic_sub(&resuV[colonne], tmp);
1153 }
1154 else //facej amont
1155 {
1156 if (fji <= 0.) R = (std::fabs(P_moins) < DMINFLOAT) ? 0. : Q_moins / P_moins;
1157 else R = (std::fabs(P_plus) < DMINFLOAT) ? 0. : Q_plus / P_plus;
1158
1159 R = minmod(R);
1160 R *= fji;
1161
1162 double tmp;
1163 if (!is_dirichlet_faces(facei))
1164 tmp = optimum(R, fij_low);
1165 else
1166 tmp = R;
1167
1168 Kokkos::atomic_sub(&resuV[ligne], tmp);
1169 Kokkos::atomic_add(&resuV[colonne], tmp);
1170 }
1171 }
1172 });
1173 end_gpu_timer(__KERNEL_NAME__);
1174
1175 //Pour les faces de bord
1176 //IL N'Y A RIEN a FAIRE TOUT EST FAIT DANS LA FONCTION AJOUTER_DIFFUSION
1177 //QUI EST PARFAITEMENT COMPLeTER PAR LA FONCTION AJOUTER_ANTIDIFFUSION
1178
1179 //Retour du resultat
1180 return tab_resu;
1181}
1182
1183DoubleTab&
1184Op_Conv_Muscl_New_VEF_Face::ajouter_antidiffusion_v1(const DoubleTab& tab_Kij, const DoubleTab& tab_Fij,const DoubleTab& tab_transporte, DoubleTab& tab_resu) const
1185{
1186 const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
1187 const int nb_elem_tot=domaine_VEF.nb_elem_tot();
1188 const int nfa7 = domaine_VEF.type_elem().nb_facette();
1189
1190 int nb_comp=1;
1191 if (tab_transporte.nb_dim()!=1) nb_comp=tab_transporte.dimension(1);
1192
1193 double alpha = alpha_;
1194 //Pour les faces internes
1195 CIntTabView KEL=domaine_VEF.type_elem().KEL().view_ro();
1196 CIntTabView elem_faces=domaine_VEF.elem_faces().view_ro();
1197 CDoubleArrView transporteV = static_cast<const DoubleVect&>(tab_transporte).view_ro();
1198 CIntTabView face_voisins=domaine_VEF.face_voisins().view_ro();
1199 CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
1200 CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
1201 CDoubleTabView4 Fij = tab_Fij.view_ro<4>();
1202 CIntArrView is_dirichlet_faces = is_dirichlet_faces_.view_ro();
1203 DoubleArrView resuV = static_cast<DoubleVect&>(tab_resu).view_rw();
1204 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
1205 Kokkos::RangePolicy<>(0, nb_elem_tot), KOKKOS_LAMBDA(
1206 const int elem)
1207 {
1208 for (int dim=0; dim<nb_comp; dim++)
1209 for (int fa7 = 0; fa7 < nfa7; fa7++)
1210 {
1211 int facei_loc = KEL(0, fa7);
1212 int facej_loc = KEL(1, fa7);
1213
1214 int facei = elem_faces(elem, facei_loc);
1215 int facej = elem_faces(elem, facej_loc);
1216
1217 double kij = Kij(elem, facei_loc, facej_loc);
1218 double Pi_plus = 0;
1219 double Pi_moins = 0;
1220 double Qi_plus = 0;
1221 double Qi_moins = 0;
1222 calculer_senseur(Kij, Fij, transporteV, dim, nb_comp, facei, elem_faces, face_voisins, num_fac_loc,
1223 Pi_plus,
1224 Pi_moins, Qi_plus, Qi_moins);
1225 double Pj_plus = 0;
1226 double Pj_moins = 0;
1227 double Qj_plus = 0;
1228 double Qj_moins = 0;
1229 calculer_senseur(Kij, Fij, transporteV, dim, nb_comp, facej, elem_faces, face_voisins, num_fac_loc,
1230 Pj_plus,
1231 Pj_moins, Qj_plus, Qj_moins);
1232
1233 double fij = alpha * Fij(elem, facei_loc, facej_loc, dim);
1234 double fji = alpha * Fij(elem, facej_loc, facei_loc, dim);
1235
1236 int ligne = facei * nb_comp + dim;
1237 int colonne = facej * nb_comp + dim;
1238
1239 if (kij >= 0.)
1240 {
1241 double Ri, Rj;
1242 //Face amont : facei
1243 if (fij >= 0.)
1244 {
1245 Ri = (std::fabs(Pi_plus) < DMINFLOAT) ? 0. : Qi_plus / Pi_plus;
1246 Rj = (std::fabs(Pj_moins) < DMINFLOAT) ? 0. : Qj_moins /
1247 Pj_moins;//car fji=-fij
1248 }
1249 else
1250 {
1251 Ri = (std::fabs(Pi_moins) < DMINFLOAT) ? 0. : Qi_moins / Pi_moins;
1252 Rj = (std::fabs(Pj_plus) < DMINFLOAT) ? 0. : Qj_plus /
1253 Pj_plus;//car fji=-fij
1254 }
1255
1256 if (is_dirichlet_faces(facej))
1257 Rj = DMAXFLOAT;//on n'a pas besoin de prendre le min quand il y a une face de Dirichlet
1258 double R = (Ri <= Rj) ? Ri : Rj;
1259 R = minmod(R);
1260
1261 double tmp = R * fij;
1262 Kokkos::atomic_add(&resuV[ligne], tmp);
1263 Kokkos::atomic_sub(&resuV[colonne], tmp);
1264 }
1265 else
1266 {
1267 double Ri, Rj;
1268 //Face amont : facej
1269 if (fji <= 0.)
1270 {
1271 Rj = (std::fabs(Pj_moins) < DMINFLOAT) ? 0. : Qj_moins / Pj_moins;
1272 Ri = (std::fabs(Pi_plus) < DMINFLOAT) ? 0. : Qi_plus / Pi_plus;
1273 }
1274 else
1275 {
1276 Rj = (std::fabs(Pj_plus) < DMINFLOAT) ? 0. : Qj_plus / Pj_plus;
1277 Ri = (std::fabs(Pi_moins) < DMINFLOAT) ? 0. : Qi_moins / Pi_moins;
1278 }
1279
1280 if (is_dirichlet_faces(facei))
1281 Ri = DMAXFLOAT;//on n'a pas besoin de prendre le min quand il y a une face de Dirichlet
1282 double R = (Ri <= Rj) ? Ri : Rj;
1283 R = minmod(R);
1284
1285 double tmp = R * fji;
1286 Kokkos::atomic_sub(&resuV[ligne], tmp);
1287 Kokkos::atomic_add(&resuV[colonne], tmp);
1288 }
1289 }
1290 });
1291 end_gpu_timer(__KERNEL_NAME__);
1292
1293 //Pour les faces de bord
1294 //IL N'Y A RIEN a FAIRE TOUT EST FAIT DANS LA FONCTION AJOUTER_DIFFUSION
1295 //QUI EST PARFAITEMENT COMPLeTER PAR LA FONCTION AJOUTER_ANTIDIFFUSION
1296
1297 //Retour du resultat
1298 // resuV+=antidiff;
1299 return tab_resu;
1300}
1301
1302//ATTENTION : suppose les parametres P_plus, P_moins, Q_plus, Q_moins nuls en entree
1303KOKKOS_INLINE_FUNCTION void
1304Op_Conv_Muscl_New_VEF_Face::calculer_senseur(CDoubleTabView3 Kij, CDoubleTabView4 Fij, CDoubleArrView transporteV,
1305 const int dim, const int nb_comp, const int face_i,
1306 CIntTabView elem_faces, CIntTabView face_voisins, CIntTabView num_fac_loc,
1307 double& P_plus, double& P_moins,
1308 double& Q_plus, double& Q_moins) const
1309{
1310 const int nb_faces_elem=(int)elem_faces.extent(1);
1311 for (int elem_voisin=0; elem_voisin<2; elem_voisin++)
1312 {
1313 int elem = face_voisins(face_i,elem_voisin);
1314 if (elem!=-1)
1315 {
1316 int face_i_loc = num_fac_loc(face_i,elem_voisin);
1317#ifndef TRUST_USE_GPU
1318 assert(face_i_loc>=0);
1319 assert(face_i_loc<nb_faces_elem);
1320#endif
1321 //On travaille sur les faces de "elem"
1322 for (int face_k_loc=0; face_k_loc<nb_faces_elem; face_k_loc++)
1323 {
1324 int face_k=elem_faces(elem,face_k_loc);
1325 double kik=Kij(elem,face_i_loc,face_k_loc);
1326 //
1327 //Calcul des variables intermediaires
1328 //
1329 double inci=transporteV[face_i*nb_comp+dim];
1330 double inck=transporteV[face_k*nb_comp+dim];
1331
1332 double fik_low=kik*(inck-inci);
1333 double fik_high=Fij(elem,face_i_loc,face_k_loc,dim);
1334
1335 // Codage optimise:
1336 if (kik<0)
1337 {
1338 if (fik_low>0) Q_plus+=fik_low;
1339 else Q_moins+=fik_low;
1340 }
1341 else
1342 {
1343 if (fik_high>0) P_plus+=fik_high;
1344 else P_moins+=fik_high;
1345 }
1346#ifndef TRUST_USE_GPU
1347 assert(P_plus>=0);
1348 assert(Q_plus>=0);
1349 assert(P_moins<=0);
1350 assert(Q_moins<=0);
1351#endif
1352 //
1353 //Fin du calcul des variables intermediaires
1354 //
1355 }//fin du for sur "face_k_loc"
1356 }//fin du if sur "elem!=-1"
1357 }//fin du for sur "elem_voisin"
1358}
1359
1360void Op_Conv_Muscl_New_VEF_Face::mettre_a_jour_pour_periodicite(const DoubleTab& Kij, const DoubleTab& transporte, DoubleTab& tab_resu) const
1361{
1362 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1363
1364 const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
1365 const int nb_comp = (tab_resu.nb_dim()==1) ? 1 : tab_resu.dimension(1);
1366
1367 //Faces de bord
1368 int old_centered = old_centered_;
1369 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1370 {
1371 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1372 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1373 int num2=le_bord.nb_faces();
1374
1375 if (sub_type(Periodique,la_cl.valeur()))
1376 {
1377 const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
1378 if (old_centered) Process::exit("Not available for periodicity."); // See code below commented, cause not tested and no datafile option to test it !
1379 //const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
1380 //const DoubleVect& transporteV = transporte;
1381 //const IntTab& face_voisins=domaine_VEF.face_voisins();
1382 //const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
1383 CIntArrView le_bord_num_face = le_bord.num_face().view_ro();
1384 CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
1385 DoubleArrView resu = static_cast<DoubleVect&>(tab_resu).view_rw();
1386 Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), num2, KOKKOS_LAMBDA(const int ind_face)
1387 {
1388 int facei=le_bord_num_face(ind_face);
1389 int ind_face_associee=face_associee(ind_face);
1390 int faceiAss=le_bord_num_face(ind_face_associee);
1391
1392 if (facei<faceiAss)
1393 for (int dim=0; dim<nb_comp; dim++)
1394 {
1395 int ligne=facei*nb_comp+dim;
1396 int ligneAss=faceiAss*nb_comp+dim;
1397 Kokkos::atomic_add(&resu[ligneAss], resu[ligne]);
1398 Kokkos::atomic_store(&resu[ligne], resu[ligneAss]);
1399 }
1400 /*
1401 if (old_centered)
1402 {
1403 //Pour le 1er element voisin
1404 int elem=face_voisins(facei,0);
1405 //facei_loc=-1;
1406 int facei_loc=num_fac_loc(facei,0);
1407 assert(facei_loc!=-1);
1408 double kii=Kij(elem,facei_loc,facei_loc);
1409
1410 //Pour le 2eme element voisin
1411 elem=face_voisins(facei,1);
1412 //facei_loc=-1;
1413 facei_loc=num_fac_loc(facei,1);
1414 assert(facei_loc!=-1);
1415 kii+=Kij(elem,facei_loc,facei_loc);
1416
1417 for (int dim=0; dim<nb_comp; dim++)
1418 {
1419 int ligne=facei*nb_comp+dim;
1420 int ligneAss=faceiAss*nb_comp+dim;
1421
1422 double tmp=kii*transporteV[ligneAss];
1423 resuV[ligneAss]-=tmp;
1424 resuV[ligne]-=tmp;
1425 }
1426 } */
1427
1428 });//fin du for sur "face_i"
1429 end_gpu_timer(__KERNEL_NAME__);
1430 }//fin du if sur "Periodique"
1431
1432 }//fin du for sur "n_bord"
1433}
1434
1435//Fonction qui initialise les attributs "elem_nb_faces_dirichlet_"
1436//et "elem_faces_dirichlet_"
1437//REMARQUE : "elem_nb_faces_dirichlet_" contient le nombre de faces de Dirichlet
1438//pour chaque element du maillage
1439//REMARQUE : "elem_faces_dirichlet_" le numero global des faces de Dirichlet
1440//contenu dans un element quelconque du maillage
1441void Op_Conv_Muscl_New_VEF_Face::calculer_data_pour_dirichlet()
1442{
1443 const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
1444 const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
1445
1446 const IntTab& face_sommets = domaine_VEF.face_sommets();
1447 const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
1448 const IntTab& elem_sommets = domaine_VEF.domaine().les_elems();
1449
1450 const int nb_elem_tot = domaine_VEF.nb_elem_tot();
1451 const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
1452 const int nb_som_faces = domaine_VEF.nb_som_face();
1453 const int nb_som_elem = domaine_VEF.domaine().nb_som_elem();
1454
1455 //On rajoute un flag aux elements qui ont au moins un sommet qui appartient
1456 //une face de Dirichlet ou Dirichlet_homogene
1457 if (old_centered_)
1458 {
1459 is_element_for_upwinding_.resize(nb_elem_tot);
1460 is_element_for_upwinding_=0;
1461
1462 IntTab est_un_sommet_de_bord(domaine_VEF.nb_som_tot());
1463 est_un_sommet_de_bord=0;
1464
1465 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1466 {
1467 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1468 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1469 int nb_faces_tot=le_bord.nb_faces_tot();
1470
1471 if ( (sub_type(Dirichlet,la_cl.valeur()))
1472 || (sub_type(Dirichlet_homogene,la_cl.valeur()))
1473 )
1474 for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
1475 {
1476 int face = le_bord.num_face(ind_face);
1477 for (int som_loc=0; som_loc<nb_som_faces; som_loc++)
1478 {
1479 int som=face_sommets(face,som_loc);
1480 est_un_sommet_de_bord(som)=1;
1481 }
1482 }
1483 }
1484 ToDo_Kokkos("critical");
1485 for (int elem=0; elem<nb_elem_tot; elem++)
1486 {
1487 int rang=rang_elem_non_std(elem);
1488 if (rang!=-1) is_element_for_upwinding_(elem)=1;
1489 else
1490 for (int som_loc=0; som_loc<nb_som_elem; som_loc++)
1491 {
1492 int som=elem_sommets(elem,som_loc);
1493 if (est_un_sommet_de_bord(som))
1494 is_element_for_upwinding_(elem)=1;
1495 }
1496 }
1497 }
1498 is_dirichlet_faces_.resize(domaine_VEF.nb_faces_tot());
1499 is_dirichlet_faces_=0;
1500
1501 for (int n_bord=0; n_bord<nb_bord; n_bord++)
1502 {
1503 const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
1504 const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
1505 int nb_faces_tot=le_bord.nb_faces_tot();
1506
1507 if ( (sub_type(Dirichlet,la_cl.valeur()))
1508 || (sub_type(Dirichlet_homogene,la_cl.valeur()))
1509 )
1510 for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
1511 {
1512 int face = le_bord.num_face(ind_face);
1513 is_dirichlet_faces_(face)=1;
1514 }
1515 }
1516}
1517
1518void Op_Conv_Muscl_New_VEF_Face::completer()
1519{
1520 Op_Conv_VEF_Face::completer();
1521 calculer_data_pour_dirichlet();
1522
1523 //alpha_tab.resize_array(le_dom_vef->nb_faces_tot());
1524 //alpha_tab = alpha_;
1525
1526 //beta_.resize_array(le_dom_vef->nb_faces_tot());
1527 //beta_=1.;
1528}
1529
1530void Op_Conv_Muscl_New_VEF_Face::ajouter_contribution(const DoubleTab& transporte_2, Matrice_Morse& matrice) const
1531{
1532}
1533
1534void Op_Conv_Muscl_New_VEF_Face::modifier_pour_Cl (Matrice_Morse& matrice, DoubleTab& secmem) const
1535{
1536 Op_Conv_VEF_Face::modifier_pour_Cl(matrice,secmem);
1537}
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
Champ_P1NC::calcul_gradient
static DoubleTab & calcul_gradient(const DoubleTab &, DoubleTab &, const Domaine_Cl_VEF &)
Definition Champ_P1NC.cpp:915
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::nb_som_elem
int nb_som_elem() const
Renvoie le nombre de sommets des elements geometriques constituants le domaine.
Definition Domaine.h:474
Domaine_32_64::les_elems
IntTab_t & les_elems()
Definition Domaine.h:129
Domaine_Cl_VEF
Definition Domaine_Cl_VEF.h:40
Domaine_Cl_VEF::type_elem_Cl
int type_elem_Cl(int i) const
Definition Domaine_Cl_VEF.h:62
Domaine_Cl_VEF::vecteur_face_facette_Cl
DoubleTab & vecteur_face_facette_Cl()
Definition Domaine_Cl_VEF.h:60
Domaine_Cl_VEF::normales_facettes_Cl
DoubleTab & normales_facettes_Cl()
Definition Domaine_Cl_VEF.h:56
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_VEF::rang_elem_non_std
IntVect & rang_elem_non_std()
Definition Domaine_VEF.h:86
Domaine_VEF::type_elem
const Elem_VEF_base & type_elem() const
Definition Domaine_VEF.h:75
Domaine_VEF::facette_normales
auto & facette_normales()
Definition Domaine_VEF.h:84
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::xv
double xv(int num_face, int k) const
Definition Domaine_VF.h:76
Domaine_VF::get_num_fac_loc
const IntTab & get_num_fac_loc() const
Definition Domaine_VF.h:140
Domaine_VF::face_sommets
int face_sommets(int i, int j) const
renvoie le numero du ieme sommet de la face num_face.
Definition Domaine_VF.h:583
Domaine_VF::nb_som_face
int nb_som_face() const
renvoie le nombre de sommets par face.
Definition Domaine_VF.h:494
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::premiere_face_int
int premiere_face_int() const
une face est interne ssi elle separe deux elements.
Definition Domaine_VF.h:463
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::nb_elem_tot
int nb_elem_tot() const
Definition Domaine_dis_base.h:50
Domaine_dis_base::domaine
const Domaine & domaine() const
Definition Domaine_dis_base.h:46
Domaine_dis_base::nb_som_tot
int nb_som_tot() const
Definition Domaine_dis_base.h:52
Echange_impose_base
classe Echange_impose_base: Cette condition limite sert uniquement pour l'equation d'energie.
Definition Echange_impose_base.h:41
Elem_VEF_base::KEL
const IntTab & KEL() const
Definition Elem_VEF_base.h:31
Elem_VEF_base::nb_facette
virtual int nb_facette() const =0
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
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::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
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
classe Neumann_sortie_libre Cette classe represente une frontiere ouverte sans vitesse imposee
Definition Neumann_sortie_libre.h:34
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
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::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_Muscl_New_VEF_Face
class Op_Conv_Muscl_New_VEF_Face
Definition Op_Conv_Muscl_New_VEF_Face.h:43
Op_Conv_Muscl_New_VEF_Face::calculer_flux_bords
public_for_cuda void calculer_flux_bords(const DoubleTab &, const DoubleTab &, const DoubleTab &) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:765
Op_Conv_Muscl_New_VEF_Face::ajouter_antidiffusion_v1
DoubleTab & ajouter_antidiffusion_v1(const DoubleTab &, const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:1184
Op_Conv_Muscl_New_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_Muscl_New_VEF_Face.cpp:1518
Op_Conv_Muscl_New_VEF_Face::modifier_pour_Cl
void modifier_pour_Cl(Matrice_Morse &, DoubleTab &) const override
DOES NOTHING - to override in derived classes.
Definition Op_Conv_Muscl_New_VEF_Face.cpp:1534
Op_Conv_Muscl_New_VEF_Face::remplir_fluent
void remplir_fluent() const override
Definition Op_Conv_Muscl_New_VEF_Face.cpp:569
Op_Conv_Muscl_New_VEF_Face::ajouter_antidiffusion_v2
DoubleTab & ajouter_antidiffusion_v2(const DoubleTab &, const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:1082
Op_Conv_Muscl_New_VEF_Face::calculer_flux_operateur_centre
void calculer_flux_operateur_centre(DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const DoubleTab &, const int, const DoubleTab &, const DoubleTab &) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:268
Op_Conv_Muscl_New_VEF_Face::calculer_coefficients_operateur_centre
void calculer_coefficients_operateur_centre(DoubleTab &, DoubleTab &, DoubleTab &, DoubleTab &, const int, const DoubleTab &vitesse) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:98
Op_Conv_Muscl_New_VEF_Face::ajouter_contribution
void ajouter_contribution(const DoubleTab &, Matrice_Morse &) const override
Definition Op_Conv_Muscl_New_VEF_Face.cpp:1530
Op_Conv_Muscl_New_VEF_Face::calculer_dt_stab
double calculer_dt_stab() const override
Calcul dt_stab.
Definition Op_Conv_Muscl_New_VEF_Face.cpp:659
Op_Conv_Muscl_New_VEF_Face::mettre_a_jour_pour_periodicite
void mettre_a_jour_pour_periodicite(const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:1360
Op_Conv_Muscl_New_VEF_Face::ajouter
DoubleTab & ajouter(const DoubleTab &, DoubleTab &) const override
Definition Op_Conv_Muscl_New_VEF_Face.cpp:593
Op_Conv_Muscl_New_VEF_Face::ajouter_operateur_centre
DoubleTab & ajouter_operateur_centre(const DoubleTab &, const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:864
Op_Conv_Muscl_New_VEF_Face::ajouter_diffusion
DoubleTab & ajouter_diffusion(const DoubleTab &, const DoubleTab &, const DoubleTab &, DoubleTab &) const
Definition Op_Conv_Muscl_New_VEF_Face.cpp:999
Op_Conv_VEF_Face
class Op_Conv_VEF_Face
Definition Op_Conv_VEF_Face.h:35
Op_Conv_VEF_Face::gradient_elem_
DoubleTab gradient_elem_
Definition Op_Conv_VEF_Face.h:74
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::alpha_
double alpha_
Definition Op_Conv_VEF_Face.h:65
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
class Op_Conv_VEF_base
Definition Op_Conv_VEF_base.h:37
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_Conv_base::fixer_dt_stab_conv
void fixer_dt_stab_conv(double dt)
Definition Operateur_Conv_base.h:63
Operateur_base::flux_bords_
DoubleTab flux_bords_
Definition Operateur_base.h:140
Operateur_base::calculer_dt_stab
virtual double calculer_dt_stab() const
Calcul dt_stab.
Definition Operateur_base.cpp:129
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
Process::mp_max
static double mp_max(double)
Definition Process.cpp:376
Process::exit
static void exit(int exit_code=-1)
Routine de sortie de TRUST dans une region Kokkos.
Definition Process.cpp:455
Sortie
Classe de base des flux de sortie.
Definition Sortie.h:52
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::view_wo
std::enable_if_t< is_default_exec_space< EXEC_SPACE >, View< _TYPE_, _SHAPE_ > > view_wo()
Definition TRUSTTab.h:276
TRUSTTab::nb_dim
int nb_dim() const
Definition TRUSTTab.h:199
TRUSTTab::resize
void resize(_SIZE_ n, RESIZE_OPTIONS opt=RESIZE_OPTIONS::COPY_INIT)
Definition TRUSTTab.tpp:469
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::mp_max_vect
_TYPE_ mp_max_vect(Mp_vect_options opt=VECT_REAL_ITEMS) const
Definition TRUSTVect.h:158