/[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 1517 by artak, Fri Apr 18 02:36:37 2008 UTC revision 2251 by gross, Fri Feb 6 06:50:39 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
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.__ATOL= None
68    
69        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
70            """
71            assigns values to model parameters
72    
73            @param f: volumetic sources/sinks
74            @type f: scalar value on the domain (e.g. L{Data})
75            @param g: flux sources/sinks
76            @type g: vector values on the domain (e.g. L{Data})
77            @param location_of_fixed_pressure: mask for locations where pressure is fixed
78            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
79            @param location_of_fixed_flux:  mask for locations where flux is fixed.
80            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
81            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
82                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
83                                 C{v} on the main diagonal is used.
84            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
85    
86            @note: the values of parameters which are not set by calling C{setValue} are not altered.
87            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
88                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
89                   is along the M{x_i} axis.
90            """
91            if f !=None:
92               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
93               if f.isEmpty():
94                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
95               else:
96                   if f.getRank()>0: raise ValueError,"illegal rank of f."
97               self.f=f
98            if g !=None:  
99               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
100               if g.isEmpty():
101                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
102               else:
103                 if not g.getShape()==(self.domain.getDim(),):
104                   raise ValueError,"illegal shape of g"
105               self.__g=g
106    
107            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
108            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
109    
110            if permeability!=None:
111               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
112               if perm.getRank()==0:
113                   perm=perm*util.kronecker(self.domain.getDim())
114               elif perm.getRank()==1:
115                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
116                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
117               elif perm.getRank()==2:
118                  pass
119               else:
120                  raise ValueError,"illegal rank of permeability."
121               self.__permeability=perm
122               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
123    
124    
125        def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):
126            """
127            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
128            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
129            Note that C{g} and C{f} are used, see L{setValue}.
130            
131            @param p: pressure.
132            @type p: scalar value on the domain (e.g. L{Data}).
133            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
134            @type fixed_flux: vector values on the domain (e.g. L{Data}).
135            @param tol: relative tolerance to be used.
136            @type tol: positive C{float}.
137            @return: flux
138            @rtype: L{Data}
139            @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}}
140                   for the permeability M{k_{ij}}
141            """
142            self.__pde_v.setTolerance(tol)
143            g=self.__g
144            f=self.__f
145            self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
146            if p == None:
147               self.__pde_v.setValue(Y=g)
148            else:
149               self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))
150            return self.__pde_v.getSolution(verbose=show_details)
151    
152        def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):
153            """
154            returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}
155            on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).
156            Note that C{g} is used, see L{setValue}.
157            
158            @param v: flux.
159            @type v: vector-valued on the domain (e.g. L{Data}).
160            @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.
161            @type fixed_pressure: vector values on the domain (e.g. L{Data}).
162            @param tol: relative tolerance to be used.
163            @type tol: positive C{float}.
164            @return: pressure
165            @rtype: L{Data}
166            @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}
167                   for the permeability M{k_{ij}}
168            """
169            self.__pde_v.setTolerance(tol)
170            g=self.__g
171            self.__pde_p.setValue(r=fixed_pressure)
172            if v == None:
173               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))
174            else:
175               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))
176            return self.__pde_p.getSolution(verbose=show_details)
177    
178        def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):
179            """
180            set the tolerance C{ATOL} used to terminate the solution process. It is used
181    
182            M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}
183    
184            where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.
185    
186            The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual
187    
188            M{r=Q^*(g-Qp-v)}
189    
190            the condition
191    
192            M{<(Q^*Q)^{-1} r,r> <= ATOL}
193    
194            holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
195    
196            @param atol: absolute tolerance for the pressure
197            @type atol: non-negative C{float}
198            @param rtol: relative tolerance for the pressure
199            @type rtol: non-negative C{float}
200            @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the
201            L{getPressure} method or use  the value from a previous time step.
202            @type p_ref: scalar value on the domain (e.g. L{Data}).
203            @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the
204            L{getFlux} method or use  the value from a previous time step.
205            @type v_ref: vector-valued on the domain (e.g. L{Data}).
206            @return: used absolute tolerance.
207            @rtype: positive C{float}
208            """
209            g=self.__g
210            if not v_ref == None:
211               f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)
212            else:
213               f1=util.integrate(util.length(util.interpolate(g))**2)
214            if not p_ref == None:
215               f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)
216            else:
217               f2=0
218            self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))
219            if self.__ATOL<=0:
220               raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL
221            return self.__ATOL
222            
223        def getTolerance(self):
224            """
225            returns the current tolerance.
226      
227            @return: used absolute tolerance.
228            @rtype: positive C{float}
229            """
230            if self.__ATOL==None:
231               raise ValueError,"no tolerance is defined."
232            return self.__ATOL
233    
234        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
235             """
236             solves the problem.
237    
238             The iteration is terminated if the residual norm is less then self.getTolerance().
239    
240             @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.
241             @type u0: vector value on the domain (e.g. L{Data}).
242             @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.
243             @type p0: scalar value on the domain (e.g. L{Data}).
244             @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
245             @type sub_rtol: positive-negative C{float}
246             @param verbose: if set some information on iteration progress are printed
247             @type verbose: C{bool}
248             @param show_details:  if set information on the subiteration process are printed.
249             @type show_details: C{bool}
250             @return: flux and pressure
251             @rtype: C{tuple} of L{Data}.
252    
253             @note: The problem is solved as a least squares form
254    
255             M{(I+D^*D)u+Qp=D^*f+g}
256             M{Q^*u+Q^*Qp=Q^*g}
257    
258             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
259             We eliminate the flux form the problem by setting
260    
261             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
262    
263             form the first equation. Inserted into the second equation we get
264    
265             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
266            
267             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
268             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
269             """
270             self.verbose=verbose
271             self.show_details= show_details and self.verbose
272             self.__pde_v.setTolerance(sub_rtol)
273             self.__pde_p.setTolerance(sub_rtol)
274             ATOL=self.getTolerance()
275             if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
276             #########################################################################################################################
277             #
278             #   we solve:
279             #  
280             #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )
281             #
282             #   residual is
283             #
284             #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)
285             #
286             #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC
287             #
288             #    we use (g - Qp, v) to represent the residual. not that
289             #
290             #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
291             #
292             #   while the initial residual is
293             #
294             #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC
295             #  
296             d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))
297             self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)
298             v00=self.__pde_v.getSolution(verbose=show_details)
299             if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
300             self.__pde_v.setValue(r=Data())
301             # start CG
302             r=ArithmeticTuple(d0, v00)
303             p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
304             return r[1],p
305    
306        def __Aprod_PCG(self,dp):
307              if self.show_details: print "DarcyFlux: Applying operator"
308              #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)
309              mQdp=util.tensor_mult(self.__permeability,util.grad(dp))
310              self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
311              du=self.__pde_v.getSolution(verbose=self.show_details)
312              return ArithmeticTuple(mQdp,du)
313    
314        def __inner_PCG(self,p,r):
315             a=util.tensor_mult(self.__permeability,util.grad(p))
316             f0=util.integrate(util.inner(a,r[0]))
317             f1=util.integrate(util.inner(a,r[1]))
318             # print "__inner_PCG:",f0,f1,"->",f0-f1
319             return f0-f1
320    
321        def __Msolve_PCG(self,r):
322              if self.show_details: print "DarcyFlux: Applying preconditioner"
323              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())
324              return self.__pde_p.getSolution(verbose=self.show_details)
325    
326  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
327        """       """
328        solves       solves
329    
330            -(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}
331                  u_{i,i}=0                  u_{i,i}=0
332    
333            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
334            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
335    
336        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
337    
338        typical usage:       typical usage:
339    
340              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
341              sp.setTolerance()              sp.setTolerance()
342              sp.initialize(...)              sp.initialize(...)
343              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
344        """       """
345        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
346             """
347             initialize the Stokes Problem
348    
349             @param domain: domain of the problem. The approximation order needs to be two.
350             @type domain: L{Domain}
351             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
352             """
353           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
354           self.domain=domain           self.domain=domain
355           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
356           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())
357           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
358  #         self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
359             # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
360                            
361           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
362           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
363             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
364           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
365    
366           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
367           self.__pde_proj.setReducedOrderOn()          """
368           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
369           self.__pde_proj.setValue(D=1.)  
370            @param f: external force
371        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
372            @param fixed_u_mask: mask of locations with fixed velocity.
373            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
374            @param eta: viscosity
375            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
376            @param surface_stress: normal surface stress
377            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
378            @param stress: initial stress
379        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
380            @note: All values needs to be set.
381    
382            """
383          self.eta=eta          self.eta=eta
384          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
385      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
386          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
387          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
388              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
389              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
390      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
391          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)
392            self.__f=f
393            self.__surface_stress=surface_stress
394            self.__stress=stress
395    
396    #===============================================================================================================
397         def inner_pBv(self,p,v):
398             """
399             returns inner product of element p and div(v)
400    
401             @param p: a pressure increment
402             @param v: a residual
403             @return: inner product of element p and div(v)
404             @rtype: C{float}
405             """
406             return util.integrate(-p*util.div(v))
407    
408        def B(self,arg):       def inner_p(self,p0,p1):
409           d=util.div(arg)           """
410           self.__pde_proj.setValue(Y=d)           Returns inner product of p0 and p1
411           self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
412           return self.__pde_proj.getSolution(verbose=self.show_details)           @param p0: a pressure
413             @param p1: a pressure
414        def inner(self,p0,p1):           @return: inner product of p0 and p1
415           s0=util.interpolate(p0,Function(self.domain))           @rtype: C{float}
416           s1=util.interpolate(p1,Function(self.domain))           """
417             s0=util.interpolate(p0/self.eta,Function(self.domain))
418             s1=util.interpolate(p1/self.eta,Function(self.domain))
419           return util.integrate(s0*s1)           return util.integrate(s0*s1)
420    
421        def getStress(self,u):       def norm_v(self,v):
422           mg=util.grad(u)           """
423           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
424    
425             @param v: a velovity
426             @return: norm of v
427             @rtype: non-negative C{float}
428             """
429             return util.sqrt(util.integrate(util.length(util.grad(v))))
430    
431        def solve_A(self,u,p):       def getV(self, p, v0):
432           """           """
433           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
434    
435             @param p: a pressure
436             @param v0: a initial guess for the value v to return.
437             @return: v given as M{v= A^{-1} (f-B^*p)}
438           """           """
439           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
440           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)
441           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
442                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
443             else:
444                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
445             out=self.__pde_u.getSolution(verbose=self.show_details)
446             return  out
447    
448    
449             raise NotImplementedError,"no v calculation implemented."
450    
451    
452         def norm_Bv(self,v):
453            """
454            Returns Bv (overwrite).
455    
456            @rtype: equal to the type of p
457            @note: boundary conditions on p should be zero!
458            """
459            return util.sqrt(util.integrate(util.div(v)**2))
460    
461        def solve_prec(self,p):       def solve_AinvBt(self,p):
462           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           """
463           self.__pde_prec.setValue(Y=p)           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
          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'):  
       if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*self.getTolerance())  
           self.iter+=1  
           if norm_r <= norm_b*self.getTolerance():  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
464    
465             @param p: a pressure increment
466             @return: the solution of M{Av=B^*p}
467             @note: boundary conditions on v should be zero!
468             """
469             self.__pde_u.setTolerance(self.getSubProblemTolerance())
470             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
471             out=self.__pde_u.getSolution(verbose=self.show_details)
472             return  out
473    
474         def solve_precB(self,v):
475             """
476             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
477             with accuracy L{self.getSubProblemTolerance()} (overwrite).
478    
479             @param v: velocity increment
480             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
481             @note: boundary conditions on p are zero.
482             """
483             self.__pde_prec.setValue(Y=-util.div(v))
484             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
485             return self.__pde_prec.getSolution(verbose=self.show_details)

Legend:
Removed from v.1517  
changed lines
  Added in v.2251

  ViewVC Help
Powered by ViewVC 1.1.26