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

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

  ViewVC Help
Powered by ViewVC 1.1.26