/[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 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 2625 by jfenwick, Fri Aug 21 06:30:25 2009 UTC
# Line 1  Line 1 
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
4  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
5  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
6  #  #
# Line 10  Line 10 
10  #  #
11  ########################################################  ########################################################
12    
13  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
15  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
16  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
17  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
18  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
19  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
20    
21  """  """
22  Some models for flow  Some models for flow
23    
24  @var __author__: name of author  :var __author__: name of author
25  @var __copyright__: copyrights  :var __copyright__: copyrights
26  @var __license__: licence agreement  :var __license__: licence agreement
27  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
28  @var __version__: version  :var __version__: version
29  @var __date__: date of the version  :var __date__: date of the version
30  """  """
31    
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
37  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38    
39  class DarcyFlow(object):  class DarcyFlow(object):
40      """      """
41      solves the problem      solves the problem
42    
43      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
44      M{u_{i,i} = f}      *u_{i,i} = f*
45    
46      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,
47    
48      @note: The problem is solved in a least squares formulation.      :note: The problem is solved in a least squares formulation.
49      """      """
50    
51      def __init__(self, domain,useReduced=False):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
52          """          """
53          initializes the Darcy flux problem          initializes the Darcy flux problem
54          @param domain: domain of the problem          :param domain: domain of the problem
55          @type domain: L{Domain}          :type domain: `Domain`
56        :param useReduced: uses reduced oreder on flux and pressure
57        :type useReduced: ``bool``
58        :param adaptSubTolerance: switches on automatic subtolerance selection
59        :type adaptSubTolerance: ``bool``  
60          """          """
61          self.domain=domain          self.domain=domain
62            if weight == None:
63               s=self.domain.getSize()
64               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65               # self.__l=(3.*util.longestEdge(self.domain))**2
66               # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
67            else:
68               self.__l=weight
69          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
70          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
71          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
72          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
73          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
74          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
75          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
76          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
77          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
78          self.__ATOL= None          self.setTolerance()
79            self.setAbsoluteTolerance()
80        self.__adaptSubTolerance=adaptSubTolerance
81        self.verbose=False
82        def getSolverOptionsFlux(self):
83        """
84        Returns the solver options used to solve the flux problems
85        
86        *(I+D^*D)u=F*
87        
88        :return: `SolverOptions`
89        """
90        return self.__pde_v.getSolverOptions()
91        def setSolverOptionsFlux(self, options=None):
92        """
93        Sets the solver options used to solve the flux problems
94        
95        *(I+D^*D)u=F*
96        
97        If ``options`` is not present, the options are reset to default
98        :param options: `SolverOptions`
99        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
100        """
101        return self.__pde_v.setSolverOptions(options)
102        def getSolverOptionsPressure(self):
103        """
104        Returns the solver options used to solve the pressure problems
105        
106        *(Q^*Q)p=Q^*G*
107        
108        :return: `SolverOptions`
109        """
110        return self.__pde_p.getSolverOptions()
111        def setSolverOptionsPressure(self, options=None):
112        """
113        Sets the solver options used to solve the pressure problems
114        
115        *(Q^*Q)p=Q^*G*
116        
117        If ``options`` is not present, the options are reset to default
118        :param options: `SolverOptions`
119        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
120        """
121        return self.__pde_p.setSolverOptions(options)
122    
123      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124          """          """
125          assigns values to model parameters          assigns values to model parameters
126    
127          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
128          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `Data`)
129          @param g: flux sources/sinks          :param g: flux sources/sinks
130          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `Data`)
131          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :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. L{Data})          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
133          @param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
134          @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. `Data`)
135          @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
136                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
137                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
138          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
139    
140          @note: the values of parameters which are not set by calling C{setValue} are not altered.          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
141          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
142                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
143                 is along the M{x_i} axis.                 is along the *x_i* axis.
144          """          """
145          if f !=None:          if f !=None:
146             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
147             if f.isEmpty():             if f.isEmpty():
148                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
149             else:             else:
150                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
151             self.f=f             self.__f=f
152          if g !=None:            if g !=None:
153             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
154             if g.isEmpty():             if g.isEmpty():
155               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 121  class DarcyFlow(object): Line 175  class DarcyFlow(object):
175             self.__permeability=perm             self.__permeability=perm
176             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
177    
178        def setTolerance(self,rtol=1e-4):
     def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):  
