/[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 2414 by jfenwick, Mon Mar 30 02:13:58 2009 UTC revision 2415 by gross, Wed May 13 02:48:39 2009 UTC
# Line 33  __author__="Lutz Gross, l.gross@uq.edu.a Line 33  __author__="Lutz Gross, l.gross@uq.edu.a
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE  from flows import StokesProblemCartesian
37  from pdetools import Defect, NewtonGMRES, ArithmeticTuple  from pdetools import MaxIterReached
38    
39  class PowerLaw(object):  class PowerLaw(object):
40      """      """
# Line 280  class PowerLaw(object): Line 280  class PowerLaw(object):
280           iter =0           iter =0
281           converged=False           converged=False
282           tau=eta_eff*gamma_dot           tau=eta_eff*gamma_dot
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 "PowerLaw: Start calculation of eta_eff (tolerance = %s)\nPowerLaw: initial 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)
# Line 298  class PowerLaw(object): Line 298  class PowerLaw(object):
298               d=util.Lsup(eta_eff-eta_eff_old)               d=util.Lsup(eta_eff-eta_eff_old)
299               l=util.Lsup(eta_eff)               l=util.Lsup(eta_eff)
300               iter+=1               iter+=1
301               if self.__verbose: print "step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))               if self.__verbose: print "PowerLaw: step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
302               converged= d<= rtol* l               converged= d<= rtol* l
303             if self.__verbose: print "PowerLaw: Start calculation of eta_eff finalized after %s steps."%iter
304           return eta_eff           return eta_eff
305    
306  #====================================================================================================================================  #====================================================================================================================================
307    class Rheology(object):
 class IncompressibleIsotropicKelvinFlow(Defect):  
308        """        """
309        This class implements the rheology of an isotropic Kelvin material.        General framework to implement a rheology
   
       Typical usage::  
   
           sp = IncompressibleIsotropicKelvinFlow(domain, stress=0, v=0)  
           sp.setTolerance()  
           sp.initialize(...)  
           v,p = sp.solve(v0, p0)  
   
       @note: This model has been used in the self-consistent plate-mantle model  
              proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}  
              and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:  
              I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},  
              see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract?CRETRY=1&SRETRY=0>}.  
   
310        """        """
311        def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, useJaumannStress=True, **kwargs):        def __init__(self, domain, stress=None, v=None, p=None, t=0, verbose=True):
312           """           """
313           Initializes the model.           Initializes the rheology
314    
315           @param domain: problem domain           @param domain: problem domain
316           @type domain: L{domain}           @type domain: L{Domain}
317           @param stress: initial deviatoric stress           @param stress: initial (deviatoric) stress
318             @type stress: a tensor value/field of order 2
319           @param v: initial velocity field           @param v: initial velocity field
320             @type stress: a vector value/field
321           @param p: initial pressure           @param p: initial pressure
322             @type p: a scalar value/field
323           @param t: initial time           @param t: initial time
324           @param useJaumannStress: C{True} if Jaumann stress is used           @type t: C{float}
                                   (not supported yet)  
325           """           """
          super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)  
326           self.__domain=domain           self.__domain=domain
327           self.__t=t           self.__t=t
328           self.__vol=util.integrate(1.,Function(self.__domain))           self.__verbose=verbose
          self.__useJaumannStress=useJaumannStress  
          self.__v_boundary=Vector(0,Solution(self.__domain))  
