/[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 3360 by jfenwick, Thu Nov 18 00:20:21 2010 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 *  from escript import *
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    class DarcyFlow(object):
41       """
42       solves the problem
43      
44       *u_i+k_{ij}*p_{,j} = g_i*
45       *u_{i,i} = f*
46      
47       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48      
49       :note: The problem is solved in a least squares formulation.
50       """
51      
52       def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):
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          :type adaptSubTolerance: ``bool``
61          :param solveForFlux: if True the solver solves for the flux (do not use!)
62          :type solveForFlux: ``bool``  
63          """
64          self.domain=domain
65          self.solveForFlux=solveForFlux
66          self.useReduced=useReduced
67          self.__adaptSubTolerance=adaptSubTolerance
68          self.verbose=False
69          
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    
78          self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting
79          # self.__pde_l.setSymmetryOn()
80          # if self.useReduced: self.__pde_l.setReducedOrderOn()
81          self.setTolerance()
82          self.setAbsoluteTolerance()
83          self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
84          self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
85          
86       def getSolverOptionsFlux(self):
87          """
88          Returns the solver options used to solve the flux problems
89          
90          *K^{-1} u=F*
91          
92          :return: `SolverOptions`
93          """
94          return self.__pde_k.getSolverOptions()
95          
96       def setSolverOptionsFlux(self, options=None):
97          """
98          Sets the solver options used to solve the flux problems
99          
100          *K^{-1}u=F*
101          
102          If ``options`` is not present, the options are reset to default
103          
104          :param options: `SolverOptions`
105          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
106          """
107          return self.__pde_v.setSolverOptions(options)
108        
109       def getSolverOptionsPressure(self):
110          """
111          Returns the solver options used to solve the pressure problems
112          
113          *(Q^* K Q)p=-Q^*G*
114          
115          :return: `SolverOptions`
116          """
117          return self.__pde_p.getSolverOptions()
118          
119       def setSolverOptionsPressure(self, options=None):
120          """
121          Sets the solver options used to solve the pressure problems
122          
123          *(Q^* K Q)p=-Q^*G*
124          
125          If ``options`` is not present, the options are reset to default
126          
127          :param options: `SolverOptions`
128          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
129          """
130          return self.__pde_p.setSolverOptions(options)
131    
132  class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):     def getSolverOptionsWeighting(self):
133        """        """
134        solves        Returns the solver options used to solve the pressure problems
135    
136            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        *(D K D^*)p=-D F*
                 u_{i,i}=0  
137    
138            u=0 where  fixed_u_mask>0        :return: `SolverOptions`
139            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress        """
140          return self.__pde_l.getSolverOptions()
141    
142        if surface_stress is not give 0 is assumed.     def setSolverOptionsWeighting(self, options=None):
143          """
144          Sets the solver options used to solve the pressure problems
145    
146        typical usage:        *(D K D^*)p=-D F*
147    
148              sp=StokesProblemCartesian(domain)        If ``options`` is not present, the options are reset to default
149              sp.setTolerance()  
150              sp.initialize(...)        :param options: `SolverOptions`
151              v,p=sp.solve(v0,p0)        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
152        """        """
153        def __init__(self,domain,**kwargs):        return self.__pde_l.setSolverOptions(options)
          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  
154    
155    
156  class StokesProblemCartesian(HomogeneousSaddlePointProblem):     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
157          """
158          assigns values to model parameters
159    
160          :param f: volumetic sources/sinks
161          :type f: scalar value on the domain (e.g. `Data`)
162          :param g: flux sources/sinks
163          :type g: vector values on the domain (e.g. `Data`)
164          :param location_of_fixed_pressure: mask for locations where pressure is fixed
165          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
166          :param location_of_fixed_flux:  mask for locations where flux is fixed.
167          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
168          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
169          :type permeability: scalar or tensor values on the domain (e.g. `Data`)
170    
171          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
172          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
173          or the normal component of the flux (``location_of_fixed_flux[i]>0``) if direction of the normal is along the *x_i* axis.
174    
175          """
176          if f !=None:
177         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
178         if f.isEmpty():
179            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
180         else:
181            if f.getRank()>0: raise ValueError,"illegal rank of f."
182            self.__f=f
183          if g !=None:
184         g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
185         if g.isEmpty():
186            g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         else:
188            if not g.getShape()==(self.domain.getDim(),):
189               raise ValueError,"illegal shape of g"
190            self.__g=g
191          if location_of_fixed_pressure!=None:
192               self.__pde_p.setValue(q=location_of_fixed_pressure)
193               #self.__pde_l.setValue(q=location_of_fixed_pressure)
194          if location_of_fixed_flux!=None:
195               self.__pde_k.setValue(q=location_of_fixed_flux)
196                
197          if permeability!=None:
198         perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
199         if perm.getRank()==0:
200            perm_inv=(1./perm)*util.kronecker(self.domain.getDim())
201            perm=perm*util.kronecker(self.domain.getDim())
202         elif perm.getRank()==2:
203            perm_inv=util.inverse(perm)
204         else:
205            raise ValueError,"illegal rank of permeability."
206    
207         self.__permeability=perm
208         self.__permeability_inv=perm_inv
209         self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10
210         #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10
211         if  self.solveForFlux:
212            self.__pde_k.setValue(D=self.__permeability_inv)
213         else:
214            self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
215         self.__pde_p.setValue(A=self.__permeability)
216         #self.__pde_l.setValue(D=1/self.__l)
217             #self.__pde_l.setValue(A=self.__permeability)
218    
219       def __applWeight(self, v, f=None):
220          # solves L p = f-Dv with p = 0
221          if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"
222          if f == None:
223         return -util.div(v)*self.__l
224          else:
225         return (f-util.div(v))*self.__l
226          # if f == None:
227          #      self.__pde_l.setValue(Y=-util.div(v))  
228          # else:
229          #      return (f-util.div(v))/self.__l
230          # return self.__pde_l.getSolution()
231          
232       def __getPressure(self, v, p0, g=None):
233          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
234          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
235          if g == None:
236         self.__pde_p.setValue(X=-v, r=p0)
237          else:
238         self.__pde_p.setValue(X=g-v, r=p0)      
239          p=self.__pde_p.getSolution()
240          return p
241    
242       def __Aprod_v(self,dv):
243          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
244          dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)
245          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
246          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
247          return ArithmeticTuple(a,b,-dp)
248    
249       def __Msolve_PCG_v(self,r):
250          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
251          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
252          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())
253          # self.__pde_p.getOperator().saveMM("prec.mm")
254          return self.__pde_k.getSolution()
255    
256       def __inner_PCG_v(self,v,r):
257          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
258          
259       def __Aprod_p(self,dp):
260          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
261          Gdp=util.grad(dp)
262          self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())
263          du=self.__pde_k.getSolution()
264          # self.__pde_v.getOperator().saveMM("proj.mm")
265          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
266    
267       def __getFlux(self,p, v0, f=None, g=None):
268          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
269          if f!=None:
270         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
271          self.__pde_k.setValue(r=v0)
272          g2=util.tensor_mult(self.__permeability_inv,g)
273          if p == None:
274         self.__pde_k.setValue(Y=g2)
275          else:
276         self.__pde_k.setValue(Y=g2-util.grad(p))
277          return self.__pde_k.getSolution()  
278          
279          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
280       def __Msolve_PCG_p(self,r):
281          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
282          self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())
283          # self.__pde_p.getOperator().saveMM("prec.mm")
284          return self.__pde_p.getSolution()
285            
286       def __inner_PCG_p(self,p,r):
287           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
288    
289       def __L2(self,v):
290          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
291    
292       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
293          """
294          solves the problem.
295          
296          The iteration is terminated if the residual norm is less then self.getTolerance().
297    
298          :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.
299          :type u0: vector value on the domain (e.g. `Data`).
300          :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.
301          :type p0: scalar value on the domain (e.g. `Data`).
302          :param verbose: if set some information on iteration progress are printed
303          :type verbose: ``bool``
304          :return: flux and pressure
305          :rtype: ``tuple`` of `Data`.
306    
307          :note: The problem is solved as a least squares form
308          *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
309          *G^*u+*G^* K Gp=G^*g*
310          where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
311          """
312          self.verbose=verbose
313          rtol=self.getTolerance()
314          atol=self.getAbsoluteTolerance()
315          self.setSubProblemTolerance()
316          num_corrections=0
317          converged=False
318          norm_r=None
319          
320          # Eliminate the hydrostatic pressure:
321          if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
322          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))        
323          p0=self.__pde_p.getSolution()
324          g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
325          norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5    
326    
327          p=p0*0
328          if self.solveForFlux:
329         v=u0.copy()
330          else:
331         v=self.__getFlux(p, u0, f=self.__f, g=g2)
332    
333          while not converged and norm_g2 > 0:
334         Gp=util.grad(p)
335         KGp=util.tensor_mult(self.__permeability,Gp)
336         if self.verbose:
337            def_p=g2-(v+KGp)
338            def_v=self.__f-util.div(v)
339            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
340            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
341            print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
342            print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
343            print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
344         ATOL=atol+rtol*norm_g2
345         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
346         if norm_r == None or norm_r>ATOL:
347            if num_corrections>max_num_corrections:
348               raise ValueError,"maximum number of correction steps reached."
349          
350            if self.solveForFlux:
351               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
352               v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
353                   self.__Aprod_v,
354                   v,
355                   self.__Msolve_PCG_v,
356                   self.__inner_PCG_v,
357                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
358               p=r[2]
359            else:
360               # initial residual is r=G^*(g2-KGp - v)
361               p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
362                     self.__Aprod_p,
363                     p,
364                     self.__Msolve_PCG_p,
365                     self.__inner_PCG_p,
366                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
367               v=r[1]
368            if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
369            num_corrections+=1
370         else:
371            if self.verbose: print "DarcyFlux: stopping criterium reached."
372            converged=True
373          return v,p+p0
374       def setTolerance(self,rtol=1e-4):
375          """
376          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
377    
378          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
379          
380          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
381          
382          :param rtol: relative tolerance for the pressure
383          :type rtol: non-negative ``float``
384          """
385          if rtol<0:
386         raise ValueError,"Relative tolerance needs to be non-negative."
387          self.__rtol=rtol
388       def getTolerance(self):
389          """
390          returns the relative tolerance
391          :return: current relative tolerance
392          :rtype: ``float``
393        """        """
394        solves        return self.__rtol
395    
396            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i     def setAbsoluteTolerance(self,atol=0.):
397          """
398          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
399          
400          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
401    
402    
403          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}*.
404    
405          :param atol: absolute tolerance for the pressure
406          :type atol: non-negative ``float``
407          """
408          if atol<0:
409         raise ValueError,"Absolute tolerance needs to be non-negative."
410          self.__atol=atol
411       def getAbsoluteTolerance(self):
412          """
413          returns the absolute tolerance
414          :return: current absolute tolerance
415          :rtype: ``float``
416          """
417          return self.__atol
418       def getSubProblemTolerance(self):
419          """
420          Returns a suitable subtolerance
421          :type: ``float``
422          """
423          return max(util.EPSILON**(0.5),self.getTolerance()**2)
424    
425       def setSubProblemTolerance(self):
426          """
427          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
428          """
429          if self.__adaptSubTolerance:
430         sub_tol=self.getSubProblemTolerance()
431         self.getSolverOptionsFlux().setTolerance(sub_tol)
432         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
433         self.getSolverOptionsPressure().setTolerance(sub_tol)
434         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
435         self.getSolverOptionsWeighting().setTolerance(sub_tol)
436         self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)
437         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
438    
439    
440    class DarcyFlowOld(object):
441        """
442        solves the problem
443    
444        *u_i+k_{ij}*p_{,j} = g_i*
445        *u_{i,i} = f*
446    
447        where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
448    
449        :note: The problem is solved in a least squares formulation.
450        """
451    
452        def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
453            """
454            initializes the Darcy flux problem
455            :param domain: domain of the problem
456            :type domain: `Domain`
457        :param useReduced: uses reduced oreder on flux and pressure
458        :type useReduced: ``bool``
459        :param adaptSubTolerance: switches on automatic subtolerance selection
460        :type adaptSubTolerance: ``bool``  
461            """
462            self.domain=domain
463            if weight == None:
464               s=self.domain.getSize()
465               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
466               # self.__l=(3.*util.longestEdge(self.domain))**2
467               #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
468            else:
469               self.__l=weight
470            self.__pde_v=LinearPDESystem(domain)
471            if useReduced: self.__pde_v.setReducedOrderOn()
472            self.__pde_v.setSymmetryOn()
473            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
474            self.__pde_p=LinearSinglePDE(domain)
475            self.__pde_p.setSymmetryOn()
476            if useReduced: self.__pde_p.setReducedOrderOn()
477            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
478            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
479            self.setTolerance()
480            self.setAbsoluteTolerance()
481        self.__adaptSubTolerance=adaptSubTolerance
482        self.verbose=False
483        def getSolverOptionsFlux(self):
484        """
485        Returns the solver options used to solve the flux problems
486        
487        *(I+D^*D)u=F*
488        
489        :return: `SolverOptions`
490        """
491        return self.__pde_v.getSolverOptions()
492        def setSolverOptionsFlux(self, options=None):
493        """
494        Sets the solver options used to solve the flux problems
495        
496        *(I+D^*D)u=F*
497        
498        If ``options`` is not present, the options are reset to default
499        :param options: `SolverOptions`
500        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
501        """
502        return self.__pde_v.setSolverOptions(options)
503        def getSolverOptionsPressure(self):
504        """
505        Returns the solver options used to solve the pressure problems
506        
507        *(Q^*Q)p=Q^*G*
508        
509        :return: `SolverOptions`
510        """
511        return self.__pde_p.getSolverOptions()
512        def setSolverOptionsPressure(self, options=None):
513        """
514        Sets the solver options used to solve the pressure problems
515        
516        *(Q^*Q)p=Q^*G*
517        
518        If ``options`` is not present, the options are reset to default
519        :param options: `SolverOptions`
520        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
521        """
522        return self.__pde_p.setSolverOptions(options)
523    
524        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
525            """
526            assigns values to model parameters
527    
528            :param f: volumetic sources/sinks
529            :type f: scalar value on the domain (e.g. `Data`)
530            :param g: flux sources/sinks
531            :type g: vector values on the domain (e.g. `Data`)
532            :param location_of_fixed_pressure: mask for locations where pressure is fixed
533            :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
534            :param location_of_fixed_flux:  mask for locations where flux is fixed.
535            :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
536            :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
537                                 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
538                                 ``v`` on the main diagonal is used.
539            :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
540    
541            :note: the values of parameters which are not set by calling ``setValue`` are not altered.
542            :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
543                   or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
544                   is along the *x_i* axis.
545            """
546            if f !=None:
547               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
548               if f.isEmpty():
549                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
550               else:
551                   if f.getRank()>0: raise ValueError,"illegal rank of f."
552               self.__f=f
553            if g !=None:
554               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
555               if g.isEmpty():
556                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
557               else:
558                 if not g.getShape()==(self.domain.getDim(),):
559                   raise ValueError,"illegal shape of g"
560               self.__g=g
561    
562            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
563            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
564    
565            if permeability!=None:
566               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
567               if perm.getRank()==0:
568                   perm=perm*util.kronecker(self.domain.getDim())
569               elif perm.getRank()==1:
570                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
571                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
572               elif perm.getRank()==2:
573                  pass
574               else:
575                  raise ValueError,"illegal rank of permeability."
576               self.__permeability=perm
577               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
578    
579        def setTolerance(self,rtol=1e-4):
580            """
581            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
582    
583            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
584    
585            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}*.
586    
587            :param rtol: relative tolerance for the pressure
588            :type rtol: non-negative ``float``
589            """
590            if rtol<0:
591                raise ValueError,"Relative tolerance needs to be non-negative."
592            self.__rtol=rtol
593        def getTolerance(self):
594            """
595            returns the relative tolerance
596    
597            :return: current relative tolerance
598            :rtype: ``float``
599            """
600            return self.__rtol
601    
602        def setAbsoluteTolerance(self,atol=0.):
603            """
604            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
605    
606            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
607    
608            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}*.
609    
610            :param atol: absolute tolerance for the pressure
611            :type atol: non-negative ``float``
612            """
613            if atol<0:
614                raise ValueError,"Absolute tolerance needs to be non-negative."
615            self.__atol=atol
616        def getAbsoluteTolerance(self):
617           """
618           returns the absolute tolerance
619          
620           :return: current absolute tolerance
621           :rtype: ``float``
622           """
623           return self.__atol
624        def getSubProblemTolerance(self):
625        """
626        Returns a suitable subtolerance
627        @type: ``float``
628        """
629        return max(util.EPSILON**(0.75),self.getTolerance()**2)
630        def setSubProblemTolerance(self):
631             """
632             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
633             """
634         if self.__adaptSubTolerance:
635             sub_tol=self.getSubProblemTolerance()
636                 self.getSolverOptionsFlux().setTolerance(sub_tol)
637             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
638             self.getSolverOptionsPressure().setTolerance(sub_tol)
639             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
640             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
641    
642        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
643             """
644             solves the problem.
645    
646             The iteration is terminated if the residual norm is less then self.getTolerance().
647    
648             :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.
649             :type u0: vector value on the domain (e.g. `Data`).
650             :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.
651             :type p0: scalar value on the domain (e.g. `Data`).
652             :param verbose: if set some information on iteration progress are printed
653             :type verbose: ``bool``
654             :return: flux and pressure
655             :rtype: ``tuple`` of `Data`.
656    
657             :note: The problem is solved as a least squares form
658    
659             *(I+D^*D)u+Qp=D^*f+g*
660             *Q^*u+Q^*Qp=Q^*g*
661    
662             where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
663             We eliminate the flux form the problem by setting
664    
665             *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
666    
667             form the first equation. Inserted into the second equation we get
668    
669             *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
670    
671             which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
672             PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
673             """
674             self.verbose=verbose
675             rtol=self.getTolerance()
676             atol=self.getAbsoluteTolerance()
677         self.setSubProblemTolerance()
678             num_corrections=0
679             converged=False
680             p=p0
681             norm_r=None
682             while not converged:
683                   v=self.getFlux(p, fixed_flux=u0)
684                   Qp=self.__Q(p)
685                   norm_v=self.__L2(v)
686                   norm_Qp=self.__L2(Qp)
687                   if norm_v == 0.:
688                      if norm_Qp == 0.:
689                         return v,p
690                      else:
691                        fac=norm_Qp
692                   else:
693                      if norm_Qp == 0.:
694                        fac=norm_v
695                      else:
696                        fac=2./(1./norm_v+1./norm_Qp)
697                   ATOL=(atol+rtol*fac)
698                   if self.verbose:
699                        print "DarcyFlux: L2 norm of v = %e."%norm_v
700                        print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
701                        print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
702                        print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
703                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
704                   if norm_r == None or norm_r>ATOL:
705                       if num_corrections>max_num_corrections:
706                             raise ValueError,"maximum number of correction steps reached."
707                       p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
708                       num_corrections+=1
709                   else:
710                       converged=True
711             return v,p
712        def __L2(self,v):
713             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
714    
715        def __Q(self,p):
716              return util.tensor_mult(self.__permeability,util.grad(p))
717    
718        def __Aprod(self,dp):
719              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
720              Qdp=self.__Q(dp)
721              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
722              du=self.__pde_v.getSolution()
723              # self.__pde_v.getOperator().saveMM("proj.mm")
724              return Qdp+du
725        def __inner_GMRES(self,r,s):
726             return util.integrate(util.inner(r,s))
727    
728        def __inner_PCG(self,p,r):
729             return util.integrate(util.inner(self.__Q(p), r))
730    
731        def __Msolve_PCG(self,r):
732          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
733              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
734              # self.__pde_p.getOperator().saveMM("prec.mm")
735              return self.__pde_p.getSolution()
736    
737        def getFlux(self,p=None, fixed_flux=Data()):
738            """
739            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
740            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
741            Note that ``g`` and ``f`` are used, see `setValue`.
742    
743            :param p: pressure.
744            :type p: scalar value on the domain (e.g. `Data`).
745            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
746            :type fixed_flux: vector values on the domain (e.g. `Data`).
747            :return: flux
748            :rtype: `Data`
749            :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}*
750                   for the permeability *k_{ij}*
751            """
752        self.setSubProblemTolerance()
753            g=self.__g
754            f=self.__f
755            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
756            if p == None:
757               self.__pde_v.setValue(Y=g)
758            else:
759               self.__pde_v.setValue(Y=g-self.__Q(p))
760            return self.__pde_v.getSolution()
761    
762    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
763         """
764         solves
765    
766              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
767                  u_{i,i}=0                  u_{i,i}=0
768    
769            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
770            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
771    
772        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
773    
774        typical usage:       typical usage:
775    
776              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
777              sp.setTolerance()              sp.setTolerance()
778              sp.initialize(...)              sp.initialize(...)
779              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
780        """       """
781        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
782             """
783             initialize the Stokes Problem
784    
785             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
786             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
787             with macro elements for the pressure.
788    
789             :param domain: domain of the problem.
790             :type domain: `Domain`
791             """
792           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
793           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
794           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())
795           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
796           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)      
               
797           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
798           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
799           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
800    
801           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
802           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
803         self.__pde_proj.setValue(D=1)
804           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
805    
806        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
807          self.eta=eta           """
808          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
809      self.__pde_u.setValue(A=Data())      
810          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
811          for j in range(self.domain.getDim()):       """
812              A[i,j,j,i] += 1.       return self.__pde_u.getSolverOptions()
813              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
814      self.__pde_prec.setValue(D=1/self.eta)           """
815          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.
816        
817        def B(self,arg):       :param options: new solver  options
818           d=util.div(arg)       :type options: `SolverOptions`
819           self.__pde_proj.setValue(Y=d)       """
820           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setSolverOptions(options)
821           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
822             """
823         returns the solver options used  solve the equation for pressure.
824         :rtype: `SolverOptions`
825         """
826         return self.__pde_prec.getSolverOptions()
827         def setSolverOptionsPressure(self, options=None):
828             """
829         set the solver options for solving the equation for pressure.
830         :param options: new solver  options
831         :type options: `SolverOptions`
832         """
833         self.__pde_prec.setSolverOptions(options)
834    
835         def setSolverOptionsDiv(self, options=None):
836             """
837         set the solver options for solving the equation to project the divergence of
838         the velocity onto the function space of presure.
839        
840         :param options: new solver options
841         :type options: `SolverOptions`
842         """
843         self.__pde_proj.setSolverOptions(options)
844         def getSolverOptionsDiv(self):
845             """
846         returns the solver options for solving the equation to project the divergence of
847         the velocity onto the function space of presure.
848        
849         :rtype: `SolverOptions`
850         """
851         return self.__pde_proj.getSolverOptions()
852    
853        def inner(self,p0,p1):       def updateStokesEquation(self, v, p):
854             """
855             updates the Stokes equation to consider dependencies from ``v`` and ``p``
856             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
857             """
858             pass
859         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
860            """
861            assigns new values to the model parameters.
862    
863            :param f: external force
864            :type f: `Vector` object in `FunctionSpace` `Function` or similar
865            :param fixed_u_mask: mask of locations with fixed velocity.
866            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
867            :param eta: viscosity
868            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
869            :param surface_stress: normal surface stress
870            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
871            :param stress: initial stress
872        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
873            """
874            if eta !=None:
875                k=util.kronecker(self.domain.getDim())
876                kk=util.outer(k,k)
877                self.eta=util.interpolate(eta, Function(self.domain))
878            self.__pde_prec.setValue(D=1/self.eta)
879                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
880            if restoration_factor!=None:
881                n=self.domain.getNormal()
882                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
883            if fixed_u_mask!=None:
884                self.__pde_u.setValue(q=fixed_u_mask)
885            if f!=None: self.__f=f
886            if surface_stress!=None: self.__surface_stress=surface_stress
887            if stress!=None: self.__stress=stress
888    
889         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
890            """
891            assigns values to the model parameters
892    
893            :param f: external force
894            :type f: `Vector` object in `FunctionSpace` `Function` or similar
895            :param fixed_u_mask: mask of locations with fixed velocity.
896            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
897            :param eta: viscosity
898            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
899            :param surface_stress: normal surface stress
900            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
901            :param stress: initial stress
902        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
903            """
904            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
905    
906         def Bv(self,v,tol):
907             """
908             returns inner product of element p and div(v)
909    
910             :param v: a residual
911             :return: inner product of element p and div(v)
912             :rtype: ``float``
913             """
914             self.__pde_proj.setValue(Y=-util.div(v))
915         self.getSolverOptionsDiv().setTolerance(tol)
916         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
917             out=self.__pde_proj.getSolution()
918             return out
919    
920         def inner_pBv(self,p,Bv):
921             """
922             returns inner product of element p and Bv=-div(v)
923    
924             :param p: a pressure increment
925             :param Bv: a residual
926             :return: inner product of element p and Bv=-div(v)
927             :rtype: ``float``
928             """
929             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
930    
931         def inner_p(self,p0,p1):
932             """
933             Returns inner product of p0 and p1
934    
935             :param p0: a pressure
936             :param p1: a pressure
937             :return: inner product of p0 and p1
938             :rtype: ``float``
939             """
940           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
941           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
942           return util.integrate(s0*s1)           return util.integrate(s0*s1)
943    
944        def inner_a(self,a0,a1):       def norm_v(self,v):
945           p0=util.interpolate(a0[1],Function(self.domain))           """
946           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
          alfa=(1/self.vol)*util.integrate(p0)  
          beta=(1/self.vol)*util.integrate(p1)  
      v0=util.grad(a0[0])  
      v1=util.grad(a1[0])  
          return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))  
   
   
       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)-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  
