/[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 877 - (hide annotations)
Wed Oct 25 03:06:58 2006 UTC (12 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 20538 byte(s)
backtraking in the saddlepoint problem (not perfect yet)
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 __del__(self):
136     return
137    
138     def __call__(self, input_data):
139     """
140 jgs 149 Projects input_data onto a continuous function
141 jgs 121
142 jgs 149 @param input_data: The input_data to be projected.
143 jgs 121 """
144 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
145 jgs 121 if input_data.getRank()==0:
146     self.__pde.setValue(Y = input_data)
147     out=self.__pde.getSolution()
148     elif input_data.getRank()==1:
149     for i0 in range(input_data.getShape()[0]):
150     self.__pde.setValue(Y = input_data[i0])
151     out[i0]=self.__pde.getSolution()
152     elif input_data.getRank()==2:
153     for i0 in range(input_data.getShape()[0]):
154     for i1 in range(input_data.getShape()[1]):
155     self.__pde.setValue(Y = input_data[i0,i1])
156     out[i0,i1]=self.__pde.getSolution()
157     elif input_data.getRank()==3:
158     for i0 in range(input_data.getShape()[0]):
159     for i1 in range(input_data.getShape()[1]):
160     for i2 in range(input_data.getShape()[2]):
161     self.__pde.setValue(Y = input_data[i0,i1,i2])
162     out[i0,i1,i2]=self.__pde.getSolution()
163     else:
164     for i0 in range(input_data.getShape()[0]):
165     for i1 in range(input_data.getShape()[1]):
166     for i2 in range(input_data.getShape()[2]):
167     for i3 in range(input_data.getShape()[3]):
168     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
169     out[i0,i1,i2,i3]=self.__pde.getSolution()
170     return out
171    
172 gross 525 class NoPDE:
173     """
174     solves the following problem for u:
175 jgs 121
176 gross 525 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
177    
178     with constraint
179    
180     M{u[j]=r[j]} where M{q[j]>0}
181    
182     where D, Y, r and q are given functions of rank 1.
183    
184     In the case of scalars this takes the form
185    
186     M{D*u=Y}
187    
188     with constraint
189    
190     M{u=r} where M{q>0}
191    
192     where D, Y, r and q are given scalar functions.
193    
194     The constraint is overwriting any other condition.
195    
196 gross 720 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
197     that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
198     thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
199 gross 525 """
200     def __init__(self,domain,D=None,Y=None,q=None,r=None):
201     """
202     initialize the problem
203    
204     @param domain: domain of the PDE.
205     @type domain: L{Domain}
206     @param D: coefficient of the solution.
207 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
208 gross 525 @param Y: right hand side
209 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
210 gross 525 @param q: location of constraints
211 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
212 gross 525 @param r: value of solution at locations of constraints
213 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
214 gross 525 """
215     self.__domain=domain
216     self.__D=D
217     self.__Y=Y
218     self.__q=q
219     self.__r=r
220     self.__u=None
221     self.__function_space=escript.Solution(self.__domain)
222     def setReducedOn(self):
223     """
224     sets the L{FunctionSpace} of the solution to L{ReducedSolution}
225     """
226     self.__function_space=escript.ReducedSolution(self.__domain)
227     self.__u=None
228    
229     def setReducedOff(self):
230     """
231     sets the L{FunctionSpace} of the solution to L{Solution}
232     """
233     self.__function_space=escript.Solution(self.__domain)
234     self.__u=None
235    
236     def setValue(self,D=None,Y=None,q=None,r=None):
237     """
238     assigns values to the parameters.
239    
240     @param D: coefficient of the solution.
241 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
242 gross 525 @param Y: right hand side
243 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
244 gross 525 @param q: location of constraints
245 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
246 gross 525 @param r: value of solution at locations of constraints
247 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
248 gross 525 """
249     if not D==None:
250     self.__D=D
251     self.__u=None
252     if not Y==None:
253     self.__Y=Y
254     self.__u=None
255     if not q==None:
256     self.__q=q
257     self.__u=None
258     if not r==None:
259     self.__r=r
260     self.__u=None
261    
262     def getSolution(self):
263     """
264     returns the solution
265    
266     @return: the solution of the problem
267     @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
268     """
269     if self.__u==None:
270     if self.__D==None:
271     raise ValueError,"coefficient D is undefined"
272     D=escript.Data(self.__D,self.__function_space)
273     if D.getRank()>1:
274     raise ValueError,"coefficient D must have rank 0 or 1"
275     if self.__Y==None:
276     self.__u=escript.Data(0.,D.getShape(),self.__function_space)
277     else:
278     self.__u=util.quotient(self.__Y,D)
279     if not self.__q==None:
280     q=util.wherePositive(escript.Data(self.__q,self.__function_space))
281     self.__u*=(1.-q)
282     if not self.__r==None: self.__u+=q*self.__r
283     return self.__u
284    
285 jgs 147 class Locator:
286     """
287 jgs 149 Locator provides access to the values of data objects at a given
288     spatial coordinate x.
289    
290     In fact, a Locator object finds the sample in the set of samples of a
291     given function space or domain where which is closest to the given
292     point x.
293 jgs 147 """
294    
295     def __init__(self,where,x=numarray.zeros((3,))):
296 jgs 149 """
297     Initializes a Locator to access values in Data objects on the Doamin
298     or FunctionSpace where for the sample point which
299     closest to the given point x.
300     """
301     if isinstance(where,escript.FunctionSpace):
302 jgs 147 self.__function_space=where
303 jgs 121 else:
304 jgs 149 self.__function_space=escript.ContinuousFunction(where)
305 gross 525 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
306 jgs 121
307 jgs 147 def __str__(self):
308 jgs 149 """
309     Returns the coordinates of the Locator as a string.
310     """
311 jgs 147 return "<Locator %s>"%str(self.getX())
312 jgs 121
313 jgs 147 def getFunctionSpace(self):
314 jgs 149 """
315     Returns the function space of the Locator.
316     """
317 jgs 147 return self.__function_space
318    
319 jgs 121 def getId(self):
320 jgs 149 """
321     Returns the identifier of the location.
322     """
323 jgs 121 return self.__id
324    
325 jgs 147 def getX(self):
326 jgs 149 """
327     Returns the exact coordinates of the Locator.
328     """
329 jgs 147 return self(self.getFunctionSpace().getX())
330 jgs 121
331 jgs 147 def __call__(self,data):
332 jgs 149 """
333     Returns the value of data at the Locator of a Data object otherwise
334     the object is returned.
335     """
336 jgs 147 return self.getValue(data)
337 jgs 121
338 jgs 147 def getValue(self,data):
339 jgs 149 """
340     Returns the value of data at the Locator if data is a Data object
341     otherwise the object is returned.
342     """
343     if isinstance(data,escript.Data):
344 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
345 bcumming 790 #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
346     out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
347 jgs 147 else:
348 bcumming 790 #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
349     out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
350 jgs 147 if data.getRank()==0:
351     return out[0]
352     else:
353     return out
354     else:
355     return data
356 jgs 149
357 gross 867 class SaddlePointProblem(object):
358     """
359     This implements a solver for a saddlepoint problem
360    
361 gross 877 M{f(u,p)=0}
362     M{g(u)=0}
363 gross 867
364     for u and p. The problem is solved with an inexact Uszawa scheme for p:
365    
366 gross 877 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})
367     M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
368 gross 867
369     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
370     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'
371     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
372     in fact the role of a preconditioner.
373     """
374     def __init__(self,verbose=False,*args):
375     """
376     initializes the problem
377    
378     @parm verbose: switches on the printing out some information
379     @type verbose: C{bool}
380     @note: this method may be overwritten by a particular saddle point problem
381     """
382     self.__verbose=verbose
383 gross 877 self.relaxation=1.
384 gross 867
385     def trace(self,text):
386     """
387     prints text if verbose has been set
388    
389     @parm text: a text message
390     @type text: C{str}
391     """
392     if self.__verbose: print "%s: %s"%(str(self),text)
393    
394 gross 873 def solve_f(self,u,p,tol=1.e-8):
395 gross 867 """
396     solves
397    
398     A_f du = f(u,p)
399    
400     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
401    
402     @param u: current approximation of u
403     @type u: L{escript.Data}
404     @param p: current approximation of p
405     @type p: L{escript.Data}
406 gross 873 @param tol: tolerance expected for du
407 gross 867 @type tol: C{float}
408     @return: increment du
409     @rtype: L{escript.Data}
410     @note: this method has to be overwritten by a particular saddle point problem
411     """
412     pass
413    
414 gross 873 def solve_g(self,u,tol=1.e-8):
415 gross 867 """
416     solves
417    
418     Q_g dp = g(u)
419    
420     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.
421    
422     @param u: current approximation of u
423     @type u: L{escript.Data}
424 gross 873 @param tol: tolerance expected for dp
425     @type tol: C{float}
426 gross 867 @return: increment dp
427     @rtype: L{escript.Data}
428     @note: this method has to be overwritten by a particular saddle point problem
429     """
430     pass
431    
432     def inner(self,p0,p1):
433     """
434     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
435     @return: inner product of p0 and p1
436     @rtype: C{float}
437     """
438     pass
439    
440 gross 877 subiter_max=3
441     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
442     """
443     runs the solver
444 gross 873
445 gross 877 @param u0: initial guess for C{u}
446     @type u0: L{esys.escript.Data}
447     @param p0: initial guess for C{p}
448     @type p0: L{esys.escript.Data}
449     @param tolerance: tolerance for relative error in C{u} and C{p}
450     @type tolerance: positive C{float}
451     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
452     @type tolerance_u: positive C{float}
453     @param iter_max: maximum number of iteration steps.
454     @type iter_max: C{int}
455     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
456     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
457     @type accepted_reduction: positive C{float} or C{None}
458     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
459     @type relaxation: C{float} or C{None}
460     """
461     tol=1.e-2
462     if tolerance_u==None: tolerance_u=tolerance
463     if not relaxation==None: self.relaxation=relaxation
464     if accepted_reduction ==None:
465     angle_limit=0.
466     elif accepted_reduction>=1.:
467     angle_limit=0.
468     else:
469     angle_limit=util.sqrt(1-accepted_reduction**2)
470     self.iter=0
471     u=u0
472     p=p0
473     #
474     # initialize things:
475     #
476     converged=False
477     #
478     # start loop:
479     #
480     # initial search direction is g
481     #
482     while not converged :
483     if self.iter>iter_max:
484     raise ArithmeticError("no convergence after %s steps."%self.iter)
485     f_new=self.solve_f(u,p,tol)
486     norm_f_new = util.Lsup(f_new)
487     u_new=u-f_new
488     g_new=self.solve_g(u_new,tol)
489     self.iter+=1
490     norm_g_new = util.sqrt(self.inner(g_new,g_new))
491     if norm_f_new==0. and norm_g_new==0.: return u, p
492     if self.iter>1 and not accepted_reduction==None:
493     #
494     # did we manage to reduce the norm of G? I
495     # if not we start a backtracking procedure
496     #
497     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
498     if norm_g_new > accepted_reduction * norm_g:
499     sub_iter=0
500     s=self.relaxation
501     d=g
502     g_last=g
503     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
504     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
505     dg= g_new-g_last
506     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
507     rad=self.inner(g_new,dg)/self.relaxation
508     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
509     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
510     if abs(rad) < norm_dg*norm_g_new * angle_limit:
511     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
512     break
513     r=self.relaxation
514     self.relaxation= - rad/norm_dg**2
515     s+=self.relaxation
516     #####
517     # a=g_new+self.relaxation*dg/r
518     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
519     #####
520     g_last=g_new
521     p+=self.relaxation*d
522     f_new=self.solve_f(u,p,tol)
523     u_new=u-f_new
524     g_new=self.solve_g(u_new,tol)
525     self.iter+=1
526     norm_f_new = util.Lsup(f_new)
527     norm_g_new = util.sqrt(self.inner(g_new,g_new))
528     # print " ",sub_iter," new g norm",norm_g_new
529     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
530     #
531     # can we expect reduction of g?
532     #
533     # u_last=u_new
534     sub_iter+=1
535     self.relaxation=s
536     #
537     # check for convergence:
538     #
539     norm_u_new = util.Lsup(u_new)
540     p_new=p+self.relaxation*g_new
541     norm_p_new = util.sqrt(self.inner(p_new,p_new))
542     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))
543 gross 873
544 gross 877 if self.iter>1:
545     dg2=g_new-g
546     df2=f_new-f
547     norm_dg2=util.sqrt(self.inner(dg2,dg2))
548     norm_df2=util.Lsup(df2)
549     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
550     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
551     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
552     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
553     converged=True
554     break
555     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
556     self.trace("convergence after %s steps."%self.iter)
557     return u,p
558     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
559     # tol=1.e-2
560     # iter=0
561     # converged=False
562     # u=u0*1.
563     # p=p0*1.
564     # while not converged and iter<iter_max:
565     # du=self.solve_f(u,p,tol)
566     # u-=du
567     # norm_du=util.Lsup(du)
568     # norm_u=util.Lsup(u)
569     #
570     # dp=self.relaxation*self.solve_g(u,tol)
571     # p+=dp
572     # norm_dp=util.sqrt(self.inner(dp,dp))
573     # norm_p=util.sqrt(self.inner(p,p))
574     # 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)
575     # iter+=1
576     #
577     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
578     # if converged:
579     # print "convergence after %s steps."%iter
580     # else:
581     # raise ArithmeticError("no convergence after %s steps."%iter)
582     #
583     # return u,p
584 gross 873
585 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