/[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 300 - (hide annotations)
Fri Dec 2 05:47:17 2005 UTC (14 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 6298 byte(s)
runs again
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 149 from esys.escript import *
18 jgs 153 from esys.escript.linearPDEs import *
19 jgs 149 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 gross 300 bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) \
66 jgs 102 * 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 123 # mypde.getOperator().saveMM("g.mm")
86 jgs 102
87     # Solve for Approximate Solution
88 jgs 117 tm=time()
89 jgs 149 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
90 jgs 117 tm=time()-tm
91 jgs 102
92     # Report Results
93     error=Lsup(u - u_approx)/Lsup(u)
94 jgs 117 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
95 jgs 102
96     return error
97    
98     #################################################################################################################
99    
100 jgs 115 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
101 jgs 102
102 jgs 117 if prec=="":
103     prec_id=0
104     else:
105     prec_id=eval("LinearPDE.%s"%prec)
106     solver_method_id=eval("LinearPDE.%s"%solver_method)
107 jgs 102 print "Scalar solver:"
108     recDim=array([len_x0,1.,1.])
109     # Define Computational Domain
110     numElem=int((totalNumElem/(len_x0*1.))**(1./numDim))
111     elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
112     # Set Mesh
113     if (numDim == 2):
114     mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
115     l0 = len_x0, l1 = 1.)
116     totElem=elemDim[0]*elemDim[1]
117     elif (numDim == 3):
118     mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
119     l0 = len_x0, l1 = 1., l2 = 1.)
120     totElem=elemDim[0]*elemDim[1]*elemDim[2]
121    
122     print " length of domain: ",recDim[:numDim]
123     print " requested elements: ",totalNumElem
124     print " num elements: ",totElem
125    
126     # Set Mesh Descriptors
127     meshDim = mesh.getDim()
128     contfunc = ContinuousFunction(mesh)
129     func = Function(mesh)
130     x = contfunc.getX()
131    
132     # Set Boundary Mask / pdelib Template "q" Parameter Vector
133     bndryMask = Scalar(value = 0, what = contfunc)
134     for i in range(meshDim):
135 gross 300 bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) * 1.0
136 jgs 102
137     # Set True Solution / pdelib Template "r" Parameter Vector
138     u = Scalar(value = 0, what = contfunc)
139     for j in range(meshDim):
140     u += x[j] * x[j]
141    
142     # Build the pdelib System Matrix and RHS
143 jgs 104 mypde=LinearPDE(mesh)
144     mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
145 jgs 117 mypde.setSolverMethod(solver_method_id)
146 jgs 102
147     # Solve for Approximate Solution
148 jgs 117 tm=time()
149 jgs 149 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
150 jgs 117 tm=time()-tm
151 jgs 102
152     # Report Results
153     error=Lsup(u - u_approx)/Lsup(u)
154 jgs 117 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
155 jgs 102
156     return error
157    
158     #######################################################################################
159    
160    
161     print "Test is started:"
162     print "----------------"
163     error=0.
164 jgs 121 for numDim in [2, 3]:
165 jgs 123 # for totalNumElem in [51200]:
166 jgs 117 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
167     for problem in [solveScalar,solveVector]:
168 jgs 123 #for problem in [solveVector]:
169 jgs 147 error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
170 jgs 122 error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
171     #if totalNumElem*2**numDim*numDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
172     # for solver_method in [ "PCG" ]:
173     # for prec in [ "JACOBI", "ILU0" ]:
174     # error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
175 jgs 102 print "----------------"
176     print "maximum error over all tests is ",error
177     print "----------------"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26