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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26