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

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

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

trunk/esys2/escript/py_src/pdetools.py revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC trunk/escript/py_src/pdetools.py revision 1312 by ksteube, Mon Sep 24 06:18:44 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 a some tools related to PDEs currently includes:  Provides some tools related to PDEs.
19    
20       Projector - to project a discontinuous  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  from escript import *  import escript
46  from linearPDEs import LinearPDE  import linearPDEs
47  import numarray  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:  class Projector:
130    """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""    """
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):    def __init__(self, domain, reduce = True, fast=True):
135      """      """
136      @brief Create a continuous function space projector for a domain.      Create a continuous function space projector for a domain.
137    
138      @param domain Domain of the projection.      @param domain: Domain of the projection.
139      @param reduce Flag to reduce projection order (default is True)      @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)      @param fast: Flag to use a fast method based on matrix lumping (default is true)
141      """      """
142      self.__pde = LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
143      self.__pde.setLumping(fast)      if fast:
144          self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
145      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
146      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
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      @brief 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=Data(0.,input_data.getShape(),what=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 66  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:
 class Locator:  
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          Locator provides a  access the values of data objects at a given spatial coordinate x.       where D, Y, r and q are given scalar functions.
         In fact, a Locator object finds the sample in the set of samples of a given function space or domain where which is closest  
         to the given point x.  
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,))):       def __init__(self,where,x=numarray.zeros((3,))):
309         """initializes a Locator to access values in Data objects on the Doamin or FunctionSpace where for the sample point which         """
310            closest to the given point x"""         Initializes a Locator to access values in Data objects on the Doamin
311         if isinstance(where,FunctionSpace):         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            self.__function_space=where
321         else:         else:
322            self.__function_space=ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
323         self.__id=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         """returns the coordinates of the Locator as a string"""         """
332         return "<Locator %s>"%str(self.getX())         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):       def getFunctionSpace(self):
356          """returns the function space of the Locator"""          """
357        Returns the function space of the Locator.
358        """
359          return self.__function_space          return self.__function_space
360    
361       def getId(self):       def getId(self,item=None):
362          """returns the identifier of the location"""          """
363          return self.__id      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    
      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          """returns the value of data at the Locator of a Data object otherwise the object is returned."""          """
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)          return self.getValue(data)
380    
381       def getValue(self,data):       def getValue(self,data):
382          """returns the value of data at the Locator if data is a Data object otherwise the object is returned."""          """
383          if isinstance(data,Data):      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():             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    class SolverSchemeException(Exception):
412       """
413       exceptions thrown by solvers
414       """
415       pass
416    
417    class IndefinitePreconditioner(SolverSchemeException):
418       """
419       the preconditioner is not positive definite.
420       """
421       pass
422    class MaxIterReached(SolverSchemeException):
423       """
424       maxium number of iteration steps is reached.
425       """
426       pass
427    class IterationBreakDown(SolverSchemeException):
428       """
429       iteration scheme econouters an incurable breakdown.
430       """
431       pass
432    class NegativeNorm(SolverSchemeException):
433       """
434       a norm calculation returns a negative norm.
435       """
436       pass
437    
438    def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):
439       """
440       Solver for
441    
442       M{Ax=b}
443    
444       with a symmetric and positive definite operator A (more details required!).
445       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
446    
447       The iteration is terminated if
448    
449       M{norm(r) <= tolerance * norm(b)}
450    
451       where C{norm()} defines a norm and
452    
453       M{r = b-Ax}
454    
455       the residual.
456    
457       For details on the preconditioned conjugate gradient method see the book:
458    
459       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
460       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
461       C. Romine, and H. van der Vorst.
462    
463       @param b: the right hand side of the liner system. C{b} is altered.
464       @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)
465       @param x: an initial guess for the solution
466       @type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)
467       @param Aprod: returns the value Ax
468       @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.
469       @param Msolve: solves Mx=r
470       @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.
471       @param bilinearform: inner product C{<x,r>}
472       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.
473       @param norm: norm calculation for the residual C{r=b-Ax}.
474       @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.
475       @param verbose: switches on the printing out some information
476       @type verbose: C{bool}
477       @param iter_max: maximum number of iteration steps.
478       @type iter_max: C{int}
479       @param tolerance: tolerance
480       @type tolerance: positive C{float}
481       @return: the solution apprximation and the corresponding residual
482       @rtype: C{tuple} of an S type and and an R type object.A
483       @warning: C{b} ans C{x} are altered.
484       """
485       if verbose:
486            print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance   =%e"%(iter_max,tolerance)
487       iter=0
488       normb = norm(b)
489       if normb<0: raise NegativeNorm
490    
491       b += (-1)*Aprod(x)
492       r=b
493       rhat=Msolve(r)
494       d = rhat;
495       rhat_dot_r = bilinearform(rhat, r)
496    
497       while True:
498           normr=norm(r)
499           if normr<0: raise NegativeNorm
500           if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)
501           if normr <= tolerance * normb: return x,r
502    
503           iter+=1 # k = iter = 1 first time through
504           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
505    
506           q=Aprod(d)
507           alpha = rhat_dot_r / bilinearform(d, q)
508           x += alpha * d
509           r += (-alpha) * q
510    
511           rhat=Msolve(r)
512           rhat_dot_r_new = bilinearform(rhat, r)
513           beta = rhat_dot_r_new / rhat_dot_r
514           rhat+=beta * d
515           d=rhat
516    
517           rhat_dot_r = rhat_dot_r_new
518    
519    
520    class SaddlePointProblem(object):
521       """
522       This implements a solver for a saddlepoint problem
523    
524       M{f(u,p)=0}
525       M{g(u)=0}
526    
527       for u and p. The problem is solved with an inexact Uszawa scheme for p:
528    
529       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
530       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
531    
532       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
533       A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'
534       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
535       in fact the role of a preconditioner.
536       """
537       def __init__(self,verbose=False,*args):
538           """
539           initializes the problem
540    
541           @param verbose: switches on the printing out some information
542           @type verbose: C{bool}
543           @note: this method may be overwritten by a particular saddle point problem
544           """
545           if not isinstance(verbose,bool):
546                raise TypeError("verbose needs to be of type bool.")
547           self.__verbose=verbose
548           self.relaxation=1.
549    
550       def trace(self,text):
551           """
552           prints text if verbose has been set
553    
554           @param text: a text message
555           @type text: C{str}
556           """
557           if self.__verbose: print "%s: %s"%(str(self),text)
558    
559       def solve_f(self,u,p,tol=1.e-8):
560           """
561           solves
562    
563           A_f du = f(u,p)
564    
565           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
566    
567           @param u: current approximation of u
568           @type u: L{escript.Data}
569           @param p: current approximation of p
570           @type p: L{escript.Data}
571           @param tol: tolerance expected for du
572           @type tol: C{float}
573           @return: increment du
574           @rtype: L{escript.Data}
575           @note: this method has to be overwritten by a particular saddle point problem
576           """
577           pass
578    
579       def solve_g(self,u,tol=1.e-8):
580           """
581           solves
582    
583           Q_g dp = g(u)
584    
585           with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.
586    
587           @param u: current approximation of u
588           @type u: L{escript.Data}
589           @param tol: tolerance expected for dp
590           @type tol: C{float}
591           @return: increment dp
592           @rtype: L{escript.Data}
593           @note: this method has to be overwritten by a particular saddle point problem
594           """
595           pass
596    
597       def inner(self,p0,p1):
598           """
599           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
600           @return: inner product of p0 and p1
601           @rtype: C{float}
602           """
603           pass
604    
605       subiter_max=3
606       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
607            """
608            runs the solver
609    
610            @param u0: initial guess for C{u}
611            @type u0: L{esys.escript.Data}
612            @param p0: initial guess for C{p}
613            @type p0: L{esys.escript.Data}
614            @param tolerance: tolerance for relative error in C{u} and C{p}
615            @type tolerance: positive C{float}
616            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
617            @type tolerance_u: positive C{float}
618            @param iter_max: maximum number of iteration steps.
619            @type iter_max: C{int}
620            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
621                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
622            @type accepted_reduction: positive C{float} or C{None}
623            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
624            @type relaxation: C{float} or C{None}
625            """
626            tol=1.e-2
627            if tolerance_u==None: tolerance_u=tolerance
628            if not relaxation==None: self.relaxation=relaxation
629            if accepted_reduction ==None:
630                  angle_limit=0.
631            elif accepted_reduction>=1.:
632                  angle_limit=0.
633            else:
634                  angle_limit=util.sqrt(1-accepted_reduction**2)
635            self.iter=0
636            u=u0
637            p=p0
638            #
639            #   initialize things:
640            #
641            converged=False
642            #
643            #  start loop:
644            #
645            #  initial search direction is g
646            #
647            while not converged :
648                if self.iter>iter_max:
649                    raise ArithmeticError("no convergence after %s steps."%self.iter)
650                f_new=self.solve_f(u,p,tol)
651                norm_f_new = util.Lsup(f_new)
652                u_new=u-f_new
653                g_new=self.solve_g(u_new,tol)
654                self.iter+=1
655                norm_g_new = util.sqrt(self.inner(g_new,g_new))
656                if norm_f_new==0. and norm_g_new==0.: return u, p
657                if self.iter>1 and not accepted_reduction==None:
658                   #
659                   #   did we manage to reduce the norm of G? I
660                   #   if not we start a backtracking procedure
661                   #
662                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
663                   if norm_g_new > accepted_reduction * norm_g:
664                      sub_iter=0
665                      s=self.relaxation
666                      d=g
667                      g_last=g
668                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
669                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
670                         dg= g_new-g_last
671                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
672                         rad=self.inner(g_new,dg)/self.relaxation
673                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
674                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
675                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
676                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
677                             break
678                         r=self.relaxation
679                         self.relaxation= - rad/norm_dg**2
680                         s+=self.relaxation
681                         #####
682                         # a=g_new+self.relaxation*dg/r
683                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
684                         #####
685                         g_last=g_new
686                         p+=self.relaxation*d
687                         f_new=self.solve_f(u,p,tol)
688                         u_new=u-f_new
689                         g_new=self.solve_g(u_new,tol)
690                         self.iter+=1
691                         norm_f_new = util.Lsup(f_new)
692                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
693                         # print "   ",sub_iter," new g norm",norm_g_new
694                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
695                         #
696                         #   can we expect reduction of g?
697                         #
698                         # u_last=u_new
699                         sub_iter+=1
700                      self.relaxation=s
701                #
702                #  check for convergence:
703                #
704                norm_u_new = util.Lsup(u_new)
705                p_new=p+self.relaxation*g_new
706                norm_p_new = util.sqrt(self.inner(p_new,p_new))
707                self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
708    
709                if self.iter>1:
710                   dg2=g_new-g
711                   df2=f_new-f
712                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
713                   norm_df2=util.Lsup(df2)
714                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
715                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
716                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
717                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
718                       converged=True
719                f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
720            self.trace("convergence after %s steps."%self.iter)
721            return u,p
722    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
723    #      tol=1.e-2
724    #      iter=0
725    #      converged=False
726    #      u=u0*1.
727    #      p=p0*1.
728    #      while not converged and iter<iter_max:
729    #          du=self.solve_f(u,p,tol)
730    #          u-=du
731    #          norm_du=util.Lsup(du)
732    #          norm_u=util.Lsup(u)
733    #        
734    #          dp=self.relaxation*self.solve_g(u,tol)
735    #          p+=dp
736    #          norm_dp=util.sqrt(self.inner(dp,dp))
737    #          norm_p=util.sqrt(self.inner(p,p))
738    #          print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
739    #          iter+=1
740    #
741    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
742    #      if converged:
743    #          print "convergence after %s steps."%iter
744    #      else:
745    #          raise ArithmeticError("no convergence after %s steps."%iter)
746    #
747    #      return u,p
748              
749    def MaskFromBoundaryTag(function_space,*tags):
750       """
751       create a mask on the given function space which one for samples
752       that touch the boundary tagged by tags.
753    
754       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
755    
756       @param function_space: a given function space
757       @type function_space: L{escript.FunctionSpace}
758       @param tags: boundray tags
759       @type tags: C{str}
760       @return: a mask which marks samples used by C{function_space} that are touching the
761                boundary tagged by any of the given tags.
762       @rtype: L{escript.Data} of rank 0
763       """
764       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
765       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
766       for t in tags: d.setTaggedValue(t,1.)
767       pde.setValue(y=d)
768       out=util.whereNonZero(pde.getRightHandSide())
769       if out.getFunctionSpace() == function_space:
770          return out
771       else:
772          return util.whereNonZero(util.interpolate(out,function_space))
773    

Legend:
Removed from v.147  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26