/[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 2108 by gross, Fri Nov 28 05:09:23 2008 UTC revision 3990 by caltinay, Tue Sep 25 05:03:20 2012 UTC
# Line 1  Line 1 
1  ########################################################  # -*- coding: utf-8 -*-
2    ##############################################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2012 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
21  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
22    
23  """  """
24  Some models for flow  Some models for flow
25    
26  @var __author__: name of author  :var __author__: name of author
27  @var __copyright__: copyrights  :var __copyright__: copyrights
28  @var __license__: licence agreement  :var __license__: licence agreement
29  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
30  @var __version__: version  :var __version__: version
31  @var __date__: date of the version  :var __date__: date of the version
32  """  """
33    
34  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
35    
36  from escript import *  from . import escript
37  import util  from . import util
38  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
39  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
40    
41  class DarcyFlow(object):  class DarcyFlow(object):
42      """     """
43      solves the problem     solves the problem
44      
45      M{u_i+k_{ij}*p_{,j} = g_i}     *u_i+k_{ij}*p_{,j} = g_i*
46      M{u_{i,i} = f}     *u_{i,i} = f*
47      
48       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
49      
50       :cvar EVAL: direct pressure gradient evaluation for flux
51       :cvar POST: global postprocessing of flux by solving the PDE *K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j*
52                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
53       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
54       """
55       EVAL="EVAL"
56       SIMPLE="EVAL"
57       POST="POST"
58       SMOOTH="SMOOTH"
59       def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
60          """
61          initializes the Darcy flux problem.
62    
63      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,        :param domain: domain of the problem
64          :type domain: `Domain`
65          :param useReduced: uses reduced oreder on flux and pressure
66          :type useReduced: ``bool``
67          :param solver: solver method
68          :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
69          :param verbose: if ``True`` some information on the iteration progress are printed.
70          :type verbose: ``bool``
71          :param w: weighting factor for `DarcyFlow.POST` solver
72          :type w: ``float``
73          
74          """
75          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
76              raise ValueError("unknown solver %d."%solver)
77    
78      @note: The problem is solved in a least squares formulation.        self.domain=domain
79      """        self.solver=solver
80          self.useReduced=useReduced
81          self.verbose=verbose
82          self.l=None
83          self.w=None
84        
85          self.__pde_p=LinearSinglePDE(domain)
86          self.__pde_p.setSymmetryOn()
87          if self.useReduced: self.__pde_p.setReducedOrderOn()
88    
89          if self.solver  == self.EVAL:
90             self.__pde_v=None
91             if self.verbose: print("DarcyFlow: simple solver is used.")
92    
93          elif self.solver  == self.POST:
94             if util.inf(w)<0.:
95                raise ValueError("Weighting factor must be non-negative.")
96             if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
97             self.__pde_v=LinearPDESystem(domain)
98             self.__pde_v.setSymmetryOn()
99             if self.useReduced: self.__pde_v.setReducedOrderOn()
100             self.w=w
101             x=self.domain.getX()
102             self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
103             #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
104    
105          elif self.solver  == self.SMOOTH:
106             self.__pde_v=LinearPDESystem(domain)
107             self.__pde_v.setSymmetryOn()
108             if self.useReduced: self.__pde_v.setReducedOrderOn()
109             if self.verbose: print("DarcyFlow: flux smoothing is used.")
110             self.w=0
111    
112          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
113          self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
114          self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
115          self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
116          self.ref_point_id=None
117          self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
119          self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
120          self.perm_scale=1.
121        
122            
123       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124          """
125          assigns values to model parameters
126    
127      def __init__(self, domain):        :param f: volumetic sources/sinks
128          """        :type f: scalar value on the domain (e.g. `escript.Data`)
129          initializes the Darcy flux problem        :param g: flux sources/sinks
130          @param domain: domain of the problem        :type g: vector values on the domain (e.g. `escript.Data`)
131          @type domain: L{Domain}        :param location_of_fixed_pressure: mask for locations where pressure is fixed
132          """        :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
133          self.domain=domain        :param location_of_fixed_flux:  mask for locations where flux is fixed.
134          self.__pde_v=LinearPDESystem(domain)        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
135          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
136          self.__pde_v.setSymmetryOn()        :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
137          self.__pde_p=LinearSinglePDE(domain)  
138          self.__pde_p.setSymmetryOn()        :note: the values of parameters which are not set by calling ``setValue`` are not altered.
139          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        :note: at any point on the boundary of the domain the pressure
140          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               (``location_of_fixed_pressure`` >0) or the normal component of the
141                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
142                 is along the *x_i* axis.
143    
144      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):        """
145          """        if location_of_fixed_pressure!=None:
146          assigns values to model parameters             self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
147               self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
148               if not self.location_of_fixed_pressure.getTupleForGlobalDataPoint(*self.ref_point_id)[0] > 0: raise ValueError("pressure needs to be fixed at least one point.")
149               self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
150               if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
151               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
152          if location_of_fixed_flux!=None:
153              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
154              if not self.__pde_v == None:
155                  self.__pde_v.setValue(q=self.location_of_fixed_flux)
156                
157          if permeability!=None:
158        
159             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
160             self.perm_scale=util.Lsup(util.length(perm))
161             if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
162             perm=perm*(1./self.perm_scale)
163            
164             if perm.getRank()==0:
165    
166          @param f: volumetic sources/sinks              perm_inv=(1./perm)
167          @type f: scalar value on the domain (e.g. L{Data})              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
168          @param g: flux sources/sinks              perm=perm*util.kronecker(self.domain.getDim())
169          @type g: vector values on the domain (e.g. L{Data})          
170          @param location_of_fixed_pressure: mask for locations where pressure is fixed          
171          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})           elif perm.getRank()==2:
172          @param location_of_fixed_flux:  mask for locations where flux is fixed.              perm_inv=util.inverse(perm)
173          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})           else:
174          @param permeability: permeability tensor. If scalar C{s} is given the tensor with              raise ValueError("illegal rank of permeability.")
175                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with          
176                               C{v} on the main diagonal is used.           self.__permeability=perm
177          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})           self.__permeability_inv=perm_inv
178        
179          @note: the values of parameters which are not set by calling C{setValue} are not altered.           #====================
180          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)           self.__pde_p.setValue(A=self.__permeability)
181                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal           if self.solver  == self.EVAL:
182                 is along the M{x_i} axis.                pass # no extra work required
183          """           elif self.solver  == self.POST:
184          if f !=None:                k=util.kronecker(self.domain.getDim())
185             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
186             if f.isEmpty():                #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
187                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
188             else:           elif self.solver  == self.SMOOTH:
189                 if f.getRank()>0: raise ValueError,"illegal rank of f."              self.__pde_v.setValue(D=self.__permeability_inv)
190             self.f=f  
191          if g !=None:          if g != None:
192             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
193             if g.isEmpty():          if g.isEmpty():
194               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195             else:          else:
196               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
197                 raise ValueError,"illegal shape of g"          self.__g=g
198             self.__g=g          self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
199            self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
200          if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)        if f !=None:
201          if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
202             if f.isEmpty():      
203          if permeability!=None:               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))           else:
205             if perm.getRank()==0:               if f.getRank()>0: raise ValueError("illegal rank of f.")
206                 perm=perm*util.kronecker(self.domain.getDim())           self.__f=f
            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))  
