/[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 1387 - (hide annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 8 months ago) by trankine
Original Path: temp/escript/py_src/pdetools.py
File MIME type: text/x-python
File size: 29931 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
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 gross 351 class TimeIntegrationManager:
52     """
53     a simple mechanism to manage time dependend values.
54    
55 gross 720 typical usage is::
56 gross 351
57 gross 720 dt=0.1 # time increment
58     tm=TimeIntegrationManager(inital_value,p=1)
59     while t<1.
60     v_guess=tm.extrapolate(dt) # extrapolate to t+dt
61     v=...
62     tm.checkin(dt,v)
63     t+=dt
64 gross 351
65 gross 720 @note: currently only p=1 is supported.
66 gross 351 """
67     def __init__(self,*inital_values,**kwargs):
68     """
69     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
70     """
71     if kwargs.has_key("p"):
72     self.__p=kwargs["p"]
73     else:
74     self.__p=1
75     if kwargs.has_key("time"):
76     self.__t=kwargs["time"]
77     else:
78     self.__t=0.
79     self.__v_mem=[inital_values]
80     self.__order=0
81     self.__dt_mem=[]
82     self.__num_val=len(inital_values)
83    
84     def getTime(self):
85     return self.__t
86 gross 396 def getValue(self):
87 gross 409 out=self.__v_mem[0]
88     if len(out)==1:
89     return out[0]
90     else:
91     return out
92    
93 gross 351 def checkin(self,dt,*values):
94     """
95     adds new values to the manager. the p+1 last value get lost
96     """
97     o=min(self.__order+1,self.__p)
98     self.__order=min(self.__order+1,self.__p)
99     v_mem_new=[values]
100     dt_mem_new=[dt]
101     for i in range(o-1):
102     v_mem_new.append(self.__v_mem[i])
103     dt_mem_new.append(self.__dt_mem[i])
104     v_mem_new.append(self.__v_mem[o-1])
105     self.__order=o
106     self.__v_mem=v_mem_new
107     self.__dt_mem=dt_mem_new
108     self.__t+=dt
109    
110     def extrapolate(self,dt):
111     """
112     extrapolates to dt forward in time.
113     """
114     if self.__order==0:
115     out=self.__v_mem[0]
116     else:
117     out=[]
118     for i in range(self.__num_val):
119     out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
120    
121     if len(out)==0:
122     return None
123     elif len(out)==1:
124     return out[0]
125     else:
126     return out
127 gross 396
128 gross 351
129 jgs 121 class Projector:
130 jgs 149 """
131     The Projector is a factory which projects a discontiuous function onto a
132     continuous function on the a given domain.
133     """
134 jgs 121 def __init__(self, domain, reduce = True, fast=True):
135     """
136 jgs 149 Create a continuous function space projector for a domain.
137 jgs 121
138 jgs 149 @param domain: Domain of the projection.
139     @param reduce: Flag to reduce projection order (default is True)
140     @param fast: Flag to use a fast method based on matrix lumping (default is true)
141 jgs 121 """
142 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
143 jgs 148 if fast:
144 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
145 jgs 121 self.__pde.setSymmetryOn()
146     self.__pde.setReducedOrderTo(reduce)
147     self.__pde.setValue(D = 1.)
148 ksteube 1312 return
149 jgs 121
150     def __call__(self, input_data):
151     """
152 jgs 149 Projects input_data onto a continuous function
153 jgs 121
154 jgs 149 @param input_data: The input_data to be projected.
155 jgs 121 """
156 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
157 gross 1122 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
158 jgs 121 if input_data.getRank()==0:
159     self.__pde.setValue(Y = input_data)
160     out=self.__pde.getSolution()
161     elif input_data.getRank()==1:
162     for i0 in range(input_data.getShape()[0]):
163     self.__pde.setValue(Y = input_data[i0])
164     out[i0]=self.__pde.getSolution()
165     elif input_data.getRank()==2:
166     for i0 in range(input_data.getShape()[0]):
167     for i1 in range(input_data.getShape()[1]):
168     self.__pde.setValue(Y = input_data[i0,i1])
169     out[i0,i1]=self.__pde.getSolution()
170     elif input_data.getRank()==3:
171     for i0 in range(input_data.getShape()[0]):
172     for i1 in range(input_data.getShape()[1]):
173     for i2 in range(input_data.getShape()[2]):
174     self.__pde.setValue(Y = input_data[i0,i1,i2])
175     out[i0,i1,i2]=self.__pde.getSolution()
176     else:
177     for i0 in range(input_data.getShape()[0]):
178     for i1 in range(input_data.getShape()[1]):
179     for i2 in range(input_data.getShape()[2]):
180     for i3 in range(input_data.getShape()[3]):
181     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
182     out[i0,i1,i2,i3]=self.__pde.getSolution()
183     return out
184    
185 gross 525 class NoPDE:
186     """
187     solves the following problem for u:
188 jgs 121
189 gross 525 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
190    
191     with constraint
192    
193     M{u[j]=r[j]} where M{q[j]>0}
194    
195     where D, Y, r and q are given functions of rank 1.
196    
197     In the case of scalars this takes the form
198    
199     M{D*u=Y}
200    
201     with constraint
202    
203     M{u=r} where M{q>0}
204    
205     where D, Y, r and q are given scalar functions.
206    
207     The constraint is overwriting any other condition.
208    
209 gross 720 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
210     that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
211     thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
212 gross 525 """
213     def __init__(self,domain,D=None,Y=None,q=None,r=None):
214     """
215     initialize the problem
216    
217     @param domain: domain of the PDE.
218     @type domain: L{Domain}
219     @param D: coefficient of the solution.
220 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
221 gross 525 @param Y: right hand side
222 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
223 gross 525 @param q: location of constraints
224 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
225 gross 525 @param r: value of solution at locations of constraints
226 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
227 gross 525 """
228     self.__domain=domain
229     self.__D=D
230     self.__Y=Y
231     self.__q=q
232     self.__r=r
233     self.__u=None
234     self.__function_space=escript.Solution(self.__domain)
235     def setReducedOn(self):
236     """
237     sets the L{FunctionSpace} of the solution to L{ReducedSolution}
238     """
239     self.__function_space=escript.ReducedSolution(self.__domain)
240     self.__u=None
241    
242     def setReducedOff(self):
243     """
244     sets the L{FunctionSpace} of the solution to L{Solution}
245     """
246     self.__function_space=escript.Solution(self.__domain)
247     self.__u=None
248    
249     def setValue(self,D=None,Y=None,q=None,r=None):
250     """
251     assigns values to the parameters.
252    
253     @param D: coefficient of the solution.
254 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
255 gross 525 @param Y: right hand side
256 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
257 gross 525 @param q: location of constraints
258 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
259 gross 525 @param r: value of solution at locations of constraints
260 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
261 gross 525 """
262     if not D==None:
263     self.__D=D
264     self.__u=None
265     if not Y==None:
266     self.__Y=Y
267     self.__u=None
268     if not q==None:
269     self.__q=q
270     self.__u=None
271     if not r==None:
272     self.__r=r
273     self.__u=None
274    
275     def getSolution(self):
276     """
277     returns the solution
278    
279     @return: the solution of the problem
280     @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
281     """
282     if self.__u==None:
283     if self.__D==None:
284     raise ValueError,"coefficient D is undefined"
285     D=escript.Data(self.__D,self.__function_space)
286     if D.getRank()>1:
287     raise ValueError,"coefficient D must have rank 0 or 1"
288     if self.__Y==None:
289     self.__u=escript.Data(0.,D.getShape(),self.__function_space)
290     else:
291     self.__u=util.quotient(self.__Y,D)
292     if not self.__q==None:
293     q=util.wherePositive(escript.Data(self.__q,self.__function_space))
294     self.__u*=(1.-q)
295     if not self.__r==None: self.__u+=q*self.__r
296     return self.__u
297    
298 jgs 147 class Locator:
299     """
300 jgs 149 Locator provides access to the values of data objects at a given
301     spatial coordinate x.
302    
303     In fact, a Locator object finds the sample in the set of samples of a
304     given function space or domain where which is closest to the given
305     point x.
306 jgs 147 """
307    
308     def __init__(self,where,x=numarray.zeros((3,))):
309 jgs 149 """
310     Initializes a Locator to access values in Data objects on the Doamin
311     or FunctionSpace where for the sample point which
312     closest to the given point x.
313 gross 880
314     @param where: function space
315     @type where: L{escript.FunctionSpace}
316     @param x: coefficient of the solution.
317     @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
318 jgs 149 """
319     if isinstance(where,escript.FunctionSpace):
320 jgs 147 self.__function_space=where
321 jgs 121 else:
322 jgs 149 self.__function_space=escript.ContinuousFunction(where)
323 gross 880 if isinstance(x, list):
324     self.__id=[]
325     for p in x:
326 gross 921 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
327 gross 880 else:
328 gross 921 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
329 jgs 121
330 jgs 147 def __str__(self):
331 jgs 149 """
332     Returns the coordinates of the Locator as a string.
333     """
334 gross 880 x=self.getX()
335     if instance(x,list):
336     out="["
337     first=True
338     for xx in x:
339     if not first:
340     out+=","
341     else:
342     first=False
343     out+=str(xx)
344     out+="]>"
345     else:
346     out=str(x)
347     return out
348 jgs 121
349 gross 880 def getX(self):
350     """
351     Returns the exact coordinates of the Locator.
352     """
353     return self(self.getFunctionSpace().getX())
354    
355 jgs 147 def getFunctionSpace(self):
356 jgs 149 """
357     Returns the function space of the Locator.
358     """
359 jgs 147 return self.__function_space
360    
361 gross 880 def getId(self,item=None):
362 jgs 149 """
363     Returns the identifier of the location.
364     """
365 gross 880 if item == None:
366     return self.__id
367     else:
368     if isinstance(self.__id,list):
369     return self.__id[item]
370     else:
371     return self.__id
372 jgs 121
373    
374 jgs 147 def __call__(self,data):
375 jgs 149 """
376     Returns the value of data at the Locator of a Data object otherwise
377     the object is returned.
378     """
379 jgs 147 return self.getValue(data)
380 jgs 121
381 jgs 147 def getValue(self,data):
382 jgs 149 """
383     Returns the value of data at the Locator if data is a Data object
384     otherwise the object is returned.
385     """
386     if isinstance(data,escript.Data):
387 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
388 gross 880 dat=data
389 jgs 147 else:
390 gross 880 dat=data.interpolate(self.getFunctionSpace())
391     id=self.getId()
392     r=data.getRank()
393     if isinstance(id,list):
394     out=[]
395     for i in id:
396 gross 921 o=data.getValueOfGlobalDataPoint(*i)
397 gross 880 if data.getRank()==0:
398     out.append(o[0])
399     else:
400     out.append(o)
401     return out
402 jgs 147 else:
403 gross 921 out=data.getValueOfGlobalDataPoint(*id)
404 gross 880 if data.getRank()==0:
405     return out[0]
406     else:
407     return out
408 jgs 147 else:
409     return data
410 jgs 149
411 ksteube 1312 class SolverSchemeException(Exception):
412     """
413     exceptions thrown by solvers
414     """
415     pass
416    
417     class IndefinitePreconditioner(SolverSchemeException):
418     """
419     the preconditioner is not positive definite.
420     """
421     pass
422     class MaxIterReached(SolverSchemeException):
423     """
424     maxium number of iteration steps is reached.
425     """
426     pass
427     class IterationBreakDown(SolverSchemeException):
428     """
429     iteration scheme econouters an incurable breakdown.
430     """
431     pass
432     class NegativeNorm(SolverSchemeException):
433     """
434     a norm calculation returns a negative norm.
435     """
436     pass
437    
438 gross 1330 class IterationHistory(object):
439 ksteube 1312 """
440 gross 1330 The IterationHistory class is used to define a stopping criterium. It keeps track of the
441     residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
442     a given tolerance.
443     """
444     def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
445     """
446     Initialization
447    
448     @param tolerance: tolerance
449     @type tolerance: positive C{float}
450     @param verbose: switches on the printing out some information
451     @type verbose: C{bool}
452     """
453     if not tolerance>0.:
454     raise ValueError,"tolerance needs to be positive."
455     self.tolerance=tolerance
456     self.verbose=verbose
457     self.history=[]
458     def stoppingcriterium(self,norm_r,r,x):
459     """
460     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.
461    
462    
463     @param norm_r: current residual norm
464     @type norm_r: non-negative C{float}
465     @param r: current residual (not used)
466     @param x: current solution approximation (not used)
467     @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
468     @rtype: C{bool}
469    
470     """
471     self.history.append(norm_r)
472     if self.verbose: print "iter: %s: inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])
473     return self.history[-1]<=self.tolerance * self.history[0]
474    
475     def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
476     """
477 ksteube 1312 Solver for
478    
479     M{Ax=b}
480    
481     with a symmetric and positive definite operator A (more details required!).
482     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
483    
484 gross 1330 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
485 ksteube 1312
486     For details on the preconditioned conjugate gradient method see the book:
487    
488     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
489     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
490     C. Romine, and H. van der Vorst.
491    
492     @param b: the right hand side of the liner system. C{b} is altered.
493 gross 1330 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
494 ksteube 1312 @param Aprod: returns the value Ax
495 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}.
496 ksteube 1312 @param Msolve: solves Mx=r
497 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
498     type like argument C{x}.
499 ksteube 1312 @param bilinearform: inner product C{<x,r>}
500 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}.
501     @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.
502     @type stoppingcriterium: function that returns C{True} or C{False}
503     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
504     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
505 ksteube 1312 @param iter_max: maximum number of iteration steps.
506     @type iter_max: C{int}
507 gross 1330 @return: the solution approximation and the corresponding residual
508     @rtype: C{tuple}
509     @warning: C{b} and C{x} are altered.
510 ksteube 1312 """
511     iter=0
512 gross 1330 if x==None:
513     x=0*b
514     else:
515     b += (-1)*Aprod(x)
516 ksteube 1312 r=b
517     rhat=Msolve(r)
518 gross 1330 d = rhat
519 ksteube 1312 rhat_dot_r = bilinearform(rhat, r)
520 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521 ksteube 1312
522 gross 1330 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
523     iter+=1
524 ksteube 1312 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
525    
526     q=Aprod(d)
527     alpha = rhat_dot_r / bilinearform(d, q)
528     x += alpha * d
529     r += (-alpha) * q
530    
531     rhat=Msolve(r)
532     rhat_dot_r_new = bilinearform(rhat, r)
533     beta = rhat_dot_r_new / rhat_dot_r
534     rhat+=beta * d
535     d=rhat
536    
537     rhat_dot_r = rhat_dot_r_new
538 gross 1330 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
539 ksteube 1312
540 gross 1330 return x,r
541 ksteube 1312
542 gross 1331 class ArithmeticTuple(object):
543     """
544     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
545    
546     example of usage:
547    
548     from esys.escript import Data
549     from numarray import array
550     a=Data(...)
551     b=array([1.,4.])
552     x=ArithmeticTuple(a,b)
553     y=5.*x
554    
555     """
556     def __init__(self,*args):
557     """
558     initialize object with elements args.
559    
560     @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
561     """
562     self.__items=list(args)
563    
564     def __len__(self):
565     """
566     number of items
567    
568     @return: number of items
569     @rtype: C{int}
570     """
571     return len(self.__items)
572    
573     def __getitem__(self,index):
574     """
575     get an item
576    
577     @param index: item to be returned
578     @type index: C{int}
579     @return: item with index C{index}
580     """
581     return self.__items.__getitem__(index)
582    
583     def __mul__(self,other):
584     """
585     scaling from the right
586    
587     @param other: scaling factor
588     @type other: C{float}
589     @return: itemwise self*other
590     @rtype: L{ArithmeticTuple}
591     """
592     out=[]
593     for i in range(len(self)):
594     out.append(self[i]*other)
595     return ArithmeticTuple(*tuple(out))
596    
597     def __rmul__(self,other):
598     """
599     scaling from the left
600    
601     @param other: scaling factor
602     @type other: C{float}
603     @return: itemwise other*self
604     @rtype: L{ArithmeticTuple}
605     """
606     out=[]
607     for i in range(len(self)):
608     out.append(other*self[i])
609     return ArithmeticTuple(*tuple(out))
610    
611     def __iadd__(self,other):
612     """
613     in-place add of other to self
614    
615     @param other: increment
616     @type other: C{ArithmeticTuple}
617     """
618     if len(self) != len(other):
619     raise ValueError,"tuple length must match."
620     for i in range(len(self)):
621     self.__items[i]+=other[i]
622     return self
623    
624 gross 867 class SaddlePointProblem(object):
625     """
626     This implements a solver for a saddlepoint problem
627    
628 gross 877 M{f(u,p)=0}
629     M{g(u)=0}
630 gross 867
631     for u and p. The problem is solved with an inexact Uszawa scheme for p:
632    
633 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
634 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
635 gross 867
636     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
637     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'
638     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
639     in fact the role of a preconditioner.
640     """
641     def __init__(self,verbose=False,*args):
642     """
643     initializes the problem
644    
645 ksteube 990 @param verbose: switches on the printing out some information
646 gross 867 @type verbose: C{bool}
647     @note: this method may be overwritten by a particular saddle point problem
648     """
649 gross 1107 if not isinstance(verbose,bool):
650     raise TypeError("verbose needs to be of type bool.")
651 gross 1106 self.__verbose=verbose
652 gross 877 self.relaxation=1.
653 gross 867
654     def trace(self,text):
655     """
656     prints text if verbose has been set
657    
658 ksteube 990 @param text: a text message
659 gross 867 @type text: C{str}
660     """
661     if self.__verbose: print "%s: %s"%(str(self),text)
662    
663 gross 873 def solve_f(self,u,p,tol=1.e-8):
664 gross 867 """
665     solves
666    
667     A_f du = f(u,p)
668    
669     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
670    
671     @param u: current approximation of u
672     @type u: L{escript.Data}
673     @param p: current approximation of p
674     @type p: L{escript.Data}
675 gross 873 @param tol: tolerance expected for du
676 gross 867 @type tol: C{float}
677     @return: increment du
678     @rtype: L{escript.Data}
679     @note: this method has to be overwritten by a particular saddle point problem
680     """
681     pass
682    
683 gross 873 def solve_g(self,u,tol=1.e-8):
684 gross 867 """
685     solves
686    
687     Q_g dp = g(u)
688    
689     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.
690    
691     @param u: current approximation of u
692     @type u: L{escript.Data}
693 gross 873 @param tol: tolerance expected for dp
694     @type tol: C{float}
695 gross 867 @return: increment dp
696     @rtype: L{escript.Data}
697     @note: this method has to be overwritten by a particular saddle point problem
698     """
699     pass
700    
701     def inner(self,p0,p1):
702     """
703     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
704     @return: inner product of p0 and p1
705     @rtype: C{float}
706     """
707     pass
708    
709 gross 877 subiter_max=3
710     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
711     """
712     runs the solver
713 gross 873
714 gross 877 @param u0: initial guess for C{u}
715     @type u0: L{esys.escript.Data}
716     @param p0: initial guess for C{p}
717     @type p0: L{esys.escript.Data}
718     @param tolerance: tolerance for relative error in C{u} and C{p}
719     @type tolerance: positive C{float}
720     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
721     @type tolerance_u: positive C{float}
722     @param iter_max: maximum number of iteration steps.
723     @type iter_max: C{int}
724     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
725     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
726     @type accepted_reduction: positive C{float} or C{None}
727     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
728     @type relaxation: C{float} or C{None}
729     """
730     tol=1.e-2
731     if tolerance_u==None: tolerance_u=tolerance
732     if not relaxation==None: self.relaxation=relaxation
733     if accepted_reduction ==None:
734     angle_limit=0.
735     elif accepted_reduction>=1.:
736     angle_limit=0.
737     else:
738     angle_limit=util.sqrt(1-accepted_reduction**2)
739     self.iter=0
740     u=u0
741     p=p0
742     #
743     # initialize things:
744     #
745     converged=False
746     #
747     # start loop:
748     #
749     # initial search direction is g
750     #
751     while not converged :
752     if self.iter>iter_max:
753     raise ArithmeticError("no convergence after %s steps."%self.iter)
754     f_new=self.solve_f(u,p,tol)
755     norm_f_new = util.Lsup(f_new)
756     u_new=u-f_new
757     g_new=self.solve_g(u_new,tol)
758     self.iter+=1
759     norm_g_new = util.sqrt(self.inner(g_new,g_new))
760     if norm_f_new==0. and norm_g_new==0.: return u, p
761     if self.iter>1 and not accepted_reduction==None:
762     #
763     # did we manage to reduce the norm of G? I
764     # if not we start a backtracking procedure
765     #
766     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
767     if norm_g_new > accepted_reduction * norm_g:
768     sub_iter=0
769     s=self.relaxation
770     d=g
771     g_last=g
772     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
773     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
774     dg= g_new-g_last
775     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
776     rad=self.inner(g_new,dg)/self.relaxation
777     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
778     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
779     if abs(rad) < norm_dg*norm_g_new * angle_limit:
780     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
781     break
782     r=self.relaxation
783     self.relaxation= - rad/norm_dg**2
784     s+=self.relaxation
785     #####
786     # a=g_new+self.relaxation*dg/r
787     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
788     #####
789     g_last=g_new
790     p+=self.relaxation*d
791     f_new=self.solve_f(u,p,tol)
792     u_new=u-f_new
793     g_new=self.solve_g(u_new,tol)
794     self.iter+=1
795     norm_f_new = util.Lsup(f_new)
796     norm_g_new = util.sqrt(self.inner(g_new,g_new))
797     # print " ",sub_iter," new g norm",norm_g_new
798     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
799     #
800     # can we expect reduction of g?
801     #
802     # u_last=u_new
803     sub_iter+=1
804     self.relaxation=s
805     #
806     # check for convergence:
807     #
808     norm_u_new = util.Lsup(u_new)
809     p_new=p+self.relaxation*g_new
810     norm_p_new = util.sqrt(self.inner(p_new,p_new))
811 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))
812 gross 873
813 gross 877 if self.iter>1:
814     dg2=g_new-g
815     df2=f_new-f
816     norm_dg2=util.sqrt(self.inner(dg2,dg2))
817     norm_df2=util.Lsup(df2)
818     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
819     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
820     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
821     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
822     converged=True
823     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
824     self.trace("convergence after %s steps."%self.iter)
825     return u,p
826     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
827     # tol=1.e-2
828     # iter=0
829     # converged=False
830     # u=u0*1.
831     # p=p0*1.
832     # while not converged and iter<iter_max:
833     # du=self.solve_f(u,p,tol)
834     # u-=du
835     # norm_du=util.Lsup(du)
836     # norm_u=util.Lsup(u)
837     #
838     # dp=self.relaxation*self.solve_g(u,tol)
839     # p+=dp
840     # norm_dp=util.sqrt(self.inner(dp,dp))
841     # norm_p=util.sqrt(self.inner(p,p))
842     # 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)
843     # iter+=1
844     #
845     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
846     # if converged:
847     # print "convergence after %s steps."%iter
848     # else:
849     # raise ArithmeticError("no convergence after %s steps."%iter)
850     #
851     # return u,p
852 gross 873
853 ksteube 1312 def MaskFromBoundaryTag(function_space,*tags):
854     """
855     create a mask on the given function space which one for samples
856     that touch the boundary tagged by tags.
857    
858     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
859    
860     @param function_space: a given function space
861     @type function_space: L{escript.FunctionSpace}
862     @param tags: boundray tags
863     @type tags: C{str}
864     @return: a mask which marks samples used by C{function_space} that are touching the
865     boundary tagged by any of the given tags.
866     @rtype: L{escript.Data} of rank 0
867     """
868     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
869     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
870     for t in tags: d.setTaggedValue(t,1.)
871     pde.setValue(y=d)
872     out=util.whereNonZero(pde.getRightHandSide())
873     if out.getFunctionSpace() == function_space:
874     return out
875     else:
876     return util.whereNonZero(util.interpolate(out,function_space))
877    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26