Télécharger la plaquette Cast3M


fr en

Télécharger solve.eso

Retour à la liste

Numérotation des lignes : 

solve

  1. C SOLVE     SOURCE    CB215821  26/05/04    21:15:05     12531          
  2. c*****************************************************************************
  3. c       sous rouine utilisee dans l'ecoulement plastique de la loi
  4. c       de GURSON
  5. c**********************************************************************
  6.         subroutine solve(dphik,iter,wrkgur)
  7. c       this subroutine solve the non linear equation which solution is
  8. c       the phi increment
  9. c
  10. c---    variables
  11.         IMPLICIT INTEGER(I-N)
  12.         IMPLICIT REAL*8 (A-H,O-Z)
  13.  
  14.  
  15. C       Une SUBROUTINE nommee "split" est entree dans la norme FORTRAN 2023
  16. C       GCC ne veut plus compiler cette source car le nombre d'argument de SPLIT
  17. C       ne correspond pas au modele INTRINSIC ==> Ajout de EXTERNAL SPLIT
  18.         EXTERNAL SPLIT
  19.  
  20.         segment wrkgur
  21.          real*8 sigbar, sy0,phi0,rho0,g,b,h
  22.          real*8 epn,phin,sqrtj2,rho,sig(6)
  23.          real*8 e(7),dt
  24.          real*8 conv,tol1,tol2
  25.         endsegment
  26. *        common/prop/sigbar,sy0,phi0,rho0,g,b,h
  27. *        common/state/epn,phin,sqrtj2,rho,sig(6)
  28. *       common /toler/ conv,tol1,tol2
  29. c
  30. c***    initial data
  31.         p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
  32.         phimax=1.d0-rho/rho0*(1.d0-phi0)
  33.         if ( phimax .lt. 0.d0) phimax=0.d0
  34. c
  35. c***    find the first dphi for which syp>0
  36.         call ziro(dphik,wrkgur)
  37.         phik=phin+dphik
  38. c
  39.         if ( phik .lt. phi0) then
  40.           call split(dphik,wrkgur)
  41.           return
  42.         endif
  43. c
  44. c---    initialize the search
  45.         call func(fk,dphik,wrkgur)
  46.         iter=1
  47. c
  48.         if ( abs(fk) .le. conv) return
  49. c
  50.         if ( fk .gt. 0.d0) then
  51.            return
  52.         endif
  53. c
  54. c*******perform iterations until find fk>0
  55.         if ( phimax .gt. 0.d0) then
  56. c
  57. c---    compute the constant used to determine dphi
  58.         c=sqrt(abs(phimax-phin)/phimax/h/b)*2.d0*sigbar/3.d0*
  59.      &        (h*(1-phimax)+3.*g)
  60.         c=c/sy0/2.d0
  61. c---    iterates until converges
  62.         do 100 i=1,50
  63. c
  64. c       check for convergence
  65.           if (abs(fk).le. conv) return
  66. c
  67. c       check if fk>0 and then exit loop
  68.           if (fk .gt. 0.d0) goto 200
  69. c
  70. c       check if split is needed
  71.           if ( phik .lt. phi0) then
  72.             call split(dphik,wrkgur)
  73.             return
  74.           endif
  75. c
  76. c---    keep track of the previous step
  77.                 fj=fk
  78.                 dphij=dphik
  79. c
  80.                 iter=iter+1
  81. c
  82. c---    find new value of dphik using the interpolation function
  83.                 d=fk-c*sqrt(1.d0/abs(phimax-phik))
  84.                 phik=phimax+sign(1.d0,p)*(c/d)**2
  85.                 dphik=phik-phin
  86.                 call func(fk,dphik,wrkgur)
  87. c
  88.   100   continue
  89. c
  90. c*******case when phimax is negative
  91.         else
  92. c
  93. c---    initialize data for interpolation function
  94.         p0=b*(1-rho0/rho/(1-phi0))
  95.         c2=(h+3.d0*g)*sqrt(phin*sigbar/1.5d0/h
  96.      &      /sinh(1.5d0*p0/sigbar))
  97.      &      /1.d1
  98. c
  99. c---    iterates til converges
  100.         do 150 i=1,50
  101. c
  102. c       check for convergence
  103.           if (abs(fk).le. conv) return
  104. c
  105. c       check if fk>0 and then exit loop
  106.           if (fk .gt. 0.d0) goto 200
  107. c
  108. c       check if split is needed
  109.           if ( phik .lt. phi0) then
  110.             call split(dphik,wrkgur)
  111.             return
  112.           endif
  113. c
  114. c---    keep track of the previous step
  115.           fj=fk
  116.           dphij=dphik
  117.           iter=iter+1
  118. c
  119. c---    find new value of dphik using the interpolation function
  120.           d=fk-c2*sqrt(1.d0/phik)
  121.           phik=sign(1.d0,p)*(c2/d)**2
  122.           dphik=phik-phin
  123.           call func(fk,dphik,wrkgur)
  124. c
  125.   150   continue
  126.         endif
  127. c
  128. c---    the previous algorithm didn't converge
  129.         return
  130.   200   continue
  131. c
  132. c*******find zero using Ridders method
  133. c       see Numerical recipee in Fortran p351
  134.         do 300 i=1,200
  135.            iter=iter+1
  136. c
  137.            half=(dphik+dphij)/2
  138.            call func(fm,half,wrkgur)
  139.            s=sqrt(abs(fm**2-fk*fj))
  140.            dum=half+(half-dphij)*(sign(1.d0,fj-fk)*fm/s)
  141. c
  142. c---    new interval
  143.            call func(fdum,dum,wrkgur)
  144.            if (fdum .lt. 0.d0) then
  145.               dphij=dum
  146.               fj=fdum
  147.            else
  148.               dphik=dum
  149.               fk=fdum
  150.            endif
  151. c
  152.            if ( phin+dphij .lt. phi0) then
  153.              call split(dphik,wrkgur)
  154.              return
  155.            endif
  156. c
  157. c---    check for convergence
  158.         if (abs(fdum) .lt. conv) then
  159.               dphik=dum
  160.               if ( (phin+dphik) .le. phi0) then
  161.                  call split(dphik,wrkgur)
  162.                  return
  163.               endif
  164.               return
  165.         endif
  166. c
  167.   300   continue
  168.         return
  169. c
  170.         end
  171. c*************************************************************
  172.  
  173.  
  174.  
  175.  
  176.  
  177.  
  178.  
  179.  
  180.  
  181.  
  182.  

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