/[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 153 - (show annotations)
Tue Oct 25 01:51:20 2005 UTC (13 years, 5 months ago) by jgs
File MIME type: text/x-python
File size: 6306 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

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.escript.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 # mypde.getOperator().saveMM("g.mm")
86
87 # Solve for Approximate Solution
88 tm=time()
89 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
90 tm=time()-tm
91
92 # Report Results
93 error=Lsup(u - u_approx)/Lsup(u)
94 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
95
96 return error
97
98 #################################################################################################################
99
100 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
101
102 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 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 bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 1.0
136
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 mypde=LinearPDE(mesh)
144 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
145 mypde.setSolverMethod(solver_method_id)
146
147 # Solve for Approximate Solution
148 tm=time()
149 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
150 tm=time()-tm
151
152 # Report Results
153 error=Lsup(u - u_approx)/Lsup(u)
154 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
155
156 return error
157
158 #######################################################################################
159
160
161 print "Test is started:"
162 print "----------------"
163 error=0.
164 for numDim in [2, 3]:
165 # for totalNumElem in [51200]:
166 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
167 for problem in [solveScalar,solveVector]:
168 #for problem in [solveVector]:
169 error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
170 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 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