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

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

Legend:
Removed from v.2208  
changed lines
  Added in v.3501

  ViewVC Help
Powered by ViewVC 1.1.26