/[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 1414 - (show annotations)
Thu Feb 14 10:01:43 2008 UTC (11 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 33536 byte(s)
a first verion of a Stokes solver
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 HomogeneousSaddlePointProblem(object):
625 """
626 This provides a framwork for solving homogeneous saddle point problem of the form
627
628 Av+B^*p=f
629 Bv =0
630
631 for the unknowns v and p and given operators A and B and given right hand side f.
632 B^* is the adjoint operator of B is the given inner product.
633
634 """
635 def __init__(self,**kwargs):
636 self.setTolerance()
637 self.setToleranceReductionFactor()
638
639 def initialize(self):
640 """
641 initialize the problem (overwrite)
642 """
643 pass
644 def B(self,v):
645 """
646 returns Bv (overwrite)
647 @rtype: equal to the type of p
648
649 @note: boundary conditions on p should be zero!
650 """
651 pass
652
653 def inner(self,p0,p1):
654 """
655 returns inner product of two element p0 and p1 (overwrite)
656
657 @type p0: equal to the type of p
658 @type p1: equal to the type of p
659 @rtype: C{float}
660
661 @rtype: equal to the type of p
662 """
663 pass
664
665 def solve_A(self,u,p):
666 """
667 solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
668
669 @rtype: equal to the type of v
670 @note: boundary conditions on v should be zero!
671 """
672 pass
673
674 def solve_prec(self,p):
675 """
676 provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
677
678 @rtype: equal to the type of p
679 """
680 pass
681
682 def stoppingcriterium(self,Bv,v,p):
683 """
684 returns a True if iteration is terminated. (overwrite)
685
686 @rtype: C{bool}
687 """
688 pass
689
690 def __inner(self,p,r):
691 return self.inner(p,r[1])
692
693 def __stoppingcriterium(self,norm_r,r,p):
694 return self.stoppingcriterium(r[1],r[0],p)
695
696 def setTolerance(self,tolerance=1.e-8):
697 self.__tol=tolerance
698 def getTolerance(self):
699 return self.__tol
700 def setToleranceReductionFactor(self,reduction=0.01):
701 self.__reduction=reduction
702 def getSubProblemTolerance(self):
703 return self.__reduction*self.getTolerance()
704
705 def solve(self,v,p,max_iter=20, verbose=False, show_details=False):
706 """
707 solves the saddle point problem using initial guesses v and p.
708
709 @param max_iter: maximum number of iteration steps.
710 """
711 self.verbose=verbose
712 self.show_details=show_details and self.verbose
713
714 # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
715
716 self.__z=v+self.solve_A(v,p*0)
717 Bz=self.B(self.__z)
718 #
719 # solve BA^-1B^*p = Bz
720 #
721 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
722 #
723 # with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
724 # A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
725 #
726 self.iter=0
727 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
728 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
729 return r[0],p
730
731 def __Msolve(self,r):
732 return self.solve_prec(r[1])
733
734 def __Aprod(self,p):
735 # return BA^-1B*p
736 #solve Av =-B^*p as Av =f-Az-B^*p
737 v=self.solve_A(self.__z,p)
738 return ArithmeticTuple(v, self.B(v))
739
740
741 class SaddlePointProblem(object):
742 """
743 This implements a solver for a saddlepoint problem
744
745 M{f(u,p)=0}
746 M{g(u)=0}
747
748 for u and p. The problem is solved with an inexact Uszawa scheme for p:
749
750 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
751 M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
752
753 where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
754 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'
755 Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
756 in fact the role of a preconditioner.
757 """
758 def __init__(self,verbose=False,*args):
759 """
760 initializes the problem
761
762 @param verbose: switches on the printing out some information
763 @type verbose: C{bool}
764 @note: this method may be overwritten by a particular saddle point problem
765 """
766 if not isinstance(verbose,bool):
767 raise TypeError("verbose needs to be of type bool.")
768 self.__verbose=verbose
769 self.relaxation=1.
770
771 def trace(self,text):
772 """
773 prints text if verbose has been set
774
775 @param text: a text message
776 @type text: C{str}
777 """
778 if self.__verbose: print "%s: %s"%(str(self),text)
779
780 def solve_f(self,u,p,tol=1.e-8):
781 """
782 solves
783
784 A_f du = f(u,p)
785
786 with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
787
788 @param u: current approximation of u
789 @type u: L{escript.Data}
790 @param p: current approximation of p
791 @type p: L{escript.Data}
792 @param tol: tolerance expected for du
793 @type tol: C{float}
794 @return: increment du
795 @rtype: L{escript.Data}
796 @note: this method has to be overwritten by a particular saddle point problem
797 """
798 pass
799
800 def solve_g(self,u,tol=1.e-8):
801 """
802 solves
803
804 Q_g dp = g(u)
805
806 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.
807
808 @param u: current approximation of u
809 @type u: L{escript.Data}
810 @param tol: tolerance expected for dp
811 @type tol: C{float}
812 @return: increment dp
813 @rtype: L{escript.Data}
814 @note: this method has to be overwritten by a particular saddle point problem
815 """
816 pass
817
818 def inner(self,p0,p1):
819 """
820 inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
821 @return: inner product of p0 and p1
822 @rtype: C{float}
823 """
824 pass
825
826 subiter_max=3
827 def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
828 """
829 runs the solver
830
831 @param u0: initial guess for C{u}
832 @type u0: L{esys.escript.Data}
833 @param p0: initial guess for C{p}
834 @type p0: L{esys.escript.Data}
835 @param tolerance: tolerance for relative error in C{u} and C{p}
836 @type tolerance: positive C{float}
837 @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
838 @type tolerance_u: positive C{float}
839 @param iter_max: maximum number of iteration steps.
840 @type iter_max: C{int}
841 @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
842 relaxation factor. If C{accepted_reduction=None} no backtracking is used.
843 @type accepted_reduction: positive C{float} or C{None}
844 @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
845 @type relaxation: C{float} or C{None}
846 """
847 tol=1.e-2
848 if tolerance_u==None: tolerance_u=tolerance
849 if not relaxation==None: self.relaxation=relaxation
850 if accepted_reduction ==None:
851 angle_limit=0.
852 elif accepted_reduction>=1.:
853 angle_limit=0.
854 else:
855 angle_limit=util.sqrt(1-accepted_reduction**2)
856 self.iter=0
857 u=u0
858 p=p0
859 #
860 # initialize things:
861 #
862 converged=False
863 #
864 # start loop:
865 #
866 # initial search direction is g
867 #
868 while not converged :
869 if self.iter>iter_max:
870 raise ArithmeticError("no convergence after %s steps."%self.iter)
871 f_new=self.solve_f(u,p,tol)
872 norm_f_new = util.Lsup(f_new)
873 u_new=u-f_new
874 g_new=self.solve_g(u_new,tol)
875 self.iter+=1
876 norm_g_new = util.sqrt(self.inner(g_new,g_new))
877 if norm_f_new==0. and norm_g_new==0.: return u, p
878 if self.iter>1 and not accepted_reduction==None:
879 #
880 # did we manage to reduce the norm of G? I
881 # if not we start a backtracking procedure
882 #
883 # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
884 if norm_g_new > accepted_reduction * norm_g:
885 sub_iter=0
886 s=self.relaxation
887 d=g
888 g_last=g
889 self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
890 while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
891 dg= g_new-g_last
892 norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
893 rad=self.inner(g_new,dg)/self.relaxation
894 # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
895 # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
896 if abs(rad) < norm_dg*norm_g_new * angle_limit:
897 if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
898 break
899 r=self.relaxation
900 self.relaxation= - rad/norm_dg**2
901 s+=self.relaxation
902 #####
903 # a=g_new+self.relaxation*dg/r
904 # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
905 #####
906 g_last=g_new
907 p+=self.relaxation*d
908 f_new=self.solve_f(u,p,tol)
909 u_new=u-f_new
910 g_new=self.solve_g(u_new,tol)
911 self.iter+=1
912 norm_f_new = util.Lsup(f_new)
913 norm_g_new = util.sqrt(self.inner(g_new,g_new))
914 # print " ",sub_iter," new g norm",norm_g_new
915 self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
916 #
917 # can we expect reduction of g?
918 #
919 # u_last=u_new
920 sub_iter+=1
921 self.relaxation=s
922 #
923 # check for convergence:
924 #
925 norm_u_new = util.Lsup(u_new)
926 p_new=p+self.relaxation*g_new
927 norm_p_new = util.sqrt(self.inner(p_new,p_new))
928 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))
929
930 if self.iter>1:
931 dg2=g_new-g
932 df2=f_new-f
933 norm_dg2=util.sqrt(self.inner(dg2,dg2))
934 norm_df2=util.Lsup(df2)
935 # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
936 tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
937 tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
938 if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
939 converged=True
940 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
941 self.trace("convergence after %s steps."%self.iter)
942 return u,p
943 # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
944 # tol=1.e-2
945 # iter=0
946 # converged=False
947 # u=u0*1.
948 # p=p0*1.
949 # while not converged and iter<iter_max:
950 # du=self.solve_f(u,p,tol)
951 # u-=du
952 # norm_du=util.Lsup(du)
953 # norm_u=util.Lsup(u)
954 #
955 # dp=self.relaxation*self.solve_g(u,tol)
956 # p+=dp
957 # norm_dp=util.sqrt(self.inner(dp,dp))
958 # norm_p=util.sqrt(self.inner(p,p))
959 # 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)
960 # iter+=1
961 #
962 # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
963 # if converged:
964 # print "convergence after %s steps."%iter
965 # else:
966 # raise ArithmeticError("no convergence after %s steps."%iter)
967 #
968 # return u,p
969
970 def MaskFromBoundaryTag(function_space,*tags):
971 """
972 create a mask on the given function space which one for samples
973 that touch the boundary tagged by tags.
974
975 usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
976
977 @param function_space: a given function space
978 @type function_space: L{escript.FunctionSpace}
979 @param tags: boundray tags
980 @type tags: C{str}
981 @return: a mask which marks samples used by C{function_space} that are touching the
982 boundary tagged by any of the given tags.
983 @rtype: L{escript.Data} of rank 0
984 """
985 pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
986 d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
987 for t in tags: d.setTaggedValue(t,1.)
988 pde.setValue(y=d)
989 out=util.whereNonZero(pde.getRightHandSide())
990 if out.getFunctionSpace() == function_space:
991 return out
992 else:
993 return util.whereNonZero(util.interpolate(out,function_space))
994
995

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26