/[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 2123 - (show annotations)
Wed Dec 3 03:26:02 2008 UTC (10 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 15979 byte(s)
typo fixed
1 ########################################################
2 #
3 # Copyright (c) 2003-2008 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 ########################################################
12
13 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14 Earth Systems Science Computational Center (ESSCC)
15 http://www.uq.edu.au/esscc
16 Primary Business: Queensland, Australia"""
17 __license__="""Licensed under the Open Software License version 3.0
18 http://www.opensource.org/licenses/osl-3.0.php"""
19 __url__="http://www.uq.edu.au/esscc/escript-finley"
20
21 """
22 Some models for flow
23
24 @var __author__: name of author
25 @var __copyright__: copyrights
26 @var __license__: licence agreement
27 @var __url__: url entry point on documentation
28 @var __version__: version
29 @var __date__: date of the version
30 """
31
32 __author__="Lutz Gross, l.gross@uq.edu.au"
33
34 from escript import *
35 import util
36 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG
38
39 class DarcyFlow(object):
40 """
41 solves the problem
42
43 M{u_i+k_{ij}*p_{,j} = g_i}
44 M{u_{i,i} = f}
45
46 where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47
48 @note: The problem is solved in a least squares formulation.
49 """
50
51 def __init__(self, domain):
52 """
53 initializes the Darcy flux problem
54 @param domain: domain of the problem
55 @type domain: L{Domain}
56 """
57 self.domain=domain
58 self.__pde_v=LinearPDESystem(domain)
59 self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
60 self.__pde_v.setSymmetryOn()
61 self.__pde_p=LinearSinglePDE(domain)
62 self.__pde_p.setSymmetryOn()
63 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
64 self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
65
66 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
67 """
68 assigns values to model parameters
69
70 @param f: volumetic sources/sinks
71 @type f: scalar value on the domain (e.g. L{Data})
72 @param g: flux sources/sinks
73 @type g: vector values on the domain (e.g. L{Data})
74 @param location_of_fixed_pressure: mask for locations where pressure is fixed
75 @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
76 @param location_of_fixed_flux: mask for locations where flux is fixed.
77 @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
78 @param permeability: permeability tensor. If scalar C{s} is given the tensor with
79 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
80 C{v} on the main diagonal is used.
81 @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
82
83 @note: the values of parameters which are not set by calling C{setValue} are not altered.
84 @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
85 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
86 is along the M{x_i} axis.
87 """
88 if f !=None:
89 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
90 if f.isEmpty():
91 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
92 else:
93 if f.getRank()>0: raise ValueError,"illegal rank of f."
94 self.f=f
95 if g !=None:
96 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
97 if g.isEmpty():
98 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
99 else:
100 if not g.getShape()==(self.domain.getDim(),):
101 raise ValueError,"illegal shape of g"
102 self.__g=g
103
104 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
105 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
106
107 if permeability!=None:
108 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
109 if perm.getRank()==0:
110 perm=perm*util.kronecker(self.domain.getDim())
111 elif perm.getRank()==1:
112 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
113 for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
114 elif perm.getRank()==2:
115 pass
116 else:
117 raise ValueError,"illegal rank of permeability."
118 self.__permeability=perm
119 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
120
121
122 def getFlux(self,p, fixed_flux=Data(),tol=1.e-8):
123 """
124 returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
125 on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
126 Note that C{g} and C{f} are used, L{setValue}.
127
128 @param p: pressure.
129 @type p: scalar value on the domain (e.g. L{Data}).
130 @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
131 @type fixed_flux: vector values on the domain (e.g. L{Data}).
132 @param tol: relative tolerance to be used.
133 @type tol: positive float.
134 @return: flux
135 @rtype: L{Data}
136 @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
137 for the permeability M{k_{ij}}
138 """
139 self.__pde_v.setTolerance(tol)
140 self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
141 return self.__pde_v.getSolution()
142
143 def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
144 """
145 solves the problem.
146
147 The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where
148 C{|q|} denotes the norm of the right hand side (see escript user's guide for details).
149
150 @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
151 @type u0: vector value on the domain (e.g. L{Data}).
152 @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
153 @type p0: scalar value on the domain (e.g. L{Data}).
154 @param atol: absolute tolerance for the pressure
155 @type atol: non-negative C{float}
156 @param rtol: relative tolerance for the pressure
157 @type rtol: non-negative C{float}
158 @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
159 @type sub_rtol: positive-negative C{float}
160 @param verbose: if set some information on iteration progress are printed
161 @type verbose: C{bool}
162 @param show_details: if set information on the subiteration process are printed.
163 @type show_details: C{bool}
164 @return: flux and pressure
165 @rtype: C{tuple} of L{Data}.
166
167 @note: The problem is solved as a least squares form
168
169 M{(I+D^*D)u+Qp=D^*f+g}
170 M{Q^*u+Q^*Qp=Q^*g}
171
172 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
173 We eliminate the flux form the problem by setting
174
175 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
176
177 form the first equation. Inserted into the second equation we get
178
179 M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0 on location_of_fixed_pressure
180
181 which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
182 PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
183 """
184 self.verbose=verbose
185 self.show_details= show_details and self.verbose
186 self.__pde_v.setTolerance(sub_rtol)
187 self.__pde_p.setTolerance(sub_rtol)
188 p2=p0*self.__pde_p.getCoefficient("q")
189 u2=u0*self.__pde_v.getCoefficient("q")
190 g=self.__g-u2-util.tensor_mult(self.__permeability,util.grad(p2))
191 f=self.__f-util.div(u2)
192 self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())
193 dv=self.__pde_v.getSolution(verbose=show_details)
194 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))
195 self.__pde_p.setValue(r=Data())
196 dp=self.__pde_p.getSolution(verbose=self.show_details)
197 norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))
198 if norm_rhs<0:
199 raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."
200 ATOL=util.sqrt(norm_rhs)*rtol +atol
201 if not ATOL>0:
202 raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)
203 rhs=ArithmeticTuple(g,dv)
204 dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0-p2, verbose=self.verbose, initial_guess=True)
205 return u2+r[1],p2+dp
206
207 def __Aprod_PCG(self,p):
208 if self.show_details: print "DarcyFlux: Applying operator"
209 Qp=util.tensor_mult(self.__permeability,util.grad(p))
210 self.__pde_v.setValue(Y=Qp,X=Data())
211 w=self.__pde_v.getSolution(verbose=self.show_details)
212 return ArithmeticTuple(Qp,w)
213
214 def __inner_PCG(self,p,r):
215 a=util.tensor_mult(self.__permeability,util.grad(p))
216 return util.integrate(util.inner(a,r[0]-r[1]))
217
218 def __Msolve_PCG(self,r):
219 if self.show_details: print "DarcyFlux: Applying preconditioner"
220 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]))
221 return self.__pde_p.getSolution(verbose=self.show_details)
222
223 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
224 """
225 solves
226
227 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
228 u_{i,i}=0
229
230 u=0 where fixed_u_mask>0
231 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
232
233 if surface_stress is not given 0 is assumed.
234
235 typical usage:
236
237 sp=StokesProblemCartesian(domain)
238 sp.setTolerance()
239 sp.initialize(...)
240 v,p=sp.solve(v0,p0)
241 """
242 def __init__(self,domain,**kwargs):
243 """
244 initialize the Stokes Problem
245
246 @param domain: domain of the problem. The approximation order needs to be two.
247 @type domain: L{Domain}
248 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
249 """
250 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
251 self.domain=domain
252 self.vol=util.integrate(1.,Function(self.domain))
253 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
254 self.__pde_u.setSymmetryOn()
255 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
256 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
257
258 self.__pde_prec=LinearPDE(domain)
259 self.__pde_prec.setReducedOrderOn()
260 self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
261 self.__pde_prec.setSymmetryOn()
262
263 self.__pde_proj=LinearPDE(domain)
264 self.__pde_proj.setReducedOrderOn()
265 self.__pde_proj.setSymmetryOn()
266 self.__pde_proj.setValue(D=1.)
267
268 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
269 """
270 assigns values to the model parameters
271
272 @param f: external force
273 @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
274 @param fixed_u_mask: mask of locations with fixed velocity.
275 @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
276 @param eta: viscosity
277 @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
278 @param surface_stress: normal surface stress
279 @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
280 @param stress: initial stress
281 @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
282 @note: All values needs to be set.
283
284 """
285 self.eta=eta
286 A =self.__pde_u.createCoefficient("A")
287 self.__pde_u.setValue(A=Data())
288 for i in range(self.domain.getDim()):
289 for j in range(self.domain.getDim()):
290 A[i,j,j,i] += 1.
291 A[i,j,i,j] += 1.
292 self.__pde_prec.setValue(D=1/self.eta)
293 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
294 self.__stress=stress
295
296 def B(self,v):
297 """
298 returns div(v)
299 @rtype: equal to the type of p
300
301 @note: boundary conditions on p should be zero!
302 """
303 if self.show_details: print "apply divergence:"
304 self.__pde_proj.setValue(Y=-util.div(v))
305 self.__pde_proj.setTolerance(self.getSubProblemTolerance())
306 return self.__pde_proj.getSolution(verbose=self.show_details)
307
308 def inner_pBv(self,p,Bv):
309 """
310 returns inner product of element p and Bv (overwrite)
311
312 @type p: equal to the type of p
313 @type Bv: equal to the type of result of operator B
314 @rtype: C{float}
315
316 @rtype: equal to the type of p
317 """
318 s0=util.interpolate(p,Function(self.domain))
319 s1=util.interpolate(Bv,Function(self.domain))
320 return util.integrate(s0*s1)
321
322 def inner_p(self,p0,p1):
323 """
324 returns inner product of element p0 and p1 (overwrite)
325
326 @type p0: equal to the type of p
327 @type p1: equal to the type of p
328 @rtype: C{float}
329
330 @rtype: equal to the type of p
331 """
332 s0=util.interpolate(p0/self.eta,Function(self.domain))
333 s1=util.interpolate(p1/self.eta,Function(self.domain))
334 return util.integrate(s0*s1)
335
336 def inner_v(self,v0,v1):
337 """
338 returns inner product of two element v0 and v1 (overwrite)
339
340 @type v0: equal to the type of v
341 @type v1: equal to the type of v
342 @rtype: C{float}
343
344 @rtype: equal to the type of v
345 """
346 gv0=util.grad(v0)
347 gv1=util.grad(v1)
348 return util.integrate(util.inner(gv0,gv1))
349
350 def solve_A(self,u,p):
351 """
352 solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
353 """
354 if self.show_details: print "solve for velocity:"
355 self.__pde_u.setTolerance(self.getSubProblemTolerance())
356 if self.__stress.isEmpty():
357 self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
358 else:
359 self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
360 out=self.__pde_u.getSolution(verbose=self.show_details)
361 return out
362
363 def solve_prec(self,p):
364 if self.show_details: print "apply preconditioner:"
365 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
366 self.__pde_prec.setValue(Y=p)
367 q=self.__pde_prec.getSolution(verbose=self.show_details)
368 return q

  ViewVC Help
Powered by ViewVC 1.1.26