/[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 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC revision 3452 by caltinay, Tue Jan 25 01:53:57 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      Represents and solves the problem     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 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      M{u_i+k_{ij}*p_{,j} = g_i}     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      M{u_{i,i} = f}     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    
     where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents  
     the permeability.  
499    
500      @note: The problem is solved in a least squares formulation.  class DarcyFlowOld(object):
501      """      """
502        solves the problem
503    
504      def __init__(self, domain):      *u_i+k_{ij}*p_{,j} = g_i*
505          """      *u_{i,i} = f*
506          Initializes the Darcy flux problem.  
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          @param domain: domain of the problem      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
513          @type domain: L{Domain}          """
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          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          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          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=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          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          if useReduced: self.__pde_p.setReducedOrderOn()
537          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          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):      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: volumetric 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 value 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 value 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          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
597                               tensor with C{s} on the main diagonal is used. If                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
598                               vector C{v} is given the tensor with C{v} on the                               ``v`` on the main diagonal is used.
599                               main diagonal is used.          :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
600          @type permeability: scalar, vector or tensor values on the domain, e.g.  
601                              L{Data}          :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          @note: the values of parameters which are not set by calling                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
604                 C{setValue} are not altered                 is along the *x_i* axis.
         @note: at any point on the boundary of the domain the pressure  
                (C{location_of_fixed_pressure}) >0 or the normal component of  
                the flux (C{location_of_fixed_flux[i]}) >0 if the direction of  
                the normal is along the M{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 125  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):
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 getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):      def setAbsoluteTolerance(self,atol=0.):
663          """          """
664          Returns the flux for a given pressure C{p}.          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             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          The flux is equal to C{fixed_flux} on locations where           where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
723          C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}           We eliminate the flux form the problem by setting
         and C{f} are used.  
   
         @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 by  
                            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 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}}  
         """  
         self.__pde_v.setTolerance(tol)  
         self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)  
         return self.__pde_v.getSolution(verbose=show_details)  
   
     def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):  
         """  
         Solves the problem.  
   
         The iteration is terminated if the error in the pressure is less than  
         M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand  
         side (see escript user's guide for details).  
   
         @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.  
         @type u0: vector value on the domain, e.g. L{Data}  
         @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.  
         @type p0: scalar value on the domain, e.g. L{Data}  
         @param atol: absolute tolerance for the pressure  
         @type atol: non-negative C{float}  
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{float}  
         @param sub_rtol: tolerance to be used in the sub iteration. It is  
                          recommended that M{sub_rtol<rtol*5.e-3}  
         @type sub_rtol: positive-negative C{float}  
         @param verbose: if True information on iteration progress is printed  
         @type verbose: C{bool}  
         @param show_details: if True information on the sub-iteration process  
                              is printed  
         @type show_details: C{bool}  
         @return: flux and pressure  
         @rtype: C{tuple} of L{Data}  
   
         @note: the problem is solved in a least squares formulation:  
   
         M{(I+D^*D)u+Qp=D^*f+g}  
   
         M{Q^*u+Q^*Qp=Q^*g}  
   
         where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the  
         permeability M{k_{ij}}. We eliminate the flux from the problem by  
         setting  
   
         M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}  
   
         from the first equation. Inserted into the second equation we get  
   
         M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}  
         on C{location_of_fixed_pressure}  
   
         which is solved using the PCG method (precondition is M{Q^*Q}).  
         In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}  
         need to be solved using a sub-iteration scheme.  
         """  
         self.verbose=verbose  
         self.show_details= show_details and self.verbose  
         self.__pde_v.setTolerance(sub_rtol)  
         self.__pde_p.setTolerance(sub_rtol)  
         u2=u0*self.__pde_v.getCoefficient("q")  
         #  
         # first the reference velocity is calculated from  
         #  
         #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)  
         #  
         self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)  
         u_ref=self.__pde_v.getSolution(verbose=show_details)  
         if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)  
         self.__pde_v.setValue(r=Data())  
         #  
         #   and then we calculate a reference pressure  
         #  
         #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)  
         #  
         self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)  
         p_ref=self.__pde_p.getSolution(verbose=self.show_details)  
         if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)  
         self.__pde_p.setValue(r=Data())  
         #  
         #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref  
         #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref  
         #  
         #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du  
         #  
         #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref  
         #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref  
         #  
         #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)  
         #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)  
         #  
         #  we use the (du,Qdp) to represent the resudual  
         #  Q^*Q is a preconditioner  
         #  
         #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>  
         #  
         Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))  
         norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))  
         ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)  
         if not ATOL>0:  
             raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)  
         if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)  
         #  
         #   caclulate the initial residual  
         #  
         self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())  
         du=self.__pde_v.getSolution(verbose=show_details)  
         r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)  
         dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
         util.saveVTK("d.vtu",p=dp,p_ref=p_ref)  
         return u_ref+r[1],dp  
   
     def __Aprod_PCG(self,p):  
         if self.show_details: print "DarcyFlux: Applying operator"  
         Qp=util.tensor_mult(self.__permeability,util.grad(p))  
         self.__pde_v.setValue(Y=Qp,X=Data())  
         w=self.__pde_v.getSolution(verbose=self.show_details)  
         return ArithmeticTuple(-Qp,w)  
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):      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))
         out=-util.integrate(util.inner(a,r[0]+r[1]))  
         return out  
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]))            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  class StokesProblemCartesian(HomogeneousSaddlePointProblem):      def getFlux(self,p=None, fixed_flux=escript.Data()):
796        """          """
797        Represents and solves the problem          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        M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
821         """
822         solves
823    
824        M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
825                    u_{i,i}=0
826    
827        M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}            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
829    
830        If C{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           Initializes the Stokes Problem.           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. The approximation order needs           :param domain: domain of the problem.
848                          to be two.           :type domain: `Domain`
          @type domain: L{Domain}  
          @warning: The approximation order needs to be two otherwise you may  
                    see oscillations in the pressure.  
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           Assigns values to the model parameters.       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           @param f: external force       def updateStokesEquation(self, v, p):
912           @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar           """
913           @param fixed_u_mask: mask of locations with fixed velocity           updates the Stokes equation to consider dependencies from ``v`` and ``p``
914           @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
915                               or similar           """
916           @param eta: viscosity           pass
917           @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
918           @param surface_stress: normal surface stress          """
919           @type surface_stress: L{Vector} object on L{FunctionSpace},          assigns new values to the model parameters.
920                                 L{FunctionOnBoundary} or similar  
921           @param stress: initial stress          :param f: external force
922           @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or          :type f: `Vector` object in `FunctionSpace` `Function` or similar
923                         similar          :param fixed_u_mask: mask of locations with fixed velocity.
924           @note: All values need to be set.          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
925           """          :param eta: viscosity
926           self.eta=eta          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
927           A =self.__pde_u.createCoefficient("A")          :param surface_stress: normal surface stress
928           self.__pde_u.setValue(A=Data())          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
929           for i in range(self.domain.getDim()):          :param stress: initial stress
930               for j in range(self.domain.getDim()):      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
931                   A[i,j,j,i] += 1.          """
932                   A[i,j,i,j] += 1.          if eta !=None:
933           self.__pde_prec.setValue(D=1/self.eta)              k=util.kronecker(self.domain.getDim())
934           self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)              kk=util.outer(k,k)
935           self.__stress=stress              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 B(self,v):       def Bv(self,v,tol):
965           """           """
966           Returns M{div(v)}.           returns inner product of element p and div(v)
          @return: M{div(v)}  
          @rtype: equal to the type of p  
967    
968           @note: Boundary conditions on p should be zero!           :param v: a residual
969             :return: inner product of element p and div(v)
970             :rtype: ``float``
971           """           """
972           if self.show_details: print "apply divergence:"           self.__pde_proj.setValue(Y=-util.div(v))
973           self.__pde_proj.setValue(Y=-util.div(v))       self.getSolverOptionsDiv().setTolerance(tol)
974           self.__pde_proj.setTolerance(self.getSubProblemTolerance())       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
975           return self.__pde_proj.getSolution(verbose=self.show_details)           out=self.__pde_proj.getSolution()
976             return out
977    
978        def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
979           """           """
980           Returns inner product of element p and Bv (overwrite).           returns inner product of element p and Bv=-div(v)
981    
982           @type p: equal to the type of p           :param p: a pressure increment
983           @type Bv: equal to the type of result of operator B           :param Bv: a residual
984           @return: inner product of p and Bv           :return: inner product of element p and Bv=-div(v)
985           @rtype: equal to the type of p           :rtype: ``float``
986           """           """
987           s0=util.interpolate(p,Function(self.domain))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
          s1=util.interpolate(Bv,Function(self.domain))  
          return util.integrate(s0*s1)  
988    
989        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
990           """           """
991           Returns inner product of element p0 and p1 (overwrite).           Returns inner product of p0 and p1
992    
993           @type p0: equal to the type of p           :param p0: a pressure
994           @type p1: equal to the type of p           :param p1: a pressure
995           @return: inner product of p0 and p1           :return: inner product of p0 and p1
996           @rtype: equal to the type of p           :rtype: ``float``
997           """           """
998           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
999           s1=util.interpolate(p1/self.eta,Function(self.domain))           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 elements v0 and v1 (overwrite).           returns the norm of v
1005    
1006           @type v0: equal to the type of v           :param v: a velovity
1007           @type v1: equal to the type of v           :return: norm of v
1008           @return: inner product of v0 and v1           :rtype: non-negative ``float``
          @rtype: equal to the type of v  
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 M{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 norm_Bv(self,Bv):
1033            """
1034            Returns Bv (overwrite).
1035    
1036        def solve_prec(self,p):          :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           Applies the preconditioner.           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           if self.show_details: print "apply preconditioner:"           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1050           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           out=self.__pde_u.getSolution()
1051           self.__pde_prec.setValue(Y=p)           return  out
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
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.2169  
changed lines
  Added in v.3452

  ViewVC Help
Powered by ViewVC 1.1.26