/[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 4446 - (hide annotations)
Tue Jun 11 04:00:15 2013 UTC (6 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 67323 byte(s)
This should fix the py3 error...

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26