179          """          """
180          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
         on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
         Note that C{g} and C{f} are used, see L{setValue}.  
           
         @param p: pressure.  
         @type p: scalar value on the domain (e.g. L{Data}).  
         @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.  
         @type fixed_flux: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: flux  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
         """  
         self.__pde_v.setTolerance(tol)  
         g=self.__g  
         f=self.__f  
         self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
         else:  
            self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))  
         return self.__pde_v.getSolution(verbose=show_details)  
181    
182      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
183    
184            where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
185    
186            :param rtol: relative tolerance for the pressure
187            :type rtol: non-negative ``float``
188          """          """
189          returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
190          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
191          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
192                def getTolerance(self):
         @param v: flux.  
         @type v: vector-valued on the domain (e.g. L{Data}).  
         @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.  
         @type fixed_pressure: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: pressure  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
193          """          """
194          self.__pde_v.setTolerance(tol)          returns the relative tolerance
         g=self.__g  
         self.__pde_p.setValue(r=fixed_pressure)  
         if v == None:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))  
         else:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))  
         return self.__pde_p.getSolution(verbose=show_details)  
195    
196      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          :return: current relative tolerance
197            :rtype: ``float``
198          """          """
199          set the tolerance C{ATOL} used to terminate the solution process. It is used          return self.__rtol
   
         M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}  
   
         where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.  
200    
201          The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual      def setAbsoluteTolerance(self,atol=0.):
   
         M{r=Q^*(g-Qp-v)}  
   
         the condition  
   
         M{<(Q^*Q)^{-1} r,r> <= ATOL}  
   
         holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}  
   
         @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 p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the  
         L{getPressure} method or use  the value from a previous time step.  
         @type p_ref: scalar value on the domain (e.g. L{Data}).  
         @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the  
         L{getFlux} method or use  the value from a previous time step.  
         @type v_ref: vector-valued on the domain (e.g. L{Data}).  
         @return: used absolute tolerance.  
         @rtype: positive C{float}  
202          """          """
203          g=self.__g          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
204          if not v_ref == None:  
205             f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
206          else:  
207             f1=util.integrate(util.length(util.interpolate(g))**2)          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
208          if not p_ref == None:  
209             f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)          :param atol: absolute tolerance for the pressure
210          else:          :type atol: non-negative ``float``
211             f2=0          """
212          self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))          if atol<0:
213          if self.__ATOL<=0:              raise ValueError,"Absolute tolerance needs to be non-negative."
214             raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL          self.__atol=atol
215          return self.__ATOL      def getAbsoluteTolerance(self):
216                   """
217      def getTolerance(self):         returns the absolute tolerance
218          """        
219          returns the current tolerance.         :return: current absolute tolerance
220             :rtype: ``float``
221          @return: used absolute tolerance.         """
222          @rtype: positive C{float}         return self.__atol
223          """      def getSubProblemTolerance(self):
224          if self.__ATOL==None:      """
225             raise ValueError,"no tolerance is defined."      Returns a suitable subtolerance
226          return self.__ATOL      @type: ``float``
227        """
228        return max(util.EPSILON**(0.75),self.getTolerance()**2)
229        def setSubProblemTolerance(self):
230             """
231             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
232             """
233         if self.__adaptSubTolerance:
234             sub_tol=self.getSubProblemTolerance()
235                 self.getSolverOptionsFlux().setTolerance(sub_tol)
236             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
237             self.getSolverOptionsPressure().setTolerance(sub_tol)
238             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
239             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
240    
241      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
242           """           """
243           solves the problem.           solves the problem.
244    
245           The iteration is terminated if the residual norm is less then self.getTolerance().           The iteration is terminated if the residual norm is less then self.getTolerance().
246    
247           @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.           :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
248           @type u0: vector value on the domain (e.g. L{Data}).           :type u0: vector value on the domain (e.g. `Data`).
249           @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.           :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
250           @type p0: scalar value on the domain (e.g. L{Data}).           :type p0: scalar value on the domain (e.g. `Data`).
251           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}           :param verbose: if set some information on iteration progress are printed
252           @type sub_rtol: positive-negative C{float}           :type verbose: ``bool``
253           @param verbose: if set some information on iteration progress are printed           :return: flux and pressure
254           @type verbose: C{bool}           :rtype: ``tuple`` of `Data`.
255           @param show_details:  if set information on the subiteration process are printed.  
256           @type show_details: C{bool}           :note: The problem is solved as a least squares form
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
   
          @note: The problem is solved as a least squares form  
