/[escript]/trunk/finley/test/python/damage.py
ViewVC logotype

Diff of /trunk/finley/test/python/damage.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3192 by gross, Mon Sep 13 10:30:24 2010 UTC revision 3193 by gross, Tue Sep 21 06:56:44 2010 UTC
# Line 49  RHO = 2800*U.kg/U.m**3 Line 49  RHO = 2800*U.kg/U.m**3
49  G = 10*U.m/U.sec**2     *0  G = 10*U.m/U.sec**2     *0
50  DT=1.*U.sec  DT=1.*U.sec
51  VMAX=-0.01*U.m/U.sec  VMAX=-0.01*U.m/U.sec
52    DT_MAX=0.1*U.sec
53    
54  T_END=100000*U.yr                       # end time  T_END=100000*U.yr                       # end time
55    
# Line 140  dt_old=None Line 141  dt_old=None
141  while t<T_END:  while t<T_END:
142    
143      print "start time step %d"%(n+1,)      print "start time step %d"%(n+1,)
     if n>1:  
        print" eps :"  
        print sqI2**2  
        print 2*mu*(sqI2-gamma/(2*mu))*eps_e[1,1]-gamma*(sqI2-lame/gamma*I1)*sqI2  
144    
145      I1=trace(eps_e)      I1=trace(eps_e)
146      sqI2=length(eps_e)      sqI2=length(eps_e)
147      xi=safeDiv(I1,sqI2)      xi=safeDiv(I1,sqI2)
148      i_xi=safeDiv(sqI2,I1)      i_xi=safeDiv(sqI2,I1)
149      print "\txi = [ %e, %e]"%(inf(xi),sup(xi))      print "\txi = [ %e, %e]"%(inf(xi),sup(xi))
150    
151        if n>1:
152           print" error:",Lsup((2*mu*sqI2-gamma*I1)*eps_e[1,1]-(gamma*sqI2-lame*I1)*sqI2)/Lsup(length(sigma)*sqI2)
153           print" error:",Lsup(sigma[1,1])/Lsup(length(sigma))
154            
155      # update damage parameter:      # update damage parameter:
156      m=wherePositive(xi-XI_0)      m=wherePositive(xi-XI_0)
# Line 213  while t<T_END: Line 214  while t<T_END:
214            error=fac*0.5*dt**2            error=fac*0.5*dt**2
215            print "\testimated local error for time step size %e is %e"%(dt,error)            print "\testimated local error for time step size %e is %e"%(dt,error)
216            dt_new=sqrt(2./3.*ODE_TOL*2/fac)            dt_new=sqrt(2./3.*ODE_TOL*2/fac)
217            dt_new=min(max(dt_new,dt/5),dt*5) # avid rapit changes            dt_new=min(max(dt_new,dt/5),dt*5,DT_MAX) # avid rapit changes
218            print "\tINFO: new time step size %e"%dt_new            print "\tINFO: new time step size %e"%dt_new
219      dt, dt_old=dt_new, dt      dt, dt_old=dt_new, dt

Legend:
Removed from v.3192  
changed lines
  Added in v.3193

  ViewVC Help
Powered by ViewVC 1.1.26