/[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 2349 by gross, Mon Mar 30 08:14:23 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"""
# Line 21  __url__="https://launchpad.net/escript-f Line 21  __url__="https://launchpad.net/escript-f
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, GMRES  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          self.__l=util.longestEdge(self.domain)**2          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()
# Line 67  class DarcyFlow(object): Line 77  class DarcyFlow(object):
77          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
78          self.setTolerance()          self.setTolerance()
79          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
80          self.setSubProblemTolerance()      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"))
# Line 126  class DarcyFlow(object): Line 177  class DarcyFlow(object):
177    
178      def setTolerance(self,rtol=1e-4):      def setTolerance(self,rtol=1e-4):
179          """          """
180          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
181    
182          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
183    
184          where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          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          :param rtol: relative tolerance for the pressure
187          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
188          """          """
189          if rtol<0:          if rtol<0:
190              raise ValueError,"Relative tolerance needs to be non-negative."              raise ValueError,"Relative tolerance needs to be non-negative."
# Line 142  class DarcyFlow(object): Line 193  class DarcyFlow(object):
193          """          """
194          returns the relative tolerance          returns the relative tolerance
195    
196          @return: current relative tolerance          :return: current relative tolerance
197          @rtype: C{float}          :rtype: ``float``
198          """          """
199          return self.__rtol          return self.__rtol
200    
201      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
202          """          """
203          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
204    
205          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
206    
207          where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          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    
209          @param atol: absolute tolerance for the pressure          :param atol: absolute tolerance for the pressure
210          @type atol: non-negative C{float}          :type atol: non-negative ``float``
211          """          """
212          if atol<0:          if atol<0:
213              raise ValueError,"Absolute tolerance needs to be non-negative."              raise ValueError,"Absolute tolerance needs to be non-negative."
# Line 165  class DarcyFlow(object): Line 216  class DarcyFlow(object):
216         """         """
217         returns the absolute tolerance         returns the absolute tolerance
218                
219         @return: current absolute tolerance         :return: current absolute tolerance
220         @rtype: C{float}         :rtype: ``float``
221         """         """
222         return self.__atol         return self.__atol
   
     def setSubProblemTolerance(self,rtol=None):  
          """  
          Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present  
          C{self.getTolerance()**2} is used.  
   
          @param rtol: relative tolerence  
          @type rtol: positive C{float}  
          """  
          if rtol == None:  
               if self.getTolerance()<=0.:  
                   raise ValueError,"A positive relative tolerance must be set."  
               self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)  
          else:  
              if rtol<=0:  
                  raise ValueError,"sub-problem tolerance must be positive."  
              self.__sub_tol=max(util.EPSILON**(0.75),rtol)  
   
223      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
224           """      """
225           Returns the subproblem reduction factor.      Returns a suitable subtolerance
226        @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           @return: subproblem reduction factor      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
          @rtype: C{float}  
          """  
          return self.__sub_tol  
   
     def solve(self,u0,p0, max_iter=100, verbose=False, show_details=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 verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
252           @type verbose: C{bool}           :type verbose: ``bool``
253           @param show_details:  if set information on the subiteration process are printed.           :return: flux and pressure
254           @type show_details: C{bool}           :rtype: ``tuple`` of `Data`.
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
255    
256           @note: The problem is solved as a least squares form           :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 or True           self.verbose=verbose
          self.show_details= show_details and self.verbose  
274           rtol=self.getTolerance()           rtol=self.getTolerance()
275           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
276           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()       self.setSubProblemTolerance()
   
277           num_corrections=0           num_corrections=0
278           converged=False           converged=False
279           p=p0           p=p0
280           norm_r=None           norm_r=None
281           while not converged:           while not converged:
282                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)                 v=self.getFlux(p, fixed_flux=u0)
283                 Qp=self.__Q(p)                 Qp=self.__Q(p)
284                 norm_v=self.__L2(v)                 norm_v=self.__L2(v)
285                 norm_Qp=self.__L2(Qp)                 norm_Qp=self.__L2(Qp)
# Line 259  class DarcyFlow(object): Line 297  class DarcyFlow(object):
297                 if self.verbose:                 if self.verbose:
298                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
299                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
300                        print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
301                        print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
302                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
303                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
304                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
305                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
306                     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.1*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                     num_corrections+=1                     num_corrections+=1
308                 else:                 else:
309                     converged=True                     converged=True
310           return v,p           return v,p
 #  
 #                
 #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp  
 #               #===========================================================================  
 #               norm_r_hat=self.__L2(r_hat)  
 #               norm_v=self.__L2(v)  
 #               norm_g=self.__L2(g)  
 #               norm_gv=self.__L2(g-v)  
 #               norm_Qp=self.__L2(Qp)  
 #               norm_gQp=self.__L2(g-Qp)  
 #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))  
 #               fac=min(norm_v,norm_Qp,norm_gv)  
 #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))  
 #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat  
 #               if r!=None:  
 #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat  
 #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()  
 #                   self.setSubProblemTolerance(sub_tol)  
 #                   print "subtol_new=",self.getSubProblemTolerance()  
 #               print "norm_v = ",norm_v  
 #               print "norm_gv = ",norm_gv  
 #               print "norm_Qp = ",norm_Qp  
 #               print "norm_gQp = ",norm_gQp  
 #               print "norm_g = ",norm_g  
 #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)  
 #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)  
 #               if fac == 0:  
 #                   if self.verbose: print "DarcyFlux: trivial case!"  
 #                   return v,p  
 #               #===============================================================================  
 #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))  
 #               # norm_Qp=self.__L2(Qp)  
 #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))  
 #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp  
 #  
 #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))  
 #               if self.verbose:  
 #                   print "DarcyFlux: residual = %e"%norm_r_hat  
 #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
 #               if norm_r_hat <= ATOL:  
 #                   print "DarcyFlux: iteration finalized."  
 #                   converged=True  
 #               else:  
 #                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)  
 #                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)  
 #                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
 #               print "norm_r =",norm_r  
 #         return v,p  
311      def __L2(self,v):      def __L2(self,v):
312           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
313    
# Line 323  class DarcyFlow(object): Line 315  class DarcyFlow(object):
315            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
316    
317      def __Aprod(self,dp):      def __Aprod(self,dp):
318            self.__pde_v.setTolerance(self.getSubProblemTolerance())            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
           if self.show_details: print "DarcyFlux: Applying operator"  
319            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
320            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
321            du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)            du=self.__pde_v.getSolution()
322              # self.__pde_v.getOperator().saveMM("proj.mm")
323            return Qdp+du            return Qdp+du
324      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
325           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 336  class DarcyFlow(object): Line 328  class DarcyFlow(object):
328           return util.integrate(util.inner(self.__Q(p), r))           return util.integrate(util.inner(self.__Q(p), r))
329    
330      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
331            self.__pde_p.setTolerance(self.getSubProblemTolerance())        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
           if self.show_details: print "DarcyFlux: Applying preconditioner"  
