/[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

trunk/esys2/escript/py_src/pdetools.py revision 146 by jgs, Fri Jul 29 01:44:29 2005 UTC trunk/escript/py_src/pdetools.py revision 1107 by gross, Thu Apr 19 02:14:18 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    
3  """  """
4  provides a some tools related to PDEs currently includes:  Provides some tools related to PDEs.
5    
6       Projector - to project a discontinuous  Currently includes:
7        - Projector - to project a discontinuous
8        - Locator - to trace values in data objects at a certain location
9        - TimeIntegrationManager - to handel extraplotion in time
10        - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
11    
12    @var __author__: name of author
13    @var __copyright__: copyrights
14    @var __license__: licence agreement
15    @var __url__: url entry point on documentation
16    @var __version__: version
17    @var __date__: date of the version
18    """
19    
20    __author__="Lutz Gross, l.gross@uq.edu.au"
21    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
22                        http://www.access.edu.au
23                    Primary Business: Queensland, Australia"""
24    __license__="""Licensed under the Open Software License version 3.0
25                 http://www.opensource.org/licenses/osl-3.0.php"""
26    __url__="http://www.iservo.edu.au/esys"
27    __version__="$Revision$"
28    __date__="$Date$"
29    
 """  
30    
31  from escript import *  import escript
32  from linearPDEs import LinearPDE  import linearPDEs
33  import numarray  import numarray
34    import util
35    
36    class TimeIntegrationManager:
37      """
38      a simple mechanism to manage time dependend values.
39    
40      typical usage is::
41    
42         dt=0.1 # time increment
43         tm=TimeIntegrationManager(inital_value,p=1)
44         while t<1.
45             v_guess=tm.extrapolate(dt) # extrapolate to t+dt
46             v=...
47             tm.checkin(dt,v)
48             t+=dt
49    
50      @note: currently only p=1 is supported.
51      """
52      def __init__(self,*inital_values,**kwargs):
53         """
54         sets up the value manager where inital_value is the initial value and p is order used for extrapolation
55         """
56         if kwargs.has_key("p"):
57                self.__p=kwargs["p"]
58         else:
59                self.__p=1
60         if kwargs.has_key("time"):
61                self.__t=kwargs["time"]
62         else:
63                self.__t=0.
64         self.__v_mem=[inital_values]
65         self.__order=0
66         self.__dt_mem=[]
67         self.__num_val=len(inital_values)
68    
69      def getTime(self):
70          return self.__t
71      def getValue(self):
72          out=self.__v_mem[0]
73          if len(out)==1:
74              return out[0]
75          else:
76              return out
77    
78      def checkin(self,dt,*values):
79          """
80          adds new values to the manager. the p+1 last value get lost
81          """
82          o=min(self.__order+1,self.__p)
83          self.__order=min(self.__order+1,self.__p)
84          v_mem_new=[values]
85          dt_mem_new=[dt]
86          for i in range(o-1):
87             v_mem_new.append(self.__v_mem[i])
88             dt_mem_new.append(self.__dt_mem[i])
89          v_mem_new.append(self.__v_mem[o-1])
90          self.__order=o
91          self.__v_mem=v_mem_new
92          self.__dt_mem=dt_mem_new
93          self.__t+=dt
94    
95      def extrapolate(self,dt):
96          """
97          extrapolates to dt forward in time.
98          """
99          if self.__order==0:
100             out=self.__v_mem[0]
101          else:
102            out=[]
103            for i in range(self.__num_val):
104               out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
105    
106          if len(out)==0:
107             return None
108          elif len(out)==1:
109             return out[0]
110          else:
111             return out
112    
113    
114  class Projector:  class Projector:
115    """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""    """
116      The Projector is a factory which projects a discontiuous function onto a
117      continuous function on the a given domain.
118      """
119    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce = True, fast=True):
120      """      """
121      @brief Create a continuous function space projector for a domain.      Create a continuous function space projector for a domain.
122    
123      @param domain Domain of the projection.      @param domain: Domain of the projection.
124      @param reduce Flag to reduce projection order (default is True)      @param reduce: Flag to reduce projection order (default is True)
125      @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)
126      """      """
127      self.__pde = LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
128      self.__pde.setLumping(fast)      if fast:
129          self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
130      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
131      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
132      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
133      return      return
134    
   def __del__(self):  
     return  
   
135    def __call__(self, input_data):    def __call__(self, input_data):
136      """      """
137      @brief projects input_data onto a continuous function      Projects input_data onto a continuous function
138    
139      @param input_data  The input_data to be projected.      @param input_data: The input_data to be projected.
140      """      """
141      out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
142      if input_data.getRank()==0:      if input_data.getRank()==0:
143          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
144          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 66  class Projector: Line 166  class Projector:
166                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
167      return out      return out
168    
169    class NoPDE:
170         """
171         solves the following problem for u:
172    
173         M{kronecker[i,j]*D[j]*u[j]=Y[i]}
174    
175         with constraint
176    
177         M{u[j]=r[j]}  where M{q[j]>0}
178    
179         where D, Y, r and q are given functions of rank 1.
180    
181         In the case of scalars this takes the form
182    
183         M{D*u=Y}
184    
185         with constraint
186    
187         M{u=r}  where M{q>0}
188    
189         where D, Y, r and q are given scalar functions.
190    
191         The constraint is overwriting any other condition.
192    
193         @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
194                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
195                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
196         """
197         def __init__(self,domain,D=None,Y=None,q=None,r=None):
198             """
199             initialize the problem
200    
201             @param domain: domain of the PDE.
202             @type domain: L{Domain}
203             @param D: coefficient of the solution.
204             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
205             @param Y: right hand side
206             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207             @param q: location of constraints
208             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209             @param r: value of solution at locations of constraints
210             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211             """
212             self.__domain=domain
213             self.__D=D
214             self.__Y=Y
215             self.__q=q
216             self.__r=r
217             self.__u=None
218             self.__function_space=escript.Solution(self.__domain)
219         def setReducedOn(self):
220             """
221             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
222             """
223             self.__function_space=escript.ReducedSolution(self.__domain)
224             self.__u=None
225    
226         def setReducedOff(self):
227             """
228             sets the L{FunctionSpace} of the solution to L{Solution}
229             """
230             self.__function_space=escript.Solution(self.__domain)
231             self.__u=None
232            
233         def setValue(self,D=None,Y=None,q=None,r=None):
234             """
235             assigns values to the parameters.
236    
237             @param D: coefficient of the solution.
238             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
239             @param Y: right hand side
240             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241             @param q: location of constraints
242             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243             @param r: value of solution at locations of constraints
244             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245             """
246             if not D==None:
247                self.__D=D
248                self.__u=None
249             if not Y==None:
250                self.__Y=Y
251                self.__u=None
252             if not q==None:
253                self.__q=q
254                self.__u=None
255             if not r==None:
256                self.__r=r
257                self.__u=None
258    
259         def getSolution(self):
260             """
261             returns the solution
262            
263             @return: the solution of the problem
264             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
265             """
266             if self.__u==None:
267                if self.__D==None:
268                   raise ValueError,"coefficient D is undefined"
269                D=escript.Data(self.__D,self.__function_space)
270                if D.getRank()>1:
271                   raise ValueError,"coefficient D must have rank 0 or 1"
272                if self.__Y==None:
273                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
274                else:
275                   self.__u=util.quotient(self.__Y,D)
276                if not self.__q==None:
277                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
278                    self.__u*=(1.-q)
279                    if not self.__r==None: self.__u+=q*self.__r
280             return self.__u
281                
282    class Locator:
283         """
284         Locator provides access to the values of data objects at a given
285         spatial coordinate x.  
286        
287         In fact, a Locator object finds the sample in the set of samples of a
288         given function space or domain where which is closest to the given
289         point x.
290         """
291    
292         def __init__(self,where,x=numarray.zeros((3,))):
293           """
294           Initializes a Locator to access values in Data objects on the Doamin
295           or FunctionSpace where for the sample point which
296           closest to the given point x.
297    
298           @param where: function space
299           @type where: L{escript.FunctionSpace}
300           @param x: coefficient of the solution.
301           @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
302           """
303           if isinstance(where,escript.FunctionSpace):
304              self.__function_space=where
305           else:
306              self.__function_space=escript.ContinuousFunction(where)
307           if isinstance(x, list):
308               self.__id=[]
309               for p in x:
310                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
311           else:
312               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
313    
314  class Location:       def __str__(self):
315       """Location provides a factory to access the values of data objects at a given spatial coordinate x.         """
316          In fact, a Location object finds the sample in the set of samples of a given function space or domain which is closest         Returns the coordinates of the Locator as a string.
317          to the given point x. """         """
318           x=self.getX()
319       def __init__(self,x,where):         if instance(x,list):
320         """initializes a Location to access values in Data objects on the Doamin or FunctionSpace where for the sample point which            out="["
321            closest to the given point x"""            first=True
322         if isinstance(where,FunctionSpace):            for xx in x:
323            self.__where=where              if not first:
324                    out+=","
325                else:
326                    first=False
327                out+=str(xx)
328              out+="]>"
329         else:         else:
330            self.__where=ContinuousFunction(where)            out=str(x)
331         self.__id=length(x-self.__where.getX()).minarg()         return out
332    
      def getValue(self,data):  
         """returns the value of data at the Location at a numarray object"""  
         if isinstance(data,Data):  
            return data  
         else:  
            if not data.getFunctionSpace()==self.getFunctionSpace():  
              raise ValueError,"function space of data obejct does not match function space of Location"  
            else:  
              return data.getValue(self.getId())  
333       def getX(self):       def getX(self):
334          """returns the exact coordinates of the Location"""          """
335          return self.getValue(self.getFunctionSpace().getX())      Returns the exact coordinates of the Locator.
336        """
337       def getId(self):          return self(self.getFunctionSpace().getX())
         """returns the identifier of the location"""  
         return self.__id  
338    
339       def getFunctionSpace(self):       def getFunctionSpace(self):
340          """returns the function space of the Location"""          """
341        Returns the function space of the Locator.
342        """
343          return self.__function_space          return self.__function_space
344    
345       def __str__(self):       def getId(self,item=None):
346         """returns the coordinates of the Location as a string"""          """
347         return str(self.getX())      Returns the identifier of the location.
348        """
349            if item == None:
350               return self.__id
351            else:
352               if isinstance(self.__id,list):
353                  return self.__id[item]
354               else:
355                  return self.__id
356    
 def testProjector(domain):  
       """runs a few test of the Projector factory and returns the largest error plus a description of the test this error occured"""  
       error_max=0.  
       error_text=""  
       x=ContinuousFunction(domain).getX()  
       for f in [True,False]:  
          p=Projector(domain,reduce=False,fast=f)  
          for r in range(5):  
             text="range %s , fast=%s"%(r,f)  
             if r==0:  
                td_ref=x[0]  
             elif r==1:  
                td_ref=x  
             elif r==2:  
                td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])  
             elif r==3:  
                td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])  
             elif r==3:  
                td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])  
             td=p(td_ref.interpolate(Function(domain)))  
             er=Lsup(td-td_ref)/Lsup(td_ref)  
             print text," error = ",er  
             if er>error_max:  
                  error_max=er  
                  error_text=text  
       return error_max,error_text  
   
   
 # this should be removed later  
 if __name__=="__main__":  
     from esys.finley import Rectangle  
     txt=testProjector(Rectangle(56,61))  
     print "test Projector: ",txt  
