/[escript]/trunk/escript/py_src/rheologies.py
ViewVC logotype

Diff of /trunk/escript/py_src/rheologies.py

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

revision 2300 by gross, Wed Mar 11 08:17:57 2009 UTC revision 2301 by gross, Thu Mar 12 01:35:26 2009 UTC
# Line 1  Line 1 
   
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2008 by University of Queensland
# Line 67  class PowerLaw(object): Line 66  class PowerLaw(object):
66           self.__tau_t=[1. for i in xrange(self.__numMaterials)]           self.__tau_t=[1. for i in xrange(self.__numMaterials)]
67           self.__power=[1. for i in xrange(self.__numMaterials)]           self.__power=[1. for i in xrange(self.__numMaterials)]
68           self.__tau_Y=None           self.__tau_Y=None
69           self.__friction=0           self.__friction=None
70           self.__mu=None           self.__mu=None
71           self.__verbose=verbose           self.__verbose=verbose
72           self.setEtaTolerance()           self.setEtaTolerance()
# Line 109  class PowerLaw(object): Line 108  class PowerLaw(object):
108           """           """
109           return self.__rtol           return self.__rtol
110      #===========================================================================      #===========================================================================
111      def setDruckerPragerLaw(self,tau_Y=None,friction=0):      def setDruckerPragerLaw(self,tau_Y=None,friction=None):
112            """            """
113            Sets the parameters for the Drucker-Prager model.            Sets the parameters for the Drucker-Prager model.
114    
# Line 225  class PowerLaw(object): Line 224  class PowerLaw(object):
224                 self.setPowerLaw(id=i, eta_N=eta_N[i],tau_t=tau_t[i],power=power[i])                 self.setPowerLaw(id=i, eta_N=eta_N[i],tau_t=tau_t[i],power=power[i])
225    
226      #===========================================================================      #===========================================================================
227      def getEtaEff(self,gamma_dot, eta0=None, pressure=0,dt=None, iter_max=10):      def getEtaEff(self,gamma_dot, eta0=None, pressure=None,dt=None, iter_max=10):
228           """           """
229           returns the effective viscosity eta_eff such that           returns the effective viscosity eta_eff such that
230    
# Line 242  class PowerLaw(object): Line 241  class PowerLaw(object):
241           @type iter_max: C{int}           @type iter_max: C{int}
242           @return: effective viscosity.           @return: effective viscosity.
243           """           """
244             SMALL=1./(util.DBLE_MAX/100.)
245           numMaterial=self.getNumMaterials()           numMaterial=self.getNumMaterials()
246           s=[1.-1./p for p in self.getPower() ]           s=[1.-1./p for p in self.getPower() ]
247           eta_N=self.getEtaN()           eta_N=self.getEtaN()
# Line 265  class PowerLaw(object): Line 265  class PowerLaw(object):
265           if mu !=None and dt == None:           if mu !=None and dt == None:
266               raise ValueError,"Time stepsize dt must be given."               raise ValueError,"Time stepsize dt must be given."
267           if dt !=None:           if dt !=None:
268               if dt<=0: raise ValueEror,"time step size must be positive."               if dt<=0: raise ValueError,"time step size must be positive."
269           if tau_Y==None:           if tau_Y==None and fric==None:
270               eta_max=None               eta_max=None
271           else:           else:
272                if util.inf(pressure)*fric<0:              if fric == None:
273                     raise ValueError,"pressure needs to be non-negative."                  eta_max=tau_Y/(gamma_dot+SMALL*util.whereZero(gamma_dot))
274                eta_max=(tau_Y+fric*pressure)/(gamma_dot+util.whereZero(gamma_dot)/util.DBLE_MAX)              else:
275                    if tau_Y==None: tau_Y==0
276                    if util.inf(fric)<=0:
277                        raise ValueError,"if friction present it needs to be positive."
278                    eta_max=fric*util.clip(tau_Y/fric+pressure,minval=0)/(gamma_dot+SMALL*util.whereZero(gamma_dot))
279           rtol=self.getEtaTolerance()           rtol=self.getEtaTolerance()
280           iter =0           iter =0
281           converged=False           converged=False
# Line 279  class PowerLaw(object): Line 283  class PowerLaw(object):
283           if self.__verbose: print "Start calculation of eta_eff (tolerance = %s)\ninitial max eta_eff = %s, tau = %s."%(rtol,util.Lsup(eta_eff),util.Lsup(tau))           if self.__verbose: print "Start calculation of eta_eff (tolerance = %s)\ninitial max eta_eff = %s, tau = %s."%(rtol,util.Lsup(eta_eff),util.Lsup(tau))
284           while not converged:           while not converged:
285               if iter>max(iter_max,1):               if iter>max(iter_max,1):
286                  raise RunTimeError,"tolerance not reached after %s steps."%max(iter_max,1)                  raise RuntimeError,"tolerance not reached after %s steps."%max(iter_max,1)
287               #===========================================               #===========================================
288               theta=0. # =1/eta               theta=0. # =1/eta
289               omega=0. # = tau*theta'= eta'*tau/eta**2               omega=0. # = tau*theta'= eta'*tau/eta**2
Legend:
Removed from v.2300  
changed lines
  Added in v.2301
  ViewVC Help
Powered by ViewVC 1.1.26