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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1809 - (show 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
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
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 __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 """
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 from pdetools import HomogeneousSaddlePointProblem,Projector
39
40 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 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
66
67 # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)
68 # self.__pde_proj.setReducedOrderOn()
69 # self.__pde_proj.setSymmetryOn()
70 # 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 # 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
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 # return util.integrate(1/self.__inv_eta**2*p0*p1)
95
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 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
139
140 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 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
166
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 self.__pde_prec.setValue(D=1/self.eta)
185 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 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 def getStress(self,u):
209 mg=util.grad(u)
210 return 2.*self.eta*util.symmetric(mg)
211 def getEtaEffective(self):
212 return self.eta
213
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 self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))
220 return self.__pde_u.getSolution(verbose=self.show_details)
221
222
223 def solve_prec(self,p):
224 #proj=Projector(domain=self.domain, reduce = True, fast=False)
225 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
226 self.__pde_prec.setValue(Y=p)
227 q=self.__pde_prec.getSolution(verbose=self.show_details)
228 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 q0=util.interpolate(q,Function(self.domain))
236 print util.inf(q*q0),util.sup(q*q0)
237 q-=(1/self.vol)*util.integrate(q0)
238 print util.inf(q*q0),util.sup(q*q0)
239 return q
240
241 def stoppingcriterium(self,Bv,v,p):
242 n_r=util.sqrt(self.inner(Bv,Bv))
243 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 self.iter+=1
248 if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
249 if self.verbose: print "PCG terminated after %s steps."%self.iter
250 return True
251 else:
252 return False
253 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 self.iter+=1
258
259 if norm_r <= norm_b*TOL:
260 if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
261 return True
262 else:
263 return False
264
265

  ViewVC Help
Powered by ViewVC 1.1.26