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

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

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

revision 2549 by jfenwick, Mon Jul 20 06:43:47 2009 UTC revision 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-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 21  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40    
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, weight=None, useReduced=False, adaptSubTolerance=True):      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      :param useReduced: uses reduced oreder on flux and pressure
762      @type useReduced: C{bool}      :type useReduced: ``bool``
763      @param adaptSubTolerance: switches on automatic subtolerance selection      :param adaptSubTolerance: switches on automatic subtolerance selection
764      @type adaptSubTolerance: C{bool}          :type adaptSubTolerance: ``bool``  
765          """          """
766          self.domain=domain          self.domain=domain
767          if weight == None:          if weight == None:
768             s=self.domain.getSize()             s=self.domain.getSize()
769             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
770             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l2=(3.*util.longestEdge(self.domain))**2
771             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
772          else:          else:
773             self.__l=weight             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=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
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.setTolerance()          self.setTolerance()
784          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
785      self.__adaptSubTolerance=adaptSubTolerance      self.__adaptSubTolerance=adaptSubTolerance
# Line 83  class DarcyFlow(object): Line 788  class DarcyFlow(object):
788      """      """
789      Returns the solver options used to solve the flux problems      Returns the solver options used to solve the flux problems
790            
791      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
792            
793      @return: L{SolverOptions}      :return: `SolverOptions`
794      """      """
795      return self.__pde_v.getSolverOptions()      return self.__pde_v.getSolverOptions()
796      def setSolverOptionsFlux(self, options=None):      def setSolverOptionsFlux(self, options=None):
797      """      """
798      Sets the solver options used to solve the flux problems      Sets the solver options used to solve the flux problems
799            
800      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
801            
802      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
803      @param options: L{SolverOptions}      :param options: `SolverOptions`
804      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
805      """      """
806      return self.__pde_v.setSolverOptions(options)      return self.__pde_v.setSolverOptions(options)
807      def getSolverOptionsPressure(self):      def getSolverOptionsPressure(self):
808      """      """
809      Returns the solver options used to solve the pressure problems      Returns the solver options used to solve the pressure problems
810            
811      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
812            
813      @return: L{SolverOptions}      :return: `SolverOptions`
814      """      """
815      return self.__pde_p.getSolverOptions()      return self.__pde_p.getSolverOptions()
816      def setSolverOptionsPressure(self, options=None):      def setSolverOptionsPressure(self, options=None):
817      """      """
818      Sets the solver options used to solve the pressure problems      Sets the solver options used to solve the pressure problems
819            
820      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
821            
822      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
823      @param options: L{SolverOptions}      :param options: `SolverOptions`
824      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
825      """      """
826      return self.__pde_p.setSolverOptions(options)      return self.__pde_p.setSolverOptions(options)
827    
# Line 124  class DarcyFlow(object): Line 829  class DarcyFlow(object):
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 177  class DarcyFlow(object): Line 882  class DarcyFlow(object):
882    
883      def setTolerance(self,rtol=1e-4):      def setTolerance(self,rtol=1e-4):
884          """          """
885          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
886    
887          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
888    
889          where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
890    
891          @param rtol: relative tolerance for the pressure          :param rtol: relative tolerance for the pressure
892          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
893          """          """
894          if rtol<0:          if rtol<0:
895              raise ValueError,"Relative tolerance needs to be non-negative."              raise ValueError,"Relative tolerance needs to be non-negative."
# Line 193  class DarcyFlow(object): Line 898  class DarcyFlow(object):
898          """          """
899          returns the relative tolerance          returns the relative tolerance
900    
901          @return: current relative tolerance          :return: current relative tolerance
902          @rtype: C{float}          :rtype: ``float``
903          """          """
904          return self.__rtol          return self.__rtol
905    
906      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
907          """          """
908          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
909    
910          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
911    
912          where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
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``
916          """          """
917          if atol<0:          if atol<0:
918              raise ValueError,"Absolute tolerance needs to be non-negative."              raise ValueError,"Absolute tolerance needs to be non-negative."
# Line 216  class DarcyFlow(object): Line 921  class DarcyFlow(object):
921         """         """
922         returns the absolute tolerance         returns the absolute tolerance
923                
924         @return: current absolute tolerance         :return: current absolute tolerance
925         @rtype: C{float}         :rtype: ``float``
926         """         """
927         return self.__atol         return self.__atol
928      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
929      """      """
930      Returns a suitable subtolerance      Returns a suitable subtolerance
931      @type: C{float}      @type: ``float``
932      """      """
933      return max(util.EPSILON**(0.75),self.getTolerance()**2)      return max(util.EPSILON**(0.75),self.getTolerance()**2)
934      def setSubProblemTolerance(self):      def setSubProblemTolerance(self):
# Line 244  class DarcyFlow(object): Line 949  class DarcyFlow(object):
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 verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
957           @type verbose: C{bool}           :type verbose: ``bool``
958           @return: flux and pressure           :return: flux and pressure
959           @rtype: C{tuple} of L{Data}.           :rtype: ``tuple`` of `escript.Data`.
960    
961           @note: The problem is solved as a least squares form           :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           rtol=self.getTolerance()           rtol=self.getTolerance()
980           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
981       self.setSubProblemTolerance()       self.setSubProblemTolerance()
       
