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

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

  ViewVC Help
Powered by ViewVC 1.1.26