/[escript]/temp/escript/py_src/pdetools.py
ViewVC logotype

Annotation of /temp/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1330 - (hide annotations)
Mon Oct 22 04:54:49 2007 UTC (11 years, 10 months ago) by gross
Original Path: trunk/escript/py_src/pdetools.py
File MIME type: text/x-python
File size: 27987 byte(s)
more flexible stopping criterium for the PCG
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 867 class SaddlePointProblem(object):
543     """
544     This implements a solver for a saddlepoint problem
545    
546 gross 877 M{f(u,p)=0}
547     M{g(u)=0}
548 gross 867
549     for u and p. The problem is solved with an inexact Uszawa scheme for p:
550    
551 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
552 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
553 gross 867
554     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
555     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'
556     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
557     in fact the role of a preconditioner.
558     """
559     def __init__(self,verbose=False,*args):
560     """
561     initializes the problem
562    
563 ksteube 990 @param verbose: switches on the printing out some information
564 gross 867 @type verbose: C{bool}
565     @note: this method may be overwritten by a particular saddle point problem
566     """
567 gross 1107 if not isinstance(verbose,bool):
568     raise TypeError("verbose needs to be of type bool.")
569 gross 1106 self.__verbose=verbose
570 gross 877 self.relaxation=1.
571 gross 867
572     def trace(self,text):
573     """
574     prints text if verbose has been set
575    
576 ksteube 990 @param text: a text message
577 gross 867 @type text: C{str}
578     """
579     if self.__verbose: print "%s: %s"%(str(self),text)
580    
581 gross 873 def solve_f(self,u,p,tol=1.e-8):
582 gross 867 """
583     solves
584    
585     A_f du = f(u,p)
586    
587     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
588    
589     @param u: current approximation of u
590     @type u: L{escript.Data}
591     @param p: current approximation of p
592     @type p: L{escript.Data}
593 gross 873 @param tol: tolerance expected for du
594 gross 867 @type tol: C{float}
595     @return: increment du
596     @rtype: L{escript.Data}
597     @note: this method has to be overwritten by a particular saddle point problem
598     """
599     pass
600    
601 gross 873 def solve_g(self,u,tol=1.e-8):
602 gross 867 """
603     solves
604    
605     Q_g dp = g(u)
606    
607     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.
608    
609     @param u: current approximation of u
610     @type u: L{escript.Data}
611 gross 873 @param tol: tolerance expected for dp
612     @type tol: C{float}
613 gross 867 @return: increment dp
614     @rtype: L{escript.Data}
615     @note: this method has to be overwritten by a particular saddle point problem
616     """
617     pass
618    
619     def inner(self,p0,p1):
620     """
621     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
622     @return: inner product of p0 and p1
623     @rtype: C{float}
624     """
625     pass
626    
627 gross 877 subiter_max=3
628     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
629     """
630     runs the solver
631 gross 873
632 gross 877 @param u0: initial guess for C{u}
633     @type u0: L{esys.escript.Data}
634     @param p0: initial guess for C{p}
635     @type p0: L{esys.escript.Data}
636     @param tolerance: tolerance for relative error in C{u} and C{p}
637     @type tolerance: positive C{float}
638     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
639     @type tolerance_u: positive C{float}
640     @param iter_max: maximum number of iteration steps.
641     @type iter_max: C{int}
642     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
643     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
644     @type accepted_reduction: positive C{float} or C{None}
645     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
646     @type relaxation: C{float} or C{None}
647     """
648     tol=1.e-2
649     if tolerance_u==None: tolerance_u=tolerance
650     if not relaxation==None: self.relaxation=relaxation
651     if accepted_reduction ==None:
652     angle_limit=0.
653     elif accepted_reduction>=1.:
654     angle_limit=0.
655     else:
656     angle_limit=util.sqrt(1-accepted_reduction**2)
657     self.iter=0
658     u=u0
659     p=p0
660     #
661     # initialize things:
662     #
663     converged=False
664     #
665     # start loop:
666     #
667     # initial search direction is g
668     #
669     while not converged :
670     if self.iter>iter_max:
671     raise ArithmeticError("no convergence after %s steps."%self.iter)
672     f_new=self.solve_f(u,p,tol)
673     norm_f_new = util.Lsup(f_new)
674     u_new=u-f_new
675     g_new=self.solve_g(u_new,tol)
676     self.iter+=1
677     norm_g_new = util.sqrt(self.inner(g_new,g_new))
678     if norm_f_new==0. and norm_g_new==0.: return u, p
679     if self.iter>1 and not accepted_reduction==None:
680     #
681     # did we manage to reduce the norm of G? I
682     # if not we start a backtracking procedure
683     #
684     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
685     if norm_g_new > accepted_reduction * norm_g:
686     sub_iter=0
687     s=self.relaxation
688     d=g
689     g_last=g
690     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
691     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
692     dg= g_new-g_last
693     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
694     rad=self.inner(g_new,dg)/self.relaxation
695     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
696     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
697     if abs(rad) < norm_dg*norm_g_new * angle_limit:
698     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
699     break
700     r=self.relaxation
701     self.relaxation= - rad/norm_dg**2
702     s+=self.relaxation
703     #####
704     # a=g_new+self.relaxation*dg/r
705     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
706     #####
707     g_last=g_new
708     p+=self.relaxation*d
709     f_new=self.solve_f(u,p,tol)
710     u_new=u-f_new
711     g_new=self.solve_g(u_new,tol)
712     self.iter+=1
713     norm_f_new = util.Lsup(f_new)
714     norm_g_new = util.sqrt(self.inner(g_new,g_new))
715     # print " ",sub_iter," new g norm",norm_g_new
716     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
717     #
718     # can we expect reduction of g?
719     #
720     # u_last=u_new
721     sub_iter+=1
722     self.relaxation=s
723     #
724     # check for convergence:
725     #
726     norm_u_new = util.Lsup(u_new)
727     p_new=p+self.relaxation*g_new
728     norm_p_new = util.sqrt(self.inner(p_new,p_new))
729 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))
730 gross 873
731 gross 877 if self.iter>1:
732     dg2=g_new-g
733     df2=f_new-f
734     norm_dg2=util.sqrt(self.inner(dg2,dg2))
735     norm_df2=util.Lsup(df2)
736     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
737     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
738     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
739     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
740     converged=True
741     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
742     self.trace("convergence after %s steps."%self.iter)
743     return u,p
744     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
745     # tol=1.e-2
746     # iter=0
747     # converged=False
748     # u=u0*1.
749     # p=p0*1.
750     # while not converged and iter<iter_max:
751     # du=self.solve_f(u,p,tol)
752     # u-=du
753     # norm_du=util.Lsup(du)
754     # norm_u=util.Lsup(u)
755     #
756     # dp=self.relaxation*self.solve_g(u,tol)
757     # p+=dp
758     # norm_dp=util.sqrt(self.inner(dp,dp))
759     # norm_p=util.sqrt(self.inner(p,p))
760     # 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)
761     # iter+=1
762     #
763     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
764     # if converged:
765     # print "convergence after %s steps."%iter
766     # else:
767     # raise ArithmeticError("no convergence after %s steps."%iter)
768     #
769     # return u,p
770 gross 873
771 ksteube 1312 def MaskFromBoundaryTag(function_space,*tags):
772     """
773     create a mask on the given function space which one for samples
774     that touch the boundary tagged by tags.
775    
776     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
777    
778     @param function_space: a given function space
779     @type function_space: L{escript.FunctionSpace}
780     @param tags: boundray tags
781     @type tags: C{str}
782     @return: a mask which marks samples used by C{function_space} that are touching the
783     boundary tagged by any of the given tags.
784     @rtype: L{escript.Data} of rank 0
785     """
786     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
787     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
788     for t in tags: d.setTaggedValue(t,1.)
789     pde.setValue(y=d)
790     out=util.whereNonZero(pde.getRightHandSide())
791     if out.getFunctionSpace() == function_space:
792     return out
793     else:
794     return util.whereNonZero(util.interpolate(out,function_space))
795    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26