/[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 3975 - (show annotations)
Thu Sep 20 01:54:06 2012 UTC (7 years, 1 month ago) by caltinay
File MIME type: text/x-python
File size: 64550 byte(s)
Merged symbolic branch into trunk. Curious what daniel and spartacus have to
say...

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26