/[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 2169 - (hide annotations)
Wed Dec 17 03:08:58 2008 UTC (10 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 18141 byte(s)
Assorted spelling, grammar, whitespace and copy/paste error fixes (Part 2).
All epydoc warnings for these files have been fixed.
This commit should be a no-op.

1 ksteube 1809 ########################################################
2 gross 1414 #
3 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
4     # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1414 #
7 ksteube 1809 # 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 gross 1414 #
11 ksteube 1809 ########################################################
12 gross 1414
13 ksteube 1809 __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 gross 1414 """
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 gross 2100 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37 gross 2156 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm
38 gross 1414
39 gross 2100 class DarcyFlow(object):
40     """
41 caltinay 2169 Represents and solves the problem
42 gross 1659
43 gross 2100 M{u_i+k_{ij}*p_{,j} = g_i}
44 caltinay 2169
45 gross 2100 M{u_{i,i} = f}
46 gross 1659
47 caltinay 2169 where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents
48     the permeability.
49 gross 1659
50 gross 2100 @note: The problem is solved in a least squares formulation.
51     """
52 gross 1659
53 gross 2100 def __init__(self, domain):
54     """
55 caltinay 2169 Initializes the Darcy flux problem.
56    
57 gross 2100 @param domain: domain of the problem
58     @type domain: L{Domain}
59     """
60     self.domain=domain
61     self.__pde_v=LinearPDESystem(domain)
62     self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
63     self.__pde_v.setSymmetryOn()
64     self.__pde_p=LinearSinglePDE(domain)
65     self.__pde_p.setSymmetryOn()
66     self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
67     self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
68 gross 1659
69 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
70     """
71 caltinay 2169 Assigns values to model parameters.
72 gross 1659
73 caltinay 2169 @param f: volumetric sources/sinks
74     @type f: scalar value on the domain, e.g. L{Data}
75 gross 2100 @param g: flux sources/sinks
76 caltinay 2169 @type g: vector value on the domain, e.g. L{Data}
77 gross 2100 @param location_of_fixed_pressure: mask for locations where pressure is fixed
78 caltinay 2169 @type location_of_fixed_pressure: scalar value on the domain, e.g. L{Data}
79     @param location_of_fixed_flux: mask for locations where flux is fixed
80     @type location_of_fixed_flux: vector value on the domain (e.g. L{Data})
81     @param permeability: permeability tensor. If scalar C{s} is given the
82     tensor with C{s} on the main diagonal is used. If
83     vector C{v} is given the tensor with C{v} on the
84     main diagonal is used.
85     @type permeability: scalar, vector or tensor values on the domain, e.g.
86     L{Data}
87 gross 1659
88 caltinay 2169 @note: the values of parameters which are not set by calling
89     C{setValue} are not altered
90     @note: at any point on the boundary of the domain the pressure
91     (C{location_of_fixed_pressure}) >0 or the normal component of
92     the flux (C{location_of_fixed_flux[i]}) >0 if the direction of
93     the normal is along the M{x_i} axis.
94 gross 2100 """
95 caltinay 2169 if f !=None:
96 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
97     if f.isEmpty():
98     f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
99     else:
100     if f.getRank()>0: raise ValueError,"illegal rank of f."
101     self.f=f
102 caltinay 2169 if g !=None:
103 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
104     if g.isEmpty():
105     g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
106     else:
107     if not g.getShape()==(self.domain.getDim(),):
108     raise ValueError,"illegal shape of g"
109     self.__g=g
110 gross 1659
111 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
112     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
113 gross 1659
114 gross 2100 if permeability!=None:
115     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
116     if perm.getRank()==0:
117     perm=perm*util.kronecker(self.domain.getDim())
118     elif perm.getRank()==1:
119     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
120     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
121     elif perm.getRank()==2:
122     pass
123     else:
124     raise ValueError,"illegal rank of permeability."
125     self.__permeability=perm
126     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
127 gross 1659
128    
129 gross 2156 def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
130 gross 2100 """
131 caltinay 2169 Returns the flux for a given pressure C{p}.
132    
133     The flux is equal to C{fixed_flux} on locations where
134     C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}
135     and C{f} are used.
136    
137     @param p: pressure
138     @type p: scalar value on the domain, e.g. L{Data}
139     @param fixed_flux: flux on the locations of the domain marked by
140     C{location_of_fixed_flux}
141     @type fixed_flux: vector values on the domain, e.g. L{Data}
142     @param tol: relative tolerance to be used
143     @type tol: positive float
144     @return: flux
145 gross 2100 @rtype: L{Data}
146 caltinay 2169 @note: the method uses the least squares solution
147     M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator
148     and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
149 gross 2100 """
150     self.__pde_v.setTolerance(tol)
151 gross 2123 self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
152 gross 2156 return self.__pde_v.getSolution(verbose=show_details)
153 gross 1659
154 caltinay 2169 def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
155     """
156     Solves the problem.
157 gross 1659
158 caltinay 2169 The iteration is terminated if the error in the pressure is less than
159     M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand
160     side (see escript user's guide for details).
161 gross 1659
162 caltinay 2169 @param u0: initial guess for the flux. At locations in the domain
163     marked by C{location_of_fixed_flux} the value of C{u0} is
164     kept unchanged.
165     @type u0: vector value on the domain, e.g. L{Data}
166     @param p0: initial guess for the pressure. At locations in the domain
167     marked by C{location_of_fixed_pressure} the value of C{p0}
168     is kept unchanged.
169     @type p0: scalar value on the domain, e.g. L{Data}
170     @param atol: absolute tolerance for the pressure
171     @type atol: non-negative C{float}
172     @param rtol: relative tolerance for the pressure
173     @type rtol: non-negative C{float}
174     @param sub_rtol: tolerance to be used in the sub iteration. It is
175     recommended that M{sub_rtol<rtol*5.e-3}
176     @type sub_rtol: positive-negative C{float}
177     @param verbose: if True information on iteration progress is printed
178     @type verbose: C{bool}
179     @param show_details: if True information on the sub-iteration process
180     is printed
181     @type show_details: C{bool}
182     @return: flux and pressure
183     @rtype: C{tuple} of L{Data}
184 gross 1659
185 caltinay 2169 @note: the problem is solved in a least squares formulation:
186 gross 1659
187 caltinay 2169 M{(I+D^*D)u+Qp=D^*f+g}
188 gross 2100
189 caltinay 2169 M{Q^*u+Q^*Qp=Q^*g}
190 gross 2100
191 caltinay 2169 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the
192     permeability M{k_{ij}}. We eliminate the flux from the problem by
193     setting
194    
195     M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}
196    
197     from the first equation. Inserted into the second equation we get
198    
199     M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}
200     on C{location_of_fixed_pressure}
201    
202     which is solved using the PCG method (precondition is M{Q^*Q}).
203     In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}
204     need to be solved using a sub-iteration scheme.
205     """
206     self.verbose=verbose
207     self.show_details= show_details and self.verbose
208     self.__pde_v.setTolerance(sub_rtol)
209     self.__pde_p.setTolerance(sub_rtol)
210     u2=u0*self.__pde_v.getCoefficient("q")
211     #
212     # first the reference velocity is calculated from
213     #
214     # (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
215     #
216     self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
217     u_ref=self.__pde_v.getSolution(verbose=show_details)
218     if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
219     self.__pde_v.setValue(r=Data())
220     #
221     # and then we calculate a reference pressure
222     #
223     # Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
224     #
225     self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
226     p_ref=self.__pde_p.getSolution(verbose=self.show_details)
227     if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
228     self.__pde_p.setValue(r=Data())
229     #
230     # (I+D^*D)du + Qdp = - Qp_ref u=du+u_ref
231     # Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0 p=dp+pref
232     #
233     # du= -(I+D^*D)^(-1} Q(p_ref+dp) u = u_ref+du
234     #
235     # => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
236     # or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
237     #
238     # r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
239     # with du=-(I+D^*D)^(-1} Q(p_ref+dp)
240     #
241     # we use the (du,Qdp) to represent the resudual
242     # Q^*Q is a preconditioner
243     #
244     # <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
245     #
246     Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
247     norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
248     ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
249     if not ATOL>0:
250     raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)
251     if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
252     #
253     # caclulate the initial residual
254     #
255     self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
256     du=self.__pde_v.getSolution(verbose=show_details)
257     r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
258     dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
259     util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
260     return u_ref+r[1],dp
261    
262 gross 2100 def __Aprod_PCG(self,p):
263 caltinay 2169 if self.show_details: print "DarcyFlux: Applying operator"
264     Qp=util.tensor_mult(self.__permeability,util.grad(p))
265     self.__pde_v.setValue(Y=Qp,X=Data())
266     w=self.__pde_v.getSolution(verbose=self.show_details)
267     return ArithmeticTuple(-Qp,w)
268 gross 2100
269     def __inner_PCG(self,p,r):
270 caltinay 2169 a=util.tensor_mult(self.__permeability,util.grad(p))
271     out=-util.integrate(util.inner(a,r[0]+r[1]))
272     return out
273 gross 2100
274     def __Msolve_PCG(self,r):
275 caltinay 2169 if self.show_details: print "DarcyFlux: Applying preconditioner"
276     self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))
277     return self.__pde_p.getSolution(verbose=self.show_details)
278 gross 2100
279 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
280     """
281 caltinay 2169 Represents and solves the problem
282 gross 1414
283 caltinay 2169 M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}
284 gross 1414
285 caltinay 2169 M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0
286 gross 1414
287 caltinay 2169 M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}
288 gross 1414
289 caltinay 2169 If C{surface_stress} is not given 0 is assumed.
290 gross 1414
291 caltinay 2169 Typical usage::
292    
293     sp = StokesProblemCartesian(domain)
294     sp.setTolerance()
295     sp.initialize(...)
296     v,p = sp.solve(v0,p0)
297 gross 1414 """
298     def __init__(self,domain,**kwargs):
299 gross 2100 """
300 caltinay 2169 Initializes the Stokes Problem.
301 gross 2100
302 caltinay 2169 @param domain: domain of the problem. The approximation order needs
303     to be two.
304 gross 2100 @type domain: L{Domain}
305 caltinay 2169 @warning: The approximation order needs to be two otherwise you may
306     see oscillations in the pressure.
307 gross 2100 """
308 gross 1414 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
309     self.domain=domain
310     self.vol=util.integrate(1.,Function(self.domain))
311     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
312     self.__pde_u.setSymmetryOn()
313 gross 2100 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
314     # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
315 caltinay 2169
316 gross 1414 self.__pde_prec=LinearPDE(domain)
317     self.__pde_prec.setReducedOrderOn()
318 gross 2156 # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
319 gross 1414 self.__pde_prec.setSymmetryOn()
320    
321     self.__pde_proj=LinearPDE(domain)
322     self.__pde_proj.setReducedOrderOn()
323     self.__pde_proj.setSymmetryOn()
324     self.__pde_proj.setValue(D=1.)
325    
326 gross 2100 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
327 caltinay 2169 """
328     Assigns values to the model parameters.
329 gross 2100
330 caltinay 2169 @param f: external force
331     @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
332     @param fixed_u_mask: mask of locations with fixed velocity
333     @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}
334     or similar
335     @param eta: viscosity
336     @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar
337     @param surface_stress: normal surface stress
338     @type surface_stress: L{Vector} object on L{FunctionSpace},
339     L{FunctionOnBoundary} or similar
340     @param stress: initial stress
341     @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or
342     similar
343     @note: All values need to be set.
344     """
345     self.eta=eta
346     A =self.__pde_u.createCoefficient("A")
347     self.__pde_u.setValue(A=Data())
348     for i in range(self.domain.getDim()):
349     for j in range(self.domain.getDim()):
350     A[i,j,j,i] += 1.
351     A[i,j,i,j] += 1.
352     self.__pde_prec.setValue(D=1/self.eta)
353     self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
354     self.__stress=stress
355 gross 1414
356 gross 2100 def B(self,v):
357 caltinay 2169 """
358     Returns M{div(v)}.
359     @return: M{div(v)}
360     @rtype: equal to the type of p
361 gross 1414
362 caltinay 2169 @note: Boundary conditions on p should be zero!
363     """
364     if self.show_details: print "apply divergence:"
365     self.__pde_proj.setValue(Y=-util.div(v))
366     self.__pde_proj.setTolerance(self.getSubProblemTolerance())
367     return self.__pde_proj.getSolution(verbose=self.show_details)
368 gross 2100
369     def inner_pBv(self,p,Bv):
370     """
371 caltinay 2169 Returns inner product of element p and Bv (overwrite).
372    
373 gross 2100 @type p: equal to the type of p
374     @type Bv: equal to the type of result of operator B
375 caltinay 2169 @return: inner product of p and Bv
376 gross 2100 @rtype: equal to the type of p
377     """
378     s0=util.interpolate(p,Function(self.domain))
379     s1=util.interpolate(Bv,Function(self.domain))
380 gross 1414 return util.integrate(s0*s1)
381    
382 gross 2100 def inner_p(self,p0,p1):
383     """
384 caltinay 2169 Returns inner product of element p0 and p1 (overwrite).
385    
386 gross 2100 @type p0: equal to the type of p
387     @type p1: equal to the type of p
388 caltinay 2169 @return: inner product of p0 and p1
389 gross 2100 @rtype: equal to the type of p
390     """
391     s0=util.interpolate(p0/self.eta,Function(self.domain))
392     s1=util.interpolate(p1/self.eta,Function(self.domain))
393     return util.integrate(s0*s1)
394 artak 1550
395 gross 2100 def inner_v(self,v0,v1):
396     """
397 caltinay 2169 Returns inner product of two elements v0 and v1 (overwrite).
398    
399 gross 2100 @type v0: equal to the type of v
400     @type v1: equal to the type of v
401 caltinay 2169 @return: inner product of v0 and v1
402 gross 2100 @rtype: equal to the type of v
403     """
404 caltinay 2169 gv0=util.grad(v0)
405     gv1=util.grad(v1)
406 gross 2100 return util.integrate(util.inner(gv0,gv1))
407    
408 gross 1414 def solve_A(self,u,p):
409     """
410 caltinay 2169 Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).
411 gross 1414 """
412 gross 2100 if self.show_details: print "solve for velocity:"
413 gross 1414 self.__pde_u.setTolerance(self.getSubProblemTolerance())
414 gross 2100 if self.__stress.isEmpty():
415     self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
416     else:
417     self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
418     out=self.__pde_u.getSolution(verbose=self.show_details)
419 caltinay 2169 return out
420 gross 1414
421     def solve_prec(self,p):
422 caltinay 2169 """
423     Applies the preconditioner.
424     """
425 gross 2100 if self.show_details: print "apply preconditioner:"
426 gross 1414 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
427     self.__pde_prec.setValue(Y=p)
428     q=self.__pde_prec.getSolution(verbose=self.show_details)
429 artak 1554 return q
430 caltinay 2169

  ViewVC Help
Powered by ViewVC 1.1.26