257    
258           M{(I+D^*D)u+Qp=D^*f+g}           *(I+D^*D)u+Qp=D^*f+g*
259           M{Q^*u+Q^*Qp=Q^*g}           *Q^*u+Q^*Qp=Q^*g*
260    
261           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
262           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
263    
264           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
265    
266           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
267    
268           M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure           *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
269            
270           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
271           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
272           """           """
273           self.verbose=verbose           self.verbose=verbose
274           self.show_details= show_details and self.verbose           rtol=self.getTolerance()
275           self.__pde_v.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
276           self.__pde_p.setTolerance(sub_rtol)       self.setSubProblemTolerance()
277           ATOL=self.getTolerance()      
278           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL           num_corrections=0
279           #########################################################################################################################           converged=False
280           #           p=p0
281           #   we solve:           norm_r=None
282           #             while not converged:
283           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )                 v=self.getFlux(p, fixed_flux=u0)
284           #                 Qp=self.__Q(p)
285           #   residual is                 norm_v=self.__L2(v)
286           #                 norm_Qp=self.__L2(Qp)
287           #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)                 if norm_v == 0.:
288           #                    if norm_Qp == 0.:
289           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                       return v,p
290           #                    else:
291           #    we use (g - Qp, v) to represent the residual. not that                      fac=norm_Qp
292           #                 else:
293           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                    if norm_Qp == 0.:
294           #                      fac=norm_v
295           #   while the initial residual is                    else:
296           #                      fac=2./(1./norm_v+1./norm_Qp)
297           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                 ATOL=(atol+rtol*fac)
298           #                   if self.verbose:
299           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                      print "DarcyFlux: L2 norm of v = %e."%norm_v
300           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
301           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
302           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
303           self.__pde_v.setValue(r=Data())                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
304           # start CG                 if norm_r == None or norm_r>ATOL:
305           r=ArithmeticTuple(d0, v00)                     if num_corrections>max_num_corrections:
306           p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                           raise ValueError,"maximum number of correction steps reached."
307           return r[1],p                     p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
308                       num_corrections+=1
309      def __Aprod_PCG(self,dp):                 else:
310            if self.show_details: print "DarcyFlux: Applying operator"                     converged=True
311            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)           return v,p
312            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))      def __L2(self,v):
313            self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
314            du=self.__pde_v.getSolution(verbose=self.show_details)  
315            return ArithmeticTuple(mQdp,du)      def __Q(self,p):
316              return util.tensor_mult(self.__permeability,util.grad(p))
317    
318        def __Aprod(self,dp):
319              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
320              Qdp=self.__Q(dp)
321              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
322              du=self.__pde_v.getSolution()
323              # self.__pde_v.getOperator().saveMM("proj.mm")
324              return Qdp+du
325        def __inner_GMRES(self,r,s):
326             return util.integrate(util.inner(r,s))
327    
328      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
329           a=util.tensor_mult(self.__permeability,util.grad(p))           return util.integrate(util.inner(self.__Q(p), r))
          f0=util.integrate(util.inner(a,r[0]))  
          f1=util.integrate(util.inner(a,r[1]))  
          # print "__inner_PCG:",f0,f1,"->",f0-f1  
          return f0-f1  
