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

revision 1562 by gross, Wed May 21 13:04:40 2008 UTC revision 2208 by gross, Mon Jan 12 06:37:07 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
#
#######################################################
10  #  #
11    ########################################################
12
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
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"
http://www.access.edu.au
__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,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:
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:
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             #
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)
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):
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
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
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
# Line 58  class StokesProblemCartesian(Homogeneous Line 343  class StokesProblemCartesian(Homogeneous
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             """
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)           self.__pde_proj=LinearPDE(domain)
# Line 74  class StokesProblemCartesian(Homogeneous Line 368  class StokesProblemCartesian(Homogeneous
368           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
369           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
370
372            """
373            assigns values to the model parameters
374
375            @param f: external force
376            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
378            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
379            @param eta: viscosity
380            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
381            @param surface_stress: normal surface stress
382            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
383            @param stress: initial stress
384        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
385            @note: All values needs to be set.
386
387            """
388          self.eta=eta          self.eta=eta
389          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
390      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
391          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
392          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
# Line 84  class StokesProblemCartesian(Homogeneous Line 394  class StokesProblemCartesian(Homogeneous
394              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
395      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
397            self.__stress=stress
398
399        def B(self,arg):        def B(self,v):
400           d=util.div(arg)          """
401           self.__pde_proj.setValue(Y=d)          returns div(v)
402           self.__pde_proj.setTolerance(self.getSubProblemTolerance())          @rtype: equal to the type of p
403           return self.__pde_proj.getSolution(verbose=self.show_details)
404            @note: boundary conditions on p should be zero!
405        def inner(self,p0,p1):          """
406           s0=util.interpolate(p0,Function(self.domain))          if self.show_details: print "apply divergence:"
407           s1=util.interpolate(p1,Function(self.domain))          self.__pde_proj.setValue(Y=-util.div(v))
408            self.__pde_proj.setTolerance(self.getSubProblemTolerance())
409            return self.__pde_proj.getSolution(verbose=self.show_details)
410
411          def inner_pBv(self,p,Bv):
412             """
413             returns inner product of element p and Bv  (overwrite)
414
415             @type p: equal to the type of p
416             @type Bv: equal to the type of result of operator B
417             @rtype: C{float}
418
419             @rtype: equal to the type of p
420             """
421             s0=util.interpolate(p,Function(self.domain))
422             s1=util.interpolate(Bv,Function(self.domain))
423           return util.integrate(s0*s1)           return util.integrate(s0*s1)
424
425        def inner_a(self,a0,a1):        def inner_p(self,p0,p1):
426           p0=util.interpolate(a0[1],Function(self.domain))           """
427           p1=util.interpolate(a1[1],Function(self.domain))           returns inner product of element p0 and p1  (overwrite)
428           alfa=(1/self.vol)*util.integrate(p0)
429           beta=(1/self.vol)*util.integrate(p1)           @type p0: equal to the type of p
430       v0=util.grad(a0[0])           @type p1: equal to the type of p
#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))
432
433             @rtype: equal to the type of p
434             """
435             s0=util.interpolate(p0/self.eta,Function(self.domain))
436             s1=util.interpolate(p1/self.eta,Function(self.domain))
437             return util.integrate(s0*s1)
438
439        def getStress(self,u):        def inner_v(self,v0,v1):
441           return 2.*self.eta*util.symmetric(mg)           returns inner product of two element v0 and v1  (overwrite)
442
443             @type v0: equal to the type of v
444             @type v1: equal to the type of v
445             @rtype: C{float}
446
447             @rtype: equal to the type of v
448             """
451             return util.integrate(util.inner(gv0,gv1))
452
453        def solve_A(self,u,p):        def solve_A(self,u,p):
454           """           """
456           """           """
457             if self.show_details: print "solve for velocity:"
458           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
459           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           if self.__stress.isEmpty():
461             else:
463             out=self.__pde_u.getSolution(verbose=self.show_details)
464             return  out
465
466        def solve_prec(self,p):        def solve_prec(self,p):
467       #proj=Projector(domain=self.domain, reduce = True, fast=False)           if self.show_details: print "apply preconditioner:"
self.__pde_prec.setTolerance(self.getSubProblemTolerance())
self.__pde_prec.setValue(Y=p)
q=self.__pde_prec.getSolution(verbose=self.show_details)
return q

def solve_prec1(self,p):
#proj=Projector(domain=self.domain, reduce = True, fast=False)
468           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
469           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
470           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)
471           return q           return q

def stoppingcriterium(self,Bv,v,p):
n_r=util.sqrt(self.inner(Bv,Bv))
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.2208