/[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 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 3500 by gross, Mon Apr 11 03:09:06 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, NegativeNorm  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, useVPIteration=True, *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          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
61          :type useVPIteration: ``bool``    
62          """
63          self.domain=domain
64          self.useVPIteration=useVPIteration
65          self.useReduced=useReduced
66          self.verbose=False
67    
68          if self.useVPIteration:
69         # this does not work yet
70             self.__pde_k=LinearPDESystem(domain)
71             self.__pde_k.setSymmetryOn()
72             if self.useReduced: self.__pde_k.setReducedOrderOn()
73      
74             self.__pde_p=LinearSinglePDE(domain)
75             self.__pde_p.setSymmetryOn()
76             if self.useReduced: self.__pde_p.setReducedOrderOn()
77          else:
78             self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
79             self.__pde_k.setSymmetryOff()
80            
81             if self.useReduced: self.__pde_k.setReducedOrderOn()
82             C=self.__pde_k.createCoefficient("C")
83             B=self.__pde_k.createCoefficient("B")
84             for i in range(self.domain.getDim()):
85                B[i,i,self.domain.getDim()]=-1
86                C[self.domain.getDim(),i,i]=1
87                C[i,self.domain.getDim(),i]=-0.5
88                B[self.domain.getDim(),i,i]=0.5
89             self.__pde_k.setValue(C=C, B=B)
90          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
91          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
92          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
93          self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))
94          self.scale=1.
95          self.setTolerance()
96          
97    
98       def __L2(self,v):
99             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))  
100       def __inner_GMRES(self,r,s):
101             return util.integrate(util.inner(r,s))
102            
103       def __Aprod_GMRES(self,p):
104          self.__pde_k.setValue(Y=-0.5*util.grad(p), X=-p*util.kronecker(self.__pde_k.getDomain()) )
105          du=self.__pde_k.getSolution()
106          self.__pde_p.setValue(Y=-util.div(du), X=0.5*(-du+util.tensor_mult(self.__permeability,util.grad(p))))
107          return self.__pde_p.getSolution()
108            
109       def getSolverOptionsFlux(self):
110          """
111          Returns the solver options used to solve the flux problems
112          
113          *K^{-1} u=F*
114          
115          :return: `SolverOptions`
116          """
117          return self.__pde_k.getSolverOptions()      
118          
119       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
120          """
121          assigns values to model parameters
122    
123          :param f: volumetic sources/sinks
124          :type f: scalar value on the domain (e.g. `escript.Data`)
125          :param g: flux sources/sinks
126          :type g: vector values on the domain (e.g. `escript.Data`)
127          :param location_of_fixed_pressure: mask for locations where pressure is fixed
128          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
129          :param location_of_fixed_flux:  mask for locations where flux is fixed.
130          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
131          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
132          :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
133    
134          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
135          :note: at any point on the boundary of the domain the pressure
136                 (``location_of_fixed_pressure`` >0) or the normal component of the
137                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
138                 is along the *x_i* axis.
139    
140          """
141          if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure
142          if location_of_fixed_flux!=None: self.location_of_fixed_flux=location_of_fixed_flux
143          
144          if self.useVPIteration:
145             if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)
146             if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux)
147          else:
148         if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:
149            q=self.__pde_k.createCoefficient("q")
150            q[self.domain.getDim()]=self.location_of_fixed_pressure
151            q[:self.domain.getDim()]=self.location_of_fixed_flux
152            self.__pde_k.setValue(q=q)
153                
154          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
155          if permeability!=None:
156         perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
157             V=util.vol(self.domain)
158         if perm.getRank()==0:
159            perm_inv=(1./perm)
160                self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
161            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
162            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
163         elif perm.getRank()==2:
164            perm_inv=util.inverse(perm)
165                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
166            perm_inv*=(1./self.scale)
167            perm=perm*self.scale
168         else:
169            raise ValueError,"illegal rank of permeability."
170            
171         self.__permeability=perm
172         self.__permeability_inv=perm_inv
173         if self.useVPIteration:
174                self.__pde_k.setValue(D=0.5*self.__permeability_inv)
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176             else:
177                D=self.__pde_k.createCoefficient("D")
178                A=self.__pde_k.createCoefficient("A")
179                D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv
180                A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability
181                self.__pde_k.setValue(A=A, D=D)
182          if g != None:
183        g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
184        if g.isEmpty():
185              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
186        else:
187            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
188            self.__g=g
189          if f !=None:
190         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
191         if f.isEmpty():
192          
193              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
194         else:
195             if f.getRank()>0: raise ValueError,"illegal rank of f."
196             self.__f=f
197            
198       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10, iter_restart=20):
199          """
200          solves the problem.
201          
202          The iteration is terminated if the residual norm is less then self.getTolerance().
203    
204          :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.
205          :type u0: vector value on the domain (e.g. `escript.Data`).
206          :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.
207          :type p0: scalar value on the domain (e.g. `escript.Data`).
208          :param verbose: if set some information on iteration progress are printed
209          :type verbose: ``bool``
210          :return: flux and pressure
211          :rtype: ``tuple`` of `escript.Data`.
212    
213          """
214          u0_b=u0*self.location_of_fixed_flux
215          p0_b=p0*self.location_of_fixed_pressure/self.scale
216          f=self.__f-util.div(u0_b)
217          g=self.__g-u0_b - util.tensor_mult(self.__permeability,util.grad(p0_b))
218          self.verbose=verbose
219          if self.useVPIteration:
220            # get u:
221                rtol=self.getTolerance()
222            self.__pde_k.setValue(Y=0.5*util.tensor_mult(self.__permeability_inv,g),X=escript.Data())
223            du=self.__pde_k.getSolution()
224            self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))
225            p=GMRES(self.__pde_p.getSolution(),
226                    self.__Aprod_GMRES,
227                p0_b*0,
228                self.__inner_GMRES,
229                atol=0,
230                rtol=rtol,
231                iter_max=max_iter,
232                iter_restart=iter_restart, verbose=self.verbose,P_R=None)
233                self.__pde_k.setValue(Y=0.5*( util.tensor_mult(self.__permeability_inv,g) + util.grad(p)) ,
234                                      X=p*util.kronecker(self.__pde_k.getDomain()))
235                #self.__pde_k.setValue(Y=0.5*util.grad(p), 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,useReduced=False):      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          if useReduced: self.__pde_v.setReducedOrderOn()          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=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)))
773          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
774          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
775          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
776          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
777          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
778          self.__ATOL= None          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 121  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):
     def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):  
         """  
         returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}  
         on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
         Note that C{g} and C{f} are used, see L{setValue}.  
           
         @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 C{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}}  
879          """          """
880          self.__pde_v.setTolerance(tol)          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
         g=self.__g  
         f=self.__f  
         self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
         else:  
            self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))  
         return self.__pde_v.getSolution(verbose=show_details)  
881    
882      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
883    
884            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 pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
890          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
891          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
892                def getTolerance(self):
         @param v: flux.  
         @type v: vector-valued on the domain (e.g. L{Data}).  
         @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.  
         @type fixed_pressure: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: pressure  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where 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
         g=self.__g  
         self.__pde_p.setValue(r=fixed_pressure)  
         if v == None:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))  
         else:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))  
         return self.__pde_p.getSolution(verbose=show_details)  
895    
896      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          :return: current relative tolerance
897            :rtype: ``float``
898          """          """
899          set the tolerance C{ATOL} used to terminate the solution process. It is used          return self.__rtol
900    
901          M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}      def setAbsoluteTolerance(self,atol=0.):
902            """
903          where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
   
         The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual  
   
         M{r=Q^*(g-Qp-v)}  
   
         the condition  
904    
905          M{<(Q^*Q)^{-1} r,r> <= ATOL}          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
906    
907          holds. M{D} is the M{div} operator 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}*.
908    
909          @param atol: absolute tolerance for the pressure          :param atol: absolute tolerance for the pressure
910          @type atol: non-negative C{float}          :type atol: non-negative ``float``
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{float}  
         @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the  
         L{getPressure} method or use  the value from a previous time step.  
         @type p_ref: scalar value on the domain (e.g. L{Data}).  
         @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the  
         L{getFlux} method or use  the value from a previous time step.  
         @type v_ref: vector-valued on the domain (e.g. L{Data}).  
         @return: used absolute tolerance.  
         @rtype: positive C{float}  
         """  
         g=self.__g  
         if not v_ref == None:  
            f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)  
         else:  
            f1=util.integrate(util.length(util.interpolate(g))**2)  
         if not p_ref == None:  
            f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)  
         else:  
            f2=0  
         self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))  
         if self.__ATOL<=0:  
            raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL  
         return self.__ATOL  
           
     def getTolerance(self):  
