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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (11 years, 11 months ago) by ksteube
Original Path: trunk/escript/py_src/pdetools.py
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 #
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 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 class SaddlePointProblem(object):
521 """
522 This implements a solver for a saddlepoint problem
523
524 M{f(u,p)=0}
525 M{g(u)=0}
526
527 for u and p. The problem is solved with an inexact Uszawa scheme for p:
528
529 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
530 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
531
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 @param verbose: switches on the printing out some information
542 @type verbose: C{bool}
543 @note: this method may be overwritten by a particular saddle point problem
544 """
545 if not isinstance(verbose,bool):
546 raise TypeError("verbose needs to be of type bool.")
547 self.__verbose=verbose
548 self.relaxation=1.
549
550 def trace(self,text):
551 """
552 prints text if verbose has been set
553
554 @param text: a text message
555 @type text: C{str}
556 """
557 if self.__verbose: print "%s: %s"%(str(self),text)
558
559 def solve_f(self,u,p,tol=1.e-8):
560 """
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 @param tol: tolerance expected for du
572 @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 def solve_g(self,u,tol=1.e-8):
580 """
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 @param tol: tolerance expected for dp
590 @type tol: C{float}
591 @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 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
610 @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 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
709 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
749 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