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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 117 - (show annotations)
Fri Apr 1 05:48:57 2005 UTC (15 years, 2 months ago) by jgs
File MIME type: text/x-python
File size: 5985 byte(s)
*** empty log message ***

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 time import time
21
22 from numarray import *
23
24 # these values are currently fixed:
25 len_x0=1.
26 alpha=10.
27 tm=0
28
29 #############################################################################################################3
30 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
31
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 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 mypde=LinearPDE(mesh)
83 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
84 mypde.setSolverMethod(solver_method_id)
85
86 # Solve for Approximate Solution
87 tm=time()
88 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
89 tm=time()-tm
90
91 # Report Results
92 error=Lsup(u - u_approx)/Lsup(u)
93 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
94
95 return error
96
97 #################################################################################################################
98
99 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
100
101 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 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 mypde=LinearPDE(mesh)
143 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
144 mypde.setSolverMethod(solver_method_id)
145
146 # Solve for Approximate Solution
147 tm=time()
148 u_approx = mypde.getSolution(preconditioner=prec_id,iter_max=10000)
149 tm=time()-tm
150
151 # Report Results
152 error=Lsup(u - u_approx)/Lsup(u)
153 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
154
155 return error
156
157 #######################################################################################
158
159
160 print "Test is started:"
161 print "----------------"
162 error=0.
163 for numDim in [2,3]:
164 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 error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
170 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