1 : C TRAVER SOURCE CHAT 05/01/13 03:45:31 5004 2 : SUBROUTINE TRAVER(DIAM,EPSILO,NSAUV,MLREE6,TMIN,NDIM 3 : $ ,ITY1,ITYG,IEL1,IELL,NOEL1,NOUN1,DIFL,DIFT,DTREEL,IZUN 4 : $ ,IZSH,TDEP,XDEP2,NPOS,ITER,IPARPO,KSAUV,DTSTOC,DTCUMU 5 : $ ,ICHGZ,Z,JFACE,JREBO,XIREB,XNREB,TARI,XARI2) 6 : 7 : ************************************************************************* 8 : *** SP 'TRAVER' : fait avancer la particule par convection + diffusion 9 : *** entre la face de depart et la face de sortie de l'element traverse 10 : *** avec un eventuel recalibrage du pas de tps associé au dernier saut 11 : *** afin d'arriver sur une face. 12 : *** 13 : *** APPELES 1 = aucun 14 : *** APPELES 2 = 'SAUTPA','LIEUPT', 'TRJTIN','CHOINT','CALIDT','MODCAL' 15 : *** 16 : *** E = 'DIAM' "longueur caracteristique" de l'element traversé 17 : *** 'EPSILO' marge relative acceptée position % element 18 : *** 'NSAUV' taille de la liste des tps de sauvegarde 19 : *** 'MLREE6' liste des tps de sauvegarde 20 : *** 'TMIN' instant de depart du lacher de la particule 21 : *** 'NDIM' dimension de l'espace 22 : *** 'ITY1' entier caracterisant le type de l'element 23 : *** 'ITYG' entier caracterisant la geometrie de l'element 24 : *** 'IEL1' n° global de l'element considéré 25 : *** 'IELL' n° local (ds sous-maillage) de l'element considéré 26 : *** 'NOEL1' nombre de noeuds de l'element traversé 27 : *** 'NOUN1' nombre de flux de l'element traversé 28 : *** 'DIFL' coeff dispersion longitudinal ds element considere 29 : *** 'DIFT' coeff dispersion transversal ds element considere 30 : *** 'DTREEL' pas de tps réel pour le calcul avancée particule 31 : *** 'IZUN' segmt content flux aux faces % ss-maillage de 'IELL' 32 : *** 'IZSH' segment content coord reelles noeuds, fcts forme et base 33 : *** 'TDEP' tps reel départ particule 34 : *** 'XDEP2' position réelle de départ particule 35 : *** 36 : *** E/S = 'NPOS' taille maximale des tableaux du segment 'IPARPO' 37 : *** 'ITER' indice des tableaux de 'IPARPO' pour la sauvegarde 38 : *** 'IPARPO' segment ou sont sauvegardés pts trajectoire particule 39 : *** 'KSAUV' indice liste des tps de sauvegarde considéré 40 : *** 'DTSTOC' pas de tps de sauvegarde considéré 41 : *** 'DTCUMU' cumul des pas de tps entre deux sauvegardes 42 : *** 'ICHGZ' vaut 1 si saut précédent effectif, 0 sinon 43 : *** 'Z' vecteur aleatoire pour le saut diffusif 44 : *** 'JFACE' n° local de la face traversée de l'element considéré 45 : *** 'JREBO' n° local face impermeable ou se trouve particule, -1 sinon 46 : *** 'XIREB' pt d'impact sur la face impermeable 47 : *** 'XNREB' vecteur normal à la face impermeable 48 : *** 49 : *** S = 'TARI' tps reel ecoulé jusqu'à interface de sortie de l'element 50 : *** 'XARI2' coord reelles arrivee particule (interface de sortie) 51 : *** 52 : *** Rq : critere d'acceptabilité position % elemt -> EPSILO (utilisé 53 : *** aussi pour critere acceptabilité dvlpt limité (calibrage). 54 : *** 55 : *** Auteur Cyril Nou 56 : ************************************************************************* 57 : 58 : 59 : *** declaration des segments et tableaux utilisés dans sp 'TRAVER' 60 : IMPLICIT INTEGER(I-N) 61 : IMPLICIT REAL*8 (A-H,O-Z) 62 : -INC SMLREEL 63 : POINTEUR MLREE6.MLREEL 64 : SEGMENT IZSH 65 : REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) 66 : ENDSEGMENT 67 : SEGMENT IZUN 68 : REAL*8 UN(I1,I2,I3) 69 : ENDSEGMENT 70 : SEGMENT IPARPO 71 : INTEGER NAPAR(NPOS),NUMP(NPOS) 72 : REAL*8 CREF(NDIM,NPOS),TPAR(NPOS) 73 : ENDSEGMENT 74 : DIMENSION XDEP(3),XDEP2(3),XCON(3),XCON2(3),XARI(3),XARI2(3) 75 : DIMENSION XCON3(3),XARI3(3) 76 : DIMENSION UREEL(3),U(3),UZERO(3),XJAC(3,3),Z(3) 77 : DIMENSION XINTER(3),XNORMA(3),XIREB(3),XNREB(3),XC2(3),XA2(3) 78 : DIMENSION IFACE(6),XINT(3,6),XN(3,6) 79 : *** 'ICALIB' vaut 1 si calibrage applicable, 0 sinon 80 : ICALIB=0 81 : *** 'DTTEST' pas de tps utilisé pour traversee (hors calibrage) 82 : DTTEST=DTREEL 83 : NCOMPT=0 83 : CALL INITD(UZERO,3,0.D0) 84 : *********************************************************** 85 : *** AVANCEE DE LA PARTICULE APRES PAS DE TEMPS 'DTTEST' *** 86 : *********************************************************** 87 : 10 LTEST=1 88 : *** saut convectif + diffusif dans la maille 88 : CALL SAUTPA(EPSILO,ICHGZ,JREBO,XIREB,XNREB,NDIM,ITY1 89 : $ ,ITYG,NOEL1,NOUN1,DIAM,DIFL,DIFT,TDEP,XDEP2,UN(1,1,IELL) 90 : $ ,DTTEST,UREEL,U,TARI,XCON2,XARI2,IZSH,Z,IREBCO,IREBDI,NTEST) 91 : *** cas ou Jacobien nul lors approximation vitesse efmh 92 : IF (NTEST.EQ.0) THEN 93 : IEL1=-1 94 : RETURN 95 : ENDIF 96 : *** localisation particule (à 0 pres car intersection inconnue) 96 : CALL LIEUPT(0.D0,NDIM,ITYG,XARI2,IZSH,ITEST) 97 : 98 : **************************************************************** 99 : *** CAS OU PARTICULE TJS DS MAILLE APRES PAS DE TPS 'DTTEST' *** 100 : **************************************************************** 101 : 102 : IF (ITEST.EQ.1) THEN 102 : CALL TRJTIN(EPSILO,NDIM,IEL1,DTTEST,TARI,XARI2,NSAUV,MLREE6,TMIN 103 : $ ,NPOS,ITER,IPARPO,KSAUV,DTSTOC,DTCUMU,TDEP,XDEP2,ICALIB,ICHGZ) 104 : *** la particule n'est plus sur une face impermeable 104 : CALL LIEUPT(-EPSILO,NDIM,ITYG,XARI2,IZSH,ITEST) 105 : IF (ITEST.EQ.1) THEN 106 : JREBO=-1 107 : IREBCO=0 108 : IREBDI=0 109 : DO 50 I=1,3 110 : XIREB(I)=0.D0 111 : XNREB(I)=0.D0 112 : 50 CONTINUE 113 : DTTEST=DTREEL 114 : ENDIF 115 : NCOMPT=NCOMPT+1 116 : GOTO 10 117 : ENDIF 118 : 119 : ************************************************************* 120 : *** CAS OU PARTICULE HORS MAILLE APRES PAS DE TPS 'DTTEST *** 121 : ************************************************************* 122 : 123 : *** cas ou 1er saut nul (conv compensée par diff au saut précédent) 124 : IF ((DTTEST.LT.(EPSILO*EPSILO*DTREEL)).AND.(NCOMPT.EQ.0)) THEN 125 : TARI=TDEP 126 : DO 20 I=1,NDIM 127 : XARI2(I)=XDEP2(I) 128 : 20 CONTINUE 129 : RETURN 130 : ENDIF 131 : *** division pas de tps pour avoir au moins un saut ds la maille 132 : IF (ICALIB.EQ.0) THEN 133 : DTTEST=1.D-1*DTTEST 134 : ICHGZ=0 135 : GOTO 10 136 : ENDIF 137 : *** determination intersection utilisée pour calibrage particule 138 : 30 CONTINUE 138 : *30 CALL CHOINT(NDIM,ITYG,XDEP2,XCON2 139 : * $ ,XARI2,IZSH,XINTER,XNORMA,JTEST,JFACE) 140 : *** on postule pas de traversee au depart 141 : JTEST=0 142 : 142 : CALL FACTRA(NDIM,ITYG,XCON2,XARI2,IZSH,IFACE,XINT,XN,NBFAC) 143 : 35 CONTINUE 144 : IF (NBFAC.NE.0) THEN 145 : *** cas ou face(s) traversee(s) lors du saut diffusif 146 : JFAC1=JFACE 146 : CALL INTPRO(NDIM,XARI2,NBFAC,IFACE,XINT,XN,JFAC1,XINTER,XNORMA) 146 : CALL CALIDT(EPSILO,NDIM,XINTER,XNORMA,XIREB,XNREB,IREBCO,IREBDI 147 : $ ,XDEP2,TDEP,UREEL,U,DTTEST,LTEST,DTNEW,XCON3,XARI3,TARI,KTEST) 148 : JTEST=2 149 : IF(KTEST.EQ.0)THEN 150 : *** la face trouvée ne convient pas on l'enleve de la liste 151 : IF(NBFAC.GT.1)THEN 152 : IJ1=0 153 : DO 37 IJ=1,NBFAC-1 154 : IF( JFACE.EQ.IFACE(IJ))IJ1=1 155 : IFACE(IJ)=IFACE(IJ+IJ1) 156 : DO 36 II=1,NDIM 157 : XINT(II,IJ)=XINT(II,IJ+IJ1) 158 : XN(II,IJ)=XN(II,IJ+IJ1) 159 : 36 CONTINUE 160 : 37 CONTINUE 161 : ENDIF 162 : NBFAC=NBFAC-1 163 : JTEST=0 164 : GO TO 35 165 : ELSE 165 : CALL RSETD(XCON2,XCON3,NDIM) 165 : CALL RSETD(XARI2,XARI3,NDIM) 166 : JFACE=JFAC1 167 : ENDIF 168 : ELSE 168 : CALL FACTRA(NDIM,ITYG,XDEP2,XCON2,IZSH,IFACE,XINT,XN,NBFAC) 169 : IF (NBFAC.NE.0) THEN 170 : *** cas ou face(s) traversee(s) lors du saut convectif 170 : CALL INTPRO(NDIM,XCON2,NBFAC 171 : $ ,IFACE,XINT,XN,JFACE,XINTER,XNORMA) 171 : CALL CALIDT(EPSILO,NDIM,XINTER,XNORMA,XIREB,XNREB,IREBCO,IREBDI 172 : $,XDEP2,TDEP,UREEL,UZERO,DTTEST,LTEST,DTNEW,XCON2,XARI2,TARI,KTEST) 173 : JTEST=1 174 : ENDIF 175 : ENDIF 176 : IF (JTEST.EQ.0) THEN 177 : IEL1=-1 178 : RETURN 179 : ENDIF 180 : *** recalibrage particule sur face avec memes vecteurs vitesses 181 : ******************************************************************* 182 : *** ANALYSE DU CAS OU CALIBRAGE A REUSSI (0 <= DTNEW <= DTTEST) *** 183 : ******************************************************************* 184 : 185 : IF (KTEST.EQ.1) THEN 185 : CALL LIEUPT(EPSILO,NDIM,ITYG,XARI2,IZSH,ITEST) 186 : IF (ITEST.EQ.1) THEN 186 : CALL TRJTIN(EPSILO,NDIM,IEL1,DTNEW,TARI,XARI2 187 : $ ,NSAUV,MLREE6,TMIN,NPOS,ITER,IPARPO,KSAUV 188 : $ ,DTSTOC,DTCUMU,TDEP,XDEP2,ICALIB,ICHGZ) 189 : TARI=TDEP 190 : DO 40 I=1,NDIM 191 : XARI2(I)=XDEP2(I) 192 : 40 CONTINUE 193 : NCOMPT=NCOMPT+1 194 : RETURN 195 : ELSEIF (ITEST.EQ.0) THEN 196 : DTTEST=DTNEW 197 : LTEST=0 198 : GOTO 30 199 : ENDIF 200 : 201 : ************************************************************* 202 : *** TRAITEMENT DU CAS OU CALIBRAGE ECHOUE (DT INCOHERENT) *** 203 : ************************************************************* 204 : 205 : ELSEIF (KTEST.EQ.0) THEN 206 : IF ((IREBCO.EQ.1).OR.(IREBDI.EQ.1)) THEN 207 : IEL1=-1 208 : RETURN 209 : ENDIF 210 : *** duplication en cas d'echec du calibrage convection seule 211 : LTEST=1 212 : DO 80 I=1,NDIM 213 : XC2(I)=XCON2(I) 214 : XA2(I)=XARI2(I) 215 : 80 CONTINUE 216 : TABIS=TARI 216 : CALL MODCAL(EPSILO,NDIM,ITYG,XINTER,XNORMA 217 : $ ,XIREB,XNREB,IREBCO,IREBDI,XDEP2,TDEP,UREEL,U 218 : $ ,IZSH,DTTEST,XCON2,XARI2,TARI,DTNEW2,KTEST1,KTEST2,LTEST) 219 : *** cas ou calibrage convection seule non pertinent 220 : IF (KTEST1.EQ.0) THEN 221 : DO 90 I=1,NDIM 222 : XCON2(I)=XC2(I) 223 : XARI2(I)=XA2(I) 224 : 90 CONTINUE 225 : TARI=TABIS 225 : CALL TRJTIN(EPSILO,NDIM,IEL1,DTNEW,TARI,XARI2 226 : $ ,NSAUV,MLREE6,TMIN,NPOS,ITER,IPARPO,KSAUV 227 : $ ,DTSTOC,DTCUMU,TDEP,XDEP2,ICALIB,ICHGZ) 228 : TARI=TDEP 229 : DO 70 I=1,NDIM 230 : XARI2(I)=XDEP2(I) 231 : 70 CONTINUE 232 : NCOMPT=NCOMPT+1 233 : RETURN 234 : ENDIF 235 : *** cas ou calibrage convection seule pertinent 236 : IF (KTEST2.EQ.1) THEN 236 : CALL TRJTIN(EPSILO,NDIM,IEL1,DTNEW2,TARI,XARI2 237 : $ ,NSAUV,MLREE6,TMIN,NPOS,ITER,IPARPO,KSAUV 238 : $ ,DTSTOC,DTCUMU,TDEP,XDEP2,ICALIB,ICHGZ) 239 : *** la particule n'est plus sur une face impermeable 239 : CALL LIEUPT(-EPSILO,NDIM,ITYG,XARI2,IZSH,ITEST) 240 : IF (ITEST.EQ.1) THEN 241 : JREBO=-1 242 : IREBCO=0 243 : IREBDI=0 244 : DO 60 I=1,3 245 : XIREB(I)=0.D0 246 : XNREB(I)=0.D0 247 : 60 CONTINUE 248 : ENDIF 249 : NCOMPT=NCOMPT+1 250 : GOTO 10 251 : ELSEIF (KTEST2.EQ.0) THEN 252 : DTTEST=DTNEW2 253 : GOTO 30 254 : ENDIF 255 : ENDIF 256 : 257 : RETURN 258 : END 259 : 260 : 261 : 262 : 263 : 264 : 265 :
© Cast3M 2003 - All rights reserved.
Disclaimer