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

Legend:
Removed from v.2349  
changed lines
  Added in v.3500

  ViewVC Help
Powered by ViewVC 1.1.26