/[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 1331 - (show annotations)
Tue Oct 23 00:42:15 2007 UTC (11 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 29931 byte(s)
ArithmeticTuple added. It is useful for certain applications of 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 ArithmeticTuple(object):
543 """
544 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
545
546 example of usage:
547
548 from esys.escript import Data
549 from numarray import array
550 a=Data(...)
551 b=array([1.,4.])
552 x=ArithmeticTuple(a,b)
553 y=5.*x
554
555 """
556 def __init__(self,*args):
557 """
558 initialize object with elements args.
559
560 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
561 """
562 self.__items=list(args)
563
564 def __len__(self):
565 """
566 number of items
567
568 @return: number of items
569 @rtype: C{int}
570 """
571 return len(self.__items)
572
573 def __getitem__(self,index):
574 """
575 get an item
576
577 @param index: item to be returned
578 @type index: C{int}
579 @return: item with index C{index}
580 """
581 return self.__items.__getitem__(index)
582
583 def __mul__(self,other):
584 """
585 scaling from the right
586
587 @param other: scaling factor
588 @type other: C{float}
589 @return: itemwise self*other
590 @rtype: L{ArithmeticTuple}
591 """
592 out=[]
593 for i in range(len(self)):
594 out.append(self[i]*other)
595 return ArithmeticTuple(*tuple(out))
596
597 def __rmul__(self,other):
598 """
599 scaling from the left
600
601 @param other: scaling factor
602 @type other: C{float}
603 @return: itemwise other*self
604 @rtype: L{ArithmeticTuple}
605 """
606 out=[]
607 for i in range(len(self)):
608 out.append(other*self[i])
609 return ArithmeticTuple(*tuple(out))
610
611 def __iadd__(self,other):
612 """
613 in-place add of other to self
614
615 @param other: increment
616 @type other: C{ArithmeticTuple}
617 """
618 if len(self) != len(other):
619 raise ValueError,"tuple length must match."
620 for i in range(len(self)):
621 self.__items[i]+=other[i]
622 return self
623
624 class SaddlePointProblem(object):
625 """
626 This implements a solver for a saddlepoint problem
627
628 M{f(u,p)=0}
629 M{g(u)=0}
630
631 for u and p. The problem is solved with an inexact Uszawa scheme for p:
632
633 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
634 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
635
636 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
637 A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
638 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
639 in fact the role of a preconditioner.
640 """
641 def __init__(self,verbose=False,*args):
642 """
643 initializes the problem
644
645 @param verbose: switches on the printing out some information
646 @type verbose: C{bool}
647 @note: this method may be overwritten by a particular saddle point problem
648 """
649 if not isinstance(verbose,bool):
650 raise TypeError("verbose needs to be of type bool.")
651 self.__verbose=verbose
652 self.relaxation=1.
653
654 def trace(self,text):
655 """
656 prints text if verbose has been set
657
658 @param text: a text message
659 @type text: C{str}
660 """
661 if self.__verbose: print "%s: %s"%(str(self),text)
662
663 def solve_f(self,u,p,tol=1.e-8):
664 """
665 solves
666
667 A_f du = f(u,p)
668
669 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
670
671 @param u: current approximation of u
672 @type u: L{escript.Data}
673 @param p: current approximation of p
674 @type p: L{escript.Data}
675 @param tol: tolerance expected for du
676 @type tol: C{float}
677 @return: increment du
678 @rtype: L{escript.Data}
679 @note: this method has to be overwritten by a particular saddle point problem
680 """
681 pass
682
683 def solve_g(self,u,tol=1.e-8):
684 """
685 solves
686
687 Q_g dp = g(u)
688
689 with Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p.
690
691 @param u: current approximation of u
692 @type u: L{escript.Data}
693 @param tol: tolerance expected for dp
694 @type tol: C{float}
695 @return: increment dp
696 @rtype: L{escript.Data}
697 @note: this method has to be overwritten by a particular saddle point problem
698 """
699 pass
700
701 def inner(self,p0,p1):
702 """
703 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
704 @return: inner product of p0 and p1
705 @rtype: C{float}
706 """
707 pass
708
709 subiter_max=3
710 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
711 """
712 runs the solver
713
714 @param u0: initial guess for C{u}
715 @type u0: L{esys.escript.Data}
716 @param p0: initial guess for C{p}
717 @type p0: L{esys.escript.Data}
718 @param tolerance: tolerance for relative error in C{u} and C{p}
719 @type tolerance: positive C{float}
720 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
721 @type tolerance_u: positive C{float}
722 @param iter_max: maximum number of iteration steps.
723 @type iter_max: C{int}
724 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
725 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
726 @type accepted_reduction: positive C{float} or C{None}
727 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
728 @type relaxation: C{float} or C{None}
729 """
730 tol=1.e-2
731 if tolerance_u==None: tolerance_u=tolerance
732 if not relaxation==None: self.relaxation=relaxation
733 if accepted_reduction ==None:
734 angle_limit=0.
735 elif accepted_reduction>=1.:
736 angle_limit=0.
737 else:
738 angle_limit=util.sqrt(1-accepted_reduction**2)
739 self.iter=0
740 u=u0
741 p=p0
742 #
743 # initialize things:
744 #
745 converged=False
746 #
747 # start loop:
748 #
749 # initial search direction is g
750 #
751 while not converged :
752 if self.iter>iter_max:
753 raise ArithmeticError("no convergence after %s steps."%self.iter)
754 f_new=self.solve_f(u,p,tol)
755 norm_f_new = util.Lsup(f_new)
756 u_new=u-f_new
757 g_new=self.solve_g(u_new,tol)
758 self.iter+=1
759 norm_g_new = util.sqrt(self.inner(g_new,g_new))
760 if norm_f_new==0. and norm_g_new==0.: return u, p
761 if self.iter>1 and not accepted_reduction==None:
762 #
763 # did we manage to reduce the norm of G? I
764 # if not we start a backtracking procedure
765 #
766 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
767 if norm_g_new > accepted_reduction * norm_g:
768 sub_iter=0
769 s=self.relaxation
770 d=g
771 g_last=g
772 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
773 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
774 dg= g_new-g_last
775 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
776 rad=self.inner(g_new,dg)/self.relaxation
777 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
778 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
779 if abs(rad) < norm_dg*norm_g_new * angle_limit:
780 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
781 break
782 r=self.relaxation
783 self.relaxation= - rad/norm_dg**2
784 s+=self.relaxation
785 #####
786 # a=g_new+self.relaxation*dg/r
787 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
788 #####
789 g_last=g_new
790 p+=self.relaxation*d
791 f_new=self.solve_f(u,p,tol)
792 u_new=u-f_new
793 g_new=self.solve_g(u_new,tol)
794 self.iter+=1
795 norm_f_new = util.Lsup(f_new)
796 norm_g_new = util.sqrt(self.inner(g_new,g_new))
797 # print " ",sub_iter," new g norm",norm_g_new
798 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
799 #
800 # can we expect reduction of g?
801 #
802 # u_last=u_new
803 sub_iter+=1
804 self.relaxation=s
805 #
806 # check for convergence:
807 #
808 norm_u_new = util.Lsup(u_new)
809 p_new=p+self.relaxation*g_new
810 norm_p_new = util.sqrt(self.inner(p_new,p_new))
811 self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
812
813 if self.iter>1:
814 dg2=g_new-g
815 df2=f_new-f
816 norm_dg2=util.sqrt(self.inner(dg2,dg2))
817 norm_df2=util.Lsup(df2)
818 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
819 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
820 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
821 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
822 converged=True
823 f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
824 self.trace("convergence after %s steps."%self.iter)
825 return u,p
826 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
827 # tol=1.e-2
828 # iter=0
829 # converged=False
830 # u=u0*1.
831 # p=p0*1.
832 # while not converged and iter<iter_max:
833 # du=self.solve_f(u,p,tol)
834 # u-=du
835 # norm_du=util.Lsup(du)
836 # norm_u=util.Lsup(u)
837 #
838 # dp=self.relaxation*self.solve_g(u,tol)
839 # p+=dp
840 # norm_dp=util.sqrt(self.inner(dp,dp))
841 # norm_p=util.sqrt(self.inner(p,p))
842 # print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
843 # iter+=1
844 #
845 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
846 # if converged:
847 # print "convergence after %s steps."%iter
848 # else:
849 # raise ArithmeticError("no convergence after %s steps."%iter)
850 #
851 # return u,p
852
853 def MaskFromBoundaryTag(function_space,*tags):
854 """
855 create a mask on the given function space which one for samples
856 that touch the boundary tagged by tags.
857
858 usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
859
860 @param function_space: a given function space
861 @type function_space: L{escript.FunctionSpace}
862 @param tags: boundray tags
863 @type tags: C{str}
864 @return: a mask which marks samples used by C{function_space} that are touching the
865 boundary tagged by any of the given tags.
866 @rtype: L{escript.Data} of rank 0
867 """
868 pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
869 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
870 for t in tags: d.setTaggedValue(t,1.)
871 pde.setValue(y=d)
872 out=util.whereNonZero(pde.getRightHandSide())
873 if out.getFunctionSpace() == function_space:
874 return out
875 else:
876 return util.whereNonZero(util.interpolate(out,function_space))
877

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26