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

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()
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:
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
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              """
570              else:
572
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:
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:
"""
638
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:
654           else:           else:
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            #            #
717                 v0=v_last                 v0=v_last
718            else:            else:
721                  else:
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)