/[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 147 - (hide annotations)
Fri Aug 12 01:45:47 2005 UTC (14 years, 8 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/SolveTest.py
File MIME type: text/x-python
File size: 6372 byte(s)
erge of development branch dev-02 back to main trunk on 2005-08-12

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26