/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Contents of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3675 - (show annotations)
Thu Nov 17 00:53:38 2011 UTC (8 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 64068 byte(s)
pasowrap joins the trunk.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26