Esope source files

   1 : C SOLVE     SOURCE    CHAT      05/01/13    03:21:55     5004           
   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 :         segment wrkgur
  14 :          real*8 sigbar, sy0,phi0,rho0,g,b,h
  15 :          real*8 epn,phin,sqrtj2,rho,sig(6)
  16 :          real*8 e(7),dt
  17 :          real*8 conv,tol1,tol2
  18 :         endsegment
  19 : *        common/prop/sigbar,sy0,phi0,rho0,g,b,h
  20 : *        common/state/epn,phin,sqrtj2,rho,sig(6)
  21 : *       common /toler/ conv,tol1,tol2
  22 : c
  23 : c***    initial data
  24 :         p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
  25 :         phimax=1.d0-rho/rho0*(1.d0-phi0)
  26 :         if ( phimax .lt. 0.d0) phimax=0.d0
  27 : c
  28 : c***    find the first dphi for which syp>0
  28 :         call ziro(dphik,wrkgur)
  29 :         phik=phin+dphik
  30 : c
  31 :         if ( phik .lt. phi0) then
  31 :           call split(dphik,wrkgur)
  32 :           return
  33 :         endif
  34 : c
  35 : c---    initialize the search
  35 :         call func(fk,dphik,wrkgur)
  36 :         iter=1
  37 : c
  38 :         if ( abs(fk) .le. conv) return 
  39 : c
  40 :         if ( fk .gt. 0.d0) then
  41 :            return
  42 :         endif
  43 : c
  44 : c*******perform iterations until find fk>0
  45 :         if ( phimax .gt. 0.d0) then
  46 : c
  47 : c---    compute the constant used to determine dphi
  48 :         c=sqrt(abs(phimax-phin)/phimax/h/b)*2.d0*sigbar/3.d0*
  49 :      &        (h*(1-phimax)+3.*g)
  50 :         c=c/sy0/2.d0
  51 : c---    iterates until converges
  52 :         do 100 i=1,50
  53 : c
  54 : c       check for convergence
  55 :           if (abs(fk).le. conv) return
  56 : c
  57 : c       check if fk>0 and then exit loop
  58 :           if (fk .gt. 0.d0) goto 200
  59 : c
  60 : c       check if split is needed
  61 :           if ( phik .lt. phi0) then
  61 :             call split(dphik,wrkgur)
  62 :             return
  63 :           endif
  64 : c
  65 : c---    keep track of the previous step
  66 :                 fj=fk
  67 :                 dphij=dphik
  68 : c
  69 :                 iter=iter+1
  70 : c
  71 : c---    find new value of dphik using the interpolation function
  72 :                 d=fk-c*sqrt(1.d0/abs(phimax-phik))
  73 :                 phik=phimax+sign(1.d0,p)*(c/d)**2
  74 :                 dphik=phik-phin
  74 :                 call func(fk,dphik,wrkgur)
  75 : c
  76 :   100   continue
  77 : c
  78 : c*******case when phimax is negative
  79 :         else
  80 : c
  81 : c---    initialize data for interpolation function
  82 :         p0=b*(1-rho0/rho/(1-phi0))
  83 :         c2=(h+3.d0*g)*sqrt(phin*sigbar/1.5d0/h
  84 :      &      /sinh(1.5d0*p0/sigbar))
  85 :      &      /1.d1
  86 : c
  87 : c---    iterates til converges  
  88 :         do 150 i=1,50
  89 : c
  90 : c       check for convergence
  91 :           if (abs(fk).le. conv) return
  92 : c
  93 : c       check if fk>0 and then exit loop
  94 :           if (fk .gt. 0.d0) goto 200
  95 : c
  96 : c       check if split is needed
  97 :           if ( phik .lt. phi0) then
  97 :             call split(dphik,wrkgur)
  98 :             return
  99 :           endif
 100 : c
 101 : c---    keep track of the previous step
 102 :           fj=fk
 103 :           dphij=dphik
 104 :           iter=iter+1
 105 : c
 106 : c---    find new value of dphik using the interpolation function
 107 :           d=fk-c2*sqrt(1.d0/phik)
 108 :           phik=sign(1.d0,p)*(c2/d)**2
 109 :           dphik=phik-phin
 109 :           call func(fk,dphik,wrkgur)
 110 : c
 111 :   150   continue
 112 :         endif
 113 : c
 114 : c---    the previous algorithm didn't converge
 115 :         return 
 116 :   200   continue 
 117 : c
 118 : c*******find zero using Ridders method
 119 : c       see Numerical recipee in Fortran p351 
 120 :         do 300 i=1,200
 121 :            iter=iter+1
 122 : c
 123 :            half=(dphik+dphij)/2
 123 :            call func(fm,half,wrkgur)
 124 :            s=sqrt(abs(fm**2-fk*fj))
 125 :            dum=half+(half-dphij)*(sign(1.d0,fj-fk)*fm/s)
 126 : c
 127 : c---    new interval
 127 :            call func(fdum,dum,wrkgur)
 128 :            if (fdum .lt. 0.d0) then
 129 :               dphij=dum
 130 :               fj=fdum
 131 :            else
 132 :               dphik=dum
 133 :               fk=fdum
 134 :            endif
 135 : c
 136 :            if ( phin+dphij .lt. phi0) then 
 136 :              call split(dphik,wrkgur) 
 137 :              return 
 138 :            endif 
 139 : c
 140 : c---    check for convergence
 141 :         if (abs(fdum) .lt. conv) then
 142 :               dphik=dum
 143 :               if ( (phin+dphik) .le. phi0) then 
 143 :                  call split(dphik,wrkgur)
 144 :                  return 
 145 :               endif 
 146 :               return
 147 :         endif
 148 : c
 149 :   300   continue
 150 :         return
 151 : c
 152 :         end
 153 : c*************************************************************
 154 : 
 155 : 
 156 : 
 157 : 
 158 : 
 159 : 
 160 : 
 161 :  
 162 :