test/python/SolveTest.py

# $Id$

"""General test environment to test the solvers for scalar and vector equations 

   test parameters are 

     numDim = spatial dimension 
     totalNumElem = number of func in each direction 
     problem = solveScalar,solveVector

     solver_method = true/false 
     len_x0 = length of the domain in x0 direction (number of func in x0 is round(totalNumElem*len_x0) )
     alpha = a parameter of the PDE (not well defined yet)

"""

__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"""
from esys.escript import *
from esys.escript.linearPDEs import *
import esys.finley as pdelib
from time import time

from numarray import *

# these values are currently fixed:
len_x0=1.
alpha=10.
tm=0

#############################################################################################################3
def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
    
    if prec=="":
        prec_id=0
    else:
       prec_id=eval("LinearPDE.%s"%prec)
    solver_method_id=eval("LinearPDE.%s"%solver_method)
    print "Vector solver:"
    recDim=array([len_x0,1.,1.])
    # Define Computational Domain
    numElem=int((totalNumElem/(len_x0*1.))**(1./numDim))
    elemDim = array([int(len_x0*numElem), numElem, numElem],Int)

    # Set Mesh
    if (numDim == 2):
        mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
                         l0 = len_x0, l1 = 1.)
        totElem=elemDim[0]*elemDim[1]
    elif (numDim == 3):
        mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
                     l0 = len_x0, l1 = 1., l2 = 1.)
        totElem=elemDim[0]*elemDim[1]*elemDim[2]

    print "  length of domain: ",recDim[:numDim]
    print "  requested elements: ",totalNumElem
    print "  num elements:  ",totElem
    # Set Mesh Descriptors
    meshDim = mesh.getDim()
    contfunc = ContinuousFunction(mesh)
    func = Function(mesh)
    x = contfunc.getX()

    # Set Boundary Mask / pdelib Template "q" Parameter Vector
    bndryMask = Vector(value = 0, what = contfunc)
    for i in range(meshDim):
        bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) \
                * ones((numDim,))

    # Set True Solution / pdelib Template "r" Parameter Vector
    u = Vector(value = 0, what = contfunc)
    for i in range(meshDim):
        for j in range(meshDim - 1):
            u[i] += x[(i + j + 1) % meshDim]**2
    # Set pdelib Template "A" Parameter Tensor
    A = Tensor4(value = 0, what = func)
    for i in range(meshDim):
      for j in range(meshDim):
        A[i,i,j,j] += 1.
        A[i,j,j,i] += alpha 
        A[i,j,i,j] += alpha 

    # Build the pdelib System Matrix and RHS
    mypde=LinearPDE(mesh)
    mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
    mypde.setSolverMethod(solver_method_id)
    # mypde.getOperator().saveMM("g.mm")

    # Solve for Approximate Solution
    tm=time()
    u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
    tm=time()-tm

    # Report Results
    error=Lsup(u - u_approx)/Lsup(u)
    print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
  
    return error

#################################################################################################################

def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):

    if prec=="":
        prec_id=0
    else:
       prec_id=eval("LinearPDE.%s"%prec)
    solver_method_id=eval("LinearPDE.%s"%solver_method)
    print "Scalar solver:"
    recDim=array([len_x0,1.,1.])
    # Define Computational Domain
    numElem=int((totalNumElem/(len_x0*1.))**(1./numDim))
    elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
    # Set Mesh
    if (numDim == 2):
        mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
                         l0 = len_x0, l1 = 1.)
        totElem=elemDim[0]*elemDim[1]
    elif (numDim == 3):
        mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
                     l0 = len_x0, l1 = 1., l2 = 1.)
        totElem=elemDim[0]*elemDim[1]*elemDim[2]

    print "  length of domain: ",recDim[:numDim]
    print "  requested elements: ",totalNumElem
    print "  num elements:  ",totElem

    # Set Mesh Descriptors
    meshDim = mesh.getDim()
    contfunc = ContinuousFunction(mesh)
    func = Function(mesh)
    x = contfunc.getX()

    # Set Boundary Mask / pdelib Template "q" Parameter Vector
    bndryMask = Scalar(value = 0, what = contfunc)
    for i in range(meshDim):
        bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) * 1.0

    # Set True Solution / pdelib Template "r" Parameter Vector
    u = Scalar(value = 0, what = contfunc)
    for j in range(meshDim):
        u += x[j] * x[j]

    # Build the pdelib System Matrix and RHS
    mypde=LinearPDE(mesh)
    mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
    mypde.setSolverMethod(solver_method_id)

    # Solve for Approximate Solution
    tm=time()
    u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
    tm=time()-tm

    # Report Results
    error=Lsup(u - u_approx)/Lsup(u)
    print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
  
    return error

#######################################################################################


print "Test is started:"
print "----------------"
error=0.
for numDim in [2, 3]:
   # for totalNumElem in [51200]:
   for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
      for problem in [solveScalar,solveVector]:
      #for problem in [solveVector]:
         error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
         error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
         #if totalNumElem*2**numDim*numDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
         # for solver_method in [ "PCG" ]:
         #    for prec in [ "JACOBI", "ILU0" ]:
         #       error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
print "----------------"
print "maximum error over all tests is ",error
print "----------------"
1	jgs	102	# $Id$
2
3			"""General test environment to test the solvers for scalar and vector equations
4
5			test parameters are
6
7			numDim = spatial dimension
8			totalNumElem = number of func in each direction
9			problem = solveScalar,solveVector
10
11			solver_method = true/false
12			len_x0 = length of the domain in x0 direction (number of func in x0 is round(totalNumElem*len_x0) )
13			alpha = a parameter of the PDE (not well defined yet)
14
15			"""
16
17	elspeth	617	__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
18			http://www.access.edu.au
19			Primary Business: Queensland, Australia"""
20			__license__="""Licensed under the Open Software License version 3.0
21			http://www.opensource.org/licenses/osl-3.0.php"""
22	jgs	149	from esys.escript import *
23	jgs	153	from esys.escript.linearPDEs import *
24	jgs	149	import esys.finley as pdelib
25	jgs	117	from time import time
26	jgs	104
27	jgs	102	from numarray import *
28
29			# these values are currently fixed:
30			len_x0=1.
31	jgs	117	alpha=10.
32			tm=0
33	jgs	102
34			#############################################################################################################3
35	jgs	115	def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
36	jgs	117
37			if prec=="":
38			prec_id=0
39			else:
40			prec_id=eval("LinearPDE.%s"%prec)
41			solver_method_id=eval("LinearPDE.%s"%solver_method)
42	jgs	102	print "Vector solver:"
43			recDim=array([len_x0,1.,1.])
44			# Define Computational Domain
45			numElem=int((totalNumElem/(len_x01.))*(1./numDim))
46			elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
47
48			# Set Mesh
49			if (numDim == 2):
50			mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
51			l0 = len_x0, l1 = 1.)
52			totElem=elemDim[0]*elemDim[1]
53			elif (numDim == 3):
54			mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
55			l0 = len_x0, l1 = 1., l2 = 1.)
56			totElem=elemDim[0]elemDim[1]elemDim[2]
57
58			print " length of domain: ",recDim[:numDim]
59			print " requested elements: ",totalNumElem
60			print " num elements: ",totElem
61			# Set Mesh Descriptors
62			meshDim = mesh.getDim()
63			contfunc = ContinuousFunction(mesh)
64			func = Function(mesh)
65			x = contfunc.getX()
66
67			# Set Boundary Mask / pdelib Template "q" Parameter Vector
68			bndryMask = Vector(value = 0, what = contfunc)
69			for i in range(meshDim):
70	gross	300	bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) \
71	jgs	102	* ones((numDim,))
72
73			# Set True Solution / pdelib Template "r" Parameter Vector
74			u = Vector(value = 0, what = contfunc)
75			for i in range(meshDim):
76			for j in range(meshDim - 1):
77			u[i] += x[(i + j + 1) % meshDim]**2
78			# Set pdelib Template "A" Parameter Tensor
79			A = Tensor4(value = 0, what = func)
80			for i in range(meshDim):
81			for j in range(meshDim):
82			A[i,i,j,j] += 1.
83			A[i,j,j,i] += alpha
84			A[i,j,i,j] += alpha
85
86			# Build the pdelib System Matrix and RHS
87	jgs	104	mypde=LinearPDE(mesh)
88			mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
89	jgs	117	mypde.setSolverMethod(solver_method_id)
90	jgs	123	# mypde.getOperator().saveMM("g.mm")
91	jgs	102
92			# Solve for Approximate Solution
93	jgs	117	tm=time()
94	jgs	149	u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
95	jgs	117	tm=time()-tm
96	jgs	102
97			# Report Results
98			error=Lsup(u - u_approx)/Lsup(u)
99	jgs	117	print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
100	jgs	102
101			return error
102
103			#################################################################################################################
104
105	jgs	115	def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
106	jgs	102
107	jgs	117	if prec=="":
108			prec_id=0
109			else:
110			prec_id=eval("LinearPDE.%s"%prec)
111			solver_method_id=eval("LinearPDE.%s"%solver_method)
112	jgs	102	print "Scalar solver:"
113			recDim=array([len_x0,1.,1.])
114			# Define Computational Domain
115			numElem=int((totalNumElem/(len_x01.))*(1./numDim))
116			elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
117			# Set Mesh
118			if (numDim == 2):
119			mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
120			l0 = len_x0, l1 = 1.)
121			totElem=elemDim[0]*elemDim[1]
122			elif (numDim == 3):
123			mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
124			l0 = len_x0, l1 = 1., l2 = 1.)
125			totElem=elemDim[0]elemDim[1]elemDim[2]
126
127			print " length of domain: ",recDim[:numDim]
128			print " requested elements: ",totalNumElem
129			print " num elements: ",totElem
130
131			# Set Mesh Descriptors
132			meshDim = mesh.getDim()
133			contfunc = ContinuousFunction(mesh)
134			func = Function(mesh)
135			x = contfunc.getX()
136
137			# Set Boundary Mask / pdelib Template "q" Parameter Vector
138			bndryMask = Scalar(value = 0, what = contfunc)
139			for i in range(meshDim):
140	gross	300	bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) * 1.0
141	jgs	102
142			# Set True Solution / pdelib Template "r" Parameter Vector
143			u = Scalar(value = 0, what = contfunc)
144			for j in range(meshDim):
145			u += x[j] * x[j]
146
147			# Build the pdelib System Matrix and RHS
148	jgs	104	mypde=LinearPDE(mesh)
149			mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
150	jgs	117	mypde.setSolverMethod(solver_method_id)
151	jgs	102
152			# Solve for Approximate Solution
153	jgs	117	tm=time()
154	jgs	149	u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
155	jgs	117	tm=time()-tm
156	jgs	102
157			# Report Results
158			error=Lsup(u - u_approx)/Lsup(u)
159	jgs	117	print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
160	jgs	102
161			return error
162
163			#######################################################################################
164
165
166			print "Test is started:"
167			print "----------------"
168			error=0.
169	jgs	121	for numDim in [2, 3]:
170	jgs	123	# for totalNumElem in [51200]:
171	jgs	117	for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
172			for problem in [solveScalar,solveVector]:
173	jgs	123	#for problem in [solveVector]:
174	jgs	147	error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
175	jgs	122	error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
176			#if totalNumElem2numDimnumDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
177			# for solver_method in [ "PCG" ]:
178			# for prec in [ "JACOBI", "ILU0" ]:
179			# error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
180	jgs	102	print "----------------"
181			print "maximum error over all tests is ",error
182			print "----------------"
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision