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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26