/[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 617 - (show annotations)
Wed Mar 22 02:58:17 2006 UTC (13 years, 8 months ago) by elspeth
File MIME type: text/x-python
File size: 6586 byte(s)
More copyright.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26