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

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

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

revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 2354 by gross, Wed Apr 1 04:01:58 2009 UTC
# Line 16  http://www.uq.edu.au/esscc Line 16  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
# Line 34  __author__="Lutz Gross, l.gross@uq.edu.a Line 34  __author__="Lutz Gross, l.gross@uq.edu.a
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
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}      M{u_i+k_{ij}*p_{,j} = g_i}
44      M{u_{i,i} = f}      M{u_{i,i} = f}
45    
46      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      where M{p} represents the pressure and M{u} the Darcy flux. M{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):
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: L{Domain}
56          """          """
57          self.domain=domain          self.domain=domain
58            if weight == None:
59               s=self.domain.getSize()
60               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
61            else:
62               self.__l=weight
63          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
64          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
65          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
66          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)))
67            # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
68          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
69          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
70          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
71          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
73          self.__ATOL= None          self.setTolerance()
74            self.setAbsoluteTolerance()
75            self.setSubProblemTolerance()
76    
77      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):
78          """          """
# Line 78  class DarcyFlow(object): Line 86  class DarcyFlow(object):
86          @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. L{Data})
87          @param location_of_fixed_flux:  mask for locations where flux is fixed.          @param location_of_fixed_flux:  mask for locations where flux is fixed.
88          @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. L{Data})
89          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          @param permeability: permeability tensor. If scalar C{s} is given the tensor with
90                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with
91                               C{v} on the main diagonal is used.                               C{v} on the main diagonal is used.
92          @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. L{Data})
93    
# Line 88  class DarcyFlow(object): Line 96  class DarcyFlow(object):
96                 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 (C{location_of_fixed_flux[i]>0} if direction of the normal
97                 is along the M{x_i} axis.                 is along the M{x_i} axis.
98          """          """
99          if f !=None:          if f !=None:
100             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101             if f.isEmpty():             if f.isEmpty():
102                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103             else:             else:
104                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
105             self.f=f             self.__f=f
106          if g !=None:            if g !=None:
107             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108             if g.isEmpty():             if g.isEmpty():
109               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 121  class DarcyFlow(object): Line 129  class DarcyFlow(object):
129             self.__permeability=perm             self.__permeability=perm
130             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))
131    
132        def setTolerance(self,rtol=1e-4):
     def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):  
         """  
         returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}  
         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}}  
133          """          """
134          self.__pde_v.setTolerance(tol)          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
135          g=self.__g  
136          f=self.__f          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
137          self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)  
138          if p == None:          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}}.
            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)  
139    
140      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          @param rtol: relative tolerance for the pressure
141            @type rtol: non-negative C{float}
142          """          """
143          returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
144          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
145          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
146                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}}  
147          """          """
148          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)  
149    
150      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          @return: current relative tolerance
151            @rtype: C{float}
152          """          """
153          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.  
154    
155          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.):
156            """
157          M{r=Q^*(g-Qp-v)}          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
   
         the condition  
158    
159          M{<(Q^*Q)^{-1} r,r> <= ATOL}          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
160    
161          holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}          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}}.
162    
163          @param atol: absolute tolerance for the pressure          @param atol: absolute tolerance for the pressure
164          @type atol: non-negative C{float}          @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):  
165          """          """
166          returns the current tolerance.          if atol<0:
167                  raise ValueError,"Absolute tolerance needs to be non-negative."
168          @return: used absolute tolerance.          self.__atol=atol
169          @rtype: positive C{float}      def getAbsoluteTolerance(self):
170          """         """
171          if self.__ATOL==None:         returns the absolute tolerance
172             raise ValueError,"no tolerance is defined."        
173          return self.__ATOL         @return: current absolute tolerance
174           @rtype: C{float}
175           """
176           return self.__atol
177    
178        def setSubProblemTolerance(self,rtol=None):
179             """
180             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
181             C{self.getTolerance()**2} is used.
182    
183             @param rtol: relative tolerence
184             @type rtol: positive C{float}
185             """
186             if rtol == None:
187                  if self.getTolerance()<=0.:
188                      raise ValueError,"A positive relative tolerance must be set."
189                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
190             else:
191                 if rtol<=0:
192                     raise ValueError,"sub-problem tolerance must be positive."
193                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
194    
195        def getSubProblemTolerance(self):
196             """
197             Returns the subproblem reduction factor.
198    
199      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):           @return: subproblem reduction factor
200           """           @rtype: C{float}
201             """
202             return self.__sub_tol
203    
204        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
205             """
206           solves the problem.           solves the problem.
207    
208           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().
209    
210           @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 C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
211           @type u0: vector value on the domain (e.g. L{Data}).           @type u0: vector value on the domain (e.g. L{Data}).
212           @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 C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
213           @type p0: scalar value on the domain (e.g. L{Data}).           @type p0: scalar value on the domain (e.g. L{Data}).
          @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}  
          @type sub_rtol: positive-negative C{float}  
214           @param verbose: if set some information on iteration progress are printed           @param verbose: if set some information on iteration progress are printed
215           @type verbose: C{bool}           @type verbose: C{bool}
216           @param show_details:  if set information on the subiteration process are printed.           @param show_details:  if set information on the subiteration process are printed.
217           @type show_details: C{bool}           @type show_details: C{bool}
218           @return: flux and pressure           @return: flux and pressure
219           @rtype: C{tuple} of L{Data}.           @rtype: C{tuple} of L{Data}.
220    
221           @note: The problem is solved as a least squares form           @note: The problem is solved as a least squares form
222    
223           M{(I+D^*D)u+Qp=D^*f+g}           M{(I+D^*D)u+Qp=D^*f+g}
224           M{Q^*u+Q^*Qp=Q^*g}           M{Q^*u+Q^*Qp=Q^*g}
225    
226           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
227           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
228    
229           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
230    
231           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
232    
233           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           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
234            
235           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 M{Q^*Q}). In each iteration step
236           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
237           """           """
238           self.verbose=verbose           self.verbose=verbose or True
239           self.show_details= show_details and self.verbose           self.show_details= show_details and self.verbose
240           self.__pde_v.setTolerance(sub_rtol)           rtol=self.getTolerance()
241           self.__pde_p.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
242           ATOL=self.getTolerance()           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL  
244           #########################################################################################################################           num_corrections=0
245           #           converged=False
246           #   we solve:           p=p0
247           #             norm_r=None
248           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )           while not converged:
249           #                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250           #   residual is                 Qp=self.__Q(p)
251           #                 norm_v=self.__L2(v)
252           #    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)                 norm_Qp=self.__L2(Qp)
253           #                 if norm_v == 0.:
254           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                    if norm_Qp == 0.:
255           #                       return v,p
256           #    we use (g - Qp, v) to represent the residual. not that                    else:
257           #                      fac=norm_Qp
258           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                 else:
259           #                    if norm_Qp == 0.:
260           #   while the initial residual is                      fac=norm_v
261           #                    else:
262           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                      fac=2./(1./norm_v+1./norm_Qp)
263           #                   ATOL=(atol+rtol*fac)
264           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                 if self.verbose:
265           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 norm of v = %e."%norm_v
266           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
267           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
268           self.__pde_v.setValue(r=Data())                 if norm_r == None or norm_r>ATOL:
269           # start CG                     if num_corrections>max_num_corrections:
270           r=ArithmeticTuple(d0, v00)                           raise ValueError,"maximum number of correction steps reached."
271           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)
272           return r[1],p                     num_corrections+=1
273                   else:
274                       converged=True
275             return v,p
276    #
277    #              
278    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
279    #               #===========================================================================
280    #               norm_r_hat=self.__L2(r_hat)
281    #               norm_v=self.__L2(v)
282    #               norm_g=self.__L2(g)
283    #               norm_gv=self.__L2(g-v)
284    #               norm_Qp=self.__L2(Qp)
285    #               norm_gQp=self.__L2(g-Qp)
286    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
287    #               fac=min(norm_v,norm_Qp,norm_gv)
288    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
289    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
290    #               if r!=None:
291    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
292    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
293    #                   self.setSubProblemTolerance(sub_tol)
294    #                   print "subtol_new=",self.getSubProblemTolerance()
295    #               print "norm_v = ",norm_v
296    #               print "norm_gv = ",norm_gv
297    #               print "norm_Qp = ",norm_Qp
298    #               print "norm_gQp = ",norm_gQp
299    #               print "norm_g = ",norm_g
300    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
301    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
302    #               if fac == 0:
303    #                   if self.verbose: print "DarcyFlux: trivial case!"
304    #                   return v,p
305    #               #===============================================================================
306    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
307    #               # norm_Qp=self.__L2(Qp)
308    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
309    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
310    #
311    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
312    #               if self.verbose:
313    #                   print "DarcyFlux: residual = %e"%norm_r_hat
314    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
315    #               if norm_r_hat <= ATOL:
316    #                   print "DarcyFlux: iteration finalized."
317    #                   converged=True
318    #               else:
319    #                   # 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)
320    #                   # 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)
321    #                   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)
322    #               print "norm_r =",norm_r
323    #         return v,p
324        def __L2(self,v):
325             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
326    
327      def __Aprod_PCG(self,dp):      def __Q(self,p):
328            if self.show_details: print "DarcyFlux: Applying operator"            return util.tensor_mult(self.__permeability,util.grad(p))
           #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)  
           mQdp=util.tensor_mult(self.__permeability,util.grad(dp))  
           self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())  
           du=self.__pde_v.getSolution(verbose=self.show_details)  
           return ArithmeticTuple(mQdp,du)  
329    
330        def __Aprod(self,dp):
331              self.__pde_v.setTolerance(self.getSubProblemTolerance())
332              if self.show_details: print "DarcyFlux: Applying operator"
333              Qdp=self.__Q(dp)
334              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
335              du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
336              return Qdp+du
337        def __inner_GMRES(self,r,s):
338             return util.integrate(util.inner(r,s))
339    
340      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
341           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  
342    
343      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
344              self.__pde_p.setTolerance(self.getSubProblemTolerance())
345            if self.show_details: print "DarcyFlux: Applying preconditioner"            if self.show_details: print "DarcyFlux: Applying preconditioner"
346            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())
347            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
348    
349        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
350            """
351            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
352            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
353            Note that C{g} and C{f} are used, see L{setValue}.
354    
355            @param p: pressure.
356            @type p: scalar value on the domain (e.g. L{Data}).
357            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
358            @type fixed_flux: vector values on the domain (e.g. L{Data}).
359            @param tol: relative tolerance to be used.
360            @type tol: positive C{float}.
361            @return: flux
362            @rtype: L{Data}
363            @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}}
364                   for the permeability M{k_{ij}}
365            """
366            self.__pde_v.setTolerance(self.getSubProblemTolerance())
367            g=self.__g
368            f=self.__f
369            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
370            if p == None:
371               self.__pde_v.setValue(Y=g)
372            else:
373               self.__pde_v.setValue(Y=g-self.__Q(p))
374            return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
375    
376  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
377        """       """
378        solves       solves
379    
380            -(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}
381                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 383  class StokesProblemCartesian(Homogeneous
383            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
384            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
385    
386        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
387    
388        typical usage:       typical usage:
389    
390              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
391              sp.setTolerance()              sp.setTolerance()
392              sp.initialize(...)              sp.initialize(...)
393              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
394        """       """
395        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
396           """           """
397           initialize the Stokes Problem           initialize the Stokes Problem
398    
# Line 356  class StokesProblemCartesian(Homogeneous Line 406  class StokesProblemCartesian(Homogeneous
406           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())
407           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
408           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
409           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
410                
411           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
412           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
413           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
414           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
415    
416           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
          self.__pde_proj.setReducedOrderOn()  
          self.__pde_proj.setSymmetryOn()  
          self.__pde_proj.setValue(D=1.)  
   
       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):  
417          """          """
418          assigns values to the model parameters          assigns values to the model parameters
419    
# Line 383  class StokesProblemCartesian(Homogeneous Line 428  class StokesProblemCartesian(Homogeneous
428          @param stress: initial stress          @param stress: initial stress
429      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
430          @note: All values needs to be set.          @note: All values needs to be set.
431    
432          """          """
433          self.eta=eta          self.eta=eta
434          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
435      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
436          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
437          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
438              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
439              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
440      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
441          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
442            self.__f=f
443            self.__surface_stress=surface_stress
444          self.__stress=stress          self.__stress=stress
445    
446        def B(self,v):       def inner_pBv(self,p,v):
447          """           """
448          returns div(v)           returns inner product of element p and div(v)
         @rtype: equal to the type of p  
449    
450          @note: boundary conditions on p should be zero!           @param p: a pressure increment
451          """           @param v: a residual
452          if self.show_details: print "apply divergence:"           @return: inner product of element p and div(v)
         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  
