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

revision 2454 by gross, Fri May 29 03:23:25 2009 UTC revision 2455 by jfenwick, Wed Jun 3 03:29:07 2009 UTC
# Line 41  __author__="Lutz Gross, l.gross@uq.edu.a Line 41  __author__="Lutz Gross, l.gross@uq.edu.a
41
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numpy
45  import util  import util
46  import math  import math
47
# Line 223  class NoPDE: Line 223  class NoPDE:
223           @param domain: domain of the PDE           @param domain: domain of the PDE
224           @type domain: L{Domain}           @type domain: L{Domain}
225           @param D: coefficient of the solution           @param D: coefficient of the solution
226           @type D: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
227           @param Y: right hand side           @param Y: right hand side
228           @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
229           @param q: location of constraints           @param q: location of constraints
230           @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
231           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
232           @type r: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
233           """           """
234           self.__domain=domain           self.__domain=domain
235           self.__D=D           self.__D=D
# Line 258  class NoPDE: Line 258  class NoPDE:
258           Assigns values to the parameters.           Assigns values to the parameters.
259
260           @param D: coefficient of the solution           @param D: coefficient of the solution
261           @type D: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
262           @param Y: right hand side           @param Y: right hand side
263           @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
264           @param q: location of constraints           @param q: location of constraints
265           @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
266           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
267           @type r: C{float}, C{int}, C{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
268           """           """
269           if not D==None:           if not D==None:
270              self.__D=D              self.__D=D
# Line 312  class Locator: Line 312  class Locator:
312       given function space or domain which is closest to the given point x.       given function space or domain which is closest to the given point x.
313       """       """
314
315       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
316         """         """
317         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
318         or FunctionSpace for the sample point which is closest to the given         or FunctionSpace for the sample point which is closest to the given
# Line 321  class Locator: Line 321  class Locator:
321         @param where: function space         @param where: function space
322         @type where: L{escript.FunctionSpace}         @type where: L{escript.FunctionSpace}
323         @param x: coefficient of the solution         @param x: coefficient of the solution
324         @type x: C{numarray.NumArray} or C{list} of C{numarray.NumArray}         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}
325         """         """
326         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
327            self.__function_space=where            self.__function_space=where
# Line 399  class Locator: Line 399  class Locator:
399             if isinstance(id,list):             if isinstance(id,list):
400                 out=[]                 out=[]
401                 for i in id:                 for i in id:
402                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(data.getTupleForGlobalDataPoint(*i))
403                    if data.getRank()==0:                    if data.getRank()==0:
404                       out.append(o[0])                       out.append(o[0])
405                    else:                    else:
406                       out.append(o)                       out.append(o)
407                 return out                 return out
408             else:             else:
409               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(data.getTupleForGlobalDataPoint(*id))
410               if data.getRank()==0:               if data.getRank()==0:
411                  return out[0]                  return out[0]
412               else:               else:
# Line 522  def PCG(r, Aprod, x, Msolve, bilinearfor Line 522  def PCG(r, Aprod, x, Msolve, bilinearfor
522         q=Aprod(d)         q=Aprod(d)
523         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
524         x += alpha * d         x += alpha * d
525         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
526           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
527           else:
528               r += (-alpha) * q
529         rhat=Msolve(r)         rhat=Msolve(r)
530         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
531         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 645  def NewtonGMRES(defect, x, iter_max=100, Line 647  def NewtonGMRES(defect, x, iter_max=100,
647     @type defect: L{Defect}     @type defect: L{Defect}
648     @param x: initial guess for the solution, C{x} is altered.     @param x: initial guess for the solution, C{x} is altered.
649     @type x: any object type allowing basic operations such as     @type x: any object type allowing basic operations such as
650              C{numarray.NumArray}, L{Data}              C{numpy.ndarray}, L{Data}
651     @param iter_max: maximum number of iteration steps     @param iter_max: maximum number of iteration steps
652     @type iter_max: positive C{int}     @type iter_max: positive C{int}
653     @param sub_iter_max: maximum number of inner iteration steps     @param sub_iter_max: maximum number of inner iteration steps
# Line 733  def __givapp(c,s,vin): Line 735  def __givapp(c,s,vin):
735      return vrot      return vrot
736
737  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
738     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
739     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
740     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
741     g=numarray.zeros(iter_restart,numarray.Float64)     g=numpy.zeros(iter_restart,numpy.float64)
742     v=[]     v=[]
743
744     rho=defect.norm(F0)     rho=defect.norm(F0)
# Line 777  def __FDGMRES(F0, defect, x0, atol, iter Line 779  def __FDGMRES(F0, defect, x0, atol, iter
779
780          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
781          if iter > 0 :          if iter > 0 :
782              hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numpy.zeros(iter+1,numpy.float64)
783              for i in range(iter+1) : hhat[i]=h[i,iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
784              hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
785              for i in range(iter+1) : h[i,iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
# Line 800  def __FDGMRES(F0, defect, x0, atol, iter Line 802  def __FDGMRES(F0, defect, x0, atol, iter
802     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
803     # It's time to compute x and leave.     # It's time to compute x and leave.
804     if iter > 0 :     if iter > 0 :
805       y=numarray.zeros(iter,numarray.Float64)       y=numpy.zeros(iter,numpy.float64)
806       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
807       if iter > 1 :       if iter > 1 :
808          i=iter-2          i=iter-2
809          while i>=0 :          while i>=0 :
810            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
811            i=i-1            i=i-1
812       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
813       for i in range(iter-1):       for i in range(iter-1):
# Line 895  def GMRES(r, Aprod, x, bilinearform, ato Line 897  def GMRES(r, Aprod, x, bilinearform, ato
897  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
898     iter=0     iter=0
899
900     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
901     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
902     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
903     g=numarray.zeros(iter_restart+1,numarray.Float64)     g=numpy.zeros(iter_restart+1,numpy.float64)
904     v=[]     v=[]
905
906     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
# Line 966  def _GMRESm(r, Aprod, x, bilinearform, a Line 968  def _GMRESm(r, Aprod, x, bilinearform, a
968
969     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print "GMRES: iteration stopped after %s step."%iter
970     if iter > 0 :     if iter > 0 :
971       y=numarray.zeros(iter,numarray.Float64)       y=numpy.zeros(iter,numpy.float64)
972       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
973       if iter > 1 :       if iter > 1 :
974          i=iter-2          i=iter-2
975          while i>=0 :          while i>=0 :
976            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
977            i=i-1            i=i-1
978       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
979       for i in range(iter-1):       for i in range(iter-1):
# Line 1279  class ArithmeticTuple(object): Line 1281  class ArithmeticTuple(object):
1281     Example of usage::     Example of usage::
1282
1283         from esys.escript import Data         from esys.escript import Data
1284         from numarray import array         from numpy import array
1285         a=Data(...)         a=Data(...)
1286         b=array([1.,4.])         b=array([1.,4.])
1287         x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)