Sources Esope | Cast3M

calidt
C CALIDT    SOURCE    CHAT      05/01/12    21:46:19     5004
      SUBROUTINE CALIDT(EPSILO,NDIM,XINTER,XNORMA
     $     ,XIREB,XNREB,IREBCO,IREBDI,XDEP2,TDEP,UREEL,U
     $     ,DTREEL,LTEST,DTNEW,XCON2,XARI2,TARI,ITEST)
**********************************************************************
*** SP 'CALIDT' : lorsque la particule sort de la maille, permet de
*** recalibrer le dernier pas de tps de facon à ce que trajectoire de
*** la particule s'arrete à l'interface de la maille.
***
*** APPELES 1 = 'SCVECT' (fonction)
*** APPELES 2 = aucun
***
*** E = 'EPSILO' marge relative acceptée position % element
***     'NDIM' dimension de l'espace
***     'XINTER' pt d'intersection de la face traversee considérée
***     'XNORMA' normale associée à face traversee considérée
***     'XIREB' pt d'impact sur la face impermeable
***     'XNREB' vecteur normal à la face impermeable
***     'IREBCO' vaut 1 si rebond lors de la convection, 0 sinon
***     'IREBDI' vaut 1 si rebond lors de la diffusion, 0 sinon
***     'XDEP2' coordonnees de depart de la particule (avant convec)
***     'TDEP' temps reel de depart de la particule
***     'UREEL' vitesse dans element reel au pt de depart
***     'U' vecteur vitesse "diffusive" = XL*Z = XPASS*XLU*Z
***     'DTREEL' pas de temps reel avant recalibrage
***
*** E/S = 'LTEST' vaut 1 si 1er essai calibrage, 0 si 2e essai
***
*** S = 'DTNEW' pas de tps reel recalibre (le plus petit possible)
***     'XCON2' coord reelles particule apres convection apres calibrage
***     'XARI2' coord reelles particule apres diffusion apres calibrage
***     'TARI' temps reel d'arrivee de la particule
***     'ITEST' vaut 1 si calibrage donne 0 &lt; dt &lt; dtreel, 0 sinon
***
*** Rq : pour pouvoir effectuer un dvlpt limité au 2eme ordre pour le
***      calcul de delta, nous considerons le critere 4ac/b² &lt; epsdl.
***      on a alors une dependance epsdl > Cte*sqrt(EPSILO) due au
***      reste du dl au 2eme ordre.
*** Rq : 'XZPREC' (-INC CCREEL) erreur precision calcul machine
***
*** ORIGINE = CYRIL NOU
**********************************************************************
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCREEL
      DIMENSION XINTER(3),XNORMA(3),XDEP2(3),UREEL(3),U(3)
      DIMENSION XMF(3),XXCO2(3),XXAR2(3),XCON2(3),XARI2(3)
      DIMENSION XIREB(3),XNREB(3),XMFR(3),X(3)
      ITEST=1
***   initialisation des coefficients du polynome 2nd degre à resoudre
      DCARA2=0.D0
      DO 10 I=1,NDIM
         XMF(I)=XINTER(I)-XDEP2(I)
         XMFR(I)=XIREB(I)-XDEP2(I)
         DCARA2=DCARA2+XMF(I)**2
 10   CONTINUE
      PDTNOR=SCVECT(XNREB,XNORMA,NDIM)
      COEF2A=SCVECT(UREEL,XNORMA,NDIM)
      COEF2B=SCVECT(UREEL,XNREB,NDIM)
      COEFF2=COEF2A-PDTNOR*(IREBCO+2*IREBDI)*COEF2B
      COEF1A=SCVECT(U,XNORMA,NDIM)
      COEF1B=SCVECT(U,XNREB,NDIM)
      COEFF1=COEF1A-2*PDTNOR*IREBDI*COEF1B
      COEF0A=SCVECT(XMF,XNORMA,NDIM)
      COEF0B=SCVECT(XMFR,XNREB,NDIM)
      COEFF0=-(COEF0A-2*PDTNOR*IREBDI*COEF0B)
      IF (ABS(4*COEFF2*COEFF0).GE.(100*SQRT(EPSILO)*(COEFF1**2))) THEN
***      cas ou resolution de coeff2*X*X+coeff1*X+coeff0=0
         DELTA=COEFF1*COEFF1-4.D0*COEFF2*COEFF0
         IF (DELTA.GE.(-XZPREC*(DCARA2/DTREEL))) THEN
            RDELTA=SQRT(ABS(DELTA))
            X1=(-COEFF1+RDELTA)/(2.D0*COEFF2)
            X2=(-COEFF1-RDELTA)/(2.D0*COEFF2)
         ELSE
            ITEST=0
            RETURN
         ENDIF
      ELSE
***      cas ou devlpt limité de la quantité "racine de delta"
         ORD1=-COEFF0/COEFF1
         ORD2=-COEFF2*(COEFF0**2)/(COEFF1**3)
         ORD3=-2*(COEFF2**2)*(COEFF0**3)/(COEFF1**5)
         ORD4=-5*(COEFF2**3)*(COEFF0**4)/(COEFF1**7)
         ORD5=-14*(COEFF2**4)*(COEFF0**5)/(COEFF1**9)
         X1=ORD1+ORD2+ORD3+ORD4+ORD5
         IF (ABS(COEFF2).LT.(XZPREC*(SQRT(DCARA2)/DTREEL))) THEN
            X2=X1
         ELSE
            X2=-COEFF1/COEFF2-X1
         ENDIF
      ENDIF
***   recuperation des racines > 0
      DT1=-1.D0
      DT2=-1.D0
      IF (X1.GE.0.D0) DT1=X1**2
      IF (X2.GE.0.D0) DT2=X2**2
***   recuperation du plus gd dt &lt; dtreel
C nous avons considéré un recalibrage possible avec dtnew environ= dttest,
C cas qui peut se produire si particule arrive a proximité de l'interface
C sans recalibrage. si le meme cas de figure se produit avec dtnew environ= 0,
C une marge 'EPSILO' devra etre pris en compte.
      IF ((DT1.GE.0.D0).AND.(DT1.LE.((1+EPSILO)*DTREEL))) THEN
         DTNEW=MIN(DT1,DTREEL)
         IF ((DT2.GE.0.D0).AND.(DT2.LE.((1+EPSILO)*DTNEW))) THEN
            DTNEW=MIN(DT2,DTREEL)
         ENDIF
      ELSEIF ((DT2.GE.0.D0).AND.(DT2.LE.((1+EPSILO)*DTREEL))) THEN
         DTNEW=MIN(DT2,DTREEL)
      ELSE
         ITEST=0
         RETURN
      ENDIF
***   calcul des éventuels pts d'arrivee apres recalibrage
      DO 20 I=1,NDIM
         XXCO2(I)=XDEP2(I)+DTNEW*(UREEL(I)-IREBCO*COEF2B*XNREB(I))
 20   CONTINUE
      DO 30 I=1,NDIM
         X(I)=U(I)*SQRT(DTNEW)-XIREB(I)
 30   CONTINUE
      COEFD=SCVECT(X,XNREB,NDIM)
      DO 40 I=1,NDIM
         XXAR2(I)=XXCO2(I)+U(I)*SQRT(DTNEW)-2*IREBDI*COEFD*XNREB(I)
 40   CONTINUE
***   verification non invariance du pt recalibre
      IF (LTEST.EQ.0) THEN
         XNORM2=0.D0
         DO 50 I=1,NDIM
            XNORM2=XNORM2+(XXAR2(I)-XARI2(I))**2
 50      CONTINUE
         IF (XNORM2.LT.(XZPREC*DCARA2)) THEN
            ITEST=0
            RETURN
         ENDIF
      ENDIF
***   recuperation des nouveaux pts d'arrivee efffectifs
      DO 60 I=1,NDIM
         XCON2(I)=XXCO2(I)
         XARI2(I)=XXAR2(I)
         TARI=TDEP+DTNEW
 60   CONTINUE
 
      RETURN
      END