/[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 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (10 years ago) by jfenwick
File MIME type: text/x-python
File size: 63520 byte(s)
Don't panic.
Updating copyright stamps

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26