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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26