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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26