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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26