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

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 elspeth 617 __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 jgs 149 from esys.escript import *
23 jgs 153 from esys.escript.linearPDEs import *
24 jgs 149 import esys.finley as pdelib
25 jgs 117 from time import time
26 jgs 104
27 jgs 102 from numarray import *
28    
29     # these values are currently fixed:
30     len_x0=1.
31 jgs 117 alpha=10.
32     tm=0
33 jgs 102
34     #############################################################################################################3
35 jgs 115 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
36 jgs 117
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 jgs 102 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 gross 300 bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) \
71 jgs 102 * 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 jgs 104 mypde=LinearPDE(mesh)
88     mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
89 jgs 117 mypde.setSolverMethod(solver_method_id)
90 jgs 123 # mypde.getOperator().saveMM("g.mm")
91 jgs 102
92     # Solve for Approximate Solution
93 jgs 117 tm=time()
94 jgs 149 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
95 jgs 117 tm=time()-tm
96 jgs 102
97     # Report Results
98     error=Lsup(u - u_approx)/Lsup(u)
99 jgs 117 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
100 jgs 102
101     return error
102    
103     #################################################################################################################
104    
105 jgs 115 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
106 jgs 102
107 jgs 117 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 jgs 102 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 gross 300 bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) * 1.0
141 jgs 102
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 jgs 104 mypde=LinearPDE(mesh)
149     mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
150 jgs 117 mypde.setSolverMethod(solver_method_id)
151 jgs 102
152     # Solve for Approximate Solution
153 jgs 117 tm=time()
154 jgs 149 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
155 jgs 117 tm=time()-tm
156 jgs 102
157     # Report Results
158     error=Lsup(u - u_approx)/Lsup(u)
159 jgs 117 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
160 jgs 102
161     return error
162    
163     #######################################################################################
164    
165    
166     print "Test is started:"
167     print "----------------"
168     error=0.
169 jgs 121 for numDim in [2, 3]:
170 jgs 123 # for totalNumElem in [51200]:
171 jgs 117 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
172     for problem in [solveScalar,solveVector]:
173 jgs 123 #for problem in [solveVector]:
174 jgs 147 error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
175 jgs 122 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 jgs 102 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