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

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
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
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              """
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              """
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              """
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
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.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()
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
690                 #
691                 # get a new velcocity and pressure:
692                 #
694                    v0=v
695                 else:
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            """            """
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.
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