/[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 1562 by gross, Wed May 21 13:04:40 2008 UTC revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 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,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm
38    
39    class DarcyFlow(object):
40        """
41        Represents and solves the problem
42    
43        M{u_i+k_{ij}*p_{,j} = g_i}
44    
45        M{u_{i,i} = f}
46    
47        where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents
48        the permeability.
49    
50        @note: The problem is solved in a least squares formulation.
51        """
52    
53        def __init__(self, domain):
54            """
55            Initializes the Darcy flux problem.
56    
57            @param domain: domain of the problem
58            @type domain: L{Domain}
59            """
60            self.domain=domain
61            self.__pde_v=LinearPDESystem(domain)
62            self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
63            self.__pde_v.setSymmetryOn()
64            self.__pde_p=LinearSinglePDE(domain)
65            self.__pde_p.setSymmetryOn()
66            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
67            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
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: volumetric sources/sinks
74            @type f: scalar value on the domain, e.g. L{Data}
75            @param g: flux sources/sinks
76            @type g: vector value 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 value on the domain (e.g. L{Data})
81            @param permeability: permeability tensor. If scalar C{s} is given the
82                                 tensor with C{s} on the main diagonal is used. If
83                                 vector C{v} is given the tensor with C{v} on the
84                                 main diagonal is used.
85            @type permeability: scalar, vector or tensor values on the domain, e.g.
86                                L{Data}
87    
88            @note: the values of parameters which are not set by calling
89                   C{setValue} are not altered
90            @note: at any point on the boundary of the domain the pressure
91                   (C{location_of_fixed_pressure}) >0 or the normal component of
92                   the flux (C{location_of_fixed_flux[i]}) >0 if the direction of
93                   the normal is along the M{x_i} axis.
94            """
95            if f !=None:
96               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
97               if f.isEmpty():
98                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
99               else:
100                   if f.getRank()>0: raise ValueError,"illegal rank of f."
101               self.f=f
102            if g !=None:
103               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
104               if g.isEmpty():
105                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
106               else:
107                 if not g.getShape()==(self.domain.getDim(),):
108                   raise ValueError,"illegal shape of g"
109               self.__g=g
110    
111            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
112            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
113    
114            if permeability!=None:
115               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
116               if perm.getRank()==0:
117                   perm=perm*util.kronecker(self.domain.getDim())
118               elif perm.getRank()==1:
119                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
120                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
121               elif perm.getRank()==2:
122                  pass
123               else:
124                  raise ValueError,"illegal rank of permeability."
125               self.__permeability=perm
126               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
127    
128    
129        def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
130            """
131            Returns the flux for a given pressure C{p}.
132    
133            The flux is equal to C{fixed_flux} on locations where
134            C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}
135            and C{f} are used.
136    
137            @param p: pressure
138            @type p: scalar value on the domain, e.g. L{Data}
139            @param fixed_flux: flux on the locations of the domain marked by
140                               C{location_of_fixed_flux}
141            @type fixed_flux: vector values on the domain, e.g. L{Data}
142            @param tol: relative tolerance to be used
143            @type tol: positive float
144            @return: flux
145            @rtype: L{Data}
146            @note: the method uses the least squares solution
147                   M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator
148                   and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
149            """
150            self.__pde_v.setTolerance(tol)
151            self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
152            return self.__pde_v.getSolution(verbose=show_details)
153    
154        def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
155            """
156            Solves the problem.
157    
158            The iteration is terminated if the error in the pressure is less than
159            M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand
160            side (see escript user's guide for details).
161    
162            @param u0: initial guess for the flux. At locations in the domain
163                       marked by C{location_of_fixed_flux} the value of C{u0} is
164                       kept unchanged.
165            @type u0: vector value on the domain, e.g. L{Data}
166            @param p0: initial guess for the pressure. At locations in the domain
167                       marked by C{location_of_fixed_pressure} the value of C{p0}
168                       is kept unchanged.
169            @type p0: scalar value on the domain, e.g. L{Data}
170            @param atol: absolute tolerance for the pressure
171            @type atol: non-negative C{float}
172            @param rtol: relative tolerance for the pressure
173            @type rtol: non-negative C{float}
174            @param sub_rtol: tolerance to be used in the sub iteration. It is
175                             recommended that M{sub_rtol<rtol*5.e-3}
176            @type sub_rtol: positive-negative C{float}
177            @param verbose: if True information on iteration progress is printed
178            @type verbose: C{bool}
179            @param show_details: if True information on the sub-iteration process
180                                 is printed
181            @type show_details: C{bool}
182            @return: flux and pressure
183            @rtype: C{tuple} of L{Data}
184    
185            @note: the problem is solved in a least squares formulation:
186    
187            M{(I+D^*D)u+Qp=D^*f+g}
188    
189            M{Q^*u+Q^*Qp=Q^*g}
190    
191            where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the
192            permeability M{k_{ij}}. We eliminate the flux from the problem by
193            setting
194    
195            M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}
196    
197            from the first equation. Inserted into the second equation we get
198    
199            M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}
200            on C{location_of_fixed_pressure}
201    
202            which is solved using the PCG method (precondition is M{Q^*Q}).
203            In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}
204            need to be solved using a sub-iteration scheme.
205            """
206            self.verbose=verbose
207            self.show_details= show_details and self.verbose
208            self.__pde_v.setTolerance(sub_rtol)
209            self.__pde_p.setTolerance(sub_rtol)
210            u2=u0*self.__pde_v.getCoefficient("q")
211            #
212            # first the reference velocity is calculated from
213            #
214            #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
215            #
216            self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
217            u_ref=self.__pde_v.getSolution(verbose=show_details)
218            if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
219            self.__pde_v.setValue(r=Data())
220            #
221            #   and then we calculate a reference pressure
222            #
223            #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
224            #
225            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
226            p_ref=self.__pde_p.getSolution(verbose=self.show_details)
227            if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
228            self.__pde_p.setValue(r=Data())
229            #
230            #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref
231            #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref
232            #
233            #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du
234            #
235            #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
236            #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
237            #
238            #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
239            #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)
240            #
241            #  we use the (du,Qdp) to represent the resudual
242            #  Q^*Q is a preconditioner
243            #
244            #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
245            #
246            Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
247            norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
248            ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
249            if not ATOL>0:
250                raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)
251            if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
252            #
253            #   caclulate the initial residual
254            #
255            self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
256            du=self.__pde_v.getSolution(verbose=show_details)
257            r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
258            dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
259            util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
260            return u_ref+r[1],dp
261    
262        def __Aprod_PCG(self,p):
263            if self.show_details: print "DarcyFlux: Applying operator"
264            Qp=util.tensor_mult(self.__permeability,util.grad(p))
265            self.__pde_v.setValue(Y=Qp,X=Data())
266            w=self.__pde_v.getSolution(verbose=self.show_details)
267            return ArithmeticTuple(-Qp,w)
268    
269        def __inner_PCG(self,p,r):
270            a=util.tensor_mult(self.__permeability,util.grad(p))
271            out=-util.integrate(util.inner(a,r[0]+r[1]))
272            return out
273    
274        def __Msolve_PCG(self,r):
275            if self.show_details: print "DarcyFlux: Applying preconditioner"
276            self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))
277            return self.__pde_p.getSolution(verbose=self.show_details)
278    
279  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
280        """        """
281        solves        Represents and solves the problem
282    
283          M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}
284    
285            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0
                 u_{i,i}=0  
286    
287            u=0 where  fixed_u_mask>0        M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}
           eta*(u_{i,j}+u_{j,i})*n_j=surface_stress  
288    
289        if surface_stress is not give 0 is assumed.        If C{surface_stress} is not given 0 is assumed.
290    
291        typical usage:        Typical usage::
292    
293              sp=StokesProblemCartesian(domain)            sp = StokesProblemCartesian(domain)
294              sp.setTolerance()            sp.setTolerance()
295              sp.initialize(...)            sp.initialize(...)
296              v,p=sp.solve(v0,p0)            v,p = sp.solve(v0,p0)
297        """        """
298        def __init__(self,domain,**kwargs):        def __init__(self,domain,**kwargs):
299             """
300             Initializes the Stokes Problem.
301    
302             @param domain: domain of the problem. The approximation order needs
303                            to be two.
304             @type domain: L{Domain}
305             @warning: The approximation order needs to be two otherwise you may
306                       see oscillations in the pressure.
307             """
308           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
309           self.domain=domain           self.domain=domain
310           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
311           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())
312           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
313           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
314                         # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
315    
316           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
317           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
318             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
319           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
320    
321           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
# Line 74  class StokesProblemCartesian(Homogeneous Line 323  class StokesProblemCartesian(Homogeneous
323           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
324           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
325    
326        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
327          self.eta=eta           """
328          A =self.__pde_u.createCoefficientOfGeneralPDE("A")           Assigns values to the model parameters.
329      self.__pde_u.setValue(A=Data())  
330          for i in range(self.domain.getDim()):           @param f: external force
331          for j in range(self.domain.getDim()):           @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
332              A[i,j,j,i] += 1.           @param fixed_u_mask: mask of locations with fixed velocity
333              A[i,j,i,j] += 1.           @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}
334      self.__pde_prec.setValue(D=1/self.eta)                               or similar
335          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)           @param eta: viscosity
336             @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar
337        def B(self,arg):           @param surface_stress: normal surface stress
338           d=util.div(arg)           @type surface_stress: L{Vector} object on L{FunctionSpace},
339           self.__pde_proj.setValue(Y=d)                                 L{FunctionOnBoundary} or similar
340             @param stress: initial stress
341             @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or
342                           similar
343             @note: All values need to be set.
344             """
345             self.eta=eta
346             A =self.__pde_u.createCoefficient("A")
347             self.__pde_u.setValue(A=Data())
348             for i in range(self.domain.getDim()):
349                 for j in range(self.domain.getDim()):
350                     A[i,j,j,i] += 1.
351                     A[i,j,i,j] += 1.
352             self.__pde_prec.setValue(D=1/self.eta)
353             self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
354             self.__stress=stress
355    
356          def B(self,v):
357             """
358             Returns M{div(v)}.
359             @return: M{div(v)}
360             @rtype: equal to the type of p
361    
362             @note: Boundary conditions on p should be zero!
363             """
364             if self.show_details: print "apply divergence:"
365             self.__pde_proj.setValue(Y=-util.div(v))
366           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_proj.setTolerance(self.getSubProblemTolerance())
367           return self.__pde_proj.getSolution(verbose=self.show_details)           return self.__pde_proj.getSolution(verbose=self.show_details)
368    
369        def inner(self,p0,p1):        def inner_pBv(self,p,Bv):
370           s0=util.interpolate(p0,Function(self.domain))           """
371           s1=util.interpolate(p1,Function(self.domain))           Returns inner product of element p and Bv (overwrite).
372    
373             @type p: equal to the type of p
374             @type Bv: equal to the type of result of operator B
375             @return: inner product of p and Bv
376             @rtype: equal to the type of p
377             """
378             s0=util.interpolate(p,Function(self.domain))
379             s1=util.interpolate(Bv,Function(self.domain))
380           return util.integrate(s0*s1)           return util.integrate(s0*s1)
381    
382        def inner_a(self,a0,a1):        def inner_p(self,p0,p1):
383           p0=util.interpolate(a0[1],Function(self.domain))           """
384           p1=util.interpolate(a1[1],Function(self.domain))           Returns inner product of element p0 and p1 (overwrite).
          alfa=(1/self.vol)*util.integrate(p0)  
          beta=(1/self.vol)*util.integrate(p1)  
      v0=util.grad(a0[0])  
      v1=util.grad(a1[0])  
      #print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1))  
          return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))  
