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

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

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

revision 2156 by gross, Mon Dec 15 05:09:02 2008 UTC revision 3515 by gross, Thu May 19 08:20: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  class DarcyFlow(object):  class DarcyFlow(object):
41      """     """
42      solves the problem     solves the problem
43      
44      M{u_i+k_{ij}*p_{,j} = g_i}     *u_i+k_{ij}*p_{,j} = g_i*
45      M{u_{i,i} = f}     *u_{i,i} = f*
46      
47      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,
48      
49      @note: The problem is solved in a least squares formulation.     :cvar SIMPLE: simple solver
50      """     :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52      def __init__(self, domain):     :cvar SYMSTAB: solver uses symmetric stabilization
53          """     """
54          initializes the Darcy flux problem     SIMPLE="SIMPLE"
55          @param domain: domain of the problem     POST="POST"
56          @type domain: L{Domain}     STAB="STAB"
57          """     SYMSTAB="SYMSTAB"
58          self.domain=domain     def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59          self.__pde_v=LinearPDESystem(domain)        """
60          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))        initializes the Darcy flux problem
61          self.__pde_v.setSymmetryOn()        :param domain: domain of the problem
62          self.__pde_p=LinearSinglePDE(domain)        :type domain: `Domain`
63          self.__pde_p.setSymmetryOn()        :param useReduced: uses reduced oreder on flux and pressure
64          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        :type useReduced: ``bool``
65          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        :param solver: solver method
66          :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):        :param verbose: if ``True`` some information on the iteration progress are printed.
68          """        :type verbose: ``bool``
69          assigns values to model parameters        :param w: weighting factor for `DarcyFlow.POST` solver
70          :type w: ``float``
71          
72          """
73          self.domain=domain
74          self.solver=solver
75          self.useReduced=useReduced
76          self.verbose=verbose
77          self.scale=1.
78          
79          
80          self.__pde_v=LinearPDESystem(domain)
81          self.__pde_v.setSymmetryOn()
82          if self.useReduced: self.__pde_v.setReducedOrderOn()
83          self.__pde_p=LinearSinglePDE(domain)
84          self.__pde_p.setSymmetryOn()
85          if self.useReduced: self.__pde_p.setReducedOrderOn()
86          
87          if self.solver  == self.SIMPLE:
88         if self.verbose: print "DarcyFlow: simple solver is used."
89             self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90          elif self.solver  == self.POST:
91         self.w=w
92         if util.inf(w)<0.:
93            raise ValueError,"Weighting factor must be non-negative."
94         if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95          elif self.solver  == self.STAB:
96          if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97          elif  self.solver  == self.SYMSTAB:
98          if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99          else:
100        raise ValueError,"unknown solver %s"%self.solver
101          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104          self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105          self.setTolerance()
106        
107            
108       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109          """
110          assigns values to model parameters
111    
112          @param f: volumetic sources/sinks        :param f: volumetic sources/sinks
113          @type f: scalar value on the domain (e.g. L{Data})        :type f: scalar value on the domain (e.g. `escript.Data`)
114          @param g: flux sources/sinks        :param g: flux sources/sinks
115          @type g: vector values on the domain (e.g. L{Data})        :type g: vector values on the domain (e.g. `escript.Data`)
116          @param location_of_fixed_pressure: mask for locations where pressure is fixed        :param location_of_fixed_pressure: mask for locations where pressure is fixed
117          @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`)
118          @param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
119          @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`)
120          @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 ``s`` on the main diagonal is used.
121                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with        :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122                               C{v} on the main diagonal is used.  
123          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})        :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124          :note: at any point on the boundary of the domain the pressure
125          @note: the values of parameters which are not set by calling C{setValue} are not altered.               (``location_of_fixed_pressure`` >0) or the normal component of the
126          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)               flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal               is along the *x_i* axis.
                is along the M{x_i} axis.  
         """  
         if f !=None:  
            f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            if f.isEmpty():  
                f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            else:  
                if f.getRank()>0: raise ValueError,"illegal rank of f."  
            self.f=f  
         if g !=None:    
            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
            if g.isEmpty():  
              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
            else:  
              if not g.getShape()==(self.domain.getDim(),):  
                raise ValueError,"illegal shape of g"  
            self.__g=g  
   
         if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
         if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)  
   
         if permeability!=None:  
            perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
            if perm.getRank()==0:  
                perm=perm*util.kronecker(self.domain.getDim())  
            elif perm.getRank()==1:  
                perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm  
                for i in range(self.domain.getDim()): perm[i,i]=perm2[i]  
            elif perm.getRank()==2:  
               pass  
            else:  
               raise ValueError,"illegal rank of permeability."  
            self.__permeability=perm  
            self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))  
128    
129          """
130          if location_of_fixed_pressure!=None:
131               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
134              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136          
137                
138          # pressure  is rescaled by the factor 1/self.scale
139          if permeability!=None:
140        
141         perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142             V=util.vol(self.domain)
143             l=V**(1./self.domain.getDim())
144            
145         if perm.getRank()==0:
146            perm_inv=(1./perm)
147                self.scale=util.integrate(perm_inv)/V*l
148            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150            
151            
152         elif perm.getRank()==2:
153            perm_inv=util.inverse(perm)
154                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155            perm_inv*=(1./self.scale)
156            perm=perm*self.scale
157         else:
158            raise ValueError,"illegal rank of permeability."
159            
160         self.__permeability=perm
161         self.__permeability_inv=perm_inv
162         if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163        
164         if self.solver  == self.SIMPLE:
165            self.__pde_p.setValue(A=self.__permeability)
166         elif self.solver  == self.POST:
167            self.__pde_p.setValue(A=self.__permeability)
168            k=util.kronecker(self.domain.getDim())
169            self.lamb = self.w*util.length(perm_inv)*l
170            self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171         elif self.solver  == self.STAB:
172            self.__pde_p.setValue(A=0.5*self.__permeability)
173            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174         elif  self.solver  == self.SYMSTAB:
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177    
178          if g != None:
179        g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180        if g.isEmpty():
181              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182        else:
183            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184        self.__g=g
185          if f !=None:
186         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         if f.isEmpty():      
188              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189         else:
190             if f.getRank()>0: raise ValueError,"illegal rank of f."
191         self.__f=f
192       def getSolverOptionsFlux(self):
193          """
194          Returns the solver options used to solve the flux problems
195          :return: `SolverOptions`
196          """
197          return self.__pde_v.getSolverOptions()
198          
199       def setSolverOptionsFlux(self, options=None):
200          """
201          Sets the solver options used to solve the flux problems
202          If ``options`` is not present, the options are reset to default
203          :param options: `SolverOptions`
204          """
205          return self.__pde_v.setSolverOptions(options)
206        
207       def getSolverOptionsPressure(self):
208          """
209          Returns the solver options used to solve the pressure problems
210          :return: `SolverOptions`
211          """
212          return self.__pde_p.getSolverOptions()
213          
214       def setSolverOptionsPressure(self, options=None):
215          """
216          Sets the solver options used to solve the pressure problems
217          If ``options`` is not present, the options are reset to default
218          
219          :param options: `SolverOptions`
220          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221          """
222          return self.__pde_p.setSolverOptions(options)
223          
224       def setTolerance(self,rtol=1e-4):
225          """
226          sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227          
228          :param rtol: relative tolerance for the pressure
229          :type rtol: non-negative ``float``
230          """
231          if rtol<0:
232         raise ValueError,"Relative tolerance needs to be non-negative."
233          self.__rtol=rtol
234          
235       def getTolerance(self):
236          """
237          returns the relative tolerance
238          :return: current relative tolerance
239          :rtype: ``float``
240          """
241          return self.__rtol
242          
243       def solve(self,u0,p0, max_iter=100, iter_restart=20):
244          """
245          solves the problem.
246          
247          The iteration is terminated if the residual norm is less then self.getTolerance().
248    
249          :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.
250          :type u0: vector value on the domain (e.g. `escript.Data`).
251          :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.
252          :type p0: scalar value on the domain (e.g. `escript.Data`).
253          :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254          :type max_iter: ``int``
255          :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256                               A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257                               of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258          :type iter_restart: ``int``
259          :return: flux and pressure
260          :rtype: ``tuple`` of `escript.Data`.
261    
262      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):        """
263          # rescale initial guess:
264          p0=p0/self.scale
265          if self.solver  == self.SIMPLE or self.solver  == self.POST :
266            self.__pde_p.setValue(X=self.__g ,
267                                  Y=self.__f,
268                                  y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269                                  r=p0)
270            p=self.__pde_p.getSolution()
271            u = self.getFlux(p, u0)
272          elif  self.solver  == self.STAB:
273        u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274          elif  self.solver  == self.SYMSTAB:
275        u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276        
277          if self.verbose:
278            KGp=util.tensor_mult(self.__permeability,util.grad(p))
279            def_p=self.__g-(u+KGp)
280            def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281            print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282            print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283          #rescale result
284          p=p*self.scale
285          return u,p
286          
287       def getFlux(self,p, u0=None):
288          """          """
289          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291          Note that C{g} and C{f} are used, L{setValue}.          Notice that ``g`` and ``f`` are used, see `setValue`.
292            
293          @param p: pressure.          :param p: pressure.
294          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `escript.Data`).
295          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296          @type fixed_flux: vector values on the domain (e.g. L{Data}).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297          @param tol: relative tolerance to be used.          :return: flux
298          @type tol: positive float.          :rtype: `escript.Data`
         @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}}  
