/[escript]/trunk/escript/py_src/flows.py
ViewVC logotype

Diff of /trunk/escript/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.1562  
changed lines
  Added in v.3462

  ViewVC Help
Powered by ViewVC 1.1.26