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

Annotation of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3975 - (hide annotations)
Thu Sep 20 01:54:06 2012 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/x-python
File size: 64550 byte(s)
Merged symbolic branch into trunk. Curious what daniel and spartacus have to
say...

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26