453           @rtype: C{float}           @rtype: C{float}
   
          @rtype: equal to the type of p  
454           """           """
455           s0=util.interpolate(p,Function(self.domain))           return util.integrate(-p*util.div(v))
          s1=util.interpolate(Bv,Function(self.domain))  
          return util.integrate(s0*s1)  
456    
457        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
458           """           """
459           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}  
460    
461           @rtype: equal to the type of p           @param p0: a pressure
462             @param p1: a pressure
463             @return: inner product of p0 and p1
464             @rtype: C{float}
465           """           """
466           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
467           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
468           return util.integrate(s0*s1)           return util.integrate(s0*s1)
469    
470        def inner_v(self,v0,v1):       def norm_v(self,v):
471           """           """
472           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}  
473    
474           @rtype: equal to the type of v           @param v: a velovity
475             @return: norm of v
476             @rtype: non-negative C{float}
477           """           """
478       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
479    
480        def solve_A(self,u,p):       def getV(self, p, v0):
481           """           """
482           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
483    
484             @param p: a pressure
485             @param v0: a initial guess for the value v to return.
486             @return: v given as M{v= A^{-1} (f-B^*p)}
487           """           """
          if self.show_details: print "solve for velocity:"  
488           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
489             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
490           if self.__stress.isEmpty():           if self.__stress.isEmpty():
491              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))
492           else:           else:
493              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))
494           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
495           return  out           return  out
496    
497        def solve_prec(self,p):  
498           if self.show_details: print "apply preconditioner:"           raise NotImplementedError,"no v calculation implemented."
499    
500    
501         def norm_Bv(self,v):
502            """
503            Returns Bv (overwrite).
504    
505            @rtype: equal to the type of p
506            @note: boundary conditions on p should be zero!
507            """
508            return util.sqrt(util.integrate(util.div(v)**2))
509    
510         def solve_AinvBt(self,p):
511             """
512             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
513    
514             @param p: a pressure increment
515             @return: the solution of M{Av=B^*p}
516             @note: boundary conditions on v should be zero!
517             """
518             self.__pde_u.setTolerance(self.getSubProblemTolerance())
519             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
520             out=self.__pde_u.getSolution(verbose=self.show_details)
521             return  out
522    
523         def solve_precB(self,v):
524             """
525             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
526             with accuracy L{self.getSubProblemTolerance()} (overwrite).
527    
528             @param v: velocity increment
529             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
530             @note: boundary conditions on p are zero.
531             """
532             self.__pde_prec.setValue(Y=-util.div(v))
533           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
534           self.__pde_prec.setValue(Y=p)           return self.__pde_prec.getSolution(verbose=self.show_details)
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  

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

  ViewVC Help
Powered by ViewVC 1.1.26