/[escript]/branches/symbolic_from_3470/escript/py_src/pdetools.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 64117 byte(s)
Fast forward to latest trunk revision 3788.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26