/[escript]/trunk/escript/py_src/rheologies.py
ViewVC logotype

Contents of /trunk/escript/py_src/rheologies.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2301 - (show annotations)
Thu Mar 12 01:35:26 2009 UTC (10 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 22562 byte(s)
test for power law added
1 ########################################################
2 #
3 # Copyright (c) 2003-2008 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 ########################################################
12
13 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14 Earth Systems Science Computational Center (ESSCC)
15 http://www.uq.edu.au/esscc
16 Primary Business: Queensland, Australia"""
17 __license__="""Licensed under the Open Software License version 3.0
18 http://www.opensource.org/licenses/osl-3.0.php"""
19 __url__="http://www.uq.edu.au/esscc/escript-finley"
20
21 """
22 Some models for flow
23
24 @var __author__: name of author
25 @var __copyright__: copyrights
26 @var __license__: licence agreement
27 @var __url__: url entry point on documentation
28 @var __version__: version
29 @var __date__: date of the version
30 """
31
32 __author__="Lutz Gross, l.gross@uq.edu.au"
33
34 from escript import *
35 import util
36 from linearPDEs import LinearPDE
37 from pdetools import Defect, NewtonGMRES, ArithmeticTuple
38
39 class PowerLaw(object):
40 """
41 this implements the power law for a composition of a set of materials where the viscosity eta of each material is given by a
42 power law relationship of the form
43
44 M{eta=eta_N*(tau/tau_t)**(1.-1./power)}
45
46 where tau is equivalent stress and eta_N, tau_t and power are given constant. Moreover an elastic component can be considered.
47 Moreover tau meets the Drucker-Prager type yield condition
48
49 M{tau <= tau_Y + friction * pressure}
50
51 where gamma_dot is the equivalent.
52 """
53 def __init__(self, numMaterials=1,verbose=False):
54 """
55 initializes a power law
56
57 @param numMaterials: number of materials
58 @type numMaterials: C{int}
59 @param verbose: if C{True} some informations are printed.
60 @type verbose: C{bool}
61 """
62 if numMaterials<1:
63 raise ValueError,"at least one material must be defined."
64 self.__numMaterials=numMaterials
65 self.__eta_N=[None for i in xrange(self.__numMaterials)]
66 self.__tau_t=[1. for i in xrange(self.__numMaterials)]
67 self.__power=[1. for i in xrange(self.__numMaterials)]
68 self.__tau_Y=None
69 self.__friction=None
70 self.__mu=None
71 self.__verbose=verbose
72 self.setEtaTolerance()
73 #===========================================================================
74 def getNumMaterials(self):
75 """
76 returns the numebr of materials
77
78 @return: number of materials
79 @rtype: C{int}
80 """
81 return self.__numMaterials
82 def validMaterialId(self,id=0):
83 """
84 checks if a given material id is valid
85
86 @param id: a material id
87 @type id: C{int}
88 @return: C{True} is the id is valid
89 @rtype: C{bool}
90 """
91 return 0<=id and id<self.getNumMaterials()
92 def setEtaTolerance(self,rtol=util.sqrt(util.EPSILON)):
93 """
94 sets the relative tolerance for the effectice viscosity.
95
96 @param rtol: relative tolerance
97 @type rtol: positive C{float}
98 """
99 if rtol<=0:
100 raise ValueError,"rtol needs to positive."
101 self.__rtol=rtol
102 def getEtaTolerance(self):
103 """
104 returns the relative tolerance for the effectice viscosity.
105
106 @return: relative tolerance
107 @rtype rtol: positive C{float}
108 """
109 return self.__rtol
110 #===========================================================================
111 def setDruckerPragerLaw(self,tau_Y=None,friction=None):
112 """
113 Sets the parameters for the Drucker-Prager model.
114
115 @param tau_Y: yield stress
116 @param friction: friction coefficient
117 """
118 self.__tau_Y=tau_Y
119 self.__friction=friction
120 def getFriction(self):
121 """
122 returns the friction coefficient
123
124 @return: friction coefficient
125 """
126 return self.__friction
127 def getTauY(self):
128 """
129 returns the yield stress
130
131 @return: the yield stress
132 """
133 return self.__tau_Y
134 #===========================================================================
135 def getElasticShearModulus(self):
136 """
137 returns the elastic shear modulus.
138
139 @return: elastic shear modulus
140 """
141 return self.__mu
142 def setElasticShearModulus(self,mu=None):
143 """
144 Sets the elastic shear modulus.
145
146 @param mu: elastic shear modulus
147 """
148 self.__mu=mu
149 #===========================================================================
150 def getPower(self, id=None):
151 """
152 returns the power in the power law
153
154 @param id: if present, the power for material C{id} is returned.
155 @type id: C{int}
156 @return: the list of the powers for all matrials is returned. If C{id} is present only the power for material C{id} is returned.
157 """
158 if id == None:
159 return self.__power
160 else:
161 if self.validMaterialId(id):
162 return self.__power[id]
163 else:
164 raise ValueError,"Illegal material id %s."%id
165 def getEtaN(self, id=None):
166 """
167 returns the viscosity
168
169 @param id: if present, the viscosity for material C{id} is returned.
170 @type id: C{int}
171 @return: the list of the viscosities for all matrials is returned. If C{id} is present only the viscosity for material C{id} is returned.
172 """
173 if id == None:
174 return self.__eta_N
175 else:
176 if self.validMaterialId(id):
177 return self.__eta_N[id]
178 else:
179 raise ValueError,"Illegal material id %s."%id
180 def getTauT(self, id=None):
181 """
182 returns the transition stress
183
184 @param id: if present, the transition stress for material C{id} is returned.
185 @type id: C{int}
186 @return: the list of the transition stresses for all matrials is returned. If C{id} is present only the transition stress for material C{id} is returned.
187 """
188 if id == None:
189 return self.__tau_t
190 else:
191 if self.validMaterialId(id):
192 return self.__tau_t[id]
193 else:
194 raise ValueError,"Illegal material id %s."%id
195
196 def setPowerLaw(self,eta_N, id=0, tau_t=1, power=1):
197 """
198 Sets the power-law parameters for material id
199
200 @param id: material id
201 @type id: C{int}
202 @param eta_N: viscosity for tau=tau_t
203 @param tau_t: transition stress
204 @param power: power law coefficient
205 """
206 if self.validMaterialId(id):
207 self.__eta_N[id]=eta_N
208 self.__power[id]=power
209 self.__tau_t[id]=tau_t
210 else:
211 raise ValueError,"Illegal material id %s."%id
212
213 def setPowerLaws(self,eta_N, tau_t, power):
214 """
215 Sets the parameters of the power-law for all materials.
216
217 @param eta_N: list of viscosities for tau=tau_t
218 @param tau_t: list of transition stresses
219 @param power: list of power law coefficient
220 """
221 if len(eta_N)!=self.__numMaterials or len(tau_t)!=self.__numMaterials or len(power)!=self.__numMaterials:
222 raise ValueError,"%s materials are expected."%self.__numMaterials
223 for i in xrange(self.__numMaterials):
224 self.setPowerLaw(id=i, eta_N=eta_N[i],tau_t=tau_t[i],power=power[i])
225
226 #===========================================================================
227 def getEtaEff(self,gamma_dot, eta0=None, pressure=None,dt=None, iter_max=10):
228 """
229 returns the effective viscosity eta_eff such that
230
231 M{tau=eta_eff * gamma_dot}
232
233 by solving a non-linear problem for tau.
234
235 @param gamma_dot: equivalent strain gamma_dot
236 @param eta0: initial guess for the effective viscosity (e.g from a previous time step). If not present, an initial guess is calculated.
237 @param pressure: pressure used to calculate yield condition
238 @param dt: time step size. only needed if elastic component is considered.
239 @type dt: positive C{float} if present
240 @param iter_max: maximum number of iteration steps.
241 @type iter_max: C{int}
242 @return: effective viscosity.
243 """
244 SMALL=1./(util.DBLE_MAX/100.)
245 numMaterial=self.getNumMaterials()
246 s=[1.-1./p for p in self.getPower() ]
247 eta_N=self.getEtaN()
248 tau_t=self.getTauT()
249 mu=self.getElasticShearModulus()
250 fric=self.getFriction()
251 tau_Y=self.getTauY()
252 if eta0==None:
253 theta=0.
254 for i in xrange(numMaterial):
255 inv_eta_i=0**s[i]/eta_N[i]
256 theta=theta+inv_eta_i
257 if util.inf(theta)<=0:
258 raise ValueError,"unable to set positive initial guess for eta_eff. Most likely no power law with power 1 set."
259 eta_eff=1./theta
260 else:
261 if util.inf(eta0)<=0:
262 raise ValueError,"initial guess for eta_eff is not positive."
263 eta_eff=eta0
264
265 if mu !=None and dt == None:
266 raise ValueError,"Time stepsize dt must be given."
267 if dt !=None:
268 if dt<=0: raise ValueError,"time step size must be positive."
269 if tau_Y==None and fric==None:
270 eta_max=None
271 else:
272 if fric == None:
273 eta_max=tau_Y/(gamma_dot+SMALL*util.whereZero(gamma_dot))
274 else:
275 if tau_Y==None: tau_Y==0
276 if util.inf(fric)<=0:
277 raise ValueError,"if friction present it needs to be positive."
278 eta_max=fric*util.clip(tau_Y/fric+pressure,minval=0)/(gamma_dot+SMALL*util.whereZero(gamma_dot))
279 rtol=self.getEtaTolerance()
280 iter =0
281 converged=False
282 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))
284 while not converged:
285 if iter>max(iter_max,1):
286 raise RuntimeError,"tolerance not reached after %s steps."%max(iter_max,1)
287 #===========================================
288 theta=0. # =1/eta
289 omega=0. # = tau*theta'= eta'*tau/eta**2
290 if mu !=None: theta=1./(dt*mu)
291 for i in xrange(numMaterial):
292 inv_eta_i=(tau/tau_t[i])**s[i]/eta_N[i]
293 theta=theta+inv_eta_i
294 omega=omega+s[i]*inv_eta_i
295 #===========================================
296 eta_eff, eta_eff_old=util.clip(eta_eff*(theta+omega)/(eta_eff*theta**2+omega),maxval=eta_max), eta_eff
297 tau=eta_eff*gamma_dot
298 d=util.Lsup(eta_eff-eta_eff_old)
299 l=util.Lsup(eta_eff)
300 iter+=1
301 if self.__verbose: print "step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
302 converged= d<= rtol* l
303 return eta_eff
304
305 #====================================================================================================================================
306
307 class IncompressibleIsotropicKelvinFlow(Defect):
308 """
309 This class implements the rheology of an isotropic Kelvin material.
310
311 Typical usage::
312
313 sp = IncompressibleIsotropicKelvinFlow(domain, stress=0, v=0)
314 sp.setTolerance()
315 sp.initialize(...)
316 v,p = sp.solve(v0, p0)
317
318 @note: This model has been used in the self-consistent plate-mantle model
319 proposed in U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>}
320 and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
321 I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
322 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract?CRETRY=1&SRETRY=0>}.
323
324 """
325 def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, useJaumannStress=True, **kwargs):
326 """
327 Initializes the model.
328
329 @param domain: problem domain
330 @type domain: L{domain}
331 @param stress: initial deviatoric stress
332 @param v: initial velocity field
333 @param p: initial pressure
334 @param t: initial time
335 @param useJaumannStress: C{True} if Jaumann stress is used
336 (not supported yet)
337 """
338 super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)
339 self.__domain=domain
340 self.__t=t
341 self.__vol=util.integrate(1.,Function(self.__domain))
342 self.__useJaumannStress=useJaumannStress
343 self.__v_boundary=Vector(0,Solution(self.__domain))
344 #=======================
345 # state variables:
346 #
347 if isinstance(stress,Data):
348 self.__stress=util.interpolate(stress,Function(domain))
349 else:
350 self.__stress=Data(stress,(domain.getDim(),domain.getDim()),Function(domain))
351 self.__stress-=util.trace(self.__stress)*(util.kronecker(domain)/domain.getDim())
352 if isinstance(v,Data):
353 self.__v=util.interpolate(v,Solution(domain))
354 else:
355 self.__v=Data(v,(domain.getDim(),),Solution(domain))
356 if isinstance(p,Data):
357 self.__p=util.interpolate(p,ReducedSolution(domain))
358 else:
359 self.__p=Data(p,(),ReducedSolution(domain))
360 #=======================
361 # PDE related stuff
362 self.__pde_v=LinearPDE(domain,numEquations=self.getDomain().getDim(),numSolutions=self.getDomain().getDim())
363 self.__pde_v.setSymmetryOn()
364 self.__pde_v.setSolverMethod(preconditioner=LinearPDE.RILU)
365
366 self.__pde_p=LinearPDE(domain)
367 self.__pde_p.setReducedOrderOn()
368 self.__pde_p.setSymmetryOn()
369
370 self.setTolerance()
371 self.setSmall()
372
373 def useJaumannStress(self):
374 """
375 Returns C{True} if Jaumann stress is included.
376 """
377 return self.__useJaumannStress
378
379 def setSmall(self,small=util.sqrt(util.EPSILON)):
380 """
381 Sets a small value to be used.
382
383 @param small: positive small value
384 """
385 self.__small=small
386
387 def getSmall(self):
388 """
389 Returns small value.
390 @rtype: positive float
391 """
392 return self.__small
393
394 def setTolerance(self,tol=1.e-4):
395 """
396 Sets the tolerance.
397 """
398 self.__pde_v.setTolerance(tol**2)
399 self.__pde_p.setTolerance(tol**2)
400 self.__tol=tol
401
402 def getTolerance(self):
403 """
404 Returns the set tolerance.
405 @rtype: positive float
406 """
407 return self.__tol
408
409 def getDomain(self):
410 """
411 Returns the domain.
412 """
413 return self.__domain
414
415 def getTime(self):
416 """
417 Returns current time.
418 """
419 return self.__t
420
421 def setExternals(self, F=None, f=None, q=None, v_boundary=None):
422 """
423 Sets externals.
424
425 @param F: external force
426 @param f: surface force
427 @param q: location of constraints
428 """
429 if F != None: self.__pde_v.setValue(Y=F)
430 if f != None: self.__pde_v.setValue(y=f)
431 if q != None: self.__pde_v.setValue(q=q)
432 if v_boundary != None: self.__v_boundary=v_boundary
433
434 def bilinearform(self,arg0,arg1):
435 s0=util.deviatoric(util.symmetric(util.grad(arg0[0])))
436 s1=util.deviatoric(util.symmetric(util.grad(arg1[0])))
437 # s0=util.interpolate(arg0[0],Function(self.getDomain()))
438 # s1=util.interpolate(arg1[0],Function(self.getDomain()))
439 p0=util.interpolate(arg0[1],Function(self.getDomain()))
440 p1=util.interpolate(arg1[1],Function(self.getDomain()))
441 a=util.integrate(self.__p_weight**2*util.inner(s0,s1))+util.integrate(p0*p1)
442 return a
443
444 def getEtaEff(self,strain, pressure):
445 if self.__mu==None:
446 eps=util.length(strain)*util.sqrt(2)
447 else:
448 eps=util.length(strain+self.__stress/((2*self.__dt)*self.__mu))*util.sqrt(2)
449 p=util.interpolate(pressure,eps.getFunctionSpace())
450 if self.__tau_Y!= None:
451 tmp=self.__tau_Y+self.__friction*p
452 m=util.wherePositive(eps)*util.wherePositive(tmp)
453 eta_max=m*tmp/(eps+(1-m)*util.EPSILON)+(1-m)*util.DBLE_MAX
454 else:
455 eta_max=util.DBLE_MAX
456 # initial guess:
457 tau=util.length(self.__stress)/util.sqrt(2)
458 # start the iteration:
459 cc=0
460 TOL=1e-7
461 dtau=util.DBLE_MAX
462 print "tau = ", tau, "eps =",eps
463 while cc<10 and dtau>TOL*util.Lsup(tau):
464 eta_eff2,eta_eff_dash=self.evalEtaEff(tau,return_dash=True)
465 eta_eff=util.clip(eta_eff2-eta_eff_dash*tau/(1-eta_eff_dash*eps),maxval=eta_max)
466 tau, tau2=eta_eff*eps, tau
467 dtau=util.Lsup(tau2-tau)
468 print "step ",cc,dtau, util.Lsup(tau)
469 cc+=1
470 return eta_eff
471
472 def getEtaCharacteristic(self):
473 a=0
474 for i in xrange(self.__numMaterials):
475 a=a+1./self.__eta_N[i]
476 return 1/a
477
478 def evalEtaEff(self, tau, return_dash=False):
479 a=Scalar(0,tau.getFunctionSpace()) # =1/eta
480 if return_dash: a_dash=Scalar(0,tau.getFunctionSpace()) # =(1/eta)'
481 s=util.Lsup(tau)
482 if s>0:
483 m=util.wherePositive(tau)
484 tau2=s*util.EPSILON*(1.-m)+m*tau
485 for i in xrange(self.__numMaterials):
486 eta_N=self.__eta_N[i]
487 tau_t=self.__tau_t[i]
488 if tau_t==None:
489 a+=1./eta_N
490 else:
491 power=1.-1./self.__power[i]
492 c=1./(tau_t**power*eta_N)
493 a+=c*tau2**power
494 if return_dash: a_dash+=power*c*tau2**(power-1.)
495 else:
496 for i in xrange(self.__numMaterials):
497 eta_N=self.__eta_N[i]
498 power=1.-1./self.__power[i]
499 a+=util.whereZero(power)/eta_N
500 if self.__mu!=None: a+=1./(self.__dt*self.__mu)
501 out=1/a
502 if return_dash:
503 return out,-out**2*a_dash
504 else:
505 return out
506
507 def eval(self,arg):
508 v=arg[0]
509 p=arg[1]
510 D=self.getDeviatoricStrain(v)
511 eta_eff=self.getEtaEff(D,p)
512 print "eta_eff=",eta_eff
513 # solve for dv
514 self.__pde_v.setValue(A=Data()) # save memory!
515 k3=util.kronecker(Function(self.getDomain()))
516 k3Xk3=util.outer(k3,k3)
517 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()))
518 dv=self.__pde_v.getSolution(verbose=self.__verbose)
519 print "resistep dv =",dv
520 # solve for dp
521 v2=v+dv
522 self.__pde_p.setValue(D=1/eta_eff,Y=util.div(v2))
523 dp=self.__pde_p.getSolution(verbose=self.__verbose)
524 print "resistep dp =",dp
525 return ArithmeticTuple(dv,dp)
526
527 def update(self, dt, iter_max=100, inner_iter_max=20, verbose=False):
528 """
529 Updates stress, velocity and pressure for time increment dt.
530
531 @param dt: time increment
532 @param inner_iter_max: maximum number of iteration steps in the
533 incompressible solver
534 @param verbose: prints some infos in the incompressible solver
535 """
536 self.__verbose=verbose
537 self.__dt=dt
538 tol=self.getTolerance()
539 # set the initial velocity:
540 m=util.wherePositive(self.__pde_v.getCoefficient("q"))
541 v_new=self.__v*(1-m)+self.__v_boundary*m
542 # and off we go:
543 x=ArithmeticTuple(v_new, self.__p)
544 # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic(),Function(self.__domain))**2
545 self.__p_weight=self.getEtaCharacteristic()
546 # self.__p_weight=util.interpolate(1./self.getEtaCharacteristic()**2,self.__p.getFunctionSpace())
547 atol=self.norm(x)*self.__tol
548 x_new=NewtonGMRES(self, x, iter_max=iter_max,sub_iter_max=inner_iter_max, atol=atol,rtol=0., verbose=verbose)
549 self.__v=x_new[0]
550 self.__p=x_new[1]
551 1/0
552 # self.__stress=self.getUpdatedStress(...)
553 self.__t+=dt
554 return self.__v, self.__p
555
556 #========================================================================
557
558 def getNewDeviatoricStress(self,D,eta_eff=None):
559 if eta_eff==None: eta_eff=self.evalEtaEff(self.__stress,D,self.__p)
560 s=(2*eta_eff)*D
561 if self.__mu!=None: s+=eta_eff/(self.__dt*self.__mu)*self.__last_stress
562 return s
563
564 def getDeviatoricStress(self):
565 """
566 Returns current stress.
567 """
568 return self.__stress
569
570 def getDeviatoricStrain(self,velocity=None):
571 """
572 Returns strain.
573 """
574 if velocity==None: velocity=self.getVelocity()
575 return util.deviatoric(util.symmetric(util.grad(velocity)))
576
577 def getPressure(self):
578 """
579 Returns current pressure.
580 """
581 return self.__p
582
583 def getVelocity(self):
584 """
585 Returns current velocity.
586 """
587 return self.__v
588
589 def getTau(self,stress=None):
590 """
591 Returns current second stress deviatoric invariant.
592 """
593 if stress==None: stress=self.getDeviatoricStress()
594 return util.sqrt(0.5)*util.length(stress)
595
596 def getGammaDot(self,strain=None):
597 """
598 Returns current second stress deviatoric invariant.
599 """
600 if strain==None: strain=self.getDeviatoricStrain()
601 return util.sqrt(2)*util.length(strain)
602

  ViewVC Help
Powered by ViewVC 1.1.26