299          """          """
300          self.__pde_v.setTolerance(tol)          if self.solver  == self.SIMPLE or self.solver  == self.POST  :
301          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)              KGp=util.tensor_mult(self.__permeability,util.grad(p))
302          return self.__pde_v.getSolution(verbose=show_details)              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303                if u0 == None:
304      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):             self.__pde_v.setValue(r=escript.Data())
305           """          else:
306           solves the problem.             self.__pde_v.setValue(r=u0)
307                u= self.__pde_v.getSolution()
308           The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where      elif self.solver  == self.POST:
309           C{|q|} denotes the norm of the right hand side (see escript user's guide for details).              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310                                      X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311           @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.              if u0 == None:
312           @type u0: vector value on the domain (e.g. L{Data}).             self.__pde_v.setValue(r=escript.Data())
313           @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.          else:
314           @type p0: scalar value on the domain (e.g. L{Data}).             self.__pde_v.setValue(r=u0)
315           @param atol: absolute tolerance for the pressure              u= self.__pde_v.getSolution()
316           @type atol: non-negative C{float}      elif self.solver  == self.STAB:
317           @param rtol: relative tolerance for the pressure           gp=util.grad(p)
318           @type rtol: non-negative C{float}           self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}                                 X= p * util.kronecker(self.domain.getDim()),
320           @type sub_rtol: positive-negative C{float}                                 y= - p * self.domain.getNormal())                          
321           @param verbose: if set some information on iteration progress are printed           if u0 == None:
322           @type verbose: C{bool}             self.__pde_v.setValue(r=escript.Data())
323           @param show_details:  if set information on the subiteration process are printed.           else:
324           @type show_details: C{bool}             self.__pde_v.setValue(r=u0)
325           @return: flux and pressure           u= self.__pde_v.getSolution()
326           @rtype: C{tuple} of L{Data}.      elif  self.solver  == self.SYMSTAB:
327             gp=util.grad(p)
328           @note: The problem is solved as a least squares form           self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329                                   X= escript.Data() ,
330           M{(I+D^*D)u+Qp=D^*f+g}                                 y= escript.Data() )                          
331           M{Q^*u+Q^*Qp=Q^*g}           if u0 == None:
332               self.__pde_v.setValue(r=escript.Data())
333           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           else:
334           We eliminate the flux form the problem by setting             self.__pde_v.setValue(r=u0)
335             u= self.__pde_v.getSolution()
336           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux      return u
337          
338        
339       def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340              # p0 is used as an initial guess
341          u=self.getFlux(p0, u0)  
342              self.__pde_p.setValue( Y=self.__f-util.div(u),
343                                     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344                                     y= escript.Data(),
345                                     r=escript.Data())
346    
347          dp=self.__pde_p.getSolution()
348          p=GMRES(dp,
349                  self.__STAB_Aprod,
350              p0,
351              self.__inner,
352              atol=self.__norm(p0+dp)*self.getTolerance() ,
353              rtol=0.,
354              iter_max=max_iter,
355              iter_restart=iter_restart,
356              verbose=self.verbose,P_R=None)
357                
358              u=self.getFlux(p, u0)
359              return u,p
360    
361           form the first equation. Inserted into the second equation we get     def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362              # p0 is used as an initial guess
363          u=self.getFlux(p0, u0)  
364              self.__pde_p.setValue( Y= self.__f,
365                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366                                     y=  -  util.inner(self.domain.getNormal(), u),
367                                     r=escript.Data())
368          dp=self.__pde_p.getSolution()
369          
370          p=GMRES(dp,
371                  self.__SYMSTAB_Aprod,
372              p0,
373              self.__inner,
374              atol=self.__norm(p0+dp)*self.getTolerance() ,
375              rtol=0.,
376              iter_max=max_iter,
377              iter_restart=iter_restart,
378              verbose=self.verbose,P_R=None)
379                
380              u=self.getFlux(p, u0)
381              return u,p
382    
383           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     def __L2(self,v):
384             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
385      
386       def __norm(self,r):
387             return util.sqrt(self.__inner(r,r))
388                    
389           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step     def __inner(self,r,s):
390           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           return util.integrate(util.inner(r,s), escript.Function(self.domain))
391           """          
392           self.verbose=verbose     def __STAB_Aprod(self,p):
393           self.show_details= show_details and self.verbose        gp=util.grad(p)
394           self.__pde_v.setTolerance(sub_rtol)        self.__pde_v.setValue(Y=-0.5*gp,
395           self.__pde_p.setTolerance(sub_rtol)                              X=-p*util.kronecker(self.__pde_v.getDomain()),
396           u2=u0*self.__pde_v.getCoefficient("q")                              y= p * self.domain.getNormal(),  
397           #                              r=escript.Data())
398           # first the reference velocity is calculated from        u = -self.__pde_v.getSolution()
399           #        self.__pde_p.setValue(Y=util.div(u),
400           #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)                              X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
401           #                              y=escript.Data(),
402           self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)                              r=escript.Data())
403           u_ref=self.__pde_v.getSolution(verbose=show_details)      
404           if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)        return  self.__pde_p.getSolution()
405           self.__pde_v.setValue(r=Data())    
406           #     def __SYMSTAB_Aprod(self,p):
407           #   and then we calculate a reference pressure        gp=util.grad(p)
408           #        self.__pde_v.setValue(Y=0.5*gp ,
409           #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)                              X=escript.Data(),
410           #                              y=escript.Data(),  
411           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)                              r=escript.Data())
412           p_ref=self.__pde_p.getSolution(verbose=self.show_details)        u = -self.__pde_v.getSolution()
413           if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)        self.__pde_p.setValue(Y=escript.Data(),
414           self.__pde_p.setValue(r=Data())                              X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
415           #                              y=escript.Data(),
416           #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref                              r=escript.Data())
417           #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref      
418           #        return  self.__pde_p.getSolution()
419           #      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 =%, 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)  
   
     def __inner_PCG(self,p,r):  
          a=util.tensor_mult(self.__permeability,util.grad(p))  
          out=-util.integrate(util.inner(a,r[0]+r[1]))  
          return out  
   
     def __Msolve_PCG(self,r):  
           if self.show_details: print "DarcyFlux: Applying preconditioner"  
           self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))  
           return self.__pde_p.getSolution(verbose=self.show_details)  
