/[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 2849 by gross, Thu Jan 14 05:51:28 2010 UTC revision 2850 by gross, Fri Jan 15 07:14:25 2010 UTC
# Line 518  class Rheology(object): Line 518  class Rheology(object):
518            self.setVelocity(v)            self.setVelocity(v)
519            self.setPressure(p)            self.setPressure(p)
520            self.setDeviatoricStrain()            self.setDeviatoricStrain()
521              self.setGammaDot()
522            self.setTime(t)            self.setTime(t)
523    
524        def setDeviatoricStrain(self, D=None):        def setDeviatoricStrain(self, D=None):
525            """            """
526            set deviatoric strain            set deviatoric strain
527    
528            :param D: new deviatoric strain. If D is not present the current velocity is used.            :param D: new deviatoric strain. If ``D`` is not present the current velocity is used.
529            :type D: `Data` of rank 2            :type D: `Data` of rank 2
530            """            """
531            if D==None: D=util.deviatoric(util.symmetric(util.grad(2.*self.getVelocity())))            if D==None:
532            self.__D=util.deviatoric(util.interpolate(D,Function(self.getDomain())))                self.__D=self.getDeviatoricStrain(self.getVelocity())
533              else:
534                  self.__D=util.deviatoric(util.interpolate(D,Function(self.getDomain())))
535    
536        def getDeviatoricStrain(self):        def getDeviatoricStrain(self, v=None):
537            """            """
538            Returns deviatoric strain of current velocity.            Returns deviatoric strain of current velocity or if ``v`` is present the
539              deviatoric strain of velocity ``v``:
540    
541            :return: deviatoric strain            :param v: a velocity field
542              :type v: `Data` of rank 1
543              :return: deviatoric strain of the current velocity field or if ``v`` is present the deviatoric strain of velocity ``v``
544            :rtype: `Data`  of rank 2            :rtype: `Data`  of rank 2
545            """            """
546            return self.__D            if v==None:
547                 return self.__D
548              else:
549                 return util.deviatoric(util.symmetric(util.grad(v)))
550    
551        def getTau(self):        def getTau(self):
552            """            """
# Line 549  class Rheology(object): Line 558  class Rheology(object):
558            s=self.getDeviatoricStress()            s=self.getDeviatoricStress()
559            return util.sqrt(0.5)*util.length(s)            return util.sqrt(0.5)*util.length(s)
560    
561        def getGammaDot(self):        def setGammaDot(self, gammadot=None):
562              """
563              set the second invariant of deviatoric strain rate. If ``gammadot`` is not present zero is used.
564    
565              :param gammadot: second invariant of deviatoric strain rate.
566              :type gammadot: `Data` of rank 1
567              """
568              if gammadot == None:
569                   self.__gammadot = Scalar(0.,Function(self.getDomain()))
570              else:
571                   self.__gammadot=gammadot
572              
573          def getGammaDot(self, D=None):
574            """            """
575            Returns current second invariant of deviatoric strain            Returns current second invariant of deviatoric strain rate or if ``D`` is present the second invariant of ``D``.
576    
577              :param D: deviatoric strain rate tensor
578              :type D: `Data`  of rank 0
579            :return: second invariant of deviatoric strain            :return: second invariant of deviatoric strain
580            :rtype: scalar `Data`            :rtype: scalar `Data`
581            """            """
582            s=self.getDeviatoricStrain()            if D==None:
583            return util.sqrt(2)*util.length(s)                return self.__gammadot
584              else:
585                  return util.sqrt(2.)*util.length(D)
586              
587    
588  #====================================================================================================================================  #====================================================================================================================================
589    
# Line 609  class IncompressibleIsotropicFlowCartesi Line 635  class IncompressibleIsotropicFlowCartesi
635            """            """
636            return self.__eta_eff            return self.__eta_eff
637    
      def __getDeviatoricStrain(self, v):  
           """  
           Returns deviatoric strain of velocity v:  
           """  
           return util.deviatoric(util.symmetric(util.grad(v)))  
638    
639       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
640           """           """
# Line 628  class IncompressibleIsotropicFlowCartesi Line 649  class IncompressibleIsotropicFlowCartesi
649           #           #
650           #  calculate eta_eff if we don't have one or elasticity is present.           #  calculate eta_eff if we don't have one or elasticity is present.
651           #           #
          D=self.__getDeviatoricStrain(v)  
652           if mu==None:           if mu==None:
653               gamma=util.sqrt(2.)*util.length(D)               gamma=self.getGammaDot(self.getDeviatoricStrain(v))
654           else:           else:
655               gamma=util.sqrt(2.)*util.length(D+s_last/(2*dt*mu))               gamma=self.getGammaDot(self.getDeviatoricStrain(v)+s_last/(2*dt*mu))
656    
657           self.__eta_eff_save=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=self.__eta_eff_save, iter_max=self.__eta_iter_max)           self.__eta_eff_save=self.getEtaEff(gamma, pressure=p,dt=dt, eta0=self.__eta_eff_save, iter_max=self.__eta_iter_max)
658    
# Line 693  class IncompressibleIsotropicFlowCartesi Line 713  class IncompressibleIsotropicFlowCartesi
713            #            #
714            # get a new velcocity and pressure:            # get a new velcocity and pressure:
715            #            #
716            if mask_v.isEmpty() or v_b.isEmpty():            if mask_v.isEmpty():
717                 v0=v_last                 v0=v_last
718            else:            else:
719                v0=v_b*mask_v+v_last*(1.-mask_v)                if v_b.isEmpty():
720                     v0=v_last*(1.-mask_v)
721                  else:
722                     v0=v_b*mask_v+v_last*(1.-mask_v)
723    
724            v,p=self._solve(v0,p_last,verbose=self.checkVerbose(),max_iter=iter_max,usePCG=usePCG, max_correction_steps=max_correction_steps)            v,p=self._solve(v0,p_last,verbose=self.checkVerbose(),max_iter=iter_max,usePCG=usePCG, max_correction_steps=max_correction_steps)
725            #            #
# Line 704  class IncompressibleIsotropicFlowCartesi Line 727  class IncompressibleIsotropicFlowCartesi
727            #            #
728            self.setPressure(p)            self.setPressure(p)
729            self.setVelocity(v)            self.setVelocity(v)
730            D=self.__getDeviatoricStrain(v)            self.setDeviatoricStrain(self.getDeviatoricStrain(v))
731            self.setDeviatoricStrain(D)            if mu==None:
732            self.__eta_eff = self.getEtaEff(self.getGammaDot(), pressure=p,dt=dt, eta0=self.__eta_eff_save, iter_max=self.__eta_iter_max)               D=self.getDeviatoricStrain(v)
           if mu==None:            
               stress=(2*self.__eta_eff)*D  
733            else:            else:
734                stress=(2.*self.__eta_eff)*(D+s_last/(2*dt*mu))               D=self.getDeviatoricStrain(v)+s_last/(2*dt*mu)
735            self.setDeviatoricStress(stress)            gamma=self.getGammaDot(D)
736              self.setGammaDot(gamma)
737              self.__eta_eff = self.getEtaEff(self.getGammaDot(), pressure=p,dt=dt, eta0=self.__eta_eff_save, iter_max=self.__eta_iter_max)
738              self.setDeviatoricStress(2.*self.__eta_eff*D)
739            self.setTime(self.getTime()+dt)            self.setTime(self.getTime()+dt)
740            if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed."%(self.getTime(),)            if self.checkVerbose(): print "IncompressibleIsotropicFlowCartesian: iteration on time step %s completed."%(self.getTime(),)
741            return self.getVelocity(), self.getPressure()            return self.getVelocity(), self.getPressure()

Legend:
Removed from v.2849  
changed lines
  Added in v.2850

  ViewVC Help
Powered by ViewVC 1.1.26