/[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 2251 by gross, Fri Feb 6 06:50:39 2009 UTC revision 2288 by gross, Tue Feb 24 06:11:48 2009 UTC
# 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      """      """
# Line 55  class DarcyFlow(object): Line 55  class DarcyFlow(object):
55          @type domain: L{Domain}          @type domain: L{Domain}
56          """          """
57          self.domain=domain          self.domain=domain
58            self.__l=util.longestEdge(self.domain)**2
59          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
60          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
61          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
62          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)))
63          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
64          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
65          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
66          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
67          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
68          self.__ATOL= None          self.setTolerance()
69            self.setAbsoluteTolerance()
70            self.setSubProblemTolerance()
71    
72      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):
73          """          """
# Line 78  class DarcyFlow(object): Line 81  class DarcyFlow(object):
81          @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})
82          @param location_of_fixed_flux:  mask for locations where flux is fixed.          @param location_of_fixed_flux:  mask for locations where flux is fixed.
83          @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})
84          @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
85                               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
86                               C{v} on the main diagonal is used.                               C{v} on the main diagonal is used.
87          @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})
88    
# Line 88  class DarcyFlow(object): Line 91  class DarcyFlow(object):
91                 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
92                 is along the M{x_i} axis.                 is along the M{x_i} axis.
93          """          """
94          if f !=None:          if f !=None:
95             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
96             if f.isEmpty():             if f.isEmpty():
97                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
98             else:             else:
99                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
100             self.f=f             self.__f=f
101          if g !=None:            if g !=None:
102             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
103             if g.isEmpty():             if g.isEmpty():
104               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 121  class DarcyFlow(object): Line 124  class DarcyFlow(object):
124             self.__permeability=perm             self.__permeability=perm
125             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))
126    
127        def setTolerance(self,rtol=1e-4):
     def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):  
128          """          """
129          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
130          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
131          Note that C{g} and C{f} are used, see L{setValue}.          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
132            
133          @param p: pressure.          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}}.
         @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)  
