/[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 1388 - (show annotations)
Fri Jan 11 07:45:58 2008 UTC (12 years, 5 months ago) by trankine
File MIME type: text/x-python
File size: 7006 byte(s)
And get the *(&(*&(* name right
1 #
2 # $Id$
3 #
4 #######################################################
5 #
6 # Copyright 2003-2007 by ACceSS MNRF
7 # Copyright 2007 by University of Queensland
8 #
9 # http://esscc.uq.edu.au
10 # Primary Business: Queensland, Australia
11 # Licensed under the Open Software License version 3.0
12 # http://www.opensource.org/licenses/osl-3.0.php
13 #
14 #######################################################
15 #
16
17 """General test environment to test the solvers for scalar and vector equations
18
19 test parameters are
20
21 numDim = spatial dimension
22 totalNumElem = number of func in each direction
23 problem = solveScalar,solveVector
24
25 solver_method = true/false
26 len_x0 = length of the domain in x0 direction (number of func in x0 is round(totalNumElem*len_x0) )
27 alpha = a parameter of the PDE (not well defined yet)
28
29 """
30
31 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
32 http://www.access.edu.au
33 Primary Business: Queensland, Australia"""
34 __license__="""Licensed under the Open Software License version 3.0
35 http://www.opensource.org/licenses/osl-3.0.php"""
36 from esys.escript import *
37 from esys.escript.linearPDEs import *
38 import esys.finley as pdelib
39 from time import time
40
41 from numarray import *
42
43 # these values are currently fixed:
44 len_x0=1.
45 alpha=10.
46 tm=0
47
48 #############################################################################################################3
49 def solveVector(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
50
51 if prec=="":
52 prec_id=0
53 else:
54 prec_id=eval("LinearPDE.%s"%prec)
55 solver_method_id=eval("LinearPDE.%s"%solver_method)
56 print "Vector solver:"
57 recDim=array([len_x0,1.,1.])
58 # Define Computational Domain
59 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim))
60 elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
61
62 # Set Mesh
63 if (numDim == 2):
64 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
65 l0 = len_x0, l1 = 1.)
66 totElem=elemDim[0]*elemDim[1]
67 elif (numDim == 3):
68 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
69 l0 = len_x0, l1 = 1., l2 = 1.)
70 totElem=elemDim[0]*elemDim[1]*elemDim[2]
71
72 print " length of domain: ",recDim[:numDim]
73 print " requested elements: ",totalNumElem
74 print " num elements: ",totElem
75 # Set Mesh Descriptors
76 meshDim = mesh.getDim()
77 contfunc = ContinuousFunction(mesh)
78 func = Function(mesh)
79 x = contfunc.getX()
80
81 # Set Boundary Mask / pdelib Template "q" Parameter Vector
82 bndryMask = Vector(value = 0, what = contfunc)
83 for i in range(meshDim):
84 bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) \
85 * ones((numDim,))
86
87 # Set True Solution / pdelib Template "r" Parameter Vector
88 u = Vector(value = 0, what = contfunc)
89 for i in range(meshDim):
90 for j in range(meshDim - 1):
91 u[i] += x[(i + j + 1) % meshDim]**2
92 # Set pdelib Template "A" Parameter Tensor
93 A = Tensor4(value = 0, what = func)
94 for i in range(meshDim):
95 for j in range(meshDim):
96 A[i,i,j,j] += 1.
97 A[i,j,j,i] += alpha
98 A[i,j,i,j] += alpha
99
100 # Build the pdelib System Matrix and RHS
101 mypde=LinearPDE(mesh)
102 mypde.setValue(A = A, Y = - 2 * alpha * (meshDim - 1)*ones(meshDim), q = bndryMask, r = u)
103 mypde.setSolverMethod(solver_method_id)
104 # mypde.getOperator().saveMM("g.mm")
105
106 # Solve for Approximate Solution
107 tm=time()
108 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
109 tm=time()-tm
110
111 # Report Results
112 error=Lsup(u - u_approx)/Lsup(u)
113 print "@@ Vector %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
114
115 return error
116
117 #################################################################################################################
118
119 def solveScalar(numDim, totalNumElem, len_x0, alpha, solver_method,prec):
120
121 if prec=="":
122 prec_id=0
123 else:
124 prec_id=eval("LinearPDE.%s"%prec)
125 solver_method_id=eval("LinearPDE.%s"%solver_method)
126 print "Scalar solver:"
127 recDim=array([len_x0,1.,1.])
128 # Define Computational Domain
129 numElem=int((totalNumElem/(len_x0*1.))**(1./numDim))
130 elemDim = array([int(len_x0*numElem), numElem, numElem],Int)
131 # Set Mesh
132 if (numDim == 2):
133 mesh = pdelib.Rectangle(elemDim[0], elemDim[1], 2, \
134 l0 = len_x0, l1 = 1.)
135 totElem=elemDim[0]*elemDim[1]
136 elif (numDim == 3):
137 mesh = pdelib.Brick(elemDim[0], elemDim[1], elemDim[2], 2, \
138 l0 = len_x0, l1 = 1., l2 = 1.)
139 totElem=elemDim[0]*elemDim[1]*elemDim[2]
140
141 print " length of domain: ",recDim[:numDim]
142 print " requested elements: ",totalNumElem
143 print " num elements: ",totElem
144
145 # Set Mesh Descriptors
146 meshDim = mesh.getDim()
147 contfunc = ContinuousFunction(mesh)
148 func = Function(mesh)
149 x = contfunc.getX()
150
151 # Set Boundary Mask / pdelib Template "q" Parameter Vector
152 bndryMask = Scalar(value = 0, what = contfunc)
153 for i in range(meshDim):
154 bndryMask += (whereZero(x[i]) + whereZero(x[i]-recDim[i])) * 1.0
155
156 # Set True Solution / pdelib Template "r" Parameter Vector
157 u = Scalar(value = 0, what = contfunc)
158 for j in range(meshDim):
159 u += x[j] * x[j]
160
161 # Build the pdelib System Matrix and RHS
162 mypde=LinearPDE(mesh)
163 mypde.setValue(A = identity(numDim), D = alpha, Y = alpha * u - 2 * meshDim, q = bndryMask, r = u)
164 mypde.setSolverMethod(solver_method_id)
165
166 # Solve for Approximate Solution
167 tm=time()
168 u_approx = mypde.getSolution(verbose=True,preconditioner=prec_id,iter_max=10000)
169 tm=time()-tm
170
171 # Report Results
172 error=Lsup(u - u_approx)/Lsup(u)
173 print "@@ Scalar %d : %d : %s(%s): error L^sup Norm : %e, time %e"%(mypde.getDim(),totElem,solver_method,prec,error,tm)
174
175 return error
176
177 #######################################################################################
178
179
180 print "Test is started:"
181 print "----------------"
182 error=0.
183 for numDim in [2, 3]:
184 # for totalNumElem in [51200]:
185 for totalNumElem in [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200, 102400,204800]:
186 for problem in [solveScalar,solveVector]:
187 #for problem in [solveVector]:
188 error=max([problem(numDim, totalNumElem, len_x0, alpha,"PCG",""),error])
189 error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
190 #if totalNumElem*2**numDim*numDim< 200000: error=max([problem(numDim, totalNumElem, len_x0, alpha,"DIRECT",""),error])
191 # for solver_method in [ "PCG" ]:
192 # for prec in [ "JACOBI", "ILU0" ]:
193 # error=max([problem(numDim, totalNumElem, len_x0, alpha, solver_method,prec),error])
194 print "----------------"
195 print "maximum error over all tests is ",error
196 print "----------------"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26