329           #=======================           #=======================
330             #
331           # state variables:           # state variables:
332           #           #
333           if isinstance(stress,Data):           if stress == None: stress=Tensor(0.,Function(self.__domain))
334              self.__stress=util.interpolate(stress,Function(domain))           if v == None: v=Vector(0.,Solution(self.__domain))
335           else:           if p == None: p=Vector(0.,ReducedSolution(self.__domain))
336              self.__stress=Data(stress,(domain.getDim(),domain.getDim()),Function(domain))           self.setDeviatoricStress(stress)
337           self.__stress-=util.trace(self.__stress)*(util.kronecker(domain)/domain.getDim())           self.setVelocity(v)
338           if isinstance(v,Data):           self.setPressure(p)
339              self.__v=util.interpolate(v,Solution(domain))           self.setDeviatoricStrain()
340           else:           self.setTime(t)
341              self.__v=Data(v,(domain.getDim(),),Solution(domain))           #=============================================================
342           if isinstance(p,Data):           self.setExternals(F=Data(), f=Data(), fixed_v_mask=Data(), v_boundary=Data())
343              self.__p=util.interpolate(p,ReducedSolution(domain))          
344           else:        def getDomain(self):
345              self.__p=Data(p,(),ReducedSolution(domain))            """
346           #=======================            returns the domain.
          # PDE related stuff  
          self.__pde_v=LinearPDE(domain,numEquations=self.getDomain().getDim(),numSolutions=self.getDomain().getDim())  
          self.__pde_v.setSymmetryOn()  
          self.__pde_v.setSolverMethod(preconditioner=LinearPDE.RILU)  
   
          self.__pde_p=LinearPDE(domain)  
          self.__pde_p.setReducedOrderOn()  
          self.__pde_p.setSymmetryOn()  
347    
348           self.setTolerance()            @return: the domain
349           self.setSmall()            @rtype: L{Domain}
350              """
351              return self.__domain
352    
353          def getTime(self):
354              """
355              Returns current time.
356    
357              @return: current time
358              @rtype: C{float}
359              """
360              return self.__t  
361    
362          def setExternals(self, F=None, f=None, fixed_v_mask=None, v_boundary=None):
363              """
364              sets external forces and velocity constraints
365    
366              @param F: external force
367              @type F: vector value/field
368              @param f: surface force
369              @type f: vector value/field on boundary
370              @param fixed_v_mask: location of constraints maked by positive values
371              @type fixed_v_mask: vector value/field
372              @param v_boundary: value of velocity at location of constraints
373              @type v_boundary: vector value/field
374              @note: Only changing parameters need to be specified.
375              """
376              if F != None: self.__F=F
377              if f != None: self.__f=f
378              if fixed_v_mask != None: self.__fixed_v_mask=fixed_v_mask
379              if v_boundary != None: self.__v_boundary=v_boundary
380              
381          def getForce(self):
382              """
383              Returns the external force
384    
385              @return:  external force
386              @rtype: L{Data}
387              """
388              return self.__F
389    
390          def getSurfaceForce(self):
391              """
392              Returns the surface force
393    
394              @return:  surface force
395              @rtype: L{Data}
396              """
397              return self.__f
398    
399          def getVelocityConstraint(self):
400              """
401              Returns the constraint for the velocity as a pair of the
402              mask of the location of the constraint and the values.
403    
404              @return: the locations of fixed velocity and value of velocities at these locations
405              @rtype: C{tuple} of L{Data}s
406              """
407              return self.__fixed_v_mask, self.__v_boundary      
408    
409        def useJaumannStress(self):        def checkVerbose(self):
410            """            """
411            Returns C{True} if Jaumann stress is included.            Returns True if verbose is switched on
412    
413              @return: value of verbosity flag
414              @rtype: C{bool}
415            """            """
416            return self.__useJaumannStress            return self.__verbose
417    
418        def setSmall(self,small=util.sqrt(util.EPSILON)):        def setTime(self,t=0.):
419              """
420              Updates current time.
421    
422              @param t: new time mark
423              @type t: C{float}
424              """
425              self.__t=t
426          #=======================================================================================
427          def getStress(self):
428            """            """
429            Sets a small value to be used.            Returns current stress.
430    
431            @param small: positive small value            @return: current stress
432              @rtype: L{Data} of rank 2
433            """            """
434            self.__small=small            s=self.getDeviatoricStress()
435              p=self.getPressure()
436              k=util.kronecker(self.getDomain())
437              return s-p*(k/trace(k))
438                
439          def getDeviatoricStress(self):
440              """
441              Returns current deviatoric stress.
442    
443        def getSmall(self):            @return: current deviatoric stress
444              @rtype: L{Data} of rank 2
445            """            """
446            Returns small value.            return self.__stress
447            @rtype: positive float  
448          def setDeviatoricStress(self, stress):
449            """            """
450            return self.__small            Sets the current deviatoric stress
451    
452              @param stress: new deviatoric stress
453              @type stress: L{Data} of rank 2
454              """
455              dom=self.getDomain()
456              s=util.interpolate(stress,Function(dom))
457              self.__stress=util.deviatoric(s)
458    
459          def getPressure(self):
460              """
461              Returns current pressure.
462    
463              @return: current stress
464              @rtype: scalar L{Data}
465              """
466              return self.__p
467    
468          def setPressure(self, p):
469              """
470              Sets current pressure.
471              @param p: new deviatoric stress
472              @type p: scalar L{Data}
473              """
474              self.__p=util.interpolate(p,ReducedSolution(self.getDomain()))
475    
476          def getVelocity(self):
477              """
478              Returns current velocity.
479    
480              @return: current velocity
481              @rtype: vector L{Data}
482              """
483              return self.__v
484    
485          def setVelocity(self, v):
486              """
487              Sets current velocity.
488    
489              @param v: new current velocity
490              @type v: vector L{Data}
491              """
492              self.__v=util.interpolate(v,Solution(self.getDomain()))
493    
494          def setDeviatoricStrain(self, D=None):
495              """
496              set deviatoric strain
497    
498              @param D: new deviatoric strain. If D is not present the current velocity is used.
499              @type D: L{Data} of rank 2
500              """
501              if D==None: D=util.deviatoric(util.symmetric(util.grad(2.*self.getVelocity())))
502              self.__D=util.deviatoric(util.interpolate(D,Function(self.getDomain())))
503    
504          def getDeviatoricStrain(self):
505              """
506              Returns deviatoric strain of current velocity.
507    
508              @return: deviatoric strain
509              @rtype: L{Data}  of rank 2
510              """
511              return self.__D
512    
513          def getTau(self):
514              """
515              Returns current second invariant of deviatoric stress
516    
517              @return: second invariant of deviatoric stress
518              @rtype: scalar L{Data}
519              """
520              s=self.getDeviatoricStress()
521              return util.sqrt(0.5)*util.length(s)
522    
523          def getGammaDot(self):
524              """
525              Returns current second invariant of deviatoric strain
526    
527              @return: second invariant of deviatoric strain
528              @rtype: scalar L{Data}
529              """
530              s=self.getDeviatoricStrain()
531              return util.sqrt(2)*util.length(s)
532    
533        def setTolerance(self,tol=1.e-4):        def setTolerance(self,tol=1.e-4):
534            """            """
535            Sets the tolerance.            Sets the tolerance used to terminate the iteration on a time step.
536              See the implementation of the rheology for details.
537    
538              @param tol: relative tolerance to terminate iteration on time step.
539              @type tol: positive C{float}
540            """            """
541            self.__pde_v.setTolerance(tol**2)            if tol<=0.:
542            self.__pde_p.setTolerance(tol**2)                raise ValueError,"tolerance must be non-negative."
543            self.__tol=tol            self.__tol=tol
544    
545        def getTolerance(self):        def getTolerance(self):
546            """            """
547            Returns the set tolerance.            Returns the set tolerance for terminate the iteration on a time step.
548            @rtype: positive float  
549              @rtype: positive C{float}
550            """            """
551            return self.__tol            return self.__tol
552    
553        def getDomain(self):        #=======================================================================
554          def setFlowTolerance(self, tol=1.e-4):
555            """            """
556            Returns the domain.            Sets the relative tolerance for the flow solver
557    
558              @param tol: desired relative tolerance for the flow solver
559              @type tol: positive C{float}
560              @note: Typically this method is overwritten by a subclass.
561            """            """
562            return self.__domain            pass
563          def getFlowTolerance(self):
564              """
565              Returns the relative tolerance for the flow solver
566    
567        def getTime(self):            @return: tolerance of the flow solver
568              @rtype: C{float}
569              @note: Typically this method is overwritten by a subclass.
570            """            """
571            Returns current time.            pass
572          def setFlowSubTolerance(self, tol=1.e-8):
573            """            """
574            return self.__t            Sets the relative tolerance for the subsolver of the flow solver
575    
576        def setExternals(self, F=None, f=None, q=None, v_boundary=None):            @param tol: desired relative tolerance for the subsolver
577              @type tol: positive C{float}
578              @note: Typically this method is overwritten by a subclass.
579              """
580              pass
581          def getFlowSubTolerance(self):
582            """            """
583            Sets externals.            Returns the relative tolerance for the subsolver of the flow solver
584    
585            @param F: external force            @return: tolerance of the flow subsolver
586            @param f: surface force            @rtype: C{float}
587            @param q: location of constraints            @note: Typically this method is overwritten by a subclass.
588            """            """
589            if F != None: self.__pde_v.setValue(Y=F)            pass
           if f != None: self.__pde_v.setValue(y=f)  
           if q != None: self.__pde_v.setValue(q=q)  
           if v_boundary != None: self.__v_boundary=v_boundary  
   
       def bilinearform(self,arg0,arg1):  
         s0=util.deviatoric(util.symmetric(util.grad(arg0[0])))  
         s1=util.deviatoric(util.symmetric(util.grad(arg1[0])))  
         # s0=util.interpolate(arg0[0],Function(self.getDomain()))  
         # s1=util.interpolate(arg1[0],Function(self.getDomain()))  
         p0=util.interpolate(arg0[1],Function(self.getDomain()))  
         p1=util.interpolate(arg1[1],Function(self.getDomain()))  
         a=util.integrate(self.__p_weight**2*util.inner(s0,s1))+util.integrate(p0*p1)  
         return a  
   
       def getEtaEff(self,strain, pressure):  
           if self.__mu==None:  
                 eps=util.length(strain)*util.sqrt(2)  
           else:  
                 eps=util.length(strain+self.__stress/((2*self.__dt)*self.__mu))*util.sqrt(2)  
           p=util.interpolate(pressure,eps.getFunctionSpace())  
           if self.__tau_Y!= None:  
              tmp=self.__tau_Y+self.__friction*p  
              m=util.wherePositive(eps)*util.wherePositive(tmp)  
              eta_max=m*tmp/(eps+(1-m)*util.EPSILON)+(1-m)*util.DBLE_MAX  
           else:  
              eta_max=util.DBLE_MAX  
           # initial guess:  
           tau=util.length(self.__stress)/util.sqrt(2)  
           # start the iteration:  
           cc=0  
           TOL=1e-7  
           dtau=util.DBLE_MAX  
           print "tau = ", tau, "eps =",eps  
           while cc<10 and dtau>TOL*util.Lsup(tau):  
              eta_eff2,eta_eff_dash=self.evalEtaEff(tau,return_dash=True)  
              eta_eff=util.clip(eta_eff2-eta_eff_dash*tau/(1-eta_eff_dash*eps),maxval=eta_max)  
              tau, tau2=eta_eff*eps, tau  
              dtau=util.Lsup(tau2-tau)  
              print "step ",cc,dtau, util.Lsup(tau)  
              cc+=1  
           return eta_eff  