207    
208       def getSolverOptionsFlux(self):
209          """
210          Returns the solver options used to solve the flux problems
211          :return: `SolverOptions`
212          """
213          if self.__pde_v == None:
214              return None
215          else:
216              return self.__pde_v.getSolverOptions()
217          
218       def setSolverOptionsFlux(self, options=None):
219          """
220          Sets the solver options used to solve the flux problems
221          If ``options`` is not present, the options are reset to default
222          :param options: `SolverOptions`
223          """
224          if not self.__pde_v == None:
225              self.__pde_v.setSolverOptions(options)
226        
227       def getSolverOptionsPressure(self):
228          """
229          Returns the solver options used to solve the pressure problems
230          :return: `SolverOptions`
231          """
232          return self.__pde_p.getSolverOptions()
233          
234       def setSolverOptionsPressure(self, options=None):
235          """
236          Sets the solver options used to solve the pressure problems
237          If ``options`` is not present, the options are reset to default
238          
239          :param options: `SolverOptions`
240          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
241          """
242          return self.__pde_p.setSolverOptions(options)
243          
244       def solve(self, u0, p0):
245          """
246          solves the problem.
247          
248          :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.
249          :type u0: vector value on the domain (e.g. `escript.Data`).
250          :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.
251          :type p0: scalar value on the domain (e.g. `escript.Data`).
252          :return: flux and pressure
253          :rtype: ``tuple`` of `escript.Data`.
254    
255      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8):        """
256          p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
257          if self.ref_point_id == None:
258              p_ref=0
259          else:
260              p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
261          p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
262          g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
263          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
264                                Y=self.__f * 1./self.perm_scale,
265                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
266                                r=p0 - p0_hydrostatic)
267          pp=self.__pde_p.getSolution()
268          u = self._getFlux(pp, u0)
269          return u,pp + p0_hydrostatic
270          
271       def getFlux(self,p, u0=None):
272          """          """
273          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``
274          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
275          Note that C{g} and C{f} are used, L{setValue}.          Notice that ``g`` is used, see `setValue`.
276            
277          @param p: pressure.          :param p: pressure.
278          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `escript.Data`).
279          @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``.
280          @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``
281          @param tol: relative tolerance to be used.          :return: flux
282          @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}}  
283          """          """
284          self.__pde_v.setTolerance(tol)          p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
285          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=boundary_flux)          if self.ref_point_id == None:
286          return self.__pde_v.getSolution()              p_ref=0
287            else:
288      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):              p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
289           """          p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
290           solves the problem.          return self._getFlux(p-p_hydrostatic, u0)
   
          The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where  
          C{|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 set some information on iteration progress are printed  
          @type verbose: C{bool}  
          @param show_details:  if set information on the subiteration process are printed.  
          @type show_details: C{bool}  
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
   
          @note: The problem is solved as a least squares form  
   
          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 form the problem by setting  
   
          M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux  
   
          form 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 p=p0  on 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} needs 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)  
          p2=p0*self.__pde_p.getCoefficient("q")  
          u2=u0*self.__pde_v.getCoefficient("q")  
          g=self.__g-u2-util.tensor_mult(self.__permeability,util.grad(p2))  
          f=self.__f-util.div(u2)  
          self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())  
          dv=self.__pde_v.getSolution(verbose=show_details)  
          self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))  
          self.__pde_p.setValue(r=Data())  
          dp=self.__pde_p.getSolution(verbose=self.show_details)  
          norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))  
          if norm_rhs<0:  
              raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."  
          ATOL=util.sqrt(norm_rhs)*rtol +atol  
          if not ATOL>0:  
              raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)  
          rhs=ArithmeticTuple(g,dv)  
          dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0-p2, verbose=self.verbose, initial_guess=True)  
          return u2+r[1],p2+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))  
          return util.integrate(util.inner(a,r[0]-r[1]))  
   
     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)  
291    
292       def _getFlux(self, pp, u0=None):
293            """
294            returns the flux for a given pressure ``pp`` where the flux is equal to
295            ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
296            `setValue`). Notice that ``g`` is used, see `setValue`.
297    
298            :param pp: pressure.
299            :type pp: scalar value on the domain (i.e. `escript.Data`).
300            :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
301            :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
302            :return: flux
303            :rtype: `escript.Data`
304            """
305            if self.solver  == self.EVAL:
306               u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
307            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
308                self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
309                print
310                if u0 == None:
311                   self.__pde_v.setValue(r=escript.Data())
312                else:
313                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
314                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
315                u= self.__pde_v.getSolution() * self.perm_scale
316            return u
317          
318  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
319        """       """
320        solves       solves
321    
322            -(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}
323                  u_{i,i}=0                  u_{i,i}=0
# Line 230  class StokesProblemCartesian(Homogeneous Line 325  class StokesProblemCartesian(Homogeneous
325            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
326            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
327    
328        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
329    
330        typical usage:       typical usage:
331    
332              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
333              sp.setTolerance()              sp.setTolerance()
334              sp.initialize(...)              sp.initialize(...)
335              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
336        """              sp.setStokesEquation(...) # new values for some parameters
337        def __init__(self,domain,**kwargs):              v1,p1=sp.solve(v,p)
338         """
339         def __init__(self,domain,**kwargs):
340           """           """
341           initialize the Stokes Problem           initialize the Stokes Problem
342    
343           @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
344           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
345           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
346    
347             :param domain: domain of the problem.
348             :type domain: `Domain`
349           """           """
350           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
351           self.domain=domain           self.domain=domain
352           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
353           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
354           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
355           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
356           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
357           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
358    
359           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
360           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
361             self.__pde_proj.setValue(D=1)
362           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
363    
364        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
365             """
366         returns the solver options used  solve the equation for velocity.
367        
368         :rtype: `SolverOptions`
369         """
370             return self.__pde_v.getSolverOptions()
371         def setSolverOptionsVelocity(self, options=None):
372             """
373         set the solver options for solving the equation for velocity.
374        
375         :param options: new solver  options
376         :type options: `SolverOptions`
377         """
378             self.__pde_v.setSolverOptions(options)
379         def getSolverOptionsPressure(self):
380             """
381         returns the solver options used  solve the equation for pressure.
382         :rtype: `SolverOptions`
383         """
384             return self.__pde_prec.getSolverOptions()
385         def setSolverOptionsPressure(self, options=None):
386             """
387         set the solver options for solving the equation for pressure.
388         :param options: new solver  options
389         :type options: `SolverOptions`
390         """
391             self.__pde_prec.setSolverOptions(options)
392    
393         def setSolverOptionsDiv(self, options=None):
394             """
395         set the solver options for solving the equation to project the divergence of
396         the velocity onto the function space of presure.
397        
398         :param options: new solver options
399         :type options: `SolverOptions`
400         """
401             self.__pde_proj.setSolverOptions(options)
402         def getSolverOptionsDiv(self):
403             """
404         returns the solver options for solving the equation to project the divergence of
405         the velocity onto the function space of presure.
406        
407         :rtype: `SolverOptions`
408         """
409             return self.__pde_proj.getSolverOptions()
410    
411         def updateStokesEquation(self, v, p):
412             """
413             updates the Stokes equation to consider dependencies from ``v`` and ``p``
414             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
415             """
416             pass
417         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
418          """          """
419          assigns values to the model parameters          assigns new values to the model parameters.
420    
421          @param f: external force          :param f: external force
422          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
423          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
424          @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
425          @param eta: viscosity          :param eta: viscosity
426          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
427          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
428          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
429          @param stress: initial stress          :param stress: initial stress
430      @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.  
   
431          """          """
432          self.eta=eta          if eta !=None:
433          A =self.__pde_u.createCoefficient("A")              k=util.kronecker(self.domain.getDim())
434      self.__pde_u.setValue(A=Data())              kk=util.outer(k,k)
435          for i in range(self.domain.getDim()):              self.eta=util.interpolate(eta, escript.Function(self.domain))
436          for j in range(self.domain.getDim()):              self.__pde_prec.setValue(D=1/self.eta)
437              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)))
438              A[i,j,i,j] += 1.          if restoration_factor!=None:
439      self.__pde_prec.setValue(D=1/self.eta)              n=self.domain.getNormal()
440          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))
441          self.__stress=stress          if fixed_u_mask!=None:
442                self.__pde_v.setValue(q=fixed_u_mask)
443            if f!=None: self.__f=f
444            if surface_stress!=None: self.__surface_stress=surface_stress
445            if stress!=None: self.__stress=stress
446    
447        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):
448          """          """
449          returns div(v)          assigns values to the model parameters
         @rtype: equal to the type of p  
450    
451          @note: boundary conditions on p should be zero!          :param f: external force
452            :type f: `Vector` object in `FunctionSpace` `Function` or similar
453            :param fixed_u_mask: mask of locations with fixed velocity.
454            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
455            :param eta: viscosity
456            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
457            :param surface_stress: normal surface stress
458            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
459            :param stress: initial stress
460            :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
461          """          """
462          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)  
463    
464        def inner_pBv(self,p,Bv):       def Bv(self,v,tol):
465           """           """
466           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}  
467    
468           @rtype: equal to the type of p           :param v: a residual
469             :return: inner product of element p and div(v)
470             :rtype: ``float``
471           """           """
472           s0=util.interpolate(p,Function(self.domain))           self.__pde_proj.setValue(Y=-util.div(v))
473           s1=util.interpolate(Bv,Function(self.domain))           self.getSolverOptionsDiv().setTolerance(tol)
474           return util.integrate(s0*s1)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
475             out=self.__pde_proj.getSolution()
476             return out
477    
478        def inner_p(self,p0,p1):       def inner_pBv(self,p,Bv):
479           """           """
480           returns inner product of element p0 and p1  (overwrite)           returns inner product of element p and Bv=-div(v)
481            
482           @type p0: equal to the type of p           :param p: a pressure increment
483           @type p1: equal to the type of p           :param Bv: a residual
484           @rtype: C{float}           :return: inner product of element p and Bv=-div(v)
485             :rtype: ``float``
486             """
487             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
488    
489           @rtype: equal to the type of p       def inner_p(self,p0,p1):
490           """           """
491           s0=util.interpolate(p0/self.eta,Function(self.domain))           Returns inner product of p0 and p1
492           s1=util.interpolate(p1/self.eta,Function(self.domain))  
493             :param p0: a pressure
494             :param p1: a pressure
495             :return: inner product of p0 and p1
496             :rtype: ``float``
497             """
498             s0=util.interpolate(p0, escript.Function(self.domain))
499             s1=util.interpolate(p1, escript.Function(self.domain))
500           return util.integrate(s0*s1)           return util.integrate(s0*s1)
501    
502        def inner_v(self,v0,v1):       def norm_v(self,v):
503           """           """
504           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}  
505    
506           @rtype: equal to the type of v           :param v: a velovity
507             :return: norm of v
508             :rtype: non-negative ``float``
509           """           """
510       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
511       gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
512    
513        def solve_A(self,u,p):       def getDV(self, p, v, tol):
514           """           """
515           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p
516    
517             :param p: a pressure
518             :param v: a initial guess for the value v to return.
519             :return: dv given as *Adv=(f-Av-B^*p)*
520           """           """
521           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
522           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
523             self.getSolverOptionsVelocity().setTolerance(tol)
524             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
525           if self.__stress.isEmpty():           if self.__stress.isEmpty():
526              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)))
527           else:           else:
528              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)))
529           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_v.getSolution()
530             return  out
531    
532         def norm_Bv(self,Bv):
533            """
534            Returns Bv (overwrite).
535    
536            :rtype: equal to the type of p
537            :note: boundary conditions on p should be zero!
538            """
539            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
540    
541         def solve_AinvBt(self,p, tol):
542             """
543             Solves *Av=B^*p* with accuracy `tol`
544    
545             :param p: a pressure increment
546             :return: the solution of *Av=B^*p*
547             :note: boundary conditions on v should be zero!
548             """
549             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
550             out=self.__pde_v.getSolution()
551           return  out           return  out
552    
553        def solve_prec(self,p):       def solve_prec(self,Bv, tol):
554           if self.show_details: print "apply preconditioner:"           """
555           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           applies preconditioner for for *BA^{-1}B^** to *Bv*
556           self.__pde_prec.setValue(Y=p)           with accuracy ``self.getSubProblemTolerance()``
557           q=self.__pde_prec.getSolution(verbose=self.show_details)  
558           return q           :param Bv: velocity increment
559             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
560             :note: boundary conditions on p are zero.
561             """
562             self.__pde_prec.setValue(Y=Bv)
563             self.getSolverOptionsPressure().setTolerance(tol)
564             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
565             out=self.__pde_prec.getSolution()
566             return out

Legend:
Removed from v.2108  
changed lines
  Added in v.3990

  ViewVC Help
Powered by ViewVC 1.1.26