/[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 2719 by gross, Wed Oct 14 06:38:03 2009 UTC
# Line 1  Line 1 
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
4  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
5  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
6  #  #
# Line 10  Line 10 
10  #  #
11  ########################################################  ########################################################
12    
13  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
15  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
16  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
17  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
18  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
19  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
20    
21  """  """
22  Some models for flow  Some models for flow
23    
24  @var __author__: name of author  :var __author__: name of author
25  @var __copyright__: copyrights  :var __copyright__: copyrights
26  @var __license__: licence agreement  :var __license__: licence agreement
27  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
28  @var __version__: version  :var __version__: version
29  @var __date__: date of the version  :var __date__: date of the version
30  """  """
31    
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
37  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38    
39  class DarcyFlow(object):  class DarcyFlow(object):
40      """      """
41      solves the problem      solves the problem
42    
43      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
44      M{u_{i,i} = f}      *u_{i,i} = f*
45    
46      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
47    
48      @note: The problem is solved in a least squares formulation.      :note: The problem is solved in a least squares formulation.
49      """      """
50    
51      def __init__(self, domain,useReduced=False):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
52          """          """
53          initializes the Darcy flux problem          initializes the Darcy flux problem
54          @param domain: domain of the problem          :param domain: domain of the problem
55          @type domain: L{Domain}          :type domain: `Domain`
56        :param useReduced: uses reduced oreder on flux and pressure
57        :type useReduced: ``bool``
58        :param adaptSubTolerance: switches on automatic subtolerance selection
59        :type adaptSubTolerance: ``bool``  
60          """          """
61          self.domain=domain          self.domain=domain
62            if weight == None:
63               s=self.domain.getSize()
64               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65               # self.__l=(3.*util.longestEdge(self.domain))**2
66               #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
67            else:
68               self.__l=weight
69          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
70          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
71          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
72          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
73          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
74          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
75          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
76          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
77          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
78          self.__ATOL= None          self.setTolerance()
79            self.setAbsoluteTolerance()
80        self.__adaptSubTolerance=adaptSubTolerance
81        self.verbose=False
82        def getSolverOptionsFlux(self):
83        """
84        Returns the solver options used to solve the flux problems
85        
86        *(I+D^*D)u=F*
87        
88        :return: `SolverOptions`
89        """
90        return self.__pde_v.getSolverOptions()
91        def setSolverOptionsFlux(self, options=None):
92        """
93        Sets the solver options used to solve the flux problems
94        
95        *(I+D^*D)u=F*
96        
97        If ``options`` is not present, the options are reset to default
98        :param options: `SolverOptions`
99        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
100        """
101        return self.__pde_v.setSolverOptions(options)
102        def getSolverOptionsPressure(self):
103        """
104        Returns the solver options used to solve the pressure problems
105        
106        *(Q^*Q)p=Q^*G*
107        
108        :return: `SolverOptions`
109        """
110        return self.__pde_p.getSolverOptions()
111        def setSolverOptionsPressure(self, options=None):
112        """
113        Sets the solver options used to solve the pressure problems
114        
115        *(Q^*Q)p=Q^*G*
116        
117        If ``options`` is not present, the options are reset to default
118        :param options: `SolverOptions`
119        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
120        """
121        return self.__pde_p.setSolverOptions(options)
122    
123      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124          """          """
125          assigns values to model parameters          assigns values to model parameters
126    
127          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
128          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `Data`)
129          @param g: flux sources/sinks          :param g: flux sources/sinks
130          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `Data`)
131          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
132          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
133          @param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
134          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
135          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
136                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
137                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
138          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
139    
140          @note: the values of parameters which are not set by calling C{setValue} are not altered.          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
141          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
142                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
143                 is along the M{x_i} axis.                 is along the *x_i* axis.
144          """          """
145          if f !=None:          if f !=None:
146             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
147             if f.isEmpty():             if f.isEmpty():
148                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
149             else:             else:
150                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
151             self.f=f             self.__f=f
152          if g !=None:            if g !=None:
153             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
154             if g.isEmpty():             if g.isEmpty():
155               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 121  class DarcyFlow(object): Line 175  class DarcyFlow(object):
175             self.__permeability=perm             self.__permeability=perm
176             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
177    
178        def setTolerance(self,rtol=1e-4):
     def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):  
