/[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 2100 - (show annotations)
Wed Nov 26 08:13:00 2008 UTC (10 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 13741 byte(s)
This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module. 
Some documentation for both classes has been added.
The convection code is only linear at the moment.



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**2*util.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/((2*self.__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=m*tmp/(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_dash*tau/(1-eta_eff_dash*eps),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=s*util.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**power*eta_N)
260 a+=c*tau2**power
261 if return_dash: a_dash+=power*c*tau2**(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**2*a_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_eff*D+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_boundary*m
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=(2*eta_eff)*D
329 if self.__mu!=None: s+=eta_eff/(self.__dt*self.__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

  ViewVC Help
Powered by ViewVC 1.1.26