/[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 2719 - (show annotations)
Wed Oct 14 06:38:03 2009 UTC (11 years ago) by gross
File MIME type: text/x-python
File size: 64768 byte(s)
a new Stokes solver added
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, subtol_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 subtol_max: upper bound for inner tolerance
670 :type subtol_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 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
680
681 F=defect(x)
682 fnrm=defect.norm(F)
683 stop_tol=atol + rtol*fnrm
684 subtol=subtol_max
685 if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
686 if verbose: print " tolerance = %e."%subtol
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 subtol_
695 #
696 if iter > 1:
697 rat=fnrm/fnrmo
698 subtol_old=subtol
699 subtol=gamma*rat**2
700 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
701 subtol=max(min(subtol,subtol_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."%subtol
710 try:
711 xc, sub_iter=__FDGMRES(F, defect, x, subtol*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 *A* may depend weakly on *v* and *p*.
1476 """
1477 def __init__(self, **kwargs):
1478 """
1479 initializes the saddle point problem
1480 """
1481 self.resetControlParameters()
1482 self.setTolerance()
1483 self.setAbsoluteTolerance()
1484 def resetControlParameters(self,gamma=0.85, gamma_min=1.e-8,chi_max=0.1, omega_div=0.2, omega_conv=1.1, rtol_min=1.e-7, rtol_max=0.1, chi=1., C_p=1., C_v=1., safety_factor=0.3):
1485 """
1486 sets a control parameter
1487
1488 :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1489 :type gamma: ``float``
1490 :param gamma_min: minimum value for ``gamma``.
1491 :type gamma_min: ``float``
1492 :param chi_max: maximum tolerable converegence rate.
1493 :type chi_max: ``float``
1494 :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1495 :type omega_div: ``float``
1496 :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1497 :type omega_conv: ``float``
1498 :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1499 :type rtol_min: ``float``
1500 :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1501 :type rtol_max: ``float``
1502 :param chi: initial convergence measure.
1503 :type chi: ``float``
1504 :param C_p: initial value for constant to adjust pressure tolerance
1505 :type C_p: ``float``
1506 :param C_v: initial value for constant to adjust velocity tolerance
1507 :type C_v: ``float``
1508 :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1509 :type safety_factor: ``float``
1510 """
1511 self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)
1512
1513 def setControlParameter(self,gamma=None, gamma_min=None ,chi_max=None, omega_div=None, omega_conv=None, rtol_min=None, rtol_max=None, chi=None, C_p=None, C_v=None, safety_factor=None):
1514 """
1515 sets a control parameter
1516
1517 :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1518 :type gamma: ``float``
1519 :param gamma_min: minimum value for ``gamma``.
1520 :type gamma_min: ``float``
1521 :param chi_max: maximum tolerable converegence rate.
1522 :type chi_max: ``float``
1523 :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1524 :type omega_div: ``float``
1525 :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1526 :type omega_conv: ``float``
1527 :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1528 :type rtol_min: ``float``
1529 :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1530 :type rtol_max: ``float``
1531 :param chi: initial convergence measure.
1532 :type chi: ``float``
1533 :param C_p: initial value for constant to adjust pressure tolerance
1534 :type C_p: ``float``
1535 :param C_v: initial value for constant to adjust velocity tolerance
1536 :type C_v: ``float``
1537 :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1538 :type safety_factor: ``float``
1539 """
1540 if not gamma == None:
1541 if gamma<=0 or gamma>=1:
1542 raise ValueError,"Initial gamma needs to be positive and less than 1."
1543 else:
1544 gamma=self.__gamma
1545
1546 if not gamma_min == None:
1547 if gamma_min<=0 or gamma_min>=1:
1548 raise ValueError,"gamma_min needs to be positive and less than 1."
1549 else:
1550 gamma_min = self.__gamma_min
1551
1552 if not chi_max == None:
1553 if chi_max<=0 or chi_max>=1:
1554 raise ValueError,"chi_max needs to be positive and less than 1."
1555 else:
1556 chi_max = self.__chi_max
1557
1558 if not omega_div == None:
1559 if omega_div<=0 or omega_div >=1:
1560 raise ValueError,"omega_div needs to be positive and less than 1."
1561 else:
1562 omega_div=self.__omega_div
1563
1564 if not omega_conv == None:
1565 if omega_conv<1:
1566 raise ValueError,"omega_conv needs to be greater or equal to 1."
1567 else:
1568 omega_conv=self.__omega_conv
1569
1570 if not rtol_min == None:
1571 if rtol_min<=0 or rtol_min>=1:
1572 raise ValueError,"rtol_min needs to be positive and less than 1."
1573 else:
1574 rtol_min=self.__rtol_min
1575
1576 if not rtol_max == None:
1577 if rtol_max<=0 or rtol_max>=1:
1578 raise ValueError,"rtol_max needs to be positive and less than 1."
1579 else:
1580 rtol_max=self.__rtol_max
1581
1582 if not chi == None:
1583 if chi<=0 or chi>1:
1584 raise ValueError,"chi needs to be positive and less than 1."
1585 else:
1586 chi=self.__chi
1587
1588 if not C_p == None:
1589 if C_p<1:
1590 raise ValueError,"C_p need to be greater or equal to 1."
1591 else:
1592 C_p=self.__C_p
1593
1594 if not C_v == None:
1595 if C_v<1:
1596 raise ValueError,"C_v need to be greater or equal to 1."
1597 else:
1598 C_v=self.__C_v
1599
1600 if not safety_factor == None:
1601 if safety_factor<=0 or safety_factor>1:
1602 raise ValueError,"safety_factor need to be between zero and one."
1603 else:
1604 safety_factor=self.__safety_factor
1605
1606 if gamma<gamma_min:
1607 raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)
1608 if rtol_max<=rtol_min:
1609 raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)
1610
1611 self.__gamma = gamma
1612 self.__gamma_min = gamma_min
1613 self.__chi_max = chi_max
1614 self.__omega_div = omega_div
1615 self.__omega_conv = omega_conv
1616 self.__rtol_min = rtol_min
1617 self.__rtol_max = rtol_max
1618 self.__chi = chi
1619 self.__C_p = C_p
1620 self.__C_v = C_v
1621 self.__safety_factor = safety_factor
1622
1623 #=============================================================
1624 def initialize(self):
1625 """
1626 Initializes the problem (overwrite).
1627 """
1628 pass
1629
1630 def inner_pBv(self,p,Bv):
1631 """
1632 Returns inner product of element p and Bv (overwrite).
1633
1634 :param p: a pressure increment
1635 :param Bv: a residual
1636 :return: inner product of element p and Bv
1637 :rtype: ``float``
1638 :note: used if PCG is applied.
1639 """
1640 raise NotImplementedError,"no inner product for p and Bv implemented."
1641
1642 def inner_p(self,p0,p1):
1643 """
1644 Returns inner product of p0 and p1 (overwrite).
1645
1646 :param p0: a pressure
1647 :param p1: a pressure
1648 :return: inner product of p0 and p1
1649 :rtype: ``float``
1650 """
1651 raise NotImplementedError,"no inner product for p implemented."
1652
1653 def norm_v(self,v):
1654 """
1655 Returns the norm of v (overwrite).
1656
1657 :param v: a velovity
1658 :return: norm of v
1659 :rtype: non-negative ``float``
1660 """
1661 raise NotImplementedError,"no norm of v implemented."
1662 def getDV(self, p, v, tol):
1663 """
1664 return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1665
1666 :param p: pressure
1667 :param v: pressure
1668 :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1669 :note: Only *A* may depend on *v* and *p*
1670 """
1671 raise NotImplementedError,"no dv calculation implemented."
1672
1673
1674 def Bv(self,v, tol):
1675 """
1676 Returns Bv with accuracy `tol` (overwrite)
1677
1678 :rtype: equal to the type of p
1679 :note: boundary conditions on p should be zero!
1680 """
1681 raise NotImplementedError, "no operator B implemented."
1682
1683 def norm_Bv(self,Bv):
1684 """
1685 Returns the norm of Bv (overwrite).
1686
1687 :rtype: equal to the type of p
1688 :note: boundary conditions on p should be zero!
1689 """
1690 raise NotImplementedError, "no norm of Bv implemented."
1691
1692 def solve_AinvBt(self,dp, tol):
1693 """
1694 Solves *A dv=B^*dp* with accuracy `tol`
1695
1696 :param dp: a pressure increment
1697 :return: the solution of *A dv=B^*dp*
1698 :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1699 """
1700 raise NotImplementedError,"no operator A implemented."
1701
1702 def solve_prec(self,Bv, tol):
1703 """
1704 Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
1705
1706 :rtype: equal to the type of p
1707 :note: boundary conditions on p should be zero!
1708 """
1709 raise NotImplementedError,"no preconditioner for Schur complement implemented."
1710 #=============================================================
1711 def __Aprod_PCG(self,dp):
1712 dv=self.solve_AinvBt(dp, self.__subtol)
1713 return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1714
1715 def __inner_PCG(self,p,r):
1716 return self.inner_pBv(p,r[1])
1717
1718 def __Msolve_PCG(self,r):
1719 return self.solve_prec(r[1], self.__subtol)
1720 #=============================================================
1721 def __Aprod_GMRES(self,p):
1722 return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1723 def __inner_GMRES(self,p0,p1):
1724 return self.inner_p(p0,p1)
1725
1726 #=============================================================
1727 def norm_p(self,p):
1728 """
1729 calculates the norm of ``p``
1730
1731 :param p: a pressure
1732 :return: the norm of ``p`` using the inner product for pressure
1733 :rtype: ``float``
1734 """
1735 f=self.inner_p(p,p)
1736 if f<0: raise ValueError,"negative pressure norm."
1737 return math.sqrt(f)
1738
1739 def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1740 """
1741 Solves the saddle point problem using initial guesses v and p.
1742
1743 :param v: initial guess for velocity
1744 :param p: initial guess for pressure
1745 :type v: `Data`
1746 :type p: `Data`
1747 :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1748 :param max_iter: maximum number of iteration steps per correction
1749 attempt
1750 :param verbose: if True, shows information on the progress of the
1751 saddlepoint problem solver.
1752 :param iter_restart: restart the iteration after ``iter_restart`` steps
1753 (only used if useUzaw=False)
1754 :type usePCG: ``bool``
1755 :type max_iter: ``int``
1756 :type verbose: ``bool``
1757 :type iter_restart: ``int``
1758 :rtype: ``tuple`` of `Data` objects
1759 """
1760 self.verbose=verbose
1761 rtol=self.getTolerance()
1762 atol=self.getAbsoluteTolerance()
1763 correction_step=0
1764 converged=False
1765 error=None
1766 chi=None
1767 gamma=self.__gamma
1768 C_p=self.__C_p
1769 C_v=self.__C_v
1770 while not converged:
1771 if error== None or chi == None:
1772 gamma_new=gamma/self.__omega_conv
1773 else:
1774 if chi < self.__chi_max:
1775 gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)
1776 else:
1777 gamma_new=gamma*self.__omega_div
1778 gamma=max(gamma_new, self.__gamma_min)
1779 # calculate velocity for current pressure:
1780 rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)
1781 rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)
1782 self.__subtol=rtol_p**2
1783 if self.verbose: print "HomogeneousSaddlePointProblem: step %s: gamma = %e, rtol_v= %e, rtol_p=%e"%(correction_step,gamma,rtol_v,rtol_p)
1784 if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol
1785 # calculate velocity for current pressure: A*dv= F-A*v-B*p
1786 dv1=self.getDV(p,v,rtol_v)
1787 v1=v+dv1
1788 Bv1=self.Bv(v1, self.__subtol)
1789 norm_Bv1=self.norm_Bv(Bv1)
1790 norm_dv1=self.norm_v(dv1)
1791 norm_v1=self.norm_v(v1)
1792 ATOL=norm_v1*rtol+atol
1793 if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)
1794 if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1795 if max(norm_Bv1, norm_dv1) <= ATOL:
1796 converged=True
1797 v=v1
1798 else:
1799 # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1800 if usePCG:
1801 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)
1802 v2=r[0]
1803 Bv2=r[1]
1804 else:
1805 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)
1806 dv2=self.solve_AinvBt(dp, self.__subtol)
1807 v2=v1-dv2
1808 Bv2=self.Bv(v2, self.__subtol)
1809 #
1810 # convergence indicators:
1811 #
1812 norm_v2=self.norm_v(v2)
1813 norm_dv2=self.norm_v(v2-v)
1814 error_new=max(norm_dv2, norm_Bv1)
1815 norm_Bv2=self.norm_Bv(Bv2)
1816 if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)
1817 if error !=None:
1818 chi_new=error_new/error
1819 if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e"%(correction_step,chi_new)
1820 if chi != None:
1821 gamma0=max(gamma, 1-chi/chi_new)
1822 C_p*=gamma0/gamma
1823 C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)
1824 chi=chi_new
1825 error = error_new
1826 correction_step+=1
1827 if correction_step>max_correction_steps:
1828 raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1829 v,p=v2,p+dp
1830 if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1831 return v,p
1832
1833 #========================================================================
1834 def setTolerance(self,tolerance=1.e-4):
1835 """
1836 Sets the relative tolerance for (v,p).
1837
1838 :param tolerance: tolerance to be used
1839 :type tolerance: non-negative ``float``
1840 """
1841 if tolerance<0:
1842 raise ValueError,"tolerance must be positive."
1843 self.__rtol=tolerance
1844
1845 def getTolerance(self):
1846 """
1847 Returns the relative tolerance.
1848
1849 :return: relative tolerance
1850 :rtype: ``float``
1851 """
1852 return self.__rtol
1853
1854 def setAbsoluteTolerance(self,tolerance=0.):
1855 """
1856 Sets the absolute tolerance.
1857
1858 :param tolerance: tolerance to be used
1859 :type tolerance: non-negative ``float``
1860 """
1861 if tolerance<0:
1862 raise ValueError,"tolerance must be non-negative."
1863 self.__atol=tolerance
1864
1865 def getAbsoluteTolerance(self):
1866 """
1867 Returns the absolute tolerance.
1868
1869 :return: absolute tolerance
1870 :rtype: ``float``
1871 """
1872 return self.__atol
1873
1874
1875 def MaskFromBoundaryTag(domain,*tags):
1876 """
1877 Creates a mask on the Solution(domain) function space where the value is
1878 one for samples that touch the boundary tagged by tags.
1879
1880 Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1881
1882 :param domain: domain to be used
1883 :type domain: `escript.Domain`
1884 :param tags: boundary tags
1885 :type tags: ``str``
1886 :return: a mask which marks samples that are touching the boundary tagged
1887 by any of the given tags
1888 :rtype: `escript.Data` of rank 0
1889 """
1890 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1891 d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1892 for t in tags: d.setTaggedValue(t,1.)
1893 pde.setValue(y=d)
1894 return util.whereNonZero(pde.getRightHandSide())
1895
1896 def MaskFromTag(domain,*tags):
1897 """
1898 Creates a mask on the Solution(domain) function space where the value is
1899 one for samples that touch regions tagged by tags.
1900
1901 Usage: m=MaskFromTag(domain, "ham")
1902
1903 :param domain: domain to be used
1904 :type domain: `escript.Domain`
1905 :param tags: boundary tags
1906 :type tags: ``str``
1907 :return: a mask which marks samples that are touching the boundary tagged
1908 by any of the given tags
1909 :rtype: `escript.Data` of rank 0
1910 """
1911 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1912 d=escript.Scalar(0.,escript.Function(domain))
1913 for t in tags: d.setTaggedValue(t,1.)
1914 pde.setValue(Y=d)
1915 return util.whereNonZero(pde.getRightHandSide())
1916
1917

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26