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

  ViewVC Help
Powered by ViewVC 1.1.26