1 : C CALIDT SOURCE CHAT 05/01/12 21:46:19 5004 2 : SUBROUTINE CALIDT(EPSILO,NDIM,XINTER,XNORMA 3 : $ ,XIREB,XNREB,IREBCO,IREBDI,XDEP2,TDEP,UREEL,U 4 : $ ,DTREEL,LTEST,DTNEW,XCON2,XARI2,TARI,ITEST) 5 : ********************************************************************** 6 : *** SP 'CALIDT' : lorsque la particule sort de la maille, permet de 7 : *** recalibrer le dernier pas de tps de facon à ce que trajectoire de 8 : *** la particule s'arrete à l'interface de la maille. 9 : *** 10 : *** APPELES 1 = 'SCVECT' (fonction) 11 : *** APPELES 2 = aucun 12 : *** 13 : *** E = 'EPSILO' marge relative acceptée position % element 14 : *** 'NDIM' dimension de l'espace 15 : *** 'XINTER' pt d'intersection de la face traversee considérée 16 : *** 'XNORMA' normale associée à face traversee considérée 17 : *** 'XIREB' pt d'impact sur la face impermeable 18 : *** 'XNREB' vecteur normal à la face impermeable 19 : *** 'IREBCO' vaut 1 si rebond lors de la convection, 0 sinon 20 : *** 'IREBDI' vaut 1 si rebond lors de la diffusion, 0 sinon 21 : *** 'XDEP2' coordonnees de depart de la particule (avant convec) 22 : *** 'TDEP' temps reel de depart de la particule 23 : *** 'UREEL' vitesse dans element reel au pt de depart 24 : *** 'U' vecteur vitesse "diffusive" = XL*Z = XPASS*XLU*Z 25 : *** 'DTREEL' pas de temps reel avant recalibrage 26 : *** 27 : *** E/S = 'LTEST' vaut 1 si 1er essai calibrage, 0 si 2e essai 28 : *** 29 : *** S = 'DTNEW' pas de tps reel recalibre (le plus petit possible) 30 : *** 'XCON2' coord reelles particule apres convection apres calibrage 31 : *** 'XARI2' coord reelles particule apres diffusion apres calibrage 32 : *** 'TARI' temps reel d'arrivee de la particule 33 : *** 'ITEST' vaut 1 si calibrage donne 0 < dt < dtreel, 0 sinon 34 : *** 35 : *** Rq : pour pouvoir effectuer un dvlpt limité au 2eme ordre pour le 36 : *** calcul de delta, nous considerons le critere 4ac/b² < epsdl. 37 : *** on a alors une dependance epsdl > Cte*sqrt(EPSILO) due au 38 : *** reste du dl au 2eme ordre. 39 :
Warning: preg_replace(): Compilation failed: unmatched parentheses at offset 6 in C:\xampp\htdocs\cast3m_kit_local_fr_en_2014\sources.php on line 158
40 : *** 41 : *** ORIGINE = CYRIL NOU 42 : ********************************************************************** 43 : 44 : IMPLICIT INTEGER(I-N) 45 : IMPLICIT REAL*8 (A-H,O-Z) 46 : -INC CCREEL 47 : DIMENSION XINTER(3),XNORMA(3),XDEP2(3),UREEL(3),U(3) 48 : DIMENSION XMF(3),XXCO2(3),XXAR2(3),XCON2(3),XARI2(3) 49 : DIMENSION XIREB(3),XNREB(3),XMFR(3),X(3) 50 : ITEST=1 51 : *** initialisation des coefficients du polynome 2nd degre à resoudre 52 : DCARA2=0.D0 53 : DO 10 I=1,NDIM 54 : XMF(I)=XINTER(I)-XDEP2(I) 55 : XMFR(I)=XIREB(I)-XDEP2(I) 56 : DCARA2=DCARA2+XMF(I)**2 57 : 10 CONTINUE 58 : PDTNOR=SCVECT(XNREB,XNORMA,NDIM) 59 : COEF2A=SCVECT(UREEL,XNORMA,NDIM) 60 : COEF2B=SCVECT(UREEL,XNREB,NDIM) 61 : COEFF2=COEF2A-PDTNOR*(IREBCO+2*IREBDI)*COEF2B 62 : COEF1A=SCVECT(U,XNORMA,NDIM) 63 : COEF1B=SCVECT(U,XNREB,NDIM) 64 : COEFF1=COEF1A-2*PDTNOR*IREBDI*COEF1B 65 : COEF0A=SCVECT(XMF,XNORMA,NDIM) 66 : COEF0B=SCVECT(XMFR,XNREB,NDIM) 67 : COEFF0=-(COEF0A-2*PDTNOR*IREBDI*COEF0B) 68 : IF (ABS(4*COEFF2*COEFF0).GE.(100*SQRT(EPSILO)*(COEFF1**2))) THEN 69 : *** cas ou resolution de coeff2*X*X+coeff1*X+coeff0=0 70 : DELTA=COEFF1*COEFF1-4.D0*COEFF2*COEFF0 71 : IF (DELTA.GE.(-XZPREC*(DCARA2/DTREEL))) THEN 72 : RDELTA=SQRT(ABS(DELTA)) 73 : X1=(-COEFF1+RDELTA)/(2.D0*COEFF2) 74 : X2=(-COEFF1-RDELTA)/(2.D0*COEFF2) 75 : ELSE 76 : ITEST=0 77 : RETURN 78 : ENDIF 79 : ELSE 80 : *** cas ou devlpt limité de la quantité "racine de delta" 81 : ORD1=-COEFF0/COEFF1 82 : ORD2=-COEFF2*(COEFF0**2)/(COEFF1**3) 83 : ORD3=-2*(COEFF2**2)*(COEFF0**3)/(COEFF1**5) 84 : ORD4=-5*(COEFF2**3)*(COEFF0**4)/(COEFF1**7) 85 : ORD5=-14*(COEFF2**4)*(COEFF0**5)/(COEFF1**9) 86 : X1=ORD1+ORD2+ORD3+ORD4+ORD5 87 : IF (ABS(COEFF2).LT.(XZPREC*(SQRT(DCARA2)/DTREEL))) THEN 88 : X2=X1 89 : ELSE 90 : X2=-COEFF1/COEFF2-X1 91 : ENDIF 92 : ENDIF 93 : *** recuperation des racines > 0 94 : DT1=-1.D0 95 : DT2=-1.D0 96 : IF (X1.GE.0.D0) DT1=X1**2 97 : IF (X2.GE.0.D0) DT2=X2**2 98 : *** recuperation du plus gd dt < dtreel 99 : C nous avons considéré un recalibrage possible avec dtnew environ= dttest, 100 : C cas qui peut se produire si particule arrive a proximité de l'interface 101 : C sans recalibrage. si le meme cas de figure se produit avec dtnew environ= 0, 102 : C une marge 'EPSILO' devra etre pris en compte. 103 : IF ((DT1.GE.0.D0).AND.(DT1.LE.((1+EPSILO)*DTREEL))) THEN 104 : DTNEW=MIN(DT1,DTREEL) 105 : IF ((DT2.GE.0.D0).AND.(DT2.LE.((1+EPSILO)*DTNEW))) THEN 106 : DTNEW=MIN(DT2,DTREEL) 107 : ENDIF 108 : ELSEIF ((DT2.GE.0.D0).AND.(DT2.LE.((1+EPSILO)*DTREEL))) THEN 109 : DTNEW=MIN(DT2,DTREEL) 110 : ELSE 111 : ITEST=0 112 : RETURN 113 : ENDIF 114 : *** calcul des éventuels pts d'arrivee apres recalibrage 115 : DO 20 I=1,NDIM 116 : XXCO2(I)=XDEP2(I)+DTNEW*(UREEL(I)-IREBCO*COEF2B*XNREB(I)) 117 : 20 CONTINUE 118 : DO 30 I=1,NDIM 119 : X(I)=U(I)*SQRT(DTNEW)-XIREB(I) 120 : 30 CONTINUE 121 : COEFD=SCVECT(X,XNREB,NDIM) 122 : DO 40 I=1,NDIM 123 : XXAR2(I)=XXCO2(I)+U(I)*SQRT(DTNEW)-2*IREBDI*COEFD*XNREB(I) 124 : 40 CONTINUE 125 : *** verification non invariance du pt recalibre 126 : IF (LTEST.EQ.0) THEN 127 : XNORM2=0.D0 128 : DO 50 I=1,NDIM 129 : XNORM2=XNORM2+(XXAR2(I)-XARI2(I))**2 130 : 50 CONTINUE 131 : IF (XNORM2.LT.(XZPREC*DCARA2)) THEN 132 : ITEST=0 133 : RETURN 134 : ENDIF 135 : ENDIF 136 : *** recuperation des nouveaux pts d'arrivee efffectifs 137 : DO 60 I=1,NDIM 138 : XCON2(I)=XXCO2(I) 139 : XARI2(I)=XXAR2(I) 140 : TARI=TDEP+DTNEW 141 : 60 CONTINUE 142 : 143 : RETURN 144 : END 145 : 146 :
© Cast3M 2003 - All rights reserved.
Disclaimer