179          """          """
180          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
         on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
         Note that C{g} and C{f} are used, see L{setValue}.  
           
         @param p: pressure.  
         @type p: scalar value on the domain (e.g. L{Data}).  
         @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.  
         @type fixed_flux: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: flux  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
         """  
         self.__pde_v.setTolerance(tol)  
         g=self.__g  
         f=self.__f  
         self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
         else:  
            self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))  
         return self.__pde_v.getSolution(verbose=show_details)  
181    
182      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
183    
184            where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
185    
186            :param rtol: relative tolerance for the pressure
187            :type rtol: non-negative ``float``
188          """          """
189          returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
190          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
191          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
192                def getTolerance(self):
         @param v: flux.  
         @type v: vector-valued on the domain (e.g. L{Data}).  
         @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.  
         @type fixed_pressure: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: pressure  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
193          """          """
194          self.__pde_v.setTolerance(tol)          returns the relative tolerance
         g=self.__g  
         self.__pde_p.setValue(r=fixed_pressure)  
         if v == None:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))  
         else:  
            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))  
         return self.__pde_p.getSolution(verbose=show_details)  
195    
196      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          :return: current relative tolerance
197            :rtype: ``float``
198          """          """
199          set the tolerance C{ATOL} used to terminate the solution process. It is used          return self.__rtol
   
         M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}  
   
         where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.  
200    
201          The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual      def setAbsoluteTolerance(self,atol=0.):
   
         M{r=Q^*(g-Qp-v)}  
   
         the condition  
   
         M{<(Q^*Q)^{-1} r,r> <= ATOL}  
   
         holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}  
   
         @param atol: absolute tolerance for the pressure  
         @type atol: non-negative C{float}  
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{float}  
         @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the  
         L{getPressure} method or use  the value from a previous time step.  
         @type p_ref: scalar value on the domain (e.g. L{Data}).  
         @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the  
         L{getFlux} method or use  the value from a previous time step.  
         @type v_ref: vector-valued on the domain (e.g. L{Data}).  
         @return: used absolute tolerance.  
         @rtype: positive C{float}  
         """  
         g=self.__g  
         if not v_ref == None:  
            f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)  
         else:  
            f1=util.integrate(util.length(util.interpolate(g))**2)  
         if not p_ref == None:  
            f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)  
         else:  
            f2=0  
         self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))  
         if self.__ATOL<=0:  
            raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL  
         return self.__ATOL  
           
     def getTolerance(self):  