420    
421  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
422        """       """
423        solves       solves
424    
425            -(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}
426                  u_{i,i}=0                  u_{i,i}=0
# Line 264  class StokesProblemCartesian(Homogeneous Line 428  class StokesProblemCartesian(Homogeneous
428            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
429            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
430    
431        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
432    
433        typical usage:       typical usage:
434    
435              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
436              sp.setTolerance()              sp.setTolerance()
437              sp.initialize(...)              sp.initialize(...)
438              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
439        """       """
440        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
441           """           """
442           initialize the Stokes Problem           initialize the Stokes Problem
443    
444           @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
445           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
446           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
447    
448             :param domain: domain of the problem.
449             :type domain: `Domain`
450           """           """
451           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
452           self.domain=domain           self.domain=domain
453           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
454           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
455           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
456           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
457           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
458           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
459    
460           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
461           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
462         self.__pde_proj.setValue(D=1)
463           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
464    
465        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
466             """
467         returns the solver options used  solve the equation for velocity.
468        
469         :rtype: `SolverOptions`
470         """
471         return self.__pde_v.getSolverOptions()
472         def setSolverOptionsVelocity(self, options=None):
473             """
474         set the solver options for solving the equation for velocity.
475        
476         :param options: new solver  options
477         :type options: `SolverOptions`
478         """
479             self.__pde_v.setSolverOptions(options)
480         def getSolverOptionsPressure(self):
481             """
482         returns the solver options used  solve the equation for pressure.
483         :rtype: `SolverOptions`
484         """
485         return self.__pde_prec.getSolverOptions()
486         def setSolverOptionsPressure(self, options=None):
487             """
488         set the solver options for solving the equation for pressure.
489         :param options: new solver  options
490         :type options: `SolverOptions`
491         """
492         self.__pde_prec.setSolverOptions(options)
493    
494         def setSolverOptionsDiv(self, options=None):
495             """
496         set the solver options for solving the equation to project the divergence of
497         the velocity onto the function space of presure.
498        
499         :param options: new solver options
500         :type options: `SolverOptions`
501         """
502         self.__pde_proj.setSolverOptions(options)
503         def getSolverOptionsDiv(self):
504             """
505         returns the solver options for solving the equation to project the divergence of
506         the velocity onto the function space of presure.
507        
508         :rtype: `SolverOptions`
509         """
510         return self.__pde_proj.getSolverOptions()
511    
512         def updateStokesEquation(self, v, p):
513             """
514             updates the Stokes equation to consider dependencies from ``v`` and ``p``
515             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
516             """
517             pass
518         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
519          """          """
520          assigns values to the model parameters          assigns new values to the model parameters.
521    
522          @param f: external force          :param f: external force
523          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
524          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
525          @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
526          @param eta: viscosity          :param eta: viscosity
527          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
528          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
529          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
530          @param stress: initial stress          :param stress: initial stress
531      @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.  
   
