/[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 1489 - (hide annotations)
Mon Apr 14 04:29:30 2008 UTC (11 years, 4 months ago) by artak
File MIME type: text/x-python
File size: 47751 byte(s)
Initial version of TFQMR method
1 ksteube 1312 #
2 jgs 121 # $Id$
3 ksteube 1312 #
4     #######################################################
5     #
6     # Copyright 2003-2007 by ACceSS MNRF
7     # Copyright 2007 by University of Queensland
8     #
9     # http://esscc.uq.edu.au
10     # Primary Business: Queensland, Australia
11     # Licensed under the Open Software License version 3.0
12     # http://www.opensource.org/licenses/osl-3.0.php
13     #
14     #######################################################
15     #
16 jgs 121
17     """
18 jgs 149 Provides some tools related to PDEs.
19 jgs 121
20 jgs 149 Currently includes:
21     - Projector - to project a discontinuous
22 gross 351 - Locator - to trace values in data objects at a certain location
23     - TimeIntegrationManager - to handel extraplotion in time
24 gross 867 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
25 gross 637
26     @var __author__: name of author
27     @var __copyright__: copyrights
28     @var __license__: licence agreement
29     @var __url__: url entry point on documentation
30     @var __version__: version
31     @var __date__: date of the version
32 jgs 121 """
33    
34 gross 637 __author__="Lutz Gross, l.gross@uq.edu.au"
35 elspeth 609 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
36     http://www.access.edu.au
37     Primary Business: Queensland, Australia"""
38 elspeth 614 __license__="""Licensed under the Open Software License version 3.0
39     http://www.opensource.org/licenses/osl-3.0.php"""
40 gross 637 __url__="http://www.iservo.edu.au/esys"
41     __version__="$Revision$"
42     __date__="$Date$"
43 elspeth 609
44 gross 637
45 jgs 149 import escript
46     import linearPDEs
47 jgs 121 import numarray
48 jgs 149 import util
49 ksteube 1312 import math
50 jgs 121
51 artak 1465 ##### Added by Artak
52 gross 1467 # from Numeric import zeros,Int,Float64
53 artak 1465 ###################################
54    
55    
56 gross 351 class TimeIntegrationManager:
57     """
58     a simple mechanism to manage time dependend values.
59    
60 gross 720 typical usage is::
61 gross 351
62 gross 720 dt=0.1 # time increment
63     tm=TimeIntegrationManager(inital_value,p=1)
64     while t<1.
65     v_guess=tm.extrapolate(dt) # extrapolate to t+dt
66     v=...
67     tm.checkin(dt,v)
68     t+=dt
69 gross 351
70 gross 720 @note: currently only p=1 is supported.
71 gross 351 """
72     def __init__(self,*inital_values,**kwargs):
73     """
74     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
75     """
76     if kwargs.has_key("p"):
77     self.__p=kwargs["p"]
78     else:
79     self.__p=1
80     if kwargs.has_key("time"):
81     self.__t=kwargs["time"]
82     else:
83     self.__t=0.
84     self.__v_mem=[inital_values]
85     self.__order=0
86     self.__dt_mem=[]
87     self.__num_val=len(inital_values)
88    
89     def getTime(self):
90     return self.__t
91 gross 396 def getValue(self):
92 gross 409 out=self.__v_mem[0]
93     if len(out)==1:
94     return out[0]
95     else:
96     return out
97    
98 gross 351 def checkin(self,dt,*values):
99     """
100     adds new values to the manager. the p+1 last value get lost
101     """
102     o=min(self.__order+1,self.__p)
103     self.__order=min(self.__order+1,self.__p)
104     v_mem_new=[values]
105     dt_mem_new=[dt]
106     for i in range(o-1):
107     v_mem_new.append(self.__v_mem[i])
108     dt_mem_new.append(self.__dt_mem[i])
109     v_mem_new.append(self.__v_mem[o-1])
110     self.__order=o
111     self.__v_mem=v_mem_new
112     self.__dt_mem=dt_mem_new
113     self.__t+=dt
114    
115     def extrapolate(self,dt):
116     """
117     extrapolates to dt forward in time.
118     """
119     if self.__order==0:
120     out=self.__v_mem[0]
121     else:
122     out=[]
123     for i in range(self.__num_val):
124     out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
125    
126     if len(out)==0:
127     return None
128     elif len(out)==1:
129     return out[0]
130     else:
131     return out
132 gross 396
133 gross 351
134 jgs 121 class Projector:
135 jgs 149 """
136     The Projector is a factory which projects a discontiuous function onto a
137     continuous function on the a given domain.
138     """
139 jgs 121 def __init__(self, domain, reduce = True, fast=True):
140     """
141 jgs 149 Create a continuous function space projector for a domain.
142 jgs 121
143 jgs 149 @param domain: Domain of the projection.
144     @param reduce: Flag to reduce projection order (default is True)
145     @param fast: Flag to use a fast method based on matrix lumping (default is true)
146 jgs 121 """
147 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
148 jgs 148 if fast:
149 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
150 jgs 121 self.__pde.setSymmetryOn()
151     self.__pde.setReducedOrderTo(reduce)
152     self.__pde.setValue(D = 1.)
153 ksteube 1312 return
154 jgs 121
155     def __call__(self, input_data):
156     """
157 jgs 149 Projects input_data onto a continuous function
158 jgs 121
159 jgs 149 @param input_data: The input_data to be projected.
160 jgs 121 """
161 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
162 gross 1122 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
163 jgs 121 if input_data.getRank()==0:
164     self.__pde.setValue(Y = input_data)
165     out=self.__pde.getSolution()
166     elif input_data.getRank()==1:
167     for i0 in range(input_data.getShape()[0]):
168     self.__pde.setValue(Y = input_data[i0])
169     out[i0]=self.__pde.getSolution()
170     elif input_data.getRank()==2:
171     for i0 in range(input_data.getShape()[0]):
172     for i1 in range(input_data.getShape()[1]):
173     self.__pde.setValue(Y = input_data[i0,i1])
174     out[i0,i1]=self.__pde.getSolution()
175     elif input_data.getRank()==3:
176     for i0 in range(input_data.getShape()[0]):
177     for i1 in range(input_data.getShape()[1]):
178     for i2 in range(input_data.getShape()[2]):
179     self.__pde.setValue(Y = input_data[i0,i1,i2])
180     out[i0,i1,i2]=self.__pde.getSolution()
181     else:
182     for i0 in range(input_data.getShape()[0]):
183     for i1 in range(input_data.getShape()[1]):
184     for i2 in range(input_data.getShape()[2]):
185     for i3 in range(input_data.getShape()[3]):
186     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
187     out[i0,i1,i2,i3]=self.__pde.getSolution()
188     return out
189    
190 gross 525 class NoPDE:
191     """
192     solves the following problem for u:
193 jgs 121
194 gross 525 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
195    
196     with constraint
197    
198     M{u[j]=r[j]} where M{q[j]>0}
199    
200     where D, Y, r and q are given functions of rank 1.
201    
202     In the case of scalars this takes the form
203    
204     M{D*u=Y}
205    
206     with constraint
207    
208     M{u=r} where M{q>0}
209    
210     where D, Y, r and q are given scalar functions.
211    
212     The constraint is overwriting any other condition.
213    
214 gross 720 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
215     that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
216     thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
217 gross 525 """
218     def __init__(self,domain,D=None,Y=None,q=None,r=None):
219     """
220     initialize the problem
221    
222     @param domain: domain of the PDE.
223     @type domain: L{Domain}
224     @param D: coefficient of the solution.
225 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
226 gross 525 @param Y: right hand side
227 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
228 gross 525 @param q: location of constraints
229 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
230 gross 525 @param r: value of solution at locations of constraints
231 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
232 gross 525 """
233     self.__domain=domain
234     self.__D=D
235     self.__Y=Y
236     self.__q=q
237     self.__r=r
238     self.__u=None
239     self.__function_space=escript.Solution(self.__domain)
240     def setReducedOn(self):
241     """
242     sets the L{FunctionSpace} of the solution to L{ReducedSolution}
243     """
244     self.__function_space=escript.ReducedSolution(self.__domain)
245     self.__u=None
246    
247     def setReducedOff(self):
248     """
249     sets the L{FunctionSpace} of the solution to L{Solution}
250     """
251     self.__function_space=escript.Solution(self.__domain)
252     self.__u=None
253    
254     def setValue(self,D=None,Y=None,q=None,r=None):
255     """
256     assigns values to the parameters.
257    
258     @param D: coefficient of the solution.
259 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
260 gross 525 @param Y: right hand side
261 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
262 gross 525 @param q: location of constraints
263 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
264 gross 525 @param r: value of solution at locations of constraints
265 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
266 gross 525 """
267     if not D==None:
268     self.__D=D
269     self.__u=None
270     if not Y==None:
271     self.__Y=Y
272     self.__u=None
273     if not q==None:
274     self.__q=q
275     self.__u=None
276     if not r==None:
277     self.__r=r
278     self.__u=None
279    
280     def getSolution(self):
281     """
282     returns the solution
283    
284     @return: the solution of the problem
285     @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
286     """
287     if self.__u==None:
288     if self.__D==None:
289     raise ValueError,"coefficient D is undefined"
290     D=escript.Data(self.__D,self.__function_space)
291     if D.getRank()>1:
292     raise ValueError,"coefficient D must have rank 0 or 1"
293     if self.__Y==None:
294     self.__u=escript.Data(0.,D.getShape(),self.__function_space)
295     else:
296     self.__u=util.quotient(self.__Y,D)
297     if not self.__q==None:
298     q=util.wherePositive(escript.Data(self.__q,self.__function_space))
299     self.__u*=(1.-q)
300     if not self.__r==None: self.__u+=q*self.__r
301     return self.__u
302    
303 jgs 147 class Locator:
304     """
305 jgs 149 Locator provides access to the values of data objects at a given
306     spatial coordinate x.
307    
308     In fact, a Locator object finds the sample in the set of samples of a
309     given function space or domain where which is closest to the given
310     point x.
311 jgs 147 """
312    
313     def __init__(self,where,x=numarray.zeros((3,))):
314 jgs 149 """
315     Initializes a Locator to access values in Data objects on the Doamin
316     or FunctionSpace where for the sample point which
317     closest to the given point x.
318 gross 880
319     @param where: function space
320     @type where: L{escript.FunctionSpace}
321     @param x: coefficient of the solution.
322     @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
323 jgs 149 """
324     if isinstance(where,escript.FunctionSpace):
325 jgs 147 self.__function_space=where
326 jgs 121 else:
327 jgs 149 self.__function_space=escript.ContinuousFunction(where)
328 gross 880 if isinstance(x, list):
329     self.__id=[]
330     for p in x:
331 gross 921 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
332 gross 880 else:
333 gross 921 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
334 jgs 121
335 jgs 147 def __str__(self):
336 jgs 149 """
337     Returns the coordinates of the Locator as a string.
338     """
339 gross 880 x=self.getX()
340     if instance(x,list):
341     out="["
342     first=True
343     for xx in x:
344     if not first:
345     out+=","
346     else:
347     first=False
348     out+=str(xx)
349     out+="]>"
350     else:
351     out=str(x)
352     return out
353 jgs 121
354 gross 880 def getX(self):
355     """
356     Returns the exact coordinates of the Locator.
357     """
358     return self(self.getFunctionSpace().getX())
359    
360 jgs 147 def getFunctionSpace(self):
361 jgs 149 """
362     Returns the function space of the Locator.
363     """
364 jgs 147 return self.__function_space
365    
366 gross 880 def getId(self,item=None):
367 jgs 149 """
368     Returns the identifier of the location.
369     """
370 gross 880 if item == None:
371     return self.__id
372     else:
373     if isinstance(self.__id,list):
374     return self.__id[item]
375     else:
376     return self.__id
377 jgs 121
378    
379 jgs 147 def __call__(self,data):
380 jgs 149 """
381     Returns the value of data at the Locator of a Data object otherwise
382     the object is returned.
383     """
384 jgs 147 return self.getValue(data)
385 jgs 121
386 jgs 147 def getValue(self,data):
387 jgs 149 """
388     Returns the value of data at the Locator if data is a Data object
389     otherwise the object is returned.
390     """
391     if isinstance(data,escript.Data):
392 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
393 gross 880 dat=data
394 jgs 147 else:
395 gross 880 dat=data.interpolate(self.getFunctionSpace())
396     id=self.getId()
397     r=data.getRank()
398     if isinstance(id,list):
399     out=[]
400     for i in id:
401 gross 921 o=data.getValueOfGlobalDataPoint(*i)
402 gross 880 if data.getRank()==0:
403     out.append(o[0])
404     else:
405     out.append(o)
406     return out
407 jgs 147 else:
408 gross 921 out=data.getValueOfGlobalDataPoint(*id)
409 gross 880 if data.getRank()==0:
410     return out[0]
411     else:
412     return out
413 jgs 147 else:
414     return data
415 jgs 149
416 ksteube 1312 class SolverSchemeException(Exception):
417     """
418     exceptions thrown by solvers
419     """
420     pass
421    
422     class IndefinitePreconditioner(SolverSchemeException):
423     """
424     the preconditioner is not positive definite.
425     """
426     pass
427     class MaxIterReached(SolverSchemeException):
428     """
429     maxium number of iteration steps is reached.
430     """
431     pass
432     class IterationBreakDown(SolverSchemeException):
433     """
434     iteration scheme econouters an incurable breakdown.
435     """
436     pass
437     class NegativeNorm(SolverSchemeException):
438     """
439     a norm calculation returns a negative norm.
440     """
441     pass
442    
443 gross 1330 class IterationHistory(object):
444 ksteube 1312 """
445 gross 1330 The IterationHistory class is used to define a stopping criterium. It keeps track of the
446     residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
447     a given tolerance.
448     """
449     def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
450     """
451     Initialization
452    
453     @param tolerance: tolerance
454     @type tolerance: positive C{float}
455     @param verbose: switches on the printing out some information
456     @type verbose: C{bool}
457     """
458     if not tolerance>0.:
459     raise ValueError,"tolerance needs to be positive."
460     self.tolerance=tolerance
461     self.verbose=verbose
462     self.history=[]
463     def stoppingcriterium(self,norm_r,r,x):
464     """
465     returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]} is the residual norm at the first call.
466    
467    
468     @param norm_r: current residual norm
469     @type norm_r: non-negative C{float}
470     @param r: current residual (not used)
471     @param x: current solution approximation (not used)
472     @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
473     @rtype: C{bool}
474    
475     """
476     self.history.append(norm_r)
477 artak 1489 if self.verbose: print "iter: #s: inner(rhat,r) = #e"#(len(self.history)-1, self.history[-1])
478 gross 1330 return self.history[-1]<=self.tolerance * self.history[0]
479    
480 gross 1469 def stoppingcriterium2(self,norm_r,norm_b):
481     """
482     returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483    
484    
485     @param norm_r: current residual norm
486     @type norm_r: non-negative C{float}
487     @param norm_b: norm of right hand side
488     @type norm_b: non-negative C{float}
489     @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
490     @rtype: C{bool}
491    
492     """
493     self.history.append(norm_r)
494 artak 1489 if self.verbose: print "iter: #s: norm(r) = #e"#(len(self.history)-1, self.history[-1])
495 gross 1469 return self.history[-1]<=self.tolerance * norm_b
496    
497 gross 1330 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
498     """
499 ksteube 1312 Solver for
500    
501     M{Ax=b}
502    
503     with a symmetric and positive definite operator A (more details required!).
504     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
505    
506 gross 1330 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
507 ksteube 1312
508     For details on the preconditioned conjugate gradient method see the book:
509    
510     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
511     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
512     C. Romine, and H. van der Vorst.
513    
514     @param b: the right hand side of the liner system. C{b} is altered.
515 gross 1330 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
516 ksteube 1312 @param Aprod: returns the value Ax
517 gross 1330 @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.
518 ksteube 1312 @param Msolve: solves Mx=r
519 gross 1330 @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same
520     type like argument C{x}.
521 ksteube 1312 @param bilinearform: inner product C{<x,r>}
522 gross 1330 @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.
523     @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.
524     @type stoppingcriterium: function that returns C{True} or C{False}
525     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
526     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
527 ksteube 1312 @param iter_max: maximum number of iteration steps.
528     @type iter_max: C{int}
529 gross 1330 @return: the solution approximation and the corresponding residual
530     @rtype: C{tuple}
531     @warning: C{b} and C{x} are altered.
532 ksteube 1312 """
533     iter=0
534 gross 1330 if x==None:
535     x=0*b
536     else:
537     b += (-1)*Aprod(x)
538 ksteube 1312 r=b
539     rhat=Msolve(r)
540 gross 1330 d = rhat
541 ksteube 1312 rhat_dot_r = bilinearform(rhat, r)
542 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
543 ksteube 1312
544 gross 1330 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
545     iter+=1
546 ksteube 1312 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
547    
548     q=Aprod(d)
549     alpha = rhat_dot_r / bilinearform(d, q)
550     x += alpha * d
551     r += (-alpha) * q
552    
553     rhat=Msolve(r)
554     rhat_dot_r_new = bilinearform(rhat, r)
555     beta = rhat_dot_r_new / rhat_dot_r
556     rhat+=beta * d
557     d=rhat
558    
559     rhat_dot_r = rhat_dot_r_new
560 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
561 ksteube 1312
562 gross 1330 return x,r
563 ksteube 1312
564 artak 1465
565     ############################
566     # Added by Artak
567     #################################3
568    
569     #Apply a sequence of k Givens rotations, used within gmres codes
570     # vrot=givapp(c, s, vin, k)
571 gross 1467 def givapp(c,s,vin):
572     vrot=vin # warning: vin is altered!!!!
573     if isinstance(c,float):
574     vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
575     else:
576     for i in range(len(c)):
577     w1=c[i]*vrot[i]-s[i]*vrot[i+1]
578     w2=s[i]*vrot[i]+c[i]*vrot[i+1]
579     vrot[i:i+2]=w1,w2
580 artak 1465 return vrot
581    
582 artak 1475 def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
583     m=iter_restart
584 gross 1467 iter=0
585 artak 1475 while True:
586 artak 1488 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
587 artak 1475 x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588     iter+=iter_restart
589     if stopped: break
590     return x
591    
592     def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
593     iter=0
594 gross 1467 r=Msolve(b)
595     r_dot_r = bilinearform(r, r)
596     if r_dot_r<0: raise NegativeNorm,"negative norm."
597     norm_b=math.sqrt(r_dot_r)
598 artak 1465
599     if x==None:
600 artak 1488 x=0
601 gross 1467 else:
602     r=Msolve(b-Aprod(x))
603     r_dot_r = bilinearform(r, r)
604     if r_dot_r<0: raise NegativeNorm,"negative norm."
605 artak 1465
606 artak 1475 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
607     c=numarray.zeros(iter_restart,numarray.Float64)
608     s=numarray.zeros(iter_restart,numarray.Float64)
609     g=numarray.zeros(iter_restart,numarray.Float64)
610 artak 1465 v=[]
611    
612 gross 1467 rho=math.sqrt(r_dot_r)
613 artak 1488
614 gross 1467 v.append(r/rho)
615 artak 1465 g[0]=rho
616    
617 artak 1475 while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
618 artak 1465
619     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
620    
621    
622 gross 1467 p=Msolve(Aprod(v[iter]))
623 artak 1465
624     v.append(p)
625    
626     v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
627    
628     # Modified Gram-Schmidt
629     for j in range(iter+1):
630     h[j][iter]=bilinearform(v[j],v[iter+1])
631     v[iter+1]+=(-1.)*h[j][iter]*v[j]
632    
633     h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
634     v_norm2=h[iter+1][iter]
635    
636    
637     # Reorthogonalize if needed
638     if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
639     for j in range(iter+1):
640     hr=bilinearform(v[j],v[iter+1])
641     h[j][iter]=h[j][iter]+hr #vhat
642     v[iter+1] +=(-1.)*hr*v[j]
643    
644     v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
645     h[iter+1][iter]=v_norm2
646    
647     # watch out for happy breakdown
648     if v_norm2 != 0:
649     v[iter+1]=v[iter+1]/h[iter+1][iter]
650    
651     # Form and store the information for the new Givens rotation
652     if iter > 0 :
653     hhat=[]
654     for i in range(iter+1) : hhat.append(h[i][iter])
655 gross 1467 hhat=givapp(c[0:iter],s[0:iter],hhat);
656 artak 1465 for i in range(iter+1) : h[i][iter]=hhat[i]
657    
658     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
659     if mu!=0 :
660     c[iter]=h[iter][iter]/mu
661     s[iter]=-h[iter+1][iter]/mu
662     h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
663     h[iter+1][iter]=0.0
664 gross 1467 g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
665 artak 1465
666     # Update the residual norm
667     rho=abs(g[iter+1])
668     iter+=1
669    
670     # At this point either iter > iter_max or rho < tol.
671     # It's time to compute x and leave.
672    
673     if iter > 0 :
674 gross 1467 y=numarray.zeros(iter,numarray.Float64)
675 artak 1465 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
676     if iter > 1 :
677     i=iter-2
678     while i>=0 :
679 gross 1467 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
680 artak 1465 i=i-1
681     xhat=v[iter-1]*y[iter-1]
682     for i in range(iter-1):
683     xhat += v[i]*y[i]
684     else : xhat=v[0]
685    
686     x += xhat
687 artak 1488 if iter<iter_restart-1:
688 artak 1475 stopped=True
689     else:
690     stopped=False
691 artak 1465
692 artak 1475 return x,stopped
693 artak 1481
694     def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
695    
696     #
697     # minres solves the system of linear equations Ax = b
698     # where A is a symmetric matrix (possibly indefinite or singular)
699     # and b is a given vector.
700     #
701     # "A" may be a dense or sparse matrix (preferably sparse!)
702     # or the name of a function such that
703     # y = A(x)
704     # returns the product y = Ax for any given vector x.
705     #
706     # "M" defines a positive-definite preconditioner M = C C'.
707     # "M" may be a dense or sparse matrix (preferably sparse!)
708     # or the name of a function such that
709     # solves the system My = x for any given vector x.
710     #
711     #
712 artak 1482
713 artak 1481 #------------------------------------------------------------------
714     # Set up y and v for the first Lanczos vector v1.
715     # y = beta1 P' v1, where P = C**(-1).
716     # v is really P' v1.
717     #------------------------------------------------------------------
718 artak 1482 if x==None:
719     x=0*b
720     else:
721     b += (-1)*Aprod(x)
722    
723 artak 1481 r1 = b
724     y = Msolve(b)
725     beta1 = bilinearform(b,y)
726    
727     if beta1< 0: raise NegativeNorm,"negative norm."
728    
729     # If b = 0 exactly, stop with x = 0.
730     if beta1==0: return x*0.
731    
732     if beta1> 0:
733 artak 1486 beta1 = math.sqrt(beta1)
734 artak 1481
735     #------------------------------------------------------------------
736 artak 1484 # Initialize quantities.
737 artak 1481 # ------------------------------------------------------------------
738 artak 1482 iter = 0
739     Anorm = 0
740     ynorm = 0
741 artak 1481 oldb = 0
742     beta = beta1
743     dbar = 0
744     epsln = 0
745     phibar = beta1
746     rhs1 = beta1
747     rhs2 = 0
748     rnorm = phibar
749     tnorm2 = 0
750     ynorm2 = 0
751     cs = -1
752     sn = 0
753     w = b*0.
754     w2 = b*0.
755     r2 = r1
756     eps = 0.0001
757    
758     #---------------------------------------------------------------------
759     # Main iteration loop.
760     # --------------------------------------------------------------------
761 artak 1486 while not stoppingcriterium(rnorm,Anorm*ynorm): # checks ||r|| < (||A|| ||x||) * TOL
762 artak 1481
763     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
764     iter = iter + 1
765    
766     #-----------------------------------------------------------------
767     # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
768     # The general iteration is similar to the case k = 1 with v0 = 0:
769     #
770     # p1 = Operator * v1 - beta1 * v0,
771     # alpha1 = v1'p1,
772     # q2 = p2 - alpha1 * v1,
773     # beta2^2 = q2'q2,
774     # v2 = (1/beta2) q2.
775     #
776     # Again, y = betak P vk, where P = C**(-1).
777     #-----------------------------------------------------------------
778     s = 1/beta # Normalize previous vector (in y).
779     v = s*y # v = vk if P = I
780    
781     y = Aprod(v)
782 artak 1465
783 artak 1481 if iter >= 2:
784     y = y - (beta/oldb)*r1
785    
786     alfa = bilinearform(v,y) # alphak
787     y = (- alfa/beta)*r2 + y
788     r1 = r2
789     r2 = y
790     y = Msolve(r2)
791     oldb = beta # oldb = betak
792     beta = bilinearform(r2,y) # beta = betak+1^2
793     if beta < 0: raise NegativeNorm,"negative norm."
794    
795     beta = math.sqrt( beta )
796     tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
797    
798     if iter==1: # Initialize a few things.
799     gmax = abs( alfa ) # alpha1
800     gmin = gmax # alpha1
801    
802     # Apply previous rotation Qk-1 to get
803     # [deltak epslnk+1] = [cs sn][dbark 0 ]
804     # [gbar k dbar k+1] [sn -cs][alfak betak+1].
805    
806     oldeps = epsln
807     delta = cs * dbar + sn * alfa # delta1 = 0 deltak
808     gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
809     epsln = sn * beta # epsln2 = 0 epslnk+1
810     dbar = - cs * beta # dbar 2 = beta2 dbar k+1
811    
812     # Compute the next plane rotation Qk
813    
814     gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
815     gamma = max(gamma,eps)
816     cs = gbar / gamma # ck
817     sn = beta / gamma # sk
818     phi = cs * phibar # phik
819     phibar = sn * phibar # phibark+1
820    
821     # Update x.
822    
823     denom = 1/gamma
824     w1 = w2
825     w2 = w
826     w = (v - oldeps*w1 - delta*w2) * denom
827     x = x + phi*w
828    
829     # Go round again.
830    
831     gmax = max(gmax,gamma)
832     gmin = min(gmin,gamma)
833     z = rhs1 / gamma
834     ynorm2 = z*z + ynorm2
835     rhs1 = rhs2 - delta*z
836     rhs2 = - epsln*z
837    
838     # Estimate various norms and test for convergence.
839    
840     Anorm = math.sqrt( tnorm2 )
841     ynorm = math.sqrt( ynorm2 )
842    
843     rnorm = phibar
844    
845     return x
846    
847 artak 1489
848     def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
849    
850     # TFQMR solver for linear systems
851     #
852     #
853     # initialization
854     #
855     errtol = math.sqrt(bilinearform(b,b))
856     norm_b=errtol
857     kmax = iter_max
858     error = []
859    
860     if math.sqrt(bilinearform(x,x)) != 0.0:
861     r = b - Aprod(x)
862     else:
863     r = b
864    
865     r=Msolve(r)
866    
867     u1=0
868     u2=0
869     y1=0
870     y2=0
871    
872     w = r
873     y1 = r
874     iter = 0
875     d = 0
876    
877     v = Msolve(Aprod(y1))
878     u1 = v
879    
880     theta = 0.0;
881     eta = 0.0;
882     tau = math.sqrt(bilinearform(r,r))
883     error = [ error, tau ]
884     rho = tau * tau
885     m=1
886     #
887     # TFQMR iteration
888     #
889     # while ( iter < kmax-1 ):
890    
891     while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b):
892     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
893    
894     sigma = bilinearform(r,v)
895    
896     if ( sigma == 0.0 ):
897     raise 'TFQMR breakdown, sigma=0'
898    
899    
900     alpha = rho / sigma
901    
902     for j in range(2):
903     #
904     # Compute y2 and u2 only if you have to
905     #
906     if ( j == 1 ):
907     y2 = y1 - alpha * v
908     u2 = Msolve(Aprod(y2))
909    
910     m = 2 * (iter+1) - 2 + (j+1)
911     if j==0:
912     w = w - alpha * u1
913     d = y1 + ( theta * theta * eta / alpha ) * d
914     if j==1:
915     w = w - alpha * u2
916     d = y2 + ( theta * theta * eta / alpha ) * d
917    
918     theta = math.sqrt(bilinearform(w,w))/ tau
919     c = 1.0 / math.sqrt ( 1.0 + theta * theta )
920     tau = tau * theta * c
921     eta = c * c * alpha
922     x = x + eta * d
923     #
924     # Try to terminate the iteration at each pass through the loop
925     #
926     # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
927     # error = [ error, tau ]
928     # total_iters = iter
929     # break
930    
931    
932     if ( rho == 0.0 ):
933     raise 'TFQMR breakdown, rho=0'
934    
935    
936     rhon = bilinearform(r,w)
937     beta = rhon / rho;
938     rho = rhon;
939     y1 = w + beta * y2;
940     u1 = Msolve(Aprod(y1))
941     v = u1 + beta * ( u2 + beta * v )
942     error = [ error, tau ]
943     total_iters = iter
944    
945     iter = iter + 1
946    
947     print iter
948     return x
949    
950    
951 artak 1465 #############################################
952    
953 gross 1331 class ArithmeticTuple(object):
954     """
955     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
956    
957     example of usage:
958    
959     from esys.escript import Data
960     from numarray import array
961     a=Data(...)
962     b=array([1.,4.])
963     x=ArithmeticTuple(a,b)
964     y=5.*x
965    
966     """
967     def __init__(self,*args):
968     """
969     initialize object with elements args.
970    
971     @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
972     """
973     self.__items=list(args)
974    
975     def __len__(self):
976     """
977     number of items
978    
979     @return: number of items
980     @rtype: C{int}
981     """
982     return len(self.__items)
983    
984     def __getitem__(self,index):
985     """
986     get an item
987    
988     @param index: item to be returned
989     @type index: C{int}
990     @return: item with index C{index}
991     """
992     return self.__items.__getitem__(index)
993    
994     def __mul__(self,other):
995     """
996     scaling from the right
997    
998     @param other: scaling factor
999     @type other: C{float}
1000     @return: itemwise self*other
1001     @rtype: L{ArithmeticTuple}
1002     """
1003     out=[]
1004     for i in range(len(self)):
1005     out.append(self[i]*other)
1006     return ArithmeticTuple(*tuple(out))
1007    
1008     def __rmul__(self,other):
1009     """
1010     scaling from the left
1011    
1012     @param other: scaling factor
1013     @type other: C{float}
1014     @return: itemwise other*self
1015     @rtype: L{ArithmeticTuple}
1016     """
1017     out=[]
1018     for i in range(len(self)):
1019     out.append(other*self[i])
1020     return ArithmeticTuple(*tuple(out))
1021    
1022 artak 1465 #########################
1023     # Added by Artak
1024     #########################
1025     def __div__(self,other):
1026     """
1027     dividing from the right
1028    
1029     @param other: scaling factor
1030     @type other: C{float}
1031     @return: itemwise self/other
1032     @rtype: L{ArithmeticTuple}
1033     """
1034     out=[]
1035     for i in range(len(self)):
1036     out.append(self[i]/other)
1037     return ArithmeticTuple(*tuple(out))
1038    
1039     def __rdiv__(self,other):
1040     """
1041     dividing from the left
1042    
1043     @param other: scaling factor
1044     @type other: C{float}
1045     @return: itemwise other/self
1046     @rtype: L{ArithmeticTuple}
1047     """
1048     out=[]
1049     for i in range(len(self)):
1050     out.append(other/self[i])
1051     return ArithmeticTuple(*tuple(out))
1052    
1053     ##########################################33
1054    
1055 gross 1331 def __iadd__(self,other):
1056     """
1057     in-place add of other to self
1058    
1059     @param other: increment
1060     @type other: C{ArithmeticTuple}
1061     """
1062     if len(self) != len(other):
1063     raise ValueError,"tuple length must match."
1064     for i in range(len(self)):
1065     self.__items[i]+=other[i]
1066     return self
1067    
1068 gross 1414 class HomogeneousSaddlePointProblem(object):
1069     """
1070     This provides a framwork for solving homogeneous saddle point problem of the form
1071    
1072     Av+B^*p=f
1073     Bv =0
1074    
1075     for the unknowns v and p and given operators A and B and given right hand side f.
1076     B^* is the adjoint operator of B is the given inner product.
1077    
1078     """
1079     def __init__(self,**kwargs):
1080     self.setTolerance()
1081     self.setToleranceReductionFactor()
1082    
1083     def initialize(self):
1084     """
1085     initialize the problem (overwrite)
1086     """
1087     pass
1088     def B(self,v):
1089     """
1090     returns Bv (overwrite)
1091     @rtype: equal to the type of p
1092    
1093     @note: boundary conditions on p should be zero!
1094     """
1095     pass
1096    
1097     def inner(self,p0,p1):
1098     """
1099     returns inner product of two element p0 and p1 (overwrite)
1100    
1101     @type p0: equal to the type of p
1102     @type p1: equal to the type of p
1103     @rtype: C{float}
1104    
1105     @rtype: equal to the type of p
1106     """
1107     pass
1108    
1109     def solve_A(self,u,p):
1110     """
1111     solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1112    
1113     @rtype: equal to the type of v
1114     @note: boundary conditions on v should be zero!
1115     """
1116     pass
1117    
1118     def solve_prec(self,p):
1119     """
1120     provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1121    
1122     @rtype: equal to the type of p
1123     """
1124     pass
1125    
1126     def stoppingcriterium(self,Bv,v,p):
1127     """
1128     returns a True if iteration is terminated. (overwrite)
1129    
1130     @rtype: C{bool}
1131     """
1132     pass
1133    
1134     def __inner(self,p,r):
1135     return self.inner(p,r[1])
1136    
1137 artak 1465 def __inner_p(self,p1,p2):
1138     return self.inner(p1,p2)
1139    
1140 gross 1414 def __stoppingcriterium(self,norm_r,r,p):
1141     return self.stoppingcriterium(r[1],r[0],p)
1142    
1143 gross 1467 def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1144     return self.stoppingcriterium_GMRES(norm_r,norm_b)
1145 artak 1465
1146 artak 1481 def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1147     return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1148    
1149    
1150 gross 1414 def setTolerance(self,tolerance=1.e-8):
1151     self.__tol=tolerance
1152     def getTolerance(self):
1153     return self.__tol
1154     def setToleranceReductionFactor(self,reduction=0.01):
1155     self.__reduction=reduction
1156     def getSubProblemTolerance(self):
1157     return self.__reduction*self.getTolerance()
1158    
1159 artak 1488 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=10):
1160 gross 1414 """
1161     solves the saddle point problem using initial guesses v and p.
1162    
1163     @param max_iter: maximum number of iteration steps.
1164     """
1165     self.verbose=verbose
1166     self.show_details=show_details and self.verbose
1167    
1168 gross 1469 # assume p is known: then v=A^-1(f-B^*p)
1169     # which leads to BA^-1B^*p = BA^-1f
1170    
1171 gross 1414 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1172    
1173 artak 1465
1174 gross 1414 self.__z=v+self.solve_A(v,p*0)
1175 artak 1465
1176 gross 1414 Bz=self.B(self.__z)
1177     #
1178     # solve BA^-1B^*p = Bz
1179     #
1180     # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1181     #
1182     # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1183     # A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1184     #
1185     self.iter=0
1186 artak 1465 if solver=='GMRES':
1187     if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1188 artak 1488 p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1189 gross 1467 # solve Au=f-B^*p
1190     # A(u-v)=f-B^*p-Av
1191     # u=v+(u-v)
1192 artak 1465 u=v+self.solve_A(v,p)
1193 artak 1481
1194 artak 1489 if solver=='TFQMR':
1195     if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1196     p=TFQMR(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1197     # solve Au=f-B^*p
1198     # A(u-v)=f-B^*p-Av
1199     # u=v+(u-v)
1200     u=v+self.solve_A(v,p)
1201    
1202 artak 1481 if solver=='MINRES':
1203     if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1204 artak 1486 p=MINRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1205 artak 1481 # solve Au=f-B^*p
1206     # A(u-v)=f-B^*p-Av
1207     # u=v+(u-v)
1208     u=v+self.solve_A(v,p)
1209 artak 1488
1210     if solver=='PCG':
1211 artak 1465 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1212 gross 1467 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1213 artak 1465 u=r[0]
1214 artak 1481
1215 artak 1475 print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1216 gross 1414
1217 artak 1465 return u,p
1218    
1219 gross 1414 def __Msolve(self,r):
1220     return self.solve_prec(r[1])
1221    
1222 artak 1465 def __Msolve_GMRES(self,r):
1223     return self.solve_prec(r)
1224    
1225    
1226 gross 1414 def __Aprod(self,p):
1227     # return BA^-1B*p
1228     #solve Av =-B^*p as Av =f-Az-B^*p
1229 gross 1469 v=self.solve_A(self.__z,-p)
1230     return ArithmeticTuple(v, self.B(v))
1231 gross 1414
1232 artak 1465 def __Aprod_GMRES(self,p):
1233     # return BA^-1B*p
1234     #solve Av =-B^*p as Av =f-Az-B^*p
1235 gross 1469 v=self.solve_A(self.__z,-p)
1236     return self.B(v)
1237 gross 1414
1238 gross 867 class SaddlePointProblem(object):
1239     """
1240     This implements a solver for a saddlepoint problem
1241    
1242 gross 877 M{f(u,p)=0}
1243     M{g(u)=0}
1244 gross 867
1245     for u and p. The problem is solved with an inexact Uszawa scheme for p:
1246    
1247 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1248 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1249 gross 867
1250     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1251     A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
1252     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1253     in fact the role of a preconditioner.
1254     """
1255     def __init__(self,verbose=False,*args):
1256     """
1257     initializes the problem
1258    
1259 ksteube 990 @param verbose: switches on the printing out some information
1260 gross 867 @type verbose: C{bool}
1261     @note: this method may be overwritten by a particular saddle point problem
1262     """
1263 gross 1107 if not isinstance(verbose,bool):
1264     raise TypeError("verbose needs to be of type bool.")
1265 gross 1106 self.__verbose=verbose
1266 gross 877 self.relaxation=1.
1267 gross 867
1268     def trace(self,text):
1269     """
1270     prints text if verbose has been set
1271    
1272 ksteube 990 @param text: a text message
1273 gross 867 @type text: C{str}
1274     """
1275 artak 1489 if self.__verbose: print "#s: #s"%(str(self),text)
1276 gross 867
1277 gross 873 def solve_f(self,u,p,tol=1.e-8):
1278 gross 867 """
1279     solves
1280    
1281     A_f du = f(u,p)
1282    
1283     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1284    
1285     @param u: current approximation of u
1286     @type u: L{escript.Data}
1287     @param p: current approximation of p
1288     @type p: L{escript.Data}
1289 gross 873 @param tol: tolerance expected for du
1290 gross 867 @type tol: C{float}
1291     @return: increment du
1292     @rtype: L{escript.Data}
1293     @note: this method has to be overwritten by a particular saddle point problem
1294     """
1295     pass
1296    
1297 gross 873 def solve_g(self,u,tol=1.e-8):
1298 gross 867 """
1299     solves
1300    
1301     Q_g dp = g(u)
1302    
1303     with Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p.
1304    
1305     @param u: current approximation of u
1306     @type u: L{escript.Data}
1307 gross 873 @param tol: tolerance expected for dp
1308     @type tol: C{float}
1309 gross 867 @return: increment dp
1310     @rtype: L{escript.Data}
1311     @note: this method has to be overwritten by a particular saddle point problem
1312     """
1313     pass
1314    
1315     def inner(self,p0,p1):
1316     """
1317     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1318     @return: inner product of p0 and p1
1319     @rtype: C{float}
1320     """
1321     pass
1322    
1323 gross 877 subiter_max=3
1324     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1325     """
1326     runs the solver
1327 gross 873
1328 gross 877 @param u0: initial guess for C{u}
1329     @type u0: L{esys.escript.Data}
1330     @param p0: initial guess for C{p}
1331     @type p0: L{esys.escript.Data}
1332     @param tolerance: tolerance for relative error in C{u} and C{p}
1333     @type tolerance: positive C{float}
1334     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1335     @type tolerance_u: positive C{float}
1336     @param iter_max: maximum number of iteration steps.
1337     @type iter_max: C{int}
1338     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1339     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1340     @type accepted_reduction: positive C{float} or C{None}
1341     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1342     @type relaxation: C{float} or C{None}
1343     """
1344     tol=1.e-2
1345     if tolerance_u==None: tolerance_u=tolerance
1346     if not relaxation==None: self.relaxation=relaxation
1347     if accepted_reduction ==None:
1348     angle_limit=0.
1349     elif accepted_reduction>=1.:
1350     angle_limit=0.
1351     else:
1352     angle_limit=util.sqrt(1-accepted_reduction**2)
1353     self.iter=0
1354     u=u0
1355     p=p0
1356     #
1357     # initialize things:
1358     #
1359     converged=False
1360     #
1361     # start loop:
1362     #
1363     # initial search direction is g
1364     #
1365     while not converged :
1366     if self.iter>iter_max:
1367     raise ArithmeticError("no convergence after %s steps."%self.iter)
1368     f_new=self.solve_f(u,p,tol)
1369     norm_f_new = util.Lsup(f_new)
1370     u_new=u-f_new
1371     g_new=self.solve_g(u_new,tol)
1372     self.iter+=1
1373     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1374     if norm_f_new==0. and norm_g_new==0.: return u, p
1375     if self.iter>1 and not accepted_reduction==None:
1376     #
1377     # did we manage to reduce the norm of G? I
1378     # if not we start a backtracking procedure
1379     #
1380     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1381     if norm_g_new > accepted_reduction * norm_g:
1382     sub_iter=0
1383     s=self.relaxation
1384     d=g
1385     g_last=g
1386     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1387     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1388     dg= g_new-g_last
1389     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1390     rad=self.inner(g_new,dg)/self.relaxation
1391     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1392     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1393     if abs(rad) < norm_dg*norm_g_new * angle_limit:
1394     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1395     break
1396     r=self.relaxation
1397     self.relaxation= - rad/norm_dg**2
1398     s+=self.relaxation
1399     #####
1400     # a=g_new+self.relaxation*dg/r
1401     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1402     #####
1403     g_last=g_new
1404     p+=self.relaxation*d
1405     f_new=self.solve_f(u,p,tol)
1406     u_new=u-f_new
1407     g_new=self.solve_g(u_new,tol)
1408     self.iter+=1
1409     norm_f_new = util.Lsup(f_new)
1410     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1411     # print " ",sub_iter," new g norm",norm_g_new
1412     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1413     #
1414     # can we expect reduction of g?
1415     #
1416     # u_last=u_new
1417     sub_iter+=1
1418     self.relaxation=s
1419     #
1420     # check for convergence:
1421     #
1422     norm_u_new = util.Lsup(u_new)
1423     p_new=p+self.relaxation*g_new
1424     norm_p_new = util.sqrt(self.inner(p_new,p_new))
1425 ksteube 1125 self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
1426 gross 873
1427 gross 877 if self.iter>1:
1428     dg2=g_new-g
1429     df2=f_new-f
1430     norm_dg2=util.sqrt(self.inner(dg2,dg2))
1431     norm_df2=util.Lsup(df2)
1432     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1433     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1434     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1435     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1436     converged=True
1437     f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
1438     self.trace("convergence after %s steps."%self.iter)
1439     return u,p
1440     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1441     # tol=1.e-2
1442     # iter=0
1443     # converged=False
1444     # u=u0*1.
1445     # p=p0*1.
1446     # while not converged and iter<iter_max:
1447     # du=self.solve_f(u,p,tol)
1448     # u-=du
1449     # norm_du=util.Lsup(du)
1450     # norm_u=util.Lsup(u)
1451     #
1452     # dp=self.relaxation*self.solve_g(u,tol)
1453     # p+=dp
1454     # norm_dp=util.sqrt(self.inner(dp,dp))
1455     # norm_p=util.sqrt(self.inner(p,p))
1456     # print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
1457     # iter+=1
1458     #
1459     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1460     # if converged:
1461     # print "convergence after %s steps."%iter
1462     # else:
1463     # raise ArithmeticError("no convergence after %s steps."%iter)
1464     #
1465     # return u,p
1466 gross 873
1467 ksteube 1312 def MaskFromBoundaryTag(function_space,*tags):
1468     """
1469     create a mask on the given function space which one for samples
1470     that touch the boundary tagged by tags.
1471    
1472     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1473    
1474     @param function_space: a given function space
1475     @type function_space: L{escript.FunctionSpace}
1476     @param tags: boundray tags
1477     @type tags: C{str}
1478     @return: a mask which marks samples used by C{function_space} that are touching the
1479     boundary tagged by any of the given tags.
1480     @rtype: L{escript.Data} of rank 0
1481     """
1482     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1483     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1484     for t in tags: d.setTaggedValue(t,1.)
1485     pde.setValue(y=d)
1486     out=util.whereNonZero(pde.getRightHandSide())
1487     if out.getFunctionSpace() == function_space:
1488     return out
1489     else:
1490     return util.whereNonZero(util.interpolate(out,function_space))
1491    
1492 gross 1414
1493 artak 1465

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26