/[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 2061 - (hide annotations)
Wed Nov 19 03:40:21 2008 UTC (11 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 59620 byte(s)
Fixing some warnings from epydoc.
Disabled c++ signatures in python docstrings.
Removed references to Bruce in epydoc and the users guide.
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     class IterationBreakDown(SolverSchemeException):
430     """
431     iteration scheme econouters an incurable breakdown.
432     """
433     pass
434     class NegativeNorm(SolverSchemeException):
435     """
436     a norm calculation returns a negative norm.
437     """
438     pass
439    
440 gross 1330 class IterationHistory(object):
441 ksteube 1312 """
442 gross 1330 The IterationHistory class is used to define a stopping criterium. It keeps track of the
443     residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
444     a given tolerance.
445     """
446     def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
447     """
448     Initialization
449    
450     @param tolerance: tolerance
451     @type tolerance: positive C{float}
452     @param verbose: switches on the printing out some information
453     @type verbose: C{bool}
454     """
455     if not tolerance>0.:
456     raise ValueError,"tolerance needs to be positive."
457     self.tolerance=tolerance
458     self.verbose=verbose
459     self.history=[]
460     def stoppingcriterium(self,norm_r,r,x):
461     """
462     returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]} is the residual norm at the first call.
463    
464    
465     @param norm_r: current residual norm
466     @type norm_r: non-negative C{float}
467     @param r: current residual (not used)
468     @param x: current solution approximation (not used)
469     @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
470     @rtype: C{bool}
471    
472     """
473     self.history.append(norm_r)
474 ksteube 1567 if self.verbose: print "iter: %s: inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])
475 gross 1330 return self.history[-1]<=self.tolerance * self.history[0]
476    
477 artak 1519 def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
478 gross 1469 """
479     returns True if the C{norm_r} is C{tolerance}*C{norm_b}
480    
481    
482     @param norm_r: current residual norm
483     @type norm_r: non-negative C{float}
484     @param norm_b: norm of right hand side
485     @type norm_b: non-negative C{float}
486     @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
487     @rtype: C{bool}
488    
489     """
490 artak 1519 if TOL==None:
491     TOL=self.tolerance
492 gross 1469 self.history.append(norm_r)
493 ksteube 1567 if self.verbose: print "iter: %s: norm(r) = %e"#(len(self.history)-1, self.history[-1])
494 artak 1519 return self.history[-1]<=TOL * norm_b
495 gross 1469
496 gross 1330 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
497     """
498 ksteube 1312 Solver for
499    
500     M{Ax=b}
501    
502     with a symmetric and positive definite operator A (more details required!).
503     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
504    
505 gross 1330 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
506 ksteube 1312
507     For details on the preconditioned conjugate gradient method see the book:
508    
509     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
510     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
511     C. Romine, and H. van der Vorst.
512    
513     @param b: the right hand side of the liner system. C{b} is altered.
514 gross 1330 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
515 ksteube 1312 @param Aprod: returns the value Ax
516 jfenwick 2061 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.
517     The returned object needs to be of the same type like argument C{b}.
518     @param Msolve: solves Mx=r
519     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}.
520     The returned object needs to be of the same type like argument C{x}.
521     @param bilinearform: inner product C{<x,r>}
522     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.
523     The returned value is a C{float}.
524     @param stoppingcriterium: function which returns True if a stopping criterium is meet.
525     C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current
526     norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and
527     the current solution approximation. C{stoppingcriterium} is called in each iteration step.
528 gross 1330 @type stoppingcriterium: function that returns C{True} or C{False}
529     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
530     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
531 ksteube 1312 @param iter_max: maximum number of iteration steps.
532     @type iter_max: C{int}
533 gross 1330 @return: the solution approximation and the corresponding residual
534     @rtype: C{tuple}
535     @warning: C{b} and C{x} are altered.
536 ksteube 1312 """
537     iter=0
538 gross 1330 if x==None:
539     x=0*b
540     else:
541     b += (-1)*Aprod(x)
542 ksteube 1312 r=b
543     rhat=Msolve(r)
544 gross 1330 d = rhat
545 ksteube 1312 rhat_dot_r = bilinearform(rhat, r)
546 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
547 ksteube 1312
548 gross 1330 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
549     iter+=1
550 ksteube 1312 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
551    
552     q=Aprod(d)
553     alpha = rhat_dot_r / bilinearform(d, q)
554     x += alpha * d
555     r += (-alpha) * q
556    
557     rhat=Msolve(r)
558     rhat_dot_r_new = bilinearform(rhat, r)
559     beta = rhat_dot_r_new / rhat_dot_r
560     rhat+=beta * d
561     d=rhat
562    
563     rhat_dot_r = rhat_dot_r_new
564 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
565 ksteube 1312
566 gross 1330 return x,r
567 ksteube 1312
568 gross 1878 class Defect(object):
569     """
570     defines a non-linear defect F(x) of a variable x
571     """
572     def __init__(self):
573     """
574     initialize defect
575     """
576     self.setDerivativeIncrementLength()
577 artak 1465
578 gross 1878 def bilinearform(self, x0, x1):
579     """
580     returns the inner product of x0 and x1
581     @param x0: a value for x
582     @param x1: a value for x
583     @return: the inner product of x0 and x1
584     @rtype: C{float}
585     """
586     return 0
587    
588     def norm(self,x):
589     """
590     the norm of argument C{x}
591    
592     @param x: a value for x
593     @return: norm of argument x
594     @rtype: C{float}
595     @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
596     """
597     s=self.bilinearform(x,x)
598     if s<0: raise NegativeNorm,"negative norm."
599     return math.sqrt(s)
600 artak 1465
601 gross 1878
602     def eval(self,x):
603     """
604     returns the value F of a given x
605    
606     @param x: value for which the defect C{F} is evalulated.
607     @return: value of the defect at C{x}
608     """
609     return 0
610    
611     def __call__(self,x):
612     return self.eval(x)
613    
614     def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
615     """
616     sets the relative length of the increment used to approximate the derivative of the defect
617     the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
618    
619     @param inc: relative increment length
620     @type inc: positive C{float}
621     """
622     if inc<=0: raise ValueError,"positive increment required."
623     self.__inc=inc
624    
625     def getDerivativeIncrementLength(self):
626     """
627     returns the relative increment length used to approximate the derivative of the defect
628     @return: value of the defect at C{x}
629     @rtype: positive C{float}
630     """
631     return self.__inc
632    
633     def derivative(self, F0, x0, v, v_is_normalised=True):
634     """
635     returns the directional derivative at x0 in the direction of v
636    
637     @param F0: value of this defect at x0
638     @param x0: value at which derivative is calculated.
639     @param v: direction
640     @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
641     @return: derivative of this defect at x0 in the direction of C{v}
642     @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
643     maybe oepsnew verwritten to use exact evalution.
644     """
645     normx=self.norm(x0)
646     if normx>0:
647     epsnew = self.getDerivativeIncrementLength() * normx
648     else:
649     epsnew = self.getDerivativeIncrementLength()
650     if not v_is_normalised:
651     normv=self.norm(v)
652     if normv<=0:
653     return F0*0
654     else:
655     epsnew /= normv
656     F1=self.eval(x0 + epsnew * v)
657     return (F1-F0)/epsnew
658    
659     ######################################
660     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):
661     """
662     solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
663    
664     M{norm(F(x) <= atol + rtol * norm(F(x0)}
665    
666     where M{x0} is the initial guess.
667    
668     @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
669     @type defect: L{Defect}
670     @param x: initial guess for the solution, C{x} is altered.
671     @type x: any object type allowing basic operations such as L{numarray.NumArray}, L{Data}
672     @param iter_max: maximum number of iteration steps
673     @type iter_max: positive C{int}
674     @param sub_iter_max:
675     @type sub_iter_max:
676     @param atol: absolute tolerance for the solution
677     @type atol: positive C{float}
678     @param rtol: relative tolerance for the solution
679     @type rtol: positive C{float}
680     @param gamma: tolerance safety factor for inner iteration
681     @type gamma: positive C{float}, less than 1
682     @param sub_tol_max: upper bound for inner tolerance.
683     @type sub_tol_max: positive C{float}, less than 1
684     @return: an approximation of the solution with the desired accuracy
685     @rtype: same type as the initial guess.
686     """
687     lmaxit=iter_max
688     if atol<0: raise ValueError,"atol needs to be non-negative."
689     if rtol<0: raise ValueError,"rtol needs to be non-negative."
690     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
691     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
692     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
693    
694     F=defect(x)
695     fnrm=defect.norm(F)
696     stop_tol=atol + rtol*fnrm
697     sub_tol=sub_tol_max
698     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
699     if verbose: print " tolerance = %e."%sub_tol
700     iter=1
701     #
702     # main iteration loop
703     #
704     while not fnrm<=stop_tol:
705     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
706     #
707     # adjust sub_tol_
708     #
709     if iter > 1:
710     rat=fnrm/fnrmo
711     sub_tol_old=sub_tol
712     sub_tol=gamma*rat**2
713     if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
714     sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
715     #
716     # calculate newton increment xc
717     # if iter_max in __FDGMRES is reached MaxIterReached is thrown
718     # if iter_restart -1 is returned as sub_iter
719     # if atol is reached sub_iter returns the numer of steps performed to get there
720     #
721     #
722     if verbose: print " subiteration (GMRES) is called with relative tolerance %e."%sub_tol
723     try:
724     xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
725     except MaxIterReached:
726     raise MaxIterReached,"maximum number of %s steps reached."%iter_max
727     if sub_iter<0:
728     iter+=sub_iter_max
729     else:
730     iter+=sub_iter
731     # ====
732     x+=xc
733     F=defect(x)
734     iter+=1
735     fnrmo, fnrm=fnrm, defect.norm(F)
736     if verbose: print " step %s: residual %e."%(iter,fnrm)
737     if verbose: print "NewtonGMRES: completed after %s steps."%iter
738     return x
739    
740     def __givapp(c,s,vin):
741     """
742     apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
743     @warning: C{vin} is altered.
744     """
745     vrot=vin
746 gross 1467 if isinstance(c,float):
747     vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
748     else:
749     for i in range(len(c)):
750     w1=c[i]*vrot[i]-s[i]*vrot[i+1]
751     w2=s[i]*vrot[i]+c[i]*vrot[i+1]
752     vrot[i:i+2]=w1,w2
753 artak 1465 return vrot
754    
755 gross 1878 def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
756 artak 1475 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
757     c=numarray.zeros(iter_restart,numarray.Float64)
758     s=numarray.zeros(iter_restart,numarray.Float64)
759     g=numarray.zeros(iter_restart,numarray.Float64)
760 artak 1465 v=[]
761    
762 gross 1878 rho=defect.norm(F0)
763     if rho<=0.: return x0*0
764 artak 1488
765 gross 1878 v.append(-F0/rho)
766 artak 1465 g[0]=rho
767 gross 1878 iter=0
768     while rho > atol and iter<iter_restart-1:
769 artak 1557
770 artak 1465 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
771    
772 gross 1878 p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
773 artak 1465 v.append(p)
774    
775 gross 1878 v_norm1=defect.norm(v[iter+1])
776 artak 1465
777 gross 1878 # Modified Gram-Schmidt
778 artak 1557 for j in range(iter+1):
779 gross 1878 h[j][iter]=defect.bilinearform(v[j],v[iter+1])
780     v[iter+1]-=h[j][iter]*v[j]
781 artak 1465
782 gross 1878 h[iter+1][iter]=defect.norm(v[iter+1])
783 artak 1465 v_norm2=h[iter+1][iter]
784    
785 gross 1878 # Reorthogonalize if needed
786 artak 1465 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
787 gross 1878 for j in range(iter+1):
788     hr=defect.bilinearform(v[j],v[iter+1])
789     h[j][iter]=h[j][iter]+hr
790     v[iter+1] -= hr*v[j]
791 artak 1465
792 gross 1878 v_norm2=defect.norm(v[iter+1])
793     h[iter+1][iter]=v_norm2
794     # watch out for happy breakdown
795 artak 1550 if not v_norm2 == 0:
796 gross 1878 v[iter+1]=v[iter+1]/h[iter+1][iter]
797 artak 1465
798 gross 1878 # Form and store the information for the new Givens rotation
799 artak 1465 if iter > 0 :
800 artak 1550 hhat=numarray.zeros(iter+1,numarray.Float64)
801     for i in range(iter+1) : hhat[i]=h[i][iter]
802 gross 1878 hhat=__givapp(c[0:iter],s[0:iter],hhat);
803 artak 1465 for i in range(iter+1) : h[i][iter]=hhat[i]
804    
805     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
806 artak 1557
807 artak 1465 if mu!=0 :
808     c[iter]=h[iter][iter]/mu
809     s[iter]=-h[iter+1][iter]/mu
810     h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
811     h[iter+1][iter]=0.0
812 gross 1878 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
813 artak 1465
814 gross 1878 # Update the residual norm
815 artak 1465 rho=abs(g[iter+1])
816     iter+=1
817    
818 gross 1878 # At this point either iter > iter_max or rho < tol.
819     # It's time to compute x and leave.
820 artak 1465 if iter > 0 :
821 gross 1467 y=numarray.zeros(iter,numarray.Float64)
822 artak 1465 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
823     if iter > 1 :
824     i=iter-2
825     while i>=0 :
826 gross 1467 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
827 artak 1465 i=i-1
828     xhat=v[iter-1]*y[iter-1]
829     for i in range(iter-1):
830     xhat += v[i]*y[i]
831 gross 1878 else :
832     xhat=v[0] * 0
833 artak 1557
834 artak 1488 if iter<iter_restart-1:
835 gross 1878 stopped=iter
836 artak 1475 else:
837 gross 1878 stopped=-1
838 artak 1465
839 gross 1878 return xhat,stopped
840 artak 1481
841 gross 1878 ##############################################
842     def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
843     ################################################
844     m=iter_restart
845     iter=0
846     xc=x
847     while True:
848     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
849     xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
850     if stopped: break
851     iter+=iter_restart
852     return xc
853 artak 1550
854 gross 1878 def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
855     iter=0
856     r=Msolve(b)
857     r_dot_r = bilinearform(r, r)
858     if r_dot_r<0: raise NegativeNorm,"negative norm."
859     norm_b=math.sqrt(r_dot_r)
860 artak 1519
861     if x==None:
862 gross 1878 x=0*b
863     else:
864     r=Msolve(b-Aprod(x))
865     r_dot_r = bilinearform(r, r)
866     if r_dot_r<0: raise NegativeNorm,"negative norm."
867 artak 1519
868     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
869     c=numarray.zeros(iter_restart,numarray.Float64)
870     s=numarray.zeros(iter_restart,numarray.Float64)
871     g=numarray.zeros(iter_restart,numarray.Float64)
872     v=[]
873    
874     rho=math.sqrt(r_dot_r)
875    
876     v.append(r/rho)
877     g[0]=rho
878    
879 gross 1878 while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
880 artak 1519
881     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
882    
883 gross 1878 p=Msolve(Aprod(v[iter]))
884 artak 1519
885     v.append(p)
886    
887     v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
888    
889     # Modified Gram-Schmidt
890 gross 1878 for j in range(iter+1):
891 artak 1519 h[j][iter]=bilinearform(v[j],v[iter+1])
892 gross 1878 v[iter+1]-=h[j][iter]*v[j]
893 artak 1519
894     h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
895     v_norm2=h[iter+1][iter]
896    
897     # Reorthogonalize if needed
898     if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
899 gross 1878 for j in range(iter+1):
900 artak 1519 hr=bilinearform(v[j],v[iter+1])
901 gross 1878 h[j][iter]=h[j][iter]+hr
902     v[iter+1] -= hr*v[j]
903 artak 1519
904     v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
905     h[iter+1][iter]=v_norm2
906    
907     # watch out for happy breakdown
908 gross 1878 if not v_norm2 == 0:
909 artak 1519 v[iter+1]=v[iter+1]/h[iter+1][iter]
910    
911     # Form and store the information for the new Givens rotation
912     if iter > 0 :
913 gross 1878 hhat=numarray.zeros(iter+1,numarray.Float64)
914     for i in range(iter+1) : hhat[i]=h[i][iter]
915     hhat=__givapp(c[0:iter],s[0:iter],hhat);
916 artak 1519 for i in range(iter+1) : h[i][iter]=hhat[i]
917    
918     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
919 gross 1878
920 artak 1519 if mu!=0 :
921     c[iter]=h[iter][iter]/mu
922     s[iter]=-h[iter+1][iter]/mu
923     h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
924     h[iter+1][iter]=0.0
925 gross 1878 g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
926 artak 1519
927     # Update the residual norm
928 gross 1878
929 artak 1519 rho=abs(g[iter+1])
930     iter+=1
931    
932 gross 1878 # At this point either iter > iter_max or rho < tol.
933     # It's time to compute x and leave.
934    
935 artak 1519 if iter > 0 :
936     y=numarray.zeros(iter,numarray.Float64)
937     y[iter-1] = g[iter-1] / h[iter-1][iter-1]
938     if iter > 1 :
939     i=iter-2
940     while i>=0 :
941     y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
942     i=i-1
943     xhat=v[iter-1]*y[iter-1]
944     for i in range(iter-1):
945     xhat += v[i]*y[i]
946     else : xhat=v[0]
947 gross 1878
948 artak 1519 x += xhat
949     if iter<iter_restart-1:
950     stopped=True
951     else:
952     stopped=False
953    
954     return x,stopped
955    
956 artak 1550 #################################################
957 artak 1481 def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
958 artak 1550 #################################################
959 artak 1481 #
960     # minres solves the system of linear equations Ax = b
961     # where A is a symmetric matrix (possibly indefinite or singular)
962     # and b is a given vector.
963     #
964     # "A" may be a dense or sparse matrix (preferably sparse!)
965     # or the name of a function such that
966     # y = A(x)
967     # returns the product y = Ax for any given vector x.
968     #
969     # "M" defines a positive-definite preconditioner M = C C'.
970     # "M" may be a dense or sparse matrix (preferably sparse!)
971     # or the name of a function such that
972     # solves the system My = x for any given vector x.
973     #
974     #
975 artak 1482
976 artak 1481 #------------------------------------------------------------------
977     # Set up y and v for the first Lanczos vector v1.
978     # y = beta1 P' v1, where P = C**(-1).
979     # v is really P' v1.
980     #------------------------------------------------------------------
981 artak 1482 if x==None:
982     x=0*b
983     else:
984     b += (-1)*Aprod(x)
985    
986 artak 1481 r1 = b
987     y = Msolve(b)
988 artak 1550 beta1 = bilinearform(y,b)
989 artak 1481
990     if beta1< 0: raise NegativeNorm,"negative norm."
991    
992     # If b = 0 exactly, stop with x = 0.
993     if beta1==0: return x*0.
994    
995     if beta1> 0:
996 artak 1486 beta1 = math.sqrt(beta1)
997 artak 1481
998     #------------------------------------------------------------------
999 artak 1484 # Initialize quantities.
1000 artak 1481 # ------------------------------------------------------------------
1001 artak 1482 iter = 0
1002     Anorm = 0
1003     ynorm = 0
1004 artak 1481 oldb = 0
1005     beta = beta1
1006     dbar = 0
1007     epsln = 0
1008     phibar = beta1
1009     rhs1 = beta1
1010     rhs2 = 0
1011     rnorm = phibar
1012     tnorm2 = 0
1013     ynorm2 = 0
1014     cs = -1
1015     sn = 0
1016     w = b*0.
1017     w2 = b*0.
1018     r2 = r1
1019     eps = 0.0001
1020    
1021     #---------------------------------------------------------------------
1022     # Main iteration loop.
1023     # --------------------------------------------------------------------
1024 artak 1517 while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'): # checks ||r|| < (||A|| ||x||) * TOL
1025 artak 1481
1026     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1027     iter = iter + 1
1028    
1029     #-----------------------------------------------------------------
1030     # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1031     # The general iteration is similar to the case k = 1 with v0 = 0:
1032     #
1033     # p1 = Operator * v1 - beta1 * v0,
1034     # alpha1 = v1'p1,
1035     # q2 = p2 - alpha1 * v1,
1036     # beta2^2 = q2'q2,
1037     # v2 = (1/beta2) q2.
1038     #
1039     # Again, y = betak P vk, where P = C**(-1).
1040     #-----------------------------------------------------------------
1041     s = 1/beta # Normalize previous vector (in y).
1042     v = s*y # v = vk if P = I
1043    
1044     y = Aprod(v)
1045 artak 1465
1046 artak 1481 if iter >= 2:
1047     y = y - (beta/oldb)*r1
1048    
1049     alfa = bilinearform(v,y) # alphak
1050 artak 1550 y += (- alfa/beta)*r2
1051 artak 1481 r1 = r2
1052     r2 = y
1053     y = Msolve(r2)
1054     oldb = beta # oldb = betak
1055 artak 1550 beta = bilinearform(y,r2) # beta = betak+1^2
1056 artak 1481 if beta < 0: raise NegativeNorm,"negative norm."
1057    
1058     beta = math.sqrt( beta )
1059     tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1060    
1061     if iter==1: # Initialize a few things.
1062     gmax = abs( alfa ) # alpha1
1063     gmin = gmax # alpha1
1064    
1065     # Apply previous rotation Qk-1 to get
1066     # [deltak epslnk+1] = [cs sn][dbark 0 ]
1067     # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1068    
1069     oldeps = epsln
1070     delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1071     gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1072     epsln = sn * beta # epsln2 = 0 epslnk+1
1073     dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1074    
1075     # Compute the next plane rotation Qk
1076    
1077     gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1078     gamma = max(gamma,eps)
1079     cs = gbar / gamma # ck
1080     sn = beta / gamma # sk
1081     phi = cs * phibar # phik
1082     phibar = sn * phibar # phibark+1
1083    
1084     # Update x.
1085    
1086     denom = 1/gamma
1087     w1 = w2
1088     w2 = w
1089     w = (v - oldeps*w1 - delta*w2) * denom
1090 artak 1550 x += phi*w
1091 artak 1481
1092     # Go round again.
1093    
1094     gmax = max(gmax,gamma)
1095     gmin = min(gmin,gamma)
1096     z = rhs1 / gamma
1097     ynorm2 = z*z + ynorm2
1098     rhs1 = rhs2 - delta*z
1099     rhs2 = - epsln*z
1100    
1101     # Estimate various norms and test for convergence.
1102    
1103     Anorm = math.sqrt( tnorm2 )
1104     ynorm = math.sqrt( ynorm2 )
1105    
1106     rnorm = phibar
1107    
1108     return x
1109 artak 1489
1110     def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1111    
1112     # TFQMR solver for linear systems
1113     #
1114     #
1115     # initialization
1116     #
1117     errtol = math.sqrt(bilinearform(b,b))
1118     norm_b=errtol
1119     kmax = iter_max
1120     error = []
1121    
1122     if math.sqrt(bilinearform(x,x)) != 0.0:
1123     r = b - Aprod(x)
1124     else:
1125     r = b
1126    
1127     r=Msolve(r)
1128    
1129     u1=0
1130     u2=0
1131     y1=0
1132     y2=0
1133    
1134     w = r
1135     y1 = r
1136     iter = 0
1137     d = 0
1138    
1139     v = Msolve(Aprod(y1))
1140     u1 = v
1141    
1142     theta = 0.0;
1143     eta = 0.0;
1144     tau = math.sqrt(bilinearform(r,r))
1145     error = [ error, tau ]
1146     rho = tau * tau
1147     m=1
1148     #
1149     # TFQMR iteration
1150     #
1151     # while ( iter < kmax-1 ):
1152    
1153 artak 1517 while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1154 artak 1489 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1155    
1156     sigma = bilinearform(r,v)
1157    
1158     if ( sigma == 0.0 ):
1159     raise 'TFQMR breakdown, sigma=0'
1160    
1161    
1162     alpha = rho / sigma
1163    
1164     for j in range(2):
1165     #
1166     # Compute y2 and u2 only if you have to
1167     #
1168     if ( j == 1 ):
1169     y2 = y1 - alpha * v
1170     u2 = Msolve(Aprod(y2))
1171    
1172     m = 2 * (iter+1) - 2 + (j+1)
1173     if j==0:
1174     w = w - alpha * u1
1175     d = y1 + ( theta * theta * eta / alpha ) * d
1176     if j==1:
1177     w = w - alpha * u2
1178     d = y2 + ( theta * theta * eta / alpha ) * d
1179    
1180     theta = math.sqrt(bilinearform(w,w))/ tau
1181     c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1182     tau = tau * theta * c
1183     eta = c * c * alpha
1184     x = x + eta * d
1185     #
1186     # Try to terminate the iteration at each pass through the loop
1187     #
1188     # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1189     # error = [ error, tau ]
1190     # total_iters = iter
1191     # break
1192    
1193    
1194     if ( rho == 0.0 ):
1195     raise 'TFQMR breakdown, rho=0'
1196    
1197    
1198     rhon = bilinearform(r,w)
1199     beta = rhon / rho;
1200     rho = rhon;
1201     y1 = w + beta * y2;
1202     u1 = Msolve(Aprod(y1))
1203     v = u1 + beta * ( u2 + beta * v )
1204     error = [ error, tau ]
1205     total_iters = iter
1206    
1207     iter = iter + 1
1208    
1209     return x
1210    
1211    
1212 artak 1465 #############################################
1213    
1214 gross 1331 class ArithmeticTuple(object):
1215     """
1216     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1217    
1218     example of usage:
1219    
1220     from esys.escript import Data
1221     from numarray import array
1222     a=Data(...)
1223     b=array([1.,4.])
1224     x=ArithmeticTuple(a,b)
1225     y=5.*x
1226    
1227     """
1228     def __init__(self,*args):
1229     """
1230     initialize object with elements args.
1231    
1232     @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1233     """
1234     self.__items=list(args)
1235    
1236     def __len__(self):
1237     """
1238     number of items
1239    
1240     @return: number of items
1241     @rtype: C{int}
1242     """
1243     return len(self.__items)
1244    
1245     def __getitem__(self,index):
1246     """
1247     get an item
1248    
1249     @param index: item to be returned
1250     @type index: C{int}
1251     @return: item with index C{index}
1252     """
1253     return self.__items.__getitem__(index)
1254    
1255     def __mul__(self,other):
1256     """
1257     scaling from the right
1258    
1259     @param other: scaling factor
1260     @type other: C{float}
1261     @return: itemwise self*other
1262     @rtype: L{ArithmeticTuple}
1263     """
1264     out=[]
1265 gross 1896 try:
1266     l=len(other)
1267     if l!=len(self):
1268     raise ValueError,"length of of arguments don't match."
1269     for i in range(l): out.append(self[i]*other[i])
1270     except TypeError:
1271     for i in range(len(self)): out.append(self[i]*other)
1272 gross 1331 return ArithmeticTuple(*tuple(out))
1273    
1274     def __rmul__(self,other):
1275     """
1276     scaling from the left
1277    
1278     @param other: scaling factor
1279     @type other: C{float}
1280     @return: itemwise other*self
1281     @rtype: L{ArithmeticTuple}
1282     """
1283     out=[]
1284 gross 1896 try:
1285     l=len(other)
1286     if l!=len(self):
1287     raise ValueError,"length of of arguments don't match."
1288     for i in range(l): out.append(other[i]*self[i])
1289     except TypeError:
1290     for i in range(len(self)): out.append(other*self[i])
1291 gross 1331 return ArithmeticTuple(*tuple(out))
1292    
1293 artak 1465 def __div__(self,other):
1294     """
1295     dividing from the right
1296    
1297     @param other: scaling factor
1298     @type other: C{float}
1299     @return: itemwise self/other
1300     @rtype: L{ArithmeticTuple}
1301     """
1302 gross 1896 return self*(1/other)
1303 artak 1465
1304     def __rdiv__(self,other):
1305     """
1306     dividing from the left
1307    
1308     @param other: scaling factor
1309     @type other: C{float}
1310     @return: itemwise other/self
1311     @rtype: L{ArithmeticTuple}
1312     """
1313     out=[]
1314 gross 1896 try:
1315     l=len(other)
1316     if l!=len(self):
1317     raise ValueError,"length of of arguments don't match."
1318     for i in range(l): out.append(other[i]/self[i])
1319     except TypeError:
1320     for i in range(len(self)): out.append(other/self[i])
1321 artak 1465 return ArithmeticTuple(*tuple(out))
1322    
1323 gross 1331 def __iadd__(self,other):
1324     """
1325     in-place add of other to self
1326    
1327     @param other: increment
1328     @type other: C{ArithmeticTuple}
1329     """
1330     if len(self) != len(other):
1331     raise ValueError,"tuple length must match."
1332     for i in range(len(self)):
1333     self.__items[i]+=other[i]
1334     return self
1335    
1336 artak 1550 def __add__(self,other):
1337     """
1338     add other to self
1339    
1340     @param other: increment
1341     @type other: C{ArithmeticTuple}
1342     """
1343 gross 1896 out=[]
1344     try:
1345     l=len(other)
1346     if l!=len(self):
1347     raise ValueError,"length of of arguments don't match."
1348     for i in range(l): out.append(self[i]+other[i])
1349     except TypeError:
1350     for i in range(len(self)): out.append(self[i]+other)
1351     return ArithmeticTuple(*tuple(out))
1352 artak 1550
1353     def __sub__(self,other):
1354     """
1355     subtract other from self
1356    
1357     @param other: increment
1358     @type other: C{ArithmeticTuple}
1359     """
1360 gross 1896 out=[]
1361     try:
1362     l=len(other)
1363     if l!=len(self):
1364     raise ValueError,"length of of arguments don't match."
1365     for i in range(l): out.append(self[i]-other[i])
1366     except TypeError:
1367     for i in range(len(self)): out.append(self[i]-other)
1368     return ArithmeticTuple(*tuple(out))
1369 artak 1557
1370     def __isub__(self,other):
1371     """
1372     subtract other from self
1373 artak 1550
1374 artak 1557 @param other: increment
1375     @type other: C{ArithmeticTuple}
1376     """
1377     if len(self) != len(other):
1378     raise ValueError,"tuple length must match."
1379     for i in range(len(self)):
1380     self.__items[i]-=other[i]
1381     return self
1382    
1383 artak 1550 def __neg__(self):
1384     """
1385     negate
1386    
1387     """
1388 gross 1896 out=[]
1389 artak 1550 for i in range(len(self)):
1390 gross 1896 out.append(-self[i])
1391     return ArithmeticTuple(*tuple(out))
1392 artak 1550
1393    
1394 gross 1414 class HomogeneousSaddlePointProblem(object):
1395     """
1396 gross 1878 This provides a framwork for solving linear homogeneous saddle point problem of the form
1397 gross 1414
1398     Av+B^*p=f
1399     Bv =0
1400    
1401     for the unknowns v and p and given operators A and B and given right hand side f.
1402     B^* is the adjoint operator of B is the given inner product.
1403    
1404     """
1405     def __init__(self,**kwargs):
1406     self.setTolerance()
1407     self.setToleranceReductionFactor()
1408    
1409     def initialize(self):
1410     """
1411     initialize the problem (overwrite)
1412     """
1413     pass
1414     def B(self,v):
1415     """
1416     returns Bv (overwrite)
1417     @rtype: equal to the type of p
1418    
1419     @note: boundary conditions on p should be zero!
1420     """
1421     pass
1422    
1423     def inner(self,p0,p1):
1424     """
1425     returns inner product of two element p0 and p1 (overwrite)
1426    
1427     @type p0: equal to the type of p
1428     @type p1: equal to the type of p
1429     @rtype: equal to the type of p
1430     """
1431     pass
1432    
1433     def solve_A(self,u,p):
1434     """
1435     solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1436    
1437     @rtype: equal to the type of v
1438     @note: boundary conditions on v should be zero!
1439     """
1440     pass
1441    
1442     def solve_prec(self,p):
1443     """
1444     provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1445    
1446     @rtype: equal to the type of p
1447     """
1448     pass
1449    
1450     def stoppingcriterium(self,Bv,v,p):
1451     """
1452     returns a True if iteration is terminated. (overwrite)
1453    
1454     @rtype: C{bool}
1455     """
1456     pass
1457    
1458     def __inner(self,p,r):
1459     return self.inner(p,r[1])
1460    
1461 artak 1465 def __inner_p(self,p1,p2):
1462     return self.inner(p1,p2)
1463 artak 1550
1464     def __inner_a(self,a1,a2):
1465     return self.inner_a(a1,a2)
1466 artak 1465
1467 artak 1550 def __inner_a1(self,a1,a2):
1468     return self.inner(a1[1],a2[1])
1469    
1470 gross 1414 def __stoppingcriterium(self,norm_r,r,p):
1471     return self.stoppingcriterium(r[1],r[0],p)
1472    
1473 artak 1519 def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1474     return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1475 artak 1465
1476 gross 1414 def setTolerance(self,tolerance=1.e-8):
1477     self.__tol=tolerance
1478     def getTolerance(self):
1479     return self.__tol
1480     def setToleranceReductionFactor(self,reduction=0.01):
1481     self.__reduction=reduction
1482     def getSubProblemTolerance(self):
1483     return self.__reduction*self.getTolerance()
1484    
1485 artak 1517 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1486 gross 1414 """
1487     solves the saddle point problem using initial guesses v and p.
1488    
1489     @param max_iter: maximum number of iteration steps.
1490     """
1491     self.verbose=verbose
1492     self.show_details=show_details and self.verbose
1493    
1494 gross 1469 # assume p is known: then v=A^-1(f-B^*p)
1495     # which leads to BA^-1B^*p = BA^-1f
1496    
1497 gross 1414 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1498     self.__z=v+self.solve_A(v,p*0)
1499 artak 1557 Bz=self.B(self.__z)
1500 gross 1414 #
1501     # solve BA^-1B^*p = Bz
1502     #
1503     #
1504     #
1505     self.iter=0
1506 artak 1465 if solver=='GMRES':
1507     if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1508 artak 1517 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1509 gross 1467 # solve Au=f-B^*p
1510     # A(u-v)=f-B^*p-Av
1511     # u=v+(u-v)
1512 artak 1465 u=v+self.solve_A(v,p)
1513 artak 1481
1514 artak 1489 if solver=='TFQMR':
1515 artak 1517 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1516     p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1517 artak 1489 # solve Au=f-B^*p
1518     # A(u-v)=f-B^*p-Av
1519     # u=v+(u-v)
1520     u=v+self.solve_A(v,p)
1521    
1522 artak 1481 if solver=='MINRES':
1523     if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1524 artak 1517 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1525 artak 1481 # solve Au=f-B^*p
1526     # A(u-v)=f-B^*p-Av
1527     # u=v+(u-v)
1528     u=v+self.solve_A(v,p)
1529 artak 1488
1530 artak 1550 if solver=='GMRESC':
1531     if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1532 artak 1557 p0=self.solve_prec1(Bz)
1533     #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1534     #p-=alfa
1535 artak 1550 x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)
1536     #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())
1537 artak 1557
1538 artak 1550 # solve Au=f-B^*p
1539     # A(u-v)=f-B^*p-Av
1540     # u=v+(u-v)
1541     p=x[1]
1542 artak 1557 u=v+self.solve_A(v,p)
1543 artak 1550 #p=x[1]
1544     #u=x[0]
1545    
1546 artak 1488 if solver=='PCG':
1547 gross 1552 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1548     #
1549     # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1550     # A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1551 artak 1465 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1552 gross 1467 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1553 artak 1465 u=r[0]
1554 gross 1659 # print "DIFF=",util.integrate(p)
1555 artak 1481
1556 gross 1659 # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1557 gross 1414
1558 artak 1465 return u,p
1559    
1560 gross 1414 def __Msolve(self,r):
1561     return self.solve_prec(r[1])
1562    
1563 artak 1517 def __Msolve2(self,r):
1564 artak 1550 return self.solve_prec(r*1)
1565 artak 1465
1566 artak 1550 def __Mempty(self,r):
1567     return r
1568 artak 1465
1569 artak 1550
1570 gross 1414 def __Aprod(self,p):
1571     # return BA^-1B*p
1572 gross 1552 #solve Av =B^*p as Av =f-Az-B^*(-p)
1573 gross 1469 v=self.solve_A(self.__z,-p)
1574     return ArithmeticTuple(v, self.B(v))
1575 gross 1414
1576 artak 1550 def __Anext(self,x):
1577     # return next v,p
1578     #solve Av =-B^*p as Av =f-Az-B^*p
1579    
1580     pc=x[1]
1581     v=self.solve_A(self.__z,-pc)
1582 artak 1554 p=self.solve_prec1(self.B(v))
1583 artak 1550
1584     return ArithmeticTuple(v,p)
1585    
1586    
1587 artak 1517 def __Aprod2(self,p):
1588 artak 1465 # return BA^-1B*p
1589 gross 1552 #solve Av =B^*p as Av =f-Az-B^*(-p)
1590 gross 1469 v=self.solve_A(self.__z,-p)
1591     return self.B(v)
1592 gross 1414
1593 artak 1519 def __Aprod_Newton(self,p):
1594 artak 1550 # return BA^-1B*p - Bz
1595 artak 1519 #solve Av =-B^*p as Av =f-Az-B^*p
1596     v=self.solve_A(self.__z,-p)
1597     return self.B(v-self.__z)
1598    
1599 artak 1550 def __Aprod_Newton2(self,x):
1600     # return BA^-1B*p - Bz
1601     #solve Av =-B^*p as Av =f-Az-B^*p
1602     pc=x[1]
1603     v=self.solve_A(self.__z,-pc)
1604 artak 1554 p=self.solve_prec1(self.B(v-self.__z))
1605 artak 1550 return ArithmeticTuple(v,p)
1606    
1607 gross 1878
1608     def MaskFromBoundaryTag(domain,*tags):
1609     """
1610 gross 1956 creates a mask on the Solution(domain) function space which one for samples
1611 gross 1878 that touch the boundary tagged by tags.
1612    
1613     usage: m=MaskFromBoundaryTag(domain,"left", "right")
1614    
1615     @param domain: a given domain
1616     @type domain: L{escript.Domain}
1617     @param tags: boundray tags
1618     @type tags: C{str}
1619     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1620     @rtype: L{escript.Data} of rank 0
1621     """
1622     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1623 gross 1956 d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1624 gross 1878 for t in tags: d.setTaggedValue(t,1.)
1625     pde.setValue(y=d)
1626     return util.whereNonZero(pde.getRightHandSide())
1627     #============================================================================================================================================
1628 gross 867 class SaddlePointProblem(object):
1629     """
1630     This implements a solver for a saddlepoint problem
1631    
1632 gross 877 M{f(u,p)=0}
1633     M{g(u)=0}
1634 gross 867
1635     for u and p. The problem is solved with an inexact Uszawa scheme for p:
1636    
1637 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1638 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1639 gross 867
1640     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1641     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'
1642     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1643     in fact the role of a preconditioner.
1644     """
1645     def __init__(self,verbose=False,*args):
1646     """
1647     initializes the problem
1648    
1649 ksteube 990 @param verbose: switches on the printing out some information
1650 gross 867 @type verbose: C{bool}
1651     @note: this method may be overwritten by a particular saddle point problem
1652     """
1653 gross 1659 print "SaddlePointProblem should not be used anymore!"
1654 gross 1107 if not isinstance(verbose,bool):
1655     raise TypeError("verbose needs to be of type bool.")
1656 gross 1106 self.__verbose=verbose
1657 gross 877 self.relaxation=1.
1658 gross 1878 DeprecationWarning("SaddlePointProblem should not be used anymore.")
1659 gross 867
1660     def trace(self,text):
1661     """
1662     prints text if verbose has been set
1663    
1664 ksteube 990 @param text: a text message
1665 gross 867 @type text: C{str}
1666     """
1667 ksteube 1567 if self.__verbose: print "%s: %s"%(str(self),text)
1668 gross 867
1669 gross 873 def solve_f(self,u,p,tol=1.e-8):
1670 gross 867 """
1671     solves
1672    
1673     A_f du = f(u,p)
1674    
1675     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1676    
1677     @param u: current approximation of u
1678     @type u: L{escript.Data}
1679     @param p: current approximation of p
1680     @type p: L{escript.Data}
1681 gross 873 @param tol: tolerance expected for du
1682 gross 867 @type tol: C{float}
1683     @return: increment du
1684     @rtype: L{escript.Data}
1685     @note: this method has to be overwritten by a particular saddle point problem
1686     """
1687     pass
1688    
1689 gross 873 def solve_g(self,u,tol=1.e-8):
1690 gross 867 """
1691     solves
1692    
1693     Q_g dp = g(u)
1694    
1695     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.
1696    
1697     @param u: current approximation of u
1698     @type u: L{escript.Data}
1699 gross 873 @param tol: tolerance expected for dp
1700     @type tol: C{float}
1701 gross 867 @return: increment dp
1702     @rtype: L{escript.Data}
1703     @note: this method has to be overwritten by a particular saddle point problem
1704     """
1705     pass
1706    
1707     def inner(self,p0,p1):
1708     """
1709     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1710     @return: inner product of p0 and p1
1711     @rtype: C{float}
1712     """
1713     pass
1714    
1715 gross 877 subiter_max=3
1716     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1717     """
1718     runs the solver
1719 gross 873
1720 gross 877 @param u0: initial guess for C{u}
1721     @type u0: L{esys.escript.Data}
1722     @param p0: initial guess for C{p}
1723     @type p0: L{esys.escript.Data}
1724     @param tolerance: tolerance for relative error in C{u} and C{p}
1725     @type tolerance: positive C{float}
1726     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1727     @type tolerance_u: positive C{float}
1728     @param iter_max: maximum number of iteration steps.
1729     @type iter_max: C{int}
1730     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1731     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1732     @type accepted_reduction: positive C{float} or C{None}
1733     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1734     @type relaxation: C{float} or C{None}
1735     """
1736     tol=1.e-2
1737     if tolerance_u==None: tolerance_u=tolerance
1738     if not relaxation==None: self.relaxation=relaxation
1739     if accepted_reduction ==None:
1740     angle_limit=0.
1741     elif accepted_reduction>=1.:
1742     angle_limit=0.
1743     else:
1744     angle_limit=util.sqrt(1-accepted_reduction**2)
1745     self.iter=0
1746     u=u0
1747     p=p0
1748     #
1749     # initialize things:
1750     #
1751     converged=False
1752     #
1753     # start loop:
1754     #
1755     # initial search direction is g
1756     #
1757     while not converged :
1758     if self.iter>iter_max:
1759     raise ArithmeticError("no convergence after %s steps."%self.iter)
1760     f_new=self.solve_f(u,p,tol)
1761     norm_f_new = util.Lsup(f_new)
1762     u_new=u-f_new
1763     g_new=self.solve_g(u_new,tol)
1764     self.iter+=1
1765     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1766     if norm_f_new==0. and norm_g_new==0.: return u, p
1767     if self.iter>1 and not accepted_reduction==None:
1768     #
1769     # did we manage to reduce the norm of G? I
1770     # if not we start a backtracking procedure
1771     #
1772     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1773     if norm_g_new > accepted_reduction * norm_g:
1774     sub_iter=0
1775     s=self.relaxation
1776     d=g
1777     g_last=g
1778     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1779     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1780     dg= g_new-g_last
1781     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1782     rad=self.inner(g_new,dg)/self.relaxation
1783     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1784     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1785     if abs(rad) < norm_dg*norm_g_new * angle_limit:
1786     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1787     break
1788     r=self.relaxation
1789     self.relaxation= - rad/norm_dg**2
1790     s+=self.relaxation
1791     #####
1792     # a=g_new+self.relaxation*dg/r
1793     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1794     #####
1795     g_last=g_new
1796     p+=self.relaxation*d
1797     f_new=self.solve_f(u,p,tol)
1798     u_new=u-f_new
1799     g_new=self.solve_g(u_new,tol)
1800     self.iter+=1
1801     norm_f_new = util.Lsup(f_new)
1802     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1803     # print " ",sub_iter," new g norm",norm_g_new
1804     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1805     #
1806     # can we expect reduction of g?
1807     #
1808     # u_last=u_new
1809     sub_iter+=1
1810     self.relaxation=s
1811     #
1812     # check for convergence:
1813     #
1814     norm_u_new = util.Lsup(u_new)
1815     p_new=p+self.relaxation*g_new
1816     norm_p_new = util.sqrt(self.inner(p_new,p_new))
1817 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))
1818 gross 873
1819 gross 877 if self.iter>1:
1820     dg2=g_new-g
1821     df2=f_new-f
1822     norm_dg2=util.sqrt(self.inner(dg2,dg2))
1823     norm_df2=util.Lsup(df2)
1824     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1825     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1826     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1827     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1828     converged=True
1829     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
1830     self.trace("convergence after %s steps."%self.iter)
1831     return u,p

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26