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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1673 - (hide annotations)
Thu Jul 24 22:28:50 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 9782 byte(s)
more algorithms for level set
1 gross 1414 # $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 artak 1550 from pdetools import HomogeneousSaddlePointProblem,Projector
40 gross 1414
41 gross 1659 class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):
42     """
43     solves
44    
45     -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i
46     u_{i,i}=0
47    
48     u=0 where fixed_u_mask>0
49     eta*(u_{i,j}+u_{j,i})*n_j=surface_stress
50    
51     if surface_stress is not give 0 is assumed.
52    
53     typical usage:
54    
55     sp=StokesProblemCartesian(domain)
56     sp.setTolerance()
57     sp.initialize(...)
58     v,p=sp.solve(v0,p0)
59     """
60     def __init__(self,domain,**kwargs):
61     HomogeneousSaddlePointProblem.__init__(self,**kwargs)
62     self.domain=domain
63     self.vol=util.integrate(1.,Function(self.domain))
64     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
65     self.__pde_u.setSymmetryOn()
66 gross 1661 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
67 gross 1659
68 gross 1661 # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)
69     # self.__pde_proj.setReducedOrderOn()
70     # self.__pde_proj.setSymmetryOn()
71 gross 1659 # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)
72    
73     def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):
74     self.eta=eta
75     A =self.__pde_u.createCoefficientOfGeneralPDE("A")
76     self.__pde_u.setValue(A=Data())
77     for i in range(self.domain.getDim()):
78     for j in range(self.domain.getDim()):
79     A[i,j,j,i] += 1.
80     A[i,j,i,j] += 1.
81     # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))
82     self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
83    
84 gross 1661 # self.__pde_proj.setValue(D=1/eta)
85     # self.__pde_proj.setValue(Y=1.)
86     # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))
87     self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))
88 gross 1659
89     def B(self,arg):
90     a=util.div(arg, ReducedFunction(self.domain))
91     return a-util.integrate(a)/self.vol
92    
93     def inner(self,p0,p1):
94     return util.integrate(p0*p1)
95 gross 1673 # return util.integrate(1/self.__inv_eta**2*p0*p1)
96 gross 1659
97     def getStress(self,u):
98     mg=util.grad(u)
99     return 2.*self.eta*util.symmetric(mg)
100     def getEtaEffective(self):
101     return self.eta
102    
103     def solve_A(self,u,p):
104     """
105     solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
106     """
107     self.__pde_u.setTolerance(self.getSubProblemTolerance())
108     self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))
109     return self.__pde_u.getSolution(verbose=self.show_details)
110    
111    
112     def solve_prec(self,p):
113     a=self.__inv_eta*p
114     return a-util.integrate(a)/self.vol
115    
116     def stoppingcriterium(self,Bv,v,p):
117     n_r=util.sqrt(self.inner(Bv,Bv))
118     n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
119     if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)
120     if self.iter == 0: self.__n_v=n_v;
121     self.__n_v, n_v_old =n_v, self.__n_v
122     self.iter+=1
123     if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
124     if self.verbose: print "PCG terminated after %s steps."%self.iter
125     return True
126     else:
127     return False
128 gross 1673 def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
129     if TOL==None:
130     TOL=self.getTolerance()
131     if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
132     self.iter+=1
133    
134     if norm_r <= norm_b*TOL:
135     if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
136     return True
137     else:
138     return False
139 gross 1659
140    
141 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
142     """
143     solves
144    
145     -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i
146     u_{i,i}=0
147    
148     u=0 where fixed_u_mask>0
149     eta*(u_{i,j}+u_{j,i})*n_j=surface_stress
150    
151     if surface_stress is not give 0 is assumed.
152    
153     typical usage:
154    
155     sp=StokesProblemCartesian(domain)
156     sp.setTolerance()
157     sp.initialize(...)
158     v,p=sp.solve(v0,p0)
159     """
160     def __init__(self,domain,**kwargs):
161     HomogeneousSaddlePointProblem.__init__(self,**kwargs)
162     self.domain=domain
163     self.vol=util.integrate(1.,Function(self.domain))
164     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
165     self.__pde_u.setSymmetryOn()
166 gross 1639 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
167 gross 1414
168     self.__pde_prec=LinearPDE(domain)
169     self.__pde_prec.setReducedOrderOn()
170     self.__pde_prec.setSymmetryOn()
171    
172     self.__pde_proj=LinearPDE(domain)
173     self.__pde_proj.setReducedOrderOn()
174     self.__pde_proj.setSymmetryOn()
175     self.__pde_proj.setValue(D=1.)
176    
177     def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):
178     self.eta=eta
179     A =self.__pde_u.createCoefficientOfGeneralPDE("A")
180     self.__pde_u.setValue(A=Data())
181     for i in range(self.domain.getDim()):
182     for j in range(self.domain.getDim()):
183     A[i,j,j,i] += 1.
184     A[i,j,i,j] += 1.
185 artak 1554 self.__pde_prec.setValue(D=1/self.eta)
186 gross 1414 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
187    
188     def B(self,arg):
189     d=util.div(arg)
190     self.__pde_proj.setValue(Y=d)
191     self.__pde_proj.setTolerance(self.getSubProblemTolerance())
192     return self.__pde_proj.getSolution(verbose=self.show_details)
193    
194     def inner(self,p0,p1):
195     s0=util.interpolate(p0,Function(self.domain))
196     s1=util.interpolate(p1,Function(self.domain))
197     return util.integrate(s0*s1)
198    
199 artak 1550 def inner_a(self,a0,a1):
200     p0=util.interpolate(a0[1],Function(self.domain))
201     p1=util.interpolate(a1[1],Function(self.domain))
202     alfa=(1/self.vol)*util.integrate(p0)
203     beta=(1/self.vol)*util.integrate(p1)
204     v0=util.grad(a0[0])
205     v1=util.grad(a1[0])
206     return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))
207    
208    
209 gross 1414 def getStress(self,u):
210     mg=util.grad(u)
211     return 2.*self.eta*util.symmetric(mg)
212 gross 1639 def getEtaEffective(self):
213     return self.eta
214 gross 1414
215     def solve_A(self,u,p):
216     """
217     solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
218     """
219     self.__pde_u.setTolerance(self.getSubProblemTolerance())
220 gross 1470 self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))
221 gross 1414 return self.__pde_u.getSolution(verbose=self.show_details)
222    
223 artak 1550
224 gross 1414 def solve_prec(self,p):
225 artak 1550 #proj=Projector(domain=self.domain, reduce = True, fast=False)
226 gross 1414 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
227     self.__pde_prec.setValue(Y=p)
228     q=self.__pde_prec.getSolution(verbose=self.show_details)
229 artak 1554 return q
230    
231     def solve_prec1(self,p):
232     #proj=Projector(domain=self.domain, reduce = True, fast=False)
233     self.__pde_prec.setTolerance(self.getSubProblemTolerance())
234     self.__pde_prec.setValue(Y=p)
235     q=self.__pde_prec.getSolution(verbose=self.show_details)
236 artak 1550 q0=util.interpolate(q,Function(self.domain))
237 gross 1562 print util.inf(q*q0),util.sup(q*q0)
238 artak 1550 q-=(1/self.vol)*util.integrate(q0)
239 gross 1562 print util.inf(q*q0),util.sup(q*q0)
240 gross 1414 return q
241 artak 1550
242 gross 1414 def stoppingcriterium(self,Bv,v,p):
243     n_r=util.sqrt(self.inner(Bv,Bv))
244 gross 1552 n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
245     if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)
246     if self.iter == 0: self.__n_v=n_v;
247     self.__n_v, n_v_old =n_v, self.__n_v
248 gross 1414 self.iter+=1
249 gross 1552 if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
250 gross 1414 if self.verbose: print "PCG terminated after %s steps."%self.iter
251     return True
252     else:
253     return False
254 artak 1519 def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
255     if TOL==None:
256     TOL=self.getTolerance()
257     if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
258 artak 1465 self.iter+=1
259 artak 1519
260     if norm_r <= norm_b*TOL:
261 artak 1517 if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
262 artak 1465 return True
263     else:
264     return False
265 artak 1481
266 artak 1517

  ViewVC Help
Powered by ViewVC 1.1.26