/[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 1122 - (hide annotations)
Tue May 1 03:21:04 2007 UTC (12 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 21624 byte(s)
As the LinearPDE uses coefficients as reduced even if they are handed
over as full the projector runs into a problem when reduced and full
arguments are used in the same projector. Now the projector resets the
coefficients befor starting the projection.


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