/[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 1557 - (hide annotations)
Mon May 19 04:44:27 2008 UTC (11 years, 9 months ago) by artak
File MIME type: text/x-python
File size: 57307 byte(s)
some changes to ArithmenticTuple class
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 artak 1550 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 artak 1550 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    
1048     etamax=min(etanew,etamax)
1049     etamax=max(etamax,.5*stop_tol/fnrm)
1050    
1051     return x
1052    
1053 artak 1489 def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1054    
1055     # TFQMR solver for linear systems
1056     #
1057     #
1058     # initialization
1059     #
1060     errtol = math.sqrt(bilinearform(b,b))
1061     norm_b=errtol
1062     kmax = iter_max
1063     error = []
1064    
1065     if math.sqrt(bilinearform(x,x)) != 0.0:
1066     r = b - Aprod(x)
1067     else:
1068     r = b
1069    
1070     r=Msolve(r)
1071    
1072     u1=0
1073     u2=0
1074     y1=0
1075     y2=0
1076    
1077     w = r
1078     y1 = r
1079     iter = 0
1080     d = 0
1081    
1082     v = Msolve(Aprod(y1))
1083     u1 = v
1084    
1085     theta = 0.0;
1086     eta = 0.0;
1087     tau = math.sqrt(bilinearform(r,r))
1088     error = [ error, tau ]
1089     rho = tau * tau
1090     m=1
1091     #
1092     # TFQMR iteration
1093     #
1094     # while ( iter < kmax-1 ):
1095    
1096 artak 1517 while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1097 artak 1489 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1098    
1099     sigma = bilinearform(r,v)
1100    
1101     if ( sigma == 0.0 ):
1102     raise 'TFQMR breakdown, sigma=0'
1103    
1104    
1105     alpha = rho / sigma
1106    
1107     for j in range(2):
1108     #
1109     # Compute y2 and u2 only if you have to
1110     #
1111     if ( j == 1 ):
1112     y2 = y1 - alpha * v
1113     u2 = Msolve(Aprod(y2))
1114    
1115     m = 2 * (iter+1) - 2 + (j+1)
1116     if j==0:
1117     w = w - alpha * u1
1118     d = y1 + ( theta * theta * eta / alpha ) * d
1119     if j==1:
1120     w = w - alpha * u2
1121     d = y2 + ( theta * theta * eta / alpha ) * d
1122    
1123     theta = math.sqrt(bilinearform(w,w))/ tau
1124     c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1125     tau = tau * theta * c
1126     eta = c * c * alpha
1127     x = x + eta * d
1128     #
1129     # Try to terminate the iteration at each pass through the loop
1130     #
1131     # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1132     # error = [ error, tau ]
1133     # total_iters = iter
1134     # break
1135    
1136    
1137     if ( rho == 0.0 ):
1138     raise 'TFQMR breakdown, rho=0'
1139    
1140    
1141     rhon = bilinearform(r,w)
1142     beta = rhon / rho;
1143     rho = rhon;
1144     y1 = w + beta * y2;
1145     u1 = Msolve(Aprod(y1))
1146     v = u1 + beta * ( u2 + beta * v )
1147     error = [ error, tau ]
1148     total_iters = iter
1149    
1150     iter = iter + 1
1151    
1152     return x
1153    
1154    
1155 artak 1465 #############################################
1156    
1157 gross 1331 class ArithmeticTuple(object):
1158     """
1159     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1160    
1161     example of usage:
1162    
1163     from esys.escript import Data
1164     from numarray import array
1165     a=Data(...)
1166     b=array([1.,4.])
1167     x=ArithmeticTuple(a,b)
1168     y=5.*x
1169    
1170     """
1171     def __init__(self,*args):
1172     """
1173     initialize object with elements args.
1174    
1175     @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1176     """
1177     self.__items=list(args)
1178    
1179     def __len__(self):
1180     """
1181     number of items
1182    
1183     @return: number of items
1184     @rtype: C{int}
1185     """
1186     return len(self.__items)
1187    
1188     def __getitem__(self,index):
1189     """
1190     get an item
1191    
1192     @param index: item to be returned
1193     @type index: C{int}
1194     @return: item with index C{index}
1195     """
1196     return self.__items.__getitem__(index)
1197    
1198     def __mul__(self,other):
1199     """
1200     scaling from the right
1201    
1202     @param other: scaling factor
1203     @type other: C{float}
1204     @return: itemwise self*other
1205     @rtype: L{ArithmeticTuple}
1206     """
1207     out=[]
1208 artak 1557 other=1.*other
1209     if isinstance(other,float):
1210     for i in range(len(self)):
1211 gross 1331 out.append(self[i]*other)
1212 artak 1557 else:
1213     for i in range(len(self)):
1214     out.append(self[i]*other[i])
1215 gross 1331 return ArithmeticTuple(*tuple(out))
1216    
1217     def __rmul__(self,other):
1218     """
1219     scaling from the left
1220    
1221     @param other: scaling factor
1222     @type other: C{float}
1223     @return: itemwise other*self
1224     @rtype: L{ArithmeticTuple}
1225     """
1226     out=[]
1227 artak 1557 other=1.*other
1228     if isinstance(other,float):
1229     for i in range(len(self)):
1230 gross 1331 out.append(other*self[i])
1231 artak 1557 else:
1232     for i in range(len(self)):
1233     out.append(other[i]*self[i])
1234 gross 1331 return ArithmeticTuple(*tuple(out))
1235    
1236 artak 1465 #########################
1237     # Added by Artak
1238     #########################
1239     def __div__(self,other):
1240     """
1241     dividing from the right
1242    
1243     @param other: scaling factor
1244     @type other: C{float}
1245     @return: itemwise self/other
1246     @rtype: L{ArithmeticTuple}
1247     """
1248     out=[]
1249 artak 1557 other=1.*other
1250     if isinstance(other,float):
1251     for i in range(len(self)):
1252 gross 1552 out.append(self[i]*(1./other))
1253 artak 1557 else:
1254     for i in range(len(self)):
1255     out.append(self[i]*(1./other[i]))
1256 artak 1465 return ArithmeticTuple(*tuple(out))
1257    
1258     def __rdiv__(self,other):
1259     """
1260     dividing from the left
1261    
1262     @param other: scaling factor
1263     @type other: C{float}
1264     @return: itemwise other/self
1265     @rtype: L{ArithmeticTuple}
1266     """
1267     out=[]
1268 artak 1557 other=1.*other
1269     if isinstance(other,float):
1270     for i in range(len(self)):
1271     out.append(other*(1./self[i]))
1272     else:
1273     for i in range(len(self)):
1274     out.append(other[i]*(1./self[i]))
1275 artak 1465 return ArithmeticTuple(*tuple(out))
1276    
1277     ##########################################33
1278    
1279 gross 1331 def __iadd__(self,other):
1280     """
1281     in-place add of other to self
1282    
1283     @param other: increment
1284     @type other: C{ArithmeticTuple}
1285     """
1286     if len(self) != len(other):
1287     raise ValueError,"tuple length must match."
1288     for i in range(len(self)):
1289     self.__items[i]+=other[i]
1290     return self
1291    
1292 artak 1550 def __add__(self,other):
1293     """
1294     add other to self
1295    
1296     @param other: increment
1297     @type other: C{ArithmeticTuple}
1298     """
1299 artak 1554 # if not isinstance(other,float):
1300     if len(self) != len(other):
1301     raise ValueError,"tuple length must match."
1302     for i in range(len(self)):
1303     self.__items[i]+=other[i]
1304     # else:
1305     # for i in range(len(self)):
1306     # self.__items[i]+=other
1307 artak 1550
1308     return self
1309    
1310     def __sub__(self,other):
1311     """
1312     subtract other from self
1313    
1314     @param other: increment
1315     @type other: C{ArithmeticTuple}
1316     """
1317     if len(self) != len(other):
1318     raise ValueError,"tuple length must match."
1319     for i in range(len(self)):
1320     self.__items[i]-=other[i]
1321     return self
1322 artak 1557
1323     def __isub__(self,other):
1324     """
1325     subtract other from self
1326 artak 1550
1327 artak 1557 @param other: increment
1328     @type other: C{ArithmeticTuple}
1329     """
1330     if len(self) != len(other):
1331     raise ValueError,"tuple length must match."
1332     for i in range(len(self)):
1333     self.__items[i]-=other[i]
1334     return self
1335    
1336 artak 1550 def __neg__(self):
1337     """
1338     negate
1339    
1340     """
1341     for i in range(len(self)):
1342     self.__items[i]=-self.__items[i]
1343     return self
1344    
1345    
1346 gross 1414 class HomogeneousSaddlePointProblem(object):
1347     """
1348     This provides a framwork for solving homogeneous saddle point problem of the form
1349    
1350     Av+B^*p=f
1351     Bv =0
1352    
1353     for the unknowns v and p and given operators A and B and given right hand side f.
1354     B^* is the adjoint operator of B is the given inner product.
1355    
1356     """
1357     def __init__(self,**kwargs):
1358     self.setTolerance()
1359     self.setToleranceReductionFactor()
1360    
1361     def initialize(self):
1362     """
1363     initialize the problem (overwrite)
1364     """
1365     pass
1366     def B(self,v):
1367     """
1368     returns Bv (overwrite)
1369     @rtype: equal to the type of p
1370    
1371     @note: boundary conditions on p should be zero!
1372     """
1373     pass
1374    
1375     def inner(self,p0,p1):
1376     """
1377     returns inner product of two element p0 and p1 (overwrite)
1378    
1379     @type p0: equal to the type of p
1380     @type p1: equal to the type of p
1381     @rtype: C{float}
1382    
1383     @rtype: equal to the type of p
1384     """
1385     pass
1386    
1387     def solve_A(self,u,p):
1388     """
1389     solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1390    
1391     @rtype: equal to the type of v
1392     @note: boundary conditions on v should be zero!
1393     """
1394     pass
1395    
1396     def solve_prec(self,p):
1397     """
1398     provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1399    
1400     @rtype: equal to the type of p
1401     """
1402     pass
1403    
1404     def stoppingcriterium(self,Bv,v,p):
1405     """
1406     returns a True if iteration is terminated. (overwrite)
1407    
1408     @rtype: C{bool}
1409     """
1410     pass
1411    
1412     def __inner(self,p,r):
1413     return self.inner(p,r[1])
1414    
1415 artak 1465 def __inner_p(self,p1,p2):
1416     return self.inner(p1,p2)
1417 artak 1550
1418     def __inner_a(self,a1,a2):
1419     return self.inner_a(a1,a2)
1420 artak 1465
1421 artak 1550 def __inner_a1(self,a1,a2):
1422     return self.inner(a1[1],a2[1])
1423    
1424 gross 1414 def __stoppingcriterium(self,norm_r,r,p):
1425     return self.stoppingcriterium(r[1],r[0],p)
1426    
1427 artak 1519 def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1428     return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1429 artak 1465
1430 gross 1414 def setTolerance(self,tolerance=1.e-8):
1431     self.__tol=tolerance
1432     def getTolerance(self):
1433     return self.__tol
1434     def setToleranceReductionFactor(self,reduction=0.01):
1435     self.__reduction=reduction
1436     def getSubProblemTolerance(self):
1437     return self.__reduction*self.getTolerance()
1438    
1439 artak 1517 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1440 gross 1414 """
1441     solves the saddle point problem using initial guesses v and p.
1442    
1443     @param max_iter: maximum number of iteration steps.
1444     """
1445     self.verbose=verbose
1446     self.show_details=show_details and self.verbose
1447    
1448 gross 1469 # assume p is known: then v=A^-1(f-B^*p)
1449     # which leads to BA^-1B^*p = BA^-1f
1450    
1451 gross 1414 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1452     self.__z=v+self.solve_A(v,p*0)
1453 artak 1557 Bz=self.B(self.__z)
1454 gross 1414 #
1455     # solve BA^-1B^*p = Bz
1456     #
1457     #
1458     #
1459     self.iter=0
1460 artak 1465 if solver=='GMRES':
1461     if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1462 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)
1463 gross 1467 # solve Au=f-B^*p
1464     # A(u-v)=f-B^*p-Av
1465     # u=v+(u-v)
1466 artak 1465 u=v+self.solve_A(v,p)
1467 artak 1481
1468 artak 1519 if solver=='NewtonGMRES':
1469     if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1470 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())
1471 artak 1519 # solve Au=f-B^*p
1472     # A(u-v)=f-B^*p-Av
1473     # u=v+(u-v)
1474     u=v+self.solve_A(v,p)
1475    
1476    
1477 artak 1489 if solver=='TFQMR':
1478 artak 1517 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1479     p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1480 artak 1489 # solve Au=f-B^*p
1481     # A(u-v)=f-B^*p-Av
1482     # u=v+(u-v)
1483     u=v+self.solve_A(v,p)
1484    
1485 artak 1481 if solver=='MINRES':
1486     if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1487 artak 1517 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1488 artak 1481 # solve Au=f-B^*p
1489     # A(u-v)=f-B^*p-Av
1490     # u=v+(u-v)
1491     u=v+self.solve_A(v,p)
1492 artak 1488
1493 artak 1550 if solver=='GMRESC':
1494     if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1495 artak 1557 p0=self.solve_prec1(Bz)
1496     #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1497     #p-=alfa
1498 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)
1499     #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())
1500 artak 1557
1501 artak 1550 # solve Au=f-B^*p
1502     # A(u-v)=f-B^*p-Av
1503     # u=v+(u-v)
1504     p=x[1]
1505 artak 1557 u=v+self.solve_A(v,p)
1506 artak 1550 #p=x[1]
1507     #u=x[0]
1508    
1509 artak 1488 if solver=='PCG':
1510 gross 1552 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1511     #
1512     # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1513     # A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1514 artak 1465 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1515 gross 1467 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1516 artak 1465 u=r[0]
1517 gross 1552 print "DIFF=",util.integrate(p)
1518 artak 1481
1519 artak 1475 print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1520 gross 1414
1521 artak 1465 return u,p
1522    
1523 gross 1414 def __Msolve(self,r):
1524     return self.solve_prec(r[1])
1525    
1526 artak 1517 def __Msolve2(self,r):
1527 artak 1550 return self.solve_prec(r*1)
1528 artak 1465
1529 artak 1550 def __Mempty(self,r):
1530     return r
1531 artak 1465
1532 artak 1550
1533 gross 1414 def __Aprod(self,p):
1534     # return BA^-1B*p
1535 gross 1552 #solve Av =B^*p as Av =f-Az-B^*(-p)
1536 gross 1469 v=self.solve_A(self.__z,-p)
1537     return ArithmeticTuple(v, self.B(v))
1538 gross 1414
1539 artak 1550 def __Anext(self,x):
1540     # return next v,p
1541     #solve Av =-B^*p as Av =f-Az-B^*p
1542    
1543     pc=x[1]
1544     v=self.solve_A(self.__z,-pc)
1545 artak 1554 p=self.solve_prec1(self.B(v))
1546 artak 1550
1547     return ArithmeticTuple(v,p)
1548    
1549    
1550 artak 1517 def __Aprod2(self,p):
1551 artak 1465 # return BA^-1B*p
1552 gross 1552 #solve Av =B^*p as Av =f-Az-B^*(-p)
1553 gross 1469 v=self.solve_A(self.__z,-p)
1554     return self.B(v)
1555 gross 1414
1556 artak 1519 def __Aprod_Newton(self,p):
1557 artak 1550 # return BA^-1B*p - Bz
1558 artak 1519 #solve Av =-B^*p as Av =f-Az-B^*p
1559     v=self.solve_A(self.__z,-p)
1560     return self.B(v-self.__z)
1561    
1562 artak 1550 def __Aprod_Newton2(self,x):
1563     # return BA^-1B*p - Bz
1564     #solve Av =-B^*p as Av =f-Az-B^*p
1565     pc=x[1]
1566     v=self.solve_A(self.__z,-pc)
1567 artak 1554 p=self.solve_prec1(self.B(v-self.__z))
1568 artak 1550 return ArithmeticTuple(v,p)
1569    
1570 gross 867 class SaddlePointProblem(object):
1571     """
1572     This implements a solver for a saddlepoint problem
1573    
1574 gross 877 M{f(u,p)=0}
1575     M{g(u)=0}
1576 gross 867
1577     for u and p. The problem is solved with an inexact Uszawa scheme for p:
1578    
1579 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1580 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1581 gross 867
1582     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1583     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'
1584     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1585     in fact the role of a preconditioner.
1586     """
1587     def __init__(self,verbose=False,*args):
1588     """
1589     initializes the problem
1590    
1591 ksteube 990 @param verbose: switches on the printing out some information
1592 gross 867 @type verbose: C{bool}
1593     @note: this method may be overwritten by a particular saddle point problem
1594     """
1595 gross 1107 if not isinstance(verbose,bool):
1596     raise TypeError("verbose needs to be of type bool.")
1597 gross 1106 self.__verbose=verbose
1598 gross 877 self.relaxation=1.
1599 gross 867
1600     def trace(self,text):
1601     """
1602     prints text if verbose has been set
1603    
1604 ksteube 990 @param text: a text message
1605 gross 867 @type text: C{str}
1606     """
1607 artak 1550 if self.__verbose: print "#s: #s"%(str(self),text)
1608 gross 867
1609 gross 873 def solve_f(self,u,p,tol=1.e-8):
1610 gross 867 """
1611     solves
1612    
1613     A_f du = f(u,p)
1614    
1615     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1616    
1617     @param u: current approximation of u
1618     @type u: L{escript.Data}
1619     @param p: current approximation of p
1620     @type p: L{escript.Data}
1621 gross 873 @param tol: tolerance expected for du
1622 gross 867 @type tol: C{float}
1623     @return: increment du
1624     @rtype: L{escript.Data}
1625     @note: this method has to be overwritten by a particular saddle point problem
1626     """
1627     pass
1628    
1629 gross 873 def solve_g(self,u,tol=1.e-8):
1630 gross 867 """
1631     solves
1632    
1633     Q_g dp = g(u)
1634    
1635     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.
1636    
1637     @param u: current approximation of u
1638     @type u: L{escript.Data}
1639 gross 873 @param tol: tolerance expected for dp
1640     @type tol: C{float}
1641 gross 867 @return: increment dp
1642     @rtype: L{escript.Data}
1643     @note: this method has to be overwritten by a particular saddle point problem
1644     """
1645     pass
1646    
1647     def inner(self,p0,p1):
1648     """
1649     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1650     @return: inner product of p0 and p1
1651     @rtype: C{float}
1652     """
1653     pass
1654    
1655 gross 877 subiter_max=3
1656     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1657     """
1658     runs the solver
1659 gross 873
1660 gross 877 @param u0: initial guess for C{u}
1661     @type u0: L{esys.escript.Data}
1662     @param p0: initial guess for C{p}
1663     @type p0: L{esys.escript.Data}
1664     @param tolerance: tolerance for relative error in C{u} and C{p}
1665     @type tolerance: positive C{float}
1666     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1667     @type tolerance_u: positive C{float}
1668     @param iter_max: maximum number of iteration steps.
1669     @type iter_max: C{int}
1670     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1671     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1672     @type accepted_reduction: positive C{float} or C{None}
1673     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1674     @type relaxation: C{float} or C{None}
1675     """
1676     tol=1.e-2
1677     if tolerance_u==None: tolerance_u=tolerance
1678     if not relaxation==None: self.relaxation=relaxation
1679     if accepted_reduction ==None:
1680     angle_limit=0.
1681     elif accepted_reduction>=1.:
1682     angle_limit=0.
1683     else:
1684     angle_limit=util.sqrt(1-accepted_reduction**2)
1685     self.iter=0
1686     u=u0
1687     p=p0
1688     #
1689     # initialize things:
1690     #
1691     converged=False
1692     #
1693     # start loop:
1694     #
1695     # initial search direction is g
1696     #
1697     while not converged :
1698     if self.iter>iter_max:
1699     raise ArithmeticError("no convergence after %s steps."%self.iter)
1700     f_new=self.solve_f(u,p,tol)
1701     norm_f_new = util.Lsup(f_new)
1702     u_new=u-f_new
1703     g_new=self.solve_g(u_new,tol)
1704     self.iter+=1
1705     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1706     if norm_f_new==0. and norm_g_new==0.: return u, p
1707     if self.iter>1 and not accepted_reduction==None:
1708     #
1709     # did we manage to reduce the norm of G? I
1710     # if not we start a backtracking procedure
1711     #
1712     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1713     if norm_g_new > accepted_reduction * norm_g:
1714     sub_iter=0
1715     s=self.relaxation
1716     d=g
1717     g_last=g
1718     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1719     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1720     dg= g_new-g_last
1721     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1722     rad=self.inner(g_new,dg)/self.relaxation
1723     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1724     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1725     if abs(rad) < norm_dg*norm_g_new * angle_limit:
1726     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1727     break
1728     r=self.relaxation
1729     self.relaxation= - rad/norm_dg**2
1730     s+=self.relaxation
1731     #####
1732     # a=g_new+self.relaxation*dg/r
1733     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1734     #####
1735     g_last=g_new
1736     p+=self.relaxation*d
1737     f_new=self.solve_f(u,p,tol)
1738     u_new=u-f_new
1739     g_new=self.solve_g(u_new,tol)
1740     self.iter+=1
1741     norm_f_new = util.Lsup(f_new)
1742     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1743     # print " ",sub_iter," new g norm",norm_g_new
1744     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1745     #
1746     # can we expect reduction of g?
1747     #
1748     # u_last=u_new
1749     sub_iter+=1
1750     self.relaxation=s
1751     #
1752     # check for convergence:
1753     #
1754     norm_u_new = util.Lsup(u_new)
1755     p_new=p+self.relaxation*g_new
1756     norm_p_new = util.sqrt(self.inner(p_new,p_new))
1757 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))
1758 gross 873
1759 gross 877 if self.iter>1:
1760     dg2=g_new-g
1761     df2=f_new-f
1762     norm_dg2=util.sqrt(self.inner(dg2,dg2))
1763     norm_df2=util.Lsup(df2)
1764     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1765     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1766     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1767     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1768     converged=True
1769     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
1770     self.trace("convergence after %s steps."%self.iter)
1771     return u,p
1772     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1773     # tol=1.e-2
1774     # iter=0
1775     # converged=False
1776     # u=u0*1.
1777     # p=p0*1.
1778     # while not converged and iter<iter_max:
1779     # du=self.solve_f(u,p,tol)
1780     # u-=du
1781     # norm_du=util.Lsup(du)
1782     # norm_u=util.Lsup(u)
1783     #
1784     # dp=self.relaxation*self.solve_g(u,tol)
1785     # p+=dp
1786     # norm_dp=util.sqrt(self.inner(dp,dp))
1787     # norm_p=util.sqrt(self.inner(p,p))
1788     # 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)
1789     # iter+=1
1790     #
1791     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1792     # if converged:
1793     # print "convergence after %s steps."%iter
1794     # else:
1795     # raise ArithmeticError("no convergence after %s steps."%iter)
1796     #
1797     # return u,p
1798 gross 873
1799 ksteube 1312 def MaskFromBoundaryTag(function_space,*tags):
1800     """
1801     create a mask on the given function space which one for samples
1802     that touch the boundary tagged by tags.
1803    
1804     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1805    
1806     @param function_space: a given function space
1807     @type function_space: L{escript.FunctionSpace}
1808     @param tags: boundray tags
1809     @type tags: C{str}
1810     @return: a mask which marks samples used by C{function_space} that are touching the
1811     boundary tagged by any of the given tags.
1812     @rtype: L{escript.Data} of rank 0
1813     """
1814     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1815     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1816     for t in tags: d.setTaggedValue(t,1.)
1817     pde.setValue(y=d)
1818     out=util.whereNonZero(pde.getRightHandSide())
1819     if out.getFunctionSpace() == function_space:
1820     return out
1821     else:
1822     return util.whereNonZero(util.interpolate(out,function_space))
1823    
1824 gross 1414
1825 artak 1465

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26