/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Annotation of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2100 - (hide annotations)
Wed Nov 26 08:13:00 2008 UTC (11 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 62936 byte(s)
This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module. 
Some documentation for both classes has been added.
The convection code is only linear at the moment.



1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 ksteube 1312 #
12 ksteube 1809 ########################################################
13 jgs 121
14 ksteube 1809 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 jgs 121 """
23 jgs 149 Provides some tools related to PDEs.
24 jgs 121
25 jgs 149 Currently includes:
26     - Projector - to project a discontinuous
27 gross 351 - Locator - to trace values in data objects at a certain location
28     - TimeIntegrationManager - to handel extraplotion in time
29 gross 867 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30 gross 637
31     @var __author__: name of author
32     @var __copyright__: copyrights
33     @var __license__: licence agreement
34     @var __url__: url entry point on documentation
35     @var __version__: version
36     @var __date__: date of the version
37 jgs 121 """
38    
39 gross 637 __author__="Lutz Gross, l.gross@uq.edu.au"
40 elspeth 609
41 gross 637
42 jgs 149 import escript
43     import linearPDEs
44 jgs 121 import numarray
45 jgs 149 import util
46 ksteube 1312 import math
47 jgs 121
48 artak 1465 ##### Added by Artak
49 gross 1467 # from Numeric import zeros,Int,Float64
50 artak 1465 ###################################
51    
52    
53 gross 351 class TimeIntegrationManager:
54     """
55     a simple mechanism to manage time dependend values.
56    
57 gross 720 typical usage is::
58 gross 351
59 gross 720 dt=0.1 # time increment
60     tm=TimeIntegrationManager(inital_value,p=1)
61     while t<1.
62     v_guess=tm.extrapolate(dt) # extrapolate to t+dt
63     v=...
64     tm.checkin(dt,v)
65     t+=dt
66 gross 351
67 gross 720 @note: currently only p=1 is supported.
68 gross 351 """
69     def __init__(self,*inital_values,**kwargs):
70     """
71     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
72     """
73     if kwargs.has_key("p"):
74     self.__p=kwargs["p"]
75     else:
76     self.__p=1
77     if kwargs.has_key("time"):
78     self.__t=kwargs["time"]
79     else:
80     self.__t=0.
81     self.__v_mem=[inital_values]
82     self.__order=0
83     self.__dt_mem=[]
84     self.__num_val=len(inital_values)
85    
86     def getTime(self):
87     return self.__t
88 gross 396 def getValue(self):
89 gross 409 out=self.__v_mem[0]
90     if len(out)==1:
91     return out[0]
92     else:
93     return out
94    
95 gross 351 def checkin(self,dt,*values):
96     """
97     adds new values to the manager. the p+1 last value get lost
98     """
99     o=min(self.__order+1,self.__p)
100     self.__order=min(self.__order+1,self.__p)
101     v_mem_new=[values]
102     dt_mem_new=[dt]
103     for i in range(o-1):
104     v_mem_new.append(self.__v_mem[i])
105     dt_mem_new.append(self.__dt_mem[i])
106     v_mem_new.append(self.__v_mem[o-1])
107     self.__order=o
108     self.__v_mem=v_mem_new
109     self.__dt_mem=dt_mem_new
110     self.__t+=dt
111    
112     def extrapolate(self,dt):
113     """
114     extrapolates to dt forward in time.
115     """
116     if self.__order==0:
117     out=self.__v_mem[0]
118     else:
119     out=[]
120     for i in range(self.__num_val):
121     out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
122    
123     if len(out)==0:
124     return None
125     elif len(out)==1:
126     return out[0]
127     else:
128     return out
129 gross 396
130 gross 351
131 jgs 121 class Projector:
132 jgs 149 """
133     The Projector is a factory which projects a discontiuous function onto a
134     continuous function on the a given domain.
135     """
136 jgs 121 def __init__(self, domain, reduce = True, fast=True):
137     """
138 jgs 149 Create a continuous function space projector for a domain.
139 jgs 121
140 jgs 149 @param domain: Domain of the projection.
141     @param reduce: Flag to reduce projection order (default is True)
142     @param fast: Flag to use a fast method based on matrix lumping (default is true)
143 jgs 121 """
144 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
145 jgs 148 if fast:
146 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
147 jgs 121 self.__pde.setSymmetryOn()
148     self.__pde.setReducedOrderTo(reduce)
149     self.__pde.setValue(D = 1.)
150 ksteube 1312 return
151 jgs 121
152     def __call__(self, input_data):
153     """
154 jgs 149 Projects input_data onto a continuous function
155 jgs 121
156 jgs 149 @param input_data: The input_data to be projected.
157 jgs 121 """
158 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
159 gross 1122 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
160 jgs 121 if input_data.getRank()==0:
161     self.__pde.setValue(Y = input_data)
162     out=self.__pde.getSolution()
163     elif input_data.getRank()==1:
164     for i0 in range(input_data.getShape()[0]):
165     self.__pde.setValue(Y = input_data[i0])
166     out[i0]=self.__pde.getSolution()
167     elif input_data.getRank()==2:
168     for i0 in range(input_data.getShape()[0]):
169     for i1 in range(input_data.getShape()[1]):
170     self.__pde.setValue(Y = input_data[i0,i1])
171     out[i0,i1]=self.__pde.getSolution()
172     elif input_data.getRank()==3:
173     for i0 in range(input_data.getShape()[0]):
174     for i1 in range(input_data.getShape()[1]):
175     for i2 in range(input_data.getShape()[2]):
176     self.__pde.setValue(Y = input_data[i0,i1,i2])
177     out[i0,i1,i2]=self.__pde.getSolution()
178     else:
179     for i0 in range(input_data.getShape()[0]):
180     for i1 in range(input_data.getShape()[1]):
181     for i2 in range(input_data.getShape()[2]):
182     for i3 in range(input_data.getShape()[3]):
183     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
184     out[i0,i1,i2,i3]=self.__pde.getSolution()
185     return out
186    
187 gross 525 class NoPDE:
188     """
189     solves the following problem for u:
190 jgs 121
191 gross 525 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
192    
193     with constraint
194    
195     M{u[j]=r[j]} where M{q[j]>0}
196    
197     where D, Y, r and q are given functions of rank 1.
198    
199     In the case of scalars this takes the form
200    
201     M{D*u=Y}
202    
203     with constraint
204    
205     M{u=r} where M{q>0}
206    
207     where D, Y, r and q are given scalar functions.
208    
209     The constraint is overwriting any other condition.
210    
211 gross 720 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
212     that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
213     thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
214 gross 525 """
215     def __init__(self,domain,D=None,Y=None,q=None,r=None):
216     """
217     initialize the problem
218    
219     @param domain: domain of the PDE.
220     @type domain: L{Domain}
221     @param D: coefficient of the solution.
222 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
223 gross 525 @param Y: right hand side
224 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
225 gross 525 @param q: location of constraints
226 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
227 gross 525 @param r: value of solution at locations of constraints
228 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
229 gross 525 """
230     self.__domain=domain
231     self.__D=D
232     self.__Y=Y
233     self.__q=q
234     self.__r=r
235     self.__u=None
236     self.__function_space=escript.Solution(self.__domain)
237     def setReducedOn(self):
238     """
239     sets the L{FunctionSpace} of the solution to L{ReducedSolution}
240     """
241     self.__function_space=escript.ReducedSolution(self.__domain)
242     self.__u=None
243    
244     def setReducedOff(self):
245     """
246     sets the L{FunctionSpace} of the solution to L{Solution}
247     """
248     self.__function_space=escript.Solution(self.__domain)
249     self.__u=None
250    
251     def setValue(self,D=None,Y=None,q=None,r=None):
252     """
253     assigns values to the parameters.
254    
255     @param D: coefficient of the solution.
256 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
257 gross 525 @param Y: right hand side
258 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
259 gross 525 @param q: location of constraints
260 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
261 gross 525 @param r: value of solution at locations of constraints
262 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
263 gross 525 """
264     if not D==None:
265     self.__D=D
266     self.__u=None
267     if not Y==None:
268     self.__Y=Y
269     self.__u=None
270     if not q==None:
271     self.__q=q
272     self.__u=None
273     if not r==None:
274     self.__r=r
275     self.__u=None
276    
277     def getSolution(self):
278     """
279     returns the solution
280    
281     @return: the solution of the problem
282     @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
283     """
284     if self.__u==None:
285     if self.__D==None:
286     raise ValueError,"coefficient D is undefined"
287     D=escript.Data(self.__D,self.__function_space)
288     if D.getRank()>1:
289     raise ValueError,"coefficient D must have rank 0 or 1"
290     if self.__Y==None:
291     self.__u=escript.Data(0.,D.getShape(),self.__function_space)
292     else:
293     self.__u=util.quotient(self.__Y,D)
294     if not self.__q==None:
295     q=util.wherePositive(escript.Data(self.__q,self.__function_space))
296     self.__u*=(1.-q)
297     if not self.__r==None: self.__u+=q*self.__r
298     return self.__u
299    
300 jgs 147 class Locator:
301     """
302 jgs 149 Locator provides access to the values of data objects at a given
303     spatial coordinate x.
304    
305     In fact, a Locator object finds the sample in the set of samples of a
306     given function space or domain where which is closest to the given
307     point x.
308 jgs 147 """
309    
310     def __init__(self,where,x=numarray.zeros((3,))):
311 jgs 149 """
312     Initializes a Locator to access values in Data objects on the Doamin
313     or FunctionSpace where for the sample point which
314     closest to the given point x.
315 gross 880
316     @param where: function space
317     @type where: L{escript.FunctionSpace}
318     @param x: coefficient of the solution.
319     @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
320 jgs 149 """
321     if isinstance(where,escript.FunctionSpace):
322 jgs 147 self.__function_space=where
323 jgs 121 else:
324 jgs 149 self.__function_space=escript.ContinuousFunction(where)
325 gross 880 if isinstance(x, list):
326     self.__id=[]
327     for p in x:
328 gross 921 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
329 gross 880 else:
330 gross 921 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
331 jgs 121
332 jgs 147 def __str__(self):
333 jgs 149 """
334     Returns the coordinates of the Locator as a string.
335     """
336 gross 880 x=self.getX()
337     if instance(x,list):
338     out="["
339     first=True
340     for xx in x:
341     if not first:
342     out+=","
343     else:
344     first=False
345     out+=str(xx)
346     out+="]>"
347     else:
348     out=str(x)
349     return out
350 jgs 121
351 gross 880 def getX(self):
352     """
353     Returns the exact coordinates of the Locator.
354     """
355     return self(self.getFunctionSpace().getX())
356    
357 jgs 147 def getFunctionSpace(self):
358 jgs 149 """
359     Returns the function space of the Locator.
360     """
361 jgs 147 return self.__function_space
362    
363 gross 880 def getId(self,item=None):
364 jgs 149 """
365     Returns the identifier of the location.
366     """
367 gross 880 if item == None:
368     return self.__id
369     else:
370     if isinstance(self.__id,list):
371     return self.__id[item]
372     else:
373     return self.__id
374 jgs 121
375    
376 jgs 147 def __call__(self,data):
377 jgs 149 """
378     Returns the value of data at the Locator of a Data object otherwise
379     the object is returned.
380     """
381 jgs 147 return self.getValue(data)
382 jgs 121
383 jgs 147 def getValue(self,data):
384 jgs 149 """
385     Returns the value of data at the Locator if data is a Data object
386     otherwise the object is returned.
387     """
388     if isinstance(data,escript.Data):
389 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
390 gross 880 dat=data
391 jgs 147 else:
392 gross 880 dat=data.interpolate(self.getFunctionSpace())
393     id=self.getId()
394     r=data.getRank()
395     if isinstance(id,list):
396     out=[]
397     for i in id:
398 gross 921 o=data.getValueOfGlobalDataPoint(*i)
399 gross 880 if data.getRank()==0:
400     out.append(o[0])
401     else:
402     out.append(o)
403     return out
404 jgs 147 else:
405 gross 921 out=data.getValueOfGlobalDataPoint(*id)
406 gross 880 if data.getRank()==0:
407     return out[0]
408     else:
409     return out
410 jgs 147 else:
411     return data
412 jgs 149
413 ksteube 1312 class SolverSchemeException(Exception):
414     """
415     exceptions thrown by solvers
416     """
417     pass
418    
419     class IndefinitePreconditioner(SolverSchemeException):
420     """
421     the preconditioner is not positive definite.
422     """
423     pass
424     class MaxIterReached(SolverSchemeException):
425     """
426     maxium number of iteration steps is reached.
427     """
428     pass
429 gross 2100 class CorrectionFailed(SolverSchemeException):
430     """
431     no convergence has been achieved in the solution correction scheme.
432     """
433     pass
434 ksteube 1312 class IterationBreakDown(SolverSchemeException):
435     """
436     iteration scheme econouters an incurable breakdown.
437     """
438     pass
439     class NegativeNorm(SolverSchemeException):
440     """
441     a norm calculation returns a negative norm.
442     """
443     pass
444    
445 gross 2100 def PCG(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, initial_guess=True, verbose=False):
446 ksteube 1312 """
447     Solver for
448    
449     M{Ax=b}
450    
451     with a symmetric and positive definite operator A (more details required!).
452     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
453    
454 gross 1330 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
455 ksteube 1312
456     For details on the preconditioned conjugate gradient method see the book:
457    
458     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
459     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
460     C. Romine, and H. van der Vorst.
461    
462     @param b: the right hand side of the liner system. C{b} is altered.
463 gross 1330 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
464 ksteube 1312 @param Aprod: returns the value Ax
465 jfenwick 2061 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
466     The returned object needs to be of the same type like argument C{b}.
467     @param Msolve: solves Mx=r
468     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}.
469     The returned object needs to be of the same type like argument C{x}.
470     @param bilinearform: inner product C{<x,r>}
471     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
472     The returned value is a C{float}.
473     @param stoppingcriterium: function which returns True if a stopping criterium is meet.
474     C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current
475     norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and
476     the current solution approximation. C{stoppingcriterium} is called in each iteration step.
477 gross 1330 @type stoppingcriterium: function that returns C{True} or C{False}
478     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
479     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
480 ksteube 1312 @param iter_max: maximum number of iteration steps.
481     @type iter_max: C{int}
482 gross 1330 @return: the solution approximation and the corresponding residual
483     @rtype: C{tuple}
484     @warning: C{b} and C{x} are altered.
485 ksteube 1312 """
486     iter=0
487 gross 1330 if x==None:
488     x=0*b
489 gross 2100 r=b*1.
490     initial_guess=False
491     if initial_guess:
492     r=b+(-1)*Aprod(x)
493     x=x*1.
494 gross 1330 else:
495 gross 2100 x=x*0
496     r=b*1
497 ksteube 1312 rhat=Msolve(r)
498 gross 1330 d = rhat
499 ksteube 1312 rhat_dot_r = bilinearform(rhat, r)
500 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
501 gross 2100 norm_b=math.sqrt(rhat_dot_r)
502     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_b, atol+rtol*norm_b)
503 ksteube 1312
504 gross 2100 while not math.sqrt(rhat_dot_r) <= atol+rtol*norm_b:
505 gross 1330 iter+=1
506 ksteube 1312 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
507    
508     q=Aprod(d)
509     alpha = rhat_dot_r / bilinearform(d, q)
510     x += alpha * d
511     r += (-alpha) * q
512    
513     rhat=Msolve(r)
514     rhat_dot_r_new = bilinearform(rhat, r)
515     beta = rhat_dot_r_new / rhat_dot_r
516     rhat+=beta * d
517     d=rhat
518    
519     rhat_dot_r = rhat_dot_r_new
520 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521 gross 2100 if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
522     if verbose: print "PCG: tolerance reached after %s steps."%iter
523 gross 1330 return x,r
524 ksteube 1312
525 gross 1878 class Defect(object):
526     """
527     defines a non-linear defect F(x) of a variable x
528     """
529     def __init__(self):
530     """
531     initialize defect
532     """
533     self.setDerivativeIncrementLength()
534 artak 1465
535 gross 1878 def bilinearform(self, x0, x1):
536     """
537     returns the inner product of x0 and x1
538     @param x0: a value for x
539     @param x1: a value for x
540     @return: the inner product of x0 and x1
541     @rtype: C{float}
542     """
543     return 0
544    
545     def norm(self,x):
546     """
547     the norm of argument C{x}
548    
549     @param x: a value for x
550     @return: norm of argument x
551     @rtype: C{float}
552     @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
553     """
554     s=self.bilinearform(x,x)
555     if s<0: raise NegativeNorm,"negative norm."
556     return math.sqrt(s)
557 artak 1465
558 gross 1878
559     def eval(self,x):
560     """
561     returns the value F of a given x
562    
563     @param x: value for which the defect C{F} is evalulated.
564     @return: value of the defect at C{x}
565     """
566     return 0
567    
568     def __call__(self,x):
569     return self.eval(x)
570    
571     def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
572     """
573     sets the relative length of the increment used to approximate the derivative of the defect
574     the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
575    
576     @param inc: relative increment length
577     @type inc: positive C{float}
578     """
579     if inc<=0: raise ValueError,"positive increment required."
580     self.__inc=inc
581    
582     def getDerivativeIncrementLength(self):
583     """
584     returns the relative increment length used to approximate the derivative of the defect
585     @return: value of the defect at C{x}
586     @rtype: positive C{float}
587     """
588     return self.__inc
589    
590     def derivative(self, F0, x0, v, v_is_normalised=True):
591     """
592     returns the directional derivative at x0 in the direction of v
593    
594     @param F0: value of this defect at x0
595     @param x0: value at which derivative is calculated.
596     @param v: direction
597     @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
598     @return: derivative of this defect at x0 in the direction of C{v}
599     @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
600     maybe oepsnew verwritten to use exact evalution.
601     """
602     normx=self.norm(x0)
603     if normx>0:
604     epsnew = self.getDerivativeIncrementLength() * normx
605     else:
606     epsnew = self.getDerivativeIncrementLength()
607     if not v_is_normalised:
608     normv=self.norm(v)
609     if normv<=0:
610     return F0*0
611     else:
612     epsnew /= normv
613     F1=self.eval(x0 + epsnew * v)
614     return (F1-F0)/epsnew
615    
616     ######################################
617     def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
618     """
619     solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
620    
621     M{norm(F(x) <= atol + rtol * norm(F(x0)}
622    
623     where M{x0} is the initial guess.
624    
625     @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
626     @type defect: L{Defect}
627     @param x: initial guess for the solution, C{x} is altered.
628     @type x: any object type allowing basic operations such as L{numarray.NumArray}, L{Data}
629     @param iter_max: maximum number of iteration steps
630     @type iter_max: positive C{int}
631     @param sub_iter_max:
632     @type sub_iter_max:
633     @param atol: absolute tolerance for the solution
634     @type atol: positive C{float}
635     @param rtol: relative tolerance for the solution
636     @type rtol: positive C{float}
637     @param gamma: tolerance safety factor for inner iteration
638     @type gamma: positive C{float}, less than 1
639     @param sub_tol_max: upper bound for inner tolerance.
640     @type sub_tol_max: positive C{float}, less than 1
641     @return: an approximation of the solution with the desired accuracy
642     @rtype: same type as the initial guess.
643     """
644     lmaxit=iter_max
645     if atol<0: raise ValueError,"atol needs to be non-negative."
646     if rtol<0: raise ValueError,"rtol needs to be non-negative."
647     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
648     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
649     if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max
650    
651     F=defect(x)
652     fnrm=defect.norm(F)
653     stop_tol=atol + rtol*fnrm
654     sub_tol=sub_tol_max
655     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
656     if verbose: print " tolerance = %e."%sub_tol
657     iter=1
658     #
659     # main iteration loop
660     #
661     while not fnrm<=stop_tol:
662     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
663     #
664     # adjust sub_tol_
665     #
666     if iter > 1:
667     rat=fnrm/fnrmo
668     sub_tol_old=sub_tol
669     sub_tol=gamma*rat**2
670     if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
671     sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
672     #
673     # calculate newton increment xc
674     # if iter_max in __FDGMRES is reached MaxIterReached is thrown
675     # if iter_restart -1 is returned as sub_iter
676     # if atol is reached sub_iter returns the numer of steps performed to get there
677     #
678     #
679     if verbose: print " subiteration (GMRES) is called with relative tolerance %e."%sub_tol
680     try:
681     xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
682     except MaxIterReached:
683     raise MaxIterReached,"maximum number of %s steps reached."%iter_max
684     if sub_iter<0:
685     iter+=sub_iter_max
686     else:
687     iter+=sub_iter
688     # ====
689     x+=xc
690     F=defect(x)
691     iter+=1
692     fnrmo, fnrm=fnrm, defect.norm(F)
693     if verbose: print " step %s: residual %e."%(iter,fnrm)
694     if verbose: print "NewtonGMRES: completed after %s steps."%iter
695     return x
696    
697     def __givapp(c,s,vin):
698     """
699     apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
700     @warning: C{vin} is altered.
701     """
702     vrot=vin
703 gross 1467 if isinstance(c,float):
704     vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
705     else:
706     for i in range(len(c)):
707     w1=c[i]*vrot[i]-s[i]*vrot[i+1]
708     w2=s[i]*vrot[i]+c[i]*vrot[i+1]
709     vrot[i:i+2]=w1,w2
710 artak 1465 return vrot
711    
712 gross 1878 def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
713 artak 1475 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
714     c=numarray.zeros(iter_restart,numarray.Float64)
715     s=numarray.zeros(iter_restart,numarray.Float64)
716     g=numarray.zeros(iter_restart,numarray.Float64)
717 artak 1465 v=[]
718    
719 gross 1878 rho=defect.norm(F0)
720     if rho<=0.: return x0*0
721 artak 1488
722 gross 1878 v.append(-F0/rho)
723 artak 1465 g[0]=rho
724 gross 1878 iter=0
725     while rho > atol and iter<iter_restart-1:
726 artak 1557
727 artak 1465 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
728    
729 gross 1878 p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
730 artak 1465 v.append(p)
731    
732 gross 1878 v_norm1=defect.norm(v[iter+1])
733 artak 1465
734 gross 1878 # Modified Gram-Schmidt
735 artak 1557 for j in range(iter+1):
736 gross 2100 h[j,iter]=defect.bilinearform(v[j],v[iter+1])
737     v[iter+1]-=h[j,iter]*v[j]
738 artak 1465
739 gross 2100 h[iter+1,iter]=defect.norm(v[iter+1])
740     v_norm2=h[iter+1,iter]
741 artak 1465
742 gross 1878 # Reorthogonalize if needed
743 artak 1465 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
744 gross 1878 for j in range(iter+1):
745     hr=defect.bilinearform(v[j],v[iter+1])
746 gross 2100 h[j,iter]=h[j,iter]+hr
747 gross 1878 v[iter+1] -= hr*v[j]
748 artak 1465
749 gross 1878 v_norm2=defect.norm(v[iter+1])
750 gross 2100 h[iter+1,iter]=v_norm2
751 gross 1878 # watch out for happy breakdown
752 artak 1550 if not v_norm2 == 0:
753 gross 2100 v[iter+1]=v[iter+1]/h[iter+1,iter]
754 artak 1465
755 gross 1878 # Form and store the information for the new Givens rotation
756 artak 1465 if iter > 0 :
757 artak 1550 hhat=numarray.zeros(iter+1,numarray.Float64)
758 gross 2100 for i in range(iter+1) : hhat[i]=h[i,iter]
759 gross 1878 hhat=__givapp(c[0:iter],s[0:iter],hhat);
760 gross 2100 for i in range(iter+1) : h[i,iter]=hhat[i]
761 artak 1465
762 gross 2100 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
763 artak 1557
764 artak 1465 if mu!=0 :
765 gross 2100 c[iter]=h[iter,iter]/mu
766     s[iter]=-h[iter+1,iter]/mu
767     h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
768     h[iter+1,iter]=0.0
769 gross 1878 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
770 artak 1465
771 gross 1878 # Update the residual norm
772 artak 1465 rho=abs(g[iter+1])
773     iter+=1
774    
775 gross 1878 # At this point either iter > iter_max or rho < tol.
776     # It's time to compute x and leave.
777 artak 1465 if iter > 0 :
778 gross 1467 y=numarray.zeros(iter,numarray.Float64)
779 gross 2100 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
780 artak 1465 if iter > 1 :
781     i=iter-2
782     while i>=0 :
783 gross 2100 y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
784 artak 1465 i=i-1
785     xhat=v[iter-1]*y[iter-1]
786     for i in range(iter-1):
787     xhat += v[i]*y[i]
788 gross 1878 else :
789     xhat=v[0] * 0
790 artak 1557
791 artak 1488 if iter<iter_restart-1:
792 gross 1878 stopped=iter
793 artak 1475 else:
794 gross 1878 stopped=-1
795 artak 1465
796 gross 1878 return xhat,stopped
797 artak 1481
798 gross 1878 ##############################################
799 gross 2100 def GMRES(b, Aprod, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, iter_restart=20, verbose=False):
800 gross 1878 ################################################
801     m=iter_restart
802     iter=0
803 gross 2100 if x == None:
804     xc=b*0
805     else:
806     xc=x*1
807     if rtol>0:
808     r_dot_r = bilinearform(b, b)
809     if r_dot_r<0: raise NegativeNorm,"negative norm."
810     atol2=atol+rtol*math.sqrt(r_dot_r)
811     if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
812     else:
813     atol2=atol
814     if verbose: print "GMRES: absolute tolerance = %e"%atol2
815    
816 gross 1878 while True:
817     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
818 gross 2100 xc,stopped=__GMRESm(b, Aprod, bilinearform, atol2, x=xc, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
819     iter+=iter_restart
820 gross 1878 if stopped: break
821 gross 2100 if verbose: print "GMRES: restart."
822     if verbose: print "GMRES: tolerance has reached."
823 gross 1878 return xc
824 artak 1550
825 gross 2100 def __GMRESm(b, Aprod, bilinearform, atol, x, iter_max=100, iter_restart=20, verbose=False):
826 gross 1878 iter=0
827 artak 1519
828 gross 2100 h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
829 artak 1519 c=numarray.zeros(iter_restart,numarray.Float64)
830     s=numarray.zeros(iter_restart,numarray.Float64)
831 gross 2100 g=numarray.zeros(iter_restart+1,numarray.Float64)
832 artak 1519 v=[]
833    
834 gross 2100 r=b-Aprod(x)
835     r_dot_r = bilinearform(r, r)
836     if r_dot_r<0: raise NegativeNorm,"negative norm."
837 artak 1519 rho=math.sqrt(r_dot_r)
838    
839     v.append(r/rho)
840     g[0]=rho
841    
842 gross 2100 if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
843     while not (rho<=atol or iter==iter_restart):
844 artak 1519
845     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
846    
847 gross 2100 p=Aprod(v[iter])
848 artak 1519 v.append(p)
849    
850     v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
851    
852     # Modified Gram-Schmidt
853 gross 1878 for j in range(iter+1):
854 gross 2100 h[j,iter]=bilinearform(v[j],v[iter+1])
855     v[iter+1]-=h[j,iter]*v[j]
856 artak 1519
857 gross 2100 h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
858     v_norm2=h[iter+1,iter]
859 artak 1519
860     # Reorthogonalize if needed
861     if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
862 gross 1878 for j in range(iter+1):
863 artak 1519 hr=bilinearform(v[j],v[iter+1])
864 gross 2100 h[j,iter]=h[j,iter]+hr
865 gross 1878 v[iter+1] -= hr*v[j]
866 artak 1519
867     v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
868 gross 2100 h[iter+1,iter]=v_norm2
869 artak 1519
870     # watch out for happy breakdown
871 gross 1878 if not v_norm2 == 0:
872 gross 2100 v[iter+1]=v[iter+1]/h[iter+1,iter]
873 artak 1519
874     # Form and store the information for the new Givens rotation
875 gross 2100 if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
876     mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
877 artak 1519
878     if mu!=0 :
879 gross 2100 c[iter]=h[iter,iter]/mu
880     s[iter]=-h[iter+1,iter]/mu
881     h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
882     h[iter+1,iter]=0.0
883 gross 1878 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
884 artak 1519 # Update the residual norm
885 gross 1878
886 artak 1519 rho=abs(g[iter+1])
887 gross 2100 if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
888 artak 1519 iter+=1
889    
890 gross 1878 # At this point either iter > iter_max or rho < tol.
891     # It's time to compute x and leave.
892    
893 gross 2100 if verbose: print "GMRES: iteration stopped after %s step."%iter
894 artak 1519 if iter > 0 :
895     y=numarray.zeros(iter,numarray.Float64)
896 gross 2100 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
897 artak 1519 if iter > 1 :
898     i=iter-2
899     while i>=0 :
900 gross 2100 y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
901 artak 1519 i=i-1
902     xhat=v[iter-1]*y[iter-1]
903     for i in range(iter-1):
904     xhat += v[i]*y[i]
905 gross 2100 else:
906     xhat=v[0] * 0
907 gross 1878
908 artak 1519 x += xhat
909     if iter<iter_restart-1:
910     stopped=True
911     else:
912     stopped=False
913    
914     return x,stopped
915    
916 artak 1550 #################################################
917 gross 2100 def MINRES(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
918 artak 1550 #################################################
919 artak 1481 #
920     # minres solves the system of linear equations Ax = b
921     # where A is a symmetric matrix (possibly indefinite or singular)
922     # and b is a given vector.
923     #
924     # "A" may be a dense or sparse matrix (preferably sparse!)
925     # or the name of a function such that
926     # y = A(x)
927     # returns the product y = Ax for any given vector x.
928     #
929     # "M" defines a positive-definite preconditioner M = C C'.
930     # "M" may be a dense or sparse matrix (preferably sparse!)
931     # or the name of a function such that
932     # solves the system My = x for any given vector x.
933     #
934     #
935 artak 1482
936 artak 1481 #------------------------------------------------------------------
937     # Set up y and v for the first Lanczos vector v1.
938     # y = beta1 P' v1, where P = C**(-1).
939     # v is really P' v1.
940     #------------------------------------------------------------------
941 artak 1482 if x==None:
942     x=0*b
943     else:
944     b += (-1)*Aprod(x)
945    
946 artak 1481 r1 = b
947     y = Msolve(b)
948 artak 1550 beta1 = bilinearform(y,b)
949 artak 1481
950     if beta1< 0: raise NegativeNorm,"negative norm."
951    
952     # If b = 0 exactly, stop with x = 0.
953     if beta1==0: return x*0.
954    
955     if beta1> 0:
956 artak 1486 beta1 = math.sqrt(beta1)
957 artak 1481
958     #------------------------------------------------------------------
959 artak 1484 # Initialize quantities.
960 artak 1481 # ------------------------------------------------------------------
961 artak 1482 iter = 0
962     Anorm = 0
963     ynorm = 0
964 artak 1481 oldb = 0
965     beta = beta1
966     dbar = 0
967     epsln = 0
968     phibar = beta1
969     rhs1 = beta1
970     rhs2 = 0
971     rnorm = phibar
972     tnorm2 = 0
973     ynorm2 = 0
974     cs = -1
975     sn = 0
976     w = b*0.
977     w2 = b*0.
978     r2 = r1
979     eps = 0.0001
980    
981     #---------------------------------------------------------------------
982     # Main iteration loop.
983     # --------------------------------------------------------------------
984 gross 2100 while not rnorm<=atol+rtol*Anorm*ynorm: # checks ||r|| < (||A|| ||x||) * TOL
985 artak 1481
986     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
987     iter = iter + 1
988    
989     #-----------------------------------------------------------------
990     # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
991     # The general iteration is similar to the case k = 1 with v0 = 0:
992     #
993     # p1 = Operator * v1 - beta1 * v0,
994     # alpha1 = v1'p1,
995     # q2 = p2 - alpha1 * v1,
996     # beta2^2 = q2'q2,
997     # v2 = (1/beta2) q2.
998     #
999     # Again, y = betak P vk, where P = C**(-1).
1000     #-----------------------------------------------------------------
1001     s = 1/beta # Normalize previous vector (in y).
1002     v = s*y # v = vk if P = I
1003    
1004     y = Aprod(v)
1005 artak 1465
1006 artak 1481 if iter >= 2:
1007     y = y - (beta/oldb)*r1
1008    
1009     alfa = bilinearform(v,y) # alphak
1010 artak 1550 y += (- alfa/beta)*r2
1011 artak 1481 r1 = r2
1012     r2 = y
1013     y = Msolve(r2)
1014     oldb = beta # oldb = betak
1015 artak 1550 beta = bilinearform(y,r2) # beta = betak+1^2
1016 artak 1481 if beta < 0: raise NegativeNorm,"negative norm."
1017    
1018     beta = math.sqrt( beta )
1019     tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1020    
1021     if iter==1: # Initialize a few things.
1022     gmax = abs( alfa ) # alpha1
1023     gmin = gmax # alpha1
1024    
1025     # Apply previous rotation Qk-1 to get
1026     # [deltak epslnk+1] = [cs sn][dbark 0 ]
1027     # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1028    
1029     oldeps = epsln
1030     delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1031     gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1032     epsln = sn * beta # epsln2 = 0 epslnk+1
1033     dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1034    
1035     # Compute the next plane rotation Qk
1036    
1037     gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1038     gamma = max(gamma,eps)
1039     cs = gbar / gamma # ck
1040     sn = beta / gamma # sk
1041     phi = cs * phibar # phik
1042     phibar = sn * phibar # phibark+1
1043    
1044     # Update x.
1045    
1046     denom = 1/gamma
1047     w1 = w2
1048     w2 = w
1049     w = (v - oldeps*w1 - delta*w2) * denom
1050 artak 1550 x += phi*w
1051 artak 1481
1052     # Go round again.
1053    
1054     gmax = max(gmax,gamma)
1055     gmin = min(gmin,gamma)
1056     z = rhs1 / gamma
1057     ynorm2 = z*z + ynorm2
1058     rhs1 = rhs2 - delta*z
1059     rhs2 = - epsln*z
1060    
1061     # Estimate various norms and test for convergence.
1062    
1063     Anorm = math.sqrt( tnorm2 )
1064     ynorm = math.sqrt( ynorm2 )
1065    
1066     rnorm = phibar
1067    
1068     return x
1069 artak 1489
1070 gross 2100 def TFQMR(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
1071 artak 1489
1072     # TFQMR solver for linear systems
1073     #
1074     #
1075     # initialization
1076     #
1077     errtol = math.sqrt(bilinearform(b,b))
1078     norm_b=errtol
1079     kmax = iter_max
1080     error = []
1081    
1082     if math.sqrt(bilinearform(x,x)) != 0.0:
1083     r = b - Aprod(x)
1084     else:
1085     r = b
1086    
1087     r=Msolve(r)
1088    
1089     u1=0
1090     u2=0
1091     y1=0
1092     y2=0
1093    
1094     w = r
1095     y1 = r
1096     iter = 0
1097     d = 0
1098    
1099     v = Msolve(Aprod(y1))
1100     u1 = v
1101    
1102     theta = 0.0;
1103     eta = 0.0;
1104     tau = math.sqrt(bilinearform(r,r))
1105     error = [ error, tau ]
1106     rho = tau * tau
1107     #
1108     # TFQMR iteration
1109     #
1110     # while ( iter < kmax-1 ):
1111    
1112 gross 2100 while not tau<=atol+rtol*norm_b:
1113 artak 1489 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1114    
1115     sigma = bilinearform(r,v)
1116    
1117     if ( sigma == 0.0 ):
1118     raise 'TFQMR breakdown, sigma=0'
1119    
1120    
1121     alpha = rho / sigma
1122    
1123     for j in range(2):
1124     #
1125     # Compute y2 and u2 only if you have to
1126     #
1127     if ( j == 1 ):
1128     y2 = y1 - alpha * v
1129     u2 = Msolve(Aprod(y2))
1130    
1131     m = 2 * (iter+1) - 2 + (j+1)
1132     if j==0:
1133     w = w - alpha * u1
1134     d = y1 + ( theta * theta * eta / alpha ) * d
1135     if j==1:
1136     w = w - alpha * u2
1137     d = y2 + ( theta * theta * eta / alpha ) * d
1138    
1139     theta = math.sqrt(bilinearform(w,w))/ tau
1140     c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1141     tau = tau * theta * c
1142     eta = c * c * alpha
1143     x = x + eta * d
1144     #
1145     # Try to terminate the iteration at each pass through the loop
1146     #
1147     # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1148     # error = [ error, tau ]
1149     # total_iters = iter
1150     # break
1151    
1152    
1153     if ( rho == 0.0 ):
1154     raise 'TFQMR breakdown, rho=0'
1155    
1156    
1157     rhon = bilinearform(r,w)
1158     beta = rhon / rho;
1159     rho = rhon;
1160     y1 = w + beta * y2;
1161     u1 = Msolve(Aprod(y1))
1162     v = u1 + beta * ( u2 + beta * v )
1163     error = [ error, tau ]
1164     total_iters = iter
1165    
1166     iter = iter + 1
1167    
1168     return x
1169    
1170    
1171 artak 1465 #############################################
1172    
1173 gross 1331 class ArithmeticTuple(object):
1174     """
1175     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1176    
1177     example of usage:
1178    
1179     from esys.escript import Data
1180     from numarray import array
1181     a=Data(...)
1182     b=array([1.,4.])
1183     x=ArithmeticTuple(a,b)
1184     y=5.*x
1185    
1186     """
1187     def __init__(self,*args):
1188     """
1189     initialize object with elements args.
1190    
1191     @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1192     """
1193     self.__items=list(args)
1194    
1195     def __len__(self):
1196     """
1197     number of items
1198    
1199     @return: number of items
1200     @rtype: C{int}
1201     """
1202     return len(self.__items)
1203    
1204     def __getitem__(self,index):
1205     """
1206     get an item
1207    
1208     @param index: item to be returned
1209     @type index: C{int}
1210     @return: item with index C{index}
1211     """
1212     return self.__items.__getitem__(index)
1213    
1214     def __mul__(self,other):
1215     """
1216     scaling from the right
1217    
1218     @param other: scaling factor
1219     @type other: C{float}
1220     @return: itemwise self*other
1221     @rtype: L{ArithmeticTuple}
1222     """
1223     out=[]
1224 gross 1896 try:
1225     l=len(other)
1226     if l!=len(self):
1227     raise ValueError,"length of of arguments don't match."
1228     for i in range(l): out.append(self[i]*other[i])
1229     except TypeError:
1230     for i in range(len(self)): out.append(self[i]*other)
1231 gross 1331 return ArithmeticTuple(*tuple(out))
1232    
1233     def __rmul__(self,other):
1234     """
1235     scaling from the left
1236    
1237     @param other: scaling factor
1238     @type other: C{float}
1239     @return: itemwise other*self
1240     @rtype: L{ArithmeticTuple}
1241     """
1242     out=[]
1243 gross 1896 try:
1244     l=len(other)
1245     if l!=len(self):
1246     raise ValueError,"length of of arguments don't match."
1247     for i in range(l): out.append(other[i]*self[i])
1248     except TypeError:
1249     for i in range(len(self)): out.append(other*self[i])
1250 gross 1331 return ArithmeticTuple(*tuple(out))
1251    
1252 artak 1465 def __div__(self,other):
1253     """
1254     dividing from the right
1255    
1256     @param other: scaling factor
1257     @type other: C{float}
1258     @return: itemwise self/other
1259     @rtype: L{ArithmeticTuple}
1260     """
1261 gross 1896 return self*(1/other)
1262 artak 1465
1263     def __rdiv__(self,other):
1264     """
1265     dividing from the left
1266    
1267     @param other: scaling factor
1268     @type other: C{float}
1269     @return: itemwise other/self
1270     @rtype: L{ArithmeticTuple}
1271     """
1272     out=[]
1273 gross 1896 try:
1274     l=len(other)
1275     if l!=len(self):
1276     raise ValueError,"length of of arguments don't match."
1277     for i in range(l): out.append(other[i]/self[i])
1278     except TypeError:
1279     for i in range(len(self)): out.append(other/self[i])
1280 artak 1465 return ArithmeticTuple(*tuple(out))
1281    
1282 gross 1331 def __iadd__(self,other):
1283     """
1284     in-place add of other to self
1285    
1286     @param other: increment
1287     @type other: C{ArithmeticTuple}
1288     """
1289     if len(self) != len(other):
1290     raise ValueError,"tuple length must match."
1291     for i in range(len(self)):
1292     self.__items[i]+=other[i]
1293     return self
1294    
1295 artak 1550 def __add__(self,other):
1296     """
1297     add other to self
1298    
1299     @param other: increment
1300     @type other: C{ArithmeticTuple}
1301     """
1302 gross 1896 out=[]
1303     try:
1304     l=len(other)
1305     if l!=len(self):
1306     raise ValueError,"length of of arguments don't match."
1307     for i in range(l): out.append(self[i]+other[i])
1308     except TypeError:
1309     for i in range(len(self)): out.append(self[i]+other)
1310     return ArithmeticTuple(*tuple(out))
1311 artak 1550
1312     def __sub__(self,other):
1313     """
1314     subtract other from self
1315    
1316     @param other: increment
1317     @type other: C{ArithmeticTuple}
1318     """
1319 gross 1896 out=[]
1320     try:
1321     l=len(other)
1322     if l!=len(self):
1323     raise ValueError,"length of of arguments don't match."
1324     for i in range(l): out.append(self[i]-other[i])
1325     except TypeError:
1326     for i in range(len(self)): out.append(self[i]-other)
1327     return ArithmeticTuple(*tuple(out))
1328 artak 1557
1329     def __isub__(self,other):
1330     """
1331     subtract other from self
1332 artak 1550
1333 artak 1557 @param other: increment
1334     @type other: C{ArithmeticTuple}
1335     """
1336     if len(self) != len(other):
1337     raise ValueError,"tuple length must match."
1338     for i in range(len(self)):
1339     self.__items[i]-=other[i]
1340     return self
1341    
1342 artak 1550 def __neg__(self):
1343     """
1344     negate
1345    
1346     """
1347 gross 1896 out=[]
1348 artak 1550 for i in range(len(self)):
1349 gross 1896 out.append(-self[i])
1350     return ArithmeticTuple(*tuple(out))
1351 artak 1550
1352    
1353 gross 1414 class HomogeneousSaddlePointProblem(object):
1354     """
1355 gross 1878 This provides a framwork for solving linear homogeneous saddle point problem of the form
1356 gross 1414
1357     Av+B^*p=f
1358     Bv =0
1359    
1360     for the unknowns v and p and given operators A and B and given right hand side f.
1361 gross 2100 B^* is the adjoint operator of B.
1362 gross 1414 """
1363     def __init__(self,**kwargs):
1364     self.setTolerance()
1365 gross 2100 self.setAbsoluteTolerance()
1366     self.setSubToleranceReductionFactor()
1367 gross 1414
1368 gross 2100 #=============================================================
1369 gross 1414 def initialize(self):
1370     """
1371     initialize the problem (overwrite)
1372     """
1373     pass
1374 gross 2100
1375 gross 1414 def B(self,v):
1376 gross 2100 """
1377     returns Bv (overwrite)
1378     @rtype: equal to the type of p
1379    
1380     @note: boundary conditions on p should be zero!
1381     """
1382     raise NotImplementedError,"no operator B implemented."
1383    
1384     def inner_pBv(self,p,Bv):
1385 gross 1414 """
1386 gross 2100 returns inner product of element p and Bv (overwrite)
1387    
1388     @type p: equal to the type of p
1389     @type Bv: equal to the type of result of operator B
1390     @rtype: C{float}
1391    
1392 gross 1414 @rtype: equal to the type of p
1393     """
1394 gross 2100 raise NotImplementedError,"no inner product for p implemented."
1395 gross 1414
1396 gross 2100 def inner_p(self,p0,p1):
1397 gross 1414 """
1398 gross 2100 returns inner product of element p0 and p1 (overwrite)
1399 gross 1414
1400     @type p0: equal to the type of p
1401     @type p1: equal to the type of p
1402     @rtype: equal to the type of p
1403     """
1404 gross 2100 raise NotImplementedError,"no inner product for p implemented."
1405    
1406     def inner_v(self,v0,v1):
1407     """
1408     returns inner product of two element v0 and v1 (overwrite)
1409    
1410     @type v0: equal to the type of v
1411     @type v1: equal to the type of v
1412     @rtype: C{float}
1413    
1414     @rtype: equal to the type of v
1415     """
1416     raise NotImplementedError,"no inner product for v implemented."
1417 gross 1414 pass
1418    
1419     def solve_A(self,u,p):
1420     """
1421 gross 2100 solves Av=f-Au-B^*p with accuracy self.getSubProblemTolerance() (overwrite)
1422 gross 1414
1423     @rtype: equal to the type of v
1424     @note: boundary conditions on v should be zero!
1425     """
1426 gross 2100 raise NotImplementedError,"no operator A implemented."
1427 gross 1414
1428     def solve_prec(self,p):
1429     """
1430 gross 2100 provides a preconditioner for BA^{-1}B^* with accuracy self.getSubProblemTolerance() (overwrite)
1431 gross 1414
1432     @rtype: equal to the type of p
1433 gross 2100 @note: boundary conditions on p should be zero!
1434 gross 1414 """
1435 gross 2100 raise NotImplementedError,"no preconditioner for Schur complement implemented."
1436     #=============================================================
1437     def __Aprod_PCG(self,p):
1438     # return (v,Bv) with v=A^-1B*p
1439     #solve Av =B^*p as Av =f-Az-B^*(-p)
1440     v=self.solve_A(self.__z,-p)
1441     return ArithmeticTuple(v, self.B(v))
1442 gross 1414
1443 gross 2100 def __inner_PCG(self,p,r):
1444     return self.inner_pBv(p,r[1])
1445 gross 1414
1446 gross 2100 def __Msolve_PCG(self,r):
1447     return self.solve_prec(r[1])
1448     #=============================================================
1449     def __Aprod_GMRES(self,x):
1450     # return w,q from v, p
1451     # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)
1452     # and Sq=B(v-w)
1453     v=x[0]
1454     p=x[1]
1455     w=self.solve_A(self.__z-v,-p)
1456     Bw=self.B(w-v)
1457     q=self.solve_prec(Bw)
1458     return ArithmeticTuple(w,q)
1459 gross 1414
1460 gross 2100 def __inner_GMRES(self,a1,a2):
1461     return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])
1462 artak 1465
1463 gross 2100 #=============================================================
1464     def norm(self,v,p):
1465     f=self.inner_p(p,p)+self.inner_v(v,v)
1466     if f<0:
1467     raise ValueError,"negative norm."
1468     return math.sqrt(f)
1469 artak 1550
1470 gross 2100 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):
1471     """
1472     solves the saddle point problem using initial guesses v and p.
1473 gross 1414
1474 gross 2100 @param v: initial guess for velocity
1475     @param p: initial guess for pressure
1476     @type v: L{Data}
1477     @type p: L{Data}
1478     @param useUzawa: indicate the usage of the Uzawa scheme. Otherwise the problem is solved in coupled form.
1479     @param max_iter: maximum number of iteration steps per correction attempt.
1480     @param verbose: show information on the progress of the saddlepoint problem solver.
1481     @param show_details: show details of the sub problems solves
1482     @param iter_restart: restart the iteration after C{iter_restart} steps (only used if useUzaw=False)
1483     @param max_correction_steps: maximum number of iteration steps in the attempt get |Bv| to zero.
1484     @return: new approximations for velocity and pressure
1485     @type useUzawa: C{bool}
1486     @type max_iter: C{int}
1487     @type verbose: C{bool}
1488     @type show_details: C{bool}
1489     @type iter_restart: C{int}
1490     @type max_correction_steps: C{int}
1491     @rtype: C{tuple} of L{Data} objects
1492     """
1493     self.verbose=verbose
1494     self.show_details=show_details and self.verbose
1495     #=================================================================================
1496     # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
1497     self.setSubProblemTolerance(self.getTolerance(), apply_reduction=True)
1498     self.__z=v+self.solve_A(v,p*0)
1499     Bz=self.B(self.__z)
1500     p0=self.solve_prec(Bz)
1501     # tolerances:
1502     rtol=self.getTolerance()
1503     atol=self.getAbsoluteTolerance()
1504     f0=self.norm(self.__z,p0)
1505     if not f0>0:
1506     return self.__z, p*0
1507     ATOL=rtol*f0+atol
1508     if not ATOL>0:
1509     raise ValueError,"overall absolute tolerance needs to be positive."
1510     if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
1511     # initialization
1512     self.iter=0
1513     correction_step=0
1514     converged=False
1515     # initial guess:
1516     q=p*1
1517     if not useUzawa:
1518     u=v
1519     while not converged :
1520     if useUzawa:
1521     self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1522     # assume p is known: then v=z-A^-1B^*p
1523     #
1524     # which leads to BA^-1B^*p = Bz
1525     #
1526     # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bu
1527     #
1528     # note that q=BA^-1B^*p = B (A^-1B^*p)=Bw with Aw=B^*p = f - Az - B^*(-p)
1529     #
1530     if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())
1531     q,r=PCG(ArithmeticTuple(self.__z,Bz),self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=q, verbose=self.verbose)
1532     u=r[0]
1533     Bu=r[1]
1534     else:
1535     self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1536     #
1537     # with v=dv+z
1538     #
1539     # A dv + B p = A (v-z) + Bp = 0
1540     # B dv = - Bz
1541     #
1542     # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1}) -S^{-1}]] to the problem
1543     # the right hand side is then [0, S^{-1} Bz ]
1544     #
1545     #
1546     if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
1547     x=GMRES(ArithmeticTuple(0*self.__z,p0),self.__Aprod_GMRES,self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, x=ArithmeticTuple(u-self.__z,q),iter_restart=iter_restart, verbose=self.verbose)
1548     self.setSubProblemTolerance(ATOL/f0, apply_reduction=False)
1549     u=self.__z+x[0]
1550     q=x[1]
1551     Bu=self.B(u)
1552     # now we check if |Bu| ~ 0 or more precise |Bu| <= rtol * |v|
1553     nu=self.inner_v(u,u)
1554     nz=self.inner_v(self.__z,self.__z)
1555     p2=self.solve_prec(Bu)
1556     nBu=self.inner_p(p2,p2)
1557     if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
1558     nu=math.sqrt(nu)
1559     nBu=math.sqrt(nBu)
1560     if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)
1561     QTOL=atol+nu*rtol
1562     if nBu <= QTOL:
1563     converged=True
1564     else:
1565     ATOL=QTOL/nBu*ATOL*0.3
1566     if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)
1567     converged=False
1568     correction_step+=1
1569     if correction_step>max_correction_steps:
1570     raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1571     if self.verbose: print "saddle point solver: tolerance reached."
1572     return u,q
1573 artak 1465
1574 gross 2100 #==========================================================================================================================
1575     def setTolerance(self,tolerance=1.e-4):
1576     """
1577     sets the relative tolerance for (v,p)
1578    
1579     @param tolerance: tolerance to be used
1580     @type tolerance: non-negative C{float}
1581     """
1582     if tolerance<0:
1583     raise ValueError,"tolerance must be positive."
1584     self.__rtol=tolerance
1585 gross 1414 def getTolerance(self):
1586 gross 2100 """
1587     returns the relative tolerance
1588 gross 1414
1589 gross 2100 @return: relative tolerance
1590     @rtype: C{float}
1591     """
1592     return self.__rtol
1593     def setAbsoluteTolerance(self,tolerance=0.):
1594     """
1595     sets the absolute tolerance
1596 gross 1414
1597 gross 2100 @param tolerance: tolerance to be used
1598     @type tolerance: non-negative C{float}
1599     """
1600     if tolerance<0:
1601     raise ValueError,"tolerance must be non-negative."
1602     self.__atol=tolerance
1603     def getAbsoluteTolerance(self):
1604     """
1605     returns the absolute tolerance
1606 gross 1414
1607 gross 2100 @return: absolute tolerance
1608     @rtype: C{float}
1609     """
1610     return self.__atol
1611     def setSubToleranceReductionFactor(self,reduction=None):
1612     """
1613     sets the reduction factor for the tolerance used to solve the subproblem in each iteration step.
1614     if set to C{None} the square of the problem tolerance is used.
1615    
1616     @param reduction: reduction factor.
1617     @type reduction: positive C{float} less or equal one or C{None}
1618     """
1619     if not reduction == None:
1620     if reduction<=0 or reduction>1:
1621     raise ValueError,"reduction factor (=%e) must be positive and less or equal one."%reduction
1622     self.__reduction=reduction
1623     def getSubToleranceReductionFactor(self):
1624     """
1625     returns the subproblem reduction factor
1626 gross 1469
1627 gross 2100 @return: subproblem reduction factor
1628     @rtype: C{float} or C{None}
1629     """
1630     return self.__reduction
1631 artak 1481
1632 gross 2100 def setSubProblemTolerance(self,rtol,apply_reduction=False):
1633     """
1634     sets the relative tolerance to solve the subproblem(s).
1635 artak 1489
1636 gross 2100 @param rtol: relative tolerence
1637     @type rtol: positive C{float}
1638     @param apply_reduction: if set the the reduction factor is applied.
1639     @type apply_reduction: C{bool}
1640     """
1641     if rtol<=0:
1642     raise ValueError,"tolerance must be positive."
1643     if apply_reduction:
1644     if self.__reduction == None:
1645     self.__sub_tol=max(rtol**2,util.EPSILON*10)
1646     else:
1647     self.__sub_tol=max(self.__reduction*rtol,util.EPSILON*10)
1648     else:
1649     self.__sub_tol=rtol
1650     def getSubProblemTolerance(self):
1651     """
1652     returns the subproblem reduction factor
1653 artak 1557
1654 gross 2100 @return: subproblem reduction factor
1655     @rtype: C{float}
1656     """
1657     return self.__sub_tol
1658 artak 1550
1659 gross 1878 def MaskFromBoundaryTag(domain,*tags):
1660     """
1661 gross 2100 create a mask on the Solution(domain) function space which one for samples
1662 gross 1878 that touch the boundary tagged by tags.
1663    
1664     usage: m=MaskFromBoundaryTag(domain,"left", "right")
1665    
1666     @param domain: a given domain
1667     @type domain: L{escript.Domain}
1668     @param tags: boundray tags
1669     @type tags: C{str}
1670     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1671     @rtype: L{escript.Data} of rank 0
1672     """
1673     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1674 gross 1956 d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1675 gross 1878 for t in tags: d.setTaggedValue(t,1.)
1676     pde.setValue(y=d)
1677     return util.whereNonZero(pde.getRightHandSide())
1678     #============================================================================================================================================
1679 gross 867 class SaddlePointProblem(object):
1680     """
1681     This implements a solver for a saddlepoint problem
1682    
1683 gross 877 M{f(u,p)=0}
1684     M{g(u)=0}
1685 gross 867
1686     for u and p. The problem is solved with an inexact Uszawa scheme for p:
1687    
1688 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1689 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1690 gross 867
1691     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1692     A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
1693     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1694     in fact the role of a preconditioner.
1695     """
1696     def __init__(self,verbose=False,*args):
1697     """
1698     initializes the problem
1699    
1700 ksteube 990 @param verbose: switches on the printing out some information
1701 gross 867 @type verbose: C{bool}
1702     @note: this method may be overwritten by a particular saddle point problem
1703     """
1704 gross 1659 print "SaddlePointProblem should not be used anymore!"
1705 gross 1107 if not isinstance(verbose,bool):
1706     raise TypeError("verbose needs to be of type bool.")
1707 gross 1106 self.__verbose=verbose
1708 gross 877 self.relaxation=1.
1709 gross 1878 DeprecationWarning("SaddlePointProblem should not be used anymore.")
1710 gross 867
1711     def trace(self,text):
1712     """
1713     prints text if verbose has been set
1714    
1715 ksteube 990 @param text: a text message
1716 gross 867 @type text: C{str}
1717     """
1718 ksteube 1567 if self.__verbose: print "%s: %s"%(str(self),text)
1719 gross 867
1720 gross 873 def solve_f(self,u,p,tol=1.e-8):
1721 gross 867 """
1722     solves
1723    
1724     A_f du = f(u,p)
1725    
1726     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1727    
1728     @param u: current approximation of u
1729     @type u: L{escript.Data}
1730     @param p: current approximation of p
1731     @type p: L{escript.Data}
1732 gross 873 @param tol: tolerance expected for du
1733 gross 867 @type tol: C{float}
1734     @return: increment du
1735     @rtype: L{escript.Data}
1736     @note: this method has to be overwritten by a particular saddle point problem
1737     """
1738     pass
1739    
1740 gross 873 def solve_g(self,u,tol=1.e-8):
1741 gross 867 """
1742     solves
1743    
1744     Q_g dp = g(u)
1745    
1746     with Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p.
1747    
1748     @param u: current approximation of u
1749     @type u: L{escript.Data}
1750 gross 873 @param tol: tolerance expected for dp
1751     @type tol: C{float}
1752 gross 867 @return: increment dp
1753     @rtype: L{escript.Data}
1754     @note: this method has to be overwritten by a particular saddle point problem
1755     """
1756     pass
1757    
1758     def inner(self,p0,p1):
1759     """
1760     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1761     @return: inner product of p0 and p1
1762     @rtype: C{float}
1763     """
1764     pass
1765    
1766 gross 877 subiter_max=3
1767     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1768     """
1769     runs the solver
1770 gross 873
1771 gross 877 @param u0: initial guess for C{u}
1772     @type u0: L{esys.escript.Data}
1773     @param p0: initial guess for C{p}
1774     @type p0: L{esys.escript.Data}
1775     @param tolerance: tolerance for relative error in C{u} and C{p}
1776     @type tolerance: positive C{float}
1777     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1778     @type tolerance_u: positive C{float}
1779     @param iter_max: maximum number of iteration steps.
1780     @type iter_max: C{int}
1781     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1782     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1783     @type accepted_reduction: positive C{float} or C{None}
1784     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1785     @type relaxation: C{float} or C{None}
1786     """
1787     tol=1.e-2
1788     if tolerance_u==None: tolerance_u=tolerance
1789     if not relaxation==None: self.relaxation=relaxation
1790     if accepted_reduction ==None:
1791     angle_limit=0.
1792     elif accepted_reduction>=1.:
1793     angle_limit=0.
1794     else:
1795     angle_limit=util.sqrt(1-accepted_reduction**2)
1796     self.iter=0
1797     u=u0
1798     p=p0
1799     #
1800     # initialize things:
1801     #
1802     converged=False
1803     #
1804     # start loop:
1805     #
1806     # initial search direction is g
1807     #
1808     while not converged :
1809     if self.iter>iter_max:
1810     raise ArithmeticError("no convergence after %s steps."%self.iter)
1811     f_new=self.solve_f(u,p,tol)
1812     norm_f_new = util.Lsup(f_new)
1813     u_new=u-f_new
1814     g_new=self.solve_g(u_new,tol)
1815     self.iter+=1
1816     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1817     if norm_f_new==0. and norm_g_new==0.: return u, p
1818     if self.iter>1 and not accepted_reduction==None:
1819     #
1820     # did we manage to reduce the norm of G? I
1821     # if not we start a backtracking procedure
1822     #
1823     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1824     if norm_g_new > accepted_reduction * norm_g:
1825     sub_iter=0
1826     s=self.relaxation
1827     d=g
1828     g_last=g
1829     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1830     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1831     dg= g_new-g_last
1832     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1833     rad=self.inner(g_new,dg)/self.relaxation
1834     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1835     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1836     if abs(rad) < norm_dg*norm_g_new * angle_limit:
1837     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1838     break
1839     r=self.relaxation
1840     self.relaxation= - rad/norm_dg**2
1841     s+=self.relaxation
1842     #####
1843     # a=g_new+self.relaxation*dg/r
1844     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1845     #####
1846     g_last=g_new
1847     p+=self.relaxation*d
1848     f_new=self.solve_f(u,p,tol)
1849     u_new=u-f_new
1850     g_new=self.solve_g(u_new,tol)
1851     self.iter+=1
1852     norm_f_new = util.Lsup(f_new)
1853     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1854     # print " ",sub_iter," new g norm",norm_g_new
1855     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1856     #
1857     # can we expect reduction of g?
1858     #
1859     # u_last=u_new
1860     sub_iter+=1
1861     self.relaxation=s
1862     #
1863     # check for convergence:
1864     #
1865     norm_u_new = util.Lsup(u_new)
1866     p_new=p+self.relaxation*g_new
1867     norm_p_new = util.sqrt(self.inner(p_new,p_new))
1868 artak 1550 self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
1869 gross 873
1870 gross 877 if self.iter>1:
1871     dg2=g_new-g
1872     df2=f_new-f
1873     norm_dg2=util.sqrt(self.inner(dg2,dg2))
1874     norm_df2=util.Lsup(df2)
1875     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1876     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1877     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1878     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1879     converged=True
1880     f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
1881     self.trace("convergence after %s steps."%self.iter)
1882     return u,p

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26