/[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 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 3440 by gross, Fri Jan 14 00:04:53 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __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"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40    print dir(escript)
41    
42  class DarcyFlow(object):  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 (``location_of_fixed_pressure`` >0)
163          or the normal component of the flux (``location_of_fixed_flux[i]>0``) if direction of the normal 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 then 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          :param verbose: if set some information on iteration progress are printed
298          :type verbose: ``bool``
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      solves the problem
503    
504      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
505      M{u_{i,i} = f}      *u_{i,i} = f*
506    
507      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      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.      :note: The problem is solved in a least squares formulation.
510      """      """
511    
512      def __init__(self, domain,useReduced=False):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
513          """          """
514          initializes the Darcy flux problem          initializes the Darcy flux problem
515          @param domain: domain of the problem          :param domain: domain of the problem
516          @type domain: L{Domain}          :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          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)          self.__pde_v=LinearPDESystem(domain)
531          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
532          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
533          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          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)          self.__pde_p=LinearSinglePDE(domain)
535          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
536          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
537          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
538          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
539          self.__ATOL= None          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):      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          assigns values to model parameters
587    
588          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
589          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `escript.Data`)
590          @param g: flux sources/sinks          :param g: flux sources/sinks
591          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `escript.Data`)
592          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :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. L{Data})          :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.          :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. L{Data})          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
596          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
597                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
598                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
599          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          :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 C{setValue} are not altered.          :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 (C{location_of_fixed_pressure} >0)          :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 (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
604                 is along the M{x_i} axis.                 is along the *x_i* axis.
605          """          """
606          if f !=None:          if f !=None:
607             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
608             if f.isEmpty():             if f.isEmpty():
609                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
610             else:             else:
611                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
612             self.f=f             self.__f=f
613          if g !=None:            if g !=None:
614             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
615             if g.isEmpty():             if g.isEmpty():
616               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
617             else:             else:
618               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
619                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 121  class DarcyFlow(object): Line 636  class DarcyFlow(object):
636             self.__permeability=perm             self.__permeability=perm
637             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
638    
639        def setTolerance(self,rtol=1e-4):
     def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):  
         """  
         returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}  
         on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
         Note that C{g} and C{f} are used, see L{setValue}.  
           
         @param p: pressure.  
         @type p: scalar value on the domain (e.g. L{Data}).  
         @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.  
         @type fixed_flux: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: flux  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
640          """          """
641          self.__pde_v.setTolerance(tol)          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
642          g=self.__g  
643          f=self.__f          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
         self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
         else:  
            self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))  
         return self.__pde_v.getSolution(verbose=show_details)  
644    
645      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          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          returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
651          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
652          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
653                def getTolerance(self):
         @param v: flux.  
         @type v: vector-valued on the domain (e.g. L{Data}).  
         @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.  
         @type fixed_pressure: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: pressure  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
654          """          """
655          self.__pde_v.setTolerance(tol)          returns the relative tolerance
         g=self.__g  
         self.__pde_p.setValue(r=fixed_pressure)  
         if v == None:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))  
         else:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))  
         return self.__pde_p.getSolution(verbose=show_details)  
656    
657      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          :return: current relative tolerance
658            :rtype: ``float``
659          """          """
660          set the tolerance C{ATOL} used to terminate the solution process. It is used          return self.__rtol
   
         M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}  
   
         where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.  
661    
662          The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual      def setAbsoluteTolerance(self,atol=0.):
663            """
664          M{r=Q^*(g-Qp-v)}          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
   
         the condition  
665    
666          M{<(Q^*Q)^{-1} r,r> <= ATOL}          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
667    
668          holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}          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          :param atol: absolute tolerance for the pressure
671          @type atol: non-negative C{float}          :type atol: non-negative ``float``
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{float}  
         @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the  
         L{getPressure} method or use  the value from a previous time step.  
         @type p_ref: scalar value on the domain (e.g. L{Data}).  
         @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the  
         L{getFlux} method or use  the value from a previous time step.  
         @type v_ref: vector-valued on the domain (e.g. L{Data}).  
         @return: used absolute tolerance.  
         @rtype: positive C{float}  
