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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (show annotations)
Wed Dec 15 07:08:39 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: 6002 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 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