/[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 1809 - (hide annotations)
Thu Sep 25 06:43:44 2008 UTC (10 years, 10 months ago) by ksteube
File MIME type: text/x-python
File size: 57088 byte(s)
Copyright updated in all python files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26