finley/py_src/finleybench.py

# $Id:$

#
#      COPYRIGHT ACcESS 2004 -  All Rights Reserved
#
#   This software is the property of ACcESS.  No part of this code
#   may be copied in any form or by any means without the expressed written
#   consent of ACcESS.  Copying, use or modification of this software
#   by any unauthorised person is illegal unless that
#   person has a software license agreement with ACcESS.
#

"""
some benchmarks for tetsing the finley solver. The idea is to develop a set of standart benchmarks

  * Laplace2Dorder1_?k
  * Laplace3Dorder2_?k

where ? is approximatively the number of unknowns in 1000.

@var __author__: name of author
@var __licence__: 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"
__licence__="contact: esys@access.uq.edu.au"
__url__="http://www.iservo.edu.au/esys/escript"
__version__="$Revision:$"
__date__="$Date:$"

from esys.escript import Lsup,whereZero,kronecker
from esys.escript.benchmark import BenchmarkProblem, Options, BenchmarkFilter
import esys.finley 
from esys.escript.linearPDEs import LinearPDE
import os

class FinleyFilter(BenchmarkFilter):
   """
   defines a filter for L{FinleyProblem} characteristics
   """
   TIME="t [sec]"
   ERROR="rel. error"


   def __init__(self,args=None):
      """
      sets up the filter

      @param args: list of value names to be filtered
      @type args: C{list} of L{TIME}, L{ERROR}
      """
      if args==None: args=[FinleyFilter.TIME,FinleyFilter.ERROR]
      super(FinleyFilter,self).__init__()
      self.__args=args

   def getResultNames(self):
       """
       return the names of the results produced when run() is called.
       
       @return: names the list of the names to be used when the results of the run() call are printed
       @rtype: C{list} of C{str}
       """
       return self.__args

   def __call__(self,result):
       """
       filters out the characteristic values
       
       @param result: characteristics rturned by a L{FinleyProblem} run
       @type result: C{dict}
       @return: filtered values
       @rtype: C{list} of C{str}
       """
       out=[]
       for a in self.__args:
         out.append(result[a])
       return out

class FinleyOptions(Options):
   """
   finley solver options to be handed over to paso

   """
   def __init__(self,solver_method=None,
                     preconditioner=None,
                     package=None,
                     tolerance=None,
                     verbose=False):
       self.strmap={
                      LinearPDE.DIRECT : "DIRECT",
                      LinearPDE.PCG:  "PCG",
                      LinearPDE.CR:  "CR",
                      LinearPDE.CGS: "CGS",
                      LinearPDE.BICGSTAB: "BICGSTAB",
                      LinearPDE.SSOR: "SSOR",
                      LinearPDE.ILU0: "ILU0",
                      LinearPDE.ILUT: "ILUT",
                      LinearPDE.JACOBI: "JACOBI",
                      LinearPDE.GMRES:  "GMRES",
                      LinearPDE.PRES20:  "PRES20",
                      LinearPDE.LUMPING:  "LUMPIMG",
                      LinearPDE.NO_REORDERING:  "NO_REORDERING",
                      LinearPDE.MINIMUM_FILL_IN:  "MINIMUM_FILL_IN",
                      LinearPDE.NESTED_DISSECTION: "NESTED_DISSECTION",
                      LinearPDE.SCSL:  "SCSL",
                      LinearPDE.MKL:  "MKL",
                      LinearPDE.UMFPACK: "UMFPACK",
                      LinearPDE.PASO:  "PASO"
                  }
       name=""
       if solver_method==None: 
             solver_method==LinearPDE.PRES20
       else:
             name+=self.strmap[solver_method]
       if preconditioner==None: 
             preconditioner==LinearPDE.JACOBI
       else:
             if not name=="": name+="+"
             name+=self.strmap[preconditioner]
       if package==None: 
             package==LinearPDE.PASO
       else:
             if not name=="": name+=" with "
             name+=self.strmap[package]
       if tolerance==None: 
             tolerance=1.e-8
       else:
             if not name=="": name+=", "
             name+="tol = %s"%tolerance
       self.solver_method=solver_method
       self.preconditioner=preconditioner
       self.tolerance=tolerance
       self.package=package
       self.verbose=verbose
       super(FinleyOptions,self).__init__(name=name)


class FinleyProblem(BenchmarkProblem):
   """
   The general benchmark problem for Finley
   """
   def run(self,options):
       """
       creates a domain and a PDE on this domain, solves it (with the given options) and returns the 
       elapsed time and the error.
       
       @param options: paso solver options
       @type options:  L{PasoOptions}
       @return:  elapsed time and the error of calculated solution
       @rtype: pair of C{float}
       """
       domain=self.getDomain()
       pde,u=self.getTestProblem(domain)
       pde.setTolerance(options.tolerance)
       pde.setSolverMethod(options.solver_method,options.preconditioner)
       pde.setSolverPackage(options.package)
       a=os.times()[4]
       uh=pde.getSolution(verbose=options.verbose)
       a=os.times()[4]-a
       if u==None:
          return {FinleyFilter.TIME : a , FinleyFilter.ERROR : None }
       else:
          error=Lsup(u-uh)/Lsup(u)
          return {FinleyFilter.TIME : a , FinleyFilter.ERROR : error }

   def getTestProblem(self,domain):
       """
       returns a PDEto be solved and an exact solution

       @param domain: the PDE domain
       @type domain: L{escript.Domain}
       @return: a linear PDE to be solved an a reference solution
       @rtype: L{LinearPDE},L{escript.Data}
       @remark: must be overwritten by a particular problem
       """
       raise NotImplementedError

   def getDomain(self):
       """
       returns the domain of the problem

       @return: a domain
       @rtype: L{escript.Domain}
       @remark: must be overwritten by a particular problem
       """
       raise NotImplementedError

class RegularFinleyProblem(FinleyProblem):
    """
    base class for finley problem on a rectangular mesh
    """
    def __init__(self,n=1,order=1,dim=2):
       """
       sets up a recangular mesh in finley on a unit cube/square
 
       @param n: number of elements in each spactial direction
       @type n: C{int}
       @param order: element order
       @type order: 1 or 2
       @param dim: spatial dimension
       @type n: 2 or 3
       """
       super(RegularFinleyProblem,self).__init__(name=str((order*n+1)**dim))
       self.__n=n
       self.__order=order
       self.__dim=dim

    def getDomain(self):
       """
       returns the unit square/cube with a rectangular mesh

       @return: a domain
       @rtype: L{escript.Domain}
       """
       if self.__dim==2:
          domain=esys.finley.Rectangle(n0=self.__n,n1=self.__n,order=self.__order)
       else:
          domain=esys.finley.Brick(n0=self.__n,n1=self.__n,n2=self.__n,order=self.__order)
       return domain

class LaplaceProblem(RegularFinleyProblem):
    """
    base class for the Lapalce eqaution on a rectangular mesh
    """
    def getTestProblem(self,domain):
         """
         returns a PDE and a test solution on the given domain
     
         @param doamin: a domain
         @type domain: L{escript.Domain}
         @return: the Laplace equation and a test solution
         @rtype: C{tuple} of C{LinearPDE} and C{escript.Data}
         """
         x=domain.getX()
         msk=whereZero(x[0])+whereZero(x[0]-1.)
         u=x[0]
         for i in range(1,domain.getDim()):
            msk+=whereZero(x[i])+whereZero(x[i]-1.)
            u*=(x[i]-i)
         pde=LinearPDE(domain)
         pde.setSymmetryOn() 
         pde.setValue(A=kronecker(domain),q=msk,r=u)
         return pde,u

class Laplace2DOrder1_30k(LaplaceProblem): 
    def __init__(self):
         super(Laplace2DOrder1_30k,self).__init__(n=176,order=1,dim=2)
class Laplace2DOrder2_30k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder2_30k,self).__init__(n=88,order=2,dim=2)
class Laplace2DOrder1_60k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder1_60k,self).__init__(n=248,order=1,dim=2)
class Laplace2DOrder2_60k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder2_60k,self).__init__(n=124,order=2,dim=2)
class Laplace2DOrder1_120k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder1_120k,self).__init__(n=349,order=1,dim=2)
class Laplace2DOrder2_120k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder2_120k,self).__init__(n=175,order=2,dim=2)
class Laplace2DOrder1_240k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder1_240k,self).__init__(n=492,order=1,dim=2)
class Laplace2DOrder2_240k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder2_240k,self).__init__(n=246,order=2,dim=2)
class Laplace2DOrder1_480k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder1_480k,self).__init__(n=694,order=1,dim=2)
class Laplace2DOrder2_480k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder2_480k,self).__init__(n=347,order=2,dim=2)
class Laplace2DOrder1_960k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder1_960k,self).__init__(n=978,order=1,dim=2)
class Laplace2DOrder2_960k(LaplaceProblem):
    def __init__(self):
         super(Laplace2DOrder2_960k,self).__init__(n=489,order=2,dim=2)
1	gross	381	# $Id:$
2
3			#
4			# COPYRIGHT ACcESS 2004 - All Rights Reserved
5			#
6			# This software is the property of ACcESS. No part of this code
7			# may be copied in any form or by any means without the expressed written
8			# consent of ACcESS. Copying, use or modification of this software
9			# by any unauthorised person is illegal unless that
10			# person has a software license agreement with ACcESS.
11			#
12
13			"""
14			some benchmarks for tetsing the finley solver. The idea is to develop a set of standart benchmarks
15
16			* Laplace2Dorder1_?k
17			* Laplace3Dorder2_?k
18
19			where ? is approximatively the number of unknowns in 1000.
20
21			@var __author__: name of author
22			@var __licence__: licence agreement
23			var __url__: url entry point on documentation
24			@var __version__: version
25			@var __date__: date of the version
26			"""
27
28			__author__="Lutz Gross, l.gross@uq.edu.au"
29			__licence__="contact: esys@access.uq.edu.au"
30			__url__="http://www.iservo.edu.au/esys/escript"
31			__version__="$Revision:$"
32			__date__="$Date:$"
33
34	gross	383	from esys.escript import Lsup,whereZero,kronecker
35	gross	381	from esys.escript.benchmark import BenchmarkProblem, Options, BenchmarkFilter
36			import esys.finley
37	gross	383	from esys.escript.linearPDEs import LinearPDE
38	gross	381	import os
39
40			class FinleyFilter(BenchmarkFilter):
41			"""
42			defines a filter for L{FinleyProblem} characteristics
43			"""
44			TIME="t [sec]"
45			ERROR="rel. error"
46
47
48			def __init__(self,args=None):
49			"""
50			sets up the filter
51
52			@param args: list of value names to be filtered
53			@type args: C{list} of L{TIME}, L{ERROR}
54			"""
55			if args==None: args=[FinleyFilter.TIME,FinleyFilter.ERROR]
56			super(FinleyFilter,self).__init__()
57			self.__args=args
58
59			def getResultNames(self):
60			"""
61			return the names of the results produced when run() is called.
62
63			@return: names the list of the names to be used when the results of the run() call are printed
64			@rtype: C{list} of C{str}
65			"""
66			return self.__args
67
68			def __call__(self,result):
69			"""
70			filters out the characteristic values
71
72			@param result: characteristics rturned by a L{FinleyProblem} run
73			@type result: C{dict}
74			@return: filtered values
75			@rtype: C{list} of C{str}
76			"""
77			out=[]
78			for a in self.__args:
79			out.append(result[a])
80			return out
81
82	gross	383	class FinleyOptions(Options):
83			"""
84			finley solver options to be handed over to paso
85	gross	381
86	gross	383	"""
87			def __init__(self,solver_method=None,
88			preconditioner=None,
89			package=None,
90			tolerance=None,
91			verbose=False):
92			self.strmap={
93			LinearPDE.DIRECT : "DIRECT",
94			LinearPDE.PCG: "PCG",
95			LinearPDE.CR: "CR",
96			LinearPDE.CGS: "CGS",
97			LinearPDE.BICGSTAB: "BICGSTAB",
98			LinearPDE.SSOR: "SSOR",
99			LinearPDE.ILU0: "ILU0",
100			LinearPDE.ILUT: "ILUT",
101			LinearPDE.JACOBI: "JACOBI",
102			LinearPDE.GMRES: "GMRES",
103			LinearPDE.PRES20: "PRES20",
104			LinearPDE.LUMPING: "LUMPIMG",
105			LinearPDE.NO_REORDERING: "NO_REORDERING",
106			LinearPDE.MINIMUM_FILL_IN: "MINIMUM_FILL_IN",
107			LinearPDE.NESTED_DISSECTION: "NESTED_DISSECTION",
108			LinearPDE.SCSL: "SCSL",
109			LinearPDE.MKL: "MKL",
110			LinearPDE.UMFPACK: "UMFPACK",
111			LinearPDE.PASO: "PASO"
112			}
113			name=""
114			if solver_method==None:
115			solver_method==LinearPDE.PRES20
116			else:
117			name+=self.strmap[solver_method]
118			if preconditioner==None:
119			preconditioner==LinearPDE.JACOBI
120			else:
121			if not name=="": name+="+"
122			name+=self.strmap[preconditioner]
123			if package==None:
124			package==LinearPDE.PASO
125			else:
126			if not name=="": name+=" with "
127			name+=self.strmap[package]
128			if tolerance==None:
129			tolerance=1.e-8
130			else:
131			if not name=="": name+=", "
132			name+="tol = %s"%tolerance
133			self.solver_method=solver_method
134			self.preconditioner=preconditioner
135			self.tolerance=tolerance
136			self.package=package
137			self.verbose=verbose
138			super(FinleyOptions,self).__init__(name=name)
139
140
141
142	gross	381	class FinleyProblem(BenchmarkProblem):
143			"""
144			The general benchmark problem for Finley
145			"""
146			def run(self,options):
147			"""
148			creates a domain and a PDE on this domain, solves it (with the given options) and returns the
149			elapsed time and the error.
150
151			@param options: paso solver options
152			@type options: L{PasoOptions}
153			@return: elapsed time and the error of calculated solution
154			@rtype: pair of C{float}
155			"""
156			domain=self.getDomain()
157			pde,u=self.getTestProblem(domain)
158	gross	383	pde.setTolerance(options.tolerance)
159			pde.setSolverMethod(options.solver_method,options.preconditioner)
160			pde.setSolverPackage(options.package)
161	gross	381	a=os.times()[4]
162	gross	383	uh=pde.getSolution(verbose=options.verbose)
163	gross	381	a=os.times()[4]-a
164			if u==None:
165	gross	383	return {FinleyFilter.TIME : a , FinleyFilter.ERROR : None }
166	gross	381	else:
167			error=Lsup(u-uh)/Lsup(u)
168	gross	383	return {FinleyFilter.TIME : a , FinleyFilter.ERROR : error }
169	gross	381
170			def getTestProblem(self,domain):
171			"""
172			returns a PDEto be solved and an exact solution
173
174			@param domain: the PDE domain
175			@type domain: L{escript.Domain}
176			@return: a linear PDE to be solved an a reference solution
177			@rtype: L{LinearPDE},L{escript.Data}
178			@remark: must be overwritten by a particular problem
179			"""
180			raise NotImplementedError
181
182			def getDomain(self):
183			"""
184			returns the domain of the problem
185
186			@return: a domain
187			@rtype: L{escript.Domain}
188			@remark: must be overwritten by a particular problem
189			"""
190			raise NotImplementedError
191
192			class RegularFinleyProblem(FinleyProblem):
193			"""
194			base class for finley problem on a rectangular mesh
195			"""
196			def __init__(self,n=1,order=1,dim=2):
197			"""
198			sets up a recangular mesh in finley on a unit cube/square
199
200			@param n: number of elements in each spactial direction
201			@type n: C{int}
202			@param order: element order
203			@type order: 1 or 2
204			@param dim: spatial dimension
205			@type n: 2 or 3
206			"""
207			super(RegularFinleyProblem,self).__init__(name=str((ordern+1)*dim))
208			self.__n=n
209			self.__order=order
210			self.__dim=dim
211
212			def getDomain(self):
213			"""
214			returns the unit square/cube with a rectangular mesh
215
216			@return: a domain
217			@rtype: L{escript.Domain}
218			"""
219	gross	383	if self.__dim==2:
220			domain=esys.finley.Rectangle(n0=self.__n,n1=self.__n,order=self.__order)
221	gross	381	else:
222	gross	383	domain=esys.finley.Brick(n0=self.__n,n1=self.__n,n2=self.__n,order=self.__order)
223	gross	381	return domain
224
225			class LaplaceProblem(RegularFinleyProblem):
226			"""
227			base class for the Lapalce eqaution on a rectangular mesh
228			"""
229			def getTestProblem(self,domain):
230			"""
231			returns a PDE and a test solution on the given domain
232
233			@param doamin: a domain
234			@type domain: L{escript.Domain}
235			@return: the Laplace equation and a test solution
236			@rtype: C{tuple} of C{LinearPDE} and C{escript.Data}
237			"""
238			x=domain.getX()
239			msk=whereZero(x[0])+whereZero(x[0]-1.)
240			u=x[0]
241			for i in range(1,domain.getDim()):
242			msk+=whereZero(x[i])+whereZero(x[i]-1.)
243			u*=(x[i]-i)
244	gross	383	pde=LinearPDE(domain)
245	gross	381	pde.setSymmetryOn()
246	gross	383	pde.setValue(A=kronecker(domain),q=msk,r=u)
247	gross	381	return pde,u
248
249			class Laplace2DOrder1_30k(LaplaceProblem):
250			def __init__(self):
251			super(Laplace2DOrder1_30k,self).__init__(n=176,order=1,dim=2)
252			class Laplace2DOrder2_30k(LaplaceProblem):
253			def __init__(self):
254			super(Laplace2DOrder2_30k,self).__init__(n=88,order=2,dim=2)
255			class Laplace2DOrder1_60k(LaplaceProblem):
256			def __init__(self):
257			super(Laplace2DOrder1_60k,self).__init__(n=248,order=1,dim=2)
258			class Laplace2DOrder2_60k(LaplaceProblem):
259			def __init__(self):
260			super(Laplace2DOrder2_60k,self).__init__(n=124,order=2,dim=2)
261			class Laplace2DOrder1_120k(LaplaceProblem):
262			def __init__(self):
263			super(Laplace2DOrder1_120k,self).__init__(n=349,order=1,dim=2)
264			class Laplace2DOrder2_120k(LaplaceProblem):
265			def __init__(self):
266			super(Laplace2DOrder2_120k,self).__init__(n=175,order=2,dim=2)
267			class Laplace2DOrder1_240k(LaplaceProblem):
268			def __init__(self):
269			super(Laplace2DOrder1_240k,self).__init__(n=492,order=1,dim=2)
270			class Laplace2DOrder2_240k(LaplaceProblem):
271			def __init__(self):
272			super(Laplace2DOrder2_240k,self).__init__(n=246,order=2,dim=2)
273			class Laplace2DOrder1_480k(LaplaceProblem):
274			def __init__(self):
275			super(Laplace2DOrder1_480k,self).__init__(n=694,order=1,dim=2)
276			class Laplace2DOrder2_480k(LaplaceProblem):
277			def __init__(self):
278			super(Laplace2DOrder2_480k,self).__init__(n=347,order=2,dim=2)
279			class Laplace2DOrder1_960k(LaplaceProblem):
280			def __init__(self):
281			super(Laplace2DOrder1_960k,self).__init__(n=978,order=1,dim=2)
282			class Laplace2DOrder2_960k(LaplaceProblem):
283			def __init__(self):
284			super(Laplace2DOrder2_960k,self).__init__(n=489,order=2,dim=2)