/[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 1673 by gross, Thu Jul 24 22:28:50 2008 UTC revision 3501 by gross, Wed Apr 13 04:07:53 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,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40    
41  class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):  
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        solves        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            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        :param f: volumetic sources/sinks
125                  u_{i,i}=0        :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            u=0 where  fixed_u_mask>0        """
142            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress        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        if surface_stress is not give 0 is assumed.        """
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        typical usage:        *|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):
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              sp=StokesProblemCartesian(domain)  
389              sp.setTolerance()     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
             sp.initialize(...)  
             v,p=sp.solve(v0,p0)  
390        """        """
391        def __init__(self,domain,**kwargs):        assigns values to model parameters
          HomogeneousSaddlePointProblem.__init__(self,**kwargs)  
          self.domain=domain  
          self.vol=util.integrate(1.,Function(self.domain))  
          self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())  
          self.__pde_u.setSymmetryOn()  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
          # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)  
          # self.__pde_proj.setReducedOrderOn()  
          # self.__pde_proj.setSymmetryOn()  
          # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)  
   
       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):  
         self.eta=eta  
         A =self.__pde_u.createCoefficientOfGeneralPDE("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.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)  
   
         # self.__pde_proj.setValue(D=1/eta)  
         # self.__pde_proj.setValue(Y=1.)  
         # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))  
         self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
   
       def B(self,arg):  
          a=util.div(arg, ReducedFunction(self.domain))  
          return a-util.integrate(a)/self.vol  
   
       def inner(self,p0,p1):  
          return util.integrate(p0*p1)  
          # return util.integrate(1/self.__inv_eta**2*p0*p1)  
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
       def getEtaEffective(self):  
          return self.eta  
   
       def solve_A(self,u,p):  
          """  
          solves Av=f-Au-B^*p (v=0 on fixed_u_mask)  
          """  
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
          self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))  
          return  self.__pde_u.getSolution(verbose=self.show_details)  
   
   
       def solve_prec(self,p):  
         a=self.__inv_eta*p  
         return a-util.integrate(a)/self.vol  
   
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= 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  
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    
 class StokesProblemCartesian(HomogeneousSaddlePointProblem):  
