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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26