/[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 113 - (show annotations)
Mon Feb 28 07:06:33 2005 UTC (14 years, 1 month ago) by jgs
File MIME type: text/x-python
File size: 5526 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 from esys.escript import *
18 from esys.linearPDEs import *
19 import esys.finley as pdelib
20
21 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 mypde=LinearPDE(mesh)
76 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
77 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 mypde=LinearPDE(mesh)
130 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
131 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]:
151 # for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400]:
152 for problem in [solveScalar,solveVector]:
153 # for solver_method in [ LinearPDE.PRES20, LinearPDE.PCG, LinearPDE.DIRECT, LinearPDE.BICGSTAB]:
154 for solver_method in [ LinearPDE.PCG ]:
155 error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method),error])
156 print "----------------"
157 print "maximum error over all tests is ",error
158 print "----------------"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26