/[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 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years ago) by ksteube
File MIME type: text/x-python
File size: 16158 byte(s)
Copyright updated in all python files

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 HomogeneousSaddlePointProblem,Projector
39
40 class PlateMantelModel(HomogeneousSaddlePointProblem):
41 """
42 This class implements the reology of a self-consistent plate-mantel model proposed in
43 U{Hans-Bernd Muhlhaus<emailto:h.muhlhaus@uq.edu.au>} and U{Klaus Regenauer-Lieb<mailto:klaus.regenauer-lieb@csiro.au>}:
44 I{Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection},
45 see U{doi: 10.1111/j.1365-246X.2005.02742.x<http://www3.interscience.wiley.com/journal/118661486/abstract?CRETRY=1&SRETRY=0>}.
46
47 typical usage:
48
49 sp=PlateMantelModel(domain,stress=0,v=0)
50 sp.setTolerance()
51 sp.initialize(...)
52 v,p=sp.solve(v0,p0)
53 """
54 def __init__(self, domain, stress=0, v=0, p=0, t=0, useJaumannStress=True, **kwargs):
55 """
56 initializes PlateMantelModel
57
58 @param domain: problem domain
59 @type domain: L{domain}
60 @param stress: initial stress
61 @param v: initial velocity field
62 @param t: initial time
63 @param useJaumannStress: C{True} if Jaumann stress is used
64 """
65 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
66 self.__domain=domain
67 self.__t=t
68 self.setSmall()
69 self.__vol=util.integrate(1.,Function(self.__domain))
70 self.__useJaumannStress=useJaumannStress
71 #=======================
72 # state variables:
73 #
74 if isinstance(stress,Data):
75 self.__stress=util.interpolate(stress,Function(domain))
76 else:
77 self.__stress=Data(stress,(domain.getDim(),domain.getDim()),Function(domain))
78 self.__stress-=util.trace(self.__stress)*(util.kronecker(domain)/domain.getDim())
79 if isinstance(v,Data):
80 self.__v=util.interpolate(v,Solution(domain))
81 else:
82 self.__v=Data(v,(domain.getDim(),),Solution(domain))
83 if isinstance(p,Data):
84 self.__p=util.interpolate(p,ReducedSolution(domain))
85 else:
86 self.__p=Data(p,(),ReducedSolution(domain))
87 self.__tau=util.sqrt(0.5)*util.length(self.__stress)
88 self.__D=util.symmetric(util.grad(self.__v))
89 #=======================
90 # parameters
91 #
92 self.__mu=None
93 self.__eta_N=None
94 self.__eta_0=None
95 self.__tau_0=None
96 self.__n=None
97 self.__eta_Y=None
98 self.__tau_Y=None
99 self.__n_Y=None
100 #=======================
101 # solme value to track:
102 #
103 self.__mu_eff=None
104 self.__h_eff=None
105 self.__eta_eff=None
106 #=======================
107 # PDE related stuff
108 self.__pde_u=LinearPDE(domain,numEquations=self.getDomain().getDim(),numSolutions=self.getDomain().getDim())
109 self.__pde_u.setSymmetryOn()
110 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
111
112 self.__pde_prec=LinearPDE(domain)
113 self.__pde_prec.setReducedOrderOn()
114 self.__pde_prec.setSymmetryOn()
115
116 self.__pde_proj=LinearPDE(domain)
117 self.__pde_proj.setReducedOrderOn()
118 self.__pde_proj.setSymmetryOn()
119 self.__pde_proj.setValue(D=1.)
120
121 def useJaumannStress(self):
122 """
123 return C{True} if Jaumann stress is included.
124 """
125 return self.__useJaumannStress
126 def setSmall(self,small=util.sqrt(util.EPSILON)):
127 """
128 sets small value
129
130 @param small: positive small value
131 """
132 self.__small=small
133 def getSmall(self):
134 """
135 returns small value
136 @rtype: positive float
137 """
138 return self.__small
139 def getDomain(self):
140 """
141 returns the domain
142 """
143 return self.__domain
144 def getStress(self):
145 """
146 returns current stress
147 """
148 return self.__stress
149 def getStretching(self):
150 """
151 return stretching
152 """
153 return self.__D
154 def getTau(self):
155 """
156 returns current second stress deviatoric invariant
157 """
158 return self.__tau
159 def getPressure(self):
160 """
161 returns current pressure
162 """
163 return self.__p
164 def getVelocity(self):
165 """
166 returns current velocity
167 """
168 return self.__v
169 def getTime(self):
170 """
171 returns current time
172 """
173 return self.__t
174 def getVolume(self):
175 """
176 returns domain volume
177 """
178 return self.__vol
179
180 def getEtaEffective(self):
181 """
182 returns effective shear viscocity
183 """
184 return self.__eta_eff
185 def getMuEffective(self):
186 """
187 returns effective shear modulus
188 """
189 return self.__mu_eff
190 def getHEffective(self):
191 """
192 returns effective h
193 """
194 return self.__h_eff
195 def getMechanicalPower(self):
196 """
197 returns the locl mechanical power M{D_ij Stress_ij}
198 """
199 return util.inner(self.getStress(),self.getStretching())
200
201 def initialize(self, mu=None, eta_N=None, eta_0=None, tau_0=None, n=None, eta_Y=None, tau_Y=None, n_Y=None, F=None, q=None):
202 """
203 sets the model parameters.
204
205 @param mu: shear modulus. If C{mu==None} ...
206 @param eta_N: Newtonian viscosity.
207 @param eta_0: power law viscovity for tau=tau0
208 @param tau_0: reference tau for power low. If C{tau_0==None} constant viscosity is assumed.
209 @param n: power of power law. if C{n=None}, constant viscosity is assumed.
210 @param eta_Y: =None
211 @param tau_Y: =None
212 @param n_Y: =None
213 @param F: external force.
214 @param q: location of constraints
215 """
216 if mu != None: self.__mu=mu
217 if eta_N != None: self.__eta_N=eta_N
218 if eta_0 != None: self.__eta_0=eta_0
219 if tau_0 != None: self.__tau_0=tau_0
220 if n != None: self.__n=n
221 if eta_Y != None: self.__eta_Y=eta_Y
222 if tau_Y != None: self.__tau_Y=tau_Y
223 if n_Y != None: self.__n_Y=n_Y
224 if F != None: self.__pde_u.setValue(Y=F)
225 if q != None: self.__pde_u.setValue(q=q)
226
227 def updateEffectiveCoefficients(self, dt, tau):
228 """
229 updates the effective coefficints depending on tau and time step size dt.
230 """
231 h_vis=None
232 h_yie=None
233 # calculate eta_eff = ( 1/eta_N + 1/eta_vis + 1/eta_yie)^{-1}
234 # with eta_vis = eta_0 * (tau/tau_0) ^ (1-n)
235 # eta_yie = eta_Y * (tau/tau_Y) ^ (1-n_Y)
236 s=Scalar(0.,Function(self.getDomain()))
237 if self.__eta_N!=None:
238 s+=1/self.__eta_N
239 if self.__eta_0!=None and self.__tau_0 != None and self.__n!=None:
240 if self.__tau_0 !=None and self.__n!=None:
241 eta_vis=self.__eta_0*(tau/self.__tau_0)**(1-self.__n)
242 s+=1./(eta_vis+self.__eta_0*self.getSmall())
243 h_vis=eta_vis/(self.__n-1)
244 if self.__tau_Y!=None and self.__eta_Y!=None and self.__n_Y!=None:
245 eta_yie=self.__eta_Y*(tau/self.__tau_Y)**(1-self.__n_Y)
246 s+=1/(eta_yie+self.getSmall()*self.__eta_Y)
247 h_yie=self.__eta_Y/(self.__n_Y-1)
248 self.__eta_eff=1/s
249 # calculate eta_eff = ( 1/h_vis + 1/h_yie)^{-1}
250 # with h_vis = eta_vis/(n-1)
251 # h_yie = eta_yie/(n_Y-1)
252 if h_vis == None:
253 if h_yie==None:
254 self__h_eff=None
255 else:
256 self__h_eff=h_yie
257 else:
258 if h_yie==None:
259 self__h_eff=h_vis
260 else:
261 self__h_eff=1/((1./h_vis)+(1./h_yie))
262 # calculate mu_eff = ( 1/mu + dt/eta_eff)^{-1} = mu*eta_eff/(mu*dt+eta_eff)
263 if self.__mu == None:
264 self.__mu_eff=self.__eta_eff/dt
265 else:
266 self.__mu_eff=1./(1./self.__mu+dt/self.__eta_eff)
267
268 def update(self,dt,max_inner_iter=20, verbose=False, show_details=False, tol=10., solver="PCG"):
269 """
270 updates stress, velocity and pressure for time increment dt
271
272 @param dt: time increment
273 @param max_inner_iter: maximum number of iteration steps in the incompressible solver
274 @param verbose: prints some infos in the incompressible solve
275 @param show_details: prints some infos while solving PDEs
276 @param tol: tolerance for the time step
277 """
278 stress_last=self.getStress()
279 # we should use something like FGMRES to merge the two iterations!
280 e=10.*tol
281 # get values from last time step and use them as initial guess:
282 v=self.__v
283 stress=self.__stress
284 p=self.__p
285 tau=self.__tau
286 while e > tol: # revise stress calculation if this is used.
287 #
288 # update the effective coefficients:
289 #
290 self.updateEffectiveCoefficients(dt,tau)
291 eta_eff=self.getEtaEffective()
292 mu_eff=self.getMuEffective()
293 h_eff=self.getHEffective()
294 #
295 # create some temporary variables:
296 #
297 self.__pde_u.setValue(A=Data()) # save memory!
298 k3=util.kronecker(Function(self.getDomain()))
299 k3Xk3=util.outer(k3,k3)
300 self.__f3=mu_eff*dt
301 A=self.__f3*(util.swap_axes(k3Xk3,0,3)+util.swap_axes(k3Xk3,1,3))
302 print "mueff=",util.inf(mu_eff),util.sup(mu_eff)
303 if h_eff == None:
304 s=0
305 self.__f0=0
306 self.__f1=mu_eff/eta_eff*dt
307 else:
308 s=mu_eff*dt/(h_eff+mu_eff*dt)
309 Lsup_tau=Lsup(tau)
310 if Lsup>0:
311 self.__f0=mu_eff*s*dt/(tau+Lsup_tau*self.getSmall())**2
312 A+=util.outer((-self.__f0)*stress,stress)
313 else:
314 self.__f0=0.
315 self.__f1=mu_eff/eta_eff*(1-s)*dt
316
317 if self.useJaumannStress():
318 self.__f2=mu_eff/eta_eff*dt**2
319 sXk3=util.outer(stress,self.__f2*(-k3/2))
320
321 A+=util.swap_axes(sXk3,0,3)
322 A+=util.swap_axes(sXk3,1,3)
323 A-=util.swap_axes(sXk3,0,2)
324 A-=util.swap_axes(sXk3,1,2)
325 else:
326 self.__f2=0
327 self.__pde_u.setValue(A=A)
328 self.__pde_prec.setValue(D=1/mu_eff)
329
330 print "X f0:",util.inf(self.__f0), util.sup(self.__f0)
331 print "X f1:",util.inf(self.__f1), util.sup(self.__f1)
332 print "X f2:",util.inf(self.__f2), util.sup(self.__f2)
333 print "X f3:",util.inf(self.__f3), util.sup(self.__f3)
334
335 v_old=v
336 v,p=self.solve(v,p,max_iter=max_inner_iter, verbose=verbose, show_details=show_details, solver=solver)
337 # update stress
338 stress=stress_last+dt*self.getStressChange(v)
339 stress-=util.trace(stress)*(util.kronecker(self.getDomain())/self.getDomain().getDim())
340 tau=util.sqrt(0.5)*util.length(stress)
341 # calculate error:
342 e=util.Lsup(v_old-v)/util.Lsup(v)
343 # state variables:
344 self.__t+=dt
345 self.__v=v
346 self.__p=p
347 self.__stress=stress
348 self.__tau=tau
349 self.__D=util.symmetric(util.grad(v))
350
351 def getStressChange(self,v):
352 """
353 returns the stress change due to a given velocity field v
354 """
355 stress=self.getStress()
356 g=util.grad(v)
357 D=util.symmetric(g)
358 W=util.nonsymmetric(g)
359 U=util.nonsymmetric(util.tensor_mult(W,stress))
360 dstress=2*self.__f3*D-self.__f0*util.inner(stress,D)*stress-self.__f1*stress-2*self.__f2*U
361 return dstress
362
363 def B(self,arg):
364 """
365 div operator
366 """
367 d=util.div(arg)
368 self.__pde_proj.setValue(Y=d)
369 self.__pde_proj.setTolerance(self.getSubProblemTolerance())
370 return self.__pde_proj.getSolution(verbose=self.show_details)
371
372 def solve_prec(self,p):
373 """
374 preconditioner
375 """
376 #proj=Projector(domain=self.getDomain(), reduce = True, fast=False)
377 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
378 self.__pde_prec.setValue(Y=p)
379 q=self.__pde_prec.getSolution(verbose=self.show_details)
380 return q
381
382 def inner(self,p0,p1):
383 """
384 inner product for pressure
385 """
386 s0=util.interpolate(p0,Function(self.getDomain()))
387 s1=util.interpolate(p1,Function(self.getDomain()))
388 return util.integrate(s0*s1)
389
390 def solve_A(self,u,p):
391 """
392 solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
393 """
394 self.__pde_u.setTolerance(self.getSubProblemTolerance())
395 self.__pde_u.setValue(X=-self.getStressChange(u)-p*util.kronecker(self.getDomain()))
396 return self.__pde_u.getSolution(verbose=self.show_details)
397
398 ##############Added by Artak ##########################
399 def inner_a(self,a0,a1):
400 p0=util.interpolate(a0[1],Function(self.domain))
401 p1=util.interpolate(a1[1],Function(self.domain))
402 alfa=(1/self.vol)*util.integrate(p0)
403 beta=(1/self.vol)*util.integrate(p1)
404 v0=util.grad(a0[0])
405 v1=util.grad(a1[0])
406 return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))
407
408 def stoppingcriterium(self,Bv,v,p):
409 n_r=util.sqrt(self.inner(Bv,Bv))
410 n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
411 if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)
412 if self.iter == 0: self.__n_v=n_v;
413 self.__n_v, n_v_old =n_v, self.__n_v
414 self.iter+=1
415 if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
416 if self.verbose: print "PCG terminated after %s steps."%self.iter
417 return True
418 else:
419 return False
420
421 def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
422 if TOL==None:
423 TOL=self.getTolerance()
424 if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
425 self.iter+=1
426
427 if norm_r <= norm_b*TOL:
428 if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
429 return True
430 else:
431 return False
432
433 ######################################################

  ViewVC Help
Powered by ViewVC 1.1.26