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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26