/[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 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC revision 2344 by jfenwick, Mon Mar 30 02:13: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      Represents and 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      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
     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):      def __init__(self, domain,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            self.__l=util.longestEdge(self.domain)**2
59          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
60          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          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=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()
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.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          """          """
74          Assigns values to model parameters.          assigns values to model parameters
75    
76          @param f: volumetric sources/sinks          @param f: volumetic sources/sinks
77          @type f: scalar value on the domain, e.g. L{Data}          @type f: scalar value on the domain (e.g. L{Data})
78          @param g: flux sources/sinks          @param g: flux sources/sinks
79          @type g: vector value on the domain, e.g. L{Data}          @type g: vector values on the domain (e.g. L{Data})
80          @param location_of_fixed_pressure: mask for locations where pressure is fixed          @param location_of_fixed_pressure: mask for locations where pressure is fixed
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 value 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          @param permeability: permeability tensor. If scalar C{s} is given the tensor with
85                               tensor with C{s} on the main diagonal is used. If                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with
86                               vector C{v} is given the tensor with C{v} on the                               C{v} on the main diagonal is used.
87                               main diagonal is used.          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
88          @type permeability: scalar, vector or tensor values on the domain, e.g.  
89                              L{Data}          @note: the values of parameters which are not set by calling C{setValue} are not altered.
90            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
91          @note: the values of parameters which are not set by calling                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
92                 C{setValue} are not altered                 is along the M{x_i} axis.
         @note: at any point on the boundary of the domain the pressure  
                (C{location_of_fixed_pressure}) >0 or the normal component of  
                the flux (C{location_of_fixed_flux[i]}) >0 if the direction of  
                the normal 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"))
# Line 98  class DarcyFlow(object): Line 97  class DarcyFlow(object):
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():
# Line 125  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):
128            """
129            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
130    
131            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
132    
133            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}}.
134    
135      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):          @param rtol: relative tolerance for the pressure
136            @type rtol: non-negative C{float}
137            """
138            if rtol<0:
139                raise ValueError,"Relative tolerance needs to be non-negative."
140            self.__rtol=rtol
141        def getTolerance(self):
142          """          """
143          Returns the flux for a given pressure C{p}.          returns the relative tolerance
144    
145          The flux is equal to C{fixed_flux} on locations where          @return: current relative tolerance
146          C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}          @rtype: C{float}
         and C{f} are used.  
   
         @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 by  
                            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 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}}  
147          """          """
148          self.__pde_v.setTolerance(tol)          return self.__rtol
         self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)  
         return self.__pde_v.getSolution(verbose=show_details)  
