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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2156 - (hide annotations)
Mon Dec 15 05:09:02 2008 UTC (10 years, 8 months ago) by gross
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 17615 byte(s)
some modifications to the iterative solver to make them easier to use. 
There are also improved versions of the Darcy flux solver and the incompressible solver.


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     solves the problem
42 gross 1659
43 gross 2100 M{u_i+k_{ij}*p_{,j} = g_i}
44     M{u_{i,i} = f}
45 gross 1659
46 gross 2100 where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47 gross 1659
48 gross 2100 @note: The problem is solved in a least squares formulation.
49     """
50 gross 1659
51 gross 2100 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 gross 1659
66 gross 2100 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 gross 1659
70 gross 2100 @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 gross 1659
83 gross 2100 @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 gross 1659
104 gross 2100 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 gross 1659
107 gross 2100 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 gross 1659
121    
122 gross 2156 def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
123 gross 2100 """
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 gross 2123 self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
141 gross 2156 return self.__pde_v.getSolution(verbose=show_details)
142 gross 1659
143 gross 2100 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 gross 1659
150 gross 2100 @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 gross 1659
169 gross 2100 M{(I+D^*D)u+Qp=D^*f+g}
170     M{Q^*u+Q^*Qp=Q^*g}
171 gross 1659
172 gross 2100 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 gross 1659
175 gross 2100 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     u2=u0*self.__pde_v.getCoefficient("q")
189 gross 2156 #
190     # first the reference velocity is calculated from
191     #
192     # (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
193     #
194     self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
195     u_ref=self.__pde_v.getSolution(verbose=show_details)
196     if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
197     self.__pde_v.setValue(r=Data())
198     #
199     # and then we calculate a reference pressure
200     #
201     # Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
202     #
203     self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
204     p_ref=self.__pde_p.getSolution(verbose=self.show_details)
205     if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
206 gross 2100 self.__pde_p.setValue(r=Data())
207 gross 2156 #
208     # (I+D^*D)du + Qdp = - Qp_ref u=du+u_ref
209     # Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0 p=dp+pref
210     #
211     # du= -(I+D^*D)^(-1} Q(p_ref+dp) u = u_ref+du
212     #
213     # => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
214     # or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
215     #
216     # r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
217     # with du=-(I+D^*D)^(-1} Q(p_ref+dp)
218     #
219     # we use the (du,Qdp) to represent the resudual
220     # Q^*Q is a preconditioner
221     #
222     # <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
223     #
224     Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
225     norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
226     ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
227 gross 2100 if not ATOL>0:
228 gross 2156 raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, norm_rhs, atol)
229     if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
230     #
231     # caclulate the initial residual
232     #
233     self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
234     du=self.__pde_v.getSolution(verbose=show_details)
235     r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
236     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)
237     util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
238     return u_ref+r[1],dp
239 gross 2100
240     def __Aprod_PCG(self,p):
241     if self.show_details: print "DarcyFlux: Applying operator"
242     Qp=util.tensor_mult(self.__permeability,util.grad(p))
243     self.__pde_v.setValue(Y=Qp,X=Data())
244     w=self.__pde_v.getSolution(verbose=self.show_details)
245 gross 2156 return ArithmeticTuple(-Qp,w)
246 gross 2100
247     def __inner_PCG(self,p,r):
248     a=util.tensor_mult(self.__permeability,util.grad(p))
249 gross 2156 out=-util.integrate(util.inner(a,r[0]+r[1]))
250     return out
251 gross 2100
252     def __Msolve_PCG(self,r):
253     if self.show_details: print "DarcyFlux: Applying preconditioner"
254 gross 2156 self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))
255 gross 2100 return self.__pde_p.getSolution(verbose=self.show_details)
256    
257 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
258     """
259     solves
260    
261 gross 2100 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
262 gross 1414 u_{i,i}=0
263    
264     u=0 where fixed_u_mask>0
265 gross 2100 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
266 gross 1414
267 gross 2100 if surface_stress is not given 0 is assumed.
268 gross 1414
269     typical usage:
270    
271     sp=StokesProblemCartesian(domain)
272     sp.setTolerance()
273     sp.initialize(...)
274     v,p=sp.solve(v0,p0)
275     """
276     def __init__(self,domain,**kwargs):
277 gross 2100 """
278     initialize the Stokes Problem
279    
280     @param domain: domain of the problem. The approximation order needs to be two.
281     @type domain: L{Domain}
282     @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
283     """
284 gross 1414 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
285     self.domain=domain
286     self.vol=util.integrate(1.,Function(self.domain))
287     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
288     self.__pde_u.setSymmetryOn()
289 gross 2100 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
290     # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
291 gross 1414
292     self.__pde_prec=LinearPDE(domain)
293     self.__pde_prec.setReducedOrderOn()
294 gross 2156 # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
295 gross 1414 self.__pde_prec.setSymmetryOn()
296    
297     self.__pde_proj=LinearPDE(domain)
298     self.__pde_proj.setReducedOrderOn()
299     self.__pde_proj.setSymmetryOn()
300     self.__pde_proj.setValue(D=1.)
301    
302 gross 2100 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
303     """
304     assigns values to the model parameters
305    
306     @param f: external force
307     @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
308     @param fixed_u_mask: mask of locations with fixed velocity.
309     @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
310     @param eta: viscosity
311     @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
312     @param surface_stress: normal surface stress
313     @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
314     @param stress: initial stress
315     @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
316     @note: All values needs to be set.
317    
318     """
319 gross 1414 self.eta=eta
320 gross 1841 A =self.__pde_u.createCoefficient("A")
321 gross 1414 self.__pde_u.setValue(A=Data())
322     for i in range(self.domain.getDim()):
323     for j in range(self.domain.getDim()):
324     A[i,j,j,i] += 1.
325     A[i,j,i,j] += 1.
326 artak 1554 self.__pde_prec.setValue(D=1/self.eta)
327 gross 1414 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
328 gross 2100 self.__stress=stress
329 gross 1414
330 gross 2100 def B(self,v):
331     """
332     returns div(v)
333     @rtype: equal to the type of p
334 gross 1414
335 gross 2100 @note: boundary conditions on p should be zero!
336     """
337     if self.show_details: print "apply divergence:"
338     self.__pde_proj.setValue(Y=-util.div(v))
339     self.__pde_proj.setTolerance(self.getSubProblemTolerance())
340     return self.__pde_proj.getSolution(verbose=self.show_details)
341    
342     def inner_pBv(self,p,Bv):
343     """
344     returns inner product of element p and Bv (overwrite)
345    
346     @type p: equal to the type of p
347     @type Bv: equal to the type of result of operator B
348     @rtype: C{float}
349    
350     @rtype: equal to the type of p
351     """
352     s0=util.interpolate(p,Function(self.domain))
353     s1=util.interpolate(Bv,Function(self.domain))
354 gross 1414 return util.integrate(s0*s1)
355    
356 gross 2100 def inner_p(self,p0,p1):
357     """
358     returns inner product of element p0 and p1 (overwrite)
359    
360     @type p0: equal to the type of p
361     @type p1: equal to the type of p
362     @rtype: C{float}
363 artak 1550
364 gross 2100 @rtype: equal to the type of p
365     """
366     s0=util.interpolate(p0/self.eta,Function(self.domain))
367     s1=util.interpolate(p1/self.eta,Function(self.domain))
368     return util.integrate(s0*s1)
369 artak 1550
370 gross 2100 def inner_v(self,v0,v1):
371     """
372     returns inner product of two element v0 and v1 (overwrite)
373    
374     @type v0: equal to the type of v
375     @type v1: equal to the type of v
376     @rtype: C{float}
377 gross 1414
378 gross 2100 @rtype: equal to the type of v
379     """
380     gv0=util.grad(v0)
381     gv1=util.grad(v1)
382     return util.integrate(util.inner(gv0,gv1))
383    
384 gross 1414 def solve_A(self,u,p):
385     """
386     solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
387     """
388 gross 2100 if self.show_details: print "solve for velocity:"
389 gross 1414 self.__pde_u.setTolerance(self.getSubProblemTolerance())
390 gross 2100 if self.__stress.isEmpty():
391     self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
392     else:
393     self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
394     out=self.__pde_u.getSolution(verbose=self.show_details)
395     return out
396 gross 1414
397     def solve_prec(self,p):
398 gross 2100 if self.show_details: print "apply preconditioner:"
399 gross 1414 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
400     self.__pde_prec.setValue(Y=p)
401     q=self.__pde_prec.getSolution(verbose=self.show_details)
402 artak 1554 return q

  ViewVC Help
Powered by ViewVC 1.1.26