385    
386             @type p0: equal to the type of p
387             @type p1: equal to the type of p
388             @return: inner product of p0 and p1
389             @rtype: equal to the type of p
390             """
391             s0=util.interpolate(p0/self.eta,Function(self.domain))
392             s1=util.interpolate(p1/self.eta,Function(self.domain))
393             return util.integrate(s0*s1)
394    
395        def getStress(self,u):        def inner_v(self,v0,v1):
396           mg=util.grad(u)           """
397           return 2.*self.eta*util.symmetric(mg)           Returns inner product of two elements v0 and v1 (overwrite).
398    
399             @type v0: equal to the type of v
400             @type v1: equal to the type of v
401             @return: inner product of v0 and v1
402             @rtype: equal to the type of v
403             """
404             gv0=util.grad(v0)
405             gv1=util.grad(v1)
406             return util.integrate(util.inner(gv0,gv1))
407    
408        def solve_A(self,u,p):        def solve_A(self,u,p):
409           """           """
410           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).
411           """           """
412             if self.show_details: print "solve for velocity:"
413           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
414           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           if self.__stress.isEmpty():
415           return  self.__pde_u.getSolution(verbose=self.show_details)              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
416             else:
417                self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
418             out=self.__pde_u.getSolution(verbose=self.show_details)
419             return out
420    
421        def solve_prec(self,p):        def solve_prec(self,p):
422       #proj=Projector(domain=self.domain, reduce = True, fast=False)           """
423           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           Applies the preconditioner.
424           self.__pde_prec.setValue(Y=p)           """
425           q=self.__pde_prec.getSolution(verbose=self.show_details)           if self.show_details: print "apply preconditioner:"
          return q  
   
       def solve_prec1(self,p):  
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
426           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
427           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
428           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
      q0=util.interpolate(q,Function(self.domain))  
          print util.inf(q*q0),util.sup(q*q0)  
          q-=(1/self.vol)*util.integrate(q0)  
          print util.inf(q*q0),util.sup(q*q0)  
429           return q           return q
430    
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           print abs(n_v_old-self.__n_v), n_v, self.getTolerance()  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= 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  
   
   

Legend:
Removed from v.1562  
changed lines
  Added in v.2169

  ViewVC Help
Powered by ViewVC 1.1.26