911          """          """
912          returns the current tolerance.          if atol<0:
913                  raise ValueError,"Absolute tolerance needs to be non-negative."
914          @return: used absolute tolerance.          self.__atol=atol
915          @rtype: positive C{float}      def getAbsoluteTolerance(self):
916          """         """
917          if self.__ATOL==None:         returns the absolute tolerance
918             raise ValueError,"no tolerance is defined."        
919          return self.__ATOL         :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, show_details=False, sub_rtol=1.e-8):      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
942           """           """
943           solves the problem.           solves the problem.
944    
945           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().
946    
947           @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.
948           @type u0: vector value on the domain (e.g. L{Data}).           :type u0: vector value on the domain (e.g. `escript.Data`).
949           @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.
950           @type p0: scalar value on the domain (e.g. L{Data}).           :type p0: scalar value on the domain (e.g. `escript.Data`).
951           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}           :param verbose: if set some information on iteration progress are printed
952           @type sub_rtol: positive-negative C{float}           :type verbose: ``bool``
953           @param verbose: if set some information on iteration progress are printed           :return: flux and pressure
954           @type verbose: C{bool}           :rtype: ``tuple`` of `escript.Data`.
955           @param show_details:  if set information on the subiteration process are printed.  
956           @type show_details: C{bool}           :note: The problem is solved as a least squares form
          @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           ATOL=self.getTolerance()           num_corrections=0
978           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL           converged=False
979           #########################################################################################################################           p=p0
980           #           norm_r=None
981           #   we solve:           while not converged:
982           #                   v=self.getFlux(p, fixed_flux=u0)
983           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )                 Qp=self.__Q(p)
984           #                 norm_v=self.__L2(v)
985           #   residual is                 norm_Qp=self.__L2(Qp)
986           #                 if norm_v == 0.:
987           #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)                    if norm_Qp == 0.:
988           #                       return v,p
989           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                    else:
990           #                      fac=norm_Qp
991           #    we use (g - Qp, v) to represent the residual. not that                 else:
992           #                    if norm_Qp == 0.:
993           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                      fac=norm_v
994           #                    else:
995           #   while the initial residual is                      fac=2./(1./norm_v+1./norm_Qp)
996           #                 ATOL=(atol+rtol*fac)
997           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                 if self.verbose:
998           #                        print "DarcyFlux: L2 norm of v = %e."%norm_v
999           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
1000           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
1001           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
1002           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
1003           self.__pde_v.setValue(r=Data())                 if norm_r == None or norm_r>ATOL:
1004           # start CG                     if num_corrections>max_num_corrections:
1005           r=ArithmeticTuple(d0, v00)                           raise ValueError,"maximum number of correction steps reached."
1006           p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=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)
1007           return r[1],p                     num_corrections+=1
1008                   else:
1009      def __Aprod_PCG(self,dp):                     converged=True
1010            if self.show_details: print "DarcyFlux: Applying operator"           return v,p
1011            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)      def __L2(self,v):
1012            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1013            self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())  
1014            du=self.__pde_v.getSolution(verbose=self.show_details)      def __Q(self,p):
1015            return ArithmeticTuple(mQdp,du)            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))
          f0=util.integrate(util.inner(a,r[0]))  
          f1=util.integrate(util.inner(a,r[1]))  
          # print "__inner_PCG:",f0,f1,"->",f0-f1  
          return f0-f1  
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]), r=Data())            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 333  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         def inner_p(self,p0,p1):
1229             """
1230             Returns inner product of p0 and p1
1231    
1232           @rtype: equal to the type of p           :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/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
1238           s1=util.interpolate(p1/self.eta,Function(self.domain))           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))
1250       gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
1251    
1252        def solve_A(self,u,p):       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.2208  
changed lines
  Added in v.3500

  ViewVC Help
Powered by ViewVC 1.1.26