/[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 1489 - (show annotations)
Mon Apr 14 04:29:30 2008 UTC (11 years, 8 months ago) by artak
File MIME type: text/x-python
File size: 47751 byte(s)
Initial version of TFQMR method
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 ##### Added by Artak
52 # from Numeric import zeros,Int,Float64
53 ###################################
54
55
56 class TimeIntegrationManager:
57 """
58 a simple mechanism to manage time dependend values.
59
60 typical usage is::
61
62 dt=0.1 # time increment
63 tm=TimeIntegrationManager(inital_value,p=1)
64 while t<1.
65 v_guess=tm.extrapolate(dt) # extrapolate to t+dt
66 v=...
67 tm.checkin(dt,v)
68 t+=dt
69
70 @note: currently only p=1 is supported.
71 """
72 def __init__(self,*inital_values,**kwargs):
73 """
74 sets up the value manager where inital_value is the initial value and p is order used for extrapolation
75 """
76 if kwargs.has_key("p"):
77 self.__p=kwargs["p"]
78 else:
79 self.__p=1
80 if kwargs.has_key("time"):
81 self.__t=kwargs["time"]
82 else:
83 self.__t=0.
84 self.__v_mem=[inital_values]
85 self.__order=0
86 self.__dt_mem=[]
87 self.__num_val=len(inital_values)
88
89 def getTime(self):
90 return self.__t
91 def getValue(self):
92 out=self.__v_mem[0]
93 if len(out)==1:
94 return out[0]
95 else:
96 return out
97
98 def checkin(self,dt,*values):
99 """
100 adds new values to the manager. the p+1 last value get lost
101 """
102 o=min(self.__order+1,self.__p)
103 self.__order=min(self.__order+1,self.__p)
104 v_mem_new=[values]
105 dt_mem_new=[dt]
106 for i in range(o-1):
107 v_mem_new.append(self.__v_mem[i])
108 dt_mem_new.append(self.__dt_mem[i])
109 v_mem_new.append(self.__v_mem[o-1])
110 self.__order=o
111 self.__v_mem=v_mem_new
112 self.__dt_mem=dt_mem_new
113 self.__t+=dt
114
115 def extrapolate(self,dt):
116 """
117 extrapolates to dt forward in time.
118 """
119 if self.__order==0:
120 out=self.__v_mem[0]
121 else:
122 out=[]
123 for i in range(self.__num_val):
124 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
125
126 if len(out)==0:
127 return None
128 elif len(out)==1:
129 return out[0]
130 else:
131 return out
132
133
134 class Projector:
135 """
136 The Projector is a factory which projects a discontiuous function onto a
137 continuous function on the a given domain.
138 """
139 def __init__(self, domain, reduce = True, fast=True):
140 """
141 Create a continuous function space projector for a domain.
142
143 @param domain: Domain of the projection.
144 @param reduce: Flag to reduce projection order (default is True)
145 @param fast: Flag to use a fast method based on matrix lumping (default is true)
146 """
147 self.__pde = linearPDEs.LinearPDE(domain)
148 if fast:
149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
150 self.__pde.setSymmetryOn()
151 self.__pde.setReducedOrderTo(reduce)
152 self.__pde.setValue(D = 1.)
153 return
154
155 def __call__(self, input_data):
156 """
157 Projects input_data onto a continuous function
158
159 @param input_data: The input_data to be projected.
160 """
161 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
162 self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
163 if input_data.getRank()==0:
164 self.__pde.setValue(Y = input_data)
165 out=self.__pde.getSolution()
166 elif input_data.getRank()==1:
167 for i0 in range(input_data.getShape()[0]):
168 self.__pde.setValue(Y = input_data[i0])
169 out[i0]=self.__pde.getSolution()
170 elif input_data.getRank()==2:
171 for i0 in range(input_data.getShape()[0]):
172 for i1 in range(input_data.getShape()[1]):
173 self.__pde.setValue(Y = input_data[i0,i1])
174 out[i0,i1]=self.__pde.getSolution()
175 elif input_data.getRank()==3:
176 for i0 in range(input_data.getShape()[0]):
177 for i1 in range(input_data.getShape()[1]):
178 for i2 in range(input_data.getShape()[2]):
179 self.__pde.setValue(Y = input_data[i0,i1,i2])
180 out[i0,i1,i2]=self.__pde.getSolution()
181 else:
182 for i0 in range(input_data.getShape()[0]):
183 for i1 in range(input_data.getShape()[1]):
184 for i2 in range(input_data.getShape()[2]):
185 for i3 in range(input_data.getShape()[3]):
186 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
187 out[i0,i1,i2,i3]=self.__pde.getSolution()
188 return out
189
190 class NoPDE:
191 """
192 solves the following problem for u:
193
194 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
195
196 with constraint
197
198 M{u[j]=r[j]} where M{q[j]>0}
199
200 where D, Y, r and q are given functions of rank 1.
201
202 In the case of scalars this takes the form
203
204 M{D*u=Y}
205
206 with constraint
207
208 M{u=r} where M{q>0}
209
210 where D, Y, r and q are given scalar functions.
211
212 The constraint is overwriting any other condition.
213
214 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
215 that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
216 thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
217 """
218 def __init__(self,domain,D=None,Y=None,q=None,r=None):
219 """
220 initialize the problem
221
222 @param domain: domain of the PDE.
223 @type domain: L{Domain}
224 @param D: coefficient of the solution.
225 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
226 @param Y: right hand side
227 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
228 @param q: location of constraints
229 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
230 @param r: value of solution at locations of constraints
231 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
232 """
233 self.__domain=domain
234 self.__D=D
235 self.__Y=Y
236 self.__q=q
237 self.__r=r
238 self.__u=None
239 self.__function_space=escript.Solution(self.__domain)
240 def setReducedOn(self):
241 """
242 sets the L{FunctionSpace} of the solution to L{ReducedSolution}
243 """
244 self.__function_space=escript.ReducedSolution(self.__domain)
245 self.__u=None
246
247 def setReducedOff(self):
248 """
249 sets the L{FunctionSpace} of the solution to L{Solution}
250 """
251 self.__function_space=escript.Solution(self.__domain)
252 self.__u=None
253
254 def setValue(self,D=None,Y=None,q=None,r=None):
255 """
256 assigns values to the parameters.
257
258 @param D: coefficient of the solution.
259 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
260 @param Y: right hand side
261 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
262 @param q: location of constraints
263 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
264 @param r: value of solution at locations of constraints
265 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
266 """
267 if not D==None:
268 self.__D=D
269 self.__u=None
270 if not Y==None:
271 self.__Y=Y
272 self.__u=None
273 if not q==None:
274 self.__q=q
275 self.__u=None
276 if not r==None:
277 self.__r=r
278 self.__u=None
279
280 def getSolution(self):
281 """
282 returns the solution
283
284 @return: the solution of the problem
285 @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
286 """
287 if self.__u==None:
288 if self.__D==None:
289 raise ValueError,"coefficient D is undefined"
290 D=escript.Data(self.__D,self.__function_space)
291 if D.getRank()>1:
292 raise ValueError,"coefficient D must have rank 0 or 1"
293 if self.__Y==None:
294 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
295 else:
296 self.__u=util.quotient(self.__Y,D)
297 if not self.__q==None:
298 q=util.wherePositive(escript.Data(self.__q,self.__function_space))
299 self.__u*=(1.-q)
300 if not self.__r==None: self.__u+=q*self.__r
301 return self.__u
302
303 class Locator:
304 """
305 Locator provides access to the values of data objects at a given
306 spatial coordinate x.
307
308 In fact, a Locator object finds the sample in the set of samples of a
309 given function space or domain where which is closest to the given
310 point x.
311 """
312
313 def __init__(self,where,x=numarray.zeros((3,))):
314 """
315 Initializes a Locator to access values in Data objects on the Doamin
316 or FunctionSpace where for the sample point which
317 closest to the given point x.
318
319 @param where: function space
320 @type where: L{escript.FunctionSpace}
321 @param x: coefficient of the solution.
322 @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
323 """
324 if isinstance(where,escript.FunctionSpace):
325 self.__function_space=where
326 else:
327 self.__function_space=escript.ContinuousFunction(where)
328 if isinstance(x, list):
329 self.__id=[]
330 for p in x:
331 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
332 else:
333 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
334
335 def __str__(self):
336 """
337 Returns the coordinates of the Locator as a string.
338 """
339 x=self.getX()
340 if instance(x,list):
341 out="["
342 first=True
343 for xx in x:
344 if not first:
345 out+=","
346 else:
347 first=False
348 out+=str(xx)
349 out+="]>"
350 else:
351 out=str(x)
352 return out
353
354 def getX(self):
355 """
356 Returns the exact coordinates of the Locator.
357 """
358 return self(self.getFunctionSpace().getX())
359
360 def getFunctionSpace(self):
361 """
362 Returns the function space of the Locator.
363 """
364 return self.__function_space
365
366 def getId(self,item=None):
367 """
368 Returns the identifier of the location.
369 """
370 if item == None:
371 return self.__id
372 else:
373 if isinstance(self.__id,list):
374 return self.__id[item]
375 else:
376 return self.__id
377
378
379 def __call__(self,data):
380 """
381 Returns the value of data at the Locator of a Data object otherwise
382 the object is returned.
383 """
384 return self.getValue(data)
385
386 def getValue(self,data):
387 """
388 Returns the value of data at the Locator if data is a Data object
389 otherwise the object is returned.
390 """
391 if isinstance(data,escript.Data):
392 if data.getFunctionSpace()==self.getFunctionSpace():
393 dat=data
394 else:
395 dat=data.interpolate(self.getFunctionSpace())
396 id=self.getId()
397 r=data.getRank()
398 if isinstance(id,list):
399 out=[]
400 for i in id:
401 o=data.getValueOfGlobalDataPoint(*i)
402 if data.getRank()==0:
403 out.append(o[0])
404 else:
405 out.append(o)
406 return out
407 else:
408 out=data.getValueOfGlobalDataPoint(*id)
409 if data.getRank()==0:
410 return out[0]
411 else:
412 return out
413 else:
414 return data
415
416 class SolverSchemeException(Exception):
417 """
418 exceptions thrown by solvers
419 """
420 pass
421
422 class IndefinitePreconditioner(SolverSchemeException):
423 """
424 the preconditioner is not positive definite.
425 """
426 pass
427 class MaxIterReached(SolverSchemeException):
428 """
429 maxium number of iteration steps is reached.
430 """
431 pass
432 class IterationBreakDown(SolverSchemeException):
433 """
434 iteration scheme econouters an incurable breakdown.
435 """
436 pass
437 class NegativeNorm(SolverSchemeException):
438 """
439 a norm calculation returns a negative norm.
440 """
441 pass
442
443 class IterationHistory(object):
444 """
445 The IterationHistory class is used to define a stopping criterium. It keeps track of the
446 residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
447 a given tolerance.
448 """
449 def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
450 """
451 Initialization
452
453 @param tolerance: tolerance
454 @type tolerance: positive C{float}
455 @param verbose: switches on the printing out some information
456 @type verbose: C{bool}
457 """
458 if not tolerance>0.:
459 raise ValueError,"tolerance needs to be positive."
460 self.tolerance=tolerance
461 self.verbose=verbose
462 self.history=[]
463 def stoppingcriterium(self,norm_r,r,x):
464 """
465 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.
466
467
468 @param norm_r: current residual norm
469 @type norm_r: non-negative C{float}
470 @param r: current residual (not used)
471 @param x: current solution approximation (not used)
472 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
473 @rtype: C{bool}
474
475 """
476 self.history.append(norm_r)
477 if self.verbose: print "iter: #s: inner(rhat,r) = #e"#(len(self.history)-1, self.history[-1])
478 return self.history[-1]<=self.tolerance * self.history[0]
479
480 def stoppingcriterium2(self,norm_r,norm_b):
481 """
482 returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483
484
485 @param norm_r: current residual norm
486 @type norm_r: non-negative C{float}
487 @param norm_b: norm of right hand side
488 @type norm_b: non-negative C{float}
489 @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
490 @rtype: C{bool}
491
492 """
493 self.history.append(norm_r)
494 if self.verbose: print "iter: #s: norm(r) = #e"#(len(self.history)-1, self.history[-1])
495 return self.history[-1]<=self.tolerance * norm_b
496
497 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
498 """
499 Solver for
500
501 M{Ax=b}
502
503 with a symmetric and positive definite operator A (more details required!).
504 It uses the conjugate gradient method with preconditioner M providing an approximation of A.
505
506 The iteration is terminated if the C{stoppingcriterium} function return C{True}.
507
508 For details on the preconditioned conjugate gradient method see the book:
509
510 Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
511 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
512 C. Romine, and H. van der Vorst.
513
514 @param b: the right hand side of the liner system. C{b} is altered.
515 @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
516 @param Aprod: returns the value Ax
517 @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}.
518 @param Msolve: solves Mx=r
519 @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
520 type like argument C{x}.
521 @param bilinearform: inner product C{<x,r>}
522 @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}.
523 @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.
524 @type stoppingcriterium: function that returns C{True} or C{False}
525 @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
526 @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
527 @param iter_max: maximum number of iteration steps.
528 @type iter_max: C{int}
529 @return: the solution approximation and the corresponding residual
530 @rtype: C{tuple}
531 @warning: C{b} and C{x} are altered.
532 """
533 iter=0
534 if x==None:
535 x=0*b
536 else:
537 b += (-1)*Aprod(x)
538 r=b
539 rhat=Msolve(r)
540 d = rhat
541 rhat_dot_r = bilinearform(rhat, r)
542 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
543
544 while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
545 iter+=1
546 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
547
548 q=Aprod(d)
549 alpha = rhat_dot_r / bilinearform(d, q)
550 x += alpha * d
551 r += (-alpha) * q
552
553 rhat=Msolve(r)
554 rhat_dot_r_new = bilinearform(rhat, r)
555 beta = rhat_dot_r_new / rhat_dot_r
556 rhat+=beta * d
557 d=rhat
558
559 rhat_dot_r = rhat_dot_r_new
560 if rhat_dot_r<0: raise NegativeNorm,"negative norm."
561
562 return x,r
563
564
565 ############################
566 # Added by Artak
567 #################################3
568
569 #Apply a sequence of k Givens rotations, used within gmres codes
570 # vrot=givapp(c, s, vin, k)
571 def givapp(c,s,vin):
572 vrot=vin # warning: vin is altered!!!!
573 if isinstance(c,float):
574 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
575 else:
576 for i in range(len(c)):
577 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
578 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
579 vrot[i:i+2]=w1,w2
580 return vrot
581
582 def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
583 m=iter_restart
584 iter=0
585 while True:
586 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
587 x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588 iter+=iter_restart
589 if stopped: break
590 return x
591
592 def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
593 iter=0
594 r=Msolve(b)
595 r_dot_r = bilinearform(r, r)
596 if r_dot_r<0: raise NegativeNorm,"negative norm."
597 norm_b=math.sqrt(r_dot_r)
598
599 if x==None:
600 x=0
601 else:
602 r=Msolve(b-Aprod(x))
603 r_dot_r = bilinearform(r, r)
604 if r_dot_r<0: raise NegativeNorm,"negative norm."
605
606 h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
607 c=numarray.zeros(iter_restart,numarray.Float64)
608 s=numarray.zeros(iter_restart,numarray.Float64)
609 g=numarray.zeros(iter_restart,numarray.Float64)
610 v=[]
611
612 rho=math.sqrt(r_dot_r)
613
614 v.append(r/rho)
615 g[0]=rho
616
617 while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
618
619 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
620
621
622 p=Msolve(Aprod(v[iter]))
623
624 v.append(p)
625
626 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
627
628 # Modified Gram-Schmidt
629 for j in range(iter+1):
630 h[j][iter]=bilinearform(v[j],v[iter+1])
631 v[iter+1]+=(-1.)*h[j][iter]*v[j]
632
633 h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
634 v_norm2=h[iter+1][iter]
635
636
637 # Reorthogonalize if needed
638 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
639 for j in range(iter+1):
640 hr=bilinearform(v[j],v[iter+1])
641 h[j][iter]=h[j][iter]+hr #vhat
642 v[iter+1] +=(-1.)*hr*v[j]
643
644 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
645 h[iter+1][iter]=v_norm2
646
647 # watch out for happy breakdown
648 if v_norm2 != 0:
649 v[iter+1]=v[iter+1]/h[iter+1][iter]
650
651 # Form and store the information for the new Givens rotation
652 if iter > 0 :
653 hhat=[]
654 for i in range(iter+1) : hhat.append(h[i][iter])
655 hhat=givapp(c[0:iter],s[0:iter],hhat);
656 for i in range(iter+1) : h[i][iter]=hhat[i]
657
658 mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
659 if mu!=0 :
660 c[iter]=h[iter][iter]/mu
661 s[iter]=-h[iter+1][iter]/mu
662 h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
663 h[iter+1][iter]=0.0
664 g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
665
666 # Update the residual norm
667 rho=abs(g[iter+1])
668 iter+=1
669
670 # At this point either iter > iter_max or rho < tol.
671 # It's time to compute x and leave.
672
673 if iter > 0 :
674 y=numarray.zeros(iter,numarray.Float64)
675 y[iter-1] = g[iter-1] / h[iter-1][iter-1]
676 if iter > 1 :
677 i=iter-2
678 while i>=0 :
679 y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
680 i=i-1
681 xhat=v[iter-1]*y[iter-1]
682 for i in range(iter-1):
683 xhat += v[i]*y[i]
684 else : xhat=v[0]
685
686 x += xhat
687 if iter<iter_restart-1:
688 stopped=True
689 else:
690 stopped=False
691
692 return x,stopped
693
694 def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
695
696 #
697 # minres solves the system of linear equations Ax = b
698 # where A is a symmetric matrix (possibly indefinite or singular)
699 # and b is a given vector.
700 #
701 # "A" may be a dense or sparse matrix (preferably sparse!)
702 # or the name of a function such that
703 # y = A(x)
704 # returns the product y = Ax for any given vector x.
705 #
706 # "M" defines a positive-definite preconditioner M = C C'.
707 # "M" may be a dense or sparse matrix (preferably sparse!)
708 # or the name of a function such that
709 # solves the system My = x for any given vector x.
710 #
711 #
712
713 #------------------------------------------------------------------
714 # Set up y and v for the first Lanczos vector v1.
715 # y = beta1 P' v1, where P = C**(-1).
716 # v is really P' v1.
717 #------------------------------------------------------------------
718 if x==None:
719 x=0*b
720 else:
721 b += (-1)*Aprod(x)
722
723 r1 = b
724 y = Msolve(b)
725 beta1 = bilinearform(b,y)
726
727 if beta1< 0: raise NegativeNorm,"negative norm."
728
729 # If b = 0 exactly, stop with x = 0.
730 if beta1==0: return x*0.
731
732 if beta1> 0:
733 beta1 = math.sqrt(beta1)
734
735 #------------------------------------------------------------------
736 # Initialize quantities.
737 # ------------------------------------------------------------------
738 iter = 0
739 Anorm = 0
740 ynorm = 0
741 oldb = 0
742 beta = beta1
743 dbar = 0
744 epsln = 0
745 phibar = beta1
746 rhs1 = beta1
747 rhs2 = 0
748 rnorm = phibar
749 tnorm2 = 0
750 ynorm2 = 0
751 cs = -1
752 sn = 0
753 w = b*0.
754 w2 = b*0.
755 r2 = r1
756 eps = 0.0001
757
758 #---------------------------------------------------------------------
759 # Main iteration loop.
760 # --------------------------------------------------------------------
761 while not stoppingcriterium(rnorm,Anorm*ynorm): # checks ||r|| < (||A|| ||x||) * TOL
762
763 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
764 iter = iter + 1
765
766 #-----------------------------------------------------------------
767 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
768 # The general iteration is similar to the case k = 1 with v0 = 0:
769 #
770 # p1 = Operator * v1 - beta1 * v0,
771 # alpha1 = v1'p1,
772 # q2 = p2 - alpha1 * v1,
773 # beta2^2 = q2'q2,
774 # v2 = (1/beta2) q2.
775 #
776 # Again, y = betak P vk, where P = C**(-1).
777 #-----------------------------------------------------------------
778 s = 1/beta # Normalize previous vector (in y).
779 v = s*y # v = vk if P = I
780
781 y = Aprod(v)
782
783 if iter >= 2:
784 y = y - (beta/oldb)*r1
785
786 alfa = bilinearform(v,y) # alphak
787 y = (- alfa/beta)*r2 + y
788 r1 = r2
789 r2 = y
790 y = Msolve(r2)
791 oldb = beta # oldb = betak
792 beta = bilinearform(r2,y) # beta = betak+1^2
793 if beta < 0: raise NegativeNorm,"negative norm."
794
795 beta = math.sqrt( beta )
796 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
797
798 if iter==1: # Initialize a few things.
799 gmax = abs( alfa ) # alpha1
800 gmin = gmax # alpha1
801
802 # Apply previous rotation Qk-1 to get
803 # [deltak epslnk+1] = [cs sn][dbark 0 ]
804 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
805
806 oldeps = epsln
807 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
808 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
809 epsln = sn * beta # epsln2 = 0 epslnk+1
810 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
811
812 # Compute the next plane rotation Qk
813
814 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
815 gamma = max(gamma,eps)
816 cs = gbar / gamma # ck
817 sn = beta / gamma # sk
818 phi = cs * phibar # phik
819 phibar = sn * phibar # phibark+1
820
821 # Update x.
822
823 denom = 1/gamma
824 w1 = w2
825 w2 = w
826 w = (v - oldeps*w1 - delta*w2) * denom
827 x = x + phi*w
828
829 # Go round again.
830
831 gmax = max(gmax,gamma)
832 gmin = min(gmin,gamma)
833 z = rhs1 / gamma
834 ynorm2 = z*z + ynorm2
835 rhs1 = rhs2 - delta*z
836 rhs2 = - epsln*z
837
838 # Estimate various norms and test for convergence.
839
840 Anorm = math.sqrt( tnorm2 )
841 ynorm = math.sqrt( ynorm2 )
842
843 rnorm = phibar
844
845 return x
846
847
848 def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
849
850 # TFQMR solver for linear systems
851 #
852 #
853 # initialization
854 #
855 errtol = math.sqrt(bilinearform(b,b))
856 norm_b=errtol
857 kmax = iter_max
858 error = []
859
860 if math.sqrt(bilinearform(x,x)) != 0.0:
861 r = b - Aprod(x)
862 else:
863 r = b
864
865 r=Msolve(r)
866
867 u1=0
868 u2=0
869 y1=0
870 y2=0
871
872 w = r
873 y1 = r
874 iter = 0
875 d = 0
876
877 v = Msolve(Aprod(y1))
878 u1 = v
879
880 theta = 0.0;
881 eta = 0.0;
882 tau = math.sqrt(bilinearform(r,r))
883 error = [ error, tau ]
884 rho = tau * tau
885 m=1
886 #
887 # TFQMR iteration
888 #
889 # while ( iter < kmax-1 ):
890
891 while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b):
892 if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
893
894 sigma = bilinearform(r,v)
895
896 if ( sigma == 0.0 ):
897 raise 'TFQMR breakdown, sigma=0'
898
899
900 alpha = rho / sigma
901
902 for j in range(2):
903 #
904 # Compute y2 and u2 only if you have to
905 #
906 if ( j == 1 ):
907 y2 = y1 - alpha * v
908 u2 = Msolve(Aprod(y2))
909
910 m = 2 * (iter+1) - 2 + (j+1)
911 if j==0:
912 w = w - alpha * u1
913 d = y1 + ( theta * theta * eta / alpha ) * d
914 if j==1:
915 w = w - alpha * u2
916 d = y2 + ( theta * theta * eta / alpha ) * d
917
918 theta = math.sqrt(bilinearform(w,w))/ tau
919 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
920 tau = tau * theta * c
921 eta = c * c * alpha
922 x = x + eta * d
923 #
924 # Try to terminate the iteration at each pass through the loop
925 #
926 # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
927 # error = [ error, tau ]
928 # total_iters = iter
929 # break
930
931
932 if ( rho == 0.0 ):
933 raise 'TFQMR breakdown, rho=0'
934
935
936 rhon = bilinearform(r,w)
937 beta = rhon / rho;
938 rho = rhon;
939 y1 = w + beta * y2;
940 u1 = Msolve(Aprod(y1))
941 v = u1 + beta * ( u2 + beta * v )
942 error = [ error, tau ]
943 total_iters = iter
944
945 iter = iter + 1
946
947 print iter
948 return x
949
950
951 #############################################
952
953 class ArithmeticTuple(object):
954 """
955 tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
956
957 example of usage:
958
959 from esys.escript import Data
960 from numarray import array
961 a=Data(...)
962 b=array([1.,4.])
963 x=ArithmeticTuple(a,b)
964 y=5.*x
965
966 """
967 def __init__(self,*args):
968 """
969 initialize object with elements args.
970
971 @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
972 """
973 self.__items=list(args)
974
975 def __len__(self):
976 """
977 number of items
978
979 @return: number of items
980 @rtype: C{int}
981 """
982 return len(self.__items)
983
984 def __getitem__(self,index):
985 """
986 get an item
987
988 @param index: item to be returned
989 @type index: C{int}
990 @return: item with index C{index}
991 """
992 return self.__items.__getitem__(index)
993
994 def __mul__(self,other):
995 """
996 scaling from the right
997
998 @param other: scaling factor
999 @type other: C{float}
1000 @return: itemwise self*other
1001 @rtype: L{ArithmeticTuple}
1002 """
1003 out=[]
1004 for i in range(len(self)):
1005 out.append(self[i]*other)
1006 return ArithmeticTuple(*tuple(out))
1007
1008 def __rmul__(self,other):
1009 """
1010 scaling from the left
1011
1012 @param other: scaling factor
1013 @type other: C{float}
1014 @return: itemwise other*self
1015 @rtype: L{ArithmeticTuple}
1016 """
1017 out=[]
1018 for i in range(len(self)):
1019 out.append(other*self[i])
1020 return ArithmeticTuple(*tuple(out))
1021
1022 #########################
1023 # Added by Artak
1024 #########################
1025 def __div__(self,other):
1026 """
1027 dividing from the right
1028
1029 @param other: scaling factor
1030 @type other: C{float}
1031 @return: itemwise self/other
1032 @rtype: L{ArithmeticTuple}
1033 """
1034 out=[]
1035 for i in range(len(self)):
1036 out.append(self[i]/other)
1037 return ArithmeticTuple(*tuple(out))
1038
1039 def __rdiv__(self,other):
1040 """
1041 dividing from the left
1042
1043 @param other: scaling factor
1044 @type other: C{float}
1045 @return: itemwise other/self
1046 @rtype: L{ArithmeticTuple}
1047 """
1048 out=[]
1049 for i in range(len(self)):
1050 out.append(other/self[i])
1051 return ArithmeticTuple(*tuple(out))
1052
1053 ##########################################33
1054
1055 def __iadd__(self,other):
1056 """
1057 in-place add of other to self
1058
1059 @param other: increment
1060 @type other: C{ArithmeticTuple}
1061 """
1062 if len(self) != len(other):
1063 raise ValueError,"tuple length must match."
1064 for i in range(len(self)):
1065 self.__items[i]+=other[i]
1066 return self
1067
1068 class HomogeneousSaddlePointProblem(object):
1069 """
1070 This provides a framwork for solving homogeneous saddle point problem of the form
1071
1072 Av+B^*p=f
1073 Bv =0
1074
1075 for the unknowns v and p and given operators A and B and given right hand side f.
1076 B^* is the adjoint operator of B is the given inner product.
1077
1078 """
1079 def __init__(self,**kwargs):
1080 self.setTolerance()
1081 self.setToleranceReductionFactor()
1082
1083 def initialize(self):
1084 """
1085 initialize the problem (overwrite)
1086 """
1087 pass
1088 def B(self,v):
1089 """
1090 returns Bv (overwrite)
1091 @rtype: equal to the type of p
1092
1093 @note: boundary conditions on p should be zero!
1094 """
1095 pass
1096
1097 def inner(self,p0,p1):
1098 """
1099 returns inner product of two element p0 and p1 (overwrite)
1100
1101 @type p0: equal to the type of p
1102 @type p1: equal to the type of p
1103 @rtype: C{float}
1104
1105 @rtype: equal to the type of p
1106 """
1107 pass
1108
1109 def solve_A(self,u,p):
1110 """
1111 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1112
1113 @rtype: equal to the type of v
1114 @note: boundary conditions on v should be zero!
1115 """
1116 pass
1117
1118 def solve_prec(self,p):
1119 """
1120 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1121
1122 @rtype: equal to the type of p
1123 """
1124 pass
1125
1126 def stoppingcriterium(self,Bv,v,p):
1127 """
1128 returns a True if iteration is terminated. (overwrite)
1129
1130 @rtype: C{bool}
1131 """
1132 pass
1133
1134 def __inner(self,p,r):
1135 return self.inner(p,r[1])
1136
1137 def __inner_p(self,p1,p2):
1138 return self.inner(p1,p2)
1139
1140 def __stoppingcriterium(self,norm_r,r,p):
1141 return self.stoppingcriterium(r[1],r[0],p)
1142
1143 def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1144 return self.stoppingcriterium_GMRES(norm_r,norm_b)
1145
1146 def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1147 return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1148
1149
1150 def setTolerance(self,tolerance=1.e-8):
1151 self.__tol=tolerance
1152 def getTolerance(self):
1153 return self.__tol
1154 def setToleranceReductionFactor(self,reduction=0.01):
1155 self.__reduction=reduction
1156 def getSubProblemTolerance(self):
1157 return self.__reduction*self.getTolerance()
1158
1159 def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=10):
1160 """
1161 solves the saddle point problem using initial guesses v and p.
1162
1163 @param max_iter: maximum number of iteration steps.
1164 """
1165 self.verbose=verbose
1166 self.show_details=show_details and self.verbose
1167
1168 # assume p is known: then v=A^-1(f-B^*p)
1169 # which leads to BA^-1B^*p = BA^-1f
1170
1171 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
1172
1173
1174 self.__z=v+self.solve_A(v,p*0)
1175
1176 Bz=self.B(self.__z)
1177 #
1178 # solve BA^-1B^*p = Bz
1179 #
1180 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1181 #
1182 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1183 # A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1184 #
1185 self.iter=0
1186 if solver=='GMRES':
1187 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1188 p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1189 # solve Au=f-B^*p
1190 # A(u-v)=f-B^*p-Av
1191 # u=v+(u-v)
1192 u=v+self.solve_A(v,p)
1193
1194 if solver=='TFQMR':
1195 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1196 p=TFQMR(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1197 # solve Au=f-B^*p
1198 # A(u-v)=f-B^*p-Av
1199 # u=v+(u-v)
1200 u=v+self.solve_A(v,p)
1201
1202 if solver=='MINRES':
1203 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1204 p=MINRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1205 # solve Au=f-B^*p
1206 # A(u-v)=f-B^*p-Av
1207 # u=v+(u-v)
1208 u=v+self.solve_A(v,p)
1209
1210 if solver=='PCG':
1211 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1212 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1213 u=r[0]
1214
1215 print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1216
1217 return u,p
1218
1219 def __Msolve(self,r):
1220 return self.solve_prec(r[1])
1221
1222 def __Msolve_GMRES(self,r):
1223 return self.solve_prec(r)
1224
1225
1226 def __Aprod(self,p):
1227 # return BA^-1B*p
1228 #solve Av =-B^*p as Av =f-Az-B^*p
1229 v=self.solve_A(self.__z,-p)
1230 return ArithmeticTuple(v, self.B(v))
1231
1232 def __Aprod_GMRES(self,p):
1233 # return BA^-1B*p
1234 #solve Av =-B^*p as Av =f-Az-B^*p
1235 v=self.solve_A(self.__z,-p)
1236 return self.B(v)
1237
1238 class SaddlePointProblem(object):
1239 """
1240 This implements a solver for a saddlepoint problem
1241
1242 M{f(u,p)=0}
1243 M{g(u)=0}
1244
1245 for u and p. The problem is solved with an inexact Uszawa scheme for p:
1246
1247 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1248 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
1249
1250 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
1251 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'
1252 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1253 in fact the role of a preconditioner.
1254 """
1255 def __init__(self,verbose=False,*args):
1256 """
1257 initializes the problem
1258
1259 @param verbose: switches on the printing out some information
1260 @type verbose: C{bool}
1261 @note: this method may be overwritten by a particular saddle point problem
1262 """
1263 if not isinstance(verbose,bool):
1264 raise TypeError("verbose needs to be of type bool.")
1265 self.__verbose=verbose
1266 self.relaxation=1.
1267
1268 def trace(self,text):
1269 """
1270 prints text if verbose has been set
1271
1272 @param text: a text message
1273 @type text: C{str}
1274 """
1275 if self.__verbose: print "#s: #s"%(str(self),text)
1276
1277 def solve_f(self,u,p,tol=1.e-8):
1278 """
1279 solves
1280
1281 A_f du = f(u,p)
1282
1283 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
1284
1285 @param u: current approximation of u
1286 @type u: L{escript.Data}
1287 @param p: current approximation of p
1288 @type p: L{escript.Data}
1289 @param tol: tolerance expected for du
1290 @type tol: C{float}
1291 @return: increment du
1292 @rtype: L{escript.Data}
1293 @note: this method has to be overwritten by a particular saddle point problem
1294 """
1295 pass
1296
1297 def solve_g(self,u,tol=1.e-8):
1298 """
1299 solves
1300
1301 Q_g dp = g(u)
1302
1303 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.
1304
1305 @param u: current approximation of u
1306 @type u: L{escript.Data}
1307 @param tol: tolerance expected for dp
1308 @type tol: C{float}
1309 @return: increment dp
1310 @rtype: L{escript.Data}
1311 @note: this method has to be overwritten by a particular saddle point problem
1312 """
1313 pass
1314
1315 def inner(self,p0,p1):
1316 """
1317 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
1318 @return: inner product of p0 and p1
1319 @rtype: C{float}
1320 """
1321 pass
1322
1323 subiter_max=3
1324 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1325 """
1326 runs the solver
1327
1328 @param u0: initial guess for C{u}
1329 @type u0: L{esys.escript.Data}
1330 @param p0: initial guess for C{p}
1331 @type p0: L{esys.escript.Data}
1332 @param tolerance: tolerance for relative error in C{u} and C{p}
1333 @type tolerance: positive C{float}
1334 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
1335 @type tolerance_u: positive C{float}
1336 @param iter_max: maximum number of iteration steps.
1337 @type iter_max: C{int}
1338 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
1339 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
1340 @type accepted_reduction: positive C{float} or C{None}
1341 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
1342 @type relaxation: C{float} or C{None}
1343 """
1344 tol=1.e-2
1345 if tolerance_u==None: tolerance_u=tolerance
1346 if not relaxation==None: self.relaxation=relaxation
1347 if accepted_reduction ==None:
1348 angle_limit=0.
1349 elif accepted_reduction>=1.:
1350 angle_limit=0.
1351 else:
1352 angle_limit=util.sqrt(1-accepted_reduction**2)
1353 self.iter=0
1354 u=u0
1355 p=p0
1356 #
1357 # initialize things:
1358 #
1359 converged=False
1360 #
1361 # start loop:
1362 #
1363 # initial search direction is g
1364 #
1365 while not converged :
1366 if self.iter>iter_max:
1367 raise ArithmeticError("no convergence after %s steps."%self.iter)
1368 f_new=self.solve_f(u,p,tol)
1369 norm_f_new = util.Lsup(f_new)
1370 u_new=u-f_new
1371 g_new=self.solve_g(u_new,tol)
1372 self.iter+=1
1373 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1374 if norm_f_new==0. and norm_g_new==0.: return u, p
1375 if self.iter>1 and not accepted_reduction==None:
1376 #
1377 # did we manage to reduce the norm of G? I
1378 # if not we start a backtracking procedure
1379 #
1380 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1381 if norm_g_new > accepted_reduction * norm_g:
1382 sub_iter=0
1383 s=self.relaxation
1384 d=g
1385 g_last=g
1386 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1387 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
1388 dg= g_new-g_last
1389 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1390 rad=self.inner(g_new,dg)/self.relaxation
1391 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1392 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1393 if abs(rad) < norm_dg*norm_g_new * angle_limit:
1394 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
1395 break
1396 r=self.relaxation
1397 self.relaxation= - rad/norm_dg**2
1398 s+=self.relaxation
1399 #####
1400 # a=g_new+self.relaxation*dg/r
1401 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1402 #####
1403 g_last=g_new
1404 p+=self.relaxation*d
1405 f_new=self.solve_f(u,p,tol)
1406 u_new=u-f_new
1407 g_new=self.solve_g(u_new,tol)
1408 self.iter+=1
1409 norm_f_new = util.Lsup(f_new)
1410 norm_g_new = util.sqrt(self.inner(g_new,g_new))
1411 # print " ",sub_iter," new g norm",norm_g_new
1412 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1413 #
1414 # can we expect reduction of g?
1415 #
1416 # u_last=u_new
1417 sub_iter+=1
1418 self.relaxation=s
1419 #
1420 # check for convergence:
1421 #
1422 norm_u_new = util.Lsup(u_new)
1423 p_new=p+self.relaxation*g_new
1424 norm_p_new = util.sqrt(self.inner(p_new,p_new))
1425 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))
1426
1427 if self.iter>1:
1428 dg2=g_new-g
1429 df2=f_new-f
1430 norm_dg2=util.sqrt(self.inner(dg2,dg2))
1431 norm_df2=util.Lsup(df2)
1432 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1433 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1434 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1435 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1436 converged=True
1437 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
1438 self.trace("convergence after %s steps."%self.iter)
1439 return u,p
1440 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
1441 # tol=1.e-2
1442 # iter=0
1443 # converged=False
1444 # u=u0*1.
1445 # p=p0*1.
1446 # while not converged and iter<iter_max:
1447 # du=self.solve_f(u,p,tol)
1448 # u-=du
1449 # norm_du=util.Lsup(du)
1450 # norm_u=util.Lsup(u)
1451 #
1452 # dp=self.relaxation*self.solve_g(u,tol)
1453 # p+=dp
1454 # norm_dp=util.sqrt(self.inner(dp,dp))
1455 # norm_p=util.sqrt(self.inner(p,p))
1456 # 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)
1457 # iter+=1
1458 #
1459 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
1460 # if converged:
1461 # print "convergence after %s steps."%iter
1462 # else:
1463 # raise ArithmeticError("no convergence after %s steps."%iter)
1464 #
1465 # return u,p
1466
1467 def MaskFromBoundaryTag(function_space,*tags):
1468 """
1469 create a mask on the given function space which one for samples
1470 that touch the boundary tagged by tags.
1471
1472 usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1473
1474 @param function_space: a given function space
1475 @type function_space: L{escript.FunctionSpace}
1476 @param tags: boundray tags
1477 @type tags: C{str}
1478 @return: a mask which marks samples used by C{function_space} that are touching the
1479 boundary tagged by any of the given tags.
1480 @rtype: L{escript.Data} of rank 0
1481 """
1482 pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1483 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1484 for t in tags: d.setTaggedValue(t,1.)
1485 pde.setValue(y=d)
1486 out=util.whereNonZero(pde.getRightHandSide())
1487 if out.getFunctionSpace() == function_space:
1488 return out
1489 else:
1490 return util.whereNonZero(util.interpolate(out,function_space))
1491
1492
1493

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26