357    
358         def __call__(self,data):
359            """
360        Returns the value of data at the Locator of a Data object otherwise
361        the object is returned.
362        """
363            return self.getValue(data)
364    
365         def getValue(self,data):
366            """
367        Returns the value of data at the Locator if data is a Data object
368        otherwise the object is returned.
369        """
370            if isinstance(data,escript.Data):
371               if data.getFunctionSpace()==self.getFunctionSpace():
372                 dat=data
373               else:
374                 dat=data.interpolate(self.getFunctionSpace())
375               id=self.getId()
376               r=data.getRank()
377               if isinstance(id,list):
378                   out=[]
379                   for i in id:
380                      o=data.getValueOfGlobalDataPoint(*i)
381                      if data.getRank()==0:
382                         out.append(o[0])
383                      else:
384                         out.append(o)
385                   return out
386               else:
387                 out=data.getValueOfGlobalDataPoint(*id)
388                 if data.getRank()==0:
389                    return out[0]
390                 else:
391                    return out
392            else:
393               return data
394    
395    class SaddlePointProblem(object):
396       """
397       This implements a solver for a saddlepoint problem
398    
399       M{f(u,p)=0}
400       M{g(u)=0}
401    
402       for u and p. The problem is solved with an inexact Uszawa scheme for p:
403    
404       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
405       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
406    
407       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
408       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'
409       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
410       in fact the role of a preconditioner.
411       """
412       def __init__(self,verbose=False,*args):
413           """
414           initializes the problem
415    
416           @param verbose: switches on the printing out some information
417           @type verbose: C{bool}
418           @note: this method may be overwritten by a particular saddle point problem
419           """
420           if not isinstance(verbose,bool):
421                raise TypeError("verbose needs to be of type bool.")
422           self.__verbose=verbose
423           self.relaxation=1.
424    
425       def trace(self,text):
426           """
427           prints text if verbose has been set
428    
429           @param text: a text message
430           @type text: C{str}
431           """
432           if self.__verbose: print "%s: %s"%(str(self),text)
433    
434       def solve_f(self,u,p,tol=1.e-8):
435           """
436           solves
437    
438           A_f du = f(u,p)
439    
440           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
441    
442           @param u: current approximation of u
443           @type u: L{escript.Data}
444           @param p: current approximation of p
445           @type p: L{escript.Data}
446           @param tol: tolerance expected for du
447           @type tol: C{float}
448           @return: increment du
449           @rtype: L{escript.Data}
450           @note: this method has to be overwritten by a particular saddle point problem
451           """
452           pass
453    
454       def solve_g(self,u,tol=1.e-8):
455           """
456           solves
457    
458           Q_g dp = g(u)
459    
460           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.
461    
462           @param u: current approximation of u
463           @type u: L{escript.Data}
464           @param tol: tolerance expected for dp
465           @type tol: C{float}
466           @return: increment dp
467           @rtype: L{escript.Data}
468           @note: this method has to be overwritten by a particular saddle point problem
469           """
470           pass
471    
472       def inner(self,p0,p1):
473           """
474           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
475           @return: inner product of p0 and p1
476           @rtype: C{float}
477           """
478           pass
479    
480       subiter_max=3
481       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
482            """
483            runs the solver
484    
485            @param u0: initial guess for C{u}
486            @type u0: L{esys.escript.Data}
487            @param p0: initial guess for C{p}
488            @type p0: L{esys.escript.Data}
489            @param tolerance: tolerance for relative error in C{u} and C{p}
490            @type tolerance: positive C{float}
491            @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
492            @type tolerance_u: positive C{float}
493            @param iter_max: maximum number of iteration steps.
494            @type iter_max: C{int}
495            @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the
496                                       relaxation factor. If C{accepted_reduction=None} no backtracking is used.
497            @type accepted_reduction: positive C{float} or C{None}
498            @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
499            @type relaxation: C{float} or C{None}
500            """
501            tol=1.e-2
502            if tolerance_u==None: tolerance_u=tolerance
503            if not relaxation==None: self.relaxation=relaxation
504            if accepted_reduction ==None:
505                  angle_limit=0.
506            elif accepted_reduction>=1.:
507                  angle_limit=0.
508            else:
509                  angle_limit=util.sqrt(1-accepted_reduction**2)
510            self.iter=0
511            u=u0
512            p=p0
513            #
514            #   initialize things:
515            #
516            converged=False
517            #
518            #  start loop:
519            #
520            #  initial search direction is g
521            #
522            while not converged :
523                if self.iter>iter_max:
524                    raise ArithmeticError("no convergence after %s steps."%self.iter)
525                f_new=self.solve_f(u,p,tol)
526                norm_f_new = util.Lsup(f_new)
527                u_new=u-f_new
528                g_new=self.solve_g(u_new,tol)
529                self.iter+=1
530                norm_g_new = util.sqrt(self.inner(g_new,g_new))
531                if norm_f_new==0. and norm_g_new==0.: return u, p
532                if self.iter>1 and not accepted_reduction==None:
533                   #
534                   #   did we manage to reduce the norm of G? I
535                   #   if not we start a backtracking procedure
536                   #
537                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
538                   if norm_g_new > accepted_reduction * norm_g:
539                      sub_iter=0
540                      s=self.relaxation
541                      d=g
542                      g_last=g
543                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
544                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
545                         dg= g_new-g_last
546                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
547                         rad=self.inner(g_new,dg)/self.relaxation
548                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
549                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
550                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
551                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
552                             break
553                         r=self.relaxation
554                         self.relaxation= - rad/norm_dg**2
555                         s+=self.relaxation
556                         #####
557                         # a=g_new+self.relaxation*dg/r
558                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
559                         #####
560                         g_last=g_new
561                         p+=self.relaxation*d
562                         f_new=self.solve_f(u,p,tol)
563                         u_new=u-f_new
564                         g_new=self.solve_g(u_new,tol)
565                         self.iter+=1
566                         norm_f_new = util.Lsup(f_new)
567                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
568                         # print "   ",sub_iter," new g norm",norm_g_new
569                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
570                         #
571                         #   can we expect reduction of g?
572                         #
573                         # u_last=u_new
574                         sub_iter+=1
575                      self.relaxation=s
576                #
577                #  check for convergence:
578                #
579                norm_u_new = util.Lsup(u_new)
580                p_new=p+self.relaxation*g_new
581                norm_p_new = util.sqrt(self.inner(p_new,p_new))
582                self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))
583    
584                if self.iter>1:
585                   dg2=g_new-g
586                   df2=f_new-f
587                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
588                   norm_df2=util.Lsup(df2)
589                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
590                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
591                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
592                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
593                       converged=True
594                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
595            self.trace("convergence after %s steps."%self.iter)
596            return u,p
597    #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
598    #      tol=1.e-2
599    #      iter=0
600    #      converged=False
601    #      u=u0*1.
602    #      p=p0*1.
603    #      while not converged and iter<iter_max:
604    #          du=self.solve_f(u,p,tol)
605    #          u-=du
606    #          norm_du=util.Lsup(du)
607    #          norm_u=util.Lsup(u)
608    #        
609    #          dp=self.relaxation*self.solve_g(u,tol)
610    #          p+=dp
611    #          norm_dp=util.sqrt(self.inner(dp,dp))
612    #          norm_p=util.sqrt(self.inner(p,p))
613    #          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)
614    #          iter+=1
615    #
616    #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
617    #      if converged:
618    #          print "convergence after %s steps."%iter
619    #      else:
620    #          raise ArithmeticError("no convergence after %s steps."%iter)
621    #
622    #      return u,p
623              
624    # vim: expandtab shiftwidth=4:

Legend:
Removed from v.146  
changed lines
  Added in v.1107

  ViewVC Help
Powered by ViewVC 1.1.26