134    
135      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          @param rtol: relative tolerance for the pressure
136            @type rtol: non-negative C{float}
137          """          """
138          returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
139          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
140          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
141                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}}  
142          """          """
143          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)  
144    
145      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          @return: current relative tolerance
146            @rtype: C{float}
147          """          """
148          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.  
149    
150          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.):
151            """
152          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  
153    
154          M{<(Q^*Q)^{-1} r,r> <= ATOL}          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
155    
156          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}}.
157    
158          @param atol: absolute tolerance for the pressure          @param atol: absolute tolerance for the pressure
159          @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):  
160          """          """
161          returns the current tolerance.          if atol<0:
162                  raise ValueError,"Absolute tolerance needs to be non-negative."
163          @return: used absolute tolerance.          self.__atol=atol
164          @rtype: positive C{float}      def getAbsoluteTolerance(self):
165          """         """
166          if self.__ATOL==None:         returns the absolute tolerance
167             raise ValueError,"no tolerance is defined."        
168          return self.__ATOL         @return: current absolute tolerance
169           @rtype: C{float}
170           """
171           return self.__atol
172    
173        def setSubProblemTolerance(self,rtol=None):
174             """
175             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
176             C{self.getTolerance()**2} is used.
177    
178             @param rtol: relative tolerence
179             @type rtol: positive C{float}
180             """
181             if rtol == None:
182                  if self.getTolerance()<=0.:
183                      raise ValueError,"A positive relative tolerance must be set."
184                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
185             else:
186                 if rtol<=0:
187                     raise ValueError,"sub-problem tolerance must be positive."
188                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
189    
190        def getSubProblemTolerance(self):
191             """
192             Returns the subproblem reduction factor.
193    
194             @return: subproblem reduction factor
195             @rtype: C{float}
196             """
197             return self.__sub_tol
198    
199      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, show_details=False, max_num_corrections=10):
200           """           """
201           solves the problem.           solves the problem.
202    
203           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().
204    
205           @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.
206           @type u0: vector value on the domain (e.g. L{Data}).           @type u0: vector value on the domain (e.g. L{Data}).
207           @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.
208           @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}  
209           @param verbose: if set some information on iteration progress are printed           @param verbose: if set some information on iteration progress are printed
210           @type verbose: C{bool}           @type verbose: C{bool}
211           @param show_details:  if set information on the subiteration process are printed.           @param show_details:  if set information on the subiteration process are printed.
212           @type show_details: C{bool}           @type show_details: C{bool}
213           @return: flux and pressure           @return: flux and pressure
214           @rtype: C{tuple} of L{Data}.           @rtype: C{tuple} of L{Data}.
215    
216           @note: The problem is solved as a least squares form           @note: The problem is solved as a least squares form
217    
218           M{(I+D^*D)u+Qp=D^*f+g}           M{(I+D^*D)u+Qp=D^*f+g}
219           M{Q^*u+Q^*Qp=Q^*g}           M{Q^*u+Q^*Qp=Q^*g}
220    
221           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}}.
222           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
223    
224           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
225    
226           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
227    
228           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
229            
230           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
231           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.
232           """           """
233           self.verbose=verbose           self.verbose=verbose or True
234           self.show_details= show_details and self.verbose           self.show_details= show_details and self.verbose
235           self.__pde_v.setTolerance(sub_rtol)           rtol=self.getTolerance()
236           self.__pde_p.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
237           ATOL=self.getTolerance()           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
238           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL  
239           #########################################################################################################################           num_corrections=0
240           #           converged=False
241           #   we solve:           p=p0
242           #             norm_r=None
243           #      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:
244           #                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
245           #   residual is                 Qp=self.__Q(p)
246           #                 norm_v=self.__L2(v)
247           #    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)
248           #                 if norm_v == 0.:
249           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                    if norm_Qp == 0.:
250           #                       return v,p
251           #    we use (g - Qp, v) to represent the residual. not that                    else:
252           #                      fac=norm_Qp
253           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                 else:
254           #                    if norm_Qp == 0.:
255           #   while the initial residual is                      fac=norm_v
256           #                    else:
257           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                      fac=2./(1./norm_v+1./norm_Qp)
258           #                   ATOL=(atol+rtol*fac)
259           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                 if self.verbose:
260           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 norm of v = %e."%norm_v
261           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
262           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
263           self.__pde_v.setValue(r=Data())                 if norm_r == None or norm_r>ATOL:
264           # start CG                     if num_corrections>max_num_corrections:
265           r=ArithmeticTuple(d0, v00)                           raise ValueError,"maximum number of correction steps reached."
266           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)
267           return r[1],p                     num_corrections+=1
268                   else:
269                       converged=True
270             return v,p
271    #
272    #              
273    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
274    #               #===========================================================================
275    #               norm_r_hat=self.__L2(r_hat)
276    #               norm_v=self.__L2(v)
277    #               norm_g=self.__L2(g)
278    #               norm_gv=self.__L2(g-v)
279    #               norm_Qp=self.__L2(Qp)
280    #               norm_gQp=self.__L2(g-Qp)
281    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
282    #               fac=min(norm_v,norm_Qp,norm_gv)
283    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
284    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
285    #               if r!=None:
286    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
287    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
288    #                   self.setSubProblemTolerance(sub_tol)
289    #                   print "subtol_new=",self.getSubProblemTolerance()
290    #               print "norm_v = ",norm_v
291    #               print "norm_gv = ",norm_gv
292    #               print "norm_Qp = ",norm_Qp
293    #               print "norm_gQp = ",norm_gQp
294    #               print "norm_g = ",norm_g
295    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
296    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
297    #               if fac == 0:
298    #                   if self.verbose: print "DarcyFlux: trivial case!"
299    #                   return v,p
300    #               #===============================================================================
301    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
302    #               # norm_Qp=self.__L2(Qp)
303    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
304    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
305    #
306    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
307    #               if self.verbose:
308    #                   print "DarcyFlux: residual = %e"%norm_r_hat
309    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
310    #               if norm_r_hat <= ATOL:
311    #                   print "DarcyFlux: iteration finalized."
312    #                   converged=True
313    #               else:
314    #                   # 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)
315    #                   # 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)
316    #                   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)
317    #               print "norm_r =",norm_r
318    #         return v,p
319        def __L2(self,v):
320             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
321    
322        def __Q(self,p):
323              return util.tensor_mult(self.__permeability,util.grad(p))
324    
325      def __Aprod_PCG(self,dp):      def __Aprod(self,dp):
326              self.__pde_v.setTolerance(self.getSubProblemTolerance())
327            if self.show_details: print "DarcyFlux: Applying operator"            if self.show_details: print "DarcyFlux: Applying operator"
328            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)            Qdp=self.__Q(dp)
329            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
           self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())  
