1 : C SUPGEF SOURCE CHAT 05/01/13 03:27:12 5004 2 : SUBROUTINE SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 3 : & NES,NP,NPG,IAXI,XCOOR, 4 : & LE,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP, 5 : & TN,IPADT,UN,IPADU,NPTD,ANUK, 6 : & WT,WS,HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, 7 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 8 : C-------------------------------------------------------------------- 9 : C Calcul des fonctions tests associées à la formulation Eléments 10 : C Finis Petrov-Galerkin afin de stabiliser les termes de convection. 11 : C On en profite pour calculer des temps caractéristiques. 12 : C-------------------------------------------------------------------- 13 : C Cette subroutine étant écrite en fortran pur, certains arguments 14 : C servent uniquement au dimensionnement des tableaux. 15 : C Afin de doper les calculs, on decoupe la boucle sur les éléments 16 : C par paquets afin de bénéficier à fond de la vectorisation. 17 : C-------------------------------------------------------------------- 18 : C 19 : C-------------------------- 20 : C Paramètre Entrée/Sortie : 21 : C-------------------------- 22 : C 23 : C /S FN : Fonction de forme associé à la transformation géométrique 24 : C /S GR : Gradient de FN dans l'élément de référence 25 : C /S PG : Poids d'intégration associé à chaque point de Gauss 26 : C /S XYZ : Coordonnées des sommets de l'élément 27 : C /S AJ : Jacobien 28 : C /S HR : Gradient des fonctions de forme dans l'élément courant 29 : C /S PGSQ : Produit du poids d'intégration par detJ 30 : C /S RPG : Rayon associé à chaque point de Gauss (cas axisymétrique) 31 : C 32 : C E/ NES : Dimension en espace du problème traité (2 en 2D et 3 en 3D) 33 : C E/ NP : Nombre de points support de ddl pour l'élément considéré 34 : C E/ NPG : Nombre de points de Gauss 35 : C E/ IAXI : Précise l'axe de symétrie dans le cas axisymétrique 36 : C (en fait l'axe de symétrie est toujours y -> IAXI=2) 37 : C E/ XCOOR : Coordonnées de l'ensemble des points (voir SMCOORD.INC) 38 : C Petrov-Galerkin choisi pour stabiliser la convection 39 : C Petrov-Galerkin choisi pour stabiliser la convection 40 : C /S HK : Gradient de FN dans l'élément courant 41 : C /S PGSK : Produit du poids d'intégration par detJ 42 : C /S RPGK : Rayon associé à chaque point de Gauss (cas axisymétrique) 43 : C /S AIRE : Volume de chaque élément 44 : C /S AJK : Jacobien 45 : C /S UIL : Vecteur vitesse aux points de Gauss 46 : C /S DUIL : Gradient du champ de vitesse 47 : C (la divergence est obtenue par sommation) 48 : C E/ KDEB : Indice du début de la fenetre 49 : C E/ KFIN : Indice de fin de la fenetre 50 : C E/ LRV : Dimension de la fenetre 51 : C E/ TN : Inconnue transportée au pas de temps précédent 52 : C E/ IPADT : Correspondance numérotation locale et globale pour les 53 : C points SOMMET : J=LECT(I) : le point numéro I est rangé 54 : C en Jème position pour le CHAMP TN 55 : C E/ UN : Champ de vitesse transportant aux sommets des éléments 56 : C E/ IPADU : Correspondance numérotation locale et globale pour les 57 : C points SOMMET : J=LECT(I) : le point numéro I est rangé 58 : C en Jème position pour le CHAMP UN 59 : C E/ NPTD : Nombre de points sous-tendant le spg associé à UN 60 : C E/ ANUK : Champ contenant le coefficient de diffusion (diffusivité 61 : C pour l'énergie, rapport viscosité sur densité pour qdm) 62 : C E/ IDCEN : Contient le schéma associé au terme convectif 63 : C (1=CENTREE 2=SUPGDC 3=SUPG 4=TVISQUEUX 5=CNG) 64 : C E/ DTM1 : Pas de temps imposé 65 : C E/S DT : Pas de temps 66 : C /S DTT1 : Temps caractéristique associé à la diffusion/convection 67 : C /S DTT2 : Temps caractéristique associé à la diffusion 68 : C /S DIAEL : Longueur carac. de l'élément le plus pénalisant pour dt 69 : C /S NUEL : Numéro de l'élément le plus pénalisant pour dt 70 : C 71 : C Variable locale 72 : C 73 : C DDUIL : gradient de chaque composante de la vitesse au pt Gauss 74 : C pour un élément. 75 : C 76 : C RTNL : valeur de TN au point de gauss 77 : C 78 : C-------------------------------------------------------------------- 79 : IMPLICIT INTEGER(I-N) 80 : IMPLICIT REAL*8 (A-H,O-Z) 81 : -INC CCOPTIO 82 : C?-INC CCREEL 83 : DIMENSION IPADT(*),IPADU(*),LE(NP,*) 84 : 85 : DIMENSION XCOOR(*) 86 : DIMENSION UN(NPTD,IDIM),TN(*) 87 : DIMENSION ANUK(LRV) 88 : 89 : DIMENSION FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG),XYZ(IDIM,NP) 90 : DIMENSION HR(IDIM,NP,NPG),PGSQ(NPG),RPG(NPG),AJ(IDIM,IDIM,NPG) 91 : DIMENSION WT(LRV,NP,NPG),WS(LRV,NP,NPG),HK(LRV,IDIM,NP,NPG) 92 : DIMENSION PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV) 93 : DIMENSION AJK(LRV,IDIM,IDIM,NPG) 94 : C 95 : DIMENSION UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG) 96 : DIMENSION UPIL(3),GRAD(3) 97 : DATA ACCT,ACC2/1.,1./ , ipas/0/ 98 : C 99 : REAL*8 DDUIL(3,3),RTNL 100 : C 101 : C 102 : C- Pour les éléments de la fenetre KDEB:KFIN initialisations des 103 : C- données associées à l'élément fini utilisé : fonctions de forme, 104 : C- gradients des fonctions de forme, aire de l'élément, poids des 105 : C- points de Gauss et rayon des points de Gauss pour le cas axi. 106 : C 107 : C? write(6,*)'UN(NPTD,IDIM) ',nptd,idim 108 : C? write(6,1002)un 109 : XPETI=1.D-30 109 : CALL INITD(UPIL,3,0.D0) 110 : 111 : DO 5020 K=KDEB,KFIN 112 : KP = K-KDEB+1 113 : 114 : DO 20 I=1,NP 115 : J = LE(I,K) 116 : DO 10 N=1,IDIM 117 : XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N) 118 : 10 CONTINUE 119 : 20 CONTINUE 120 : 120 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, 121 : & NES,IDIM,NP,NPG,IAXI,AIRE(KP),AJ,SGN) 122 : 123 : DO 8934 N=1,IDIM 124 : DO 8934 M=1,IDIM 125 : DO 8934 L=1,NPG 126 : AJK(KP,M,N,L)=AJ(M,N,L) 127 : 8934 CONTINUE 128 : 128 : CALL CALJTR(GR,XYZ,NES,IDIM,NP,NPG,AJ) 129 : 130 : C 131 : DO 5030 L=1,NPG 132 : 133 : DO 40 N=1,IDIM 134 : DO 30 I=1,NP 135 : HK(KP,N,I,L) = HR(N,I,L) 136 : 30 CONTINUE 137 : 40 CONTINUE 138 : 139 : PGSK(KP,L) = PGSQ(L) 140 : RPGK(KP,L) = RPG(L) 141 : C 142 : C- Calcul en chaque élément, pour chaque point de Gauss 143 : C- UIL : Vitesse 144 : C 145 : DO 230 N=1,IDIM 146 : DUIL(KP,N,L) = XPETI 147 : UIL(KP,N,L) = XPETI 148 : DO 210 I=1,NP 149 : NF = IPADU(LE(I,K)) 150 : UIL(KP,N,L) = UIL(KP,N,L) + UN(NF,N)*FN(I,L) 151 : DUIL(KP,N,L) = DUIL(KP,N,L) + UN(NF,N)*HK(KP,N,I,L) 152 : 210 CONTINUE 153 : 230 CONTINUE 154 : C 155 : C Dans le cas d'une formulation conservative du transport, 156 : C La vitesse de transport réelle est d(uT)/d(T). Cela revient 157 : C à rajouter le terme T (dU/dT) à la vitesse U. 158 : C Le terme dU/dT peut s'écrite pour chaque composante : 159 : C dx/dT dU/dx + dy/dT dU/dy 160 : C Lorsque dx/dT ou dy/dT est infini, cela signifie que les variables 161 : C sont indépendantes de T. Le terme correspondant doit alors être 162 : C éliminé du caclul. 163 : C 164 : IF (IDCEN .EQ. 2 .OR. IDCEN .EQ. 3) THEN 165 : IF (IKOMP .EQ. 1 .OR. IKOMP .EQ. 2) THEN 166 : C 167 : C Evaluation du gradient de TN. 168 : C 169 : DO 171 N=1,IDIM 170 : GRAD(N)=0.D0 171 : DO 172 I=1,NP 172 : NT = IPADT(LE(I,K)) 173 : GRAD(N) = GRAD(N) + TN(NT) *HK(KP,N,I,L) 174 : 172 CONTINUE 175 : IF (ABS(GRAD(N)) .LT. 1.D-10) THEN 176 : C La variable x ou y est indépendante de T. On met le 177 : C gradient à l'infini pour le déconnecter des calculs 178 : C suivants. 179 : GRAD(N) = 1.D+100 180 : ENDIF 181 : 171 CONTINUE 182 : * WRITE(*,*) 'GRADT', GRAD(1),GRAD(2) 183 : C 184 : C On évalue le gradient de chaque composante de la vitesse. 185 : C 186 : DO 231 N=1,IDIM 187 : DO 231 M=1,IDIM 188 : DDUIL(N,M) = XPETI 189 : DO 211 I=1,NP 190 : NF = IPADU(LE(I,K)) 191 : DDUIL(N,M) = DDUIL(N,M) + UN(NF,N)*HK(KP,M,I,L) 192 : 211 CONTINUE 193 : 231 CONTINUE 194 : C 195 : C Correction en axisymetrique 196 : C 197 : C? IF (IAXI .EQ. 1) THEN 198 : C? DO 212 I=1,NP 199 : C? NF = IPADU(LE(I,K)) 200 : C? DDUIL(1,1)=DDUIL(1,1)+(UN(NF,1)*FN(I,L)/RPGK(KP,L)) 201 : C? 212 CONTINUE 202 : C? ENDIF 203 : C 204 : C On évalue la variable scalaire TN au point de gauss L 205 : C 206 : RTNL=0.D0 207 : DO 173 I=1,NP 208 : NT = IPADT(LE(I,K)) 209 : RTNL = RTNL + TN(NT) * FN(I,L) 210 : 173 CONTINUE 211 : * WRITE(*,*) 'RTNL', RTNL 212 : C 213 : C Modification de la vitesse U 214 : C 215 : IF (IKOMP .EQ. 1) THEN 216 : c decentrement de UN*TN en UN + TN*(dUN/dTN). 217 : DO 174 N=1,IDIM 218 : DO 174 M=1,IDIM 219 : UIL(KP,N,L) = UIL(KP,N,L)+RTNL*(DDUIL(N,M)/GRAD(M)) 220 : 174 CONTINUE 221 : ELSE 222 : c décentrement de JN=UN*TN en dJN/dTN 223 : DO 175 N=1,IDIM 224 : DO 175 M=1,IDIM 225 : UIL(KP,N,L) = DDUIL(N,M)/GRAD(M) 226 : 175 CONTINUE 227 : ENDIF 228 : C 229 : ENDIF 230 : ENDIF 231 : C 232 : C- Calcul pour chaque point de Gauss de chaque élément de : 233 : C- /L UML : Norme de la vitesse aux points de Gauss 234 : C- BM : module d'un temps caractéristique associé à la convection 235 : C- XMB : caractéristique géométrique 1/2 He 236 : C 237 : 238 : C 239 : BMI=0.D0 240 : UL=0.D0 241 : DO 310 N=1,IDIM 242 : UHAT=0.D0 243 : DO 311 M=1,IDIM 244 : UHAT=UHAT+AJ(M,N,L)*UIL(KP,M,L) 245 : 311 CONTINUE 246 : 247 : BMI=BMI+UHAT*UHAT 248 : UL=UL+UIL(KP,N,L)*UIL(KP,N,L) 249 : 310 CONTINUE 250 : BM=SQRT(BMI) + XPETI 251 : UML=SQRT(UL) 252 : XMB=UML/BM 253 : 254 : 255 : C 256 : C 257 : C- Calcul en chaque élément, pour chaque point de Gauss de 258 : C- GRAD : vecteur unitaire porté par le gradient du champ scalaire 259 : C- UP : projection de la vitesse sur la direction donnée par GRAD 260 : C- UPIL : vecteur UP*GRAD aux points de Gauss 261 : C- pour l'option SUPGDC 262 : C 263 : IF (IDCEN.EQ.2) THEN 264 : DO 170 N=1,IDIM 265 : GRAD(N)=0.D0 266 : DO 170 I=1,NP 267 : NT = IPADT(LE(I,K)) 268 : GRAD(N) = GRAD(N) + TN(NT) *HK(KP,N,I,L) 269 : 170 CONTINUE 270 : 271 : AX=0.D0 272 : DO 2301 M=1,IDIM 273 : AX = AX + GRAD(M)*GRAD(M) 274 : 2301 CONTINUE 275 : AX = SQRT(AX) + XPETI 276 : 277 : UPL=0.D0 278 : DO 2302 N=1,IDIM 279 : GRAD(N) = GRAD(N) / AX 280 : UPL = UPL + GRAD(N) * UIL(KP,N,L) 281 : 2302 CONTINUE 282 : 283 : DO 2303 N=1,IDIM 284 : UPIL(N) = GRAD(N) * UPL 285 : 2303 CONTINUE 286 : 287 : C 288 : BPI=0.D0 289 : DO 410 N=1,IDIM 290 : UPHAT=0.D0 291 : DO 411 M=1,IDIM 292 : UPHAT=UPHAT+AJ(M,N,L)*UPIL(M) 293 : 411 CONTINUE 294 : 295 : BPI=BPI+UPHAT*UPHAT 296 : 410 CONTINUE 297 : BP=SQRT(BPI) + XPETI 298 : XPB=UPL/BP 299 : 300 : ENDIF 301 : C 302 : C----------------------------- 303 : C- DECENTREMENT suivant IDCEN 304 : C----------------------------- 305 : C On calcule dans chaque cas TO1 et TO2 ainsi que le tenseur 306 : C associé à la viscosité numérique afin d'évaluer le pas de 307 : C temps de stabilité pour les schémas explicites. 308 : C 309 : C---------- 310 : C CENTREE : 311 : C---------- 312 : IF (IDCEN.EQ.1) THEN 313 : SI1 = 0.D0 314 : SI2 = 0.D0 315 : TO1 = 0.D0 316 : TO2 = 0.D0 317 : C--------- 318 : C SUPGDC : Base théorique dans : A New FE formulation for computational 319 : C fluid dynamics, II Beyond SUPG, HUGHES et al., in Comp.Meth.Appl.Mech. 320 : C Eng., vol 54, pp 341-355 (1986). 321 : C--------- 322 : ELSEIF (IDCEN.EQ.2) THEN 323 : SI1 = 1.D0 324 : SI2 = 1.D0 325 : C 326 : C- Approximation "doublement asymptotique" basée sur la vitesse moyenne 327 : C- HMK : Distance basé sur la vitesse moyenne 328 : C- ALFA : Peclet de maille basé sur la vitesse moyenne 329 : C 330 : ALFA = UML * XMB / (ANUK(KP)+XPETI) / 3.D0 331 : AKSI = MIN(1.D0,ALFA) 332 : CCT = AKSI / BM * ACCT * CMD 333 : C 334 : C- Approximation "doublement asymptotique" basée sur la projection de 335 : C- la vitesse sur le gradient du champ scalaire 336 : C- HMK : Distance basée sur la vitesse projetée 337 : C- ALFA : Peclet de maille basé sur la vitesse projetée 338 : C 339 : ALFA = UPL * XPB / (ANUK(KP)+XPETI) / 3.D0 340 : AKSI = MIN(1.D0,ALFA) 341 : CCP = AKSI / BP 342 : C 343 : CPT = CCP - CCT 344 : CC2 = MAX(0.D0,CPT) * ACC2 * CMD 345 : C 346 : TO1 = CCT 347 : TO2 = CC2 348 : C------- 349 : C SUPG : 350 : C------- 351 : ELSEIF (IDCEN.EQ.3) THEN 352 : SI1 = 1.D0 353 : SI2 = 1.D0 354 : C 355 : C- Approximation "doublement asymptotique" basée sur la vitesse moyenne 356 : C- HMK : Distance basé sur la vitesse moyenne 357 : C- ALFA : Peclet de maille basé sur la vitesse moyenne 358 : C 359 : ALFA = UML * XMB / (ANUK(KP)+XPETI) / 3.D0 360 : AKSI = MIN(1.D0,ALFA) 361 : CCT = AKSI / BM * ACCT * CMD 362 : C 363 : TO1 = CCT 364 : TO2 = 0.D0 365 : C------------------- 366 : C Tenseur Visqueux : 367 : C------------------- 368 : ELSEIF (IDCEN.EQ.4) THEN 369 : SI1 =-1.D0 370 : SI2 = 1.D0 371 : DT19 = DTM1 * 0.5D0 372 : TO1 = DT19 373 : TO2 = 0.D0 374 : C----------------------------- 375 : C Crank Nicholson généralisé : 376 : C----------------------------- 377 : ELSEIF (IDCEN.EQ.5) THEN 378 : SI1 =-1.D0 379 : SI2 = 1.D0 380 : DT19 = DTM1/6.D0 381 : TO1 = DT19 382 : TO2 = 0.D0 383 : ENDIF 384 : C 385 : C--------------------------- 386 : C Pas de temps de stabilité 387 : C--------------------------- 388 : C 389 : C La viscosité utilisée dans l'évaluation des pas de temps de stabilité 390 : C explicite ajoute à la viscosité physique la viscosité numérique : 391 : C DT1 : Correspond à une CFL de 1 et un Peclet de maille de 1 392 : C CFL=1 -> dx=Vdt | 393 : C |-> dt=2D/V2 394 : C Pe=1 -> dx=2D/V | 395 : C DT2 : Correspond au pas de temps de stabilité lié à la diffusion 396 : C FOU=1 -> dt=0.5 dx2/D 397 : C 398 : IF (IDCEN.NE.5) THEN 399 : DT0 = DT 400 : C DT1 = 2.D0 / 401 : C & ( UIL(KP,1,L)*UIL(KP,1,L)/(ANUK(KP)+XPETI) 402 : C & + UIL(KP,2,L)*UIL(KP,2,L)/(ANUK(KP)+XPETI) ) 403 : C DT2 = 0.5D0 / 404 : C & ( ANUK(KP)*AL2(KP) 405 : C & + ANUK(KP)*AH2(KP) ) 406 : DT1 = 2.D0 * ANUK(KP) / (UML * UML) 407 : DT2 = DT1 408 : C 409 : IF (DT1.LT.DT) DT=DT1 410 : IF (DT2.LT.DT) DT=DT2 411 : IF (DT.NE.DT0) THEN 412 : DTT1 = DT1 413 : DTT2 = DT2 414 : DIAEL = XMB * 2.D0 415 : NUEL = K 416 : ENDIF 417 : ENDIF 418 : C 419 : C--------------------------------------------------- 420 : C Fonction test pour la formulation Petrov-Galerkin 421 : C--------------------------------------------------- 422 : C Ce qui est diffusion numérique en explicite se transforme en 423 : C ajoutant de la viscosité numérique (Tenseurs visqueux et CNG). 424 : C WS : Fonction test pour la partie explicite 425 : C 426 : IF(IDIM.EQ.2)THEN 427 : DO 2050 I=1,NP 428 : WT(KP,I,L) = FN(I,L) + SI1 * 429 : & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L)) 430 : & +TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L))) 431 : WS(KP,I,L) = FN(I,L) + SI2 * 432 : & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L)) 433 : & +TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L))) 434 : 2050 CONTINUE 435 : C 436 : ELSEIF(IDIM.EQ.3)THEN 437 : DO 2051 I=1,NP 438 : WT(KP,I,L) = FN(I,L) + SI1 * 439 : & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L)+ 440 : & UIL(KP,3,L)*HK(KP,3,I,L) ) 441 : & +TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L)+ 442 : & UPIL(3)*HK(KP,3,I,L) )) 443 : 444 : WS(KP,I,L) = FN(I,L) + SI2 * 445 : & (TO1*(UIL(KP,1,L)*HK(KP,1,I,L)+UIL(KP,2,L)*HK(KP,2,I,L)+ 446 : & UIL(KP,3,L)*HK(KP,3,I,L)) 447 : & + TO2*(UPIL(1)*HK(KP,1,I,L)+UPIL(2)*HK(KP,2,I,L)+ 448 : & UPIL(3)*HK(KP,3,I,L) )) 449 : 2051 CONTINUE 450 : ENDIF 451 : C 452 : C- Si on est en conservatif, on rétablit les valeurs de la vitesse 453 : C de transport qui ont été modifiées car elles sont utilisées 454 : C dans d'autres subroutines ultérieures. 455 : C 456 : IF (IKOMP .EQ. 1) THEN 457 : DO 235 N=1,IDIM 458 : UIL(KP,N,L) = XPETI 459 : DO 215 I=1,NP 460 : NF = IPADU(LE(I,K)) 461 : UIL(KP,N,L) = UIL(KP,N,L) + UN(NF,N)*FN(I,L) 462 : 215 CONTINUE 463 : 235 CONTINUE 464 : ENDIF 465 : C 466 : 5030 CONTINUE 467 : 468 : 5020 CONTINUE 469 : 470 : RETURN 471 : 1001 FORMAT(20(1X,I5)) 472 : 1002 FORMAT(10(1X,1PE11.4)) 473 : END 474 : 475 : 476 : 477 :
© Cast3M 2003 - All rights reserved.
Disclaimer