672          """          """
673          g=self.__g          if atol<0:
674          if not v_ref == None:              raise ValueError,"Absolute tolerance needs to be non-negative."
675             f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)          self.__atol=atol
676          else:      def getAbsoluteTolerance(self):
677             f1=util.integrate(util.length(util.interpolate(g))**2)         """
678          if not p_ref == None:         returns the absolute tolerance
679             f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)        
680          else:         :return: current absolute tolerance
681             f2=0         :rtype: ``float``
682          self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))         """
683          if self.__ATOL<=0:         return self.__atol
684             raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL      def getSubProblemTolerance(self):
685          return self.__ATOL      """
686                Returns a suitable subtolerance
687      def getTolerance(self):      @type: ``float``
688          """      """
689          returns the current tolerance.      return max(util.EPSILON**(0.75),self.getTolerance()**2)
690          def setSubProblemTolerance(self):
691          @return: used absolute tolerance.           """
692          @rtype: positive C{float}           Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
693          """           """
694          if self.__ATOL==None:       if self.__adaptSubTolerance:
695             raise ValueError,"no tolerance is defined."           sub_tol=self.getSubProblemTolerance()
696          return self.__ATOL               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, show_details=False, sub_rtol=1.e-8):      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
703           """           """
704           solves the problem.           solves the problem.
705    
706           The iteration is terminated if the residual norm is less then self.getTolerance().           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 C{location_of_fixed_flux} the value of C{u0} is kept unchanged.           :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. L{Data}).           :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 C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.           :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. L{Data}).           :type p0: scalar value on the domain (e.g. `escript.Data`).
712           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}           :param verbose: if set some information on iteration progress are printed
713           @type sub_rtol: positive-negative C{float}           :type verbose: ``bool``
714           @param verbose: if set some information on iteration progress are printed           :return: flux and pressure
715           @type verbose: C{bool}           :rtype: ``tuple`` of `escript.Data`.
716           @param show_details:  if set information on the subiteration process are printed.  
717           @type show_details: C{bool}           :note: The problem is solved as a least squares form
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
   
          @note: The problem is solved as a least squares form  
718    
719           M{(I+D^*D)u+Qp=D^*f+g}           *(I+D^*D)u+Qp=D^*f+g*
720           M{Q^*u+Q^*Qp=Q^*g}           *Q^*u+Q^*Qp=Q^*g*
721    
722           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           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           We eliminate the flux form the problem by setting
724    
725           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           *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           form the first equation. Inserted into the second equation we get
728    
729           M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure           *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 M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
732           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
733           """           """
734           self.verbose=verbose           self.verbose=verbose
735           self.show_details= show_details and self.verbose           rtol=self.getTolerance()
736           self.__pde_v.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
737           self.__pde_p.setTolerance(sub_rtol)       self.setSubProblemTolerance()
738           ATOL=self.getTolerance()           num_corrections=0
739           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL           converged=False
740           #########################################################################################################################           p=p0
741           #           norm_r=None
742           #   we solve:           while not converged:
743           #                   v=self.getFlux(p, fixed_flux=u0)
744           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )                 Qp=self.__Q(p)
745           #                 norm_v=self.__L2(v)
746           #   residual is                 norm_Qp=self.__L2(Qp)
747           #                 if norm_v == 0.:
748           #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)                    if norm_Qp == 0.:
749           #                       return v,p
750           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                    else:
751           #                      fac=norm_Qp
752           #    we use (g - Qp, v) to represent the residual. not that                 else:
753           #                    if norm_Qp == 0.:
754           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                      fac=norm_v
755           #                    else:
756           #   while the initial residual is                      fac=2./(1./norm_v+1./norm_Qp)
757           #                 ATOL=(atol+rtol*fac)
758           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                 if self.verbose:
759           #                        print "DarcyFlux: L2 norm of v = %e."%norm_v
760           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
761           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
762           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
763           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
764           self.__pde_v.setValue(r=Data())                 if norm_r == None or norm_r>ATOL:
765           # start CG                     if num_corrections>max_num_corrections:
766           r=ArithmeticTuple(d0, v00)                           raise ValueError,"maximum number of correction steps reached."
767           p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                     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           return r[1],p                     num_corrections+=1
769                   else:
770      def __Aprod_PCG(self,dp):                     converged=True
771            if self.show_details: print "DarcyFlux: Applying operator"           return v,p
772            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)      def __L2(self,v):
773            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
774            self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())  
775            du=self.__pde_v.getSolution(verbose=self.show_details)      def __Q(self,p):
776            return ArithmeticTuple(mQdp,du)            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):      def __inner_PCG(self,p,r):
788           a=util.tensor_mult(self.__permeability,util.grad(p))           return util.integrate(util.inner(self.__Q(p), r))
          f0=util.integrate(util.inner(a,r[0]))  
          f1=util.integrate(util.inner(a,r[1]))  
          # print "__inner_PCG:",f0,f1,"->",f0-f1  
          return f0-f1  
