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)

"""

from esys.escript import *
from esys.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 += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \
                * 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)

    # Solve for Approximate Solution
    tm=time()
    u_approx = mypde.getSolution(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 += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 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(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 [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
      for problem in [solveScalar,solveVector]:
         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	jgs	104	from esys.escript import *
18			from esys.linearPDEs import *
19			import esys.finley as pdelib
20	jgs	117	from time import time
21	jgs	104
22	jgs	102	from numarray import *
23
24			# these values are currently fixed:
25			len_x0=1.
26	jgs	117	alpha=10.
27			tm=0
28	jgs	102
29			#############################################################################################################3
30	jgs	115	def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
31	jgs	117
32			if prec=="":
33			prec_id=0
34			else:
35			prec_id=eval("LinearPDE.%s"%prec)
36			solver_method_id=eval("LinearPDE.%s"%solver_method)
37	jgs	102	print "Vector solver:"
38			recDim=array([len_x0,1.,1.])
39			# Define Computational Domain
40			numElem=int((totalNumElem/(len_x01.))*(1./numDim))
41			elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
42
43			# Set Mesh
44			if (numDim == 2):
45			mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
46			l0 = len_x0, l1 = 1.)
47			totElem=elemDim[0]*elemDim[1]
48			elif (numDim == 3):
49			mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
50			l0 = len_x0, l1 = 1., l2 = 1.)
51			totElem=elemDim[0]elemDim[1]elemDim[2]
52
53			print " length of domain: ",recDim[:numDim]
54			print " requested elements: ",totalNumElem
55			print " num elements: ",totElem
56			# Set Mesh Descriptors
57			meshDim = mesh.getDim()
58			contfunc = ContinuousFunction(mesh)
59			func = Function(mesh)
60			x = contfunc.getX()
61
62			# Set Boundary Mask / pdelib Template "q" Parameter Vector
63			bndryMask = Vector(value = 0, what = contfunc)
64			for i in range(meshDim):
65			bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \
66			* ones((numDim,))
67
68			# Set True Solution / pdelib Template "r" Parameter Vector
69			u = Vector(value = 0, what = contfunc)
70			for i in range(meshDim):
71			for j in range(meshDim - 1):
72			u[i] += x[(i + j + 1) % meshDim]**2
73			# Set pdelib Template "A" Parameter Tensor
74			A = Tensor4(value = 0, what = func)
75			for i in range(meshDim):
76			for j in range(meshDim):
77			A[i,i,j,j] += 1.
78			A[i,j,j,i] += alpha
79			A[i,j,i,j] += alpha
80
81			# Build the pdelib System Matrix and RHS
82	jgs	104	mypde=LinearPDE(mesh)
83			mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
84	jgs	117	mypde.setSolverMethod(solver_method_id)
85	jgs	102
86			# Solve for Approximate Solution
87	jgs	117	tm=time()
88			u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
89			tm=time()-tm
90	jgs	102
91			# Report Results
92			error=Lsup(u - u_approx)/Lsup(u)
93	jgs	117	print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
94	jgs	102
95			return error
96
97			#################################################################################################################
98
99	jgs	115	def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
100	jgs	102
101	jgs	117	if prec=="":
102			prec_id=0
103			else:
104			prec_id=eval("LinearPDE.%s"%prec)
105			solver_method_id=eval("LinearPDE.%s"%solver_method)
106	jgs	102	print "Scalar solver:"
107			recDim=array([len_x0,1.,1.])
108			# Define Computational Domain
109			numElem=int((totalNumElem/(len_x01.))*(1./numDim))
110			elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
111			# Set Mesh
112			if (numDim == 2):
113			mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
114			l0 = len_x0, l1 = 1.)
115			totElem=elemDim[0]*elemDim[1]
116			elif (numDim == 3):
117			mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
118			l0 = len_x0, l1 = 1., l2 = 1.)
119			totElem=elemDim[0]elemDim[1]elemDim[2]
120
121			print " length of domain: ",recDim[:numDim]
122			print " requested elements: ",totalNumElem
123			print " num elements: ",totElem
124
125			# Set Mesh Descriptors
126			meshDim = mesh.getDim()
127			contfunc = ContinuousFunction(mesh)
128			func = Function(mesh)
129			x = contfunc.getX()
130
131			# Set Boundary Mask / pdelib Template "q" Parameter Vector
132			bndryMask = Scalar(value = 0, what = contfunc)
133			for i in range(meshDim):
134			bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 1.0
135
136			# Set True Solution / pdelib Template "r" Parameter Vector
137			u = Scalar(value = 0, what = contfunc)
138			for j in range(meshDim):
139			u += x[j] * x[j]
140
141			# Build the pdelib System Matrix and RHS
142	jgs	104	mypde=LinearPDE(mesh)
143			mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
144	jgs	117	mypde.setSolverMethod(solver_method_id)
145	jgs	102
146			# Solve for Approximate Solution
147	jgs	117	tm=time()
148			u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
149			tm=time()-tm
150	jgs	102
151			# Report Results
152			error=Lsup(u - u_approx)/Lsup(u)
153	jgs	117	print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
154	jgs	102
155			return error
156
157			#######################################################################################
158
159
160			print "Test is started:"
161			print "----------------"
162			error=0.
163	jgs	117	for numDim in [2,3]:
164			for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
165			for problem in [solveScalar,solveVector]:
166			if totalNumElem2numDimnumDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
167			for solver_method in [ "PCG" ]:
168			for prec in [ "JACOBI", "ILU0" ]:
169	jgs	115	error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
170	jgs	102	print "----------------"
171			print "maximum error over all tests is ",error
172			print "----------------"
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision