escript/py_src/pdetools.py

# $Id$

"""
Provides some tools related to PDEs. 

Currently includes:
    - Projector - to project a discontinuous
    - Locator - to trace values in data objects at a certain location
    - TimeIntegrationManager - to handel extraplotion in time
    - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__author__="Lutz Gross, l.gross@uq.edu.au"
__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.iservo.edu.au/esys"
__version__="$Revision$"
__date__="$Date$"


import escript
import linearPDEs
import numarray
import util

class TimeIntegrationManager:
  """
  a simple mechanism to manage time dependend values. 

  typical usage is::

     dt=0.1 # time increment
     tm=TimeIntegrationManager(inital_value,p=1)
     while t<1.
         v_guess=tm.extrapolate(dt) # extrapolate to t+dt
         v=...
         tm.checkin(dt,v)
         t+=dt

  @note: currently only p=1 is supported.
  """
  def __init__(self,*inital_values,**kwargs):
     """
     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
     """
     if kwargs.has_key("p"):
            self.__p=kwargs["p"]
     else:
            self.__p=1
     if kwargs.has_key("time"):
            self.__t=kwargs["time"]
     else:
            self.__t=0.
     self.__v_mem=[inital_values]
     self.__order=0
     self.__dt_mem=[]
     self.__num_val=len(inital_values)

  def getTime(self):
      return self.__t
  def getValue(self):
      out=self.__v_mem[0]
      if len(out)==1:
          return out[0]
      else:
          return out

  def checkin(self,dt,*values):
      """
      adds new values to the manager. the p+1 last value get lost
      """
      o=min(self.__order+1,self.__p)
      self.__order=min(self.__order+1,self.__p)
      v_mem_new=[values]
      dt_mem_new=[dt]
      for i in range(o-1):
         v_mem_new.append(self.__v_mem[i])
         dt_mem_new.append(self.__dt_mem[i])
      v_mem_new.append(self.__v_mem[o-1])
      self.__order=o
      self.__v_mem=v_mem_new
      self.__dt_mem=dt_mem_new
      self.__t+=dt

  def extrapolate(self,dt):
      """
      extrapolates to dt forward in time.
      """
      if self.__order==0:
         out=self.__v_mem[0]
      else:
        out=[]
        for i in range(self.__num_val):
           out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])

      if len(out)==0:
         return None
      elif len(out)==1:
         return out[0]
      else:
         return out
 

class Projector:
  """
  The Projector is a factory which projects a discontiuous function onto a
  continuous function on the a given domain.
  """
  def __init__(self, domain, reduce = True, fast=True):
    """
    Create a continuous function space projector for a domain.

    @param domain: Domain of the projection.
    @param reduce: Flag to reduce projection order (default is True)
    @param fast: Flag to use a fast method based on matrix lumping (default is true)
    """
    self.__pde = linearPDEs.LinearPDE(domain)
    if fast:
      self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
    self.__pde.setSymmetryOn()
    self.__pde.setReducedOrderTo(reduce)
    self.__pde.setValue(D = 1.)
    return

  def __call__(self, input_data):
    """
    Projects input_data onto a continuous function

    @param input_data: The input_data to be projected.
    """
    out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
    if input_data.getRank()==0:
        self.__pde.setValue(Y = input_data)
        out=self.__pde.getSolution()
    elif input_data.getRank()==1:
        for i0 in range(input_data.getShape()[0]):
           self.__pde.setValue(Y = input_data[i0])
           out[i0]=self.__pde.getSolution()
    elif input_data.getRank()==2:
        for i0 in range(input_data.getShape()[0]):
           for i1 in range(input_data.getShape()[1]):
              self.__pde.setValue(Y = input_data[i0,i1])
              out[i0,i1]=self.__pde.getSolution()
    elif input_data.getRank()==3:
        for i0 in range(input_data.getShape()[0]):
           for i1 in range(input_data.getShape()[1]):
              for i2 in range(input_data.getShape()[2]):
                 self.__pde.setValue(Y = input_data[i0,i1,i2])
                 out[i0,i1,i2]=self.__pde.getSolution()
    else:
        for i0 in range(input_data.getShape()[0]):
           for i1 in range(input_data.getShape()[1]):
              for i2 in range(input_data.getShape()[2]):
                 for i3 in range(input_data.getShape()[3]):
                    self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
                    out[i0,i1,i2,i3]=self.__pde.getSolution()
    return out

class NoPDE:
     """
     solves the following problem for u:

     M{kronecker[i,j]*D[j]*u[j]=Y[i]} 

     with constraint

     M{u[j]=r[j]}  where M{q[j]>0}

     where D, Y, r and q are given functions of rank 1.

     In the case of scalars this takes the form

     M{D*u=Y} 

     with constraint

     M{u=r}  where M{q>0}

     where D, Y, r and q are given scalar functions.

     The constraint is overwriting any other condition.

     @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
            that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
            thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
     """
     def __init__(self,domain,D=None,Y=None,q=None,r=None):
         """
         initialize the problem

         @param domain: domain of the PDE.
         @type domain: L{Domain}
         @param D: coefficient of the solution. 
         @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param Y: right hand side
         @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param q: location of constraints
         @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param r: value of solution at locations of constraints
         @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         """
         self.__domain=domain
         self.__D=D
         self.__Y=Y
         self.__q=q
         self.__r=r
         self.__u=None
         self.__function_space=escript.Solution(self.__domain)
     def setReducedOn(self):
         """
         sets the L{FunctionSpace} of the solution to L{ReducedSolution}
         """
         self.__function_space=escript.ReducedSolution(self.__domain)
         self.__u=None

     def setReducedOff(self):
         """
         sets the L{FunctionSpace} of the solution to L{Solution}
         """
         self.__function_space=escript.Solution(self.__domain)
         self.__u=None
         
     def setValue(self,D=None,Y=None,q=None,r=None):
         """
         assigns values to the parameters.

         @param D: coefficient of the solution. 
         @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param Y: right hand side
         @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param q: location of constraints
         @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param r: value of solution at locations of constraints
         @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         """
         if not D==None:
            self.__D=D
            self.__u=None
         if not Y==None:
            self.__Y=Y
            self.__u=None
         if not q==None:
            self.__q=q
            self.__u=None
         if not r==None:
            self.__r=r
            self.__u=None

     def getSolution(self):
         """
         returns the solution
        
         @return: the solution of the problem
         @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
         """
         if self.__u==None:
            if self.__D==None:
               raise ValueError,"coefficient D is undefined"
            D=escript.Data(self.__D,self.__function_space)
            if D.getRank()>1:
               raise ValueError,"coefficient D must have rank 0 or 1"
            if self.__Y==None:
               self.__u=escript.Data(0.,D.getShape(),self.__function_space)
            else:
               self.__u=util.quotient(self.__Y,D)
            if not self.__q==None:
                q=util.wherePositive(escript.Data(self.__q,self.__function_space))
                self.__u*=(1.-q)
                if not self.__r==None: self.__u+=q*self.__r
         return self.__u
             
class Locator:
     """
     Locator provides access to the values of data objects at a given
     spatial coordinate x.  
     
     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. 
     """

     def __init__(self,where,x=numarray.zeros((3,))):
       """
       Initializes a Locator to access values in Data objects on the Doamin 
       or FunctionSpace where for the sample point which
       closest to the given point x.

       @param where: function space
       @type where: L{escript.FunctionSpace} 
       @param x: coefficient of the solution. 
       @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
       """
       if isinstance(where,escript.FunctionSpace):
          self.__function_space=where
       else:
          self.__function_space=escript.ContinuousFunction(where)
       if isinstance(x, list):
           self.__id=[]
           for p in x:
              self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
       else:
           self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()

     def __str__(self):
       """
       Returns the coordinates of the Locator as a string.
       """
       x=self.getX()
       if instance(x,list):
          out="["
          first=True
          for xx in x:
            if not first:
                out+=","
            else:
                first=False
            out+=str(xx)
          out+="]>"
       else:
          out=str(x)
       return out

     def getX(self):
        """
    Returns the exact coordinates of the Locator.
    """
        return self(self.getFunctionSpace().getX())

     def getFunctionSpace(self):
        """
    Returns the function space of the Locator.
    """
        return self.__function_space

     def getId(self,item=None):
        """
    Returns the identifier of the location.
    """
        if item == None:
           return self.__id
        else:
           if isinstance(self.__id,list):
              return self.__id[item]
           else:
              return self.__id


     def __call__(self,data):
        """
    Returns the value of data at the Locator of a Data object otherwise 
    the object is returned.
    """
        return self.getValue(data)

     def getValue(self,data):
        """
    Returns the value of data at the Locator if data is a Data object 
    otherwise the object is returned.
    """
        if isinstance(data,escript.Data):
           if data.getFunctionSpace()==self.getFunctionSpace():
             dat=data
           else:
             dat=data.interpolate(self.getFunctionSpace())
           id=self.getId()
           r=data.getRank()
           if isinstance(id,list):
               out=[]
               for i in id:
                  o=data.getValueOfGlobalDataPoint(*i)
                  if data.getRank()==0:
                     out.append(o[0])
                  else:
                     out.append(o)
               return out
           else:
             out=data.getValueOfGlobalDataPoint(*id)
             if data.getRank()==0:
                return out[0]
             else:
                return out
        else:
           return data

class SaddlePointProblem(object):
   """
   This implements a solver for a saddlepoint problem

   M{f(u,p)=0}
   M{g(u)=0}

   for u and p. The problem is solved with an inexact Uszawa scheme for p:

   M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
   M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}

   where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of
   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'
   Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
   in fact the role of a preconditioner.
   """
   def __init__(self,verbose=False,*args):
       """
       initializes the problem

       @param verbose: switches on the printing out some information
       @type verbose: C{bool}
       @note: this method may be overwritten by a particular saddle point problem
       """
       if not isinstance(verbose,bool):
            raise TypeError("verbose needs to be of type bool.")
       self.__verbose=verbose
       self.relaxation=1.

   def trace(self,text):
       """
       prints text if verbose has been set

       @param text: a text message
       @type text: C{str}
       """
       if self.__verbose: print "%s: %s"%(str(self),text)

   def solve_f(self,u,p,tol=1.e-8):
       """
       solves 

       A_f du = f(u,p) 

       with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.

       @param u: current approximation of u
       @type u: L{escript.Data}
       @param p: current approximation of p
       @type p: L{escript.Data}
       @param tol: tolerance expected for du
       @type tol: C{float}
       @return: increment du
       @rtype: L{escript.Data}
       @note: this method has to be overwritten by a particular saddle point problem
       """
       pass

   def solve_g(self,u,tol=1.e-8):
       """
       solves 

       Q_g dp = g(u) 

       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.

       @param u: current approximation of u
       @type u: L{escript.Data}
       @param tol: tolerance expected for dp
       @type tol: C{float}
       @return: increment dp
       @rtype: L{escript.Data}
       @note: this method has to be overwritten by a particular saddle point problem
       """
       pass

   def inner(self,p0,p1):
       """
       inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
       @return: inner product of p0 and p1
       @rtype: C{float}
       """
       pass

   subiter_max=3
   def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
        """
        runs the solver

        @param u0: initial guess for C{u}
        @type u0: L{esys.escript.Data}
        @param p0: initial guess for C{p}
        @type p0: L{esys.escript.Data}
        @param tolerance: tolerance for relative error in C{u} and C{p}
        @type tolerance: positive C{float}
        @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
        @type tolerance_u: positive C{float}
        @param iter_max: maximum number of iteration steps.
        @type iter_max: C{int}
        @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the 
                                   relaxation factor. If C{accepted_reduction=None} no backtracking is used.
        @type accepted_reduction: positive C{float} or C{None}
        @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
        @type relaxation: C{float} or C{None}
        """
        tol=1.e-2
        if tolerance_u==None: tolerance_u=tolerance
        if not relaxation==None: self.relaxation=relaxation
        if accepted_reduction ==None:
              angle_limit=0.
        elif accepted_reduction>=1.:
              angle_limit=0.
        else:
              angle_limit=util.sqrt(1-accepted_reduction**2)
        self.iter=0
        u=u0
        p=p0
        #
        #   initialize things:
        #
        converged=False
        #
        #  start loop:
        #
        #  initial search direction is g
        #
        while not converged :
            if self.iter>iter_max:
                raise ArithmeticError("no convergence after %s steps."%self.iter)
            f_new=self.solve_f(u,p,tol)
            norm_f_new = util.Lsup(f_new)
            u_new=u-f_new
            g_new=self.solve_g(u_new,tol)
            self.iter+=1
            norm_g_new = util.sqrt(self.inner(g_new,g_new))
            if norm_f_new==0. and norm_g_new==0.: return u, p
            if self.iter>1 and not accepted_reduction==None:
               #
               #   did we manage to reduce the norm of G? I
               #   if not we start a backtracking procedure
               #
               # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
               if norm_g_new > accepted_reduction * norm_g:
                  sub_iter=0
                  s=self.relaxation
                  d=g
                  g_last=g
                  self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
                  while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
                     dg= g_new-g_last
                     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
                     rad=self.inner(g_new,dg)/self.relaxation
                     # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
                     # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
                     if abs(rad) < norm_dg*norm_g_new * angle_limit:
                         if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
                         break
                     r=self.relaxation
                     self.relaxation= - rad/norm_dg**2
                     s+=self.relaxation
                     #####
                     # a=g_new+self.relaxation*dg/r
                     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
                     #####
                     g_last=g_new
                     p+=self.relaxation*d
                     f_new=self.solve_f(u,p,tol)
                     u_new=u-f_new
                     g_new=self.solve_g(u_new,tol)
                     self.iter+=1
                     norm_f_new = util.Lsup(f_new)
                     norm_g_new = util.sqrt(self.inner(g_new,g_new))
                     # print "   ",sub_iter," new g norm",norm_g_new
                     self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
                     #
                     #   can we expect reduction of g?
                     #
                     # u_last=u_new
                     sub_iter+=1
                  self.relaxation=s
            #
            #  check for convergence:
            #
            norm_u_new = util.Lsup(u_new)
            p_new=p+self.relaxation*g_new
            norm_p_new = util.sqrt(self.inner(p_new,p_new))
            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))

            if self.iter>1:
               dg2=g_new-g
               df2=f_new-f
               norm_dg2=util.sqrt(self.inner(dg2,dg2))
               norm_df2=util.Lsup(df2)
               # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
               tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
               tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
               if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
                   converged=True
            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
        self.trace("convergence after %s steps."%self.iter)
        return u,p
#   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
#      tol=1.e-2
#      iter=0
#      converged=False
#      u=u0*1.
#      p=p0*1.
#      while not converged and iter<iter_max:
#          du=self.solve_f(u,p,tol)
#          u-=du
#          norm_du=util.Lsup(du)
#          norm_u=util.Lsup(u)
#        
#          dp=self.relaxation*self.solve_g(u,tol)
#          p+=dp
#          norm_dp=util.sqrt(self.inner(dp,dp))
#          norm_p=util.sqrt(self.inner(p,p))
#          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)
#          iter+=1
#
#          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
#      if converged:
#          print "convergence after %s steps."%iter
#      else:
#          raise ArithmeticError("no convergence after %s steps."%iter)
#
#      return u,p
          
# vim: expandtab shiftwidth=4:
1	jgs	121	# $Id$
2
3			"""
4	jgs	149	Provides some tools related to PDEs.
5	jgs	121
6	jgs	149	Currently includes:
7			- Projector - to project a discontinuous
8	gross	351	- Locator - to trace values in data objects at a certain location
9			- TimeIntegrationManager - to handel extraplotion in time
10	gross	867	- SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
11	gross	637
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	jgs	121	"""
19
20	gross	637	__author__="Lutz Gross, l.gross@uq.edu.au"
21	elspeth	609	__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
22			http://www.access.edu.au
23			Primary Business: Queensland, Australia"""
24	elspeth	614	__license__="""Licensed under the Open Software License version 3.0
25			http://www.opensource.org/licenses/osl-3.0.php"""
26	gross	637	__url__="http://www.iservo.edu.au/esys"
27			__version__="$Revision$"
28			__date__="$Date$"
29	elspeth	609
30	gross	637
31	jgs	149	import escript
32			import linearPDEs
33	jgs	121	import numarray
34	jgs	149	import util
35	jgs	121
36	gross	351	class TimeIntegrationManager:
37			"""
38			a simple mechanism to manage time dependend values.
39
40	gross	720	typical usage is::
41	gross	351
42	gross	720	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	gross	351
50	gross	720	@note: currently only p=1 is supported.
51	gross	351	"""
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	gross	396	def getValue(self):
72	gross	409	out=self.__v_mem[0]
73			if len(out)==1:
74			return out[0]
75			else:
76			return out
77
78	gross	351	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	gross	396
113	gross	351
114	jgs	121	class Projector:
115	jgs	149	"""
116			The Projector is a factory which projects a discontiuous function onto a
117			continuous function on the a given domain.
118			"""
119	jgs	121	def __init__(self, domain, reduce = True, fast=True):
120			"""
121	jgs	149	Create a continuous function space projector for a domain.
122	jgs	121
123	jgs	149	@param domain: Domain of the projection.
124			@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)
126	jgs	121	"""
127	jgs	149	self.__pde = linearPDEs.LinearPDE(domain)
128	jgs	148	if fast:
129	jgs	149	self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
130	jgs	121	self.__pde.setSymmetryOn()
131			self.__pde.setReducedOrderTo(reduce)
132			self.__pde.setValue(D = 1.)
133			return
134
135			def __call__(self, input_data):
136			"""
137	jgs	149	Projects input_data onto a continuous function
138	jgs	121
139	jgs	149	@param input_data: The input_data to be projected.
140	jgs	121	"""
141	gross	525	out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
142	jgs	121	if input_data.getRank()==0:
143			self.__pde.setValue(Y = input_data)
144			out=self.__pde.getSolution()
145			elif input_data.getRank()==1:
146			for i0 in range(input_data.getShape()[0]):
147			self.__pde.setValue(Y = input_data[i0])
148			out[i0]=self.__pde.getSolution()
149			elif input_data.getRank()==2:
150			for i0 in range(input_data.getShape()[0]):
151			for i1 in range(input_data.getShape()[1]):
152			self.__pde.setValue(Y = input_data[i0,i1])
153			out[i0,i1]=self.__pde.getSolution()
154			elif input_data.getRank()==3:
155			for i0 in range(input_data.getShape()[0]):
156			for i1 in range(input_data.getShape()[1]):
157			for i2 in range(input_data.getShape()[2]):
158			self.__pde.setValue(Y = input_data[i0,i1,i2])
159			out[i0,i1,i2]=self.__pde.getSolution()
160			else:
161			for i0 in range(input_data.getShape()[0]):
162			for i1 in range(input_data.getShape()[1]):
163			for i2 in range(input_data.getShape()[2]):
164			for i3 in range(input_data.getShape()[3]):
165			self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
166			out[i0,i1,i2,i3]=self.__pde.getSolution()
167			return out
168
169	gross	525	class NoPDE:
170			"""
171			solves the following problem for u:
172	jgs	121
173	gross	525	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	gross	720	@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	gross	525	"""
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	gross	720	@type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
205	gross	525	@param Y: right hand side
206	gross	720	@type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207	gross	525	@param q: location of constraints
208	gross	720	@type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209	gross	525	@param r: value of solution at locations of constraints
210	gross	720	@type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211	gross	525	"""
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	gross	720	@type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
239	gross	525	@param Y: right hand side
240	gross	720	@type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241	gross	525	@param q: location of constraints
242	gross	720	@type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243	gross	525	@param r: value of solution at locations of constraints
244	gross	720	@type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245	gross	525	"""
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	jgs	147	class Locator:
283			"""
284	jgs	149	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	jgs	147	"""
291
292			def __init__(self,where,x=numarray.zeros((3,))):
293	jgs	149	"""
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	gross	880
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	jgs	149	"""
303			if isinstance(where,escript.FunctionSpace):
304	jgs	147	self.__function_space=where
305	jgs	121	else:
306	jgs	149	self.__function_space=escript.ContinuousFunction(where)
307	gross	880	if isinstance(x, list):
308			self.__id=[]
309			for p in x:
310	gross	921	self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
311	gross	880	else:
312	gross	921	self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
313	jgs	121
314	jgs	147	def __str__(self):
315	jgs	149	"""
316			Returns the coordinates of the Locator as a string.
317			"""
318	gross	880	x=self.getX()
319			if instance(x,list):
320			out="["
321			first=True
322			for xx in x:
323			if not first:
324			out+=","
325			else:
326			first=False
327			out+=str(xx)
328			out+="]>"
329			else:
330			out=str(x)
331			return out
332	jgs	121
333	gross	880	def getX(self):
334			"""
335			Returns the exact coordinates of the Locator.
336			"""
337			return self(self.getFunctionSpace().getX())
338
339	jgs	147	def getFunctionSpace(self):
340	jgs	149	"""
341			Returns the function space of the Locator.
342			"""
343	jgs	147	return self.__function_space
344
345	gross	880	def getId(self,item=None):
346	jgs	149	"""
347			Returns the identifier of the location.
348			"""
349	gross	880	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	jgs	121
357
358	jgs	147	def __call__(self,data):
359	jgs	149	"""
360			Returns the value of data at the Locator of a Data object otherwise
361			the object is returned.
362			"""
363	jgs	147	return self.getValue(data)
364	jgs	121
365	jgs	147	def getValue(self,data):
366	jgs	149	"""
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	jgs	147	if data.getFunctionSpace()==self.getFunctionSpace():
372	gross	880	dat=data
373	jgs	147	else:
374	gross	880	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	gross	921	o=data.getValueOfGlobalDataPoint(*i)
381	gross	880	if data.getRank()==0:
382			out.append(o[0])
383			else:
384			out.append(o)
385			return out
386	jgs	147	else:
387	gross	921	out=data.getValueOfGlobalDataPoint(*id)
388	gross	880	if data.getRank()==0:
389			return out[0]
390			else:
391			return out
392	jgs	147	else:
393			return data
394	jgs	149
395	gross	867	class SaddlePointProblem(object):
396			"""
397			This implements a solver for a saddlepoint problem
398
399	gross	877	M{f(u,p)=0}
400			M{g(u)=0}
401	gross	867
402			for u and p. The problem is solved with an inexact Uszawa scheme for p:
403
404	ksteube	990	M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
405	gross	877	M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
406	gross	867
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	ksteube	990	@param verbose: switches on the printing out some information
417	gross	867	@type verbose: C{bool}
418			@note: this method may be overwritten by a particular saddle point problem
419			"""
420	gross	1107	if not isinstance(verbose,bool):
421			raise TypeError("verbose needs to be of type bool.")
422	gross	1106	self.__verbose=verbose
423	gross	877	self.relaxation=1.
424	gross	867
425			def trace(self,text):
426			"""
427			prints text if verbose has been set
428
429	ksteube	990	@param text: a text message
430	gross	867	@type text: C{str}
431			"""
432			if self.__verbose: print "%s: %s"%(str(self),text)
433
434	gross	873	def solve_f(self,u,p,tol=1.e-8):
435	gross	867	"""
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	gross	873	@param tol: tolerance expected for du
447	gross	867	@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	gross	873	def solve_g(self,u,tol=1.e-8):
455	gross	867	"""
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	gross	873	@param tol: tolerance expected for dp
465			@type tol: C{float}
466	gross	867	@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	gross	877	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	gross	873
485	gross	877	@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_dgnorm_g_new angle_limit
549			# print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
550			if abs(rad) < norm_dgnorm_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	gross	873
584	gross	877	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=tolerancenorm_dg2/(norm_gabs(self.relaxation))*norm_p_new
591			tol_eq_f=tolerance_unorm_df2/norm_fnorm_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 <= tolerancenorm_u) and (norm_dp <= tolerancenorm_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	gross	873
624	jgs	149	# vim: expandtab shiftwidth=4:
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision