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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1673 - (show annotations)
Thu Jul 24 22:28:50 2008 UTC (11 years, 1 month ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 9782 byte(s)
more algorithms for level set
1 # $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 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 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
67
68 # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)
69 # self.__pde_proj.setReducedOrderOn()
70 # self.__pde_proj.setSymmetryOn()
71 # 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 # 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
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 # return util.integrate(1/self.__inv_eta**2*p0*p1)
96
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 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
140
141 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 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
167
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 self.__pde_prec.setValue(D=1/self.eta)
186 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 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 def getStress(self,u):
210 mg=util.grad(u)
211 return 2.*self.eta*util.symmetric(mg)
212 def getEtaEffective(self):
213 return self.eta
214
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 self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))
221 return self.__pde_u.getSolution(verbose=self.show_details)
222
223
224 def solve_prec(self,p):
225 #proj=Projector(domain=self.domain, reduce = True, fast=False)
226 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
227 self.__pde_prec.setValue(Y=p)
228 q=self.__pde_prec.getSolution(verbose=self.show_details)
229 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 q0=util.interpolate(q,Function(self.domain))
237 print util.inf(q*q0),util.sup(q*q0)
238 q-=(1/self.vol)*util.integrate(q0)
239 print util.inf(q*q0),util.sup(q*q0)
240 return q
241
242 def stoppingcriterium(self,Bv,v,p):
243 n_r=util.sqrt(self.inner(Bv,Bv))
244 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 self.iter+=1
249 if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
250 if self.verbose: print "PCG terminated after %s steps."%self.iter
251 return True
252 else:
253 return False
254 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 self.iter+=1
259
260 if norm_r <= norm_b*TOL:
261 if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
262 return True
263 else:
264 return False
265
266

  ViewVC Help
Powered by ViewVC 1.1.26