332            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), 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, iter_max = 100000)            # self.__pde_p.getOperator().saveMM("prec.mm")
334              return self.__pde_p.getSolution()
335    
336      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):      def getFlux(self,p=None, fixed_flux=Data()):
337          """          """
338          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
339          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
340          Note that C{g} and C{f} are used, see L{setValue}.          Note that ``g`` and ``f`` are used, see `setValue`.
341    
342          @param p: pressure.          :param p: pressure.
343          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `Data`).
344          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          :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. L{Data}).          :type fixed_flux: vector values on the domain (e.g. `Data`).
346          @param tol: relative tolerance to be used.          :return: flux
347          @type tol: positive C{float}.          :rtype: `Data`
348          @return: flux          :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          @rtype: L{Data}                 for the permeability *k_{ij}*
         @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}}  
350          """          """
351          self.__pde_v.setTolerance(self.getSubProblemTolerance())      self.setSubProblemTolerance()
352          g=self.__g          g=self.__g
353          f=self.__f          f=self.__f
354          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
# Line 366  class DarcyFlow(object): Line 356  class DarcyFlow(object):
356             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
357          else:          else:
358             self.__pde_v.setValue(Y=g-self.__Q(p))             self.__pde_v.setValue(Y=g-self.__Q(p))
359          return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)          return self.__pde_v.getSolution()
360    
361  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
362       """       """
# Line 391  class StokesProblemCartesian(Homogeneous Line 381  class StokesProblemCartesian(Homogeneous
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       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):           self.__pde_proj=LinearPDE(domain)
401             self.__pde_proj.setReducedOrderOn()
402         self.__pde_proj.setValue(D=1)
403             self.__pde_proj.setSymmetryOn()
404    
405         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
502          @note: All values needs to be set.          """
503            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         """  
         self.eta=eta  
         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)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
