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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2549 by jfenwick, Mon Jul 20 06:43:47 2009 UTC revision 3452 by caltinay, Tue Jan 25 01:53:57 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 21  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40    print dir(escript)
41    
42  class DarcyFlow(object):  class DarcyFlow(object):
43       """
44       solves the problem
45      
46       *u_i+k_{ij}*p_{,j} = g_i*
47       *u_{i,i} = f*
48      
49       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
50      
51       :note: The problem is solved in a least squares formulation.
52       """
53      
54       def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
55          """
56          initializes the Darcy flux problem
57          :param domain: domain of the problem
58          :type domain: `Domain`
59          :param useReduced: uses reduced oreder on flux and pressure
60          :type useReduced: ``bool``
61          :param adaptSubTolerance: switches on automatic subtolerance selection
62          :type adaptSubTolerance: ``bool``
63          :param solveForFlux: if True the solver solves for the flux (do not use!)
64          :type solveForFlux: ``bool``  
65          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
66          :type useVPIteration: ``bool``    
67          """
68          self.domain=domain
69          self.useVPIteration=useVPIteration
70          self.useReduced=useReduced
71          self.weighting_scale=weighting_scale
72          if self.useVPIteration:
73             self.solveForFlux=solveForFlux
74             self.__adaptSubTolerance=adaptSubTolerance
75             self.verbose=False
76            
77             self.__pde_k=LinearPDESystem(domain)
78             self.__pde_k.setSymmetryOn()
79             if self.useReduced: self.__pde_k.setReducedOrderOn()
80      
81             self.__pde_p=LinearSinglePDE(domain)
82             self.__pde_p.setSymmetryOn()
83             if self.useReduced: self.__pde_p.setReducedOrderOn()
84             self.setTolerance()
85             self.setAbsoluteTolerance()
86          else:
87             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
88             self.__pde_k.setSymmetryOn()
89             if self.useReduced: self.__pde_k.setReducedOrderOn()
90             C=self.__pde_k.createCoefficient("C")
91             B=self.__pde_k.createCoefficient("B")
92             for i in range(self.domain.getDim()):
93                C[i,self.domain.getDim(),i]=1
94                B[self.domain.getDim(),i,i]=1
95             self.__pde_k.setValue(C=C, B=B)
96          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
97          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
98          
99       def getSolverOptionsFlux(self):
100          """
101          Returns the solver options used to solve the flux problems
102          
103          *K^{-1} u=F*
104          
105          :return: `SolverOptions`
106          """
107          return self.__pde_k.getSolverOptions()
108          
109       def setSolverOptionsFlux(self, options=None):
110          """
111          Sets the solver options used to solve the flux problems
112          
113          *K^{-1}u=F*
114          
115          If ``options`` is not present, the options are reset to default
116          
117          :param options: `SolverOptions`
118          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
119          """
120          return self.__pde_v.setSolverOptions(options)
121        
122       def getSolverOptionsPressure(self):
123          """
124          Returns the solver options used to solve the pressure problems
125          
126          *(Q^* K Q)p=-Q^*G*
127          
128          :return: `SolverOptions`
129          """
130          return self.__pde_p.getSolverOptions()
131          
132       def setSolverOptionsPressure(self, options=None):
133          """
134          Sets the solver options used to solve the pressure problems
135          
136          *(Q^* K Q)p=-Q^*G*
137          
138          If ``options`` is not present, the options are reset to default
139          
140          :param options: `SolverOptions`
141          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
142          """
143          return self.__pde_p.setSolverOptions(options)
144    
145    
146       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
147          """
148          assigns values to model parameters
149    
150          :param f: volumetic sources/sinks
151          :type f: scalar value on the domain (e.g. `escript.Data`)
152          :param g: flux sources/sinks
153          :type g: vector values on the domain (e.g. `escript.Data`)
154          :param location_of_fixed_pressure: mask for locations where pressure is fixed
155          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
156          :param location_of_fixed_flux:  mask for locations where flux is fixed.
157          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
158          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
159          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
160    
161          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
162          :note: at any point on the boundary of the domain the pressure
163                 (``location_of_fixed_pressure`` >0) or the normal component of the
164                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
165                 is along the *x_i* axis.
166    
167          """
168          if self.useVPIteration:
169             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
170             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)
171          else:
172             q=self.__pde_k.getCoefficient("q")
173             if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
174             if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
175             if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
176             self.__pde_k.setValue(q=q)
177                
178          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
179          if permeability!=None:
180         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
181             V=util.vol(self.domain)
182         if perm.getRank()==0:
183            perm_inv=(1./perm)
184                if self.useVPIteration:
185                  self.scale=1.
186                else:
187                  self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
188    
189            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
190            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
191         elif perm.getRank()==2:
192            perm_inv=util.inverse(perm)
193                if self.useVPIteration:
194                  self.scale=1.
195                else:
196                  self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
197              perm_inv*=(1./self.scale)
198              perm=perm*self.scale
199         else:
200            raise ValueError,"illegal rank of permeability."
201    
202         self.__permeability=perm
203         self.__permeability_inv=perm_inv
204    
205         self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
206             if self.useVPIteration:
207            if  self.solveForFlux:
208               self.__pde_k.setValue(D=self.__permeability_inv)
209            else:
210               self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
211            self.__pde_p.setValue(A=self.__permeability)
212             else:
213                D=self.__pde_k.createCoefficient("D")
214                A=self.__pde_k.createCoefficient("A")
215                D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
216                for i in range(self.domain.getDim()):
217                   for j in range(self.domain.getDim()):
218                     A[i,i,j,j]=self.__l2
219                A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
220                self.__pde_k.setValue(A=A, D=D)
221          if g !=None:
222         g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
223         if g.isEmpty():
224              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
225         else:
226            if not g.getShape()==(self.domain.getDim(),):
227                  raise ValueError,"illegal shape of g"
228            self.__g=g
229          elif permeability!=None:
230                 X
231          if f !=None:
232         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
233         if f.isEmpty():
234              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
235         else:
236             if f.getRank()>0: raise ValueError,"illegal rank of f."
237             self.__f=f
238    
239          
240       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
241          """
242          solves the problem.
243          
244          The iteration is terminated if the residual norm is less then self.getTolerance().
245    
246          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
247          :type u0: vector value on the domain (e.g. `escript.Data`).
248          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
249          :type p0: scalar value on the domain (e.g. `escript.Data`).
250          :param verbose: if set some information on iteration progress are printed
251          :type verbose: ``bool``
252          :return: flux and pressure
253          :rtype: ``tuple`` of `escript.Data`.
254    
255          :note: The problem is solved as a least squares form
256                 *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
257                 *G^*u+*G^* K Gp=G^*g*
258                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
259          """
260          self.verbose=verbose
261          if self.useVPIteration:
262              return self.__solveVP(u0,p0,max_iter,max_num_corrections)
263          else:
264              X=self.__pde_k.createCoefficient("X")
265              Y=self.__pde_k.createCoefficient("Y")
266              Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
267              rtmp=self.__f * self.__l2 * self.scale
268              for i in range(self.domain.getDim()): X[i,i]=rtmp
269              X[self.domain.getDim(),:]=self.__g*self.scale
270              r=self.__pde_k.createCoefficient("r")
271              r[:self.domain.getDim()]=u0*self.scale
272              r[self.domain.getDim()]=p0
273              self.__pde_k.setValue(X=X, Y=Y, r=r)
274              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
275              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
276              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
277              U=self.__pde_k.getSolution()
278              # self.__pde_k.getOperator().saveMM("k.mm")
279              u=U[:self.domain.getDim()]*(1./self.scale)
280              p=U[self.domain.getDim()]
281              if self.verbose:
282            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
283            def_p=self.__g-(u+KGp)
284            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
285            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
286            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
287              return u,p
288    
289       def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
290          """
291          solves the problem.
292          
293          The iteration is terminated if the residual norm is less than self.getTolerance().
294    
295          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
296          :type u0: vector value on the domain (e.g. `escript.Data`).
297          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
298          :type p0: scalar value on the domain (e.g. `escript.Data`).
299          :return: flux and pressure
300          :rtype: ``tuple`` of `escript.Data`.
301    
302          :note: The problem is solved as a least squares form
303                 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
304                 *G^*u+*G^* K Gp=G^*g*
305                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
306          """
307          rtol=self.getTolerance()
308          atol=self.getAbsoluteTolerance()
309          self.setSubProblemTolerance()
310          num_corrections=0
311          converged=False
312          norm_r=None
313          
314          # Eliminate the hydrostatic pressure:
315          if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
316          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))        
317          p0=self.__pde_p.getSolution()
318          g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
319          norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5    
320    
321          p=p0*0
322          if self.solveForFlux:
323         v=u0.copy()
324          else:
325         v=self.__getFlux(p, u0, f=self.__f, g=g2)
326    
327          while not converged and norm_g2 > 0:
328         Gp=util.grad(p)
329         KGp=util.tensor_mult(self.__permeability,Gp)
330         if self.verbose:
331            def_p=g2-(v+KGp)
332            def_v=self.__f-util.div(v)
333            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
334            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
335            print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
336            print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
337            print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
338         ATOL=atol+rtol*norm_g2
339         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
340         if norm_r == None or norm_r>ATOL:
341            if num_corrections>max_num_corrections:
342               raise ValueError,"maximum number of correction steps reached."
343          
344            if self.solveForFlux:
345               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
346               v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
347                   self.__Aprod_v,
348                   v,
349                   self.__Msolve_PCG_v,
350                   self.__inner_PCG_v,
351                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
352               p=r[2]
353            else:
354               # initial residual is r=G^*(g2-KGp - v)
355               p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
356                     self.__Aprod_p,
357                     p,
358                     self.__Msolve_PCG_p,
359                     self.__inner_PCG_p,
360                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
361               v=r[1]
362            if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
363            num_corrections+=1
364         else:
365            if self.verbose: print "DarcyFlux: stopping criterium reached."
366            converged=True
367          return v,p+p0
368    
369       def __applWeight(self, v, f=None):
370          # solves L p = f-Dv with p = 0
371          if self.verbose: print "DarcyFlux: Applying weighting operator"
372          if f == None:
373         return -util.div(v)*self.__l2
374          else:
375         return (f-util.div(v))*self.__l2
376       def __getPressure(self, v, p0, g=None):
377          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
378          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
379          if g == None:
380         self.__pde_p.setValue(X=-v, r=p0)
381          else:
382         self.__pde_p.setValue(X=g-v, r=p0)      
383          p=self.__pde_p.getSolution()
384          return p
385    
386       def __Aprod_v(self,dv):
387          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
388          dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
389          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
390          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
391          return ArithmeticTuple(a,b,-dp)
392    
393       def __Msolve_PCG_v(self,r):
394          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
395          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
396          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
397          return self.__pde_k.getSolution()
398    
399       def __inner_PCG_v(self,v,r):
400          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
401          
402       def __Aprod_p(self,dp):
403          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
404          Gdp=util.grad(dp)
405          self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
406          du=self.__pde_k.getSolution()
407          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
408    
409       def __getFlux(self,p, v0, f=None, g=None):
410          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
411          if f!=None:
412         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
413          self.__pde_k.setValue(r=v0)
414          g2=util.tensor_mult(self.__permeability_inv,g)
415          if p == None:
416         self.__pde_k.setValue(Y=g2)
417          else:
418         self.__pde_k.setValue(Y=g2-util.grad(p))
419          return self.__pde_k.getSolution()  
420          
421          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
422       def __Msolve_PCG_p(self,r):
423          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
424          self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
425          return self.__pde_p.getSolution()
426            
427       def __inner_PCG_p(self,p,r):
428           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
429    
430       def __L2(self,v):
431          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
432    
433       def __L2_r(self,v):
434          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
435    
436       def setTolerance(self,rtol=1e-4):
437          """
438          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
439    
440          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
441          
442          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
443          
444          :param rtol: relative tolerance for the pressure
445          :type rtol: non-negative ``float``
446          """
447          if rtol<0:
448         raise ValueError,"Relative tolerance needs to be non-negative."
449          self.__rtol=rtol
450       def getTolerance(self):
451          """
452          returns the relative tolerance
453          :return: current relative tolerance
454          :rtype: ``float``
455          """
456          return self.__rtol
457    
458       def setAbsoluteTolerance(self,atol=0.):
459          """
460          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
461          
462          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
463    
464    
465          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
466    
467          :param atol: absolute tolerance for the pressure
468          :type atol: non-negative ``float``
469          """
470          if atol<0:
471         raise ValueError,"Absolute tolerance needs to be non-negative."
472          self.__atol=atol
473       def getAbsoluteTolerance(self):
474          """
475          returns the absolute tolerance
476          :return: current absolute tolerance
477          :rtype: ``float``
478          """
479          return self.__atol
480       def getSubProblemTolerance(self):
481          """
482          Returns a suitable subtolerance
483          :type: ``float``
484          """
485          return max(util.EPSILON**(0.5),self.getTolerance()**2)
486    
487       def setSubProblemTolerance(self):
488          """
489          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
490          """
491          if self.__adaptSubTolerance:
492         sub_tol=self.getSubProblemTolerance()
493         self.getSolverOptionsFlux().setTolerance(sub_tol)
494         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
495         self.getSolverOptionsPressure().setTolerance(sub_tol)
496         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
497         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
498    
499    
500    class DarcyFlowOld(object):
501      """      """
502      solves the problem      solves the problem
503    
504      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
505      M{u_{i,i} = f}      *u_{i,i} = f*
506    
507      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
508    
509      @note: The problem is solved in a least squares formulation.      :note: The problem is solved in a least squares formulation.
510      """      """
511    
512      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
513          """          """
514          initializes the Darcy flux problem          initializes the Darcy flux problem
515          @param domain: domain of the problem          :param domain: domain of the problem
516          @type domain: L{Domain}          :type domain: `Domain`
517      @param useReduced: uses reduced oreder on flux and pressure      :param useReduced: uses reduced oreder on flux and pressure
518      @type useReduced: C{bool}      :type useReduced: ``bool``
519      @param adaptSubTolerance: switches on automatic subtolerance selection      :param adaptSubTolerance: switches on automatic subtolerance selection
520      @type adaptSubTolerance: C{bool}          :type adaptSubTolerance: ``bool``  
521          """          """
522          self.domain=domain          self.domain=domain
523          if weight == None:          if weight == None:
524             s=self.domain.getSize()             s=self.domain.getSize()
525             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
526             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l2=(3.*util.longestEdge(self.domain))**2
527             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
528          else:          else:
529             self.__l=weight             self.__l2=weight
530          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
531          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
532          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
533          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
534          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
535          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
536          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
537          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
538          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
539          self.setTolerance()          self.setTolerance()
540          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
541      self.__adaptSubTolerance=adaptSubTolerance      self.__adaptSubTolerance=adaptSubTolerance
# Line 83  class DarcyFlow(object): Line 544  class DarcyFlow(object):
544      """      """
545      Returns the solver options used to solve the flux problems      Returns the solver options used to solve the flux problems
546            
547      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
548            
549      @return: L{SolverOptions}      :return: `SolverOptions`
550      """      """
551      return self.__pde_v.getSolverOptions()      return self.__pde_v.getSolverOptions()
552      def setSolverOptionsFlux(self, options=None):      def setSolverOptionsFlux(self, options=None):
553      """      """
554      Sets the solver options used to solve the flux problems      Sets the solver options used to solve the flux problems
555            
556      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
557            
558      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
559      @param options: L{SolverOptions}      :param options: `SolverOptions`
560      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
561      """      """
562      return self.__pde_v.setSolverOptions(options)      return self.__pde_v.setSolverOptions(options)
563      def getSolverOptionsPressure(self):      def getSolverOptionsPressure(self):
564      """      """
565      Returns the solver options used to solve the pressure problems      Returns the solver options used to solve the pressure problems
566            
567      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
568            
569      @return: L{SolverOptions}      :return: `SolverOptions`
570      """      """
571      return self.__pde_p.getSolverOptions()      return self.__pde_p.getSolverOptions()
572      def setSolverOptionsPressure(self, options=None):      def setSolverOptionsPressure(self, options=None):
573      """      """
574      Sets the solver options used to solve the pressure problems      Sets the solver options used to solve the pressure problems
575            
576      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
577            
578      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
579      @param options: L{SolverOptions}      :param options: `SolverOptions`
580      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
581      """      """
582      return self.__pde_p.setSolverOptions(options)      return self.__pde_p.setSolverOptions(options)
583    
# Line 124  class DarcyFlow(object): Line 585  class DarcyFlow(object):
585          """          """
586          assigns values to model parameters          assigns values to model parameters
587    
588          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
589          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `escript.Data`)
590          @param g: flux sources/sinks          :param g: flux sources/sinks
591          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `escript.Data`)
592          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
593          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
594          @param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
595          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
596          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
597                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
598                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
599          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
600    
601          @note: the values of parameters which are not set by calling C{setValue} are not altered.          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
602          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
603                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
604                 is along the M{x_i} axis.                 is along the *x_i* axis.
605          """          """
606          if f !=None:          if f !=None:
607             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
608             if f.isEmpty():             if f.isEmpty():
609                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
610             else:             else:
611                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
612             self.__f=f             self.__f=f
613          if g !=None:          if g !=None:
614             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
615             if g.isEmpty():             if g.isEmpty():
616               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
617             else:             else:
618               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
619                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 177  class DarcyFlow(object): Line 638  class DarcyFlow(object):
638    
639      def setTolerance(self,rtol=1e-4):      def setTolerance(self,rtol=1e-4):
640          """          """
641          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
642    
643          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
644    
645          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}}.          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
646    
647          @param rtol: relative tolerance for the pressure          :param rtol: relative tolerance for the pressure
648          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
649          """          """
650          if rtol<0:          if rtol<0:
651              raise ValueError,"Relative tolerance needs to be non-negative."              raise ValueError,"Relative tolerance needs to be non-negative."
# Line 193  class DarcyFlow(object): Line 654  class DarcyFlow(object):
654          """          """
655          returns the relative tolerance          returns the relative tolerance
656    
657          @return: current relative tolerance          :return: current relative tolerance
658          @rtype: C{float}          :rtype: ``float``
659          """          """
660          return self.__rtol          return self.__rtol
661    
662      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
663          """          """
664          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
665    
666          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
667    
668          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}}.          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
669    
670          @param atol: absolute tolerance for the pressure          :param atol: absolute tolerance for the pressure
671          @type atol: non-negative C{float}          :type atol: non-negative ``float``
672          """          """
673          if atol<0:          if atol<0:
674              raise ValueError,"Absolute tolerance needs to be non-negative."              raise ValueError,"Absolute tolerance needs to be non-negative."
# Line 216  class DarcyFlow(object): Line 677  class DarcyFlow(object):
677         """         """
678         returns the absolute tolerance         returns the absolute tolerance
679                
680         @return: current absolute tolerance         :return: current absolute tolerance
681         @rtype: C{float}         :rtype: ``float``
682         """         """
683         return self.__atol         return self.__atol
684      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
685      """      """
686      Returns a suitable subtolerance      Returns a suitable subtolerance
687      @type: C{float}      @type: ``float``
688      """      """
689      return max(util.EPSILON**(0.75),self.getTolerance()**2)      return max(util.EPSILON**(0.75),self.getTolerance()**2)
690      def setSubProblemTolerance(self):      def setSubProblemTolerance(self):
# Line 244  class DarcyFlow(object): Line 705  class DarcyFlow(object):
705    
706           The iteration is terminated if the residual norm is less then self.getTolerance().           The iteration is terminated if the residual norm is less then self.getTolerance().
707    
708           @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.           :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
709           @type u0: vector value on the domain (e.g. L{Data}).           :type u0: vector value on the domain (e.g. `escript.Data`).
710           @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.           :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
711           @type p0: scalar value on the domain (e.g. L{Data}).           :type p0: scalar value on the domain (e.g. `escript.Data`).
712           @param verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
713           @type verbose: C{bool}           :type verbose: ``bool``
714           @return: flux and pressure           :return: flux and pressure
715           @rtype: C{tuple} of L{Data}.           :rtype: ``tuple`` of `escript.Data`.
716    
717           @note: The problem is solved as a least squares form           :note: The problem is solved as a least squares form
718    
719           M{(I+D^*D)u+Qp=D^*f+g}           *(I+D^*D)u+Qp=D^*f+g*
720           M{Q^*u+Q^*Qp=Q^*g}           *Q^*u+Q^*Qp=Q^*g*
721    
722           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
723           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
724    
725           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
726    
727           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
728    
729           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           *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
730    
731           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
732           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
733           """           """
734           self.verbose=verbose           self.verbose=verbose
735           rtol=self.getTolerance()           rtol=self.getTolerance()
736           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
737       self.setSubProblemTolerance()       self.setSubProblemTolerance()
       
