zctsca
C ZCTSCA SOURCE CHAT 05/01/13 04:21:34 5004 SUBROUTINE ZCTSCA C C VERSION VECTORISEE C C Les éléments sont groupés en paquets de LRV éléments, LRV étant C la longueur des registres vectoriels de la machine cible, i.e C 64 sur Cray, 128 ou 256 sur IBM 3090VF. On promène une fenêtre C de longueur LRV sur la boucle générale de longueur NEL. C C! Attention : IAXI et RPG ne sont pas utilisées... & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI, & IPADL,IKOMP,IKAS,IPADU,IPADD, & ALFE,IND1,UN,INDU,NPTU,TN,QE,IKS, & HRN,G,ALT,SGT, & VOLU,COTE,NELZ,IDCEN, & DTM1,DT,DTT1,DTT2,NUEL,DIAEL) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C C CE SP DISCRETISE UNE EQUATION GENERALE DE TRANSPORT-DIFFUSION AVEC C SOURCE. C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 C LES OPERATEURS SONT "SOUS-INTEGRES" C C C APPELE PAR YTSCAL C C C*********************************************************************** C -INC PPARAM -INC CCOPTIO -INC CCVQUA4 -INC CCREEL C C Longueur des registres vectoriels de la machine cible C On prend 64 pour ne pas augmenter la taille des tableaux C nécessaires à la vectorisation. C PARAMETER(LRV=64) DIMENSION UN(NPTU,IES),HRN(NPTD),TN(NPTD) DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*) DIMENSION ALFE(*),ALT(*) DIMENSION IPADL(*),LE(NP,*),IPADU(*),IPADD(*) DIMENSION HR(NEL,NP,IES),RPG(*),DRR(NP,NEL) DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) DIMENSION AIRE(LRV) DIMENSION AL2(LRV),AH2(LRV),AP2(LRV) DIMENSION AL(LRV),AH(LRV),AP(LRV) DIMENSION ALFT(LRV),QT(LRV) DIMENSION XMB(LRV),XMH(LRV) DIMENSION CFM(LRV) C DIMENSION CF1(LRV),CF2(LRV),CF3(LRV) DIMENSION DR(LRV,9) DIMENSION UM(LRV),UP(LRV) DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9) DIMENSION TETAC(LRV,9),TETAD(LRV,9) DIMENSION UMI(LRV,3),UPI(LRV,3) DIMENSION CXT(LRV),CYT(LRV),CXY(LRV) DIMENSION DXT(LRV),DYT(LRV),DXY(LRV) DIMENSION CZT(LRV),CXZ(LRV),CYZ(LRV) DIMENSION DZT(LRV),DXZ(LRV),DYZ(LRV) DIMENSION SBF(LRV,9) DIMENSION GRADT(LRV,3) DIMENSION BM(LRV),BP(LRV) REAL*8 G(*) SAVE IPAS,QGGT,Q1,Q2,Q3 DATA CD/1.D0/ DATA IPAS/0/ C************************************************************************ C C INITIALISATIONS DIVERSES C C WRITE(IOIMP,*)' debut YCTSCL',IES ZERMA=1.D-20 NK=K0 C ******** C * 2D * C ******** IF(IES.EQ.3)GO TO 10 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DIFFERENCES TRIANGLE / QUADRANGLE QUA4=0.D0 IF(NP.EQ.4)QUA4=1.D0 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Calcul du nombre de paquets de LRV éléments C NNN=MOD(NEL,LRV) IF(NNN.EQ.0) NPACK=NEL/LRV IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV KPACKD=1 KPACKF=NPACK C C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** C DO 7001 KPACK=KPACKD,KPACKF C C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= C C 1.D0 Calcul des limites du paquet courant. KDEB=1+(KPACK-1)*LRV KFIN=MIN(NEL,KDEB+LRV-1) C DO 7002 K=KDEB,KFIN NK=K+K0 NK1=(1-IND1)*(NK-1)+1 7002 CONTINUE IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. & (IKOMP.EQ.1.AND.IKAS.EQ.6))THEN DO 7003 K=KDEB,KFIN NK=K+K0 7003 CONTINUE ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Initialisation des UMI avant accumulation DO 7005 K=KDEB,KFIN 7005 CONTINUE DO 7006 I=1,NP C*IBMDIR* PREFER VECTOR DO 70061 K=KDEB,KFIN NF=IPADL(LE(I,K)) NFD=IPADD(LE(I,K)) NFU=IPADU(LE(I,K)) NFU=(1-INDU)*(NFU-1)+1 70061 CONTINUE 7006 CONTINUE C C Initialisation de la variable d'accumulation SBF au terme source C C WRITE(IOIMP,*)' IKomp,ikas=',IKomp,ikas C WRITE(IOIMP,*)' IKS,IND1,INDU=',IKS,IND1,INDU IF(IKOMP.EQ.0)THEN DO 70021 K=KDEB,KFIN NK=K+K0 NKS=(1-IKS)*(NK-1)+1 70021 CONTINUE DO 70062 I=1,NP C*IBMDIR* PREFER VECTOR DO 70162 K=KDEB,KFIN 70162 CONTINUE 70062 CONTINUE ELSEIF(IKOMP.EQ.1)THEN DO 70023 K=KDEB,KFIN NK=K+K0 NKS=(1-IKS)*(NK-1)+1 70023 CONTINUE DO 70064 I=1,NP C*IBMDIR* PREFER VECTOR DO 70164 K=KDEB,KFIN 70164 CONTINUE 70064 CONTINUE ENDIF DO 7007 K=KDEB,KFIN A=SQRT(A)+XPETIT 7007 CONTINUE DO 70074 K=KDEB,KFIN 70074 CONTINUE C C AVANT DECENTREMENT C IF(IKOMP.EQ.0)THEN C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008 IF(IDCEN.EQ.1)THEN DO 70081 K=KDEB,KFIN 70081 CONTINUE ELSEIF(IDCEN.EQ.2)THEN DO 7008 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 CPT=CCP-CCT CC2=0.D0 IF(CPT.GE.0.D0)CC2=CPT 7008 CONTINUE ELSEIF(IDCEN.EQ.3)THEN DO 7018 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 7018 CONTINUE ELSEIF(IDCEN.EQ.4)THEN DT19=DTM1*0.5D0 DO 7009 K=KDEB,KFIN 7009 CONTINUE ENDIF ELSEIF(IKOMP.EQ.1)THEN C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008 IF(IDCEN.EQ.1)THEN DO 71081 K=KDEB,KFIN 71081 CONTINUE ELSEIF(IDCEN.EQ.2)THEN DO 7108 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 CPT=CCP-CCT CC2=0.D0 IF(CPT.GE.0.D0)CC2=CPT CC2 = CC2*1.3D0 7108 CONTINUE ELSEIF(IDCEN.EQ.3)THEN DO 71018 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 71018 CONTINUE ELSEIF(IDCEN.EQ.4)THEN DT19=DTM1*0.5D0 DO 71009 K=KDEB,KFIN 71009 CONTINUE ENDIF ENDIF C *********************** C C AVANT CALCUL DT C IF(IKOMP.EQ.0)THEN C*IBMDIR* PREFER SCALAR DO 7010 K=KDEB,KFIN DT0=DT DT1=2.D0/ DT2=0.5D0/ IF(DT1.LT.DT)DT=DT1 IF(DT2.LT.DT)DT=DT2 C IF(DT3.LT.DT)DT=DT3 IF(DT.NE.DT0) THEN DTT1=DT1 DTT2=DT2 C DTT3=0.D0 NUEL=K END IF 7010 CONTINUE ELSEIF(IKOMP.EQ.1)THEN C*IBMDIR* PREFER SCALAR DO 7110 K=KDEB,KFIN DT0=DT DT1=2.D0/ DT2=0.5D0/ IF(DT1.LT.DT)DT=DT1 IF(DT2.LT.DT)DT=DT2 C IF(DT3.LT.DT)DT=DT3 IF(DT.NE.DT0) THEN DTT1=DT1 DTT2=DT2 C DTT3=0.D0 NUEL=K END IF 7110 CONTINUE ENDIF C Le coeur du calcul ... IF(IKOMP.EQ.0)THEN DO 7014 I=1,NP DO 70141 J= 1,NP C*IBMDIR* PREFER VECTOR DO 70142 K=KDEB,KFIN & *VGGT(J,I)*QUA4 ) C &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) ) C &+ (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4 & *VGGT(J,I)*QUA4 ) C &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) ) C &+ ALFT(KP)/12.D0*VGGT(J,I)*QUA4 C? V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2))*DR(KP,I) 70142 CONTINUE 70141 CONTINUE 7014 CONTINUE ELSEIF(IKOMP.EQ.1)THEN DO 7015 I=1,NP DO 70151 J= 1,NP C*IBMDIR* PREFER VECTOR DO 70152 K=KDEB,KFIN &+ HR(K,I,1)*HR(K,J,2) &+ HR(K,I,2)*HR(K,J,1) & *VGGT(J,I)*QUA4 ) C &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) ) C &+ (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4 & *VGGT(J,I)*QUA4 ) C &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) ) C &+ ALFT(KP)/12.D0*VGGT(J,I)*QUA4 70152 CONTINUE 70151 CONTINUE 7015 CONTINUE ENDIF C C Fin de l'accumulation dans SBF. C On ajoute ces incréments G. C DO 7017 I=1,NP C*IBMDIR* PREFER VECTOR DO 70171 K=KDEB,KFIN NF=IPADL(LE(I,K)) 70171 CONTINUE 7017 CONTINUE 7001 CONTINUE C WRITE(IOIMP,*)' G DANS YCTSCL ' C WRITE(IOIMP,1984)(M,G(M),M=1,NPTD) 1984 FORMAT(7(1X,I4,2X,1PE11.4)) C CALL ARRET(0) IPAS=1 RETURN C ******** C * 3D * C ******** 10 CONTINUE CUB8=0.D0 IF(NP.EQ.8)CUB8=1.D0 C C Calcul du nombre de paquets de LRV éléments C NNN=MOD(NEL,LRV) IF(NNN.EQ.0) NPACK=NEL/LRV IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV KPACKD=1 KPACKF=NPACK C C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** C DO 8001 KPACK=KPACKD,KPACKF C C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= C C 1.D0 Calcul des limites du paquet courant. KDEB=1+(KPACK-1)*LRV KFIN=MIN(NEL,KDEB+LRV-1) C DO 8002 K=KDEB,KFIN NK=K+K0 NK1=(1-IND1)*(NK-1)+1 C CF1(KP)=AL(KP)*AH(KP)/AP(KP) C CF2(KP)=AL(KP)*AP(KP)/AH(KP) C CF3(KP)=AP(KP)*AH(KP)/AL(KP) 8002 CONTINUE IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. & (IKOMP.EQ.1.AND.IKAS.EQ.6))THEN DO 8003 K=KDEB,KFIN NK=K+K0 8003 CONTINUE ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Initialisation des UMI avant accumulation DO 8005 K=KDEB,KFIN 8005 CONTINUE DO 8006 I=1,NP C*IBMDIR* PREFER VECTOR DO 80061 K=KDEB,KFIN NF=IPADL(LE(I,K)) NFD=IPADD(LE(I,K)) NFU=IPADU(LE(I,K)) NFU=(1-INDU)*(NFU-1)+1 80061 CONTINUE 8006 CONTINUE C C Initialisation de la variable d'accumulation SBF au terme source C M C WRITE(IOIMP,*)' IKomp,ikas=',IKomp,ikas C WRITE(IOIMP,*)' IKS,IND1,INDU=',IKS,IND1,INDU IF(IKOMP.EQ.0)THEN DO 80021 K=KDEB,KFIN NK=K+K0 NKS=(1-IKS)*(NK-1)+1 80021 CONTINUE DO 80062 I=1,NP C*IBMDIR* PREFER VECTOR DO 80162 K=KDEB,KFIN 80162 CONTINUE 80062 CONTINUE ELSEIF(IKOMP.EQ.1)THEN DO 80023 K=KDEB,KFIN NK=K+K0 NKS=(1-IKS)*(NK-1)+1 80023 CONTINUE DO 80064 I=1,NP C*IBMDIR* PREFER VECTOR DO 80164 K=KDEB,KFIN 80164 CONTINUE 80064 CONTINUE ENDIF DO 8007 K=KDEB,KFIN A=SQRT(A)+XPETIT 8007 CONTINUE DO 80074 K=KDEB,KFIN 80074 CONTINUE C C AVANT DECENTREMENT C IF(IKOMP.EQ.0)THEN C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008 IF(IDCEN.EQ.1)THEN DO 80081 K=KDEB,KFIN 80081 CONTINUE ELSEIF(IDCEN.EQ.2)THEN DO 8008 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 CPT=CCP-CCT CC2=0.D0 IF(CPT.GE.0.D0)CC2=CPT 8008 CONTINUE ELSEIF(IDCEN.EQ.3)THEN DO 8018 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 8018 CONTINUE ELSEIF(IDCEN.EQ.4)THEN DT19=DTM1*0.5D0 DO 8009 K=KDEB,KFIN 8009 CONTINUE ENDIF ELSEIF(IKOMP.EQ.1)THEN C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008 IF(IDCEN.EQ.1)THEN DO 81081 K=KDEB,KFIN 81081 CONTINUE ELSEIF(IDCEN.EQ.2)THEN DO 81008 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 CPT=CCP-CCT CC2=0.D0 IF(CPT.GE.0.D0)CC2=CPT CC2 = CC2*1.3D0 81008 CONTINUE ELSEIF(IDCEN.EQ.3)THEN DO 81018 K=KDEB,KFIN AKSI=ALFA/3.D0 IF(ALFA.GT.3.D0)AKSI=1.D0 81018 CONTINUE ELSEIF(IDCEN.EQ.4)THEN DT19=DTM1*0.5D0 DO 81009 K=KDEB,KFIN 81009 CONTINUE ENDIF ENDIF C *********************** C C AVANT CALCUL DT C IF(IKOMP.EQ.0)THEN C*IBMDIR* PREFER SCALAR DO 8010 K=KDEB,KFIN DT0=DT DT1=2.D0/ DT2=0.5D0/ IF(DT1.LT.DT)DT=DT1 IF(DT2.LT.DT)DT=DT2 C IF(DT3.LT.DT)DT=DT3 IF(DT.NE.DT0) THEN DTT1=DT1 DTT2=DT2 C DTT3=0.D0 NUEL=K END IF 8010 CONTINUE ELSEIF(IKOMP.EQ.1)THEN C*IBMDIR* PREFER SCALAR DO 8110 K=KDEB,KFIN DT0=DT DT1=2.D0/ DT2=0.5D0/ IF(DT1.LT.DT)DT=DT1 IF(DT2.LT.DT)DT=DT2 C IF(DT3.LT.DT)DT=DT3 IF(DT.NE.DT0) THEN DTT1=DT1 DTT2=DT2 C DTT3=0.D0 NUEL=K ENDIF 8110 CONTINUE ENDIF C Le coeur du calcul ... IF(IKOMP.EQ.0)THEN DO 8014 I=1,NP DO 80141 J= 1,NP C*IBMDIR* PREFER VECTOR DO 80142 K=KDEB,KFIN C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 C &+ (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 C &+ ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 80142 CONTINUE 80141 CONTINUE 8014 CONTINUE ELSEIF(IKOMP.EQ.1)THEN DO 8015 I=1,NP DO 80151 J= 1,NP C*IBMDIR* PREFER VECTOR DO 80152 K=KDEB,KFIN C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 C &+ (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 C &+ ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 80152 CONTINUE 80151 CONTINUE 8015 CONTINUE ENDIF C C Fin de l'accumulation dans SBF. C On ajoute ces incréments G. C DO 8017 I=1,NP C*IBMDIR* PREFER VECTOR DO 80171 K=KDEB,KFIN NF=IPADL(LE(I,K)) 80171 CONTINUE 8017 CONTINUE 8001 CONTINUE C WRITE(IOIMP,*)' G DANS YCTSCL ' C WRITE(IOIMP,1984)(M,G(M),M=1,NPTD) C CALL ARRET(0) IPAS=1 RETURN 1002 FORMAT(10(1X,1PE11.4)) END
© Cast3M 2003 - Tous droits réservés.
Mentions légales