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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26