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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 351 by gross, Tue Dec 13 09:12:15 2005 UTC revision 1330 by gross, Mon Oct 22 04:54:49 2007 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $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.  Provides some tools related to PDEs.
# Line 7  Currently includes: Line 21  Currently includes:
21      - Projector - to project a discontinuous      - Projector - to project a discontinuous
22      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
23      - TimeIntegrationManager - to handel extraplotion in time      - 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  import escript
46  import linearPDEs  import linearPDEs
47  import numarray  import numarray
48  import util  import util
49    import math
50    
51  class TimeIntegrationManager:  class TimeIntegrationManager:
52    """    """
53    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
54    
55    typical usage is:    typical usage is::
56    
57    dt=0.1 # time increment       dt=0.1 # time increment
58    tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
59    while t<1.       while t<1.
60        v_guess=tm.extrapolate(dt) # extrapolate to t+dt           v_guess=tm.extrapolate(dt) # extrapolate to t+dt
61        v=...           v=...
62        tm.checkin(dt,v)           tm.checkin(dt,v)
63        t+=dt           t+=dt
64    
65    @remark: currently only p=1 is supported.    @note: currently only p=1 is supported.
66    """    """
67    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
68       """       """
# Line 49  class TimeIntegrationManager: Line 83  class TimeIntegrationManager:
83    
84    def getTime(self):    def getTime(self):
85        return self.__t        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):    def checkin(self,dt,*values):
94        """        """
# Line 84  class TimeIntegrationManager: Line 124  class TimeIntegrationManager:
124           return out[0]           return out[0]
125        else:        else:
126           return out           return out
127    
128    
129  class Projector:  class Projector:
130    """    """
# Line 106  class Projector: Line 147  class Projector:
147      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
148      return      return
149    
   def __del__(self):  
     return  
   
150    def __call__(self, input_data):    def __call__(self, input_data):
151      """      """
152      Projects input_data onto a continuous function      Projects input_data onto a continuous function
153    
154      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
155      """      """
156      out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))      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:      if input_data.getRank()==0:
159          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
160          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 143  class Projector: Line 182  class Projector:
182                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
183      return out      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:  class Locator:
299       """       """
300       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given
# Line 159  class Locator: Line 310  class Locator:
310         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
311         or FunctionSpace where for the sample point which         or FunctionSpace where for the sample point which
312         closest to the given point x.         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):         if isinstance(where,escript.FunctionSpace):
320            self.__function_space=where            self.__function_space=where
321         else:         else:
322            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
323         self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()         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):       def __str__(self):
331         """         """
332         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
333         """         """
334         return "<Locator %s>"%str(self.getX())         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):       def getFunctionSpace(self):
356          """          """
# Line 178  class Locator: Line 358  class Locator:
358      """      """
359          return self.__function_space          return self.__function_space
360    
361       def getId(self):       def getId(self,item=None):
362          """          """
363      Returns the identifier of the location.      Returns the identifier of the location.
364      """      """
365          return self.__id          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    
      def getX(self):  
         """  
     Returns the exact coordinates of the Locator.  
     """  
         return self(self.getFunctionSpace().getX())  
373    
374       def __call__(self,data):       def __call__(self,data):
375          """          """
# Line 204  class Locator: Line 385  class Locator:
385      """      """
386          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
387             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
388               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
389             else:             else:
390               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
391             if data.getRank()==0:             id=self.getId()
392                return out[0]             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:             else:
403                return out               out=data.getValueOfGlobalDataPoint(*id)
404                 if data.getRank()==0:
405                    return out[0]
406                 else:
407                    return out
408          else:          else:
409             return data             return data
410    
411  # vim: expandtab shiftwidth=4:  class SolverSchemeException(Exception):
412       """
413       exceptions thrown by solvers
414       """
415       pass
416    
417    class IndefinitePreconditioner(SolverSchemeException):
418       """
419       the preconditioner is not positive definite.
420       """
421       pass
422    class MaxIterReached(SolverSchemeException):
423       """
424       maxium number of iteration steps is reached.
425       """
426       pass
427    class IterationBreakDown(SolverSchemeException):
428       """
429       iteration scheme econouters an incurable breakdown.
430       """
431       pass
432    class NegativeNorm(SolverSchemeException):
433       """
434       a norm calculation returns a negative norm.
435       """
436       pass
437    
438    class IterationHistory(object):
439       """
440       The IterationHistory class is used to define a stopping criterium. It keeps track of the
441       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
442       a given tolerance.
443       """
444       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
445          """
446          Initialization
447    
448          @param tolerance: tolerance
449          @type tolerance: positive C{float}
450          @param verbose: switches on the printing out some information
451          @type verbose: C{bool}
452          """
453          if not tolerance>0.:
454              raise ValueError,"tolerance needs to be positive."
455          self.tolerance=tolerance
456          self.verbose=verbose
457          self.history=[]
458       def stoppingcriterium(self,norm_r,r,x):
459           """
460           returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.
461    
462          
463           @param norm_r: current residual norm
464           @type norm_r: non-negative C{float}
465           @param r: current residual (not used)
466           @param x: current solution approximation (not used)
467           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
468           @rtype: C{bool}
469    
470           """
471           self.history.append(norm_r)
472           if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])
473           return self.history[-1]<=self.tolerance * self.history[0]
474    
475    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
476       """
477       Solver for
478    
479       M{Ax=b}
480    
481       with a symmetric and positive definite operator A (more details required!).
482       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
483    
484       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
485    
486       For details on the preconditioned conjugate gradient method see the book:
487    
488       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
489       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
490       C. Romine, and H. van der Vorst.
491    
492       @param b: the right hand side of the liner system. C{b} is altered.
493       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
494       @param Aprod: returns the value Ax
495       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.
496       @param Msolve: solves Mx=r
497       @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same
498    type like argument C{x}.
499       @param bilinearform: inner product C{<x,r>}
500       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.
501       @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.
502       @type stoppingcriterium: function that returns C{True} or C{False}
503       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
504       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
505       @param iter_max: maximum number of iteration steps.
506       @type iter_max: C{int}
507       @return: the solution approximation and the corresponding residual
508       @rtype: C{tuple}
509       @warning: C{b} and C{x} are altered.
510       """
511       iter=0
512       if x==None:
513          x=0*b
514       else:
515          b += (-1)*Aprod(x)
516       r=b
517       rhat=Msolve(r)
518       d = rhat
519       rhat_dot_r = bilinearform(rhat, r)
520       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521    
522       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
523           iter+=1
524           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
525    
526           q=Aprod(d)
527           alpha = rhat_dot_r / bilinearform(d, q)
528           x += alpha * d
529           r += (-alpha) * q
530    
531           rhat=Msolve(r)
532           rhat_dot_r_new = bilinearform(rhat, r)
533           beta = rhat_dot_r_new / rhat_dot_r
534           rhat+=beta * d
535           d=rhat
536    
537           rhat_dot_r = rhat_dot_r_new
538           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
539    
540       return x,r
541    
542    class SaddlePointProblem(object):
543       """
544       This implements a solver for a saddlepoint problem
545    
546       M{f(u,p)=0}
547       M{g(u)=0}
548    
549       for u and p. The problem is solved with an inexact Uszawa scheme for p:
550    
551       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
552       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
553    
554       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
555       A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
556       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
557       in fact the role of a preconditioner.
558       """
559       def __init__(self,verbose=False,*args):
560           """
561           initializes the problem
562    
563           @param verbose: switches on the printing out some information
564           @type verbose: C{bool}
565           @note: this method may be overwritten by a particular saddle point problem
566           """
567           if not isinstance(verbose,bool):
568                raise TypeError("verbose needs to be of type bool.")
569           self.__verbose=verbose
570           self.relaxation=1.
571    
572       def trace(self,text):
573           """
574           prints text if verbose has been set
575    
576           @param text: a text message
577           @type text: C{str}
578           """
579           if self.__verbose: print "%s: %s"%(str(self),text)
580    
581       def solve_f(self,u,p,tol=1.e-8):
582           """
583           solves
584    
585           A_f du = f(u,p)
586    
587           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
588    
589           @param u: current approximation of u
590           @type u: L{escript.Data}
591           @param p: current approximation of p
592           @type p: L{escript.Data}
593           @param tol: tolerance expected for du
594           @type tol: C{float}
595           @return: increment du
596           @rtype: L{escript.Data}
597           @note: this method has to be overwritten by a particular saddle point problem
598           """
599           pass
600    
601       def solve_g(self,u,tol=1.e-8):
602           """
603           solves
604    
605           Q_g dp = g(u)
606    
607           with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.
608    
609           @param u: current approximation of u
610           @type u: L{escript.Data}
611           @param tol: tolerance expected for dp
612           @type tol: C{float}
613           @return: increment dp
614           @rtype: L{escript.Data}
615           @note: this method has to be overwritten by a particular saddle point problem
616           """
617           pass
618    
619       def inner(self,p0,p1):
620           """
621           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
622           @return: inner product of p0 and p1
623           @rtype: C{float}
624           """
625           pass
626    
627       subiter_max=3
628       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
629            """
630            runs the solver
631    
632            @param u0: initial guess for C{u}
633            @type u0: L{esys.escript.Data}
634            @param p0: initial guess for C{p}
635            @type p0: L{esys.escript.Data}
636            @param tolerance: tolerance for relative error in C{u} and C{p}
637            @type tolerance: positive C{float}
638            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
639            @type tolerance_u: positive C{float}
640            @param iter_max: maximum number of iteration steps.
641            @type iter_max: C{int}
642            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
643                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
644            @type accepted_reduction: positive C{float} or C{None}
645            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
646            @type relaxation: C{float} or C{None}
647            """
648            tol=1.e-2
649            if tolerance_u==None: tolerance_u=tolerance
650            if not relaxation==None: self.relaxation=relaxation
651            if accepted_reduction ==None:
652                  angle_limit=0.
653            elif accepted_reduction>=1.:
654                  angle_limit=0.
655            else:
656                  angle_limit=util.sqrt(1-accepted_reduction**2)
657            self.iter=0
658            u=u0
659            p=p0
660            #
661            #   initialize things:
662            #
663            converged=False
664            #
665            #  start loop:
666            #
667            #  initial search direction is g
668            #
669            while not converged :
670                if self.iter>iter_max:
671                    raise ArithmeticError("no convergence after %s steps."%self.iter)
672                f_new=self.solve_f(u,p,tol)
673                norm_f_new = util.Lsup(f_new)
674                u_new=u-f_new
675                g_new=self.solve_g(u_new,tol)
676                self.iter+=1
677                norm_g_new = util.sqrt(self.inner(g_new,g_new))
678                if norm_f_new==0. and norm_g_new==0.: return u, p
679                if self.iter>1 and not accepted_reduction==None:
680                   #
681                   #   did we manage to reduce the norm of G? I
682                   #   if not we start a backtracking procedure
683                   #
684                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
685                   if norm_g_new > accepted_reduction * norm_g:
686                      sub_iter=0
687                      s=self.relaxation
688                      d=g
689                      g_last=g
690                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
691                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
692                         dg= g_new-g_last
693                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
694                         rad=self.inner(g_new,dg)/self.relaxation
695                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
696                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
697                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
698                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
699                             break
700                         r=self.relaxation
701                         self.relaxation= - rad/norm_dg**2
702                         s+=self.relaxation
703                         #####
704                         # a=g_new+self.relaxation*dg/r
705                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
706                         #####
707                         g_last=g_new
708                         p+=self.relaxation*d
709                         f_new=self.solve_f(u,p,tol)
710                         u_new=u-f_new
711                         g_new=self.solve_g(u_new,tol)
712                         self.iter+=1
713                         norm_f_new = util.Lsup(f_new)
714                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
715                         # print "   ",sub_iter," new g norm",norm_g_new
716                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
717                         #
718                         #   can we expect reduction of g?
719                         #
720                         # u_last=u_new
721                         sub_iter+=1
722                      self.relaxation=s
723                #
724                #  check for convergence:
725                #
726                norm_u_new = util.Lsup(u_new)
727                p_new=p+self.relaxation*g_new
728                norm_p_new = util.sqrt(self.inner(p_new,p_new))
729                self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
730    
731                if self.iter>1:
732                   dg2=g_new-g
733                   df2=f_new-f
734                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
735                   norm_df2=util.Lsup(df2)
736                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
737                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
738                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
739                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
740                       converged=True
741                f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
742            self.trace("convergence after %s steps."%self.iter)
743            return u,p
744    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
745    #      tol=1.e-2
746    #      iter=0
747    #      converged=False
748    #      u=u0*1.
749    #      p=p0*1.
750    #      while not converged and iter<iter_max:
751    #          du=self.solve_f(u,p,tol)
752    #          u-=du
753    #          norm_du=util.Lsup(du)
754    #          norm_u=util.Lsup(u)
755    #        
756    #          dp=self.relaxation*self.solve_g(u,tol)
757    #          p+=dp
758    #          norm_dp=util.sqrt(self.inner(dp,dp))
759    #          norm_p=util.sqrt(self.inner(p,p))
760    #          print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
761    #          iter+=1
762    #
763    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
764    #      if converged:
765    #          print "convergence after %s steps."%iter
766    #      else:
767    #          raise ArithmeticError("no convergence after %s steps."%iter)
768    #
769    #      return u,p
770              
771    def MaskFromBoundaryTag(function_space,*tags):
772       """
773       create a mask on the given function space which one for samples
774       that touch the boundary tagged by tags.
775    
776       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
777    
778       @param function_space: a given function space
779       @type function_space: L{escript.FunctionSpace}
780       @param tags: boundray tags
781       @type tags: C{str}
782       @return: a mask which marks samples used by C{function_space} that are touching the
783                boundary tagged by any of the given tags.
784       @rtype: L{escript.Data} of rank 0
785       """
786       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
787       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
788       for t in tags: d.setTaggedValue(t,1.)
789       pde.setValue(y=d)
790       out=util.whereNonZero(pde.getRightHandSide())
791       if out.getFunctionSpace() == function_space:
792          return out
793       else:
794          return util.whereNonZero(util.interpolate(out,function_space))
795    

Legend:
Removed from v.351  
changed lines
  Added in v.1330

  ViewVC Help
Powered by ViewVC 1.1.26