/[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 1107 - (hide annotations)
Thu Apr 19 02:14:18 2007 UTC (12 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 21552 byte(s)
Small bug in ReadMeh call fixed.
Test for typ of verbose argument in pdetools.SaddlePointSolver added.


1 jgs 121 # $Id$
2    
3     """
4 jgs 149 Provides some tools related to PDEs.
5 jgs 121
6 jgs 149 Currently includes:
7     - Projector - to project a discontinuous
8 gross 351 - Locator - to trace values in data objects at a certain location
9     - TimeIntegrationManager - to handel extraplotion in time
10 gross 867 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
11 gross 637
12     @var __author__: name of author
13     @var __copyright__: copyrights
14     @var __license__: licence agreement
15     @var __url__: url entry point on documentation
16     @var __version__: version
17     @var __date__: date of the version
18 jgs 121 """
19    
20 gross 637 __author__="Lutz Gross, l.gross@uq.edu.au"
21 elspeth 609 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
22     http://www.access.edu.au
23     Primary Business: Queensland, Australia"""
24 elspeth 614 __license__="""Licensed under the Open Software License version 3.0
25     http://www.opensource.org/licenses/osl-3.0.php"""
26 gross 637 __url__="http://www.iservo.edu.au/esys"
27     __version__="$Revision$"
28     __date__="$Date$"
29 elspeth 609
30 gross 637
31 jgs 149 import escript
32     import linearPDEs
33 jgs 121 import numarray
34 jgs 149 import util
35 jgs 121
36 gross 351 class TimeIntegrationManager:
37     """
38     a simple mechanism to manage time dependend values.
39    
40 gross 720 typical usage is::
41 gross 351
42 gross 720 dt=0.1 # time increment
43     tm=TimeIntegrationManager(inital_value,p=1)
44     while t<1.
45     v_guess=tm.extrapolate(dt) # extrapolate to t+dt
46     v=...
47     tm.checkin(dt,v)
48     t+=dt
49 gross 351
50 gross 720 @note: currently only p=1 is supported.
51 gross 351 """
52     def __init__(self,*inital_values,**kwargs):
53     """
54     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
55     """
56     if kwargs.has_key("p"):
57     self.__p=kwargs["p"]
58     else:
59     self.__p=1
60     if kwargs.has_key("time"):
61     self.__t=kwargs["time"]
62     else:
63     self.__t=0.
64     self.__v_mem=[inital_values]
65     self.__order=0
66     self.__dt_mem=[]
67     self.__num_val=len(inital_values)
68    
69     def getTime(self):
70     return self.__t
71 gross 396 def getValue(self):
72 gross 409 out=self.__v_mem[0]
73     if len(out)==1:
74     return out[0]
75     else:
76     return out
77    
78 gross 351 def checkin(self,dt,*values):
79     """
80     adds new values to the manager. the p+1 last value get lost
81     """
82     o=min(self.__order+1,self.__p)
83     self.__order=min(self.__order+1,self.__p)
84     v_mem_new=[values]
85     dt_mem_new=[dt]
86     for i in range(o-1):
87     v_mem_new.append(self.__v_mem[i])
88     dt_mem_new.append(self.__dt_mem[i])
89     v_mem_new.append(self.__v_mem[o-1])
90     self.__order=o
91     self.__v_mem=v_mem_new
92     self.__dt_mem=dt_mem_new
93     self.__t+=dt
94    
95     def extrapolate(self,dt):
96     """
97     extrapolates to dt forward in time.
98     """
99     if self.__order==0:
100     out=self.__v_mem[0]
101     else:
102     out=[]
103     for i in range(self.__num_val):
104     out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
105    
106     if len(out)==0:
107     return None
108     elif len(out)==1:
109     return out[0]
110     else:
111     return out
112 gross 396
113 gross 351
114 jgs 121 class Projector:
115 jgs 149 """
116     The Projector is a factory which projects a discontiuous function onto a
117     continuous function on the a given domain.
118     """
119 jgs 121 def __init__(self, domain, reduce = True, fast=True):
120     """
121 jgs 149 Create a continuous function space projector for a domain.
122 jgs 121
123 jgs 149 @param domain: Domain of the projection.
124     @param reduce: Flag to reduce projection order (default is True)
125     @param fast: Flag to use a fast method based on matrix lumping (default is true)
126 jgs 121 """
127 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
128 jgs 148 if fast:
129 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
130 jgs 121 self.__pde.setSymmetryOn()
131     self.__pde.setReducedOrderTo(reduce)
132     self.__pde.setValue(D = 1.)
133     return
134    
135     def __call__(self, input_data):
136     """
137 jgs 149 Projects input_data onto a continuous function
138 jgs 121
139 jgs 149 @param input_data: The input_data to be projected.
140 jgs 121 """
141 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
142 jgs 121 if input_data.getRank()==0:
143     self.__pde.setValue(Y = input_data)
144     out=self.__pde.getSolution()
145     elif input_data.getRank()==1:
146     for i0 in range(input_data.getShape()[0]):
147     self.__pde.setValue(Y = input_data[i0])
148     out[i0]=self.__pde.getSolution()
149     elif input_data.getRank()==2:
150     for i0 in range(input_data.getShape()[0]):
151     for i1 in range(input_data.getShape()[1]):
152     self.__pde.setValue(Y = input_data[i0,i1])
153     out[i0,i1]=self.__pde.getSolution()
154     elif input_data.getRank()==3:
155     for i0 in range(input_data.getShape()[0]):
156     for i1 in range(input_data.getShape()[1]):
157     for i2 in range(input_data.getShape()[2]):
158     self.__pde.setValue(Y = input_data[i0,i1,i2])
159     out[i0,i1,i2]=self.__pde.getSolution()
160     else:
161     for i0 in range(input_data.getShape()[0]):
162     for i1 in range(input_data.getShape()[1]):
163     for i2 in range(input_data.getShape()[2]):
164     for i3 in range(input_data.getShape()[3]):
165     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
166     out[i0,i1,i2,i3]=self.__pde.getSolution()
167     return out
168    
169 gross 525 class NoPDE:
170     """
171     solves the following problem for u:
172 jgs 121
173 gross 525 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
174    
175     with constraint
176    
177     M{u[j]=r[j]} where M{q[j]>0}
178    
179     where D, Y, r and q are given functions of rank 1.
180    
181     In the case of scalars this takes the form
182    
183     M{D*u=Y}
184    
185     with constraint
186    
187     M{u=r} where M{q>0}
188    
189     where D, Y, r and q are given scalar functions.
190    
191     The constraint is overwriting any other condition.
192    
193 gross 720 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
194     that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
195     thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
196 gross 525 """
197     def __init__(self,domain,D=None,Y=None,q=None,r=None):
198     """
199     initialize the problem
200    
201     @param domain: domain of the PDE.
202     @type domain: L{Domain}
203     @param D: coefficient of the solution.
204 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
205 gross 525 @param Y: right hand side
206 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207 gross 525 @param q: location of constraints
208 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209 gross 525 @param r: value of solution at locations of constraints
210 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211 gross 525 """
212     self.__domain=domain
213     self.__D=D
214     self.__Y=Y
215     self.__q=q
216     self.__r=r
217     self.__u=None
218     self.__function_space=escript.Solution(self.__domain)
219     def setReducedOn(self):
220     """
221     sets the L{FunctionSpace} of the solution to L{ReducedSolution}
222     """
223     self.__function_space=escript.ReducedSolution(self.__domain)
224     self.__u=None
225    
226     def setReducedOff(self):
227     """
228     sets the L{FunctionSpace} of the solution to L{Solution}
229     """
230     self.__function_space=escript.Solution(self.__domain)
231     self.__u=None
232    
233     def setValue(self,D=None,Y=None,q=None,r=None):
234     """
235     assigns values to the parameters.
236    
237     @param D: coefficient of the solution.
238 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
239 gross 525 @param Y: right hand side
240 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241 gross 525 @param q: location of constraints
242 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243 gross 525 @param r: value of solution at locations of constraints
244 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245 gross 525 """
246     if not D==None:
247     self.__D=D
248     self.__u=None
249     if not Y==None:
250     self.__Y=Y
251     self.__u=None
252     if not q==None:
253     self.__q=q
254     self.__u=None
255     if not r==None:
256     self.__r=r
257     self.__u=None
258    
259     def getSolution(self):
260     """
261     returns the solution
262    
263     @return: the solution of the problem
264     @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
265     """
266     if self.__u==None:
267     if self.__D==None:
268     raise ValueError,"coefficient D is undefined"
269     D=escript.Data(self.__D,self.__function_space)
270     if D.getRank()>1:
271     raise ValueError,"coefficient D must have rank 0 or 1"
272     if self.__Y==None:
273     self.__u=escript.Data(0.,D.getShape(),self.__function_space)
274     else:
275     self.__u=util.quotient(self.__Y,D)
276     if not self.__q==None:
277     q=util.wherePositive(escript.Data(self.__q,self.__function_space))
278     self.__u*=(1.-q)
279     if not self.__r==None: self.__u+=q*self.__r
280     return self.__u
281    
282 jgs 147 class Locator:
283     """
284 jgs 149 Locator provides access to the values of data objects at a given
285     spatial coordinate x.
286    
287     In fact, a Locator object finds the sample in the set of samples of a
288     given function space or domain where which is closest to the given
289     point x.
290 jgs 147 """
291    
292     def __init__(self,where,x=numarray.zeros((3,))):
293 jgs 149 """
294     Initializes a Locator to access values in Data objects on the Doamin
295     or FunctionSpace where for the sample point which
296     closest to the given point x.
297 gross 880
298     @param where: function space
299     @type where: L{escript.FunctionSpace}
300     @param x: coefficient of the solution.
301     @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
302 jgs 149 """
303     if isinstance(where,escript.FunctionSpace):
304 jgs 147 self.__function_space=where
305 jgs 121 else:
306 jgs 149 self.__function_space=escript.ContinuousFunction(where)
307 gross 880 if isinstance(x, list):
308     self.__id=[]
309     for p in x:
310 gross 921 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
311 gross 880 else:
312 gross 921 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
313 jgs 121
314 jgs 147 def __str__(self):
315 jgs 149 """
316     Returns the coordinates of the Locator as a string.
317     """
318 gross 880 x=self.getX()
319     if instance(x,list):
320     out="["
321     first=True
322     for xx in x:
323     if not first:
324     out+=","
325     else:
326     first=False
327     out+=str(xx)
328     out+="]>"
329     else:
330     out=str(x)
331     return out
332 jgs 121
333 gross 880 def getX(self):
334     """
335     Returns the exact coordinates of the Locator.
336     """
337     return self(self.getFunctionSpace().getX())
338    
339 jgs 147 def getFunctionSpace(self):
340 jgs 149 """
341     Returns the function space of the Locator.
342     """
343 jgs 147 return self.__function_space
344    
345 gross 880 def getId(self,item=None):
346 jgs 149 """
347     Returns the identifier of the location.
348     """
349 gross 880 if item == None:
350     return self.__id
351     else:
352     if isinstance(self.__id,list):
353     return self.__id[item]
354     else:
355     return self.__id
356 jgs 121
357    
358 jgs 147 def __call__(self,data):
359 jgs 149 """
360     Returns the value of data at the Locator of a Data object otherwise
361     the object is returned.
362     """
363 jgs 147 return self.getValue(data)
364 jgs 121
365 jgs 147 def getValue(self,data):
366 jgs 149 """
367     Returns the value of data at the Locator if data is a Data object
368     otherwise the object is returned.
369     """
370     if isinstance(data,escript.Data):
371 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
372 gross 880 dat=data
373 jgs 147 else:
374 gross 880 dat=data.interpolate(self.getFunctionSpace())
375     id=self.getId()
376     r=data.getRank()
377     if isinstance(id,list):
378     out=[]
379     for i in id:
380 gross 921 o=data.getValueOfGlobalDataPoint(*i)
381 gross 880 if data.getRank()==0:
382     out.append(o[0])
383     else:
384     out.append(o)
385     return out
386 jgs 147 else:
387 gross 921 out=data.getValueOfGlobalDataPoint(*id)
388 gross 880 if data.getRank()==0:
389     return out[0]
390     else:
391     return out
392 jgs 147 else:
393     return data
394 jgs 149
395 gross 867 class SaddlePointProblem(object):
396     """
397     This implements a solver for a saddlepoint problem
398    
399 gross 877 M{f(u,p)=0}
400     M{g(u)=0}
401 gross 867
402     for u and p. The problem is solved with an inexact Uszawa scheme for p:
403    
404 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
405 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
406 gross 867
407     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
408     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'
409     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
410     in fact the role of a preconditioner.
411     """
412     def __init__(self,verbose=False,*args):
413     """
414     initializes the problem
415    
416 ksteube 990 @param verbose: switches on the printing out some information
417 gross 867 @type verbose: C{bool}
418     @note: this method may be overwritten by a particular saddle point problem
419     """
420 gross 1107 if not isinstance(verbose,bool):
421     raise TypeError("verbose needs to be of type bool.")
422 gross 1106 self.__verbose=verbose
423 gross 877 self.relaxation=1.
424 gross 867
425     def trace(self,text):
426     """
427     prints text if verbose has been set
428    
429 ksteube 990 @param text: a text message
430 gross 867 @type text: C{str}
431     """
432     if self.__verbose: print "%s: %s"%(str(self),text)
433    
434 gross 873 def solve_f(self,u,p,tol=1.e-8):
435 gross 867 """
436     solves
437    
438     A_f du = f(u,p)
439    
440     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
441    
442     @param u: current approximation of u
443     @type u: L{escript.Data}
444     @param p: current approximation of p
445     @type p: L{escript.Data}
446 gross 873 @param tol: tolerance expected for du
447 gross 867 @type tol: C{float}
448     @return: increment du
449     @rtype: L{escript.Data}
450     @note: this method has to be overwritten by a particular saddle point problem
451     """
452     pass
453    
454 gross 873 def solve_g(self,u,tol=1.e-8):
455 gross 867 """
456     solves
457    
458     Q_g dp = g(u)
459    
460     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.
461    
462     @param u: current approximation of u
463     @type u: L{escript.Data}
464 gross 873 @param tol: tolerance expected for dp
465     @type tol: C{float}
466 gross 867 @return: increment dp
467     @rtype: L{escript.Data}
468     @note: this method has to be overwritten by a particular saddle point problem
469     """
470     pass
471    
472     def inner(self,p0,p1):
473     """
474     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
475     @return: inner product of p0 and p1
476     @rtype: C{float}
477     """
478     pass
479    
480 gross 877 subiter_max=3
481     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
482     """
483     runs the solver
484 gross 873
485 gross 877 @param u0: initial guess for C{u}
486     @type u0: L{esys.escript.Data}
487     @param p0: initial guess for C{p}
488     @type p0: L{esys.escript.Data}
489     @param tolerance: tolerance for relative error in C{u} and C{p}
490     @type tolerance: positive C{float}
491     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
492     @type tolerance_u: positive C{float}
493     @param iter_max: maximum number of iteration steps.
494     @type iter_max: C{int}
495     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
496     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
497     @type accepted_reduction: positive C{float} or C{None}
498     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
499     @type relaxation: C{float} or C{None}
500     """
501     tol=1.e-2
502     if tolerance_u==None: tolerance_u=tolerance
503     if not relaxation==None: self.relaxation=relaxation
504     if accepted_reduction ==None:
505     angle_limit=0.
506     elif accepted_reduction>=1.:
507     angle_limit=0.
508     else:
509     angle_limit=util.sqrt(1-accepted_reduction**2)
510     self.iter=0
511     u=u0
512     p=p0
513     #
514     # initialize things:
515     #
516     converged=False
517     #
518     # start loop:
519     #
520     # initial search direction is g
521     #
522     while not converged :
523     if self.iter>iter_max:
524     raise ArithmeticError("no convergence after %s steps."%self.iter)
525     f_new=self.solve_f(u,p,tol)
526     norm_f_new = util.Lsup(f_new)
527     u_new=u-f_new
528     g_new=self.solve_g(u_new,tol)
529     self.iter+=1
530     norm_g_new = util.sqrt(self.inner(g_new,g_new))
531     if norm_f_new==0. and norm_g_new==0.: return u, p
532     if self.iter>1 and not accepted_reduction==None:
533     #
534     # did we manage to reduce the norm of G? I
535     # if not we start a backtracking procedure
536     #
537     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
538     if norm_g_new > accepted_reduction * norm_g:
539     sub_iter=0
540     s=self.relaxation
541     d=g
542     g_last=g
543     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
544     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
545     dg= g_new-g_last
546     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
547     rad=self.inner(g_new,dg)/self.relaxation
548     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
549     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
550     if abs(rad) < norm_dg*norm_g_new * angle_limit:
551     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
552     break
553     r=self.relaxation
554     self.relaxation= - rad/norm_dg**2
555     s+=self.relaxation
556     #####
557     # a=g_new+self.relaxation*dg/r
558     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
559     #####
560     g_last=g_new
561     p+=self.relaxation*d
562     f_new=self.solve_f(u,p,tol)
563     u_new=u-f_new
564     g_new=self.solve_g(u_new,tol)
565     self.iter+=1
566     norm_f_new = util.Lsup(f_new)
567     norm_g_new = util.sqrt(self.inner(g_new,g_new))
568     # print " ",sub_iter," new g norm",norm_g_new
569     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
570     #
571     # can we expect reduction of g?
572     #
573     # u_last=u_new
574     sub_iter+=1
575     self.relaxation=s
576     #
577     # check for convergence:
578     #
579     norm_u_new = util.Lsup(u_new)
580     p_new=p+self.relaxation*g_new
581     norm_p_new = util.sqrt(self.inner(p_new,p_new))
582     self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))
583 gross 873
584 gross 877 if self.iter>1:
585     dg2=g_new-g
586     df2=f_new-f
587     norm_dg2=util.sqrt(self.inner(dg2,dg2))
588     norm_df2=util.Lsup(df2)
589     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
590     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
591     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
592     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
593     converged=True
594     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
595     self.trace("convergence after %s steps."%self.iter)
596     return u,p
597     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
598     # tol=1.e-2
599     # iter=0
600     # converged=False
601     # u=u0*1.
602     # p=p0*1.
603     # while not converged and iter<iter_max:
604     # du=self.solve_f(u,p,tol)
605     # u-=du
606     # norm_du=util.Lsup(du)
607     # norm_u=util.Lsup(u)
608     #
609     # dp=self.relaxation*self.solve_g(u,tol)
610     # p+=dp
611     # norm_dp=util.sqrt(self.inner(dp,dp))
612     # norm_p=util.sqrt(self.inner(p,p))
613     # 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)
614     # iter+=1
615     #
616     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
617     # if converged:
618     # print "convergence after %s steps."%iter
619     # else:
620     # raise ArithmeticError("no convergence after %s steps."%iter)
621     #
622     # return u,p
623 gross 873
624 jgs 149 # vim: expandtab shiftwidth=4:

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26