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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26