947    
948             :param v: a velovity
949             :return: norm of v
950             :rtype: non-negative ``float``
951             """
952             return util.sqrt(util.integrate(util.length(util.grad(v))**2))
953    
954    
955         def getDV(self, p, v, tol):
956             """
957             return the value for v for a given p (overwrite)
958    
959             :param p: a pressure
960             :param v: a initial guess for the value v to return.
961             :return: dv given as *Adv=(f-Av-B^*p)*
962             """
963             self.updateStokesEquation(v,p)
964             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
965         self.getSolverOptionsVelocity().setTolerance(tol)
966         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
967             if self.__stress.isEmpty():
968                self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
969             else:
970                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
971             out=self.__pde_u.getSolution()
972             return  out
973    
974         def norm_Bv(self,Bv):
975            """
976            Returns Bv (overwrite).
977    
978            :rtype: equal to the type of p
979            :note: boundary conditions on p should be zero!
980            """
981            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
982    
983         def solve_AinvBt(self,p, tol):
984             """
985             Solves *Av=B^*p* with accuracy `tol`
986    
987             :param p: a pressure increment
988             :return: the solution of *Av=B^*p*
989             :note: boundary conditions on v should be zero!
990             """
991             self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
992             out=self.__pde_u.getSolution()
993             return  out
994    
995         def solve_prec(self,Bv, tol):
996             """
997             applies preconditioner for for *BA^{-1}B^** to *Bv*
998             with accuracy `self.getSubProblemTolerance()`
999    
1000             :param Bv: velocity increment
1001             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1002             :note: boundary conditions on p are zero.
1003             """
1004             self.__pde_prec.setValue(Y=Bv)
1005         self.getSolverOptionsPressure().setTolerance(tol)
1006         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1007             out=self.__pde_prec.getSolution()
1008             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26