/[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 2386 - (show annotations)
Wed Apr 15 03:54:25 2009 UTC (10 years, 4 months ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 19964 byte(s)
some minor fixes
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__="https://launchpad.net/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, NegativeNorm, GMRES
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, weight=None, useReduced=False):
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 if weight == None:
59 s=self.domain.getSize()
60 self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
61 else:
62 self.__l=weight
63 self.__pde_v=LinearPDESystem(domain)
64 if useReduced: self.__pde_v.setReducedOrderOn()
65 self.__pde_v.setSymmetryOn()
66 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
67 # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
68 self.__pde_p=LinearSinglePDE(domain)
69 self.__pde_p.setSymmetryOn()
70 if useReduced: self.__pde_p.setReducedOrderOn()
71 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72 self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
73 self.setTolerance()
74 self.setAbsoluteTolerance()
75 self.setSubProblemTolerance()
76
77 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
78 """
79 assigns values to model parameters
80
81 @param f: volumetic sources/sinks
82 @type f: scalar value on the domain (e.g. L{Data})
83 @param g: flux sources/sinks
84 @type g: vector values on the domain (e.g. L{Data})
85 @param location_of_fixed_pressure: mask for locations where pressure is fixed
86 @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
87 @param location_of_fixed_flux: mask for locations where flux is fixed.
88 @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
89 @param permeability: permeability tensor. If scalar C{s} is given the tensor with
90 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
91 C{v} on the main diagonal is used.
92 @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
93
94 @note: the values of parameters which are not set by calling C{setValue} are not altered.
95 @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
96 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
97 is along the M{x_i} axis.
98 """
99 if f !=None:
100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101 if f.isEmpty():
102 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103 else:
104 if f.getRank()>0: raise ValueError,"illegal rank of f."
105 self.__f=f
106 if g !=None:
107 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108 if g.isEmpty():
109 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
110 else:
111 if not g.getShape()==(self.domain.getDim(),):
112 raise ValueError,"illegal shape of g"
113 self.__g=g
114
115 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
116 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
117
118 if permeability!=None:
119 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
120 if perm.getRank()==0:
121 perm=perm*util.kronecker(self.domain.getDim())
122 elif perm.getRank()==1:
123 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
124 for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
125 elif perm.getRank()==2:
126 pass
127 else:
128 raise ValueError,"illegal rank of permeability."
129 self.__permeability=perm
130 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
131
132 def setTolerance(self,rtol=1e-4):
133 """
134 sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
135
136 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
137
138 where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
139
140 @param rtol: relative tolerance for the pressure
141 @type rtol: non-negative C{float}
142 """
143 if rtol<0:
144 raise ValueError,"Relative tolerance needs to be non-negative."
145 self.__rtol=rtol
146 def getTolerance(self):
147 """
148 returns the relative tolerance
149
150 @return: current relative tolerance
151 @rtype: C{float}
152 """
153 return self.__rtol
154
155 def setAbsoluteTolerance(self,atol=0.):
156 """
157 sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
158
159 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
160
161 where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
162
163 @param atol: absolute tolerance for the pressure
164 @type atol: non-negative C{float}
165 """
166 if atol<0:
167 raise ValueError,"Absolute tolerance needs to be non-negative."
168 self.__atol=atol
169 def getAbsoluteTolerance(self):
170 """
171 returns the absolute tolerance
172
173 @return: current absolute tolerance
174 @rtype: C{float}
175 """
176 return self.__atol
177
178 def setSubProblemTolerance(self,rtol=None):
179 """
180 Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
181 C{self.getTolerance()**2} is used.
182
183 @param rtol: relative tolerence
184 @type rtol: positive C{float}
185 """
186 if rtol == None:
187 if self.getTolerance()<=0.:
188 raise ValueError,"A positive relative tolerance must be set."
189 self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
190 else:
191 if rtol<=0:
192 raise ValueError,"sub-problem tolerance must be positive."
193 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
194
195 def getSubProblemTolerance(self):
196 """
197 Returns the subproblem reduction factor.
198
199 @return: subproblem reduction factor
200 @rtype: C{float}
201 """
202 return self.__sub_tol
203
204 def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
205 """
206 solves the problem.
207
208 The iteration is terminated if the residual norm is less then self.getTolerance().
209
210 @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.
211 @type u0: vector value on the domain (e.g. L{Data}).
212 @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.
213 @type p0: scalar value on the domain (e.g. L{Data}).
214 @param verbose: if set some information on iteration progress are printed
215 @type verbose: C{bool}
216 @param show_details: if set information on the subiteration process are printed.
217 @type show_details: C{bool}
218 @return: flux and pressure
219 @rtype: C{tuple} of L{Data}.
220
221 @note: The problem is solved as a least squares form
222
223 M{(I+D^*D)u+Qp=D^*f+g}
224 M{Q^*u+Q^*Qp=Q^*g}
225
226 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
227 We eliminate the flux form the problem by setting
228
229 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
230
231 form the first equation. Inserted into the second equation we get
232
233 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
234
235 which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
236 PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
237 """
238 self.verbose=verbose
239 self.show_details= show_details and self.verbose
240 rtol=self.getTolerance()
241 atol=self.getAbsoluteTolerance()
242 if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243
244 num_corrections=0
245 converged=False
246 p=p0
247 norm_r=None
248 while not converged:
249 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250 Qp=self.__Q(p)
251 norm_v=self.__L2(v)
252 norm_Qp=self.__L2(Qp)
253 if norm_v == 0.:
254 if norm_Qp == 0.:
255 return v,p
256 else:
257 fac=norm_Qp
258 else:
259 if norm_Qp == 0.:
260 fac=norm_v
261 else:
262 fac=2./(1./norm_v+1./norm_Qp)
263 ATOL=(atol+rtol*fac)
264 if self.verbose:
265 print "DarcyFlux: L2 norm of v = %e."%norm_v
266 print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
267 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
268 if norm_r == None or norm_r>ATOL:
269 if num_corrections>max_num_corrections:
270 raise ValueError,"maximum number of correction steps reached."
271 p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
272 num_corrections+=1
273 else:
274 converged=True
275 return v,p
276 def __L2(self,v):
277 return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
278
279 def __Q(self,p):
280 return util.tensor_mult(self.__permeability,util.grad(p))
281
282 def __Aprod(self,dp):
283 self.__pde_v.setTolerance(self.getSubProblemTolerance())
284 if self.show_details: print "DarcyFlux: Applying operator"
285 Qdp=self.__Q(dp)
286 self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
287 du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
288 # self.__pde_v.getOperator().saveMM("proj.mm")
289 return Qdp+du
290 def __inner_GMRES(self,r,s):
291 return util.integrate(util.inner(r,s))
292
293 def __inner_PCG(self,p,r):
294 return util.integrate(util.inner(self.__Q(p), r))
295
296 def __Msolve_PCG(self,r):
297 self.__pde_p.setTolerance(self.getSubProblemTolerance())
298 if self.show_details: print "DarcyFlux: Applying preconditioner"
299 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
300 # self.__pde_p.getOperator().saveMM("prec.mm")
301 return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
302
303 def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
304 """
305 returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
306 on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
307 Note that C{g} and C{f} are used, see L{setValue}.
308
309 @param p: pressure.
310 @type p: scalar value on the domain (e.g. L{Data}).
311 @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
312 @type fixed_flux: vector values on the domain (e.g. L{Data}).
313 @param tol: relative tolerance to be used.
314 @type tol: positive C{float}.
315 @return: flux
316 @rtype: L{Data}
317 @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}}
318 for the permeability M{k_{ij}}
319 """
320 self.__pde_v.setTolerance(self.getSubProblemTolerance())
321 g=self.__g
322 f=self.__f
323 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
324 if p == None:
325 self.__pde_v.setValue(Y=g)
326 else:
327 self.__pde_v.setValue(Y=g-self.__Q(p))
328 return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
329
330 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
331 """
332 solves
333
334 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
335 u_{i,i}=0
336
337 u=0 where fixed_u_mask>0
338 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
339
340 if surface_stress is not given 0 is assumed.
341
342 typical usage:
343
344 sp=StokesProblemCartesian(domain)
345 sp.setTolerance()
346 sp.initialize(...)
347 v,p=sp.solve(v0,p0)
348 """
349 def __init__(self,domain,**kwargs):
350 """
351 initialize the Stokes Problem
352
353 @param domain: domain of the problem. The approximation order needs to be two.
354 @type domain: L{Domain}
355 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
356 """
357 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
358 self.domain=domain
359 self.vol=util.integrate(1.,Function(self.domain))
360 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
361 self.__pde_u.setSymmetryOn()
362 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
363 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
364
365 self.__pde_prec=LinearPDE(domain)
366 self.__pde_prec.setReducedOrderOn()
367 # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
368 self.__pde_prec.setSymmetryOn()
369
370 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
371 """
372 assigns values to the model parameters
373
374 @param f: external force
375 @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
376 @param fixed_u_mask: mask of locations with fixed velocity.
377 @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
378 @param eta: viscosity
379 @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
380 @param surface_stress: normal surface stress
381 @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
382 @param stress: initial stress
383 @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
384 @note: All values needs to be set.
385
386 """
387 self.eta=eta
388 A =self.__pde_u.createCoefficient("A")
389 self.__pde_u.setValue(A=Data())
390 for i in range(self.domain.getDim()):
391 for j in range(self.domain.getDim()):
392 A[i,j,j,i] += 1.
393 A[i,j,i,j] += 1.
394 self.__pde_prec.setValue(D=1/self.eta)
395 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
396 self.__f=f
397 self.__surface_stress=surface_stress
398 self.__stress=stress
399
400 def inner_pBv(self,p,v):
401 """
402 returns inner product of element p and div(v)
403
404 @param p: a pressure increment
405 @param v: a residual
406 @return: inner product of element p and div(v)
407 @rtype: C{float}
408 """
409 return util.integrate(-p*util.div(v))
410
411 def inner_p(self,p0,p1):
412 """
413 Returns inner product of p0 and p1
414
415 @param p0: a pressure
416 @param p1: a pressure
417 @return: inner product of p0 and p1
418 @rtype: C{float}
419 """
420 s0=util.interpolate(p0/self.eta,Function(self.domain))
421 s1=util.interpolate(p1/self.eta,Function(self.domain))
422 return util.integrate(s0*s1)
423
424 def norm_v(self,v):
425 """
426 returns the norm of v
427
428 @param v: a velovity
429 @return: norm of v
430 @rtype: non-negative C{float}
431 """
432 return util.sqrt(util.integrate(util.length(util.grad(v))))
433
434 def getV(self, p, v0):
435 """
436 return the value for v for a given p (overwrite)
437
438 @param p: a pressure
439 @param v0: a initial guess for the value v to return.
440 @return: v given as M{v= A^{-1} (f-B^*p)}
441 """
442 self.__pde_u.setTolerance(self.getSubProblemTolerance())
443 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
444 if self.__stress.isEmpty():
445 self.__pde_u.setValue(X=p*util.kronecker(self.domain))
446 else:
447 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
448 out=self.__pde_u.getSolution(verbose=self.show_details)
449 return out
450
451
452 raise NotImplementedError,"no v calculation implemented."
453
454
455 def norm_Bv(self,v):
456 """
457 Returns Bv (overwrite).
458
459 @rtype: equal to the type of p
460 @note: boundary conditions on p should be zero!
461 """
462 return util.sqrt(util.integrate(util.div(v)**2))
463
464 def solve_AinvBt(self,p):
465 """
466 Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
467
468 @param p: a pressure increment
469 @return: the solution of M{Av=B^*p}
470 @note: boundary conditions on v should be zero!
471 """
472 self.__pde_u.setTolerance(self.getSubProblemTolerance())
473 self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
474 out=self.__pde_u.getSolution(verbose=self.show_details)
475 return out
476
477 def solve_precB(self,v):
478 """
479 applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
480 with accuracy L{self.getSubProblemTolerance()} (overwrite).
481
482 @param v: velocity increment
483 @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
484 @note: boundary conditions on p are zero.
485 """
486 self.__pde_prec.setValue(Y=-util.div(v))
487 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
488 return self.__pde_prec.getSolution(verbose=self.show_details)

  ViewVC Help
Powered by ViewVC 1.1.26