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

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

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

revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 2793 by gross, Tue Dec 1 06:10:10 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.  
   
         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  
200    
201          M{r=Q^*(g-Qp-v)}      def setAbsoluteTolerance(self,atol=0.):
   
         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()           num_corrections=0
278           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL           converged=False
279           #########################################################################################################################           p=p0
280           #           norm_r=None
281           #   we solve:           while not converged:
282           #                   v=self.getFlux(p, fixed_flux=u0)
283           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )                 Qp=self.__Q(p)
284           #                 norm_v=self.__L2(v)
285           #   residual is                 norm_Qp=self.__L2(Qp)
286           #                 if norm_v == 0.:
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_Qp == 0.:
288           #                       return v,p
289           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                    else:
290           #                      fac=norm_Qp
291           #    we use (g - Qp, v) to represent the residual. not that                 else:
292           #                    if norm_Qp == 0.:
293           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                      fac=norm_v
294           #                    else:
295           #   while the initial residual is                      fac=2./(1./norm_v+1./norm_Qp)
296           #                 ATOL=(atol+rtol*fac)
297           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                 if self.verbose:
298           #                        print "DarcyFlux: L2 norm of v = %e."%norm_v
299           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
300           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
301           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
302           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
303           self.__pde_v.setValue(r=Data())                 if norm_r == None or norm_r>ATOL:
304           # start CG                     if num_corrections>max_num_corrections:
305           r=ArithmeticTuple(d0, v00)                           raise ValueError,"maximum number of correction steps reached."
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)                     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)
307           return r[1],p                     num_corrections+=1
308                   else:
309      def __Aprod_PCG(self,dp):                     converged=True
310            if self.show_details: print "DarcyFlux: Applying operator"           return v,p
311            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)      def __L2(self,v):
312            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
313            self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())  
314            du=self.__pde_v.getSolution(verbose=self.show_details)      def __Q(self,p):
315            return ArithmeticTuple(mQdp,du)            return util.tensor_mult(self.__permeability,util.grad(p))
316    
317        def __Aprod(self,dp):
318              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
319              Qdp=self.__Q(dp)
320              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
321              du=self.__pde_v.getSolution()
322              # self.__pde_v.getOperator().saveMM("proj.mm")
323              return Qdp+du
324        def __inner_GMRES(self,r,s):
325             return util.integrate(util.inner(r,s))
326    
327      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
328           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  
329    
330      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
331            if self.show_details: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
332            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())
333            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
334              return self.__pde_p.getSolution()
335    
336        def getFlux(self,p=None, fixed_flux=Data()):
337            """
338            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
339            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
340            Note that ``g`` and ``f`` are used, see `setValue`.
341    
342            :param p: pressure.
343            :type p: scalar value on the domain (e.g. `Data`).
344            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
345            :type fixed_flux: vector values on the domain (e.g. `Data`).
346            :return: flux
347            :rtype: `Data`
348            :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}*
349                   for the permeability *k_{ij}*
350            """
351        self.setSubProblemTolerance()
352            g=self.__g
353            f=self.__f
354            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
355            if p == None:
356               self.__pde_v.setValue(Y=g)
357            else:
358               self.__pde_v.setValue(Y=g-self.__Q(p))
359            return self.__pde_v.getSolution()
360    
361  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
362        """       """
363        solves       solves
364    
365            -(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}
366                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 368  class StokesProblemCartesian(Homogeneous
368            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
369            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
370    
371        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
372    
373        typical usage:       typical usage:
374    
375              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
376              sp.setTolerance()              sp.setTolerance()
377              sp.initialize(...)              sp.initialize(...)
378              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
379        """       """
380        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
381           """           """
382           initialize the Stokes Problem           initialize the Stokes Problem
383    
384           @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
385           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
386           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
387    
388             :param domain: domain of the problem.
389             :type domain: `Domain`
390           """           """
391           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
392           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
393           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())
394           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
395           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
396           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
397           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
398           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
399    
400           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
401           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
402         self.__pde_proj.setValue(D=1)
403           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
404    
405        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
406             """
407         returns the solver options used  solve the equation for velocity.
408        
409         :rtype: `SolverOptions`
410         """
411         return self.__pde_u.getSolverOptions()
412         def setSolverOptionsVelocity(self, options=None):
413             """
414         set the solver options for solving the equation for velocity.
415        
416         :param options: new solver  options
417         :type options: `SolverOptions`
418         """
419             self.__pde_u.setSolverOptions(options)
420         def getSolverOptionsPressure(self):
421             """
422         returns the solver options used  solve the equation for pressure.
423         :rtype: `SolverOptions`
424         """
425         return self.__pde_prec.getSolverOptions()
426         def setSolverOptionsPressure(self, options=None):
427             """
428         set the solver options for solving the equation for pressure.
429         :param options: new solver  options
430         :type options: `SolverOptions`
431         """
432         self.__pde_prec.setSolverOptions(options)
433    
434         def setSolverOptionsDiv(self, options=None):
435             """
436         set the solver options for solving the equation to project the divergence of
437         the velocity onto the function space of presure.
438        
439         :param options: new solver options
440         :type options: `SolverOptions`
441         """
442         self.__pde_proj.setSolverOptions(options)
443         def getSolverOptionsDiv(self):
444             """
445         returns the solver options for solving the equation to project the divergence of
446         the velocity onto the function space of presure.
447        
448         :rtype: `SolverOptions`
449         """
450         return self.__pde_proj.getSolverOptions()
451    
452         def updateStokesEquation(self, v, p):
453             """
454             updates the Stokes equation to consider dependencies from ``v`` and ``p``
455             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
456             """
457             pass
458         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
459            """
460            assigns new values to the model parameters.
461    
462            :param f: external force
463            :type f: `Vector` object in `FunctionSpace` `Function` or similar
464            :param fixed_u_mask: mask of locations with fixed velocity.
465            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
466            :param eta: viscosity
467            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
468            :param surface_stress: normal surface stress
469            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
470            :param stress: initial stress
471        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
472            """
473            if eta !=None:
474                k=util.kronecker(self.domain.getDim())
475                kk=util.outer(k,k)
476                self.eta=util.interpolate(eta, Function(self.domain))
477            self.__pde_prec.setValue(D=1/self.eta)
478                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
479            if restoration_factor!=None:
480                n=self.domain.getNormal()
481                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
482            if fixed_u_mask!=None:
483                self.__pde_u.setValue(q=fixed_u_mask)
484            if f!=None: self.__f=f
485            if surface_stress!=None: self.__surface_stress=surface_stress
486            if stress!=None: self.__stress=stress
487    
488         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
489          """          """
490          assigns values to the model parameters          assigns values to the model parameters
491    
492          @param f: external force          :param f: external force
493          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
494          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
495          @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
496          @param eta: viscosity          :param eta: viscosity
497          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
498          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
499          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
500          @param stress: initial stress          :param stress: initial stress
501      @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.  
   
502          """          """
503          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)  
         self.__stress=stress  
   
       def B(self,v):  
         """  
         returns div(v)  
         @rtype: equal to the type of p  
   
         @note: boundary conditions on p should be zero!  
         """  
         if self.show_details: print "apply divergence:"  
         self.__pde_proj.setValue(Y=-util.div(v))  
         self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
         return self.__pde_proj.getSolution(verbose=self.show_details)  
   
       def inner_pBv(self,p,Bv):  
          """  
          returns inner product of element p and Bv  (overwrite)  
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @rtype: C{float}  
504    
505           @rtype: equal to the type of p       def Bv(self,v,tol):
506           """           """
507           s0=util.interpolate(p,Function(self.domain))           returns inner product of element p and div(v)
508           s1=util.interpolate(Bv,Function(self.domain))  
509           return util.integrate(s0*s1)           :param v: a residual
510             :return: inner product of element p and div(v)
511             :rtype: ``float``
512             """
513             self.__pde_proj.setValue(Y=-util.div(v))
514         self.getSolverOptionsDiv().setTolerance(tol)
515         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
516             out=self.__pde_proj.getSolution()
517             return out
518    
519         def inner_pBv(self,p,Bv):
520             """
521             returns inner product of element p and Bv=-div(v)
522    
523        def inner_p(self,p0,p1):           :param p: a pressure increment
524             :param Bv: a residual
525             :return: inner product of element p and Bv=-div(v)
526             :rtype: ``float``
527           """           """
528           returns inner product of element p0 and p1  (overwrite)           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
529    
530           @rtype: equal to the type of p       def inner_p(self,p0,p1):
531           """           """
532           s0=util.interpolate(p0/self.eta,Function(self.domain))           Returns inner product of p0 and p1
533           s1=util.interpolate(p1/self.eta,Function(self.domain))  
534             :param p0: a pressure
535             :param p1: a pressure
536             :return: inner product of p0 and p1
537             :rtype: ``float``
538             """
539             s0=util.interpolate(p0,Function(self.domain))
540             s1=util.interpolate(p1,Function(self.domain))
541           return util.integrate(s0*s1)           return util.integrate(s0*s1)
542    
543        def inner_v(self,v0,v1):       def norm_v(self,v):
544           """           """
545           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}  
546    
547           @rtype: equal to the type of v           :param v: a velovity
548             :return: norm of v
549             :rtype: non-negative ``float``
550           """           """
551       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
552    
553        def solve_A(self,u,p):  
554         def getDV(self, p, v, tol):
555           """           """
556           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
557    
558             :param p: a pressure
559             :param v: a initial guess for the value v to return.
560             :return: dv given as *Adv=(f-Av-B^*p)*
561           """           """
562           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
563           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
564         self.getSolverOptionsVelocity().setTolerance(tol)
565         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
566           if self.__stress.isEmpty():           if self.__stress.isEmpty():
567              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
568           else:           else:
569              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
570           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
571           return  out           return  out
572    
573        def solve_prec(self,p):       def norm_Bv(self,Bv):
574           if self.show_details: print "apply preconditioner:"          """
575           self.__pde_prec.setTolerance(self.getSubProblemTolerance())          Returns Bv (overwrite).
576           self.__pde_prec.setValue(Y=p)  
577           q=self.__pde_prec.getSolution(verbose=self.show_details)          :rtype: equal to the type of p
578           return q          :note: boundary conditions on p should be zero!
579            """
580            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
581    
582         def solve_AinvBt(self,p, tol):
583             """
584             Solves *Av=B^*p* with accuracy `tol`
585    
586             :param p: a pressure increment
587             :return: the solution of *Av=B^*p*
588             :note: boundary conditions on v should be zero!
589             """
590             self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
591             out=self.__pde_u.getSolution()
592             return  out
593    
594         def solve_prec(self,Bv, tol):
595             """
596             applies preconditioner for for *BA^{-1}B^** to *Bv*
597             with accuracy `self.getSubProblemTolerance()`
598    
599             :param Bv: velocity increment
600             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
601             :note: boundary conditions on p are zero.
602             """
603             self.__pde_prec.setValue(Y=Bv)
604         self.getSolverOptionsPressure().setTolerance(tol)
605         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
606             out=self.__pde_prec.getSolution()
607             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26