330            du=self.__pde_v.getSolution(verbose=self.show_details)            du=self.__pde_v.getSolution(verbose=self.show_details)
331            return ArithmeticTuple(mQdp,du)            return Qdp+du
332        def __inner_GMRES(self,r,s):
333             return util.integrate(util.inner(r,s))
334    
335      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
336           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  
337    
338      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
339              self.__pde_p.setTolerance(self.getSubProblemTolerance())
340            if self.show_details: print "DarcyFlux: Applying preconditioner"            if self.show_details: print "DarcyFlux: Applying preconditioner"
341            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())
342            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution(verbose=self.show_details)
343    
344        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
345            """
346            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
347            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
348            Note that C{g} and C{f} are used, see L{setValue}.
349    
350            @param p: pressure.
351            @type p: scalar value on the domain (e.g. L{Data}).
352            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
353            @type fixed_flux: vector values on the domain (e.g. L{Data}).
354            @param tol: relative tolerance to be used.
355            @type tol: positive C{float}.
356            @return: flux
357            @rtype: L{Data}
358            @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}}
359                   for the permeability M{k_{ij}}
360            """
361            self.__pde_v.setTolerance(self.getSubProblemTolerance())
362            g=self.__g
363            f=self.__f
364            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
365            if p == None:
366               self.__pde_v.setValue(Y=g)
367            else:
368               self.__pde_v.setValue(Y=g-self.__Q(p))
369            return self.__pde_v.getSolution(verbose=show_details)
370    
371  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
372       """       """
373       solves       solves
374    
375            -(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}
376                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 378  class StokesProblemCartesian(Homogeneous
378            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
379            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
380    
381       if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
382    
383       typical usage:       typical usage:
384    
# Line 357  class StokesProblemCartesian(Homogeneous Line 402  class StokesProblemCartesian(Homogeneous
402           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
403           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
404           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
405                
406           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
407           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
408           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
# Line 378  class StokesProblemCartesian(Homogeneous Line 423  class StokesProblemCartesian(Homogeneous
423          @param stress: initial stress          @param stress: initial stress
424      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
425          @note: All values needs to be set.          @note: All values needs to be set.
426    
427          """          """
428          self.eta=eta          self.eta=eta
429          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
430      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
431          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
432          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
433              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
434              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
435      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
436          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
437          self.__f=f          self.__f=f
438          self.__surface_stress=surface_stress          self.__surface_stress=surface_stress
439          self.__stress=stress          self.__stress=stress
440    
 #===============================================================================================================  
441       def inner_pBv(self,p,v):       def inner_pBv(self,p,v):
442           """           """
443           returns inner product of element p and div(v)           returns inner product of element p and div(v)
# Line 433  class StokesProblemCartesian(Homogeneous Line 477  class StokesProblemCartesian(Homogeneous
477           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
478    
479           @param p: a pressure           @param p: a pressure
480           @param v0: a initial guess for the value v to return.           @param v0: a initial guess for the value v to return.
481           @return: v given as M{v= A^{-1} (f-B^*p)}           @return: v given as M{v= A^{-1} (f-B^*p)}
482           """           """
483           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
# Line 448  class StokesProblemCartesian(Homogeneous Line 492  class StokesProblemCartesian(Homogeneous
492    
493           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no v calculation implemented."
494    
495    
496       def norm_Bv(self,v):       def norm_Bv(self,v):
497          """          """
498          Returns Bv (overwrite).          Returns Bv (overwrite).
# Line 463  class StokesProblemCartesian(Homogeneous Line 507  class StokesProblemCartesian(Homogeneous
507           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
508    
509           @param p: a pressure increment           @param p: a pressure increment
510           @return: the solution of M{Av=B^*p}           @return: the solution of M{Av=B^*p}
511           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
512           """           """
513           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())

Legend:
Removed from v.2251  
changed lines
  Added in v.2288

  ViewVC Help
Powered by ViewVC 1.1.26