738           num_corrections=0           num_corrections=0
739           converged=False           converged=False
740           p=p0           p=p0
# Line 297  class DarcyFlow(object): Line 757  class DarcyFlow(object):
757                 ATOL=(atol+rtol*fac)                 ATOL=(atol+rtol*fac)
758                 if self.verbose:                 if self.verbose:
759                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
760                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
761                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
762                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
763                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
764                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
765                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
766                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
767                     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)                     p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.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)
768                     num_corrections+=1                     num_corrections+=1
769                 else:                 else:
770                     converged=True                     converged=True
771           return v,p           return v,p
772      def __L2(self,v):      def __L2(self,v):
773           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
774    
775      def __Q(self,p):      def __Q(self,p):
776            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
# Line 318  class DarcyFlow(object): Line 778  class DarcyFlow(object):
778      def __Aprod(self,dp):      def __Aprod(self,dp):
779            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
780            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
781            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
782            du=self.__pde_v.getSolution()            du=self.__pde_v.getSolution()
           # self.__pde_v.getOperator().saveMM("proj.mm")  
783            return Qdp+du            return Qdp+du
784      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
785           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 330  class DarcyFlow(object): Line 789  class DarcyFlow(object):
789    
790      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
791        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
792            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
           # self.__pde_p.getOperator().saveMM("prec.mm")  
