/[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 1519 by artak, Tue Apr 22 03:45:36 2008 UTC revision 3500 by gross, Mon Apr 11 03:09:06 2011 UTC
# Line 1  Line 1 
1  # $Id:$  # -*- coding: utf-8 -*-
2    ########################################################
3  #  #
4  #######################################################  # Copyright (c) 2003-2010 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7  #  #
8  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
9  #  # Licensed under the Open Software License version 3.0
10  #                http://esscc.uq.edu.au  # http://www.opensource.org/licenses/osl-3.0.php
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
11  #  #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem  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):
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    
 class StokesProblemCartesian(HomogeneousSaddlePointProblem):  
406        """        """
407        solves        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
742    
743            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i      *u_i+k_{ij}*p_{,j} = g_i*
744        *u_{i,i} = f*
745    
746        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.
749        """
750    
751        def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
752            """
753            initializes the Darcy flux problem
754            :param domain: domain of the problem
755            :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
762            if weight == None:
763               s=self.domain.getSize()
764               self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
765               # self.__l2=(3.*util.longestEdge(self.domain))**2
766               #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
767            else:
768               self.__l2=weight
769            self.__pde_v=LinearPDESystem(domain)
770            if useReduced: self.__pde_v.setReducedOrderOn()
771            self.__pde_v.setSymmetryOn()
772            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)
774            self.__pde_p.setSymmetryOn()
775            if useReduced: self.__pde_p.setReducedOrderOn()
776            self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
777            self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
778            self.setTolerance()
779            self.setAbsoluteTolerance()
780        self.__adaptSubTolerance=adaptSubTolerance
781        self.verbose=False
782        def getSolverOptionsFlux(self):
783        """
784        Returns the solver options used to solve the flux problems
785        
786        *(I+D^*D)u=F*
787        
788        :return: `SolverOptions`
789        """
790        return self.__pde_v.getSolverOptions()
791        def setSolverOptionsFlux(self, options=None):
792        """
793        Sets the solver options used to solve the flux problems
794        
795        *(I+D^*D)u=F*
796        
797        If ``options`` is not present, the options are reset to default
798        :param options: `SolverOptions`
799        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
800        """
801        return self.__pde_v.setSolverOptions(options)
802        def getSolverOptionsPressure(self):
803        """
804        Returns the solver options used to solve the pressure problems
805        
806        *(Q^*Q)p=Q^*G*
807        
808        :return: `SolverOptions`
809        """
810        return self.__pde_p.getSolverOptions()
811        def setSolverOptionsPressure(self, options=None):
812        """
813        Sets the solver options used to solve the pressure problems
814        
815        *(Q^*Q)p=Q^*G*
816        
817        If ``options`` is not present, the options are reset to default
818        :param options: `SolverOptions`
819        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
820        """
821        return self.__pde_p.setSolverOptions(options)
822    
823        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
824            """
825            assigns values to model parameters
826    
827            :param f: volumetic sources/sinks
828            :type f: scalar value on the domain (e.g. `escript.Data`)
829            :param g: flux sources/sinks
830            :type g: vector values on the domain (e.g. `escript.Data`)
831            :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. `escript.Data`)
833            :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. `escript.Data`)
835            :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
836                                 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
837                                 ``v`` on the main diagonal is used.
838            :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 ``setValue`` are not altered.
841            :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 (``location_of_fixed_flux[i]>0`` if direction of the normal
843                   is along the *x_i* axis.
844            """
845            if f !=None:
846               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
847               if f.isEmpty():
848                   f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
849               else:
850                   if f.getRank()>0: raise ValueError,"illegal rank of f."
851               self.__f=f
852            if g !=None:
853               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
854               if g.isEmpty():
855                 g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
856               else:
857                 if not g.getShape()==(self.domain.getDim(),):
858                   raise ValueError,"illegal shape of g"
859               self.__g=g
860    
861            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
862            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
863    
864            if permeability!=None:
865               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
866               if perm.getRank()==0:
867                   perm=perm*util.kronecker(self.domain.getDim())
868               elif perm.getRank()==1:
869                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
870                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
871               elif perm.getRank()==2:
872                  pass
873               else:
874                  raise ValueError,"illegal rank of permeability."
875               self.__permeability=perm
876               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
877    
878        def setTolerance(self,rtol=1e-4):
879            """
880            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
881    
882            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
883    
884            where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
885    
886            :param rtol: relative tolerance for the pressure
887            :type rtol: non-negative ``float``
888            """
889            if rtol<0:
890                raise ValueError,"Relative tolerance needs to be non-negative."
891            self.__rtol=rtol
892        def getTolerance(self):
893            """
894            returns the relative tolerance
895    
896            :return: current relative tolerance
897            :rtype: ``float``
898            """
899            return self.__rtol
900    
901        def setAbsoluteTolerance(self,atol=0.):
902            """
903            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
904    
905            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
906    
907            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
910            :type atol: non-negative ``float``
911            """
912            if atol<0:
913                raise ValueError,"Absolute tolerance needs to be non-negative."
914            self.__atol=atol
915        def getAbsoluteTolerance(self):
916           """
917           returns the absolute tolerance
918          
919           :return: current absolute tolerance
920           :rtype: ``float``
921           """
922           return self.__atol
923        def getSubProblemTolerance(self):
924        """
925        Returns a suitable subtolerance
926        @type: ``float``
927        """
928        return max(util.EPSILON**(0.75),self.getTolerance()**2)
929        def setSubProblemTolerance(self):
930             """
931             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
932             """
933         if self.__adaptSubTolerance:
934             sub_tol=self.getSubProblemTolerance()
935                 self.getSolverOptionsFlux().setTolerance(sub_tol)
936             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
937             self.getSolverOptionsPressure().setTolerance(sub_tol)
938             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
939             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
940    
941        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
942             """
943             solves the problem.
944    
945             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 ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
948             :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 ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
950             :type p0: scalar value on the domain (e.g. `escript.Data`).
951             :param verbose: if set some information on iteration progress are printed
952             :type verbose: ``bool``
953             :return: flux and pressure
954             :rtype: ``tuple`` of `escript.Data`.
955    
956             :note: The problem is solved as a least squares form
957    
958             *(I+D^*D)u+Qp=D^*f+g*
959             *Q^*u+Q^*Qp=Q^*g*
960    
961             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
963    
964             *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
967    
968             *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 *Q^*Q*). In each iteration step
971             PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
972             """
973             self.verbose=verbose
974             rtol=self.getTolerance()
975             atol=self.getAbsoluteTolerance()
976         self.setSubProblemTolerance()
977             num_corrections=0
978             converged=False
979             p=p0
980             norm_r=None
981             while not converged:
982                   v=self.getFlux(p, fixed_flux=u0)
983                   Qp=self.__Q(p)
984                   norm_v=self.__L2(v)
985                   norm_Qp=self.__L2(Qp)
986                   if norm_v == 0.:
987                      if norm_Qp == 0.:
988                         return v,p
989                      else:
990                        fac=norm_Qp
991                   else:
992                      if norm_Qp == 0.:
993                        fac=norm_v
994                      else:
995                        fac=2./(1./norm_v+1./norm_Qp)
996                   ATOL=(atol+rtol*fac)
997                   if self.verbose:
998                        print "DarcyFlux: L2 norm of v = %e."%norm_v
999                        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
1003                   if norm_r == None or norm_r>ATOL:
1004                       if num_corrections>max_num_corrections:
1005                             raise ValueError,"maximum number of correction steps reached."
1006                       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
1008                   else:
1009                       converged=True
1010             return v,p
1011        def __L2(self,v):
1012             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1013    
1014        def __Q(self,p):
1015              return util.tensor_mult(self.__permeability,util.grad(p))
1016    
1017        def __Aprod(self,dp):
1018              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
1019              Qdp=self.__Q(dp)
1020              self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
1021              du=self.__pde_v.getSolution()
1022              return Qdp+du
1023        def __inner_GMRES(self,r,s):
1024             return util.integrate(util.inner(r,s))
1025    
1026        def __inner_PCG(self,p,r):
1027             return util.integrate(util.inner(self.__Q(p), r))
1028    
1029        def __Msolve_PCG(self,r):
1030          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
1031              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
1032              return self.__pde_p.getSolution()
1033    
1034        def getFlux(self,p=None, fixed_flux=escript.Data()):
1035            """
1036            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
1037            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
1038            Note that ``g`` and ``f`` are used, see `setValue`.
1039    
1040            :param p: pressure.
1041            :type p: scalar value on the domain (e.g. `escript.Data`).
1042            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
1043            :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
1044            :return: flux
1045            :rtype: `escript.Data`
1046            :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
1047                   for the permeability *k_{ij}*
1048            """
1049        self.setSubProblemTolerance()
1050            g=self.__g
1051            f=self.__f
1052            self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
1053            if p == None:
1054               self.__pde_v.setValue(Y=g)
1055            else:
1056               self.__pde_v.setValue(Y=g-self.__Q(p))
1057            return self.__pde_v.getSolution()
1058    
1059    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
1060         """
1061         solves
1062    
1063              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
1064                  u_{i,i}=0                  u_{i,i}=0
1065    
1066            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
1067            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
1068    
1069        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
1070    
1071        typical usage:       typical usage:
1072    
1073              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
1074              sp.setTolerance()              sp.setTolerance()
1075              sp.initialize(...)              sp.initialize(...)
1076              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
1077        """       """
1078        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
1079             """
1080             initialize the Stokes Problem
1081    
1082             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
1083             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
1084             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(preconditioner=LinearPDE.ILU0)      
               
1094           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
1095           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
1096           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
1097    
1098           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
1099           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
1100         self.__pde_proj.setValue(D=1)
1101           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
1102    
1103        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
1104          self.eta=eta           """
1105          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
1106      self.__pde_u.setValue(A=Data())      
1107          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
1108          for j in range(self.domain.getDim()):       """
1109              A[i,j,j,i] += 1.       return self.__pde_u.getSolverOptions()
1110              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
1111      self.__pde_prec.setValue(D=1./self.eta)           """
1112          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)       set the solver options for solving the equation for velocity.
1113        
1114        def B(self,arg):       :param options: new solver  options
1115           d=util.div(arg)       :type options: `SolverOptions`
1116           self.__pde_proj.setValue(Y=d)       """
1117           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setSolverOptions(options)
1118           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
1119             """
1120        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
1121           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
1122           s1=util.interpolate(p1,Function(self.domain))       """
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
1189    
1190            :param f: external force
1191            :type f: `Vector` object in `FunctionSpace` `Function` or similar
1192            :param fixed_u_mask: mask of locations with fixed velocity.
1193            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
1194            :param eta: viscosity
1195            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
1196            :param surface_stress: normal surface stress
1197            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
1198            :param stress: initial stress
1199        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
1200            """
1201            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
1202    
1203         def Bv(self,v,tol):
1204             """
1205             returns inner product of element p and div(v)
1206    
1207             :param v: a residual
1208             :return: inner product of element p and div(v)
1209             :rtype: ``float``
1210             """
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(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
1227    
1228         def inner_p(self,p0,p1):
1229             """
1230             Returns inner product of p0 and p1
1231    
1232             :param p0: a pressure
1233             :param p1: a pressure
1234             :return: inner product of p0 and p1
1235             :rtype: ``float``
1236             """
1237             s0=util.interpolate(p0, escript.Function(self.domain))
1238             s1=util.interpolate(p1, escript.Function(self.domain))
1239           return util.integrate(s0*s1)           return util.integrate(s0*s1)
1240    
1241        def getStress(self,u):       def norm_v(self,v):
1242           mg=util.grad(u)           """
1243           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
1244    
1245        def solve_A(self,u,p):           :param v: a velovity
1246           """           :return: norm of v
1247           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           :rtype: non-negative ``float``
1248           """           """
1249           self.__pde_u.setTolerance(self.getSubProblemTolerance())           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
          self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))  
          return  self.__pde_u.getSolution(verbose=self.show_details)  
   
       def solve_prec(self,p):  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.Lsup(v)  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)  
           self.iter+=1  
           if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
1250    
1251    
1252         def getDV(self, p, v, tol):
1253             """
1254             return the value for v for a given p (overwrite)
1255    
1256             :param p: a pressure
1257             :param v: a initial guess for the value v to return.
1258             :return: dv given as *Adv=(f-Av-B^*p)*
1259             """
1260             self.updateStokesEquation(v,p)
1261             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():
1265                self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1266             else:
1267                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()
1269             return  out
1270    
1271         def norm_Bv(self,Bv):
1272            """
1273            Returns Bv (overwrite).
1274    
1275            :rtype: equal to the type of p
1276            :note: boundary conditions on p should be zero!
1277            """
1278            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1279    
1280         def solve_AinvBt(self,p, tol):
1281             """
1282             Solves *Av=B^*p* with accuracy `tol`
1283    
1284             :param p: a pressure increment
1285             :return: the solution of *Av=B^*p*
1286             :note: boundary conditions on v should be zero!
1287             """
1288             self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1289             out=self.__pde_u.getSolution()
1290             return  out
1291    
1292         def solve_prec(self,Bv, tol):
1293             """
1294             applies preconditioner for for *BA^{-1}B^** to *Bv*
1295             with accuracy `self.getSubProblemTolerance()`
1296    
1297             :param Bv: velocity increment
1298             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1299             :note: boundary conditions on p are zero.
1300             """
1301             self.__pde_prec.setValue(Y=Bv)
1302         self.getSolverOptionsPressure().setTolerance(tol)
1303         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1304             out=self.__pde_prec.getSolution()
1305             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26