149    
150      def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def setAbsoluteTolerance(self,atol=0.):
151          """          """
152          Solves the problem.          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
153    
154            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
155    
156            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    
         The iteration is terminated if the error in the pressure is less than  
         M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand  
         side (see escript user's guide for details).  
   
         @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.  
         @type u0: vector value on the domain, e.g. L{Data}  
         @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.  
         @type p0: scalar value on the domain, e.g. L{Data}  
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}
160          @param rtol: relative tolerance for the pressure          """
161          @type rtol: non-negative C{float}          if atol<0:
162          @param sub_rtol: tolerance to be used in the sub iteration. It is              raise ValueError,"Absolute tolerance needs to be non-negative."
163                           recommended that M{sub_rtol<rtol*5.e-3}          self.__atol=atol
164          @type sub_rtol: positive-negative C{float}      def getAbsoluteTolerance(self):
165          @param verbose: if True information on iteration progress is printed         """
166          @type verbose: C{bool}         returns the absolute tolerance
167          @param show_details: if True information on the sub-iteration process        
168                               is printed         @return: current absolute tolerance
169          @type show_details: C{bool}         @rtype: C{float}
170          @return: flux and pressure         """
171          @rtype: C{tuple} of L{Data}         return self.__atol
172    
173          @note: the problem is solved in a least squares formulation:      def setSubProblemTolerance(self,rtol=None):
174             """
175          M{(I+D^*D)u+Qp=D^*f+g}           Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
176             C{self.getTolerance()**2} is used.
177          M{Q^*u+Q^*Qp=Q^*g}  
178             @param rtol: relative tolerence
179          where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the           @type rtol: positive C{float}
180          permeability M{k_{ij}}. We eliminate the flux from the problem by           """
181          setting           if rtol == None:
182                  if self.getTolerance()<=0.:
183          M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}                    raise ValueError,"A positive relative tolerance must be set."
184                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
185          from the first equation. Inserted into the second equation we get           else:
186                 if rtol<=0:
187          M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}                   raise ValueError,"sub-problem tolerance must be positive."
188          on C{location_of_fixed_pressure}               self.__sub_tol=max(util.EPSILON**(0.75),rtol)
189    
190          which is solved using the PCG method (precondition is M{Q^*Q}).      def getSubProblemTolerance(self):
191          In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}           """
192          need to be solved using a sub-iteration scheme.           Returns the subproblem reduction factor.
193          """  
194          self.verbose=verbose           @return: subproblem reduction factor
195          self.show_details= show_details and self.verbose           @rtype: C{float}
196          self.__pde_v.setTolerance(sub_rtol)           """
197          self.__pde_p.setTolerance(sub_rtol)           return self.__sub_tol
198          u2=u0*self.__pde_v.getCoefficient("q")  
199          #      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
200          # first the reference velocity is calculated from           """
201          #           solves the problem.
202          #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)  
203          #           The iteration is terminated if the residual norm is less then self.getTolerance().
204          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)  
205          u_ref=self.__pde_v.getSolution(verbose=show_details)           @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          if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)           @type u0: vector value on the domain (e.g. L{Data}).
207          self.__pde_v.setValue(r=Data())           @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}).
209          #   and then we calculate a reference pressure           @param verbose: if set some information on iteration progress are printed
210          #           @type verbose: C{bool}
211          #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)           @param show_details:  if set information on the subiteration process are printed.
212          #           @type show_details: C{bool}
213          self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)           @return: flux and pressure
214          p_ref=self.__pde_p.getSolution(verbose=self.show_details)           @rtype: C{tuple} of L{Data}.
215          if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)  
216          self.__pde_p.setValue(r=Data())           @note: The problem is solved as a least squares form
217          #  
218          #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref           M{(I+D^*D)u+Qp=D^*f+g}
219          #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref           M{Q^*u+Q^*Qp=Q^*g}
220          #  
221          #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du           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
223          #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref  
224          #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
225          #  
226          #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)           form the first equation. Inserted into the second equation we get
227          #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)  
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
229          #  we use the (du,Qdp) to represent the resudual  
230          #  Q^*Q is a preconditioner           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.
232          #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>           """
233          #           self.verbose=verbose or True
234          Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))           self.show_details= show_details and self.verbose
235          norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))           rtol=self.getTolerance()
236          ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)           atol=self.getAbsoluteTolerance()
237          if not ATOL>0:           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
238              raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)  
239          if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)           num_corrections=0
240          #           converged=False
241          #   caclulate the initial residual           p=p0
242          #           norm_r=None
243          self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())           while not converged:
244          du=self.__pde_v.getSolution(verbose=show_details)                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
245          r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)                 Qp=self.__Q(p)
246          dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                 norm_v=self.__L2(v)
247          util.saveVTK("d.vtu",p=dp,p_ref=p_ref)                 norm_Qp=self.__L2(Qp)
248          return u_ref+r[1],dp                 if norm_v == 0.:
249                      if norm_Qp == 0.:
250      def __Aprod_PCG(self,p):                       return v,p
251          if self.show_details: print "DarcyFlux: Applying operator"                    else:
252          Qp=util.tensor_mult(self.__permeability,util.grad(p))                      fac=norm_Qp
253          self.__pde_v.setValue(Y=Qp,X=Data())                 else:
254          w=self.__pde_v.getSolution(verbose=self.show_details)                    if norm_Qp == 0.:
255          return ArithmeticTuple(-Qp,w)                      fac=norm_v
256                      else:
257                        fac=2./(1./norm_v+1./norm_Qp)
258                   ATOL=(atol+rtol*fac)
259                   if self.verbose:
260                        print "DarcyFlux: L2 norm of v = %e."%norm_v
261                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
262                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
263                   if norm_r == None or norm_r>ATOL:
264                       if num_corrections>max_num_corrections:
265                             raise ValueError,"maximum number of correction steps reached."
266                       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                       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(self,dp):
326              self.__pde_v.setTolerance(self.getSubProblemTolerance())
327              if self.show_details: print "DarcyFlux: Applying operator"
328              Qdp=self.__Q(dp)
329              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
330              du=self.__pde_v.getSolution(verbose=self.show_details)
331              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))
         out=-util.integrate(util.inner(a,r[0]+r[1]))  
         return out  
