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

  ViewVC Help
Powered by ViewVC 1.1.26