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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26