/[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 2344 - (hide annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 64339 byte(s)
Change __url__ to launchpad site

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26