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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26