982           num_corrections=0           num_corrections=0
983           converged=False           converged=False
984           p=p0           p=p0
# Line 297  class DarcyFlow(object): Line 1001  class DarcyFlow(object):
1001                 ATOL=(atol+rtol*fac)                 ATOL=(atol+rtol*fac)
1002                 if self.verbose:                 if self.verbose:
1003                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
1004                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
1005                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
1006                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
1007                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
1008                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
1009                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
1010                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
1011                     p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                     p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1012                     num_corrections+=1                     num_corrections+=1
1013                 else:                 else:
1014                     converged=True                     converged=True
1015           return v,p           return v,p
1016      def __L2(self,v):      def __L2(self,v):
1017           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1018    
1019      def __Q(self,p):      def __Q(self,p):
1020            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
# Line 318  class DarcyFlow(object): Line 1022  class DarcyFlow(object):
1022      def __Aprod(self,dp):      def __Aprod(self,dp):
1023            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
1024            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
1025            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
1026            du=self.__pde_v.getSolution()            du=self.__pde_v.getSolution()
           # self.__pde_v.getOperator().saveMM("proj.mm")  
1027            return Qdp+du            return Qdp+du
1028      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
1029           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 330  class DarcyFlow(object): Line 1033  class DarcyFlow(object):
1033    
1034      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
1035        if self.getSolverOptionsPressure().isVerbose(): 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), Y=Data(), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
           # self.__pde_p.getOperator().saveMM("prec.mm")  
1037            return self.__pde_p.getSolution()            return self.__pde_p.getSolution()
1038    
1039      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
1040          """          """
1041          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
1042          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
1043          Note that C{g} and C{f} are used, see L{setValue}.          Note that ``g`` and ``f`` are used, see `setValue`.
1044    
1045          @param p: pressure.          :param p: pressure.
1046          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `escript.Data`).
1047          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
1048          @type fixed_flux: vector values on the domain (e.g. L{Data}).          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
1049          @param tol: relative tolerance to be used.          :return: flux
1050          @type tol: positive C{float}.          :rtype: `escript.Data`
1051          @return: flux          :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
1052          @rtype: L{Data}                 for the permeability *k_{ij}*
         @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
