1 : C ZCTSCL SOURCE CHAT 05/01/13 04:21:41 5004 2 : SUBROUTINE ZCTSCL 3 : C 4 : C VERSION VECTORISEE 5 : C 6 : C Les éléments sont groupés en paquets de LRV éléments, LRV étant 7 : C la longueur des registres vectoriels de la machine cible, i.e 8 : C 64 sur Cray, 128 ou 256 sur IBM 3090VF. On promène une fenêtre 9 : C de longueur LRV sur la boucle générale de longueur NEL. 10 : C 11 : & (HR,RPG,DRR,LE,NEL,K0,IES,NP,IAXI, 12 : & IPADI,IPADU,IPADF,IKOMP,IKAS, 13 : & ALFE,IND1,UN,INDU,NPTU,TN,NPTD,QE,IKS, 14 : & HRN,G,NPTI,ALT,SGT, 15 : & VOLU,COTE,NELZ,IDCEN,IPG, 16 : & DTM1,DT,DTT1,DTT2,NUEL,DIAEL,FN) 17 : 18 : IMPLICIT INTEGER(I-N) 19 : IMPLICIT REAL*8 (A-H,O-Z) 20 : 21 : C*********************************************************************** 22 : C 23 : C CE SP DISCRETISE UNE EQUATION GENERALE DE TRANSPORT-DIFFUSION AVEC 24 : C SOURCE. 25 : C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI 26 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 27 : C LES OPERATEURS SONT "SOUS-INTEGRES" 28 : C 29 : C 30 : C APPELE PAR YTSCAL 31 : C 32 : C 33 : C*********************************************************************** 34 : 35 : -INC CCVQUA4 36 : -INC CCREEL 37 : C 38 : C Longueur des registres vectoriels de la machine cible 39 : C On prend 64 pour ne pas augmenter la taille des tableaux 40 : C nécessaires à la vectorisation. 41 : C 42 : PARAMETER(LRV=64) 43 : 44 : DIMENSION UN(NPTU,IES),HRN(NPTI),TN(NPTD) 45 : DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*) 46 : DIMENSION ALFE(*),ALT(*) 47 : 48 : DIMENSION IPADI(*),IPADU(*),IPADF(*),LE(NP,1) 49 : DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL) 50 : 51 : DIMENSION BF(9,9) 52 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) 53 : 54 : DIMENSION AIRE(LRV) 55 : DIMENSION AL(LRV),AH(LRV),AP(LRV) 56 : DIMENSION ALFT(LRV),QT(LRV) 57 : DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9) 58 : DIMENSION TETAC(LRV,9),TETAD(LRV,9),TETA(LRV,9) 59 : DIMENSION UMI(LRV,3) 60 : DIMENSION SBF(LRV,9) 61 : DIMENSION WT(LRV,9),CHGLD(LRV),CHGLP(LRV) 62 : 63 : REAL*8 G(NPTI),FN(NP,*) 64 : 65 : SAVE IPAS,QGGT,Q1,Q2,Q3 66 : 67 : DATA CD/1.D0/ 68 : 69 : DATA IPAS/0/ 70 : C************************************************************************ 71 : C 72 : C INITIALISATIONS DIVERSES 73 : C 74 : ZERMA=XPETIT 75 : 76 : NK=K0 77 : C ******** 78 : C * 2D * 79 : C ******** 80 : 81 : 82 : IF(IES.EQ.3)GO TO 10 83 : 84 : IAX1=0 85 : IAX2=0 86 : IF(IAXI.EQ.1)IAX2=1 87 : IF(IAXI.EQ.2)IAX1=1 88 : 89 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 90 : 91 : C DIFFERENCES TRIANGLE / QUADRANGLE 92 : QUA4=0.D0 93 : IF(NP.EQ.4)QUA4=1.D0 94 : C 95 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 96 : 97 : C 98 : C Calcul du nombre de paquets de LRV éléments 99 : C 100 : NNN=MOD(NEL,LRV) 101 : IF(NNN.EQ.0) NPACK=NEL/LRV 102 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 103 : KPACKD=1 104 : KPACKF=NPACK 105 : C 106 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 107 : C 108 : DO 7001 KPACK=KPACKD,KPACKF 109 : C 110 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 111 : C 112 : C 1. Calcul des limites du paquet courant. 113 : KDEB=1+(KPACK-1)*LRV 114 : KFIN=MIN(NEL,KDEB+LRV-1) 115 : C 116 : DO 7002 K=KDEB,KFIN 117 : KP=K-KDEB+1 118 : NK=K+K0 119 : NK1=(1-IND1)*(NK-1)+1 120 : ALFT(KP)=ALFE(NK1)+ZERMA 121 : AIRE(KP)=VOLU(NK) 122 : AL(KP)=COTE(NK,1) 123 : AH(KP)=COTE(NK,2) 124 : 7002 CONTINUE 125 : 126 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 127 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 128 : DO 7003 K=KDEB,KFIN 129 : KP=K-KDEB+1 130 : NK=K+K0 131 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 132 : 7003 CONTINUE 133 : ENDIF 134 : 135 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 136 : 137 : DO 7006 I=1,NP 138 : DO 7006 K=KDEB,KFIN 139 : KP=K-KDEB+1 140 : NF=IPADF(LE(I,K)) 141 : NU=IPADU(LE(I,K)) 142 : NI=IPADI(LE(I,K)) 143 : NFU=(1-INDU)*(NU-1)+1 144 : UIX(KP,I)=UN(NFU,1) 145 : UIY(KP,I)=UN(NFU,2) 146 : TETAC(KP,I)=HRN(NI) 147 : TETAD(KP,I)=TN(NF) 148 : 7006 CONTINUE 149 : 149 : CALL KSUPG1(WT,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV, 150 : &HRN,IPADI,UN,ALFT,NPTU,NEL,NP,DRR,HR,FN, 151 : &AIRE,AL,AH,AP,IDCEN,IPADU,LE,QUA4,IKOMP, 152 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 153 : C 154 : C Initialisation de la variable d'accumulation SBF au terme source 155 : C 156 : 157 : IF(IKOMP.EQ.0)THEN 158 : 159 : DO 70021 K=KDEB,KFIN 160 : KP=K-KDEB+1 161 : NK=K+K0 162 : NKS=(1-IKS)*(NK-1)+1 163 : QT(KP)=QE(NKS) 164 : 70021 CONTINUE 165 : 166 : IF(IPG.EQ.0)THEN 167 : DO 70062 I=1,NP 168 : DO 70062 K=KDEB,KFIN 169 : KP=K-KDEB+1 170 : SBF(KP,I)=-QT(KP)*DRR(I,K) 171 : 70062 CONTINUE 172 : ELSE 173 : DO 71062 I=1,NP 174 : DO 71062 K=KDEB,KFIN 175 : KP=K-KDEB+1 176 : SBF(KP,I)=-QT(KP)*WT(KP,I) 177 : 71062 CONTINUE 178 : ENDIF 179 : 180 : ELSEIF(IKOMP.EQ.1)THEN 181 : 182 : DO 70023 K=KDEB,KFIN 183 : KP=K-KDEB+1 184 : NK=K+K0 185 : NKS=(1-IKS)*(NK-1)+1 186 : QT(KP)=QE(NKS) 187 : 70023 CONTINUE 188 : 189 : IF(IPG.EQ.0)THEN 190 : DO 70064 I=1,NP 191 : DO 70064 K=KDEB,KFIN 192 : KP=K-KDEB+1 193 : SBF(KP,I)=-QT(KP)*DRR(I,K) 194 : 70064 CONTINUE 195 : ELSE 196 : DO 71064 I=1,NP 197 : DO 71064 K=KDEB,KFIN 198 : KP=K-KDEB+1 199 : SBF(KP,I)=-QT(KP)*WT(KP,I) 200 : 71064 CONTINUE 201 : ENDIF 202 : 203 : ENDIF 204 : 205 : C Le coeur du calcul ... 206 : 207 : IF(IKOMP.EQ.0)THEN 208 : 209 : DO 7014 I=1,NP 210 : DO 7014 J= 1,NP 211 : DO 7014 K=KDEB,KFIN 212 : KP=K-KDEB+1 213 : ZVGG=AIRE(KP)*CHGLP(KP)*VGGT(J,I) 214 : 215 : ZVGT=AIRE(KP)*( 216 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 217 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 218 : &+ CHGLD(KP)*VGGT(J,I) ) 219 : 220 : V2=(UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2))*WT(KP,I) 221 : 222 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 223 : 224 : 7014 CONTINUE 225 : 226 : ELSEIF(IKOMP.EQ.1)THEN 227 : 228 : DO 7015 I=1,NP 229 : DO 7015 J= 1,NP 230 : DO 7015 K=KDEB,KFIN 231 : KP=K-KDEB+1 232 : ZVGG=AIRE(KP)*CHGLP(KP)*VGGT(J,I) 233 : 234 : ZVGT=AIRE(KP)*( 235 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 236 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 237 : &+ CHGLD(KP)*VGGT(J,I) ) 238 : 239 : V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2))*WT(KP,I) 240 : 241 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 242 : 243 : 7015 CONTINUE 244 : 245 : ENDIF 246 : C 247 : C Fin de l'accumulation dans SBF. 248 : C On ajoute ces incréments G. 249 : C 250 : DO 7017 I=1,NP 251 : DO 7017 K=KDEB,KFIN 252 : KP=K-KDEB+1 253 : NF=IPADI(LE(I,K)) 254 : G(NF) = G(NF)+SBF(KP,I) 255 : 7017 CONTINUE 256 : 257 : 7001 CONTINUE 258 : 259 : C WRITE(6,*)' G DANS YCTSCL ' 260 : C WRITE(6,1984)(M,G(M),M=1,NPTD) 261 : 1984 FORMAT(7(1X,I4,2X,1PE11.4)) 262 : 262 : C CALL ARRET(0) 263 : IPAS=1 264 : RETURN 265 : 266 : C ******** 267 : C * 3D * 268 : C ******** 269 : 270 : 10 CONTINUE 271 : 271 : IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES) 272 : CUB8=0.D0 273 : IF(NP.EQ.8)CUB8=1.D0 274 : 275 : C 276 : C Calcul du nombre de paquets de LRV éléments 277 : C 278 : NNN=MOD(NEL,LRV) 279 : IF(NNN.EQ.0) NPACK=NEL/LRV 280 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 281 : KPACKD=1 282 : KPACKF=NPACK 283 : C 284 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 285 : C 286 : DO 8001 KPACK=KPACKD,KPACKF 287 : C 288 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 289 : C 290 : C 1. Calcul des limites du paquet courant. 291 : KDEB=1+(KPACK-1)*LRV 292 : KFIN=MIN(NEL,KDEB+LRV-1) 293 : C 294 : DO 8002 K=KDEB,KFIN 295 : KP=K-KDEB+1 296 : NK=K+K0 297 : NK1=(1-IND1)*(NK-1)+1 298 : ALFT(KP)=ALFE(NK1)+ZERMA 299 : AIRE(KP)=VOLU(NK) 300 : AL(KP)=COTE(NK,1) 301 : AH(KP)=COTE(NK,2) 302 : AP(KP)=COTE(NK,3) 303 : 8002 CONTINUE 304 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 305 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 306 : DO 8003 K=KDEB,KFIN 307 : KP=K-KDEB+1 308 : NK=K+K0 309 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 310 : 8003 CONTINUE 311 : ENDIF 312 : 313 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 314 : 315 : DO 8006 I=1,NP 316 : DO 8006 K=KDEB,KFIN 317 : KP=K-KDEB+1 318 : NF=IPADF(LE(I,K)) 319 : NU=IPADU(LE(I,K)) 320 : NI=IPADI(LE(I,K)) 321 : NFU=(1-INDU)*(NU-1)+1 322 : UIX(KP,I)=UN(NFU,1) 323 : UIY(KP,I)=UN(NFU,2) 324 : UIZ(KP,I)=UN(NFU,3) 325 : TETAC(KP,I)=HRN(NI) 326 : TETAD(KP,I)=TN(NF) 327 : 8006 CONTINUE 328 : 328 : CALL KSUPG1(WT,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV, 329 : &HRN,IPADI,UN,ALFT,NPTU,NEL,NP,DRR,HR,FN, 330 : &AIRE,AL,AH,AP,IDCEN,IPADU,LE,CUB8,IKOMP, 331 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 332 : 333 : C 334 : C Initialisation de la variable d'accumulation SBF au terme source 335 : C M 336 : IF(IKOMP.EQ.0)THEN 337 : 338 : DO 80021 K=KDEB,KFIN 339 : KP=K-KDEB+1 340 : NK=K+K0 341 : NKS=(1-IKS)*(NK-1)+1 342 : QT(KP)=QE(NKS) 343 : 80021 CONTINUE 344 : 345 : IF(IPG.EQ.0)THEN 346 : DO 80062 I=1,NP 347 : DO 80062 K=KDEB,KFIN 348 : KP=K-KDEB+1 349 : SBF(KP,I)=-QT(KP)*DRR(I,K) 350 : 80062 CONTINUE 351 : ELSE 352 : DO 81062 I=1,NP 353 : DO 81062 K=KDEB,KFIN 354 : KP=K-KDEB+1 355 : SBF(KP,I)=-QT(KP)*WT(KP,I) 356 : 81062 CONTINUE 357 : ENDIF 358 : 359 : ELSEIF(IKOMP.EQ.1)THEN 360 : 361 : DO 80023 K=KDEB,KFIN 362 : KP=K-KDEB+1 363 : NK=K+K0 364 : NKS=(1-IKS)*(NK-1)+1 365 : QT(KP)=QE(NKS) 366 : 80023 CONTINUE 367 : 368 : IF(IPG.EQ.0)THEN 369 : DO 80064 I=1,NP 370 : DO 80064 K=KDEB,KFIN 371 : KP=K-KDEB+1 372 : SBF(KP,I)=-QT(KP)*DRR(I,K) 373 : 80064 CONTINUE 374 : ELSE 375 : DO 81064 I=1,NP 376 : DO 81064 K=KDEB,KFIN 377 : KP=K-KDEB+1 378 : SBF(KP,I)=-QT(KP)*WT(KP,I) 379 : 81064 CONTINUE 380 : ENDIF 381 : 382 : ENDIF 383 : 384 : C Le coeur du calcul ... 385 : 386 : IF(IKOMP.EQ.0)THEN 387 : 388 : DO 8014 I=1,NP 389 : DO 8014 J= 1,NP 390 : DO 8014 K=KDEB,KFIN 391 : KP=K-KDEB+1 392 : 393 : ZVGG=AIRE(KP)*CHGLP(KP)*QGGT(J,I) 394 : 395 : ZVGT=AIRE(KP)*( 396 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 397 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 398 : &+ HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 399 : &+ CHGLD(KP)*QGGT(J,I) 400 : 401 : V2=UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2)+UMI(KP,3)*HR(K,J,3) 402 : 403 : SBF(KP,I)=SBF(KP,I) 404 : & +TETAC(KP,J)*(ZVGG+V2*WT(KP,I))+TETAD(KP,J)*ZVGT 405 : 406 : 8014 CONTINUE 407 : 408 : ELSEIF(IKOMP.EQ.1)THEN 409 : 410 : DO 8015 I=1,NP 411 : DO 8015 J= 1,NP 412 : DO 8015 K=KDEB,KFIN 413 : 414 : KP=K-KDEB+1 415 : 416 : ZVGG=AIRE(KP)*CHGLP(KP)*QGGT(J,I) 417 : 418 : ZVGT=AIRE(KP)*( 419 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 420 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 421 : &+ HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 422 : &+ CHGLD(KP)*QGGT(J,I) 423 : 424 : V2=UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+UIZ(KP,J)*HR(K,J,3) 425 : 426 : SBF(KP,I)=SBF(KP,I) 427 : & +TETAC(KP,J)*(ZVGG+V2*WT(KP,I))+TETAD(KP,J)*ZVGT 428 : 429 : 8015 CONTINUE 430 : 431 : ENDIF 432 : 433 : C 434 : C Fin de l'accumulation dans SBF. 435 : C On ajoute ces incréments G. 436 : C 437 : DO 8017 I=1,NP 438 : DO 8017 K=KDEB,KFIN 439 : KP=K-KDEB+1 440 : NF=IPADI(LE(I,K)) 441 : G(NF) = G(NF)+SBF(KP,I) 442 : 8017 CONTINUE 443 : 444 : 8001 CONTINUE 445 : 446 : C WRITE(6,*)' G DANS YCTSCL ' 447 : C WRITE(6,1984)(M,G(M),M=1,NPTD) 448 : 448 : C CALL ARRET(0) 449 : IPAS=1 450 : RETURN 451 : 1002 FORMAT(10(1X,1PE11.4)) 452 : END 453 : 454 : 455 : 456 : 457 : 458 : 459 : 460 :
© Cast3M 2003 - All rights reserved.
Disclaimer