/[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 1519 by artak, Tue Apr 22 03:45:36 2008 UTC revision 2264 by gross, Wed Feb 11 06:48:28 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__="http://www.uq.edu.au/esscc/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,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            self.__pde_v=LinearPDESystem(domain)
59            if useReduced: self.__pde_v.setReducedOrderOn()
60            self.__pde_v.setSymmetryOn()
61            self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
62            self.__pde_p=LinearSinglePDE(domain)
63            self.__pde_p.setSymmetryOn()
64            if useReduced: self.__pde_p.setReducedOrderOn()
65            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
66            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
67            self.setTolerance()
68            self.setAbsoluteTolerance()
69            self.setSubProblemTolerance()
70    
71        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
72            """
73            assigns values to model parameters
74    
75            @param f: volumetic sources/sinks
76            @type f: scalar value on the domain (e.g. L{Data})
77            @param g: flux sources/sinks
78            @type g: vector values on the domain (e.g. L{Data})
79            @param location_of_fixed_pressure: mask for locations where pressure is fixed
80            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
81            @param location_of_fixed_flux:  mask for locations where flux is fixed.
82            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
83            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
84                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
85                                 C{v} on the main diagonal is used.
86            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
87    
88            @note: the values of parameters which are not set by calling C{setValue} are not altered.
89            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
90                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
91                   is along the M{x_i} axis.
92            """
93            if f !=None:
94               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
95               if f.isEmpty():
96                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
97               else:
98                   if f.getRank()>0: raise ValueError,"illegal rank of f."
99               self.f=f
100            if g !=None:
101               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
102               if g.isEmpty():
103                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
104               else:
105                 if not g.getShape()==(self.domain.getDim(),):
106                   raise ValueError,"illegal shape of g"
107               self.__g=g
108    
109            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
110            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
111    
112            if permeability!=None:
113               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
114               if perm.getRank()==0:
115                   perm=perm*util.kronecker(self.domain.getDim())
116               elif perm.getRank()==1:
117                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
118                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
119               elif perm.getRank()==2:
120                  pass
121               else:
122                  raise ValueError,"illegal rank of permeability."
123               self.__permeability=perm
124               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
125    
126        def setTolerance(self,rtol=1e-4):
127            """
128            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
129    
130            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
131    
132            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}}.
133    
134            @param rtol: relative tolerance for the pressure
135            @type rtol: non-negative C{float}
136            """
137            if rtol<0:
138                raise ValueError,"Relative tolerance needs to be non-negative."
139            self.__rtol=rtol
140        def getTolerance(self):
141            """
142            returns the relative tolerance
143    
144            @return: current relative tolerance
145            @rtype: C{float}
146            """
147            return self.__rtol
148    
149        def setAbsoluteTolerance(self,atol=0.):
150            """
151            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
152    
153            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
154    
155            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}}.
156    
157            @param atol: absolute tolerance for the pressure
158            @type atol: non-negative C{float}
159            """
160            if atol<0:
161                raise ValueError,"Absolute tolerance needs to be non-negative."
162            self.__atol=atol
163        def getAbsoluteTolerance(self):
164           """
165           returns the absolute tolerance
166          
167           @return: current absolute tolerance
168           @rtype: C{float}
169           """
170           return self.__atol
171    
172        def setSubProblemTolerance(self,rtol=None):
173             """
174             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
175             C{self.getTolerance()**2} is used.
176    
177             @param rtol: relative tolerence
178             @type rtol: positive C{float}
179             """
180             if rtol == None:
181                  if self.getTolerance()<=0.:
182                      raise ValueError,"A positive relative tolerance must be set."
183                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
184             else:
185                 if rtol<=0:
186                     raise ValueError,"sub-problem tolerance must be positive."
187                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
188    
189        def getSubProblemTolerance(self):
190             """
191             Returns the subproblem reduction factor.
192    
193             @return: subproblem reduction factor
194             @rtype: C{float}
195             """
196             return self.__sub_tol
197    
198        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
199             """
200             solves the problem.
201    
202             The iteration is terminated if the residual norm is less then self.getTolerance().
203    
204             @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.
205             @type u0: vector value on the domain (e.g. L{Data}).
206             @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.
207             @type p0: scalar value on the domain (e.g. L{Data}).
208             @param verbose: if set some information on iteration progress are printed
209             @type verbose: C{bool}
210             @param show_details:  if set information on the subiteration process are printed.
211             @type show_details: C{bool}
212             @return: flux and pressure
213             @rtype: C{tuple} of L{Data}.
214    
215             @note: The problem is solved as a least squares form
216    
217             M{(I+D^*D)u+Qp=D^*f+g}
218             M{Q^*u+Q^*Qp=Q^*g}
219    
220             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
221             We eliminate the flux form the problem by setting
222    
223             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
224    
225             form the first equation. Inserted into the second equation we get
226    
227             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
228    
229             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
230             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
231             """
232             self.verbose=verbose
233             self.show_details= show_details and self.verbose
234             rtol=self.getTolerance()
235             atol=self.getAbsoluteTolerance()
236             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
237    
238             num_corrections=0
239             converged=False
240             p=p0
241             norm_r=None
242             while not converged:
243                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
244                   Qp=self.__Q(p)
245                   norm_v=self.__L2(v)
246                   norm_Qp=self.__L2(Qp)
247                   if norm_v == 0.:
248                      if norm_Qp == 0.:
249                         return v,p
250                      else:
251                        fac=norm_Qp
252                   else:
253                      if norm_Qp == 0.:
254                        fac=norm_v
255                      else:
256                        fac=2./(1./norm_v+1./norm_Qp)
257                   ATOL=(atol+rtol*fac)
258                   if self.verbose:
259                        print "DarcyFlux: L2 norm of v = %e."%norm_v
260                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
261                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
262                   if norm_r == None or norm_r>ATOL:
263                       if num_corrections>max_num_corrections:
264                             raise ValueError,"maximum number of correction steps reached."
265                       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)
266                       num_corrections+=1
267                   else:
268                       converged=True
269             return v,p
270    #
271    #              
272    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
273    #               #===========================================================================
274    #               norm_r_hat=self.__L2(r_hat)
275    #               norm_v=self.__L2(v)
276    #               norm_g=self.__L2(g)
277    #               norm_gv=self.__L2(g-v)
278    #               norm_Qp=self.__L2(Qp)
279    #               norm_gQp=self.__L2(g-Qp)
280    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
281    #               fac=min(norm_v,norm_Qp,norm_gv)
282    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
283    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
284    #               if r!=None:
285    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
286    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
287    #                   self.setSubProblemTolerance(sub_tol)
288    #                   print "subtol_new=",self.getSubProblemTolerance()
289    #               print "norm_v = ",norm_v
290    #               print "norm_gv = ",norm_gv
291    #               print "norm_Qp = ",norm_Qp
292    #               print "norm_gQp = ",norm_gQp
293    #               print "norm_g = ",norm_g
294    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
295    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
296    #               if fac == 0:
297    #                   if self.verbose: print "DarcyFlux: trivial case!"
298    #                   return v,p
299    #               #===============================================================================
300    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
301    #               # norm_Qp=self.__L2(Qp)
302    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
303    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
304    #
305    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
306    #               if self.verbose:
307    #                   print "DarcyFlux: residual = %e"%norm_r_hat
308    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
309    #               if norm_r_hat <= ATOL:
310    #                   print "DarcyFlux: iteration finalized."
311    #                   converged=True
312    #               else:
313    #                   # 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)
314    #                   # 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)
315    #                   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)
316    #               print "norm_r =",norm_r
317    #         return v,p
318        def __L2(self,v):
319             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
320    
321        def __Q(self,p):
322              return util.tensor_mult(self.__permeability,util.grad(p))
323    
324        def __Aprod(self,dp):
325              self.__pde_v.setTolerance(self.getSubProblemTolerance())
326              if self.show_details: print "DarcyFlux: Applying operator"
327              Qdp=self.__Q(dp)
328              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
329              du=self.__pde_v.getSolution(verbose=self.show_details)
330              return Qdp+du
331        def __inner_GMRES(self,r,s):
332             return util.integrate(util.inner(r,s))
333    
334        def __inner_PCG(self,p,r):
335             return util.integrate(util.inner(self.__Q(p), r))
336    
337        def __Msolve_PCG(self,r):
338              self.__pde_p.setTolerance(self.getSubProblemTolerance())
339              if self.show_details: print "DarcyFlux: Applying preconditioner"
340              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
341              return self.__pde_p.getSolution(verbose=self.show_details)
342    
343    
344        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
345            """
346            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
347            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
348            Note that C{g} and C{f} are used, see L{setValue}.
349    
350            @param p: pressure.
351            @type p: scalar value on the domain (e.g. L{Data}).
352            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
353            @type fixed_flux: vector values on the domain (e.g. L{Data}).
354            @param tol: relative tolerance to be used.
355            @type tol: positive C{float}.
356            @return: flux
357            @rtype: L{Data}
358            @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
359                   for the permeability M{k_{ij}}
360            """
361            self.__pde_v.setTolerance(self.getSubProblemTolerance())
362            g=self.__g
363            f=self.__f
364            self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
365            if p == None:
366               self.__pde_v.setValue(Y=g)
367            else:
368               self.__pde_v.setValue(Y=g-self.__Q(p))
369            return self.__pde_v.getSolution(verbose=show_details)
370    
371  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
372        """       """
373        solves       solves
374    
375            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
376                  u_{i,i}=0                  u_{i,i}=0
377    
378            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
379            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
380    
381        if surface_stress is not give 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             initialize the Stokes Problem
393    
394             @param domain: domain of the problem. The approximation order needs to be two.
395             @type domain: L{Domain}
396             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
397             """
398           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
399           self.domain=domain           self.domain=domain
400           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
401           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())
402           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
403           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
404                         # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
405    
406           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
407           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
408             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
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            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
417            @param fixed_u_mask: mask of locations with fixed velocity.
418            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
419            @param eta: viscosity
420            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
421            @param surface_stress: normal surface stress
422            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
423            @param stress: initial stress
424        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
425            @note: All values needs to be set.
426    
427        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):          """
428          self.eta=eta          self.eta=eta
429          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
430      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
431          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
432          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
433              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
434              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
435      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
436          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
437            self.__f=f
438            self.__surface_stress=surface_stress
439            self.__stress=stress
440    
441         def inner_pBv(self,p,v):
442             """
443             returns inner product of element p and div(v)
444    
445             @param p: a pressure increment
446             @param v: a residual
447             @return: inner product of element p and div(v)
448             @rtype: C{float}
449             """
450             return util.integrate(-p*util.div(v))
451    
452         def inner_p(self,p0,p1):
453             """
454             Returns inner product of p0 and p1
455    
456        def B(self,arg):           @param p0: a pressure
457           d=util.div(arg)           @param p1: a pressure
458           self.__pde_proj.setValue(Y=d)           @return: inner product of p0 and p1
459           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           @rtype: C{float}
460           return self.__pde_proj.getSolution(verbose=self.show_details)           """
461             s0=util.interpolate(p0/self.eta,Function(self.domain))
462        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))  
463           return util.integrate(s0*s1)           return util.integrate(s0*s1)
464    
465        def getStress(self,u):       def norm_v(self,v):
466           mg=util.grad(u)           """
467           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
468    
469             @param v: a velovity
470             @return: norm of v
471             @rtype: non-negative C{float}
472             """
473             return util.sqrt(util.integrate(util.length(util.grad(v))))
474    
475        def solve_A(self,u,p):       def getV(self, p, v0):
476           """           """
477           solves 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           """           """
483           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
484           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
485           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
486                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
487             else:
488                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
489             out=self.__pde_u.getSolution(verbose=self.show_details)
490             return  out
491    
       def solve_prec(self,p):  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          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 stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
492    
493             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             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             self.__pde_u.setTolerance(self.getSubProblemTolerance())
514             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
515             out=self.__pde_u.getSolution(verbose=self.show_details)
516             return  out
517    
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.1519  
changed lines
  Added in v.2264

  ViewVC Help
Powered by ViewVC 1.1.26