escript/py_src/rheologies.py


########################################################
#
# Copyright (c) 2003-2008 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2008 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.uq.edu.au/esscc/escript-finley"

"""
Some models for flow

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__author__="Lutz Gross, l.gross@uq.edu.au"

from escript import *
import util
from linearPDEs import LinearPDE
from pdetools import Defect, NewtonGMRES, ArithmeticTuple

class IncompressibleIsotropicKelvinFlow(Defect):
      """
      This class implements the reology of an isotropic kelvin material. 

      @note: this model has been used in the self-consistent plate-mantel 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>}.

      typical usage:

            sp=PlateMantelModel(domain,stress=0,v=0)
            sp.setTolerance()
            sp.initialize(...)
            v,p=sp.solve(v0,p0)
      """
      def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, useJaumannStress=True, **kwargs):
         """
         initializes PlateMantelModel 

         @param domain: problem domain
         @type domain: L{domain}
         @param stress: initial deviatoric stress
         @param v: initial velocity field
         @param p: initial pressure
         @param t: initial time
         @param useJaumannStress: C{True} if Jaumann stress is used (not supported yet)
         """
         if numMaterials<1:
            raise ValueError,"at least one material must be defined."
         super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)
         self.__domain=domain
         self.__t=t
         self.__vol=util.integrate(1.,Function(self.__domain))
         self.__useJaumannStress=useJaumannStress
         self.__numMaterials=numMaterials
         self.__eta_N=[None for i in xrange(self.__numMaterials)]
         self.__tau_t=[None for i in xrange(self.__numMaterials)]
         self.__power=[1 for i in xrange(self.__numMaterials)]
         self.__tau_Y=None
         self.__friction=None
         self.__mu=None
         self.__v_boundary=Vector(0,Solution(self.__domain))
         #=======================
         # state variables:
         #
         if isinstance(stress,Data):
            self.__stress=util.interpolate(stress,Function(domain))
         else:
            self.__stress=Data(stress,(domain.getDim(),domain.getDim()),Function(domain))
         self.__stress-=util.trace(self.__stress)*(util.kronecker(domain)/domain.getDim())
         if isinstance(v,Data):
            self.__v=util.interpolate(v,Solution(domain))
         else:
            self.__v=Data(v,(domain.getDim(),),Solution(domain)) 
         if isinstance(p,Data):
            self.__p=util.interpolate(p,ReducedSolution(domain))
         else:
            self.__p=Data(p,(),ReducedSolution(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()

         self.setTolerance()
         self.setSmall()
      def useJaumannStress(self):
          """ 
          return C{True} if Jaumann stress is included.
          """
          return self.__useJaumannStress

      def setSmall(self,small=util.sqrt(util.EPSILON)):
          """
          sets small value
      
          @param small: positive small value
          """
          self.__small=small

      def getSmall(self):
          """
          returns small value
          @rtype: positive float
          """
          return self.__small

      def setTolerance(self,tol=1.e-4):
          """
          sets the tolerance
          """
          self.__pde_v.setTolerance(tol**2)
          self.__pde_p.setTolerance(tol**2)
          self.__tol=tol

      def getTolerance(self):
          """
          returns the set tolerance
          @rtype: positive float
          """
          return self.__tol

      def getDomain(self):
          """
          returns the domain
          """
          return self.__domain
      def getTime(self):
          """
          returns current time
          """
          return self.__t

      def setPowerLaw(self,id,eta_N, tau_t=None, power=1):
          """
          sets the power-law parameters for material q
          """
          if id<0 or id>=self.__numMaterials:
              raise ValueError,"Illegal material id = %s."%id
          self.__eta_N[id]=eta_N
          self.__power[id]=power
          self.__tau_t[id]=tau_t

      def setPowerLaws(self,eta_N, tau_t, power):
          """
          set the parameters of the powerlaw for all materials
          """
          if len(eta_N)!=self.__numMaterials or len(tau_t)!=self.__numMaterials or len(power)!=self.__numMaterials:
              raise ValueError,"%s materials are expected."%self.__numMaterials
          for i in xrange(self.__numMaterials):
               self.setPowerLaw(i,eta_N[i],tau_t[i],power[i])
      def setDruckerPragerLaw(self,tau_Y=None,friction=0):
          """
          set the parameters for the Drucker-Prager model
          """
          self.__tau_Y=tau_Y
          self.__friction=friction
          
      def setElasticShearModulus(self,mu=None):
          """ 
          sets the elastic shere modulus
          """
          self.__mu=mu
      def setExternals(self, F=None, f=None, q=None, v_boundary=None):
          """ 
          sets 

          @param F: external force.
          @param f: surface force
          @param q: location of constraints
          """
          if F != None: self.__pde_v.setValue(Y=F)
          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

      def getEtaCharacteristic(self):
          a=0
          for i in xrange(self.__numMaterials):
            a=a+1./self.__eta_N[i]
          return 1/a
             
      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
             
      def eval(self,arg):
         v=arg[0]
         p=arg[1]
         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)

      def update(self,dt, iter_max=100, inner_iter_max=20, verbose=False):
          """
          updates stress, velocity and pressure for time increment dt

          @param dt: time increment
          @param max_inner_iter: maximum number of iteration steps in the incompressible solver
          @param verbose: prints some infos in the incompressible solve
          @param show_details: prints some infos while solving PDEs
          @param tol: tolerance for the time step
          """
          self.__verbose=verbose
          self.__dt=dt
          tol=self.getTolerance()
          # set the initial velocity:
          m=util.wherePositive(self.__pde_v.getCoefficient("q"))
          v_new=self.__v*(1-m)+self.__v_boundary*m
          # and off we go:
          x=ArithmeticTuple(v_new, self.__p)
          # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic(),Function(self.__domain))**2
          self.__p_weight=self.getEtaCharacteristic()
          # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic()**2,self.__p.getFunctionSpace())
          atol=self.norm(x)*self.__tol
          x_new=NewtonGMRES(self, x, iter_max=iter_max,sub_iter_max=inner_iter_max, atol=atol,rtol=0., verbose=verbose)
          self.__v=x_new[0]
          self.__p=x_new[1]
          1/0
          # self.__stress=self.getUpdatedStress(...)
          self.__t+=dt
          return self.__v, self.__p

      #=========================================================================================

      def getNewDeviatoricStress(self,D,eta_eff=None):
         if eta_eff==None: eta_eff=self.evalEtaEff(self.__stress,D,self.__p)
         s=(2*eta_eff)*D
         if self.__mu!=None: s+=eta_eff/(self.__dt*self.__mu)*self.__last_stress
         return s

      def getDeviatoricStress(self):
          """
          returns current stress
          """
          return self.__stress
      def getDeviatoricStrain(self,velocity=None):
          """
          return strain
          """
          if velocity==None: velocity=self.getVelocity()
          return util.deviatoric(util.symmetric(util.grad(velocity)))

      def getPressure(self):
          """
          returns current pressure
          """
          return self.__p
      def getVelocity(self):
          """
          returns current velocity
          """
          return self.__v

      def getTau(self,stress=None):
          """
          returns current second stress deviatoric invariant
          """
          if stress==None: stress=self.getDeviatoricStress()
          return util.sqrt(0.5)*util.length(stress)

      def getGammaDot(self,strain=None):
          """
          returns current second stress deviatoric invariant
          """
          if strain==None: strain=self.getDeviatoricStrain()
          return util.sqrt(2)*util.length(strain)

1
2	########################################################
3	#
4	# Copyright (c) 2003-2008 by University of Queensland
5	# Earth Systems Science Computational Center (ESSCC)
6	# http://www.uq.edu.au/esscc
7	#
8	# Primary Business: Queensland, Australia
9	# Licensed under the Open Software License version 3.0
10	# http://www.opensource.org/licenses/osl-3.0.php
11	#
12	########################################################
13
14	__copyright__="""Copyright (c) 2003-2008 by University of Queensland
15	Earth Systems Science Computational Center (ESSCC)
16	http://www.uq.edu.au/esscc
17	Primary Business: Queensland, Australia"""
18	__license__="""Licensed under the Open Software License version 3.0
19	http://www.opensource.org/licenses/osl-3.0.php"""
20	__url__="http://www.uq.edu.au/esscc/escript-finley"
21
22	"""
23	Some models for flow
24
25	@var __author__: name of author
26	@var __copyright__: copyrights
27	@var __license__: licence agreement
28	@var __url__: url entry point on documentation
29	@var __version__: version
30	@var __date__: date of the version
31	"""
32
33	__author__="Lutz Gross, l.gross@uq.edu.au"
34
35	from escript import *
36	import util
37	from linearPDEs import LinearPDE
38	from pdetools import Defect, NewtonGMRES, ArithmeticTuple
39
40	class IncompressibleIsotropicKelvinFlow(Defect):
41	"""
42	This class implements the reology of an isotropic kelvin material.
43
44	@note: this model has been used in the self-consistent plate-mantel model proposed in
45	U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>} and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
46	I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
47	see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract?CRETRY=1&SRETRY=0>}.
48
49	typical usage:
50
51	sp=PlateMantelModel(domain,stress=0,v=0)
52	sp.setTolerance()
53	sp.initialize(...)
54	v,p=sp.solve(v0,p0)
55	"""
56	def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, useJaumannStress=True, **kwargs):
57	"""
58	initializes PlateMantelModel
59
60	@param domain: problem domain
61	@type domain: L{domain}
62	@param stress: initial deviatoric stress
63	@param v: initial velocity field
64	@param p: initial pressure
65	@param t: initial time
66	@param useJaumannStress: C{True} if Jaumann stress is used (not supported yet)
67	"""
68	if numMaterials<1:
69	raise ValueError,"at least one material must be defined."
70	super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)
71	self.__domain=domain
72	self.__t=t
73	self.__vol=util.integrate(1.,Function(self.__domain))
74	self.__useJaumannStress=useJaumannStress
75	self.__numMaterials=numMaterials
76	self.__eta_N=[None for i in xrange(self.__numMaterials)]
77	self.__tau_t=[None for i in xrange(self.__numMaterials)]
78	self.__power=[1 for i in xrange(self.__numMaterials)]
79	self.__tau_Y=None
80	self.__friction=None
81	self.__mu=None
82	self.__v_boundary=Vector(0,Solution(self.__domain))
83	#=======================
84	# state variables:
85	#
86	if isinstance(stress,Data):
87	self.__stress=util.interpolate(stress,Function(domain))
88	else:
89	self.__stress=Data(stress,(domain.getDim(),domain.getDim()),Function(domain))
90	self.__stress-=util.trace(self.__stress)*(util.kronecker(domain)/domain.getDim())
91	if isinstance(v,Data):
92	self.__v=util.interpolate(v,Solution(domain))
93	else:
94	self.__v=Data(v,(domain.getDim(),),Solution(domain))
95	if isinstance(p,Data):
96	self.__p=util.interpolate(p,ReducedSolution(domain))
97	else:
98	self.__p=Data(p,(),ReducedSolution(domain))
99	#=======================
100	# PDE related stuff
101	self.__pde_v=LinearPDE(domain,numEquations=self.getDomain().getDim(),numSolutions=self.getDomain().getDim())
102	self.__pde_v.setSymmetryOn()
103	self.__pde_v.setSolverMethod(preconditioner=LinearPDE.RILU)
104
105	self.__pde_p=LinearPDE(domain)
106	self.__pde_p.setReducedOrderOn()
107	self.__pde_p.setSymmetryOn()
108
109	self.setTolerance()
110	self.setSmall()
111	def useJaumannStress(self):
112	"""
113	return C{True} if Jaumann stress is included.
114	"""
115	return self.__useJaumannStress
116
117	def setSmall(self,small=util.sqrt(util.EPSILON)):
118	"""
119	sets small value
120
121	@param small: positive small value
122	"""
123	self.__small=small
124
125	def getSmall(self):
126	"""
127	returns small value
128	@rtype: positive float
129	"""
130	return self.__small
131
132	def setTolerance(self,tol=1.e-4):
133	"""
134	sets the tolerance
135	"""
136	self.__pde_v.setTolerance(tol**2)
137	self.__pde_p.setTolerance(tol**2)
138	self.__tol=tol
139
140	def getTolerance(self):
141	"""
142	returns the set tolerance
143	@rtype: positive float
144	"""
145	return self.__tol
146
147	def getDomain(self):
148	"""
149	returns the domain
150	"""
151	return self.__domain
152	def getTime(self):
153	"""
154	returns current time
155	"""
156	return self.__t
157
158	def setPowerLaw(self,id,eta_N, tau_t=None, power=1):
159	"""
160	sets the power-law parameters for material q
161	"""
162	if id<0 or id>=self.__numMaterials:
163	raise ValueError,"Illegal material id = %s."%id
164	self.__eta_N[id]=eta_N
165	self.__power[id]=power
166	self.__tau_t[id]=tau_t
167
168	def setPowerLaws(self,eta_N, tau_t, power):
169	"""
170	set the parameters of the powerlaw for all materials
171	"""
172	if len(eta_N)!=self.__numMaterials or len(tau_t)!=self.__numMaterials or len(power)!=self.__numMaterials:
173	raise ValueError,"%s materials are expected."%self.__numMaterials
174	for i in xrange(self.__numMaterials):
175	self.setPowerLaw(i,eta_N[i],tau_t[i],power[i])
176	def setDruckerPragerLaw(self,tau_Y=None,friction=0):
177	"""
178	set the parameters for the Drucker-Prager model
179	"""
180	self.__tau_Y=tau_Y
181	self.__friction=friction
182
183	def setElasticShearModulus(self,mu=None):
184	"""
185	sets the elastic shere modulus
186	"""
187	self.__mu=mu
188	def setExternals(self, F=None, f=None, q=None, v_boundary=None):
189	"""
190	sets
191
192	@param F: external force.
193	@param f: surface force
194	@param q: location of constraints
195	"""
196	if F != None: self.__pde_v.setValue(Y=F)
197	if f != None: self.__pde_v.setValue(y=f)
198	if q != None: self.__pde_v.setValue(q=q)
199	if v_boundary != None: self.__v_boundary=v_boundary
200
201	def bilinearform(self,arg0,arg1):
202	s0=util.deviatoric(util.symmetric(util.grad(arg0[0])))
203	s1=util.deviatoric(util.symmetric(util.grad(arg1[0])))
204	# s0=util.interpolate(arg0[0],Function(self.getDomain()))
205	# s1=util.interpolate(arg1[0],Function(self.getDomain()))
206	p0=util.interpolate(arg0[1],Function(self.getDomain()))
207	p1=util.interpolate(arg1[1],Function(self.getDomain()))
208	a=util.integrate(self.__p_weight*2util.inner(s0,s1))+util.integrate(p0*p1)
209	return a
210
211	def getEtaEff(self,strain, pressure):
212	if self.__mu==None:
213	eps=util.length(strain)*util.sqrt(2)
214	else:
215	eps=util.length(strain+self.__stress/((2self.__dt)self.__mu))*util.sqrt(2)
216	p=util.interpolate(pressure,eps.getFunctionSpace())
217	if self.__tau_Y!= None:
218	tmp=self.__tau_Y+self.__friction*p
219	m=util.wherePositive(eps)*util.wherePositive(tmp)
220	eta_max=mtmp/(eps+(1-m)util.EPSILON)+(1-m)*util.DBLE_MAX
221	else:
222	eta_max=util.DBLE_MAX
223	# initial guess:
224	tau=util.length(self.__stress)/util.sqrt(2)
225	# start the iteration:
226	cc=0
227	TOL=1e-7
228	dtau=util.DBLE_MAX
229	print "tau = ", tau, "eps =",eps
230	while cc<10 and dtau>TOL*util.Lsup(tau):
231	eta_eff2,eta_eff_dash=self.evalEtaEff(tau,return_dash=True)
232	eta_eff=util.clip(eta_eff2-eta_eff_dashtau/(1-eta_eff_dasheps),maxval=eta_max)
233	tau, tau2=eta_eff*eps, tau
234	dtau=util.Lsup(tau2-tau)
235	print "step ",cc,dtau, util.Lsup(tau)
236	cc+=1
237	return eta_eff
238
239	def getEtaCharacteristic(self):
240	a=0
241	for i in xrange(self.__numMaterials):
242	a=a+1./self.__eta_N[i]
243	return 1/a
244
245	def evalEtaEff(self, tau, return_dash=False):
246	a=Scalar(0,tau.getFunctionSpace()) # =1/eta
247	if return_dash: a_dash=Scalar(0,tau.getFunctionSpace()) # =(1/eta)'
248	s=util.Lsup(tau)
249	if s>0:
250	m=util.wherePositive(tau)
251	tau2=sutil.EPSILON(1.-m)+m*tau
252	for i in xrange(self.__numMaterials):
253	eta_N=self.__eta_N[i]
254	tau_t=self.__tau_t[i]
255	if tau_t==None:
256	a+=1./eta_N
257	else:
258	power=1.-1./self.__power[i]
259	c=1./(tau_t*powereta_N)
260	a+=ctau2*power
261	if return_dash: a_dash+=powerctau2**(power-1.)
262	else:
263	for i in xrange(self.__numMaterials):
264	eta_N=self.__eta_N[i]
265	power=1.-1./self.__power[i]
266	a+=util.whereZero(power)/eta_N
267	if self.__mu!=None: a+=1./(self.__dt*self.__mu)
268	out=1/a
269	if return_dash:
270	return out,-out*2a_dash
271	else:
272	return out
273
274	def eval(self,arg):
275	v=arg[0]
276	p=arg[1]
277	D=self.getDeviatoricStrain(v)
278	eta_eff=self.getEtaEff(D,p)
279	print "eta_eff=",eta_eff
280	# solve for dv
281	self.__pde_v.setValue(A=Data()) # save memory!
282	k3=util.kronecker(Function(self.getDomain()))
283	k3Xk3=util.outer(k3,k3)
284	self.__pde_v.setValue(A=eta_eff(util.swap_axes(k3Xk3,0,3)+util.swap_axes(k3Xk3,1,3)),X=-eta_effD+p*util.kronecker(self.getDomain()))
285	dv=self.__pde_v.getSolution(verbose=self.__verbose)
286	print "resistep dv =",dv
287	# solve for dp
288	v2=v+dv
289	self.__pde_p.setValue(D=1/eta_eff,Y=util.div(v2))
290	dp=self.__pde_p.getSolution(verbose=self.__verbose)
291	print "resistep dp =",dp
292	return ArithmeticTuple(dv,dp)
293
294	def update(self,dt, iter_max=100, inner_iter_max=20, verbose=False):
295	"""
296	updates stress, velocity and pressure for time increment dt
297
298	@param dt: time increment
299	@param max_inner_iter: maximum number of iteration steps in the incompressible solver
300	@param verbose: prints some infos in the incompressible solve
301	@param show_details: prints some infos while solving PDEs
302	@param tol: tolerance for the time step
303	"""
304	self.__verbose=verbose
305	self.__dt=dt
306	tol=self.getTolerance()
307	# set the initial velocity:
308	m=util.wherePositive(self.__pde_v.getCoefficient("q"))
309	v_new=self.__v(1-m)+self.__v_boundarym
310	# and off we go:
311	x=ArithmeticTuple(v_new, self.__p)
312	# self.__p_weight=util.interpolate(1./self.getEtaCharacteristic(),Function(self.__domain))**2
313	self.__p_weight=self.getEtaCharacteristic()
314	# self.__p_weight=util.interpolate(1./self.getEtaCharacteristic()**2,self.__p.getFunctionSpace())
315	atol=self.norm(x)*self.__tol
316	x_new=NewtonGMRES(self, x, iter_max=iter_max,sub_iter_max=inner_iter_max, atol=atol,rtol=0., verbose=verbose)
317	self.__v=x_new[0]
318	self.__p=x_new[1]
319	1/0
320	# self.__stress=self.getUpdatedStress(...)
321	self.__t+=dt
322	return self.__v, self.__p
323
324	#=========================================================================================
325
326	def getNewDeviatoricStress(self,D,eta_eff=None):
327	if eta_eff==None: eta_eff=self.evalEtaEff(self.__stress,D,self.__p)
328	s=(2eta_eff)D
329	if self.__mu!=None: s+=eta_eff/(self.__dtself.__mu)self.__last_stress
330	return s
331
332	def getDeviatoricStress(self):
333	"""
334	returns current stress
335	"""
336	return self.__stress
337	def getDeviatoricStrain(self,velocity=None):
338	"""
339	return strain
340	"""
341	if velocity==None: velocity=self.getVelocity()
342	return util.deviatoric(util.symmetric(util.grad(velocity)))
343
344	def getPressure(self):
345	"""
346	returns current pressure
347	"""
348	return self.__p
349	def getVelocity(self):
350	"""
351	returns current velocity
352	"""
353	return self.__v
354
355	def getTau(self,stress=None):
356	"""
357	returns current second stress deviatoric invariant
358	"""
359	if stress==None: stress=self.getDeviatoricStress()
360	return util.sqrt(0.5)*util.length(stress)
361
362	def getGammaDot(self,strain=None):
363	"""
364	returns current second stress deviatoric invariant
365	"""
366	if strain==None: strain=self.getDeviatoricStrain()
367	return util.sqrt(2)*util.length(strain)
368