1053          """          """
1054      self.setSubProblemTolerance()      self.setSubProblemTolerance()
1055          g=self.__g          g=self.__g
1056          f=self.__f          f=self.__f
1057          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
1058          if p == None:          if p == None:
1059             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
1060          else:          else:
# Line 380  class StokesProblemCartesian(Homogeneous Line 1080  class StokesProblemCartesian(Homogeneous
1080              sp.initialize(...)              sp.initialize(...)
1081              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
1082       """       """
1083       def __init__(self,domain,adaptSubTolerance=True, **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       @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.           with macro elements for the pressure.
1090       @type adaptSubTolerance: C{bool}  
1091           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           :param domain: domain of the problem.
1092             :type domain: `Domain`
1093           """           """
1094           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**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            
# Line 409  class StokesProblemCartesian(Homogeneous Line 1109  class StokesProblemCartesian(Homogeneous
1109           """           """
1110       returns the solver options used  solve the equation for velocity.       returns the solver options used  solve the equation for velocity.
1111            
1112       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
1113       """       """
1114       return self.__pde_u.getSolverOptions()       return self.__pde_u.getSolverOptions()
1115       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
1116           """           """
1117       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
1118            
1119       @param options: new solver  options       :param options: new solver  options
1120       @type options: L{SolverOptions}       :type options: `SolverOptions`
1121       """       """
1122           self.__pde_u.setSolverOptions(options)           self.__pde_u.setSolverOptions(options)
1123       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
1124           """           """
1125       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
1126       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
1127       """       """
1128       return self.__pde_prec.getSolverOptions()       return self.__pde_prec.getSolverOptions()
1129       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
1130           """           """
1131       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
1132       @param options: new solver  options       :param options: new solver  options
1133       @type options: L{SolverOptions}       :type options: `SolverOptions`
1134       """       """
1135       self.__pde_prec.setSolverOptions(options)       self.__pde_prec.setSolverOptions(options)
1136    
# Line 439  class StokesProblemCartesian(Homogeneous Line 1139  class StokesProblemCartesian(Homogeneous
1139       set the solver options for solving the equation to project the divergence of       set the solver options for solving the equation to project the divergence of
1140       the velocity onto the function space of presure.       the velocity onto the function space of presure.
1141            
1142       @param options: new solver options       :param options: new solver options
1143       @type options: L{SolverOptions}       :type options: `SolverOptions`
1144       """       """
1145       self.__pde_prec.setSolverOptions(options)       self.__pde_proj.setSolverOptions(options)
1146       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
1147           """           """
1148       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
1149       the velocity onto the function space of presure.       the velocity onto the function space of presure.
1150            
1151       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
1152       """       """
1153       return self.__pde_prec.getSolverOptions()       return self.__pde_proj.getSolverOptions()
1154       def setSubProblemTolerance(self):  
1155         def updateStokesEquation(self, v, p):
1156           """           """
1157       Updates the tolerance for subproblems           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       if self.adaptSubTolerance():           pass
1161               sub_tol=self.getSubProblemTolerance()       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
1162           self.getSolverOptionsDiv().setTolerance(sub_tol)          """
1163           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)          assigns new values to the model parameters.
          self.getSolverOptionsPressure().setTolerance(sub_tol)  
          self.getSolverOptionsPressure().setAbsoluteTolerance(0.)  
          self.getSolverOptionsVelocity().setTolerance(sub_tol)  
          self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)  
           
1164    
1165       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):          :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
1205          @note: All values needs to be set.          """
1206          """          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         self.eta=eta  
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
1207    
1208       def Bv(self,v):       def Bv(self,v,tol):
1209           """           """
1210           returns inner product of element p and div(v)           returns inner product of element p and div(v)
1211    
1212           @param p: a pressure increment           :param v: a residual
1213           @param v: a residual           :return: inner product of element p and div(v)
1214           @return: inner product of element p and div(v)           :rtype: ``float``
1215           @rtype: C{float}           """
1216           """           self.__pde_proj.setValue(Y=-util.div(v))
1217           self.__pde_proj.setValue(Y=-util.div(v))       self.getSolverOptionsDiv().setTolerance(tol)
1218           return self.__pde_proj.getSolution()       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
1219             out=self.__pde_proj.getSolution()
1220             return out
1221    
1222       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
1223           """           """
1224           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
1225    
1226           @param p: a pressure increment           :param p: a pressure increment
1227           @param v: a residual           :param Bv: a residual
1228           @return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
1229           @rtype: C{float}           :rtype: ``float``
1230           """           """
1231           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
1232    
1233       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
1234           """           """
1235           Returns inner product of p0 and p1           Returns inner product of p0 and p1
1236    
1237           @param p0: a pressure           :param p0: a pressure
1238           @param p1: a pressure           :param p1: a pressure
1239           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1240           @rtype: C{float}           :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 norm_v(self,v):       def norm_v(self,v):
1247           """           """
1248           returns the norm of v           returns the norm of v
1249    
1250           @param v: a velovity           :param v: a velovity
1251           @return: norm of v           :return: norm of v
1252           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1253           """           """
1254           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
1255    
1256       def getV(self, p, v0):  
1257         def getDV(self, p, v, tol):
1258           """           """
1259           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
1260    
1261           @param p: a pressure           :param p: a pressure
1262           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
1263           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
1264           """           """
1265           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.updateStokesEquation(v,p)
1266             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=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+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()           out=self.__pde_u.getSolution()
1274           return  out           return  out
1275    
# Line 560  class StokesProblemCartesian(Homogeneous Line 1277  class StokesProblemCartesian(Homogeneous
1277          """          """
1278          Returns Bv (overwrite).          Returns Bv (overwrite).
1279    
1280          @rtype: equal to the type of p          :rtype: equal to the type of p
1281          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1282          """          """
1283          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1284    
1285       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
1286           """           """
1287           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
1288    
1289           @param p: a pressure increment           :param p: a pressure increment
1290           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
1291           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
1292           """           """
1293           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1294           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
1295           return  out           return  out
1296    
1297       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
1298           """           """
1299           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
1300           with accuracy L{self.getSubProblemTolerance()}           with accuracy `self.getSubProblemTolerance()`
1301    
1302           @param v: velocity increment           :param Bv: velocity increment
1303           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1304           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
1305           """           """
1306           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
1307           return self.__pde_prec.getSolution()       self.getSolverOptionsPressure().setTolerance(tol)
1308         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1309             out=self.__pde_prec.getSolution()
1310             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26