/[escript]/trunk/escript/py_src/rheologies.py
ViewVC logotype

Annotation of /trunk/escript/py_src/rheologies.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1639 - (hide annotations)
Mon Jul 14 08:55:25 2008 UTC (11 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 14897 byte(s)


1 gross 1639 # $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 getStrainEnergy(self):
197     """
198     returns strain energy density
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     v_old=v
332     v,p=self.solve(v,p,max_iter=max_inner_iter, verbose=verbose, show_details=show_details, solver=solver)
333     # update stress
334     stress=stress_last+dt*self.getStressChange(v)
335     stress-=util.trace(stress)*(util.kronecker(self.getDomain())/self.getDomain().getDim())
336     tau=util.sqrt(0.5)*util.length(stress)
337     # calculate error:
338     e=util.Lsup(v_old-v)/util.Lsup(v)
339     # state variables:
340     self.__t+=dt
341     self.__v=v
342     self.__p=p
343     self.__stress=stress
344     self.__tau=tau
345     self.__D=util.symmetric(util.grad(v))
346    
347     def getStressChange(self,v):
348     """
349     returns the stress change due to a given velocity field v
350     """
351     stress=self.getStress()
352     g=util.grad(v)
353     D=util.symmetric(g)
354     W=util.nonsymmetric(g)
355     U=util.nonsymmetric(util.tensor_mult(W,stress))
356     dstress=2*self.__f3*D-self.__f0*util.inner(stress,D)*stress-self.__f1*stress-2*self.__f2*U
357     return dstress
358    
359     def B(self,arg):
360     """
361     div operator
362     """
363     d=util.div(arg)
364     self.__pde_proj.setValue(Y=d)
365     self.__pde_proj.setTolerance(self.getSubProblemTolerance())
366     return self.__pde_proj.getSolution(verbose=self.show_details)
367    
368     def solve_prec(self,p):
369     """
370     preconditioner
371     """
372     #proj=Projector(domain=self.getDomain(), reduce = True, fast=False)
373     self.__pde_prec.setTolerance(self.getSubProblemTolerance())
374     self.__pde_prec.setValue(Y=p)
375     q=self.__pde_prec.getSolution(verbose=self.show_details)
376     return q
377    
378     def inner(self,p0,p1):
379     """
380     inner product for pressure
381     """
382     s0=util.interpolate(p0,Function(self.getDomain()))
383     s1=util.interpolate(p1,Function(self.getDomain()))
384     return util.integrate(s0*s1)
385    
386     def solve_A(self,u,p):
387     """
388     solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
389     """
390     self.__pde_u.setTolerance(self.getSubProblemTolerance())
391     self.__pde_u.setValue(X=-self.getStressChange(u)-p*util.kronecker(self.getDomain()))
392     return self.__pde_u.getSolution(verbose=self.show_details)
393    
394    
395     def stoppingcriterium(self,Bv,v,p):
396     n_r=util.sqrt(self.inner(Bv,Bv))
397     n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
398     if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)
399     if self.iter == 0: self.__n_v=n_v;
400     self.__n_v, n_v_old =n_v, self.__n_v
401     self.iter+=1
402     if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
403     if self.verbose: print "PCG terminated after %s steps."%self.iter
404     return True
405     else:
406     return False

  ViewVC Help
Powered by ViewVC 1.1.26