# Diff of /trunk/escript/py_src/flows.py

revision 1465 by artak, Wed Apr 2 03:28:25 2008 UTC revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC
# Line 1  Line 1
1  # \$Id:\$  ########################################################
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  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
#
#######################################################
10  #  #
11    ########################################################
12
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
19    __url__="http://www.uq.edu.au/esscc/escript-finley"
20
21  """  """
22  Some models for flow  Some models for flow
# Line 24  Some models for flow Line 30  Some models for flow
30  """  """
31
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
http://www.access.edu.au
__url__="http://www.iservo.edu.au/esys"
__version__="\$Revision:\$"
__date__="\$Date:\$"
33
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
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=boundary_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")
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, verbose=self.verbose, initial_guess=False)
205             return u2+r[1],p2+dp
206
207        def __Aprod_PCG(self,p):
208              if self.show_details: print "DarcyFlux: Applying operator"
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):
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
224        """        """
225        solves        solves
226
227            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
228                  u_{i,i}=0                  u_{i,i}=0
229
231            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            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 give 0 is assumed.        if surface_stress is not given 0 is assumed.
234
235        typical usage:        typical usage:
236
# Line 58  class StokesProblemCartesian(Homogeneous Line 240  class StokesProblemCartesian(Homogeneous
240              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
241        """        """
242        def __init__(self,domain,**kwargs):        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             """
251           self.domain=domain           self.domain=domain
252           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
253           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
254           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
255           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
256             # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
257
258           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
259           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
260             self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
261           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
262
263           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
# Line 74  class StokesProblemCartesian(Homogeneous Line 265  class StokesProblemCartesian(Homogeneous
265           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
266           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
267
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
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          self.eta=eta
286          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
287      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
288          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
289          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
290              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
291              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
292      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
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 B(self,arg):        def inner_pBv(self,p,Bv):
309           d=util.div(arg)           """
310           self.__pde_proj.setValue(Y=d)           returns inner product of element p and Bv  (overwrite)
311           self.__pde_proj.setTolerance(self.getSubProblemTolerance())
312           return self.__pde_proj.getSolution(verbose=self.show_details)           @type p: equal to the type of p
313             @type Bv: equal to the type of result of operator B
314        def inner(self,p0,p1):           @rtype: C{float}
315           s0=util.interpolate(p0,Function(self.domain))
316           s1=util.interpolate(p1,Function(self.domain))           @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)           return util.integrate(s0*s1)
321
322        def getStress(self,u):        def inner_p(self,p0,p1):
324           return 2.*self.eta*util.symmetric(mg)           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             """
348             return util.integrate(util.inner(gv0,gv1))
349
350        def solve_A(self,u,p):        def solve_A(self,u,p):
351           """           """
353           """           """
354             if self.show_details: print "solve for velocity:"
355           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
356           self.__pde_u.setValue(X=-self.getStress(u)+p*util.kronecker(self.domain))           if self.__stress.isEmpty():
358             else:
360             out=self.__pde_u.getSolution(verbose=self.show_details)
361             return  out
362
363        def solve_prec(self,p):        def solve_prec(self,p):
364             if self.show_details: print "apply preconditioner:"
365           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
366           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
367           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
368           return q           return q
def stoppingcriterium(self,Bv,v,p):
n_r=util.sqrt(self.inner(Bv,Bv))
n_v=util.Lsup(v)
if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)
self.iter+=1
if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():
if self.verbose: print "PCG terminated after %s steps."%self.iter
return True
else:
return False
def stoppingcriterium_GMRES(self,rho,r):
if self.verbose: print "GMRES step %s: L2(rho) = %s, L2(b)*TOL=%s"%(self.iter,rho,r*self.getTolerance())
self.iter+=1
if rho <= r*self.getTolerance():
if self.verbose: print "GMRES terminated after %s steps."%self.iter
return True
else:
return False

Legend:
 Removed from v.1465 changed lines Added in v.2100