410        """        """
411        solves        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
747    
748        *u_i+k_{ij}*p_{,j} = g_i*
749        *u_{i,i} = f*
750    
751        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.
754        """
755    
756        def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
757            """
758            initializes the Darcy flux problem
759            :param domain: domain of the problem
760            :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
767            if weight == None:
768               s=self.domain.getSize()
769               self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
770               # self.__l2=(3.*util.longestEdge(self.domain))**2
771               #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
772            else:
773               self.__l2=weight
774            self.__pde_v=LinearPDESystem(domain)
775            if useReduced: self.__pde_v.setReducedOrderOn()
776            self.__pde_v.setSymmetryOn()
777            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)
779            self.__pde_p.setSymmetryOn()
780            if useReduced: self.__pde_p.setReducedOrderOn()
781            self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
782            self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
783            self.setTolerance()
784            self.setAbsoluteTolerance()
785        self.__adaptSubTolerance=adaptSubTolerance
786        self.verbose=False
787        def getSolverOptionsFlux(self):
788        """
789        Returns the solver options used to solve the flux problems
790        
791        *(I+D^*D)u=F*
792        
793        :return: `SolverOptions`
794        """
795        return self.__pde_v.getSolverOptions()
796        def setSolverOptionsFlux(self, options=None):
797        """
798        Sets the solver options used to solve the flux problems
799        
800        *(I+D^*D)u=F*
801        
802        If ``options`` is not present, the options are reset to default
803        :param options: `SolverOptions`
804        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
805        """
806        return self.__pde_v.setSolverOptions(options)
807        def getSolverOptionsPressure(self):
808        """
809        Returns the solver options used to solve the pressure problems
810        
811        *(Q^*Q)p=Q^*G*
812        
813        :return: `SolverOptions`
814        """
815        return self.__pde_p.getSolverOptions()
816        def setSolverOptionsPressure(self, options=None):
817        """
818        Sets the solver options used to solve the pressure problems
819        
820        *(Q^*Q)p=Q^*G*
821        
822        If ``options`` is not present, the options are reset to default
823        :param options: `SolverOptions`
824        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
825        """
826        return self.__pde_p.setSolverOptions(options)
827    
828        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
829            """
830            assigns values to model parameters
831    
832            :param f: volumetic sources/sinks
833            :type f: scalar value on the domain (e.g. `escript.Data`)
834            :param g: flux sources/sinks
835            :type g: vector values on the domain (e.g. `escript.Data`)
836            :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. `escript.Data`)
838            :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. `escript.Data`)
840            :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
841                                 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
842                                 ``v`` on the main diagonal is used.
843            :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 ``setValue`` are not altered.
846            :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 (``location_of_fixed_flux[i]>0`` if direction of the normal
848                   is along the *x_i* axis.
849            """
850            if f !=None:
851               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
852               if f.isEmpty():
853                   f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
854               else:
855                   if f.getRank()>0: raise ValueError,"illegal rank of f."
856               self.__f=f
857            if g !=None:
858               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
859               if g.isEmpty():
860                 g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
861               else:
862                 if not g.getShape()==(self.domain.getDim(),):
863                   raise ValueError,"illegal shape of g"
864               self.__g=g
865    
866            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
867            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
868    
869            if permeability!=None:
870               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
871               if perm.getRank()==0:
872                   perm=perm*util.kronecker(self.domain.getDim())
873               elif perm.getRank()==1:
874                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
875                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
876               elif perm.getRank()==2:
877                  pass
878               else:
879                  raise ValueError,"illegal rank of permeability."
880               self.__permeability=perm
881               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
882    
883        def setTolerance(self,rtol=1e-4):
884            """
885            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
886    
887            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
888    
889            where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
890    
891            :param rtol: relative tolerance for the pressure
892            :type rtol: non-negative ``float``
893            """
894            if rtol<0:
895                raise ValueError,"Relative tolerance needs to be non-negative."
896            self.__rtol=rtol
897        def getTolerance(self):
898            """
899            returns the relative tolerance
900    
901            :return: current relative tolerance
902            :rtype: ``float``
903            """
904            return self.__rtol
905    
906        def setAbsoluteTolerance(self,atol=0.):
907            """
908            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
909    
910            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
911    
912            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
915            :type atol: non-negative ``float``
916            """
917            if atol<0:
918                raise ValueError,"Absolute tolerance needs to be non-negative."
919            self.__atol=atol
920        def getAbsoluteTolerance(self):
921           """
922           returns the absolute tolerance
923          
924           :return: current absolute tolerance
925           :rtype: ``float``
926           """
927           return self.__atol
928        def getSubProblemTolerance(self):
929        """
930        Returns a suitable subtolerance
931        @type: ``float``
932        """
933        return max(util.EPSILON**(0.75),self.getTolerance()**2)
934        def setSubProblemTolerance(self):
935             """
936             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
937             """
938         if self.__adaptSubTolerance:
939             sub_tol=self.getSubProblemTolerance()
940                 self.getSolverOptionsFlux().setTolerance(sub_tol)
941             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
942             self.getSolverOptionsPressure().setTolerance(sub_tol)
943             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
944             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
945    
946        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
947             """
948             solves the problem.
949    
950             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 ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
953             :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 ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
955             :type p0: scalar value on the domain (e.g. `escript.Data`).
956             :param verbose: if set some information on iteration progress are printed
957             :type verbose: ``bool``
958             :return: flux and pressure
959             :rtype: ``tuple`` of `escript.Data`.
960    
961             :note: The problem is solved as a least squares form
962    
963             *(I+D^*D)u+Qp=D^*f+g*
964             *Q^*u+Q^*Qp=Q^*g*
965    
966             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
968    
969             *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
972    
973            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i           *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 *Q^*Q*). In each iteration step
976             PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
977             """
978             self.verbose=verbose
979             rtol=self.getTolerance()
980             atol=self.getAbsoluteTolerance()
981         self.setSubProblemTolerance()
982             num_corrections=0
983             converged=False
984             p=p0
985             norm_r=None
986             while not converged:
987                   v=self.getFlux(p, fixed_flux=u0)
988                   Qp=self.__Q(p)
989                   norm_v=self.__L2(v)
990                   norm_Qp=self.__L2(Qp)
991                   if norm_v == 0.:
992                      if norm_Qp == 0.:
993                         return v,p
994                      else:
995                        fac=norm_Qp
996                   else:
997                      if norm_Qp == 0.:
998                        fac=norm_v
999                      else:
1000                        fac=2./(1./norm_v+1./norm_Qp)
1001                   ATOL=(atol+rtol*fac)
1002                   if self.verbose:
1003                        print "DarcyFlux: L2 norm of v = %e."%norm_v
1004                        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
1008                   if norm_r == None or norm_r>ATOL:
1009                       if num_corrections>max_num_corrections:
1010                             raise ValueError,"maximum number of correction steps reached."
1011                       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
1013                   else:
1014                       converged=True
1015             return v,p
1016        def __L2(self,v):
1017             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
1018    
1019        def __Q(self,p):
1020              return util.tensor_mult(self.__permeability,util.grad(p))
1021    
1022        def __Aprod(self,dp):
1023              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
1024              Qdp=self.__Q(dp)
1025              self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
1026              du=self.__pde_v.getSolution()
1027              return Qdp+du
1028        def __inner_GMRES(self,r,s):
1029             return util.integrate(util.inner(r,s))
1030    
1031        def __inner_PCG(self,p,r):
1032             return util.integrate(util.inner(self.__Q(p), r))
1033    
1034        def __Msolve_PCG(self,r):
1035          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
1036              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
1037              return self.__pde_p.getSolution()
1038    
1039        def getFlux(self,p=None, fixed_flux=escript.Data()):
1040            """
1041            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
1042            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
1043            Note that ``g`` and ``f`` are used, see `setValue`.
1044    
1045            :param p: pressure.
1046            :type p: scalar value on the domain (e.g. `escript.Data`).
1047            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
1048            :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
1049            :return: flux
1050            :rtype: `escript.Data`
1051            :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
1052                   for the permeability *k_{ij}*
1053            """
1054        self.setSubProblemTolerance()
1055            g=self.__g
1056            f=self.__f
1057            self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
1058            if p == None:
1059               self.__pde_v.setValue(Y=g)
1060            else:
1061               self.__pde_v.setValue(Y=g-self.__Q(p))
1062            return self.__pde_v.getSolution()
1063    
1064    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
1065         """
1066         solves
1067    
1068              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
1069                  u_{i,i}=0                  u_{i,i}=0
1070    
1071            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
1072            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
1073    
1074        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
1075    
1076        typical usage:       typical usage:
1077    
1078              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
1079              sp.setTolerance()              sp.setTolerance()
1080              sp.initialize(...)              sp.initialize(...)
1081              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
1082        """       """
1083        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
1084             """
1085             initialize the Stokes Problem
1086    
1087             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
1088             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
1089             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(preconditioner=LinearPDE.ILU0)      
               
1099           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
1100           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
1101           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
1102    
1103           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
1104           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
1105         self.__pde_proj.setValue(D=1)
1106           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
1107    
1108        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
1109          self.eta=eta           """
1110          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
1111      self.__pde_u.setValue(A=Data())      
1112          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
1113          for j in range(self.domain.getDim()):       """
1114              A[i,j,j,i] += 1.       return self.__pde_u.getSolverOptions()
1115              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
1116      self.__pde_prec.setValue(D=1/self.eta)           """
1117          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.
1118        
1119        def B(self,arg):       :param options: new solver  options
1120           d=util.div(arg)       :type options: `SolverOptions`
1121           self.__pde_proj.setValue(Y=d)       """
1122           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setSolverOptions(options)
1123           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
1124             """
1125        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
1126           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
1127           s1=util.interpolate(p1,Function(self.domain))       """
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
1194    
1195            :param f: external force
1196            :type f: `Vector` object in `FunctionSpace` `Function` or similar
1197            :param fixed_u_mask: mask of locations with fixed velocity.
1198            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
1199            :param eta: viscosity
1200            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
1201            :param surface_stress: normal surface stress
1202            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
1203            :param stress: initial stress
1204        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
1205            """
1206            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
1207    
1208         def Bv(self,v,tol):
1209             """
1210             returns inner product of element p and div(v)
1211    
1212             :param v: a residual
1213             :return: inner product of element p and div(v)
1214             :rtype: ``float``
1215             """
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(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
1232    
1233         def inner_p(self,p0,p1):
1234             """
1235             Returns inner product of p0 and p1
1236    
1237             :param p0: a pressure
1238             :param p1: a pressure
1239             :return: inner product of p0 and p1
1240             :rtype: ``float``
1241             """
1242             s0=util.interpolate(p0, escript.Function(self.domain))
1243             s1=util.interpolate(p1, escript.Function(self.domain))
1244           return util.integrate(s0*s1)           return util.integrate(s0*s1)
1245    
1246        def inner_a(self,a0,a1):       def norm_v(self,v):
1247           p0=util.interpolate(a0[1],Function(self.domain))           """
1248           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
1249           alfa=(1/self.vol)*util.integrate(p0)  
1250           beta=(1/self.vol)*util.integrate(p1)           :param v: a velovity
1251       v0=util.grad(a0[0])           :return: norm of v
1252       v1=util.grad(a1[0])           :rtype: non-negative ``float``
1253           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))           """
1254             return util.sqrt(util.integrate(util.length(util.grad(v))**2))
1255    
1256        def getStress(self,u):  
1257           mg=util.grad(u)       def getDV(self, p, v, tol):
1258           return 2.*self.eta*util.symmetric(mg)           """
1259        def getEtaEffective(self):           return the value for v for a given p (overwrite)
          return self.eta  
   
       def solve_A(self,u,p):  
          """  
          solves Av=f-Au-B^*p (v=0 on fixed_u_mask)  
          """  
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
          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):  
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
   
       def solve_prec1(self,p):  
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
      q0=util.interpolate(q,Function(self.domain))  
          print util.inf(q*q0),util.sup(q*q0)  
          q-=(1/self.vol)*util.integrate(q0)  
          print util.inf(q*q0),util.sup(q*q0)  
          return q  
   
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= 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  
1260    
1261             :param p: a pressure
1262             :param v: a initial guess for the value v to return.
1263             :return: dv given as *Adv=(f-Av-B^*p)*
1264             """
1265             self.updateStokesEquation(v,p)
1266             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
1267         self.getSolverOptionsVelocity().setTolerance(tol)
1268         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
1269             if self.__stress.isEmpty():
1270                self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1271             else:
1272                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()
1274             return  out
1275    
1276         def norm_Bv(self,Bv):
1277            """
1278            Returns Bv (overwrite).
1279    
1280            :rtype: equal to the type of p
1281            :note: boundary conditions on p should be zero!
1282            """
1283            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1284    
1285         def solve_AinvBt(self,p, tol):
1286             """
1287             Solves *Av=B^*p* with accuracy `tol`
1288    
1289             :param p: a pressure increment
1290             :return: the solution of *Av=B^*p*
1291             :note: boundary conditions on v should be zero!
1292             """
1293             self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1294             out=self.__pde_u.getSolution()
1295             return  out
1296    
1297         def solve_prec(self,Bv, tol):
1298             """
1299             applies preconditioner for for *BA^{-1}B^** to *Bv*
1300             with accuracy `self.getSubProblemTolerance()`
1301    
1302             :param Bv: velocity increment
1303             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1304             :note: boundary conditions on p are zero.
1305             """
1306             self.__pde_prec.setValue(Y=Bv)
1307         self.getSolverOptionsPressure().setTolerance(tol)
1308         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1309             out=self.__pde_prec.getSolution()
1310             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26