337    
338      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
339          if self.show_details: print "DarcyFlux: Applying preconditioner"            self.__pde_p.setTolerance(self.getSubProblemTolerance())
340          self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))            if self.show_details: print "DarcyFlux: Applying preconditioner"
341          return self.__pde_p.getSolution(verbose=self.show_details)            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)
343    
344  class StokesProblemCartesian(HomogeneousSaddlePointProblem):      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
345        """          """
346        Represents and solves the problem          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        M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
372         """
373         solves
374    
375        M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
376                    u_{i,i}=0
377    
378        M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}            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
380    
381        If C{surface_stress} is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
382    
383        Typical usage::       typical usage:
384    
385            sp = StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
386            sp.setTolerance()              sp.setTolerance()
387            sp.initialize(...)              sp.initialize(...)
388            v,p = sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
389        """       """
390        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
391           """           """
392           Initializes the Stokes Problem.           initialize the Stokes Problem
393    
394           @param domain: domain of the problem. The approximation order needs           @param domain: domain of the problem. The approximation order needs to be two.
                         to be two.  
395           @type domain: L{Domain}           @type domain: L{Domain}
396           @warning: The approximation order needs to be two otherwise you may           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
                    see oscillations in the pressure.  
397           """           """
398           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
399           self.domain=domain           self.domain=domain
# Line 318  class StokesProblemCartesian(Homogeneous Line 408  class StokesProblemCartesian(Homogeneous
408           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
409           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
410    
411           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
412           self.__pde_proj.setReducedOrderOn()          """
413           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
414           self.__pde_proj.setValue(D=1.)  
415            @param f: external force
416        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
417           """          @param fixed_u_mask: mask of locations with fixed velocity.
418           Assigns values to the model parameters.          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
419            @param eta: viscosity
420           @param f: external force          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
421           @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          @param surface_stress: normal surface stress
422           @param fixed_u_mask: mask of locations with fixed velocity          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
423           @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}          @param stress: initial stress
424                               or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
425           @param eta: viscosity          @note: All values needs to be set.
426           @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar  
427           @param surface_stress: normal surface stress          """
428           @type surface_stress: L{Vector} object on L{FunctionSpace},          self.eta=eta
429                                 L{FunctionOnBoundary} or similar          A =self.__pde_u.createCoefficient("A")
430           @param stress: initial stress      self.__pde_u.setValue(A=Data())
431           @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or          for i in range(self.domain.getDim()):
432                         similar          for j in range(self.domain.getDim()):
433           @note: All values need to be set.              A[i,j,j,i] += 1.
434           """              A[i,j,i,j] += 1.
435           self.eta=eta      self.__pde_prec.setValue(D=1/self.eta)
436           A =self.__pde_u.createCoefficient("A")          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
437           self.__pde_u.setValue(A=Data())          self.__f=f
438           for i in range(self.domain.getDim()):          self.__surface_stress=surface_stress
439               for j in range(self.domain.getDim()):          self.__stress=stress
440                   A[i,j,j,i] += 1.  
441                   A[i,j,i,j] += 1.       def inner_pBv(self,p,v):
442           self.__pde_prec.setValue(D=1/self.eta)           """
443           self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)           returns inner product of element p and div(v)
444           self.__stress=stress  
445             @param p: a pressure increment
446        def B(self,v):           @param v: a residual
447           """           @return: inner product of element p and div(v)
448           Returns M{div(v)}.           @rtype: C{float}
          @return: M{div(v)}  
          @rtype: equal to the type of p  
   
          @note: Boundary conditions on p should be zero!  
          """  
          if self.show_details: print "apply divergence:"  
          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  
          @return: inner product of p and Bv  
          @rtype: equal to the type of p  
