/[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 2100 - (hide annotations)
Wed Nov 26 08:13:00 2008 UTC (10 years, 10 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 ksteube 1809
2     ########################################################
3 gross 1639 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1639 #
8 ksteube 1809 # 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 gross 1639 #
12 ksteube 1809 ########################################################
13 gross 1639
14 ksteube 1809 __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 gross 1639 """
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 gross 2100 from pdetools import Defect, NewtonGMRES, ArithmeticTuple
39 gross 1639
40 gross 2100 class IncompressibleIsotropicKelvinFlow(Defect):
41 gross 1639 """
42 gross 2100 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 gross 1639 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 gross 2100 def __init__(self, domain, stress=0, v=0, p=0, t=0, numMaterials=1, useJaumannStress=True, **kwargs):
57 gross 1639 """
58     initializes PlateMantelModel
59    
60     @param domain: problem domain
61     @type domain: L{domain}
62 gross 2100 @param stress: initial deviatoric stress
63 gross 1639 @param v: initial velocity field
64 gross 2100 @param p: initial pressure
65 gross 1639 @param t: initial time
66 gross 2100 @param useJaumannStress: C{True} if Jaumann stress is used (not supported yet)
67 gross 1639 """
68 gross 2100 if numMaterials<1:
69     raise ValueError,"at least one material must be defined."
70     super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)
71 gross 1639 self.__domain=domain
72     self.__t=t
73     self.__vol=util.integrate(1.,Function(self.__domain))
74     self.__useJaumannStress=useJaumannStress
75 gross 2100 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 gross 1639 #=======================
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 gross 2100 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 gross 1639
105 gross 2100 self.__pde_p=LinearPDE(domain)
106     self.__pde_p.setReducedOrderOn()
107     self.__pde_p.setSymmetryOn()
108 gross 1639
109 gross 2100 self.setTolerance()
110     self.setSmall()
111 gross 1639 def useJaumannStress(self):
112     """
113     return C{True} if Jaumann stress is included.
114     """
115     return self.__useJaumannStress
116 gross 2100
117 gross 1639 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 gross 2100
125 gross 1639 def getSmall(self):
126     """
127     returns small value
128     @rtype: positive float
129     """
130     return self.__small
131 gross 2100
132     def setTolerance(self,tol=1.e-4):
133 gross 1639 """
134 gross 2100 sets the tolerance
135 gross 1639 """
136 gross 2100 self.__pde_v.setTolerance(tol**2)
137     self.__pde_p.setTolerance(tol**2)
138     self.__tol=tol
139    
140     def getTolerance(self):
141 gross 1639 """
142 gross 2100 returns the set tolerance
143     @rtype: positive float
144 gross 1639 """
145 gross 2100 return self.__tol
146    
147     def getDomain(self):
148 gross 1639 """
149 gross 2100 returns the domain
150 gross 1639 """
151 gross 2100 return self.__domain
152 gross 1639 def getTime(self):
153     """
154     returns current time
155     """
156     return self.__t
157 gross 2100
158     def setPowerLaw(self,id,eta_N, tau_t=None, power=1):
159 gross 1639 """
160 gross 2100 sets the power-law parameters for material q
161 gross 1639 """
162 gross 2100 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 gross 1639
168 gross 2100 def setPowerLaws(self,eta_N, tau_t, power):
169 gross 1639 """
170 gross 2100 set the parameters of the powerlaw for all materials
171 gross 1639 """
172 gross 2100 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 gross 1639 """
178 gross 2100 set the parameters for the Drucker-Prager model
179 gross 1639 """
180 gross 2100 self.__tau_Y=tau_Y
181     self.__friction=friction
182 gross 1639
183 gross 2100 def setElasticShearModulus(self,mu=None):
184 gross 1639 """
185 gross 2100 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 gross 1639
192     @param F: external force.
193 gross 2100 @param f: surface force
194 gross 1639 @param q: location of constraints
195     """
196 gross 2100 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 gross 1639
201 gross 2100 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 gross 1639 else:
215 gross 2100 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 gross 1639 else:
222 gross 2100 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 gross 1639
239 gross 2100 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 gross 1639 """
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 gross 2100 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 gross 1639
324 gross 2100 #=========================================================================================
325 gross 1639
326 gross 2100 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 gross 1639
332 gross 2100 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 gross 1659
344 gross 2100 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 gross 1639
355 gross 2100 def getTau(self,stress=None):
356 gross 1639 """
357 gross 2100 returns current second stress deviatoric invariant
358 gross 1639 """
359 gross 2100 if stress==None: stress=self.getDeviatoricStress()
360     return util.sqrt(0.5)*util.length(stress)
361 gross 1639
362 gross 2100 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 gross 1639

  ViewVC Help
Powered by ViewVC 1.1.26