330    
331      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
332            if self.show_details: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
333            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
334            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
335              return self.__pde_p.getSolution()
336    
337        def getFlux(self,p=None, fixed_flux=Data()):
338            """
339            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
340            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
341            Note that ``g`` and ``f`` are used, see `setValue`.
342    
343            :param p: pressure.
344            :type p: scalar value on the domain (e.g. `Data`).
345            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
346            :type fixed_flux: vector values on the domain (e.g. `Data`).
347            :return: flux
348            :rtype: `Data`
349            :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
350                   for the permeability *k_{ij}*
351            """
352        self.setSubProblemTolerance()
353            g=self.__g
354            f=self.__f
355            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
356            if p == None:
357               self.__pde_v.setValue(Y=g)
358            else:
359               self.__pde_v.setValue(Y=g-self.__Q(p))
360            return self.__pde_v.getSolution()
361    
362  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
363        """       """
364        solves       solves
365    
366            -(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}
367                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 369  class StokesProblemCartesian(Homogeneous
369            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
370            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
371    
372        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
373    
374        typical usage:       typical usage:
375    
376              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
377              sp.setTolerance()              sp.setTolerance()
378              sp.initialize(...)              sp.initialize(...)
379              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
380        """       """
381        def __init__(self,domain,**kwargs):       def __init__(self,domain,adaptSubTolerance=True, **kwargs):
382           """           """
383           initialize the Stokes Problem           initialize the Stokes Problem
384    
385           @param domain: domain of the problem. The approximation order needs to be two.           :param domain: domain of the problem. The approximation order needs to be two.
386           @type domain: L{Domain}           :type domain: `Domain`
387           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.       :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
388         :type adaptSubTolerance: ``bool``
389             :warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
390           """           """
391           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
392           self.domain=domain           self.domain=domain
393           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
394           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
395           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
396           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
397           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
398           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
399           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
400    
401           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
402           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
403         self.__pde_proj.setValue(D=1)
404           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
405    
406        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
407             """
408         returns the solver options used  solve the equation for velocity.
409        
410         :rtype: `SolverOptions`
411         """
412         return self.__pde_u.getSolverOptions()
413         def setSolverOptionsVelocity(self, options=None):
414             """
415         set the solver options for solving the equation for velocity.
416        
417         :param options: new solver  options
418         :type options: `SolverOptions`
419         """
420             self.__pde_u.setSolverOptions(options)
421         def getSolverOptionsPressure(self):
422             """
423         returns the solver options used  solve the equation for pressure.
424         :rtype: `SolverOptions`
425         """
426         return self.__pde_prec.getSolverOptions()
427         def setSolverOptionsPressure(self, options=None):
428             """
429         set the solver options for solving the equation for pressure.
430         :param options: new solver  options
431         :type options: `SolverOptions`
432         """
433         self.__pde_prec.setSolverOptions(options)
434    
435         def setSolverOptionsDiv(self, options=None):
436             """
437         set the solver options for solving the equation to project the divergence of
438         the velocity onto the function space of presure.
439        
440         :param options: new solver options
441         :type options: `SolverOptions`
442         """
443         self.__pde_prec.setSolverOptions(options)
444         def getSolverOptionsDiv(self):
445             """
446         returns the solver options for solving the equation to project the divergence of
447         the velocity onto the function space of presure.
448        
449         :rtype: `SolverOptions`
450         """
451         return self.__pde_prec.getSolverOptions()
452         def setSubProblemTolerance(self):
453             """
454         Updates the tolerance for subproblems
455             """
456         if self.adaptSubTolerance():
457                 sub_tol=self.getSubProblemTolerance()
458             self.getSolverOptionsDiv().setTolerance(sub_tol)
459             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
460             self.getSolverOptionsPressure().setTolerance(sub_tol)
461             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
462             self.getSolverOptionsVelocity().setTolerance(sub_tol)
463             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
464            
465    
466         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
467          """          """
468          assigns values to the model parameters          assigns values to the model parameters
469    
470          @param f: external force          :param f: external force
471          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
472          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
473          @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
474          @param eta: viscosity          :param eta: viscosity
475          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
476          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
477          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
478          @param stress: initial stress          :param stress: initial stress
479      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
480          @note: All values needs to be set.          :note: All values needs to be set.
   
481          """          """
482          self.eta=eta          self.eta=eta
483          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
484      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
485          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
486          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
487              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
488              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
489      self.__pde_prec.setValue(D=1/self.eta)          n=self.domain.getNormal()
490          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)      self.__pde_prec.setValue(D=1/self.eta)
491            self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))
492            self.__f=f
493            self.__surface_stress=surface_stress
494          self.__stress=stress          self.__stress=stress
495    
496        def B(self,v):       def Bv(self,v):
497          """           """
498          returns div(v)           returns inner product of element p and div(v)
         @rtype: equal to the type of p  
499    
500          @note: boundary conditions on p should be zero!           :param v: a residual
501          """           :return: inner product of element p and div(v)
502          if self.show_details: print "apply divergence:"           :rtype: ``float``
503          self.__pde_proj.setValue(Y=-util.div(v))           """
504          self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_proj.setValue(Y=-util.div(v))
505          return self.__pde_proj.getSolution(verbose=self.show_details)           return self.__pde_proj.getSolution()
506    
507        def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
508           """           """
509           returns inner product of element p and Bv  (overwrite)           returns inner product of element p and Bv=-div(v)
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @rtype: C{float}  
510    
511           @rtype: equal to the type of p           :param p: a pressure increment
512             :param Bv: a residual
513             :return: inner product of element p and Bv=-div(v)
514             :rtype: ``float``
515           """           """
516           s0=util.interpolate(p,Function(self.domain))           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
          s1=util.interpolate(Bv,Function(self.domain))  
          return util.integrate(s0*s1)  
517    
518        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
519           """           """
520           returns inner product of element p0 and p1  (overwrite)           Returns inner product of p0 and p1
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
521    
522           @rtype: equal to the type of p           :param p0: a pressure
523             :param p1: a pressure
524             :return: inner product of p0 and p1
525             :rtype: ``float``
526           """           """
527           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
528           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
529           return util.integrate(s0*s1)           return util.integrate(s0*s1)
530    
531        def inner_v(self,v0,v1):       def norm_v(self,v):
532           """           """
533           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}  
534    
535           @rtype: equal to the type of v           :param v: a velovity
536             :return: norm of v
537             :rtype: non-negative ``float``
538           """           """
539       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
540    
541        def solve_A(self,u,p):       def getV(self, p, v0):
542           """           """
543           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
544    
545             :param p: a pressure
546             :param v0: a initial guess for the value v to return.
547             :return: v given as *v= A^{-1} (f-B^*p)*
548           """           """
549           if self.show_details: print "solve for velocity:"           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
550           if self.__stress.isEmpty():           if self.__stress.isEmpty():
551              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain))
552           else:           else:
553              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
554           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
555             return  out
556    
557         def norm_Bv(self,Bv):
558            """
559            Returns Bv (overwrite).
560    
561            :rtype: equal to the type of p
562            :note: boundary conditions on p should be zero!
563            """
564            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
565    
566         def solve_AinvBt(self,p):
567             """
568             Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`
569    
570             :param p: a pressure increment
571             :return: the solution of *Av=B^*p*
572             :note: boundary conditions on v should be zero!
573             """
574             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
575             out=self.__pde_u.getSolution()
576           return  out           return  out
577    
578        def solve_prec(self,p):       def solve_prec(self,Bv):
579           if self.show_details: print "apply preconditioner:"           """
580           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           applies preconditioner for for *BA^{-1}B^** to *Bv*
581           self.__pde_prec.setValue(Y=p)           with accuracy `self.getSubProblemTolerance()`
582           q=self.__pde_prec.getSolution(verbose=self.show_details)  
583           return q           :param Bv: velocity increment
584             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
585             :note: boundary conditions on p are zero.
586             """
587             self.__pde_prec.setValue(Y=Bv)
588             return self.__pde_prec.getSolution()

Legend:
Removed from v.2208  
changed lines
  Added in v.2625

  ViewVC Help
Powered by ViewVC 1.1.26