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

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

  ViewVC Help
Powered by ViewVC 1.1.26