/[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 2954 - (hide annotations)
Fri Feb 26 06:09:47 2010 UTC (10 years ago) by gross
File MIME type: text/x-python
File size: 63508 byte(s)
bug introduced by a smart programmer to the Locator has been fixed?

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26