789    
790      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
791            if self.show_details: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
792            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data())
793            return self.__pde_p.getSolution(verbose=self.show_details)            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):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
821        """       """
822        solves       solves
823    
824            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
825                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 827  class StokesProblemCartesian(Homogeneous
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-p*n_i=surface_stress +stress_{ij}n_j            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 given 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           initialize the Stokes Problem
842    
843           @param domain: domain of the problem. The approximation order needs to be two.           The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
844           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
845           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           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(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
855           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
856           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
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(),stress=Data()):       def getSolverOptionsVelocity(self):
865             """
866         returns the solver options used  solve the equation for velocity.
867        
868         :rtype: `SolverOptions`
869         """
870         return self.__pde_u.getSolverOptions()
871         def setSolverOptionsVelocity(self, options=None):
872             """
873         set the solver options for solving the equation for velocity.
874        
875         :param options: new solver  options
876         :type options: `SolverOptions`
877         """
878             self.__pde_u.setSolverOptions(options)
879         def getSolverOptionsPressure(self):
880             """
881         returns the solver options used  solve the equation for pressure.
882         :rtype: `SolverOptions`
883         """
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          assigns values to the model parameters
950    
951          @param f: external force          :param f: external force
952          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
953          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
954          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
955          @param eta: viscosity          :param eta: viscosity
956          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
957          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
958          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
959          @param stress: initial stress          :param stress: initial stress
960      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
         @note: All values needs to be set.  
   
961          """          """
962          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)  
         self.__stress=stress  
   
       def B(self,v):  
         """  
         returns div(v)  
         @rtype: equal to the type of p  
   
         @note: boundary conditions on p should be zero!  
         """  
         if self.show_details: print "apply divergence:"  
         self.__pde_proj.setValue(Y=-util.div(v))  
         self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
         return self.__pde_proj.getSolution(verbose=self.show_details)  
963    
964        def inner_pBv(self,p,Bv):       def Bv(self,v,tol):
965           """           """
966           returns inner product of element p and Bv  (overwrite)           returns inner product of element p and div(v)
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @rtype: C{float}  
967    
968           @rtype: equal to the type of p           :param v: a residual
969             :return: inner product of element p and div(v)
970             :rtype: ``float``
971           """           """
972           s0=util.interpolate(p,Function(self.domain))           self.__pde_proj.setValue(Y=-util.div(v))
973           s1=util.interpolate(Bv,Function(self.domain))       self.getSolverOptionsDiv().setTolerance(tol)
974           return util.integrate(s0*s1)       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
975             out=self.__pde_proj.getSolution()
976             return out
977    
978        def inner_p(self,p0,p1):       def inner_pBv(self,p,Bv):
979           """           """
980           returns inner product of element p0 and p1  (overwrite)           returns inner product of element p and Bv=-div(v)
981            
982           @type p0: equal to the type of p           :param p: a pressure increment
983           @type p1: equal to the type of p           :param Bv: a residual
984           @rtype: C{float}           :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           @rtype: equal to the type of p       def inner_p(self,p0,p1):
990           """           """
991           s0=util.interpolate(p0/self.eta,Function(self.domain))           Returns inner product of p0 and p1
992           s1=util.interpolate(p1/self.eta,Function(self.domain))  
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_v(self,v0,v1):       def norm_v(self,v):
1003           """           """
1004           returns inner product of two element v0 and v1  (overwrite)           returns the norm of v
           
          @type v0: equal to the type of v  
          @type v1: equal to the type of v  
          @rtype: C{float}  
1005    
1006           @rtype: equal to the type of v           :param v: a velovity
1007             :return: norm of v
1008             :rtype: non-negative ``float``
1009           """           """
1010       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
1011    
1012        def solve_A(self,u,p):  
1013         def getDV(self, p, v, tol):
1014           """           """
1015           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           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           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
1022           self.__pde_u.setTolerance(self.getSubProblemTolerance())           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():           if self.__stress.isEmpty():
1026              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
1027           else:           else:
1028              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              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(verbose=self.show_details)           out=self.__pde_u.getSolution()
1030           return  out           return  out
1031    
1032        def solve_prec(self,p):       def norm_Bv(self,Bv):
1033           if self.show_details: print "apply preconditioner:"          """
1034           self.__pde_prec.setTolerance(self.getSubProblemTolerance())          Returns Bv (overwrite).
1035           self.__pde_prec.setValue(Y=p)  
1036           q=self.__pde_prec.getSolution(verbose=self.show_details)          :rtype: equal to the type of p
1037           return q          :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.2208  
changed lines
  Added in v.3440

  ViewVC Help
Powered by ViewVC 1.1.26