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

  ViewVC Help
Powered by ViewVC 1.1.26