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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26