/[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 1467 by gross, Wed Apr 2 08:10:37 2008 UTC revision 2351 by gross, Tue Mar 31 08:26:41 2009 UTC
# Line 1  Line 1 
1  # $Id:$  ########################################################
2  #  #
3  #######################################################  # Copyright (c) 2003-2008 by University of Queensland
4    # Earth Systems Science Computational Center (ESSCC)
5    # http://www.uq.edu.au/esscc
6  #  #
7  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
8  #  # Licensed under the Open Software License version 3.0
9  #                http://esscc.uq.edu.au  # http://www.opensource.org/licenses/osl-3.0.php
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
10  #  #
11    ########################################################
12    
13    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
16    Primary Business: Queensland, Australia"""
17    __license__="""Licensed under the Open Software License version 3.0
18    http://www.opensource.org/licenses/osl-3.0.php"""
19    __url__="https://launchpad.net/escript-finley"
20    
21  """  """
22  Some models for flow  Some models for flow
# Line 24  Some models for flow Line 30  Some models for flow
30  """  """
31    
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37  from pdetools import HomogeneousSaddlePointProblem  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38    
39    class DarcyFlow(object):
40        """
41        solves the problem
42    
43        M{u_i+k_{ij}*p_{,j} = g_i}
44        M{u_{i,i} = f}
45    
46        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.
49        """
50    
51        def __init__(self, domain, weight=None, useReduced=False):
52            """
53            initializes the Darcy flux problem
54            @param domain: domain of the problem
55            @type domain: L{Domain}
56            """
57            self.domain=domain
58            if weight == None:
59               self.__l=10.*util.longestEdge(self.domain)**2
60            else:
61               self.__l=weight
62            self.__pde_v=LinearPDESystem(domain)
63            if useReduced: self.__pde_v.setReducedOrderOn()
64            self.__pde_v.setSymmetryOn()
65            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
66            # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
67            self.__pde_p=LinearSinglePDE(domain)
68            self.__pde_p.setSymmetryOn()
69            if useReduced: self.__pde_p.setReducedOrderOn()
70            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
71            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
72            self.setTolerance()
73            self.setAbsoluteTolerance()
74            self.setSubProblemTolerance()
75    
76        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
77            """
78            assigns values to model parameters
79    
80            @param f: volumetic sources/sinks
81            @type f: scalar value on the domain (e.g. L{Data})
82            @param g: flux sources/sinks
83            @type g: vector values on the domain (e.g. L{Data})
84            @param location_of_fixed_pressure: mask for locations where pressure is fixed
85            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
86            @param location_of_fixed_flux:  mask for locations where flux is fixed.
87            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
88            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
89                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
90                                 C{v} on the main diagonal is used.
91            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
92    
93            @note: the values of parameters which are not set by calling C{setValue} are not altered.
94            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
95                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
96                   is along the M{x_i} axis.
97            """
98            if f !=None:
99               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
100               if f.isEmpty():
101                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
102               else:
103                   if f.getRank()>0: raise ValueError,"illegal rank of f."
104               self.__f=f
105            if g !=None:
106               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
107               if g.isEmpty():
108                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
109               else:
110                 if not g.getShape()==(self.domain.getDim(),):
111                   raise ValueError,"illegal shape of g"
112               self.__g=g
113    
114            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
115            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
116    
117            if permeability!=None:
118               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
119               if perm.getRank()==0:
120                   perm=perm*util.kronecker(self.domain.getDim())
121               elif perm.getRank()==1:
122                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
123                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
124               elif perm.getRank()==2:
125                  pass
126               else:
127                  raise ValueError,"illegal rank of permeability."
128               self.__permeability=perm
129               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
130    
131        def setTolerance(self,rtol=1e-4):
132            """
133            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
134    
135            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
136    
137            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}}.
138    
139            @param rtol: relative tolerance for the pressure
140            @type rtol: non-negative C{float}
141            """
142            if rtol<0:
143                raise ValueError,"Relative tolerance needs to be non-negative."
144            self.__rtol=rtol
145        def getTolerance(self):
146            """
147            returns the relative tolerance
148    
149            @return: current relative tolerance
150            @rtype: C{float}
151            """
152            return self.__rtol
153    
154        def setAbsoluteTolerance(self,atol=0.):
155            """
156            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
157    
158            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
159    
160            where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
161    
162            @param atol: absolute tolerance for the pressure
163            @type atol: non-negative C{float}
164            """
165            if atol<0:
166                raise ValueError,"Absolute tolerance needs to be non-negative."
167            self.__atol=atol
168        def getAbsoluteTolerance(self):
169           """
170           returns the absolute tolerance
171          
172           @return: current absolute tolerance
173           @rtype: C{float}
174           """
175           return self.__atol
176    
177        def setSubProblemTolerance(self,rtol=None):
178             """
179             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
180             C{self.getTolerance()**2} is used.
181    
182             @param rtol: relative tolerence
183             @type rtol: positive C{float}
184             """
185             if rtol == None:
186                  if self.getTolerance()<=0.:
187                      raise ValueError,"A positive relative tolerance must be set."
188                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
189             else:
190                 if rtol<=0:
191                     raise ValueError,"sub-problem tolerance must be positive."
192                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
193    
194        def getSubProblemTolerance(self):
195             """
196             Returns the subproblem reduction factor.
197    
198             @return: subproblem reduction factor
199             @rtype: C{float}
200             """
201             return self.__sub_tol
202    
203        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
204             """
205             solves the problem.
206    
207             The iteration is terminated if the residual norm is less then self.getTolerance().
208    
209             @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
210             @type u0: vector value on the domain (e.g. L{Data}).
211             @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
212             @type p0: scalar value on the domain (e.g. L{Data}).
213             @param verbose: if set some information on iteration progress are printed
214             @type verbose: C{bool}
215             @param show_details:  if set information on the subiteration process are printed.
216             @type show_details: C{bool}
217             @return: flux and pressure
218             @rtype: C{tuple} of L{Data}.
219    
220             @note: The problem is solved as a least squares form
221    
222             M{(I+D^*D)u+Qp=D^*f+g}
223             M{Q^*u+Q^*Qp=Q^*g}
224    
225             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
226             We eliminate the flux form the problem by setting
227    
228             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
229    
230             form the first equation. Inserted into the second equation we get
231    
232             M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure
233    
234             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
235             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
236             """
237             self.verbose=verbose or True
238             self.show_details= show_details and self.verbose
239             rtol=self.getTolerance()
240             atol=self.getAbsoluteTolerance()
241             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
242    
243             num_corrections=0
244             converged=False
245             p=p0
246             norm_r=None
247             while not converged:
248                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
249                   Qp=self.__Q(p)
250                   norm_v=self.__L2(v)
251                   norm_Qp=self.__L2(Qp)
252                   if norm_v == 0.:
253                      if norm_Qp == 0.:
254                         return v,p
255                      else:
256                        fac=norm_Qp
257                   else:
258                      if norm_Qp == 0.:
259                        fac=norm_v
260                      else:
261                        fac=2./(1./norm_v+1./norm_Qp)
262                   ATOL=(atol+rtol*fac)
263                   if self.verbose:
264                        print "DarcyFlux: L2 norm of v = %e."%norm_v
265                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
266                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
267                   if norm_r == None or norm_r>ATOL:
268                       if num_corrections>max_num_corrections:
269                             raise ValueError,"maximum number of correction steps reached."
270                       p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
271                       num_corrections+=1
272                   else:
273                       converged=True
274             return v,p
275    #
276    #              
277    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
278    #               #===========================================================================
279    #               norm_r_hat=self.__L2(r_hat)
280    #               norm_v=self.__L2(v)
281    #               norm_g=self.__L2(g)
282    #               norm_gv=self.__L2(g-v)
283    #               norm_Qp=self.__L2(Qp)
284    #               norm_gQp=self.__L2(g-Qp)
285    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
286    #               fac=min(norm_v,norm_Qp,norm_gv)
287    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
288    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
289    #               if r!=None:
290    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
291    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
292    #                   self.setSubProblemTolerance(sub_tol)
293    #                   print "subtol_new=",self.getSubProblemTolerance()
294    #               print "norm_v = ",norm_v
295    #               print "norm_gv = ",norm_gv
296    #               print "norm_Qp = ",norm_Qp
297    #               print "norm_gQp = ",norm_gQp
298    #               print "norm_g = ",norm_g
299    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
300    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
301    #               if fac == 0:
302    #                   if self.verbose: print "DarcyFlux: trivial case!"
303    #                   return v,p
304    #               #===============================================================================
305    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
306    #               # norm_Qp=self.__L2(Qp)
307    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
308    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
309    #
310    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
311    #               if self.verbose:
312    #                   print "DarcyFlux: residual = %e"%norm_r_hat
313    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
314    #               if norm_r_hat <= ATOL:
315    #                   print "DarcyFlux: iteration finalized."
316    #                   converged=True
317    #               else:
318    #                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)
319    #                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)
320    #                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
321    #               print "norm_r =",norm_r
322    #         return v,p
323        def __L2(self,v):
324             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
325    
326        def __Q(self,p):
327              return util.tensor_mult(self.__permeability,util.grad(p))
328    
329        def __Aprod(self,dp):
330              self.__pde_v.setTolerance(self.getSubProblemTolerance())
331              if self.show_details: print "DarcyFlux: Applying operator"
332              Qdp=self.__Q(dp)
333              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
334              du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
335              return Qdp+du
336        def __inner_GMRES(self,r,s):
337             return util.integrate(util.inner(r,s))
338    
339        def __inner_PCG(self,p,r):
340             return util.integrate(util.inner(self.__Q(p), r))
341    
342        def __Msolve_PCG(self,r):
343              self.__pde_p.setTolerance(self.getSubProblemTolerance())
344              if self.show_details: print "DarcyFlux: Applying preconditioner"
345              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
346              return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
347    
348        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
349            """
350            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
351            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
352            Note that C{g} and C{f} are used, see L{setValue}.
353    
354            @param p: pressure.
355            @type p: scalar value on the domain (e.g. L{Data}).
356            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
357            @type fixed_flux: vector values on the domain (e.g. L{Data}).
358            @param tol: relative tolerance to be used.
359            @type tol: positive C{float}.
360            @return: flux
361            @rtype: L{Data}
362            @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
363                   for the permeability M{k_{ij}}
364            """
365            self.__pde_v.setTolerance(self.getSubProblemTolerance())
366            g=self.__g
367            f=self.__f
368            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
369            if p == None:
370               self.__pde_v.setValue(Y=g)
371            else:
372               self.__pde_v.setValue(Y=g-self.__Q(p))
373            return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
374    
375  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
376        """       """
377        solves       solves
378    
379            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
380                  u_{i,i}=0                  u_{i,i}=0
381    
382            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
383            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
384    
385        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
386    
387        typical usage:       typical usage:
388    
389              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
390              sp.setTolerance()              sp.setTolerance()
391              sp.initialize(...)              sp.initialize(...)
392              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
393        """       """
394        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
395             """
396             initialize the Stokes Problem
397    
398             @param domain: domain of the problem. The approximation order needs to be two.
399             @type domain: L{Domain}
400             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
401             """
402           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
403           self.domain=domain           self.domain=domain
404           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
405           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
406           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
407           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
408                         # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
409    
410           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
411           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
412             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
413           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
414    
415           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
416           self.__pde_proj.setReducedOrderOn()          """
417           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
418           self.__pde_proj.setValue(D=1.)  
419            @param f: external force
420            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
421            @param fixed_u_mask: mask of locations with fixed velocity.
422            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
423            @param eta: viscosity
424            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
425            @param surface_stress: normal surface stress
426            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
427            @param stress: initial stress
428        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
429            @note: All values needs to be set.
430    
431        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):          """
432          self.eta=eta          self.eta=eta
433          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
434      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
435          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
436          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
437              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
438              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
439      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
440          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
441            self.__f=f
442            self.__surface_stress=surface_stress
443            self.__stress=stress
444    
445         def inner_pBv(self,p,v):
446             """
447             returns inner product of element p and div(v)
448    
449             @param p: a pressure increment
450             @param v: a residual
451             @return: inner product of element p and div(v)
452             @rtype: C{float}
453             """
454             return util.integrate(-p*util.div(v))
455    
456         def inner_p(self,p0,p1):
457             """
458             Returns inner product of p0 and p1
459    
460        def B(self,arg):           @param p0: a pressure
461           d=util.div(arg)           @param p1: a pressure
462           self.__pde_proj.setValue(Y=d)           @return: inner product of p0 and p1
463           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           @rtype: C{float}
464           return self.__pde_proj.getSolution(verbose=self.show_details)           """
465             s0=util.interpolate(p0/self.eta,Function(self.domain))
466        def inner(self,p0,p1):           s1=util.interpolate(p1/self.eta,Function(self.domain))
          s0=util.interpolate(p0,Function(self.domain))  
          s1=util.interpolate(p1,Function(self.domain))  
467           return util.integrate(s0*s1)           return util.integrate(s0*s1)
468    
469        def getStress(self,u):       def norm_v(self,v):
470           mg=util.grad(u)           """
471           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
472    
473             @param v: a velovity
474             @return: norm of v
475             @rtype: non-negative C{float}
476             """
477             return util.sqrt(util.integrate(util.length(util.grad(v))))
478    
479         def getV(self, p, v0):
480             """
481             return the value for v for a given p (overwrite)
482    
483             @param p: a pressure
484             @param v0: a initial guess for the value v to return.
485             @return: v given as M{v= A^{-1} (f-B^*p)}
486             """
487             self.__pde_u.setTolerance(self.getSubProblemTolerance())
488             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
489             if self.__stress.isEmpty():
490                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
491             else:
492                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
493             out=self.__pde_u.getSolution(verbose=self.show_details)
494             return  out
495    
496    
497             raise NotImplementedError,"no v calculation implemented."
498    
499        def solve_A(self,u,p):  
500         def norm_Bv(self,v):
501            """
502            Returns Bv (overwrite).
503    
504            @rtype: equal to the type of p
505            @note: boundary conditions on p should be zero!
506            """
507            return util.sqrt(util.integrate(util.div(v)**2))
508    
509         def solve_AinvBt(self,p):
510           """           """
511           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
512    
513             @param p: a pressure increment
514             @return: the solution of M{Av=B^*p}
515             @note: boundary conditions on v should be zero!
516           """           """
517           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
518           self.__pde_u.setValue(X=-self.getStress(u)+p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
519           return  self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
520             return  out
521    
522        def solve_prec(self,p):       def solve_precB(self,v):
523             """
524             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
525             with accuracy L{self.getSubProblemTolerance()} (overwrite).
526    
527             @param v: velocity increment
528             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
529             @note: boundary conditions on p are zero.
530             """
531             self.__pde_prec.setValue(Y=-util.div(v))
532           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
533           self.__pde_prec.setValue(Y=p)           return self.__pde_prec.getSolution(verbose=self.show_details)
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.Lsup(v)  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)  
           self.iter+=1  
           if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium_GMRES(self,norm_r,norm_b):  
           if self.verbose: print "GMRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_b*self.getTolerance())  
           self.iter+=1  
           if norm_r <= norm_b*self.getTolerance():  
               if self.verbose: print "GMRES terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  

Legend:
Removed from v.1467  
changed lines
  Added in v.2351

  ViewVC Help
Powered by ViewVC 1.1.26