/[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 2108 by gross, Fri Nov 28 05:09:23 2008 UTC revision 3499 by gross, Fri Apr 8 00:14:35 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 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-2008 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"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
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  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40    
41    class DarcyFlowVeryNew(object):
42       """
43       solves the problem
44      
45       *u_i+k_{ij}*p_{,j} = g_i*
46       *u_{i,i} = f*
47      
48       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
49      
50       :note: The problem is solved in a stabelized formulation.
51       """
52       def __init__(self, domain, useReduced=False, *args, **kargs):
53          """
54          initializes the Darcy flux problem
55          :param domain: domain of the problem
56          :type domain: `Domain`
57          :param useReduced: uses reduced oreder on flux and pressure
58          :type useReduced: ``bool``
59          :param adaptSubTolerance: switches on automatic subtolerance selection
60          :type adaptSubTolerance: ``bool``
61          :param solveForFlux: if True the solver solves for the flux (do not use!)
62          :type solveForFlux: ``bool``  
63          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
64          :type useVPIteration: ``bool``    
65          """
66          self.domain=domain
67          useVPIteration=False
68          self.useVPIteration=useVPIteration or True
69          #self.useVPIteration=False
70          self.useReduced=useReduced
71          self.verbose=False
72    
73          if self.useVPIteration:
74         # this does not work yet
75             self.__pde_k=LinearPDESystem(domain)
76             self.__pde_k.setSymmetryOn()
77             if self.useReduced: self.__pde_k.setReducedOrderOn()
78      
79             self.__pde_p=LinearSinglePDE(domain)
80             self.__pde_p.setSymmetryOn()
81             if self.useReduced: self.__pde_p.setReducedOrderOn()
82          else:
83             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
84             self.__pde_k.setSymmetryOff()
85            
86             if self.useReduced: self.__pde_k.setReducedOrderOn()
87             C=self.__pde_k.createCoefficient("C")
88             B=self.__pde_k.createCoefficient("B")
89             for i in range(self.domain.getDim()):
90                B[i,i,self.domain.getDim()]=-1
91                C[self.domain.getDim(),i,i]=1
92                C[i,self.domain.getDim(),i]=-0.5
93                B[self.domain.getDim(),i,i]=0.5
94             self.__pde_k.setValue(C=C, B=B)
95          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
96          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
97          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
98          self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
99          self.scale=1.
100       def __L2(self,v):
101             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))  
102       def __inner_GMRES(self,r,s):
103             return util.integrate(util.inner(r,s))
104            
105       def __Aprod_GMRES(self,p):
106          self.__pde_k.setValue(Y=0.5*util.grad(p), X=p*util.kronecker(self.__pde_k.getDomain()) )
107          du=self.__pde_k.getSolution()
108          self.__pde_p.setValue(Y=util.div(du), X=0.5*(du+util.tensor_mult(self.__permeability,util.grad(p))))
109          return self.__pde_p.getSolution()
110            
111       def getSolverOptionsFlux(self):
112          """
113          Returns the solver options used to solve the flux problems
114          
115          *K^{-1} u=F*
116          
117          :return: `SolverOptions`
118          """
119          return self.__pde_k.getSolverOptions()      
120          
121       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
122          """
123          assigns values to model parameters
124    
125          :param f: volumetic sources/sinks
126          :type f: scalar value on the domain (e.g. `escript.Data`)
127          :param g: flux sources/sinks
128          :type g: vector values on the domain (e.g. `escript.Data`)
129          :param location_of_fixed_pressure: mask for locations where pressure is fixed
130          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
131          :param location_of_fixed_flux:  mask for locations where flux is fixed.
132          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
133          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
134          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
135    
136          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
137          :note: at any point on the boundary of the domain the pressure
138                 (``location_of_fixed_pressure`` >0) or the normal component of the
139                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
140                 is along the *x_i* axis.
141    
142          """
143          if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure
144          if location_of_fixed_flux!=None: self.location_of_fixed_flux=location_of_fixed_flux
145          
146          if self.useVPIteration:
147             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)
148             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux)
149          else:
150         if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:
151            q=self.__pde_k.createCoefficient("q")
152            q[self.domain.getDim()]=self.location_of_fixed_pressure
153            q[:self.domain.getDim()]=self.location_of_fixed_flux
154            self.__pde_k.setValue(q=q)
155                
156          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
157          if permeability!=None:
158         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
159             V=util.vol(self.domain)
160         if perm.getRank()==0:
161            perm_inv=(1./perm)
162                self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
163            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
164            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
165         elif perm.getRank()==2:
166            perm_inv=util.inverse(perm)
167                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
168            perm_inv*=(1./self.scale)
169            perm=perm*self.scale
170         else:
171            raise ValueError,"illegal rank of permeability."
172            
173         self.__permeability=perm
174         self.__permeability_inv=perm_inv
175         if self.useVPIteration:
176                self.__pde_k.setValue(D=0.5*self.__permeability_inv)
177            self.__pde_p.setValue(A=0.5*self.__permeability)
178             else:
179                D=self.__pde_k.createCoefficient("D")
180                A=self.__pde_k.createCoefficient("A")
181                D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv
182                A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability
183                self.__pde_k.setValue(A=A, D=D)
184          if g != None:
185        g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
186        if g.isEmpty():
187              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
188        else:
189            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
190            self.__g=g
191          if f !=None:
192         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
193         if f.isEmpty():
194          
195              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
196         else:
197             if f.getRank()>0: raise ValueError,"illegal rank of f."
198             self.__f=f
199            
200       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
201          """
202          solves the problem.
203          
204          The iteration is terminated if the residual norm is less then self.getTolerance().
205    
206          :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.
207          :type u0: vector value on the domain (e.g. `escript.Data`).
208          :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.
209          :type p0: scalar value on the domain (e.g. `escript.Data`).
210          :param verbose: if set some information on iteration progress are printed
211          :type verbose: ``bool``
212          :return: flux and pressure
213          :rtype: ``tuple`` of `escript.Data`.
214    
215          """
216          u0_b=u0*self.location_of_fixed_flux
217          p0_b=p0*self.location_of_fixed_pressure/self.scale
218          f=self.__f-util.div(u0_b)
219          g=self.__g-u0_b - util.tensor_mult(self.__permeability,util.grad(p0_b))
220          self.verbose=verbose
221          if self.useVPIteration:
222            # get u:
223            self.__pde_k.setValue(Y=0.5*util.tensor_mult(self.__permeability_inv,g),X=escript.Data())
224            du=self.__pde_k.getSolution()
225            self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))
226            p=GMRES(self.__pde_p.getSolution(),
227                    self.__Aprod_GMRES,
228                p0_b*0,
229                self.__inner_GMRES,
230                atol=0,
231                rtol=1.e-4,
232                iter_max=100,
233                iter_restart=20, verbose=self.verbose,P_R=None)
234                self.__pde_k.setValue(Y=0.5*( util.tensor_mult(self.__permeability_inv,g) + util.grad(p)) ,
235                                      X=p*util.kronecker(self.__pde_k.getDomain()))
236            u=self.__pde_k.getSolution()
237          else:
238              X=self.__pde_k.createCoefficient("X")
239              Y=self.__pde_k.createCoefficient("Y")
240              Y[:self.domain.getDim()]=0.5*util.tensor_mult(self.__permeability_inv,g)
241              Y[self.domain.getDim()]=f
242              X[self.domain.getDim(),:]=g*0.5
243              self.__pde_k.setValue(X=X, Y=Y)
244              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
245              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
246              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
247              U=self.__pde_k.getSolution()
248              u=U[:self.domain.getDim()]
249              p=U[self.domain.getDim()]
250          # self.__pde_k.getOperator().saveMM("k.mm")
251          u=u0_b+u
252          p=(p0_b+p)*self.scale
253          if self.verbose:
254            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
255            def_p=self.__g-(u+KGp)
256            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
257            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
258            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
259          return u,p
260       def setTolerance(self,rtol=1e-4):
261          """
262          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
263    
264          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
265          
266          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
267          
268          :param rtol: relative tolerance for the pressure
269          :type rtol: non-negative ``float``
270          """
271          if rtol<0:
272         raise ValueError,"Relative tolerance needs to be non-negative."
273          self.__rtol=rtol
274       def getTolerance(self):
275          """
276          returns the relative tolerance
277          :return: current relative tolerance
278          :rtype: ``float``
279          """
280          return self.__rtol
281  class DarcyFlow(object):  class DarcyFlow(object):
282       """
283       solves the problem
284      
285       *u_i+k_{ij}*p_{,j} = g_i*
286       *u_{i,i} = f*
287      
288       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
289      
290       :note: The problem is solved in a least squares formulation.
291       """
292      
293       def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
294          """
295          initializes the Darcy flux problem
296          :param domain: domain of the problem
297          :type domain: `Domain`
298          :param useReduced: uses reduced oreder on flux and pressure
299          :type useReduced: ``bool``
300          :param adaptSubTolerance: switches on automatic subtolerance selection
301          :type adaptSubTolerance: ``bool``
302          :param solveForFlux: if True the solver solves for the flux (do not use!)
303          :type solveForFlux: ``bool``  
304          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
305          :type useVPIteration: ``bool``    
306          """
307          self.domain=domain
308          self.useVPIteration=useVPIteration
309          self.useReduced=useReduced
310          self.weighting_scale=weighting_scale
311          if self.useVPIteration:
312             self.solveForFlux=solveForFlux
313             self.__adaptSubTolerance=adaptSubTolerance
314             self.verbose=False
315            
316             self.__pde_k=LinearPDESystem(domain)
317             self.__pde_k.setSymmetryOn()
318             if self.useReduced: self.__pde_k.setReducedOrderOn()
319      
320             self.__pde_p=LinearSinglePDE(domain)
321             self.__pde_p.setSymmetryOn()
322             if self.useReduced: self.__pde_p.setReducedOrderOn()
323             self.setTolerance()
324             self.setAbsoluteTolerance()
325          else:
326             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
327             self.__pde_k.setSymmetryOn()
328             if self.useReduced: self.__pde_k.setReducedOrderOn()
329             C=self.__pde_k.createCoefficient("C")
330             B=self.__pde_k.createCoefficient("B")
331             for i in range(self.domain.getDim()):
332                C[i,self.domain.getDim(),i]=1
333                B[self.domain.getDim(),i,i]=1
334             self.__pde_k.setValue(C=C, B=B)
335          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
336          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
337          
338       def getSolverOptionsFlux(self):
339          """
340          Returns the solver options used to solve the flux problems
341          
342          *K^{-1} u=F*
343          
344          :return: `SolverOptions`
345          """
346          return self.__pde_k.getSolverOptions()
347          
348       def setSolverOptionsFlux(self, options=None):
349          """
350          Sets the solver options used to solve the flux problems
351          
352          *K^{-1}u=F*
353          
354          If ``options`` is not present, the options are reset to default
355          
356          :param options: `SolverOptions`
357          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
358          """
359          return self.__pde_v.setSolverOptions(options)
360        
361       def getSolverOptionsPressure(self):
362          """
363          Returns the solver options used to solve the pressure problems
364          
365          *(Q^* K Q)p=-Q^*G*
366          
367          :return: `SolverOptions`
368          """
369          return self.__pde_p.getSolverOptions()
370          
371       def setSolverOptionsPressure(self, options=None):
372          """
373          Sets the solver options used to solve the pressure problems
374          
375          *(Q^* K Q)p=-Q^*G*
376          
377          If ``options`` is not present, the options are reset to default
378          
379          :param options: `SolverOptions`
380          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
381          """
382          return self.__pde_p.setSolverOptions(options)
383    
384    
385       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
386          """
387          assigns values to model parameters
388    
389          :param f: volumetic sources/sinks
390          :type f: scalar value on the domain (e.g. `escript.Data`)
391          :param g: flux sources/sinks
392          :type g: vector values on the domain (e.g. `escript.Data`)
393          :param location_of_fixed_pressure: mask for locations where pressure is fixed
394          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
395          :param location_of_fixed_flux:  mask for locations where flux is fixed.
396          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
397          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
398          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
399    
400          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
401          :note: at any point on the boundary of the domain the pressure
402                 (``location_of_fixed_pressure`` >0) or the normal component of the
403                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
404                 is along the *x_i* axis.
405    
406          """
407          if self.useVPIteration:
408             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
409             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)
410          else:
411             q=self.__pde_k.getCoefficient("q")
412             if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
413             if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
414             if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
415             self.__pde_k.setValue(q=q)
416                
417          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
418          if permeability!=None:
419         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
420             V=util.vol(self.domain)
421         if perm.getRank()==0:
422            perm_inv=(1./perm)
423                if self.useVPIteration:
424                  self.scale=1.
425                else:
426                  self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
427    
428            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
429            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
430         elif perm.getRank()==2:
431            perm_inv=util.inverse(perm)
432                if self.useVPIteration:
433                  self.scale=1.
434                else:
435                  self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
436              perm_inv*=(1./self.scale)
437              perm=perm*self.scale
438         else:
439            raise ValueError,"illegal rank of permeability."
440    
441         self.__permeability=perm
442         self.__permeability_inv=perm_inv
443    
444         self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
445             if self.useVPIteration:
446            if  self.solveForFlux:
447               self.__pde_k.setValue(D=self.__permeability_inv)
448            else:
449               self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
450            self.__pde_p.setValue(A=self.__permeability)
451             else:
452                D=self.__pde_k.createCoefficient("D")
453                A=self.__pde_k.createCoefficient("A")
454                D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
455                for i in range(self.domain.getDim()):
456                   for j in range(self.domain.getDim()):
457                     A[i,i,j,j]=self.__l2
458                A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
459                self.__pde_k.setValue(A=A, D=D)
460          if g !=None:
461         g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
462         if g.isEmpty():
463              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
464         else:
465            if not g.getShape()==(self.domain.getDim(),):
466                  raise ValueError,"illegal shape of g"
467            self.__g=g
468          elif permeability!=None:
469                 X
470          if f !=None:
471         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
472         if f.isEmpty():
473              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
474         else:
475             if f.getRank()>0: raise ValueError,"illegal rank of f."
476             self.__f=f
477    
478          
479       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
480          """
481          solves the problem.
482          
483          The iteration is terminated if the residual norm is less then self.getTolerance().
484    
485          :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.
486          :type u0: vector value on the domain (e.g. `escript.Data`).
487          :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.
488          :type p0: scalar value on the domain (e.g. `escript.Data`).
489          :param verbose: if set some information on iteration progress are printed
490          :type verbose: ``bool``
491          :return: flux and pressure
492          :rtype: ``tuple`` of `escript.Data`.
493    
494          :note: The problem is solved as a least squares form
495                 *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
496                 *G^*u+*G^* K Gp=G^*g*
497                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
498          """
499          self.verbose=verbose
500          if self.useVPIteration:
501              return self.__solveVP(u0,p0,max_iter,max_num_corrections)
502          else:
503              X=self.__pde_k.createCoefficient("X")
504              Y=self.__pde_k.createCoefficient("Y")
505              Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
506              rtmp=self.__f * self.__l2 * self.scale
507              for i in range(self.domain.getDim()): X[i,i]=rtmp
508              X[self.domain.getDim(),:]=self.__g*self.scale
509              r=self.__pde_k.createCoefficient("r")
510              r[:self.domain.getDim()]=u0*self.scale
511              r[self.domain.getDim()]=p0
512              self.__pde_k.setValue(X=X, Y=Y, r=r)
513              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
514              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
515              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
516              U=self.__pde_k.getSolution()
517              # self.__pde_k.getOperator().saveMM("k.mm")
518              u=U[:self.domain.getDim()]*(1./self.scale)
519              p=U[self.domain.getDim()]
520              if self.verbose:
521            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
522            def_p=self.__g-(u+KGp)
523            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
524            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
525            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
526              return u,p
527    
528       def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
529          """
530          solves the problem.
531          
532          The iteration is terminated if the residual norm is less than self.getTolerance().
533    
534          :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.
535          :type u0: vector value on the domain (e.g. `escript.Data`).
536          :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.
537          :type p0: scalar value on the domain (e.g. `escript.Data`).
538          :return: flux and pressure
539          :rtype: ``tuple`` of `escript.Data`.
540    
541          :note: The problem is solved as a least squares form
542                 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
543                 *G^*u+*G^* K Gp=G^*g*
544                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
545          """
546          rtol=self.getTolerance()
547          atol=self.getAbsoluteTolerance()
548          self.setSubProblemTolerance()
549          num_corrections=0
550          converged=False
551          norm_r=None
552          
553          # Eliminate the hydrostatic pressure:
554          if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
555          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))        
556          p0=self.__pde_p.getSolution()
557          g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
558          norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5    
559    
560          p=p0*0
561          if self.solveForFlux:
562         v=u0.copy()
563          else:
564         v=self.__getFlux(p, u0, f=self.__f, g=g2)
565    
566          while not converged and norm_g2 > 0:
567         Gp=util.grad(p)
568         KGp=util.tensor_mult(self.__permeability,Gp)
569         if self.verbose:
570            def_p=g2-(v+KGp)
571            def_v=self.__f-util.div(v)
572            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
573            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
574            print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
575            print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
576            print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
577         ATOL=atol+rtol*norm_g2
578         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
579         if norm_r == None or norm_r>ATOL:
580            if num_corrections>max_num_corrections:
581               raise ValueError,"maximum number of correction steps reached."
582          
583            if self.solveForFlux:
584               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
585               v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
586                   self.__Aprod_v,
587                   v,
588                   self.__Msolve_PCG_v,
589                   self.__inner_PCG_v,
590                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
591               p=r[2]
592            else:
593               # initial residual is r=G^*(g2-KGp - v)
594               p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
595                     self.__Aprod_p,
596                     p,
597                     self.__Msolve_PCG_p,
598                     self.__inner_PCG_p,
599                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
600               v=r[1]
601            if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
602            num_corrections+=1
603         else:
604            if self.verbose: print "DarcyFlux: stopping criterium reached."
605            converged=True
606          return v,p+p0
607    
608       def __applWeight(self, v, f=None):
609          # solves L p = f-Dv with p = 0
610          if self.verbose: print "DarcyFlux: Applying weighting operator"
611          if f == None:
612         return -util.div(v)*self.__l2
613          else:
614         return (f-util.div(v))*self.__l2
615       def __getPressure(self, v, p0, g=None):
616          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
617          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
618          if g == None:
619         self.__pde_p.setValue(X=-v, r=p0)
620          else:
621         self.__pde_p.setValue(X=g-v, r=p0)      
622          p=self.__pde_p.getSolution()
623          return p
624    
625       def __Aprod_v(self,dv):
626          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
627          dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
628          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
629          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
630          return ArithmeticTuple(a,b,-dp)
631    
632       def __Msolve_PCG_v(self,r):
633          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
634          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
635          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
636          return self.__pde_k.getSolution()
637    
638       def __inner_PCG_v(self,v,r):
639          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
640          
641       def __Aprod_p(self,dp):
642          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
643          Gdp=util.grad(dp)
644          self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
645          du=self.__pde_k.getSolution()
646          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
647    
648       def __getFlux(self,p, v0, f=None, g=None):
649          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
650          if f!=None:
651         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
652          self.__pde_k.setValue(r=v0)
653          g2=util.tensor_mult(self.__permeability_inv,g)
654          if p == None:
655         self.__pde_k.setValue(Y=g2)
656          else:
657         self.__pde_k.setValue(Y=g2-util.grad(p))
658          return self.__pde_k.getSolution()  
659          
660          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
661       def __Msolve_PCG_p(self,r):
662          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
663          self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
664          return self.__pde_p.getSolution()
665            
666       def __inner_PCG_p(self,p,r):
667           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
668    
669       def __L2(self,v):
670          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
671    
672       def __L2_r(self,v):
673          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
674    
675       def setTolerance(self,rtol=1e-4):
676          """
677          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
678    
679          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
680          
681          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
682          
683          :param rtol: relative tolerance for the pressure
684          :type rtol: non-negative ``float``
685          """
686          if rtol<0:
687         raise ValueError,"Relative tolerance needs to be non-negative."
688          self.__rtol=rtol
689       def getTolerance(self):
690          """
691          returns the relative tolerance
692          :return: current relative tolerance
693          :rtype: ``float``
694          """
695          return self.__rtol
696    
697       def setAbsoluteTolerance(self,atol=0.):
698          """
699          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
700          
701          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
702    
703    
704          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}*.
705    
706          :param atol: absolute tolerance for the pressure
707          :type atol: non-negative ``float``
708          """
709          if atol<0:
710         raise ValueError,"Absolute tolerance needs to be non-negative."
711          self.__atol=atol
712       def getAbsoluteTolerance(self):
713          """
714          returns the absolute tolerance
715          :return: current absolute tolerance
716          :rtype: ``float``
717          """
718          return self.__atol
719       def getSubProblemTolerance(self):
720          """
721          Returns a suitable subtolerance
722          :type: ``float``
723          """
724          return max(util.EPSILON**(0.5),self.getTolerance()**2)
725    
726       def setSubProblemTolerance(self):
727          """
728          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
729          """
730          if self.__adaptSubTolerance:
731         sub_tol=self.getSubProblemTolerance()
732         self.getSolverOptionsFlux().setTolerance(sub_tol)
733         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
734         self.getSolverOptionsPressure().setTolerance(sub_tol)
735         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
736         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
737    
738    
739    class DarcyFlowOld(object):
740      """      """
741      solves the problem      solves the problem
742    
743      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
744      M{u_{i,i} = f}      *u_{i,i} = f*
745    
746      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,
747    
748      @note: The problem is solved in a least squares formulation.      :note: The problem is solved in a least squares formulation.
749      """      """
750    
751      def __init__(self, domain):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
752          """          """
753          initializes the Darcy flux problem          initializes the Darcy flux problem
754          @param domain: domain of the problem          :param domain: domain of the problem
755          @type domain: L{Domain}          :type domain: `Domain`
756        :param useReduced: uses reduced oreder on flux and pressure
757        :type useReduced: ``bool``
758        :param adaptSubTolerance: switches on automatic subtolerance selection
759        :type adaptSubTolerance: ``bool``  
760          """          """
761          self.domain=domain          self.domain=domain
762            if weight == None:
763               s=self.domain.getSize()
764               self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
765               # self.__l2=(3.*util.longestEdge(self.domain))**2
766               #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
767            else:
768               self.__l2=weight
769          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
770          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          if useReduced: self.__pde_v.setReducedOrderOn()
771          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
772            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
773          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
774          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
775          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          if useReduced: self.__pde_p.setReducedOrderOn()
776          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
777            self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
778            self.setTolerance()
779            self.setAbsoluteTolerance()
780        self.__adaptSubTolerance=adaptSubTolerance
781        self.verbose=False
782        def getSolverOptionsFlux(self):
783        """
784        Returns the solver options used to solve the flux problems
785        
786        *(I+D^*D)u=F*
787        
788        :return: `SolverOptions`
789        """
790        return self.__pde_v.getSolverOptions()
791        def setSolverOptionsFlux(self, options=None):
792        """
793        Sets the solver options used to solve the flux problems
794        
795        *(I+D^*D)u=F*
796        
797        If ``options`` is not present, the options are reset to default
798        :param options: `SolverOptions`
799        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
800        """
801        return self.__pde_v.setSolverOptions(options)
802        def getSolverOptionsPressure(self):
803        """
804        Returns the solver options used to solve the pressure problems
805        
806        *(Q^*Q)p=Q^*G*
807        
808        :return: `SolverOptions`
809        """
810        return self.__pde_p.getSolverOptions()
811        def setSolverOptionsPressure(self, options=None):
812        """
813        Sets the solver options used to solve the pressure problems
814        
815        *(Q^*Q)p=Q^*G*
816        
817        If ``options`` is not present, the options are reset to default
818        :param options: `SolverOptions`
819        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
820        """
821        return self.__pde_p.setSolverOptions(options)
822    
823      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
824          """          """
825          assigns values to model parameters          assigns values to model parameters
826    
827          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
828          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `escript.Data`)
829          @param g: flux sources/sinks          :param g: flux sources/sinks
830          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `escript.Data`)
831          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
832          @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`)
833          @param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
834          @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`)
835          @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
836                               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
837                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
838          @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`)
839    
840          @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.
841          @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)
842                 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
843                 is along the M{x_i} axis.                 is along the *x_i* axis.
844          """          """
845          if f !=None:          if f !=None:
846             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
847             if f.isEmpty():             if f.isEmpty():
848                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
849             else:             else:
850                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
851             self.f=f             self.__f=f
852          if g !=None:            if g !=None:
853             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
854             if g.isEmpty():             if g.isEmpty():
855               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
856             else:             else:
857               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
858                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 118  class DarcyFlow(object): Line 875  class DarcyFlow(object):
875             self.__permeability=perm             self.__permeability=perm
876             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
877    
878        def setTolerance(self,rtol=1e-4):
879            """
880            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
881    
882            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
883    
884      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8):          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}*.
885    
886            :param rtol: relative tolerance for the pressure
887            :type rtol: non-negative ``float``
888          """          """
889          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          if rtol<0:
890          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
891          Note that C{g} and C{f} are used, L{setValue}.          self.__rtol=rtol
892                def getTolerance(self):
         @param p: pressure.  
         @type p: scalar value on the domain (e.g. L{Data}).  
         @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.  
         @type fixed_flux: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive float.  
         @return: flux  
         @rtype: L{Data}  
         @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}}  
893          """          """
894          self.__pde_v.setTolerance(tol)          returns the relative tolerance
895          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=boundary_flux)  
896          return self.__pde_v.getSolution()          :return: current relative tolerance
897            :rtype: ``float``
898            """
899            return self.__rtol
900    
901        def setAbsoluteTolerance(self,atol=0.):
902            """
903            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
904    
905      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
906           """  
907            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}*.
908    
909            :param atol: absolute tolerance for the pressure
910            :type atol: non-negative ``float``
911            """
912            if atol<0:
913                raise ValueError,"Absolute tolerance needs to be non-negative."
914            self.__atol=atol
915        def getAbsoluteTolerance(self):
916           """
917           returns the absolute tolerance
918          
919           :return: current absolute tolerance
920           :rtype: ``float``
921           """
922           return self.__atol
923        def getSubProblemTolerance(self):
924        """
925        Returns a suitable subtolerance
926        @type: ``float``
927        """
928        return max(util.EPSILON**(0.75),self.getTolerance()**2)
929        def setSubProblemTolerance(self):
930             """
931             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
932             """
933         if self.__adaptSubTolerance:
934             sub_tol=self.getSubProblemTolerance()
935                 self.getSolverOptionsFlux().setTolerance(sub_tol)
936             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
937             self.getSolverOptionsPressure().setTolerance(sub_tol)
938             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
939             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
940    
941        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
942             """
943           solves the problem.           solves the problem.
   
          The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where  
          C{|q|} denotes the norm of the right hand side (see escript user's guide for details).  
944    
945           @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.           The iteration is terminated if the residual norm is less then self.getTolerance().
946           @type u0: vector value on the domain (e.g. L{Data}).  
947           @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 u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
948           @type p0: scalar value on the domain (e.g. L{Data}).           :type u0: vector value on the domain (e.g. `escript.Data`).
949           @param atol: absolute tolerance for the pressure           :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.
950           @type atol: non-negative C{float}           :type p0: scalar value on the domain (e.g. `escript.Data`).
951           @param rtol: relative tolerance for the pressure           :param verbose: if set some information on iteration progress are printed
952           @type rtol: non-negative C{float}           :type verbose: ``bool``
953           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}           :return: flux and pressure
954           @type sub_rtol: positive-negative C{float}           :rtype: ``tuple`` of `escript.Data`.
955           @param verbose: if set some information on iteration progress are printed  
956           @type verbose: C{bool}           :note: The problem is solved as a least squares form
          @param show_details:  if set information on the subiteration process are printed.  
          @type show_details: C{bool}  
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
   
          @note: The problem is solved as a least squares form  
957    
958           M{(I+D^*D)u+Qp=D^*f+g}           *(I+D^*D)u+Qp=D^*f+g*
959           M{Q^*u+Q^*Qp=Q^*g}           *Q^*u+Q^*Qp=Q^*g*
960    
961           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}*.
962           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
963    
964           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
965    
966           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
967    
968           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
969            
970           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
971           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.
972           """           """
973           self.verbose=verbose           self.verbose=verbose
974           self.show_details= show_details and self.verbose           rtol=self.getTolerance()
975           self.__pde_v.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
976           self.__pde_p.setTolerance(sub_rtol)       self.setSubProblemTolerance()
977           p2=p0*self.__pde_p.getCoefficient("q")           num_corrections=0
978           u2=u0*self.__pde_v.getCoefficient("q")           converged=False
979           g=self.__g-u2-util.tensor_mult(self.__permeability,util.grad(p2))           p=p0
980           f=self.__f-util.div(u2)           norm_r=None
981           self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())           while not converged:
982           dv=self.__pde_v.getSolution(verbose=show_details)                 v=self.getFlux(p, fixed_flux=u0)
983           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))                 Qp=self.__Q(p)
984           self.__pde_p.setValue(r=Data())                 norm_v=self.__L2(v)
985           dp=self.__pde_p.getSolution(verbose=self.show_details)                 norm_Qp=self.__L2(Qp)
986           norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))                 if norm_v == 0.:
987           if norm_rhs<0:                    if norm_Qp == 0.:
988               raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."                       return v,p
989           ATOL=util.sqrt(norm_rhs)*rtol +atol                    else:
990           if not ATOL>0:                      fac=norm_Qp
991               raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)                 else:
992           rhs=ArithmeticTuple(g,dv)                    if norm_Qp == 0.:
993           dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0-p2, verbose=self.verbose, initial_guess=True)                      fac=norm_v
994           return u2+r[1],p2+dp                    else:
995                                fac=2./(1./norm_v+1./norm_Qp)
996      def __Aprod_PCG(self,p):                 ATOL=(atol+rtol*fac)
997            if self.show_details: print "DarcyFlux: Applying operator"                 if self.verbose:
998            Qp=util.tensor_mult(self.__permeability,util.grad(p))                      print "DarcyFlux: L2 norm of v = %e."%norm_v
999            self.__pde_v.setValue(Y=Qp,X=Data())                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
1000            w=self.__pde_v.getSolution(verbose=self.show_details)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
1001            return ArithmeticTuple(Qp,w)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
1002                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
1003                   if norm_r == None or norm_r>ATOL:
1004                       if num_corrections>max_num_corrections:
1005                             raise ValueError,"maximum number of correction steps reached."
1006                       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)
1007                       num_corrections+=1
1008                   else:
1009                       converged=True
1010             return v,p
1011        def __L2(self,v):
1012             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1013    
1014        def __Q(self,p):
1015              return util.tensor_mult(self.__permeability,util.grad(p))
1016    
1017        def __Aprod(self,dp):
1018              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
1019              Qdp=self.__Q(dp)
1020              self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
1021              du=self.__pde_v.getSolution()
1022              return Qdp+du
1023        def __inner_GMRES(self,r,s):
1024             return util.integrate(util.inner(r,s))
1025    
1026      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
1027           a=util.tensor_mult(self.__permeability,util.grad(p))           return util.integrate(util.inner(self.__Q(p), r))
          return util.integrate(util.inner(a,r[0]-r[1]))  
1028    
1029      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
1030            if self.show_details: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
1031            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]))            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
1032            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution()
1033    
1034        def getFlux(self,p=None, fixed_flux=escript.Data()):
1035            """
1036            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
1037            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
1038            Note that ``g`` and ``f`` are used, see `setValue`.
1039    
1040            :param p: pressure.
1041            :type p: scalar value on the domain (e.g. `escript.Data`).
1042            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
1043            :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
1044            :return: flux
1045            :rtype: `escript.Data`
1046            :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}*
1047                   for the permeability *k_{ij}*
1048            """
1049        self.setSubProblemTolerance()
1050            g=self.__g
1051            f=self.__f
1052            self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
1053            if p == None:
1054               self.__pde_v.setValue(Y=g)
1055            else:
1056               self.__pde_v.setValue(Y=g-self.__Q(p))
1057            return self.__pde_v.getSolution()
1058    
1059  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
1060        """       """
1061        solves       solves
1062    
1063            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
1064                  u_{i,i}=0                  u_{i,i}=0
# Line 230  class StokesProblemCartesian(Homogeneous Line 1066  class StokesProblemCartesian(Homogeneous
1066            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
1067            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
1068    
1069        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
1070    
1071        typical usage:       typical usage:
1072    
1073              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
1074              sp.setTolerance()              sp.setTolerance()
1075              sp.initialize(...)              sp.initialize(...)
1076              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
1077        """       """
1078        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
1079           """           """
1080           initialize the Stokes Problem           initialize the Stokes Problem
1081    
1082           @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
1083           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
1084           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
1085    
1086             :param domain: domain of the problem.
1087             :type domain: `Domain`
1088           """           """
1089           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
1090           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
1091           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())
1092           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
1093           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
1094           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
1095           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
1096           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
1097    
1098           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
1099           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
1100         self.__pde_proj.setValue(D=1)
1101           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
1102    
1103        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
1104             """
1105         returns the solver options used  solve the equation for velocity.
1106        
1107         :rtype: `SolverOptions`
1108         """
1109         return self.__pde_u.getSolverOptions()
1110         def setSolverOptionsVelocity(self, options=None):
1111             """
1112         set the solver options for solving the equation for velocity.
1113        
1114         :param options: new solver  options
1115         :type options: `SolverOptions`
1116         """
1117             self.__pde_u.setSolverOptions(options)
1118         def getSolverOptionsPressure(self):
1119             """
1120         returns the solver options used  solve the equation for pressure.
1121         :rtype: `SolverOptions`
1122         """
1123         return self.__pde_prec.getSolverOptions()
1124         def setSolverOptionsPressure(self, options=None):
1125             """
1126         set the solver options for solving the equation for pressure.
1127         :param options: new solver  options
1128         :type options: `SolverOptions`
1129         """
1130         self.__pde_prec.setSolverOptions(options)
1131    
1132         def setSolverOptionsDiv(self, options=None):
1133             """
1134         set the solver options for solving the equation to project the divergence of
1135         the velocity onto the function space of presure.
1136        
1137         :param options: new solver options
1138         :type options: `SolverOptions`
1139         """
1140         self.__pde_proj.setSolverOptions(options)
1141         def getSolverOptionsDiv(self):
1142             """
1143         returns the solver options for solving the equation to project the divergence of
1144         the velocity onto the function space of presure.
1145        
1146         :rtype: `SolverOptions`
1147         """
1148         return self.__pde_proj.getSolverOptions()
1149    
1150         def updateStokesEquation(self, v, p):
1151             """
1152             updates the Stokes equation to consider dependencies from ``v`` and ``p``
1153             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
1154             """
1155             pass
1156         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
1157            """
1158            assigns new values to the model parameters.
1159    
1160            :param f: external force
1161            :type f: `Vector` object in `FunctionSpace` `Function` or similar
1162            :param fixed_u_mask: mask of locations with fixed velocity.
1163            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
1164            :param eta: viscosity
1165            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
1166            :param surface_stress: normal surface stress
1167            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
1168            :param stress: initial stress
1169        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
1170            """
1171            if eta !=None:
1172                k=util.kronecker(self.domain.getDim())
1173                kk=util.outer(k,k)
1174                self.eta=util.interpolate(eta, escript.Function(self.domain))
1175            self.__pde_prec.setValue(D=1/self.eta)
1176                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
1177            if restoration_factor!=None:
1178                n=self.domain.getNormal()
1179                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
1180            if fixed_u_mask!=None:
1181                self.__pde_u.setValue(q=fixed_u_mask)
1182            if f!=None: self.__f=f
1183            if surface_stress!=None: self.__surface_stress=surface_stress
1184            if stress!=None: self.__stress=stress
1185    
1186         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
1187          """          """
1188          assigns values to the model parameters          assigns values to the model parameters
1189    
1190          @param f: external force          :param f: external force
1191          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
1192          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
1193          @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
1194          @param eta: viscosity          :param eta: viscosity
1195          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
1196          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
1197          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
1198          @param stress: initial stress          :param stress: initial stress
1199      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
         @note: All values needs to be set.  
   
1200          """          """
1201          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         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,Y=f,y=surface_stress)  
         self.__stress=stress  
   
       def B(self,v):  
         """  
         returns div(v)  
         @rtype: equal to the type of p  
   
         @note: boundary conditions on p should be zero!  
         """  
         if self.show_details: print "apply divergence:"  
         self.__pde_proj.setValue(Y=-util.div(v))  
         self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
         return self.__pde_proj.getSolution(verbose=self.show_details)  
1202    
1203        def inner_pBv(self,p,Bv):       def Bv(self,v,tol):
1204           """           """
1205           returns inner product of element p and Bv  (overwrite)           returns inner product of element p and div(v)
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @rtype: C{float}  
1206    
1207           @rtype: equal to the type of p           :param v: a residual
1208             :return: inner product of element p and div(v)
1209             :rtype: ``float``
1210           """           """
1211           s0=util.interpolate(p,Function(self.domain))           self.__pde_proj.setValue(Y=-util.div(v))
1212           s1=util.interpolate(Bv,Function(self.domain))       self.getSolverOptionsDiv().setTolerance(tol)
1213           return util.integrate(s0*s1)       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
1214             out=self.__pde_proj.getSolution()
1215             return out
1216    
1217        def inner_p(self,p0,p1):       def inner_pBv(self,p,Bv):
1218           """           """
1219           returns inner product of element p0 and p1  (overwrite)           returns inner product of element p and Bv=-div(v)
1220            
1221           @type p0: equal to the type of p           :param p: a pressure increment
1222           @type p1: equal to the type of p           :param Bv: a residual
1223           @rtype: C{float}           :return: inner product of element p and Bv=-div(v)
1224             :rtype: ``float``
1225             """
1226             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
1227    
1228           @rtype: equal to the type of p       def inner_p(self,p0,p1):
1229           """           """
1230           s0=util.interpolate(p0/self.eta,Function(self.domain))           Returns inner product of p0 and p1
1231           s1=util.interpolate(p1/self.eta,Function(self.domain))  
1232             :param p0: a pressure
1233             :param p1: a pressure
1234             :return: inner product of p0 and p1
1235             :rtype: ``float``
1236             """
1237             s0=util.interpolate(p0, escript.Function(self.domain))
1238             s1=util.interpolate(p1, escript.Function(self.domain))
1239           return util.integrate(s0*s1)           return util.integrate(s0*s1)
1240    
1241        def inner_v(self,v0,v1):       def norm_v(self,v):
1242           """           """
1243           returns inner product of two element v0 and v1  (overwrite)           returns the norm of v
           
          @type v0: equal to the type of v  
          @type v1: equal to the type of v  
          @rtype: C{float}  
1244    
1245           @rtype: equal to the type of v           :param v: a velovity
1246             :return: norm of v
1247             :rtype: non-negative ``float``
1248           """           """
1249       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
1250    
1251        def solve_A(self,u,p):  
1252         def getDV(self, p, v, tol):
1253           """           """
1254           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
1255    
1256             :param p: a pressure
1257             :param v: a initial guess for the value v to return.
1258             :return: dv given as *Adv=(f-Av-B^*p)*
1259           """           """
1260           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
1261           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
1262         self.getSolverOptionsVelocity().setTolerance(tol)
1263         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
1264           if self.__stress.isEmpty():           if self.__stress.isEmpty():
1265              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1266           else:           else:
1267              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+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)))
1268           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
1269             return  out
1270    
1271         def norm_Bv(self,Bv):
1272            """
1273            Returns Bv (overwrite).
1274    
1275            :rtype: equal to the type of p
1276            :note: boundary conditions on p should be zero!
1277            """
1278            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1279    
1280         def solve_AinvBt(self,p, tol):
1281             """
1282             Solves *Av=B^*p* with accuracy `tol`
1283    
1284             :param p: a pressure increment
1285             :return: the solution of *Av=B^*p*
1286             :note: boundary conditions on v should be zero!
1287             """
1288             self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1289             out=self.__pde_u.getSolution()
1290           return  out           return  out
1291    
1292        def solve_prec(self,p):       def solve_prec(self,Bv, tol):
1293           if self.show_details: print "apply preconditioner:"           """
1294           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           applies preconditioner for for *BA^{-1}B^** to *Bv*
1295           self.__pde_prec.setValue(Y=p)           with accuracy `self.getSubProblemTolerance()`
1296           q=self.__pde_prec.getSolution(verbose=self.show_details)  
1297           return q           :param Bv: velocity increment
1298             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1299             :note: boundary conditions on p are zero.
1300             """
1301             self.__pde_prec.setValue(Y=Bv)
1302         self.getSolverOptionsPressure().setTolerance(tol)
1303         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1304             out=self.__pde_prec.getSolution()
1305             return out

Legend:
Removed from v.2108  
changed lines
  Added in v.3499

  ViewVC Help
Powered by ViewVC 1.1.26