1 : C SAUTPA SOURCE CHAT 05/01/13 03:11:37 5004 2 : SUBROUTINE SAUTPA(EPSILO,ICHGZ,JREBO,XIREB,XNREB 3 : $ ,NDIM,ITY1,ITYG,NOEL1,NOUN1,DIAM,DIFL,DIFT,TDEP,XDEP2 4 : $ ,UELEM,DTREEL,UREEL,U,TARI,XCON2,XARI2,IZSH,Z 5 : $ ,IREBCO,IREBDI,NTEST) 6 : * 7 : ************************************************************************* 8 : *** SP 'SAUTPA' : fait avancer la particule apres un pas de tps 'DTREEL' 9 : *** par convection (analytique ou explicite) + dispersion (explicite). 10 : *** 11 : *** APPELES 1 = 'TDRAND' 12 : *** APPELES 2 = 'REEREF', 'REFREE', 'JCBIEN', 'SAUCO1', SAUCO2', 13 : *** 'SAUDIF', 'NORMVI' 14 : *** 15 : *** E = 'EPSILO' marge relative acceptée % 0 (cause machine) 16 : *** 'ICHGZ' vaut 1 si saut précédent effectif, 0 sinon 17 : *** 'JREBO' n° local face impermeable ou se trouve particule, -1 sinon 18 : *** 'XIREB' pt d'impact sur la face impermeable 19 : *** 'XNREB' vecteur normal à la face impermeable 20 : *** 'NDIM' dimension de l'espace 21 : *** 'ITY1' entier caracterisant le type de l'element 22 : *** 'ITYG' entier caracterisant la geometrie de l'element 23 : *** 'NOEL1' nombre de noeuds de l'element traversé 24 : *** 'NOUN1' nombre de flux de l'element traversé 25 : *** 'DIAM' "longueur caracteristique" de l'element traversé 26 : *** 'DIFL' coeff dispersion longitudinal ds element considere 27 : *** 'DIFT' coeff dispersion transversal ds element considere 28 : *** 'TDEP' tps reel depart avant l'avancee de la particule 29 : *** 'XDEP2' coord reelles de depart particule 30 : *** 'UELEM' valeurs des flux aux faces 31 : *** 'DTREEL' pas de tps réel pour le calcul avancée particule 32 : *** 33 : *** S = 'UREEL' vecteur "vitesse convective" ds element reel 34 : *** 'U' vecteur "vitesse diffusive" ds element reel 35 : *** 'TARI' tps reel d'arrivee apres avancee particule 36 : *** 'XCON2' coord reelles arrivee particule apres convection 37 : *** 'XARI2' coord reelles arrivee particule apres diffusion 38 : *** 'NTEST' vaut 0 si Jacobien nul pour approximation vitesse, 0 sinon 39 : *** 40 : *** E/S = 'IZSH' segment content pour l'elemt considere au pt depart : 41 : *** 'XYZL' coordonnees reelles des noeuds (E) 42 : *** 'SHP' fonctions de forme et derivees ds elemt ref (S) 43 : *** 'SHY' fonctions de base et derivees ds elemt ref (S) 44 : *** 'Z' vecteur aleatoire entre -1 et 1 pour le saut diffusif 45 : *** 'IREBCO' vaut 1 si rebond lors de la convection, 0 sinon 46 : *** 'IREBDI' vaut 1 si rebond lors de la diffusion, 0 sinon 47 : *** 48 : *** Auteur Cyril Nou 49 : ************************************************************************* 50 : * 51 : IMPLICIT INTEGER(I-N) 52 : IMPLICIT REAL*8 (A-H,O-Z) 53 : SEGMENT IZSH 54 : REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) 55 : ENDSEGMENT 56 : DIMENSION XDEP(3),XDEP2(3),XCON2(3),XARI2(3) 57 : DIMENSION UELEM(NOUN1),UREEL(3),U(3),XJAC(3,3),Z(3) 58 : DIMENSION XIREB(3),XNREB(3),XITEST(3),XNTEST(3),X(3) 59 : 60 : C write(6,*)' uelem ', (UELEM(ii),ii=1,NOUN1) 61 : ****************************** 62 : *** INITIALISATION DU SAUT *** 63 : ****************************** 64 : 65 : *** calcul position de départ de référence 65 : CALL REEREF(NDIM,ITY1,NOEL1,IZSH,XDEP2,XDEP) 66 : *** initialisation à position de départ 67 : TARI=TDEP 68 : IREBCO=0 69 : IREBDI=0 70 : JFTEST=0 71 : DO 10 I=1,NDIM 72 : XCON2(I)=XDEP2(I) 73 : XARI2(I)=XDEP2(I) 74 : 10 CONTINUE 75 : 76 : ********************** 77 : *** SAUT CONVECTIF *** 78 : ********************** 79 : 80 : *** cas ou le saut convectif est calculé explicitement 80 : CALL SAUCON(NDIM,ITY1,ITYG,NOEL1,NOUN1 81 : $ ,DIAM,UELEM,XDEP,XDEP2,DTREEL,XCON2,IZSH,UREEL,NTEST) 82 : *** cas ou Jacobien nul lors approximation vitesse 83 : IF (NTEST.EQ.0) RETURN 84 : *** calcul de la norme de la vitesse 84 : CALL NORMVI(NDIM,UREEL,XNORM,UXY) 85 : ******************************************************** 86 : *** CAS OU TRAVERSEE FACE IMPERMEABLE PDT CONVECTION *** 87 : ******************************************************** 88 : 88 : CALL LIEUPT(0.D0,NDIM,ITYG,XCON2,IZSH,ITEST) 89 : IF (ITEST.EQ.0) THEN 89 : CALL CHOINT(NDIM,ITYG,XDEP2 90 : $ ,XDEP2,XCON2,IZSH,XITEST,XNTEST,KTEST,JFTEST) 91 : IF ((KTEST.EQ.2).AND.(JFTEST.EQ.JREBO)) THEN 92 : IREBCO=1 93 : COEFC=SCVECT(UREEL,XNREB,NDIM) 94 : COEFC=COEFC*DTREEL 95 : DO 30 I=1,NDIM 96 : XCON2(I)=XCON2(I)-IREBCO*COEFC*XNREB(I) 97 : 30 CONTINUE 98 : ENDIF 99 : ENDIF 100 : 101 : ********************************* 102 : *** SAUT DIFFUSIF (EXPLICITE) *** 103 : ********************************* 104 : 105 : IF (ICHGZ.EQ.1) THEN 106 : *** on retire vecteur aleatoire car saut précédent effectif 107 : DO 20 I=1,3 107 : CALL TDRAND(Z(I)) 108 : Z(I)=(2.D0*Z(I))-1.D0 109 : 20 CONTINUE 110 : ENDIF 110 : CALL SAUDIF(NDIM,Z,XCON2,UREEL 111 : $ ,DTREEL,XNORM,UXY,DIFL,DIFT,U,XARI2) 112 : 113 : ******************************************************* 114 : *** CAS OU TRAVERSEE FACE IMPERMEABLE PDT DIFFUSION *** 115 : ******************************************************* 116 : 116 : CALL LIEUPT(0.D0,NDIM,ITYG,XARI2,IZSH,JTEST) 117 : IF (JTEST.EQ.0) THEN 117 : CALL CHOINT(NDIM,ITYG,XCON2 118 : $ ,XCON2,XARI2,IZSH,XITEST,XNTEST,KTEST,JFTEST) 119 : IF ((KTEST.EQ.2).AND.(JFTEST.EQ.JREBO)) THEN 120 : IREBDI=1 121 : DO 40 I=1,NDIM 122 : XIREB(I)=XITEST(I) 123 : XNREB(I)=XNTEST(I) 124 : 40 CONTINUE 125 : DO 50 I=1,NDIM 126 : X(I)=XARI2(I)-XIREB(I) 127 : 50 CONTINUE 128 : COEFD=SCVECT(X,XNREB,NDIM) 129 : DO 60 I=1,NDIM 130 : XARI2(I)=XARI2(I)-2*IREBDI*COEFD*XNREB(I) 131 : 60 CONTINUE 132 : ENDIF 133 : ENDIF 134 : 135 : ******************************* 136 : *** TEMPS ARRIVEE PARTICULE *** 137 : ******************************* 138 : 139 : TARI=TDEP+DTREEL 140 : RETURN 141 : END 142 : 143 :
© Cast3M 2003 - All rights reserved.
Disclaimer