590    
       def getEtaCharacteristic(self):  
           a=0  
           for i in xrange(self.__numMaterials):  
             a=a+1./self.__eta_N[i]  
           return 1/a  
591    
592        def evalEtaEff(self, tau, return_dash=False):  #====================================================================================================================================
          a=Scalar(0,tau.getFunctionSpace())  # =1/eta  
          if return_dash: a_dash=Scalar(0,tau.getFunctionSpace()) # =(1/eta)'  
          s=util.Lsup(tau)  
          if s>0:  
             m=util.wherePositive(tau)  
             tau2=s*util.EPSILON*(1.-m)+m*tau  
             for i in xrange(self.__numMaterials):  
                  eta_N=self.__eta_N[i]  
                  tau_t=self.__tau_t[i]  
                  if tau_t==None:  
                     a+=1./eta_N  
                  else:  
                     power=1.-1./self.__power[i]  
                     c=1./(tau_t**power*eta_N)  
                     a+=c*tau2**power  
                     if return_dash: a_dash+=power*c*tau2**(power-1.)  
          else:  
             for i in xrange(self.__numMaterials):  
                  eta_N=self.__eta_N[i]  
                  power=1.-1./self.__power[i]  
                  a+=util.whereZero(power)/eta_N  
          if self.__mu!=None: a+=1./(self.__dt*self.__mu)  
          out=1/a  
          if return_dash:  
              return out,-out**2*a_dash  
          else:  
              return out  
