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

Legend:
Removed from v.1551  
changed lines
  Added in v.2156

  ViewVC Help
Powered by ViewVC 1.1.26