449           """           """
450           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)  
451    
452        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
453           """           """
454           Returns inner product of element p0 and p1 (overwrite).           Returns inner product of p0 and p1
455    
456           @type p0: equal to the type of p           @param p0: a pressure
457           @type p1: equal to the type of p           @param p1: a pressure
458           @return: inner product of p0 and p1           @return: inner product of p0 and p1
459           @rtype: equal to the type of p           @rtype: C{float}
460           """           """
461           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
462           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
463           return util.integrate(s0*s1)           return util.integrate(s0*s1)
464    
465        def inner_v(self,v0,v1):       def norm_v(self,v):
466           """           """
467           Returns inner product of two elements v0 and v1 (overwrite).           returns the norm of v
468    
469           @type v0: equal to the type of v           @param v: a velovity
470           @type v1: equal to the type of v           @return: norm of v
471           @return: inner product of v0 and v1           @rtype: non-negative C{float}
          @rtype: equal to the type of v  
472           """           """
473           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))  
474    
475        def solve_A(self,u,p):       def getV(self, p, v0):
476           """           """
477           Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).           return the value for v for a given p (overwrite)
478    
479             @param p: a pressure
480             @param v0: a initial guess for the value v to return.
481             @return: v given as M{v= A^{-1} (f-B^*p)}
482           """           """
          if self.show_details: print "solve for velocity:"  
483           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
484             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
485           if self.__stress.isEmpty():           if self.__stress.isEmpty():
486              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))
487           else:           else:
488              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))
489           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
490           return out           return  out
491    
492    
493        def solve_prec(self,p):           raise NotImplementedError,"no v calculation implemented."
494    
495    
496         def norm_Bv(self,v):
497            """
498            Returns Bv (overwrite).
499    
500            @rtype: equal to the type of p
501            @note: boundary conditions on p should be zero!
502            """
503            return util.sqrt(util.integrate(util.div(v)**2))
504    
505         def solve_AinvBt(self,p):
506           """           """
507           Applies the preconditioner.           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
508    
509             @param p: a pressure increment
510             @return: the solution of M{Av=B^*p}
511             @note: boundary conditions on v should be zero!
512           """           """
513           if self.show_details: print "apply preconditioner:"           self.__pde_u.setTolerance(self.getSubProblemTolerance())
514           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
515           self.__pde_prec.setValue(Y=p)           out=self.__pde_u.getSolution(verbose=self.show_details)
516           q=self.__pde_prec.getSolution(verbose=self.show_details)           return  out
517           return q  
518         def solve_precB(self,v):
519             """
520             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
521             with accuracy L{self.getSubProblemTolerance()} (overwrite).
522    
523             @param v: velocity increment
524             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
525             @note: boundary conditions on p are zero.
526             """
527             self.__pde_prec.setValue(Y=-util.div(v))
528             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
529             return self.__pde_prec.getSolution(verbose=self.show_details)

Legend:
Removed from v.2169  
changed lines
  Added in v.2344

  ViewVC Help
Powered by ViewVC 1.1.26