/[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 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC revision 873 by gross, Mon Oct 16 04:07:33 2006 UTC
# Line 7  Currently includes: Line 7  Currently includes:
7      - Projector - to project a discontinuous      - Projector - to project a discontinuous
8      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
9      - TimeIntegrationManager - to handel extraplotion in time      - 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  @var __author__: name of author
13  @var __copyright__: copyrights  @var __copyright__: copyrights
# Line 341  class Locator: Line 342  class Locator:
342      """      """
343          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
344             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
345               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
346               #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
347             else:             else:
348               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
349               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
350             if data.getRank()==0:             if data.getRank()==0:
351                return out[0]                return out[0]
352             else:             else:
# Line 353  class Locator: Line 354  class Locator:
354          else:          else:
355             return data             return data
356    
357    class SaddlePointProblem(object):
358       """
359       This implements a solver for a saddlepoint problem
360    
361       f(u,p)=0
362       g(u)=0
363    
364       for u and p. The problem is solved with an inexact Uszawa scheme for p:
365    
366       Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})
367       Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})
368    
369       where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
370       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'
371       Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
372       in fact the role of a preconditioner.
373       """
374       def __init__(self,verbose=False,*args):
375           """
376           initializes the problem
377    
378           @parm verbose: switches on the printing out some information
379           @type verbose: C{bool}
380           @note: this method may be overwritten by a particular saddle point problem
381           """
382           self.__verbose=verbose
383    
384       def trace(self,text):
385           """
386           prints text if verbose has been set
387    
388           @parm text: a text message
389           @type text: C{str}
390           """
391           if self.__verbose: print "%s: %s"%(str(self),text)
392    
393       def solve_f(self,u,p,tol=1.e-8):
394           """
395           solves
396    
397           A_f du = f(u,p)
398    
399           with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
400    
401           @param u: current approximation of u
402           @type u: L{escript.Data}
403           @param p: current approximation of p
404           @type p: L{escript.Data}
405           @param tol: tolerance expected for du
406           @type tol: C{float}
407           @return: increment du
408           @rtype: L{escript.Data}
409           @note: this method has to be overwritten by a particular saddle point problem
410           """
411           pass
412    
413       def solve_g(self,u,tol=1.e-8):
414           """
415           solves
416    
417           Q_g dp = g(u)
418    
419           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.
420    
421           @param u: current approximation of u
422           @type u: L{escript.Data}
423           @param tol: tolerance expected for dp
424           @type tol: C{float}
425           @return: increment dp
426           @rtype: L{escript.Data}
427           @note: this method has to be overwritten by a particular saddle point problem
428           """
429           pass
430    
431       def inner(self,p0,p1):
432           """
433           inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
434           @return: inner product of p0 and p1
435           @rtype: C{float}
436           """
437           pass
438    
439       def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,relaxation=1.):
440          tol=1.e-2
441          iter=0
442          converged=False
443          u=u0*1.
444          p=p0*1.
445          while not converged and iter<iter_max:
446              du=self.solve_f(u,p,tol)
447              u-=du
448              norm_du=util.Lsup(du)
449              norm_u=util.Lsup(u)
450            
451              dp=relaxation*self.solve_g(u,tol)
452              p+=dp
453              norm_dp=util.sqrt(self.inner(dp,dp))
454              norm_p=util.sqrt(self.inner(p,p))
455              print iter,"-th step rel. errror u,p= (%s,%s)(%s,%s)"%(norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
456              iter+=1
457    
458              converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
459          if converged:
460              print "convergence after %s steps."%iter
461          else:
462              raise ArithmeticError("no convergence after %s steps."%iter)
463    
464          return u,p
465              
466  # vim: expandtab shiftwidth=4:  # vim: expandtab shiftwidth=4:

Legend:
Removed from v.782  
changed lines
  Added in v.873

  ViewVC Help
Powered by ViewVC 1.1.26