532          """          """
533          self.eta=eta          if eta !=None:
534          A =self.__pde_u.createCoefficient("A")              k=util.kronecker(self.domain.getDim())
535      self.__pde_u.setValue(A=Data())              kk=util.outer(k,k)
536          for i in range(self.domain.getDim()):              self.eta=util.interpolate(eta, escript.Function(self.domain))
537          for j in range(self.domain.getDim()):          self.__pde_prec.setValue(D=1/self.eta)
538              A[i,j,j,i] += 1.              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
539              A[i,j,i,j] += 1.          if restoration_factor!=None:
540      self.__pde_prec.setValue(D=1/self.eta)              n=self.domain.getNormal()
541          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
542          self.__stress=stress          if fixed_u_mask!=None:
543                self.__pde_v.setValue(q=fixed_u_mask)
544            if f!=None: self.__f=f
545            if surface_stress!=None: self.__surface_stress=surface_stress
546            if stress!=None: self.__stress=stress
547    
548        def B(self,v):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
549          """          """
550          returns div(v)          assigns values to the model parameters
         @rtype: equal to the type of p  
551    
552          @note: boundary conditions on p should be zero!          :param f: external force
553            :type f: `Vector` object in `FunctionSpace` `Function` or similar
554            :param fixed_u_mask: mask of locations with fixed velocity.
555            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
556            :param eta: viscosity
557            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
558            :param surface_stress: normal surface stress
559            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
560            :param stress: initial stress
561        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
562          """          """
563          if self.show_details: print "apply divergence:"          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         self.__pde_proj.setValue(Y=-util.div(v))  
         self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
         return self.__pde_proj.getSolution(verbose=self.show_details)  
