/[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 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (11 years, 9 months ago) by ksteube
File MIME type: text/x-python
File size: 26672 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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     def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):
439     """
440     Solver for
441    
442     M{Ax=b}
443    
444     with a symmetric and positive definite operator A (more details required!).
445     It uses the conjugate gradient method with preconditioner M providing an approximation of A.
446    
447     The iteration is terminated if
448    
449     M{norm(r) <= tolerance * norm(b)}
450    
451     where C{norm()} defines a norm and
452    
453     M{r = b-Ax}
454    
455     the residual.
456    
457     For details on the preconditioned conjugate gradient method see the book:
458    
459     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
460     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
461     C. Romine, and H. van der Vorst.
462    
463     @param b: the right hand side of the liner system. C{b} is altered.
464     @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)
465     @param x: an initial guess for the solution
466     @type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)
467     @param Aprod: returns the value Ax
468     @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.
469     @param Msolve: solves Mx=r
470     @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.
471     @param bilinearform: inner product C{<x,r>}
472     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.
473     @param norm: norm calculation for the residual C{r=b-Ax}.
474     @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.
475     @param verbose: switches on the printing out some information
476     @type verbose: C{bool}
477     @param iter_max: maximum number of iteration steps.
478     @type iter_max: C{int}
479     @param tolerance: tolerance
480     @type tolerance: positive C{float}
481     @return: the solution apprximation and the corresponding residual
482     @rtype: C{tuple} of an S type and and an R type object.A
483     @warning: C{b} ans C{x} are altered.
484     """
485     if verbose:
486     print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance =%e"%(iter_max,tolerance)
487     iter=0
488     normb = norm(b)
489     if normb<0: raise NegativeNorm
490    
491     b += (-1)*Aprod(x)
492     r=b
493     rhat=Msolve(r)
494     d = rhat;
495     rhat_dot_r = bilinearform(rhat, r)
496    
497     while True:
498     normr=norm(r)
499     if normr<0: raise NegativeNorm
500     if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)
501     if normr <= tolerance * normb: return x,r
502    
503     iter+=1 # k = iter = 1 first time through
504     if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
505    
506     q=Aprod(d)
507     alpha = rhat_dot_r / bilinearform(d, q)
508     x += alpha * d
509     r += (-alpha) * q
510    
511     rhat=Msolve(r)
512     rhat_dot_r_new = bilinearform(rhat, r)
513     beta = rhat_dot_r_new / rhat_dot_r
514     rhat+=beta * d
515     d=rhat
516    
517     rhat_dot_r = rhat_dot_r_new
518    
519    
520 gross 867 class SaddlePointProblem(object):
521     """
522     This implements a solver for a saddlepoint problem
523    
524 gross 877 M{f(u,p)=0}
525     M{g(u)=0}
526 gross 867
527     for u and p. The problem is solved with an inexact Uszawa scheme for p:
528    
529 ksteube 990 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
530 gross 877 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
531 gross 867
532     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
533     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'
534     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
535     in fact the role of a preconditioner.
536     """
537     def __init__(self,verbose=False,*args):
538     """
539     initializes the problem
540    
541 ksteube 990 @param verbose: switches on the printing out some information
542 gross 867 @type verbose: C{bool}
543     @note: this method may be overwritten by a particular saddle point problem
544     """
545 gross 1107 if not isinstance(verbose,bool):
546     raise TypeError("verbose needs to be of type bool.")
547 gross 1106 self.__verbose=verbose
548 gross 877 self.relaxation=1.
549 gross 867
550     def trace(self,text):
551     """
552     prints text if verbose has been set
553    
554 ksteube 990 @param text: a text message
555 gross 867 @type text: C{str}
556     """
557     if self.__verbose: print "%s: %s"%(str(self),text)
558    
559 gross 873 def solve_f(self,u,p,tol=1.e-8):
560 gross 867 """
561     solves
562    
563     A_f du = f(u,p)
564    
565     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
566    
567     @param u: current approximation of u
568     @type u: L{escript.Data}
569     @param p: current approximation of p
570     @type p: L{escript.Data}
571 gross 873 @param tol: tolerance expected for du
572 gross 867 @type tol: C{float}
573     @return: increment du
574     @rtype: L{escript.Data}
575     @note: this method has to be overwritten by a particular saddle point problem
576     """
577     pass
578    
579 gross 873 def solve_g(self,u,tol=1.e-8):
580 gross 867 """
581     solves
582    
583     Q_g dp = g(u)
584    
585     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.
586    
587     @param u: current approximation of u
588     @type u: L{escript.Data}
589 gross 873 @param tol: tolerance expected for dp
590     @type tol: C{float}
591 gross 867 @return: increment dp
592     @rtype: L{escript.Data}
593     @note: this method has to be overwritten by a particular saddle point problem
594     """
595     pass
596    
597     def inner(self,p0,p1):
598     """
599     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
600     @return: inner product of p0 and p1
601     @rtype: C{float}
602     """
603     pass
604    
605 gross 877 subiter_max=3
606     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
607     """
608     runs the solver
609 gross 873
610 gross 877 @param u0: initial guess for C{u}
611     @type u0: L{esys.escript.Data}
612     @param p0: initial guess for C{p}
613     @type p0: L{esys.escript.Data}
614     @param tolerance: tolerance for relative error in C{u} and C{p}
615     @type tolerance: positive C{float}
616     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
617     @type tolerance_u: positive C{float}
618     @param iter_max: maximum number of iteration steps.
619     @type iter_max: C{int}
620     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
621     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
622     @type accepted_reduction: positive C{float} or C{None}
623     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
624     @type relaxation: C{float} or C{None}
625     """
626     tol=1.e-2
627     if tolerance_u==None: tolerance_u=tolerance
628     if not relaxation==None: self.relaxation=relaxation
629     if accepted_reduction ==None:
630     angle_limit=0.
631     elif accepted_reduction>=1.:
632     angle_limit=0.
633     else:
634     angle_limit=util.sqrt(1-accepted_reduction**2)
635     self.iter=0
636     u=u0
637     p=p0
638     #
639     # initialize things:
640     #
641     converged=False
642     #
643     # start loop:
644     #
645     # initial search direction is g
646     #
647     while not converged :
648     if self.iter>iter_max:
649     raise ArithmeticError("no convergence after %s steps."%self.iter)
650     f_new=self.solve_f(u,p,tol)
651     norm_f_new = util.Lsup(f_new)
652     u_new=u-f_new
653     g_new=self.solve_g(u_new,tol)
654     self.iter+=1
655     norm_g_new = util.sqrt(self.inner(g_new,g_new))
656     if norm_f_new==0. and norm_g_new==0.: return u, p
657     if self.iter>1 and not accepted_reduction==None:
658     #
659     # did we manage to reduce the norm of G? I
660     # if not we start a backtracking procedure
661     #
662     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
663     if norm_g_new > accepted_reduction * norm_g:
664     sub_iter=0
665     s=self.relaxation
666     d=g
667     g_last=g
668     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
669     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
670     dg= g_new-g_last
671     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
672     rad=self.inner(g_new,dg)/self.relaxation
673     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
674     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
675     if abs(rad) < norm_dg*norm_g_new * angle_limit:
676     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
677     break
678     r=self.relaxation
679     self.relaxation= - rad/norm_dg**2
680     s+=self.relaxation
681     #####
682     # a=g_new+self.relaxation*dg/r
683     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
684     #####
685     g_last=g_new
686     p+=self.relaxation*d
687     f_new=self.solve_f(u,p,tol)
688     u_new=u-f_new
689     g_new=self.solve_g(u_new,tol)
690     self.iter+=1
691     norm_f_new = util.Lsup(f_new)
692     norm_g_new = util.sqrt(self.inner(g_new,g_new))
693     # print " ",sub_iter," new g norm",norm_g_new
694     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
695     #
696     # can we expect reduction of g?
697     #
698     # u_last=u_new
699     sub_iter+=1
700     self.relaxation=s
701     #
702     # check for convergence:
703     #
704     norm_u_new = util.Lsup(u_new)
705     p_new=p+self.relaxation*g_new
706     norm_p_new = util.sqrt(self.inner(p_new,p_new))
707 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))
708 gross 873
709 gross 877 if self.iter>1:
710     dg2=g_new-g
711     df2=f_new-f
712     norm_dg2=util.sqrt(self.inner(dg2,dg2))
713     norm_df2=util.Lsup(df2)
714     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
715     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
716     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
717     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
718     converged=True
719     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
720     self.trace("convergence after %s steps."%self.iter)
721     return u,p
722     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
723     # tol=1.e-2
724     # iter=0
725     # converged=False
726     # u=u0*1.
727     # p=p0*1.
728     # while not converged and iter<iter_max:
729     # du=self.solve_f(u,p,tol)
730     # u-=du
731     # norm_du=util.Lsup(du)
732     # norm_u=util.Lsup(u)
733     #
734     # dp=self.relaxation*self.solve_g(u,tol)
735     # p+=dp
736     # norm_dp=util.sqrt(self.inner(dp,dp))
737     # norm_p=util.sqrt(self.inner(p,p))
738     # 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)
739     # iter+=1
740     #
741     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
742     # if converged:
743     # print "convergence after %s steps."%iter
744     # else:
745     # raise ArithmeticError("no convergence after %s steps."%iter)
746     #
747     # return u,p
748 gross 873
749 ksteube 1312 def MaskFromBoundaryTag(function_space,*tags):
750     """
751     create a mask on the given function space which one for samples
752     that touch the boundary tagged by tags.
753    
754     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
755    
756     @param function_space: a given function space
757     @type function_space: L{escript.FunctionSpace}
758     @param tags: boundray tags
759     @type tags: C{str}
760     @return: a mask which marks samples used by C{function_space} that are touching the
761     boundary tagged by any of the given tags.
762     @rtype: L{escript.Data} of rank 0
763     """
764     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
765     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
766     for t in tags: d.setTaggedValue(t,1.)
767     pde.setValue(y=d)
768     out=util.whereNonZero(pde.getRightHandSide())
769     if out.getFunctionSpace() == function_space:
770     return out
771     else:
772     return util.whereNonZero(util.interpolate(out,function_space))
773    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26