/[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 1465 - (hide annotations)
Wed Apr 2 03:28:25 2008 UTC (11 years, 3 months ago) by artak
File MIME type: text/x-python
File size: 38222 byte(s)


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     from Numeric import zeros,Int,Float32,Float64
53     ###################################
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     if self.verbose: print "iter: %s: inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])
478     return self.history[-1]<=self.tolerance * self.history[0]
479    
480     def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
481     """
482 ksteube 1312 Solver for
483    
484     M{Ax=b}
485    
486     with a symmetric and positive definite operator A (more details required!).
487     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
488    
489 gross 1330 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
490 ksteube 1312
491     For details on the preconditioned conjugate gradient method see the book:
492    
493     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
494     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
495     C. Romine, and H. van der Vorst.
496    
497     @param b: the right hand side of the liner system. C{b} is altered.
498 gross 1330 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
499 ksteube 1312 @param Aprod: returns the value Ax
500 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}.
501 ksteube 1312 @param Msolve: solves Mx=r
502 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
503     type like argument C{x}.
504 ksteube 1312 @param bilinearform: inner product C{<x,r>}
505 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}.
506     @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.
507     @type stoppingcriterium: function that returns C{True} or C{False}
508     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
509     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
510 ksteube 1312 @param iter_max: maximum number of iteration steps.
511     @type iter_max: C{int}
512 gross 1330 @return: the solution approximation and the corresponding residual
513     @rtype: C{tuple}
514     @warning: C{b} and C{x} are altered.
515 ksteube 1312 """
516     iter=0
517 gross 1330 if x==None:
518     x=0*b
519     else:
520     b += (-1)*Aprod(x)
521 ksteube 1312 r=b
522     rhat=Msolve(r)
523 gross 1330 d = rhat
524 ksteube 1312 rhat_dot_r = bilinearform(rhat, r)
525 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
526 ksteube 1312
527 gross 1330 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
528     iter+=1
529 ksteube 1312 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
530    
531     q=Aprod(d)
532     alpha = rhat_dot_r / bilinearform(d, q)
533     x += alpha * d
534     r += (-alpha) * q
535    
536     rhat=Msolve(r)
537     rhat_dot_r_new = bilinearform(rhat, r)
538     beta = rhat_dot_r_new / rhat_dot_r
539     rhat+=beta * d
540     d=rhat
541    
542     rhat_dot_r = rhat_dot_r_new
543 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
544 ksteube 1312
545 gross 1330 return x,r
546 ksteube 1312
547 artak 1465
548     ############################
549     # Added by Artak
550     #################################3
551    
552     #Apply a sequence of k Givens rotations, used within gmres codes
553     # vrot=givapp(c, s, vin, k)
554     def givapp(c,s,vin,k):
555     vrot=vin
556     for i in range(k+1):
557     w1=c[i]*vrot[i]-s[i]*vrot[i+1]
558     w2=s[i]*vrot[i]+c[i]*vrot[i+1]
559     vrot[i:i+2]=w1,w2
560     return vrot
561    
562     def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
563    
564     from numarray import dot
565    
566     v0=b
567     iter=0
568     if x==None:
569     x=0*b
570     else:
571     b += (-1.)*Aprod(x)
572     r=b
573    
574     rhat=Msolve(r)
575    
576     rhat_dot_r = bilinearform(rhat, rhat)
577     norm_r=math.sqrt(rhat_dot_r)
578    
579     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
580    
581     h=zeros((iter_max,iter_max),Float32)
582     c=zeros(iter_max,Float32)
583     s=zeros(iter_max,Float32)
584     g=zeros(iter_max,Float32)
585     v=[]
586    
587     v.append(rhat/norm_r)
588     rho=norm_r
589     g[0]=rho
590    
591     while not stoppingcriterium(rho,rho,norm_r):
592    
593     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
594    
595    
596     vhat=Aprod(v[iter])
597     p=Msolve(vhat)
598    
599     v.append(p)
600    
601     v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
602    
603     # Modified Gram-Schmidt
604     for j in range(iter+1):
605     h[j][iter]=bilinearform(v[j],v[iter+1])
606     v[iter+1]+=(-1.)*h[j][iter]*v[j]
607    
608     h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
609     v_norm2=h[iter+1][iter]
610    
611    
612     # Reorthogonalize if needed
613     if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
614     for j in range(iter+1):
615     hr=bilinearform(v[j],v[iter+1])
616     h[j][iter]=h[j][iter]+hr #vhat
617     v[iter+1] +=(-1.)*hr*v[j]
618    
619     v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
620     h[iter+1][iter]=v_norm2
621    
622     # watch out for happy breakdown
623     if v_norm2 != 0:
624     v[iter+1]=v[iter+1]/h[iter+1][iter]
625    
626     # Form and store the information for the new Givens rotation
627     if iter > 0 :
628     hhat=[]
629     for i in range(iter+1) : hhat.append(h[i][iter])
630     hhat=givapp(c[0:iter],s[0:iter],hhat,iter-1);
631     for i in range(iter+1) : h[i][iter]=hhat[i]
632    
633     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
634     if mu!=0 :
635     c[iter]=h[iter][iter]/mu
636     s[iter]=-h[iter+1][iter]/mu
637     h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
638     h[iter+1][iter]=0.0
639     g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2],0)
640    
641     # Update the residual norm
642     rho=abs(g[iter+1])
643     iter+=1
644    
645     # At this point either iter > iter_max or rho < tol.
646     # It's time to compute x and leave.
647    
648     if iter > 0 :
649     y=zeros(iter,Float32)
650     y[iter-1] = g[iter-1] / h[iter-1][iter-1]
651     if iter > 1 :
652     i=iter-2
653     while i>=0 :
654     y[i] = ( g[i] - dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
655     i=i-1
656     xhat=v[iter-1]*y[iter-1]
657     for i in range(iter-1):
658     xhat += v[i]*y[i]
659     else : xhat=v[0]
660    
661     x += xhat
662    
663     return x
664    
665     #############################################
666    
667 gross 1331 class ArithmeticTuple(object):
668     """
669     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
670    
671     example of usage:
672    
673     from esys.escript import Data
674     from numarray import array
675     a=Data(...)
676     b=array([1.,4.])
677     x=ArithmeticTuple(a,b)
678     y=5.*x
679    
680     """
681     def __init__(self,*args):
682     """
683     initialize object with elements args.
684    
685     @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
686     """
687     self.__items=list(args)
688    
689     def __len__(self):
690     """
691     number of items
692    
693     @return: number of items
694     @rtype: C{int}
695     """
696     return len(self.__items)
697    
698     def __getitem__(self,index):
699     """
700     get an item
701    
702     @param index: item to be returned
703     @type index: C{int}
704     @return: item with index C{index}
705     """
706     return self.__items.__getitem__(index)
707    
708     def __mul__(self,other):
709     """
710     scaling from the right
711    
712     @param other: scaling factor
713     @type other: C{float}
714     @return: itemwise self*other
715     @rtype: L{ArithmeticTuple}
716     """
717     out=[]
718     for i in range(len(self)):
719     out.append(self[i]*other)
720     return ArithmeticTuple(*tuple(out))
721    
722     def __rmul__(self,other):
723     """
724     scaling from the left
725    
726     @param other: scaling factor
727     @type other: C{float}
728     @return: itemwise other*self
729     @rtype: L{ArithmeticTuple}
730     """
731     out=[]
732     for i in range(len(self)):
733     out.append(other*self[i])
734     return ArithmeticTuple(*tuple(out))
735    
736 artak 1465 #########################
737     # Added by Artak
738     #########################
739     def __div__(self,other):
740     """
741     dividing from the right
742    
743     @param other: scaling factor
744     @type other: C{float}
745     @return: itemwise self/other
746     @rtype: L{ArithmeticTuple}
747     """
748     out=[]
749     for i in range(len(self)):
750     out.append(self[i]/other)
751     return ArithmeticTuple(*tuple(out))
752    
753     def __rdiv__(self,other):
754     """
755     dividing from the left
756    
757     @param other: scaling factor
758     @type other: C{float}
759     @return: itemwise other/self
760     @rtype: L{ArithmeticTuple}
761     """
762     out=[]
763     for i in range(len(self)):
764     out.append(other/self[i])
765     return ArithmeticTuple(*tuple(out))
766    
767     ##########################################33
768    
769 gross 1331 def __iadd__(self,other):
770     """
771     in-place add of other to self
772    
773     @param other: increment
774     @type other: C{ArithmeticTuple}
775     """
776     if len(self) != len(other):
777     raise ValueError,"tuple length must match."
778     for i in range(len(self)):
779     self.__items[i]+=other[i]
780     return self
781    
782 gross 1414 class HomogeneousSaddlePointProblem(object):
783     """
784     This provides a framwork for solving homogeneous saddle point problem of the form
785    
786     Av+B^*p=f
787     Bv =0
788    
789     for the unknowns v and p and given operators A and B and given right hand side f.
790     B^* is the adjoint operator of B is the given inner product.
791    
792     """
793     def __init__(self,**kwargs):
794     self.setTolerance()
795     self.setToleranceReductionFactor()
796    
797     def initialize(self):
798     """
799     initialize the problem (overwrite)
800     """
801     pass
802     def B(self,v):
803     """
804     returns Bv (overwrite)
805     @rtype: equal to the type of p
806    
807     @note: boundary conditions on p should be zero!
808     """
809     pass
810    
811     def inner(self,p0,p1):
812     """
813     returns inner product of two element p0 and p1 (overwrite)
814    
815     @type p0: equal to the type of p
816     @type p1: equal to the type of p
817     @rtype: C{float}
818    
819     @rtype: equal to the type of p
820     """
821     pass
822    
823     def solve_A(self,u,p):
824     """
825     solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
826    
827     @rtype: equal to the type of v
828     @note: boundary conditions on v should be zero!
829     """
830     pass
831    
832     def solve_prec(self,p):
833     """
834     provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
835    
836     @rtype: equal to the type of p
837     """
838     pass
839    
840     def stoppingcriterium(self,Bv,v,p):
841     """
842     returns a True if iteration is terminated. (overwrite)
843    
844     @rtype: C{bool}
845     """
846     pass
847    
848     def __inner(self,p,r):
849     return self.inner(p,r[1])
850    
851 artak 1465 def __inner_p(self,p1,p2):
852     return self.inner(p1,p2)
853    
854 gross 1414 def __stoppingcriterium(self,norm_r,r,p):
855     return self.stoppingcriterium(r[1],r[0],p)
856    
857 artak 1465 def __stoppingcriterium_GMRES(self,norm_r,rho,r):
858     return self.stoppingcriterium_GMRES(rho,r)
859    
860 gross 1414 def setTolerance(self,tolerance=1.e-8):
861     self.__tol=tolerance
862     def getTolerance(self):
863     return self.__tol
864     def setToleranceReductionFactor(self,reduction=0.01):
865     self.__reduction=reduction
866     def getSubProblemTolerance(self):
867     return self.__reduction*self.getTolerance()
868    
869 artak 1465 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
870 gross 1414 """
871     solves the saddle point problem using initial guesses v and p.
872    
873     @param max_iter: maximum number of iteration steps.
874     """
875     self.verbose=verbose
876     self.show_details=show_details and self.verbose
877    
878     # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
879    
880 artak 1465
881 gross 1414 self.__z=v+self.solve_A(v,p*0)
882 artak 1465
883 gross 1414 Bz=self.B(self.__z)
884     #
885     # solve BA^-1B^*p = Bz
886     #
887     # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
888     #
889     # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
890     # A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
891     #
892     self.iter=0
893 artak 1465 if solver=='GMRES':
894     if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
895     p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1)
896     u=v+self.solve_A(v,p)
897    
898     else:
899     if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
900     p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p*1)
901     u=r[0]
902     u=v+self.solve_A(v,p) # Lutz to check !!!!!
903 gross 1414
904 artak 1465 return u,p
905    
906 gross 1414 def __Msolve(self,r):
907     return self.solve_prec(r[1])
908    
909 artak 1465 def __Msolve_GMRES(self,r):
910     return self.solve_prec(r)
911    
912    
913 gross 1414 def __Aprod(self,p):
914     # return BA^-1B*p
915     #solve Av =-B^*p as Av =f-Az-B^*p
916     v=self.solve_A(self.__z,p)
917     return ArithmeticTuple(v, self.B(v))
918    
919 artak 1465 def __Aprod_GMRES(self,p):
920     # return BA^-1B*p
921     #solve Av =-B^*p as Av =f-Az-B^*p
922     v=self.solve_A(self.__z,p)
923     return self.B(v)
924 gross 1414
925 gross 867 class SaddlePointProblem(object):
926     """
927     This implements a solver for a saddlepoint problem
928    
929 gross 877 M{f(u,p)=0}
930     M{g(u)=0}
931 gross 867
932     for u and p. The problem is solved with an inexact Uszawa scheme for p:
933    
934 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
935 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
936 gross 867
937     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
938     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'
939     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
940     in fact the role of a preconditioner.
941     """
942     def __init__(self,verbose=False,*args):
943     """
944     initializes the problem
945    
946 ksteube 990 @param verbose: switches on the printing out some information
947 gross 867 @type verbose: C{bool}
948     @note: this method may be overwritten by a particular saddle point problem
949     """
950 gross 1107 if not isinstance(verbose,bool):
951     raise TypeError("verbose needs to be of type bool.")
952 gross 1106 self.__verbose=verbose
953 gross 877 self.relaxation=1.
954 gross 867
955     def trace(self,text):
956     """
957     prints text if verbose has been set
958    
959 ksteube 990 @param text: a text message
960 gross 867 @type text: C{str}
961     """
962     if self.__verbose: print "%s: %s"%(str(self),text)
963    
964 gross 873 def solve_f(self,u,p,tol=1.e-8):
965 gross 867 """
966     solves
967    
968     A_f du = f(u,p)
969    
970     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
971    
972     @param u: current approximation of u
973     @type u: L{escript.Data}
974     @param p: current approximation of p
975     @type p: L{escript.Data}
976 gross 873 @param tol: tolerance expected for du
977 gross 867 @type tol: C{float}
978     @return: increment du
979     @rtype: L{escript.Data}
980     @note: this method has to be overwritten by a particular saddle point problem
981     """
982     pass
983    
984 gross 873 def solve_g(self,u,tol=1.e-8):
985 gross 867 """
986     solves
987    
988     Q_g dp = g(u)
989    
990     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.
991    
992     @param u: current approximation of u
993     @type u: L{escript.Data}
994 gross 873 @param tol: tolerance expected for dp
995     @type tol: C{float}
996 gross 867 @return: increment dp
997     @rtype: L{escript.Data}
998     @note: this method has to be overwritten by a particular saddle point problem
999     """
1000     pass
1001    
1002     def inner(self,p0,p1):
1003     """
1004     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1005     @return: inner product of p0 and p1
1006     @rtype: C{float}
1007     """
1008     pass
1009    
1010 gross 877 subiter_max=3
1011     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1012     """
1013     runs the solver
1014 gross 873
1015 gross 877 @param u0: initial guess for C{u}
1016     @type u0: L{esys.escript.Data}
1017     @param p0: initial guess for C{p}
1018     @type p0: L{esys.escript.Data}
1019     @param tolerance: tolerance for relative error in C{u} and C{p}
1020     @type tolerance: positive C{float}
1021     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1022     @type tolerance_u: positive C{float}
1023     @param iter_max: maximum number of iteration steps.
1024     @type iter_max: C{int}
1025     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1026     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1027     @type accepted_reduction: positive C{float} or C{None}
1028     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1029     @type relaxation: C{float} or C{None}
1030     """
1031     tol=1.e-2
1032     if tolerance_u==None: tolerance_u=tolerance
1033     if not relaxation==None: self.relaxation=relaxation
1034     if accepted_reduction ==None:
1035     angle_limit=0.
1036     elif accepted_reduction>=1.:
1037     angle_limit=0.
1038     else:
1039     angle_limit=util.sqrt(1-accepted_reduction**2)
1040     self.iter=0
1041     u=u0
1042     p=p0
1043     #
1044     # initialize things:
1045     #
1046     converged=False
1047     #
1048     # start loop:
1049     #
1050     # initial search direction is g
1051     #
1052     while not converged :
1053     if self.iter>iter_max:
1054     raise ArithmeticError("no convergence after %s steps."%self.iter)
1055     f_new=self.solve_f(u,p,tol)
1056     norm_f_new = util.Lsup(f_new)
1057     u_new=u-f_new
1058     g_new=self.solve_g(u_new,tol)
1059     self.iter+=1
1060     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1061     if norm_f_new==0. and norm_g_new==0.: return u, p
1062     if self.iter>1 and not accepted_reduction==None:
1063     #
1064     # did we manage to reduce the norm of G? I
1065     # if not we start a backtracking procedure
1066     #
1067     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1068     if norm_g_new > accepted_reduction * norm_g:
1069     sub_iter=0
1070     s=self.relaxation
1071     d=g
1072     g_last=g
1073     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1074     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1075     dg= g_new-g_last
1076     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1077     rad=self.inner(g_new,dg)/self.relaxation
1078     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1079     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1080     if abs(rad) < norm_dg*norm_g_new * angle_limit:
1081     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1082     break
1083     r=self.relaxation
1084     self.relaxation= - rad/norm_dg**2
1085     s+=self.relaxation
1086     #####
1087     # a=g_new+self.relaxation*dg/r
1088     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1089     #####
1090     g_last=g_new
1091     p+=self.relaxation*d
1092     f_new=self.solve_f(u,p,tol)
1093     u_new=u-f_new
1094     g_new=self.solve_g(u_new,tol)
1095     self.iter+=1
1096     norm_f_new = util.Lsup(f_new)
1097     norm_g_new = util.sqrt(self.inner(g_new,g_new))
1098     # print " ",sub_iter," new g norm",norm_g_new
1099     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1100     #
1101     # can we expect reduction of g?
1102     #
1103     # u_last=u_new
1104     sub_iter+=1
1105     self.relaxation=s
1106     #
1107     # check for convergence:
1108     #
1109     norm_u_new = util.Lsup(u_new)
1110     p_new=p+self.relaxation*g_new
1111     norm_p_new = util.sqrt(self.inner(p_new,p_new))
1112 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))
1113 gross 873
1114 gross 877 if self.iter>1:
1115     dg2=g_new-g
1116     df2=f_new-f
1117     norm_dg2=util.sqrt(self.inner(dg2,dg2))
1118     norm_df2=util.Lsup(df2)
1119     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1120     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1121     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1122     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1123     converged=True
1124     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
1125     self.trace("convergence after %s steps."%self.iter)
1126     return u,p
1127     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1128     # tol=1.e-2
1129     # iter=0
1130     # converged=False
1131     # u=u0*1.
1132     # p=p0*1.
1133     # while not converged and iter<iter_max:
1134     # du=self.solve_f(u,p,tol)
1135     # u-=du
1136     # norm_du=util.Lsup(du)
1137     # norm_u=util.Lsup(u)
1138     #
1139     # dp=self.relaxation*self.solve_g(u,tol)
1140     # p+=dp
1141     # norm_dp=util.sqrt(self.inner(dp,dp))
1142     # norm_p=util.sqrt(self.inner(p,p))
1143     # 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)
1144     # iter+=1
1145     #
1146     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1147     # if converged:
1148     # print "convergence after %s steps."%iter
1149     # else:
1150     # raise ArithmeticError("no convergence after %s steps."%iter)
1151     #
1152     # return u,p
1153 gross 873
1154 ksteube 1312 def MaskFromBoundaryTag(function_space,*tags):
1155     """
1156     create a mask on the given function space which one for samples
1157     that touch the boundary tagged by tags.
1158    
1159     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1160    
1161     @param function_space: a given function space
1162     @type function_space: L{escript.FunctionSpace}
1163     @param tags: boundray tags
1164     @type tags: C{str}
1165     @return: a mask which marks samples used by C{function_space} that are touching the
1166     boundary tagged by any of the given tags.
1167     @rtype: L{escript.Data} of rank 0
1168     """
1169     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1170     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1171     for t in tags: d.setTaggedValue(t,1.)
1172     pde.setValue(y=d)
1173     out=util.whereNonZero(pde.getRightHandSide())
1174     if out.getFunctionSpace() == function_space:
1175     return out
1176     else:
1177     return util.whereNonZero(util.interpolate(out,function_space))
1178    
1179 gross 1414
1180 artak 1465

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26