/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Diff of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1519 by artak, Tue Apr 22 03:45:36 2008 UTC revision 2370 by gross, Mon Apr 6 06:41:49 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               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)
64            if useReduced: self.__pde_v.setReducedOrderOn()
65            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)
69            self.__pde_p.setSymmetryOn()
70            if useReduced: self.__pde_p.setReducedOrderOn()
71            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72            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):
78            """
79            assigns values to model parameters
80    
81            @param f: volumetic sources/sinks
82            @type f: scalar value on the domain (e.g. L{Data})
83            @param g: flux sources/sinks
84            @type g: vector values on the domain (e.g. L{Data})
85            @param location_of_fixed_pressure: mask for locations where pressure is fixed
86            @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.
88            @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
90                                 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.
92            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
93    
94            @note: the values of parameters which are not set by calling C{setValue} are not altered.
95            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
96                   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.
98            """
99            if f !=None:
100               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101               if f.isEmpty():
102                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103               else:
104                   if f.getRank()>0: raise ValueError,"illegal rank of f."
105               self.__f=f
106            if g !=None:
107               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108               if g.isEmpty():
109                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
110               else:
111                 if not g.getShape()==(self.domain.getDim(),):
112                   raise ValueError,"illegal shape of g"
113               self.__g=g
114    
115            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
116            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
117    
118            if permeability!=None:
119               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
120               if perm.getRank()==0:
121                   perm=perm*util.kronecker(self.domain.getDim())
122               elif perm.getRank()==1:
123                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
124                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
125               elif perm.getRank()==2:
126                  pass
127               else:
128                  raise ValueError,"illegal rank of permeability."
129               self.__permeability=perm
130               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            @param rtol: relative tolerance for the pressure
141            @type rtol: non-negative C{float}
142            """
143            if rtol<0:
144                raise ValueError,"Relative tolerance needs to be non-negative."
145            self.__rtol=rtol
146        def getTolerance(self):
147            """
148            returns the relative tolerance
149    
150            @return: current relative tolerance
151            @rtype: C{float}
152            """
153            return self.__rtol
154    
155        def setAbsoluteTolerance(self,atol=0.):
156            """
157            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
158    
159            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 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.
207    
208             The iteration is terminated if the residual norm is less then self.getTolerance().
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.
211             @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.
213             @type p0: scalar value on the domain (e.g. L{Data}).
214             @param verbose: if set some information on iteration progress are printed
215             @type verbose: C{bool}
216             @param show_details:  if set information on the subiteration process are printed.
217             @type show_details: C{bool}
218             @return: flux and pressure
219             @rtype: C{tuple} of L{Data}.
220    
221             @note: The problem is solved as a least squares form
222    
223             M{(I+D^*D)u+Qp=D^*f+g}
224             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}}.
227             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
230    
231             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
234    
235             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.
237             """
238             self.verbose=verbose or True
239             self.show_details= show_details and self.verbose
240             rtol=self.getTolerance()
241             atol=self.getAbsoluteTolerance()
242             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243    
244             num_corrections=0
245             converged=False
246             p=p0
247             norm_r=None
248             while not converged:
249                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250                   Qp=self.__Q(p)
251                   norm_v=self.__L2(v)
252                   norm_Qp=self.__L2(Qp)
253                   if norm_v == 0.:
254                      if norm_Qp == 0.:
255                         return v,p
256                      else:
257                        fac=norm_Qp
258                   else:
259                      if norm_Qp == 0.:
260                        fac=norm_v
261                      else:
262                        fac=2./(1./norm_v+1./norm_Qp)
263                   ATOL=(atol+rtol*fac)
264                   if self.verbose:
265                        print "DarcyFlux: L2 norm of v = %e."%norm_v
266                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
267                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
268                   if norm_r == None or norm_r>ATOL:
269                       if num_corrections>max_num_corrections:
270                             raise ValueError,"maximum number of correction steps reached."
271                       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                   else:
274                       converged=True
275             return v,p
276    #
277    #              
278    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
279    #               #===========================================================================
280    #               norm_r_hat=self.__L2(r_hat)
281    #               norm_v=self.__L2(v)
282    #               norm_g=self.__L2(g)
283    #               norm_gv=self.__L2(g-v)
284    #               norm_Qp=self.__L2(Qp)
285    #               norm_gQp=self.__L2(g-Qp)
286    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
287    #               fac=min(norm_v,norm_Qp,norm_gv)
288    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
289    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
290    #               if r!=None:
291    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
292    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
293    #                   self.setSubProblemTolerance(sub_tol)
294    #                   print "subtol_new=",self.getSubProblemTolerance()
295    #               print "norm_v = ",norm_v
296    #               print "norm_gv = ",norm_gv
297    #               print "norm_Qp = ",norm_Qp
298    #               print "norm_gQp = ",norm_gQp
299    #               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):
342             return util.integrate(util.inner(self.__Q(p), r))
343    
344        def __Msolve_PCG(self,r):
345              self.__pde_p.setTolerance(self.getSubProblemTolerance())
346              if self.show_details: print "DarcyFlux: Applying preconditioner"
347              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
348              # 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            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
383                  u_{i,i}=0                  u_{i,i}=0
384    
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=surface_stress            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 give 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
400    
401             @param domain: domain of the problem. The approximation order needs to be two.
402             @type domain: L{Domain}
403             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
404             """
405           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
406           self.domain=domain           self.domain=domain
407           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
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(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
411                         # 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)
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()):
419           self.__pde_proj.setReducedOrderOn()          """
420           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
421           self.__pde_proj.setValue(D=1.)  
422            @param f: external force
423            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
424            @param fixed_u_mask: mask of locations with fixed velocity.
425            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
426            @param eta: viscosity
427            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
428            @param surface_stress: normal surface stress
429            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
430            @param stress: initial stress
431        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
432            @note: All values needs to be set.
433    
434        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):          """
435          self.eta=eta          self.eta=eta
436          A =self.__pde_u.createCoefficientOfGeneralPDE("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
447    
448         def inner_pBv(self,p,v):
449             """
450             returns inner product of element p and div(v)
451    
452             @param p: a pressure increment
453             @param v: a residual
454             @return: inner product of element p and div(v)
455             @rtype: C{float}
456             """
457             return util.integrate(-p*util.div(v))
458    
459         def inner_p(self,p0,p1):
460             """
461             Returns inner product of p0 and p1
462    
463        def B(self,arg):           @param p0: a pressure
464           d=util.div(arg)           @param p1: a pressure
465           self.__pde_proj.setValue(Y=d)           @return: inner product of p0 and p1
466           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           @rtype: C{float}
467           return self.__pde_proj.getSolution(verbose=self.show_details)           """
468             s0=util.interpolate(p0/self.eta,Function(self.domain))
469        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))  
470           return util.integrate(s0*s1)           return util.integrate(s0*s1)
471    
472        def getStress(self,u):       def norm_v(self,v):
473           mg=util.grad(u)           """
474           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
475    
476             @param v: a velovity
477             @return: norm of v
478             @rtype: non-negative C{float}
479             """
480             return util.sqrt(util.integrate(util.length(util.grad(v))))
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           """           """
490           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
491           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)
492           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
493                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
494             else:
495                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    
       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  
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)
523             return  out
524    
525         def solve_precB(self,v):
526             """
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())
536             return self.__pde_prec.getSolution(verbose=self.show_details)

Legend:
Removed from v.1519  
changed lines
  Added in v.2370

  ViewVC Help
Powered by ViewVC 1.1.26