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

  ViewVC Help
Powered by ViewVC 1.1.26