Télécharger calidt.eso

Retour à la liste

Numérotation des lignes :

  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. *** Rq : 'XZPREC' (-INC CCREEL) erreur precision calcul machine
  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.  
  147.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales