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 numarray import *

# these values are currently fixed:
len_x0=1.
alpha=0.1

#############################################################################################################3
def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):

    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)

    # Solve for Approximate Solution
    u_approx = mypde.getSolution(preconditioner=prec,iter_max=10000)

    # Report Results
    error=Lsup(u - u_approx)/Lsup(u)
    print "    error L^sup Norm : ", error
    print "    residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx))
  
    return error

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

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

    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)

    # Solve for Approximate Solution
    u_approx = mypde.getSolution(preconditioner=prec,iter_max=10000)

    # Report Results
    error=Lsup(u - u_approx)/Lsup(u)
    print "    error L^sup Norm : ", error
    print "    residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx))
  
    return error

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


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