564    
565        def inner_pBv(self,p,Bv):       def Bv(self,v,tol):
566           """           """
567           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}  
568    
569           @rtype: equal to the type of p           :param v: a residual
570             :return: inner product of element p and div(v)
571             :rtype: ``float``
572           """           """
573           s0=util.interpolate(p,Function(self.domain))           self.__pde_proj.setValue(Y=-util.div(v))
574           s1=util.interpolate(Bv,Function(self.domain))       self.getSolverOptionsDiv().setTolerance(tol)
575           return util.integrate(s0*s1)       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
576             out=self.__pde_proj.getSolution()
577             return out
578    
579        def inner_p(self,p0,p1):       def inner_pBv(self,p,Bv):
580           """           """
581           returns inner product of element p0 and p1  (overwrite)           returns inner product of element p and Bv=-div(v)
582            
583           @type p0: equal to the type of p           :param p: a pressure increment
584           @type p1: equal to the type of p           :param Bv: a residual
585           @rtype: C{float}           :return: inner product of element p and Bv=-div(v)
586             :rtype: ``float``
587             """
588             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
589    
590         def inner_p(self,p0,p1):
591             """
592             Returns inner product of p0 and p1
593    
594           @rtype: equal to the type of p           :param p0: a pressure
595             :param p1: a pressure
596             :return: inner product of p0 and p1
597             :rtype: ``float``
598           """           """
599           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
600           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
601           return util.integrate(s0*s1)           return util.integrate(s0*s1)
602    
603        def inner_v(self,v0,v1):       def norm_v(self,v):
604           """           """
605           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}  
606    
607           @rtype: equal to the type of v           :param v: a velovity
608             :return: norm of v
609             :rtype: non-negative ``float``
610           """           """
611       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
612       gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
613    
614        def solve_A(self,u,p):       def getDV(self, p, v, tol):
615           """           """
616           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
617    
618             :param p: a pressure
619             :param v: a initial guess for the value v to return.
620             :return: dv given as *Adv=(f-Av-B^*p)*
621           """           """
622           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
623           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
624         self.getSolverOptionsVelocity().setTolerance(tol)
625         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
626           if self.__stress.isEmpty():           if self.__stress.isEmpty():
627              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
628           else:           else:
629              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
630           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_v.getSolution()
631             return  out
632    
633         def norm_Bv(self,Bv):
634            """
635            Returns Bv (overwrite).
636    
637            :rtype: equal to the type of p
638            :note: boundary conditions on p should be zero!
639            """
640            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
641    
642         def solve_AinvBt(self,p, tol):
643             """
644             Solves *Av=B^*p* with accuracy `tol`
645    
646             :param p: a pressure increment
647             :return: the solution of *Av=B^*p*
648             :note: boundary conditions on v should be zero!
649             """
650             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
651             out=self.__pde_v.getSolution()
652           return  out           return  out
653    
654        def solve_prec(self,p):       def solve_prec(self,Bv, tol):
655           if self.show_details: print "apply preconditioner:"           """
656           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           applies preconditioner for for *BA^{-1}B^** to *Bv*
657           self.__pde_prec.setValue(Y=p)           with accuracy `self.getSubProblemTolerance()`
658           q=self.__pde_prec.getSolution(verbose=self.show_details)  
659           return q           :param Bv: velocity increment
660             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
661             :note: boundary conditions on p are zero.
662             """
663             self.__pde_prec.setValue(Y=Bv)
664         self.getSolverOptionsPressure().setTolerance(tol)
665         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
666             out=self.__pde_prec.getSolution()
667             return out

Legend:
Removed from v.2156  
changed lines
  Added in v.3515

  ViewVC Help
Powered by ViewVC 1.1.26