/[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 2169 - (show annotations)
Wed Dec 17 03:08:58 2008 UTC (10 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 13704 byte(s)
Assorted spelling, grammar, whitespace and copy/paste error fixes (Part 2).
All epydoc warnings for these files have been fixed.
This commit should be a no-op.

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

  ViewVC Help
Powered by ViewVC 1.1.26