202          """          """
203          returns the current tolerance.          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
204      
205          @return: used absolute tolerance.          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
206          @rtype: positive C{float}  
207          """          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 self.__ATOL==None:  
209             raise ValueError,"no tolerance is defined."          :param atol: absolute tolerance for the pressure
210          return self.__ATOL          :type atol: non-negative ``float``
211            """
212            if atol<0:
213                raise ValueError,"Absolute tolerance needs to be non-negative."
214            self.__atol=atol
215        def getAbsoluteTolerance(self):
216           """
217           returns the absolute tolerance
218          
219           :return: current absolute tolerance
220           :rtype: ``float``
221           """
222           return self.__atol
223        def getSubProblemTolerance(self):
224        """
225        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      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.           :param domain: domain of the problem. The approximation order needs to be two.
385           @type domain: L{Domain}           :type domain: `Domain`
386           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           :warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
387           """           """
388           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
389           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
390           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())
391           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
392           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
393           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
394           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
395           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
396    
397           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
398           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
399         self.__pde_proj.setValue(D=1)
400           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
401    
402        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
403             """
404         returns the solver options used  solve the equation for velocity.
405        
406         :rtype: `SolverOptions`
407         """
408         return self.__pde_u.getSolverOptions()
409         def setSolverOptionsVelocity(self, options=None):
410             """
411         set the solver options for solving the equation for velocity.
412        
413         :param options: new solver  options
414         :type options: `SolverOptions`
415         """
416             self.__pde_u.setSolverOptions(options)
417         def getSolverOptionsPressure(self):
418             """
419         returns the solver options used  solve the equation for pressure.
420         :rtype: `SolverOptions`
421         """
422         return self.__pde_prec.getSolverOptions()
423         def setSolverOptionsPressure(self, options=None):
424             """
425         set the solver options for solving the equation for pressure.
426         :param options: new solver  options
427         :type options: `SolverOptions`
428         """
429         self.__pde_prec.setSolverOptions(options)
430    
431         def setSolverOptionsDiv(self, options=None):
432             """
433         set the solver options for solving the equation to project the divergence of
434         the velocity onto the function space of presure.
435        
436         :param options: new solver options
437         :type options: `SolverOptions`
438         """
439         self.__pde_proj.setSolverOptions(options)
440         def getSolverOptionsDiv(self):
441             """
442         returns the solver options for solving the equation to project the divergence of
443         the velocity onto the function space of presure.
444        
445         :rtype: `SolverOptions`
446         """
447         return self.__pde_proj.getSolverOptions()
448    
449         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
450          """          """
451          assigns values to the model parameters          assigns values to the model parameters
452    
453          @param f: external force          :param f: external force
454          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
455          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
456          @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
457          @param eta: viscosity          :param eta: viscosity
458          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
459          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
460          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
461          @param stress: initial stress          :param stress: initial stress
462      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
463          @note: All values needs to be set.          :note: All values needs to be set.
   
464          """          """
465          self.eta=eta          self.eta=eta
466          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
467      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
468          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
469          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
470              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
471              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
472      self.__pde_prec.setValue(D=1/self.eta)          n=self.domain.getNormal()
473          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)      self.__pde_prec.setValue(D=1/self.eta)
474            self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))
475            self.__f=f
476            self.__surface_stress=surface_stress
477          self.__stress=stress          self.__stress=stress
478    
479        def B(self,v):       def Bv(self,v,tol):
480          """           """
481          returns div(v)           returns inner product of element p and div(v)
         @rtype: equal to the type of p  
482    
483          @note: boundary conditions on p should be zero!           :param v: a residual
484          """           :return: inner product of element p and div(v)
485          if self.show_details: print "apply divergence:"           :rtype: ``float``
486          self.__pde_proj.setValue(Y=-util.div(v))           """
487          self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_proj.setValue(Y=-util.div(v)) # -???
488          return self.__pde_proj.getSolution(verbose=self.show_details)       self.getSolverOptionsDiv().setTolerance(tol)
489         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
490             out=self.__pde_proj.getSolution()
491             return out
492    
493        def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
494           """           """
495           returns inner product of element p and Bv  (overwrite)           returns inner product of element p and Bv=-div(v)
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @rtype: C{float}  
496    
497           @rtype: equal to the type of p           :param p: a pressure increment
498             :param Bv: a residual
499             :return: inner product of element p and Bv=-div(v)
500             :rtype: ``float``
501           """           """
502           s0=util.interpolate(p,Function(self.domain))           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
          s1=util.interpolate(Bv,Function(self.domain))  
          return util.integrate(s0*s1)  
503    
504        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
505           """           """
506           returns inner product of element p0 and p1  (overwrite)           Returns inner product of p0 and p1
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
507    
508           @rtype: equal to the type of p           :param p0: a pressure
509             :param p1: a pressure
510             :return: inner product of p0 and p1
511             :rtype: ``float``
512           """           """
513           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
514           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
515           return util.integrate(s0*s1)           return util.integrate(s0*s1)
516    
517        def inner_v(self,v0,v1):       def norm_v(self,v):
518           """           """
519           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}  
520    
521           @rtype: equal to the type of v           :param v: a velovity
522             :return: norm of v
523             :rtype: non-negative ``float``
524           """           """
525       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))  
526    
527        def solve_A(self,u,p):       def getDV(self, p, v, tol):
528           """           """
529           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
530    
531             :param p: a pressure
532             :param v: a initial guess for the value v to return.
533             :return: dv given as *Adv=(f-Av-B^*p)*
534           """           """
535           if self.show_details: print "solve for velocity:"           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
536           self.__pde_u.setTolerance(self.getSubProblemTolerance())       self.getSolverOptionsVelocity().setTolerance(tol)
537         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
538           if self.__stress.isEmpty():           if self.__stress.isEmpty():
539              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)))
540           else:           else:
541              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)))
542           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
543             return  out
544    
545         def norm_Bv(self,Bv):
546            """
547            Returns Bv (overwrite).
548    
549            :rtype: equal to the type of p
550            :note: boundary conditions on p should be zero!
551            """
552            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
553    
554         def solve_AinvBt(self,p, tol):
555             """
556             Solves *Av=B^*p* with accuracy `tol`
557    
558             :param p: a pressure increment
559             :return: the solution of *Av=B^*p*
560             :note: boundary conditions on v should be zero!
561             """
562             self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
563             out=self.__pde_u.getSolution()
564           return  out           return  out
565    
566        def solve_prec(self,p):       def solve_prec(self,Bv, tol):
567           if self.show_details: print "apply preconditioner:"           """
568           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           applies preconditioner for for *BA^{-1}B^** to *Bv*
569           self.__pde_prec.setValue(Y=p)           with accuracy `self.getSubProblemTolerance()`
570           q=self.__pde_prec.getSolution(verbose=self.show_details)  
571           return q           :param Bv: velocity increment
572             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
573             :note: boundary conditions on p are zero.
574             """
575             self.__pde_prec.setValue(Y=Bv)
576         self.getSolverOptionsPressure().setTolerance(tol)
577         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
578             out=self.__pde_prec.getSolution()
579             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26