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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 147 - (show annotations)
Fri Aug 12 01:45:47 2005 UTC (14 years, 7 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 # $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 # 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 from time import time
24
25 from numarray import *
26
27 # these values are currently fixed:
28 len_x0=1.
29 alpha=10.
30 tm=0
31
32 #############################################################################################################3
33 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
34
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 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 mypde=LinearPDE(mesh)
86 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
87 mypde.setSolverMethod(solver_method_id)
88 # mypde.getOperator().saveMM("g.mm")
89
90 # Solve for Approximate Solution
91 tm=time()
92 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
93 tm=time()-tm
94
95 # Report Results
96 error=Lsup(u - u_approx)/Lsup(u)
97 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
98
99 return error
100
101 #################################################################################################################
102
103 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
104
105 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 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 mypde=LinearPDE(mesh)
147 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
148 mypde.setSolverMethod(solver_method_id)
149
150 # Solve for Approximate Solution
151 tm=time()
152 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
153 tm=time()-tm
154
155 # Report Results
156 error=Lsup(u - u_approx)/Lsup(u)
157 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
158
159 return error
160
161 #######################################################################################
162
163
164 print "Test is started:"
165 print "----------------"
166 error=0.
167 for numDim in [2, 3]:
168 # for totalNumElem in [51200]:
169 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
170 for problem in [solveScalar,solveVector]:
171 #for problem in [solveVector]:
172 error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
173 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 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