/[escript]/trunk/escriptcore/py_src/pdetools.py
ViewVC logotype

Contents of /trunk/escriptcore/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26