1 : C TRJAV3 SOURCE PV 09/03/13 21:17:06 6328 2 : SUBROUTINE TRJAV3(EPSILO,NSAUV,MLREE6 3 : $ ,MELVA1,MELVA2,MELVA3,IPT9,IZVIT,IZPART,IPART 4 : $ ,TMIN,TMAX,MELEME,IELTFA,IZCENT,IFACEL,IPARPO,IZSH,SORTIE) 5 : **************************************************************************** 6 : *** SP 'TRJAV3' : fait avancer une particule n° 'IPART' jsuqu'à instant 'TMAX' 7 : *** ou jusqu'à ce que particule sorte du domaine. Pour chaque intersection 8 : *** element-trajectoire, les coordonnées réelles et le tps courant réel 9 : *** associés à la particule n° 'IPART' sont stockés via 'IPARPO'. 10 : *** 11 : *** E = 'EPSILO' erreur précision calul relative acceptee (cause machine) 12 : *** 'NSAUV' taille de la liste des tps de sauvegarde 13 : *** 'MLREE6' liste des tps de sauvegarde 14 : *** 'MELVA1' segment content coeffs disper intrinseques longitudinaux 15 : *** 'MELVA2' segment content coeffs disper intrinseques transversaux 16 : *** 'MELVA3' segment content coeffs diff moleculaire 17 : *** 'IPT9' pteur sur maillage decrivant faces impermeables du domaine 18 : *** 'IZVIT' segment decrivant les vitesses (<- 'TRJVIT' OU 'TRJFLU') 19 : *** 'IZPART' position initiale en coordonnees reference de la particule 20 : *** 'IPART' n° de la particule concernee 21 : *** 'TMIN' instant de depart de la particule 22 : *** 'TMAX' instant à ne pas depasser pour la particule 23 : *** 'MELEME' pteur sur le maillage du domaine etudie 24 : *** 'IELTFA' pteur sur la table "DOMAINE.ELTFA" 25 : *** 'IZCENT' pteur sur la table "DOMAINE.CENTRE" 26 : *** 'IFACEL' pteur sur la table " DOMAINE.FACEL" 27 : *** S = 'IPARPO' position reelle et tps correspondant de la particule aux 28 : *** interfaces des elements traversees. 29 : *** 'IZSH' segmt content fcts forme,base et coord reelles des noeuds 30 : *** 'SORTIE' tableau contenant respectivement tps et coords de sortie 31 : ****************************************************************************** 32 : 33 : *** ORIGINE = SP 'TRJAVA' modifié par PATRICK MEYNIEL 12/99 34 : *** CYRIL NOU 08/00 35 : *** 36 : *********************************************************** 37 : *** DECLARATION/ACIVATION SEGMENTS UTILISES DS 'TRJAV2' *** 38 : *********************************************************** 39 : 40 : IMPLICIT INTEGER(I-N) 41 : IMPLICIT REAL*8 (A-H,O-Z) 42 : -INC CCOPTIO 43 : -INC SMCOORD 44 : -INC SMELEME 45 : -INC SMCHAML 46 : -INC SMLREEL 47 : POINTEUR IZCENT.MELEME,IELTFA.MELEME,IFACEL.MELEME 48 : POINTEUR MLREE6.MLREEL 49 : SEGMENT IZPART 50 : INTEGER NLEPA(NPART),NUMPA(NPART) 51 : REAL*8 COORPA(NDIM,NPART) 52 : ENDSEGMENT 53 : SEGMENT IPARPO 54 : INTEGER NAPAR(NPOS),NUMP(NPOS) 55 : REAL*8 CREF(NDIM,NPOS),TPAR(NPOS) 56 : ENDSEGMENT 57 : SEGMENT IZSH 58 : REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) 59 : ENDSEGMENT 60 : SEGMENT IZVIT 61 : REAL*8 TEMTRA(NVIPT) 62 : INTEGER IPUN(NBS),IDUN(NBS),IPVPT(NVIPT),IFORML 63 : ENDSEGMENT 64 : *** declaration de segment utilise dans le cas du transitoire 65 : *** 'IPUN1' = pteurs sur les 'IZUN', 'IPUMAX' sur le 'IZUMAX' 66 : SEGMENT IZVPT 67 : INTEGER IPUN1(NBS),IPUMAX 68 : ENDSEGMENT 69 : *** declaration de segment contenant les vitesses ou flux 70 : SEGMENT IZUN 71 : REAL*8 UN(I1,I2,I3) 72 : ENDSEGMENT 73 : *** declaration de segment contenant la vitesse ou flux max par element 74 : SEGMENT IZUMAX 75 : REAL*8 UMAX(NBREL) 76 : ENDSEGMENT 77 : *** variables pour positions, stockage, aleatoire et rebond 78 : DIMENSION XDEP(3),XDEP2(3),XARI(3),XARI2(3) 79 : DIMENSION XSTOC(3),Z(3),XIREB(3),XNREB(3) 80 : DIMENSION PT1(3),PT2(3),PT3(3),PT4(3),UCEN(3) 81 : *** contient respectivement tps de sortie et coords de sortie 82 : DIMENSION SORTIE(4) 83 : *** activation si necessaire segments liés au pb (maillage, vitesse,...) 84 : SEGACT IFACEL 85 : ******************************************************** 86 : *** INITIALISATION DU LACHER DE LA PARTICULE 'IPART' *** 87 : ******************************************************** 88 : 89 : *** initialisation trajectoire % position lacher particule 90 : * WRITE(6,*) '-> Trajectoire particule n°',IPART 90 : CALL INITPA(EPSILO,IPART,IZPART,TMIN,NSAUV,MLREE6 91 : $ ,MELEME,IPT9,NDIM,NPOS,ITER,IPARPO,IVPT,IEL1,XDEP2,TDEP 92 : $ ,JFACE,JREBO,XIREB,XNREB,ICHGZ,Z,KSAUV,DTSTOC,DTCUMU,IZSH, 93 : $ IELTFA) 94 : *** initialisation element initial % position lacher particule 95 : *** localisation, caractérisation élément content particule 'IPART' 95 : CALL TRJLOC(NDIM,MELEME,IZCENT,IELTFA,IZVIT,IVPT 96 : $ ,IEL1,TDEP,NOEL1,ITY1,ITYG,IZSH,IZUN,NOUN1,IELL,DIAM,IPT1) 97 : *** calcul des coefficients de dispersion long et trans 97 : CALL COEFDI(NDIM,ITY1,ITYG,NOEL1,NOUN1,DIAM,UN(1,1,IELL) 98 : $ ,IEL1,IZCENT,MELVA1,MELVA2,MELVA3,DIFL,DIFT,UCEN,IZSH) 99 : *** cas ou soit particule ne peut pas bouger (pas conv,disp et diff) 100 : *** soit Jacobien nul dans approximation de la vitesse efmh 101 : IF (IEL1.EQ.(-1)) GOTO 40 102 : *** choix du pas de tps pour la traversee de la maille 102 : CALL TRJCDT(NDIM,DIAM,IZVIT,IVPT,IEL1,DIFL,DIFT,UCEN,DTREEL) 103 : 104 : *********************************************************** 105 : *** TRAVERSEE ENTRE 2 INTERFACES DE L'ELEMENT CONSIDERE *** 106 : *********************************************************** 107 : 108 : 10 CONTINUE 109 : *** traversee de l'element 'IEL1' ('IELL') 109 : CALL TRAVER(DIAM,EPSILO,NSAUV,MLREE6,TMIN,NDIM 110 : $ ,ITY1,ITYG,IEL1,IELL,NOEL1,NOUN1,DIFL,DIFT,DTREEL 111 : $ ,IZUN,IZSH,TDEP,XDEP2,NPOS,ITER,IPARPO,KSAUV,DTSTOC 112 : $ ,DTCUMU,ICHGZ,Z,JFACE,JREBO,XIREB,XNREB,TARI,XARI2) 113 : *** cas ou echec recalibrage sur face element 114 : IF (IEL1.EQ.(-1)) GOTO 40 115 : 116 : ************************************************************ 117 : *** RECHERCHE DU PROCHAIN ELEMENT TRAVERSE PAR PARTICULE *** 118 : ************************************************************ 119 : 120 : *** teste si face de recalibrage est une face impermeable 120 : CALL TESTFA(EPSILO,NDIM,ITYG 121 : $ ,XARI2,IZSH,JFTEST,PT1,PT2,PT3,PT4) 121 : CALL LOCIMP(NDIM,JFACE,XARI2,PT1,PT2,PT3,PT4,IPT9 122 : $ ,JREBO,XIREB,XNREB,IEL1,IELTFA) 123 : 124 : *** chgt element à traverser 124 : CALL CHGELE(EPSILO,NDIM,JREBO,XNREB,TARI 125 : $ ,XARI2,DTREEL,IVPT,MELEME,IFACEL,IZCENT,IELTFA 126 : $ ,IZVIT,IZUN,IEL1,TDEP,XDEP2,JFACE,ICHGZ,IZSH) 127 : ******************************************************** 128 : *** DIFFERENTS CAS POSSIBLES APRES TRAVERSEE ELEMENT *** 129 : ******************************************************** 130 : 131 : *** cas ou il n'y a pas d'élément suivant 132 : IF (IEL1.EQ.0) THEN 133 : SORTIE(1)=TDEP 134 : SORTIE(2)=XDEP2(1) 135 : SORTIE(3)=XDEP2(2) 136 : SORTIE(4)=XDEP2(3) 137 : GOTO 20 138 : *** cas ou tps maximum imposé est dépassé 139 : ELSEIF ((TMAX-TDEP).LE.0.D0) THEN 140 : SORTIE(1)=TDEP 141 : SORTIE(2)=XDEP2(1) 142 : SORTIE(3)=XDEP2(2) 143 : SORTIE(4)=XDEP2(3) 144 : GOTO 20 145 : *** cas ou on refait la meme démarche dans élément suivant 146 : ELSE 147 : *** localisation, caractérisation élément content particule 'IPART' 147 : CALL TRJLOC(NDIM,MELEME,IZCENT,IELTFA,IZVIT,IVPT 148 : $ ,IEL1,TDEP,NOEL1,ITY1,ITYG,IZSH,IZUN,NOUN1,IELL,DIAM,IPT1) 149 : *** calcul des coefficients de dispersion long et trans 149 : CALL COEFDI(NDIM,ITY1,ITYG,NOEL1,NOUN1,DIAM,UN(1,1,IELL) 150 : $ ,IEL1,IZCENT,MELVA1,MELVA2,MELVA3,DIFL,DIFT,UCEN,IZSH) 151 : *** choix du pas de tps pour la traversee de la maille 152 : IF (IEL1.NE.(-1)) 152 : *CALL TRJCDT(NDIM,DIAM,IZVIT,IVPT,IEL1,DIFL,DIFT,UCEN,DTREEL) 153 : GOTO 10 154 : ENDIF 155 : 156 : ************************************************************************* 157 : *** REAJUSTEMENT DE 'IPARPO' (POUR LE TRACE VIA OPERATEUR 'TRTRAJEC') *** 158 : ************************************************************************* 159 : 160 : *** recuperation de variables de travail pour la sauvegarde 161 : 20 TSTOC=TPAR(ITER) 162 : IESTO=NAPAR(ITER) 163 : DO 30 I=1,NDIM 164 : XSTOC(I)=CREF(I,ITER) 165 : 30 CONTINUE 166 : *** on supprime les pts sauvegardés dont tps sauvegarde > 'TMAX' 167 : IF (TSTOC.GT.((1+EPSILO)*TMAX)) THEN 168 : ITER=ITER-1 169 : GOTO 20 170 : ENDIF 171 : *** au moins 2 pts sauvegardés car 'TRTRAJEC' => maillage 'SEG2' 171 : IF(ITER.EQ.1)CALL TRJSTO(NDIM,IESTO,NPOS,ITER,IPARPO,TSTOC,XSTOC) 172 : IF(ABS((TPAR(ITER)-TDEP)/TDEP).GT.EPSILO) THEN 173 : NPOS=ITER+1 174 : SEGADJ IPARPO 175 : NAPAR(NPOS)=IESTO 176 : TPAR(NPOS)=TDEP 177 : DO 35 III=1,IDIM 178 : CREF(III,NPOS)=XDEP2(III) 179 : 35 CONTINUE 180 : ELSE 181 : *** on reajuste la dimension du segment 'IPARPO' 182 : NPOS=ITER 183 : SEGADJ IPARPO 184 : ENDIF 185 : *** on sort du module pour traiter la particule suivante 186 : RETURN 187 : 188 : ****************************************************** 189 : *** CAS D'ERREUR LORS DE LA TRAVERSEE DE L'ELEMENT *** 190 : ****************************************************** 191 : 192 : 40 WRITE(6,*) '==================================================' 193 : WRITE(6,*) 'LA PARTICULE N°',IPART,' VA ETRE ELIMINEE POUR' 194 : WRITE(6,*) 'OPERATEUR "TRAJ" POUR UNE DES RAISONS SUIVANTES : ' 195 : WRITE(6,*) 196 : WRITE(6,*) 'a) ni convection, ni dispersion, ni diffusion ' 197 : WRITE(6,*) 'b) Jacobien nul pour approximation vitesse efmh ' 198 : WRITE(6,*) 'c) echec recalibrage particule sur face élément ' 199 : WRITE(6,*) 'd) echec localisation particule ds un élément ' 200 : WRITE(6,*) 201 : WRITE(6,*) 'OPERATEUR "TRAJ" CONTINUE POUR PARTICULE SUIVANTE ' 202 : WRITE(6,*) '==================================================' 203 : 204 : *** indication pour elimination particule 205 : SORTIE(1)=-1.D0 206 : 207 : RETURN 208 : END 209 : 210 : 211 : 212 : 213 : 214 : 215 : 216 :
© Cast3M 2003 - All rights reserved.
Disclaimer