793            return self.__pde_p.getSolution()            return self.__pde_p.getSolution()
794    
795      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
796          """          """
797          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
798          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
799          Note that C{g} and C{f} are used, see L{setValue}.          Note that ``g`` and ``f`` are used, see `setValue`.
800    
801          @param p: pressure.          :param p: pressure.
802          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `escript.Data`).
803          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
804          @type fixed_flux: vector values on the domain (e.g. L{Data}).          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
805          @param tol: relative tolerance to be used.          :return: flux
806          @type tol: positive C{float}.          :rtype: `escript.Data`
807          @return: flux          :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
808          @rtype: L{Data}                 for the permeability *k_{ij}*
         @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}}  
                for the permeability M{k_{ij}}  
809          """          """
810      self.setSubProblemTolerance()      self.setSubProblemTolerance()
811          g=self.__g          g=self.__g
812          f=self.__f          f=self.__f
813          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
814          if p == None:          if p == None:
815             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
816          else:          else:
# Line 380  class StokesProblemCartesian(Homogeneous Line 836  class StokesProblemCartesian(Homogeneous
836              sp.initialize(...)              sp.initialize(...)
837              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
838       """       """
839       def __init__(self,domain,adaptSubTolerance=True, **kwargs):       def __init__(self,domain,**kwargs):
840           """           """
841           initialize the Stokes Problem           initialize the Stokes Problem
842    
843           @param domain: domain of the problem. The approximation order needs to be two.           The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
844           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
845       @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.           with macro elements for the pressure.
846       @type adaptSubTolerance: C{bool}  
847           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           :param domain: domain of the problem.
848             :type domain: `Domain`
849           """           """
850           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
851           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
852           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())
853           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
854            
# Line 409  class StokesProblemCartesian(Homogeneous Line 865  class StokesProblemCartesian(Homogeneous
865           """           """
866       returns the solver options used  solve the equation for velocity.       returns the solver options used  solve the equation for velocity.
867            
868       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
869       """       """
870       return self.__pde_u.getSolverOptions()       return self.__pde_u.getSolverOptions()
871       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
872           """           """
873       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
874            
875       @param options: new solver  options       :param options: new solver  options
876       @type options: L{SolverOptions}       :type options: `SolverOptions`
877       """       """
878           self.__pde_u.setSolverOptions(options)           self.__pde_u.setSolverOptions(options)
879       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
880           """           """
881       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
882       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
883       """       """
884       return self.__pde_prec.getSolverOptions()       return self.__pde_prec.getSolverOptions()
885       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
886           """           """
887       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
888       @param options: new solver  options       :param options: new solver  options
889       @type options: L{SolverOptions}       :type options: `SolverOptions`
890       """       """
891       self.__pde_prec.setSolverOptions(options)       self.__pde_prec.setSolverOptions(options)
892    
# Line 439  class StokesProblemCartesian(Homogeneous Line 895  class StokesProblemCartesian(Homogeneous
895       set the solver options for solving the equation to project the divergence of       set the solver options for solving the equation to project the divergence of
896       the velocity onto the function space of presure.       the velocity onto the function space of presure.
897            
898       @param options: new solver options       :param options: new solver options
899       @type options: L{SolverOptions}       :type options: `SolverOptions`
900       """       """
901       self.__pde_prec.setSolverOptions(options)       self.__pde_proj.setSolverOptions(options)
902       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
903           """           """
904       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
905       the velocity onto the function space of presure.       the velocity onto the function space of presure.
906            
907       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
908       """       """
909       return self.__pde_prec.getSolverOptions()       return self.__pde_proj.getSolverOptions()
910       def setSubProblemTolerance(self):  
911         def updateStokesEquation(self, v, p):
912           """           """
913       Updates the tolerance for subproblems           updates the Stokes equation to consider dependencies from ``v`` and ``p``
914             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
915           """           """
916       if self.adaptSubTolerance():           pass
917               sub_tol=self.getSubProblemTolerance()       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
918           self.getSolverOptionsDiv().setTolerance(sub_tol)          """
919           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)          assigns new values to the model parameters.
920           self.getSolverOptionsPressure().setTolerance(sub_tol)  
921           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)          :param f: external force
922           self.getSolverOptionsVelocity().setTolerance(sub_tol)          :type f: `Vector` object in `FunctionSpace` `Function` or similar
923           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)          :param fixed_u_mask: mask of locations with fixed velocity.
924                    :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
925            :param eta: viscosity
926            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
927            :param surface_stress: normal surface stress
928            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
929            :param stress: initial stress
930        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
931            """
932            if eta !=None:
933                k=util.kronecker(self.domain.getDim())
934                kk=util.outer(k,k)
935                self.eta=util.interpolate(eta, escript.Function(self.domain))
936            self.__pde_prec.setValue(D=1/self.eta)
937                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
938            if restoration_factor!=None:
939                n=self.domain.getNormal()
940                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
941            if fixed_u_mask!=None:
942                self.__pde_u.setValue(q=fixed_u_mask)
943            if f!=None: self.__f=f
944            if surface_stress!=None: self.__surface_stress=surface_stress
945            if stress!=None: self.__stress=stress
946    
947       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
948          """          """
949          assigns values to the model parameters          assigns values to the model parameters
950    
951          @param f: external force          :param f: external force
952          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
953          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
954          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
955          @param eta: viscosity          :param eta: viscosity
956          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
957          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
958          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
959          @param stress: initial stress          :param stress: initial stress
960      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
961          @note: All values needs to be set.          """
962          """          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         self.eta=eta  
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
963    
964       def Bv(self,v):       def Bv(self,v,tol):
965           """           """
966           returns inner product of element p and div(v)           returns inner product of element p and div(v)
967    
968           @param p: a pressure increment           :param v: a residual
969           @param v: a residual           :return: inner product of element p and div(v)
970           @return: inner product of element p and div(v)           :rtype: ``float``
971           @rtype: C{float}           """
972           """           self.__pde_proj.setValue(Y=-util.div(v))
973           self.__pde_proj.setValue(Y=-util.div(v))       self.getSolverOptionsDiv().setTolerance(tol)
974           return self.__pde_proj.getSolution()       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
975             out=self.__pde_proj.getSolution()
976             return out
977    
978       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
979           """           """
980           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
981    
982           @param p: a pressure increment           :param p: a pressure increment
983           @param v: a residual           :param Bv: a residual
984           @return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
985           @rtype: C{float}           :rtype: ``float``
986           """           """
987           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
988    
989       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
990           """           """
991           Returns inner product of p0 and p1           Returns inner product of p0 and p1
992    
993           @param p0: a pressure           :param p0: a pressure
994           @param p1: a pressure           :param p1: a pressure
995           @return: inner product of p0 and p1           :return: inner product of p0 and p1
996           @rtype: C{float}           :rtype: ``float``
997           """           """
998           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
999           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
1000           return util.integrate(s0*s1)           return util.integrate(s0*s1)
1001    
1002       def norm_v(self,v):       def norm_v(self,v):
1003           """           """
1004           returns the norm of v           returns the norm of v
1005    
1006           @param v: a velovity           :param v: a velovity
1007           @return: norm of v           :return: norm of v
1008           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1009           """           """
1010           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
1011    
1012       def getV(self, p, v0):  
1013         def getDV(self, p, v, tol):
1014           """           """
1015           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
1016    
1017           @param p: a pressure           :param p: a pressure
1018           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
1019           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
1020           """           """
1021           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.updateStokesEquation(v,p)
1022             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
1023         self.getSolverOptionsVelocity().setTolerance(tol)
1024         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
1025           if self.__stress.isEmpty():           if self.__stress.isEmpty():
1026              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1027           else:           else:
1028              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1029           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
1030           return  out           return  out
1031    
# Line 560  class StokesProblemCartesian(Homogeneous Line 1033  class StokesProblemCartesian(Homogeneous
1033          """          """
1034          Returns Bv (overwrite).          Returns Bv (overwrite).
1035    
1036          @rtype: equal to the type of p          :rtype: equal to the type of p
1037          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1038          """          """
1039          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1040    
1041       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
1042           """           """
1043           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
1044    
1045           @param p: a pressure increment           :param p: a pressure increment
1046           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
1047           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
1048           """           """
1049           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1050           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
1051           return  out           return  out
1052    
1053       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
1054           """           """
1055           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
1056           with accuracy L{self.getSubProblemTolerance()}           with accuracy `self.getSubProblemTolerance()`
1057    
1058           @param v: velocity increment           :param Bv: velocity increment
1059           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1060           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
1061           """           """
1062           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
1063           return self.__pde_prec.getSolution()       self.getSolverOptionsPressure().setTolerance(tol)
1064         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1065             out=self.__pde_prec.getSolution()
1066             return out

Legend:
Removed from v.2549  
changed lines
  Added in v.3452

  ViewVC Help
Powered by ViewVC 1.1.26