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

Contents of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1330 - (show annotations)
Mon Oct 22 04:54:49 2007 UTC (12 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 27987 byte(s)
more flexible stopping criterium for the PCG
1 #
2 # $Id$
3 #
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
17 """
18 Provides some tools related to PDEs.
19
20 Currently includes:
21 - Projector - to project a discontinuous
22 - Locator - to trace values in data objects at a certain location
23 - TimeIntegrationManager - to handel extraplotion in time
24 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
25
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 """
33
34 __author__="Lutz Gross, l.gross@uq.edu.au"
35 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
36 http://www.access.edu.au
37 Primary Business: Queensland, Australia"""
38 __license__="""Licensed under the Open Software License version 3.0
39 http://www.opensource.org/licenses/osl-3.0.php"""
40 __url__="http://www.iservo.edu.au/esys"
41 __version__="$Revision$"
42 __date__="$Date$"
43
44
45 import escript
46 import linearPDEs
47 import numarray
48 import util
49 import math
50
51 class TimeIntegrationManager:
52 """
53 a simple mechanism to manage time dependend values.
54
55 typical usage is::
56
57 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
65 @note: currently only p=1 is supported.
66 """
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 def getValue(self):
87 out=self.__v_mem[0]
88 if len(out)==1:
89 return out[0]
90 else:
91 return out
92
93 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
128
129 class Projector:
130 """
131 The Projector is a factory which projects a discontiuous function onto a
132 continuous function on the a given domain.
133 """
134 def __init__(self, domain, reduce = True, fast=True):
135 """
136 Create a continuous function space projector for a domain.
137
138 @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 """
142 self.__pde = linearPDEs.LinearPDE(domain)
143 if fast:
144 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
145 self.__pde.setSymmetryOn()
146 self.__pde.setReducedOrderTo(reduce)
147 self.__pde.setValue(D = 1.)
148 return
149
150 def __call__(self, input_data):
151 """
152 Projects input_data onto a continuous function
153
154 @param input_data: The input_data to be projected.
155 """
156 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
157 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
158 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 class NoPDE:
186 """
187 solves the following problem for u:
188
189 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 @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 """
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 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
221 @param Y: right hand side
222 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
223 @param q: location of constraints
224 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
225 @param r: value of solution at locations of constraints
226 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
227 """
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 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
255 @param Y: right hand side
256 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
257 @param q: location of constraints
258 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
259 @param r: value of solution at locations of constraints
260 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
261 """
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 class Locator:
299 """
300 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 """
307
308 def __init__(self,where,x=numarray.zeros((3,))):
309 """
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
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 """
319 if isinstance(where,escript.FunctionSpace):
320 self.__function_space=where
321 else:
322 self.__function_space=escript.ContinuousFunction(where)
323 if isinstance(x, list):
324 self.__id=[]
325 for p in x:
326 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
327 else:
328 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
329
330 def __str__(self):
331 """
332 Returns the coordinates of the Locator as a string.
333 """
334 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
349 def getX(self):
350 """
351 Returns the exact coordinates of the Locator.
352 """
353 return self(self.getFunctionSpace().getX())
354
355 def getFunctionSpace(self):
356 """
357 Returns the function space of the Locator.
358 """
359 return self.__function_space
360
361 def getId(self,item=None):
362 """
363 Returns the identifier of the location.
364 """
365 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
373
374 def __call__(self,data):
375 """
376 Returns the value of data at the Locator of a Data object otherwise
377 the object is returned.
378 """
379 return self.getValue(data)
380
381 def getValue(self,data):
382 """
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 if data.getFunctionSpace()==self.getFunctionSpace():
388 dat=data
389 else:
390 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 o=data.getValueOfGlobalDataPoint(*i)
397 if data.getRank()==0:
398 out.append(o[0])
399 else:
400 out.append(o)
401 return out
402 else:
403 out=data.getValueOfGlobalDataPoint(*id)
404 if data.getRank()==0:
405 return out[0]
406 else:
407 return out
408 else:
409 return data
410
411 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 class IterationHistory(object):
439 """
440 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 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 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
485
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 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
494 @param Aprod: returns the value Ax
495 @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 @param Msolve: solves Mx=r
497 @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 @param bilinearform: inner product C{<x,r>}
500 @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 @param iter_max: maximum number of iteration steps.
506 @type iter_max: C{int}
507 @return: the solution approximation and the corresponding residual
508 @rtype: C{tuple}
509 @warning: C{b} and C{x} are altered.
510 """
511 iter=0
512 if x==None:
513 x=0*b
514 else:
515 b += (-1)*Aprod(x)
516 r=b
517 rhat=Msolve(r)
518 d = rhat
519 rhat_dot_r = bilinearform(rhat, r)
520 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521
522 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
523 iter+=1
524 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 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
539
540 return x,r
541
542 class SaddlePointProblem(object):
543 """
544 This implements a solver for a saddlepoint problem
545
546 M{f(u,p)=0}
547 M{g(u)=0}
548
549 for u and p. The problem is solved with an inexact Uszawa scheme for p:
550
551 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
552 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
553
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 @param verbose: switches on the printing out some information
564 @type verbose: C{bool}
565 @note: this method may be overwritten by a particular saddle point problem
566 """
567 if not isinstance(verbose,bool):
568 raise TypeError("verbose needs to be of type bool.")
569 self.__verbose=verbose
570 self.relaxation=1.
571
572 def trace(self,text):
573 """
574 prints text if verbose has been set
575
576 @param text: a text message
577 @type text: C{str}
578 """
579 if self.__verbose: print "%s: %s"%(str(self),text)
580
581 def solve_f(self,u,p,tol=1.e-8):
582 """
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 @param tol: tolerance expected for du
594 @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 def solve_g(self,u,tol=1.e-8):
602 """
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 @param tol: tolerance expected for dp
612 @type tol: C{float}
613 @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 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
632 @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 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
731 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
771 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