/[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 1639 - (hide annotations)
Mon Jul 14 08:55:25 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 57296 byte(s)


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26