/[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 104 - (hide annotations)
Fri Dec 17 07:43:12 2004 UTC (15 years, 3 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/SolveTest.py
File MIME type: text/x-python
File size: 5495 byte(s)
*** empty log message ***

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 jgs 104 from esys.escript import *
18     from esys.linearPDEs import *
19     import esys.finley as pdelib
20    
21 jgs 102 from numarray import *
22    
23     # these values are currently fixed:
24     len_x0=1.
25     alpha=0.1
26    
27     #############################################################################################################3
28     def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method):
29    
30     print "Vector solver:"
31     recDim=array([len_x0,1.,1.])
32     # Define Computational Domain
33     numElem=int((totalNumElem/(len_x0*1.))**(1./numDim))
34     elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
35    
36     # Set Mesh
37     if (numDim == 2):
38     mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
39     l0 = len_x0, l1 = 1.)
40     totElem=elemDim[0]*elemDim[1]
41     elif (numDim == 3):
42     mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
43     l0 = len_x0, l1 = 1., l2 = 1.)
44     totElem=elemDim[0]*elemDim[1]*elemDim[2]
45    
46     print " length of domain: ",recDim[:numDim]
47     print " requested elements: ",totalNumElem
48     print " num elements: ",totElem
49     # Set Mesh Descriptors
50     meshDim = mesh.getDim()
51     contfunc = ContinuousFunction(mesh)
52     func = Function(mesh)
53     x = contfunc.getX()
54    
55     # Set Boundary Mask / pdelib Template "q" Parameter Vector
56     bndryMask = Vector(value = 0, what = contfunc)
57     for i in range(meshDim):
58     bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) \
59     * ones((numDim,))
60    
61     # Set True Solution / pdelib Template "r" Parameter Vector
62     u = Vector(value = 0, what = contfunc)
63     for i in range(meshDim):
64     for j in range(meshDim - 1):
65     u[i] += x[(i + j + 1) % meshDim]**2
66     # Set pdelib Template "A" Parameter Tensor
67     A = Tensor4(value = 0, what = func)
68     for i in range(meshDim):
69     for j in range(meshDim):
70     A[i,i,j,j] += 1.
71     A[i,j,j,i] += alpha
72     A[i,j,i,j] += alpha
73    
74     # Build the pdelib System Matrix and RHS
75 jgs 104 mypde=LinearPDE(mesh)
76     mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
77 jgs 102 mypde.setSolverMethod(solver_method)
78    
79     # Solve for Approximate Solution
80     u_approx = mypde.getSolution(iter_max=10000)
81    
82     # Report Results
83     error=Lsup(u - u_approx)/Lsup(u)
84     print " error L^sup Norm : ", error
85     print " residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx))
86    
87     return error
88    
89     #################################################################################################################
90    
91     def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method):
92    
93     print "Scalar solver:"
94     recDim=array([len_x0,1.,1.])
95     # Define Computational Domain
96     numElem=int((totalNumElem/(len_x0*1.))**(1./numDim))
97     elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
98     # Set Mesh
99     if (numDim == 2):
100     mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
101     l0 = len_x0, l1 = 1.)
102     totElem=elemDim[0]*elemDim[1]
103     elif (numDim == 3):
104     mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
105     l0 = len_x0, l1 = 1., l2 = 1.)
106     totElem=elemDim[0]*elemDim[1]*elemDim[2]
107    
108     print " length of domain: ",recDim[:numDim]
109     print " requested elements: ",totalNumElem
110     print " num elements: ",totElem
111    
112     # Set Mesh Descriptors
113     meshDim = mesh.getDim()
114     contfunc = ContinuousFunction(mesh)
115     func = Function(mesh)
116     x = contfunc.getX()
117    
118     # Set Boundary Mask / pdelib Template "q" Parameter Vector
119     bndryMask = Scalar(value = 0, what = contfunc)
120     for i in range(meshDim):
121     bndryMask += (x[i].whereZero() + (x[i]-recDim[i]).whereZero()) * 1.0
122    
123     # Set True Solution / pdelib Template "r" Parameter Vector
124     u = Scalar(value = 0, what = contfunc)
125     for j in range(meshDim):
126     u += x[j] * x[j]
127    
128     # Build the pdelib System Matrix and RHS
129 jgs 104 mypde=LinearPDE(mesh)
130     mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
131 jgs 102 mypde.setSolverMethod(solver_method)
132    
133     # Solve for Approximate Solution
134     u_approx = mypde.getSolution(iter_max=10000)
135    
136     # Report Results
137     error=Lsup(u - u_approx)/Lsup(u)
138     print " error L^sup Norm : ", error
139     print " residual L^sup Norm : ", Lsup(mypde.getResidual(u_approx))
140    
141     return error
142    
143     #######################################################################################
144    
145    
146     print "Test is started:"
147     print "----------------"
148     error=0.
149     for numDim in [2,3]:
150     for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400]:
151     for problem in [solveScalar,solveVector]:
152     # for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.DIRECT, LinearPDE.BICGSTAB]:
153     for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.BICGSTAB]:
154     error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method),error])
155     print "----------------"
156     print "maximum error over all tests is ",error
157     print "----------------"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26