593    
594        def eval(self,arg):  class IncompressibleIsotropicFlowCartesian(PowerLaw,Rheology):
595           v=arg[0]        """
596           p=arg[1]        This class implements the rheology of an isotropic Kelvin material.
          D=self.getDeviatoricStrain(v)  
          eta_eff=self.getEtaEff(D,p)  
          print "eta_eff=",eta_eff  
          # solve for dv  
          self.__pde_v.setValue(A=Data()) # save memory!  
          k3=util.kronecker(Function(self.getDomain()))  
          k3Xk3=util.outer(k3,k3)  
          self.__pde_v.setValue(A=eta_eff*(util.swap_axes(k3Xk3,0,3)+util.swap_axes(k3Xk3,1,3)),X=-eta_eff*D+p*util.kronecker(self.getDomain()))  
          dv=self.__pde_v.getSolution(verbose=self.__verbose)  
          print "resistep dv =",dv  
          # solve for dp  
          v2=v+dv  
          self.__pde_p.setValue(D=1/eta_eff,Y=util.div(v2))  
          dp=self.__pde_p.getSolution(verbose=self.__verbose)  
          print "resistep dp =",dp  
          return ArithmeticTuple(dv,dp)  
597    
598        def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False):        Typical usage::
599    
600              sp = IncompressibleIsotropicFlow(domain, stress=0, v=0)
601              sp.setTolerance()
602              sp.initialize(...)
603              v,p = sp.solve()
604    
605          @note: This model has been used in the self-consistent plate-mantle model
606                 proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}
607                 and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
608                 I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
609                 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract?CRETRY=1&SRETRY=0>}.
610    
611          """
612          def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, verbose=True):
613             """
614             Initializes the model.
615    
616             @param domain: problem domain
617             @type domain: L{Domain}
618             @param stress: initial (deviatoric) stress
619             @type stress: a tensor value/field of order 2
620             @param v: initial velocity field
621             @type stress: a vector value/field
622             @param p: initial pressure
623             @type p: a scalar value/field
624             @param t: initial time
625             @type t: C{float}
626             @param numMaterials: number of materials
627             @type numMaterials: C{int}
628             @param verbose: if C{True} some informations are printed.
629             @type verbose: C{bool}        
630             """
631             PowerLaw. __init__(self, numMaterials,verbose)
632             Rheology. __init__(self, domain, stress, v, p, t, verbose)
633             self.__solver=StokesProblemCartesian(self.getDomain(),verbose=verbose)
634             self.__eta_eff=None
635             self.setTolerance()
636             self.setFlowTolerance()
637             self.setFlowSubTolerance()
638    
639          def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False, usePCG=True):
640            """            """
641            Updates stress, velocity and pressure for time increment dt.            Updates stress, velocity and pressure for time increment dt.
642    
# Line 533  class IncompressibleIsotropicKelvinFlow( Line 645  class IncompressibleIsotropicKelvinFlow(
645                                   incompressible solver                                   incompressible solver
646            @param verbose: prints some infos in the incompressible solver            @param verbose: prints some infos in the incompressible solver
647            """            """
648            self.__verbose=verbose            if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: start iteration for t = %s."%(self.getTime()+dt,)
649            self.__dt=dt            v_last=self.getVelocity()
650            tol=self.getTolerance()            s_last=self.getDeviatoricStress()
651            # set the initial velocity:            F=self.getForce()
652            m=util.wherePositive(self.__pde_v.getCoefficient("q"))            f=self.getSurfaceForce()
653            v_new=self.__v*(1-m)+self.__v_boundary*m            mask_v,v_b=self.getVelocityConstraint()
654            # and off we go:            mu=self.getElasticShearModulus()
655            x=ArithmeticTuple(v_new, self.__p)            #=========================================================================
656            # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic(),Function(self.__domain))**2            #
657            self.__p_weight=self.getEtaCharacteristic()            #   we use velocity and pressure from the last time step as initial guess:
658            # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic()**2,self.__p.getFunctionSpace())            #
659            atol=self.norm(x)*self.__tol            v=v_last
660            x_new=NewtonGMRES(self, x, iter_max=iter_max,sub_iter_max=inner_iter_max, atol=atol,rtol=0., verbose=verbose)            p=self.getPressure()
661            self.__v=x_new[0]            #
662            self.__p=x_new[1]            #  calculate eta_eff  if we don't have one or elasticity is present.
663            1/0            #
664            # self.__stress=self.getUpdatedStress(...)            if self.__eta_eff == None or  mu!=None:
665            self.__t+=dt               D=self.__getDeviatoricStrain(v)
666            return self.__v, self.__p               if mu==None:
667                     gamma=util.sqrt(2.)*util.length(D)
668        #========================================================================               else:
669                     gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
670        def getNewDeviatoricStress(self,D,eta_eff=None):               if self.__eta_eff == None:
671           if eta_eff==None: eta_eff=self.evalEtaEff(self.__stress,D,self.__p)                   eta0=None
672           s=(2*eta_eff)*D               else:
673           if self.__mu!=None: s+=eta_eff/(self.__dt*self.__mu)*self.__last_stress                    eta0=self.__eta_eff
674           return s               eta_eff=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta0, iter_max=iter_max)
675                 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been initialied."
676              else:
677                 eta_eff = self.__eta_eff
678    
679        def getDeviatoricStress(self):            iter=0
680              converged=False
681              while not converged:
682                 #
683                 #   intialize the solver
684                 #
685                 if mu==None:          
686                    stress0=Data()
687                 else:
688                    stress0=-(eta_eff/(dt*mu))*s_last
689                 self.__solver.initialize(f=F,fixed_u_mask=mask_v,eta=eta_eff,surface_stress=f,stress=stress0)
690                 #
691                 # get a new velcocity and pressure:
692                 #
693                 if mask_v.isEmpty() or v_b.isEmpty():
694                    v0=v
695                 else:
696                    v0=v_b*mask_v+v*(1.-mask_v)
697                 v,p=self.__solver.solve(v0,p,show_details=False,
698                                              verbose=self.checkVerbose(),max_iter=inner_iter_max,usePCG=usePCG)
699                 util.saveVTK("v.xml",v=v)
700                 util.saveVTK("p.xml",p=p)
701                 #
702                 #   update eta_eff:
703                 #
704                 D=self.__getDeviatoricStrain(v)
705                 print D[0,0]
706                 print D[1,1]
707                 print D[1,0]
708                 if mu==None:
709                     gamma=util.sqrt(2.)*util.length(D)
710                 else:
711                     gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))
712                 print gamma
713                 1/0
714                 eta_eff_old ,eta_eff=eta_eff, self.getEtaEff(gamma, pressure=p,dt=dt, eta0=eta_eff, iter_max=iter_max)
715                 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: eta_eff has been updated."
716                 #
717                 # check the change on eta_eff:
718                 #
719                 diff=util.Lsup(eta_eff_old-eta_eff)
720                 n=util.Lsup(eta_eff)
721                 if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: step %s: max. change in eta_eff is %s."%(iter,diff)
722                 converged = diff <= self.getTolerance()* n
723                 iter+=1
724                 if iter >= iter_max:
725                     raise MaxIterReached,"maximum number of iteration steps on time step %e reached."%(self.getTime()+dt)
726              #
727              #   finally we can update the return values:
728              #
729              self.setPressure(p)
730              self.setVelocity(v)
731              self.setDeviatoricStrain(D)
732              if mu==None:          
733                  stress=(2*eta_eff)*D
734              else:
735                  stress=(2.*eta_eff)*(D+s_last/(2*dt*mu))
736              self.setDeviatoricStress(stress)
737              self.__eta_eff = eta_eff
738              self.setTime(self.getTime()+dt)
739              if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed after %s steps."%(self.getTime(),iter)
740              return self.getVelocity(), self.getPressure()
741    
742          def __getDeviatoricStrain(self, v):
743            """            """
744            Returns current stress.            Returns deviatoric strain of velocity v:
745            """            """
746            return self.__stress            return util.deviatoric(util.symmetric(util.grad(2.*v)))
747    
748        def getDeviatoricStrain(self,velocity=None):        def setFlowTolerance(self, tol=1.e-6):
           """  
           Returns strain.  
749            """            """
750            if velocity==None: velocity=self.getVelocity()            Sets the relative tolerance for the flow solver. See L{StokesProblemCartesian.setTolerance} for details.
           return util.deviatoric(util.symmetric(util.grad(velocity)))  
751    
752        def getPressure(self):            @param tol: desired relative tolerance for the flow solver
753              @type tol: positive C{float}
754            """            """
755            Returns current pressure.            self.__solver.setTolerance(tol)
756          def getFlowTolerance(self):
757            """            """
758            return self.__p            Returns the relative tolerance for the flow solver
759    
760        def getVelocity(self):            @return: tolerance of the flow solver
761              @rtype: C{float}
762            """            """
763            Returns current velocity.            return self.__solver.getTolerance()
764          def setFlowSubTolerance(self, tol=1.e-12):
765            """            """
766            return self.__v            Sets the relative tolerance for the subsolver of the flow solver. See L{StokesProblemCartesian.setSubProblemTolerance} for details
767    
768        def getTau(self,stress=None):            @param tol: desired relative tolerance for the subsolver
769              @type tol: positive C{float}
770            """            """
771            Returns current second stress deviatoric invariant.            self.__solver.setSubProblemTolerance(tol)
772          def getFlowSubTolerance(self):
773            """            """
774            if stress==None: stress=self.getDeviatoricStress()            Returns the relative tolerance for the subsolver of the flow solver
           return util.sqrt(0.5)*util.length(stress)  
775    
776        def getGammaDot(self,strain=None):            @return: tolerance of the flow subsolver
777            """            @rtype: C{float}
           Returns current second stress deviatoric invariant.  
778            """            """
779            if strain==None: strain=self.getDeviatoricStrain()            return self.__solver.getSubProblemTolerance()
           return util.sqrt(2)*util.length(strain)  
780    

Legend:
Removed from v.2414  
changed lines
  Added in v.2415

  ViewVC Help
Powered by ViewVC 1.1.26