/[escript]/trunk/finley/test/python/SolveTest.py
ViewVC logotype

Annotation of /trunk/finley/test/python/SolveTest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 121 - (hide annotations)
Fri May 6 04:26:16 2005 UTC (14 years, 10 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/SolveTest.py
File MIME type: text/x-python
File size: 5986 byte(s)
Merge of development branch back to main trunk on 2005-05-06

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_x0*1.))**(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_x0*1.))**(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 121 for numDim in [2, 3]:
164 jgs 117 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
165     for problem in [solveScalar,solveVector]:
166     if totalNumElem*2**numDim*numDim< 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 "----------------"

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26