504    
505       def inner_pBv(self,p,v):       def Bv(self,v,tol):
506           """           """
507           returns inner product of element p and div(v)           returns inner product of element p and div(v)
508    
509           @param p: a pressure increment           :param v: a residual
510           @param v: a residual           :return: inner product of element p and div(v)
511           @return: inner product of element p and div(v)           :rtype: ``float``
512           @rtype: C{float}           """
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             :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           return util.integrate(-p*util.div(v))           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
529    
530       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
531           """           """
532           Returns inner product of p0 and p1           Returns inner product of p0 and p1
533    
534           @param p0: a pressure           :param p0: a pressure
535           @param p1: a pressure           :param p1: a pressure
536           @return: inner product of p0 and p1           :return: inner product of p0 and p1
537           @rtype: C{float}           :rtype: ``float``
538           """           """
539           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
540           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
541           return util.integrate(s0*s1)           return util.integrate(s0*s1)
542    
543       def norm_v(self,v):       def norm_v(self,v):
544           """           """
545           returns the norm of v           returns the norm of v
546    
547           @param v: a velovity           :param v: a velovity
548           @return: norm of v           :return: norm of v
549           @rtype: non-negative C{float}           :rtype: non-negative ``float``
550           """           """
551           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
552    
553       def getV(self, p, v0):  
554         def getDV(self, p, v, tol):
555           """           """
556           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
557    
558           @param p: a pressure           :param p: a pressure
559           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
560           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
561           """           """
562           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.updateStokesEquation(v,p)
563           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           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=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+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 norm_Bv(self,Bv):
          raise NotImplementedError,"no v calculation implemented."  
   
   
      def norm_Bv(self,v):  
574          """          """
575          Returns Bv (overwrite).          Returns Bv (overwrite).
576    
577          @rtype: equal to the type of p          :rtype: equal to the type of p
578          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
579          """          """
580          return util.sqrt(util.integrate(util.div(v)**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
581    
582       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
583           """           """
584           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
585    
586           @param p: a pressure increment           :param p: a pressure increment
587           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
588           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
589           """           """
590           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
591           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           out=self.__pde_u.getSolution()
          out=self.__pde_u.getSolution(verbose=self.show_details)  
592           return  out           return  out
593    
594       def solve_precB(self,v):       def solve_prec(self,Bv, tol):
595           """           """
596           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
597           with accuracy L{self.getSubProblemTolerance()} (overwrite).           with accuracy `self.getSubProblemTolerance()`
598    
599           @param v: velocity increment           :param Bv: velocity increment
600           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
601           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
602           """           """
603           self.__pde_prec.setValue(Y=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
604           self.__pde_prec.setTolerance(self.getSubProblemTolerance())       self.getSolverOptionsPressure().setTolerance(tol)
605           return self.__pde_prec.getSolution(verbose=self.show_details)       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
606             out=self.__pde_prec.getSolution()
607             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26