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

Contents of /trunk/esys2/finley/test/python/AssembleTest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (15 years, 6 months ago) by jgs
File MIME type: text/x-python
File size: 10384 byte(s)
*** empty log message ***

1 # assemble test:
2 #
3 # $Id$
4
5 import sys
6 import os
7 import unittest
8
9 esys_root=os.getenv('ESYS_ROOT')
10 sys.path.append(esys_root+'/finley/lib')
11 sys.path.append(esys_root+'/escript/lib')
12 sys.path.append(esys_root+'/escript/py_src')
13
14 from escript import *
15 from util import *
16 from linearPDEs import *
17
18 import finley
19 from math import *
20 global seed
21 num_elem=2 # number of elements in each spatial direction
22 num_equations=3 # number of equations
23 integration_order=-1 # order of the integration scheme
24 order_u=1 # solution order
25
26
27 seed=1
28 rand_factor=sqrt(2.)
29 #
30 # The test solution is represented in the basis
31 #
32 # B=[1,x_1,x_2,x_3,x_1**2,x_2**2,x_3**2,....,x_1**order_u,x_2**order_u,x_3**order_u] for dim=3
33 # or
34 # B=[1,x_1,x_2,x_1**2,x_2**2,....,x_1**order_u,x_2**order_u] for dim=2
35 #
36 # length of B is len_Basis.
37 #
38 # any test solution u with numComp components is represented by the numComp x len_Basis matrix sol.
39 # the actual solution is then given by u=matmul(sol,B).
40 #
41
42 global maxError,coeff,total_maxError,total_coeff
43 total_maxError=0
44 total_coeff=""
45 maxError=0
46 coeff=""
47 D=()
48
49
50 def randomNum():
51 global seed
52 seed+=1
53 s=rand_factor*seed
54 return s-int(s)
55
56 def generateRandom(arg):
57 if len(arg)==0:
58 return randomNum()
59 out=numarray.ones(arg,numarray.Float64) # *randomNum()
60 return out
61
62 def algebraicGrad(u):
63 out=numarray.zeros([u.shape[0],dim,len_Basis],numarray.Float64)
64 for s in range(u.shape[0]):
65 for i in range(dim):
66 for k in range(len_Basis):
67 h=0
68 for j in range(len_Basis):
69 h+=u[s,j]*D[i][k,j]
70 out[s,i,k]=h
71 return out
72
73 def algebraicDiv(u):
74 out=numarray.zeros([u.shape[0],len_Basis],numarray.Float64)
75 for s in range(u.shape[0]):
76 for k in range(len_Basis):
77 h=0
78 for i in range(dim):
79 for j in range(len_Basis):
80 h+=u[s,i,j]*D[i][k,j]
81 out[s,k]=h
82 return out
83
84 def mult4(A,u):
85 out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
86 for s in range(A.shape[0]):
87 for i in range(A.shape[1]):
88 for k in range(len_Basis):
89 h=0
90 for t in range(A.shape[2]):
91 for j in range(A.shape[3]):
92 h+=A[s,i,t,j]*u[t,j,k]
93 out[s,i,k]=h
94 return out
95
96 def mult3_2(A,u):
97 out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
98 for s in range(A.shape[0]):
99 for i in range(A.shape[1]):
100 for k in range(len_Basis):
101 h=0
102 for t in range(A.shape[2]):
103 h+=A[s,i,t]*u[t,k]
104 out[s,i,k]=h
105 return out
106
107 def mult3_1(A,u):
108 out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
109 for s in range(A.shape[0]):
110 for k in range(len_Basis):
111 h=0
112 for t in range(A.shape[1]):
113 for j in range(A.shape[2]):
114 h+=A[s,t,j]*u[t,j,k]
115 out[s,k]=h
116 return out
117
118 def mult2(A,u):
119 out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
120 for s in range(A.shape[0]):
121 for k in range(len_Basis):
122 h=0
123 for t in range(A.shape[1]):
124 h+=A[s,t]*u[t,k]
125 out[s,k]=h
126 return out
127
128
129 def eval(u,this):
130 x=this.getX()
131 if u.rank==2:
132 out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
133 for i0 in range(u.shape[0]):
134 out[i0]=u[i0,0]
135 for p in range(order_u):
136 for j in range(dim):
137 out[i0]+=u[i0,p*dim+j+1]*x[j]**(p+1)
138 else:
139 if u.shape[0]==1:
140 out=Data(value=0,shape=(u.shape[1],),what=this,expand=True)
141 for i1 in range(u.shape[1]):
142 out[i1]=u[0,i1,0]
143 for p in range(order_u):
144 for j in range(dim):
145 out[i1]+=u[0,i1,p*dim+j+1]*x[j]**(p+1)
146
147 elif u.shape[1]==1:
148 out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
149 for i0 in range(u.shape[0]):
150 out[i0]=u[i0,0,0]
151 for p in range(order_u):
152 for j in range(dim):
153 out[i0]+=u[i0,0,p*dim+j+1]*x[j]**(p+1)
154 else:
155 out=Data(value=0,shape=(u.shape[0],u.shape[1]),what=this,expand=True)
156 for i0 in range(u.shape[0]):
157 for i1 in range(u.shape[1]):
158 out[i0,i1]=u[i0,i1,0]
159 for p in range(order_u):
160 for j in range(dim):
161 out[i0,i1]+=u[i0,i1,p*dim+j+1]*x[j]**(p+1)
162 return out
163
164 def checkSystem(text,operator,u,rhs):
165 global maxError,coeff
166 error=Lsup(operator*u-rhs)
167 print "@@ "+text+" error: ",error
168 if error>=maxError:
169 maxError=error
170 coeff=text
171 #
172 #
173 def TestSystem(numEqu,numComp,mydomain,reduce):
174 elem=Function(mydomain)
175 face_elem=FunctionOnBoundary(mydomain)
176 nodes=ContinuousFunction(mydomain)
177 nrml=face_elem.getNormal()
178 #
179 # test solution:
180 #
181 u=generateRandom([numComp,len_Basis])
182 # u=numarray.zeros([numComp,len_Basis])
183 # u[0,0]=1
184 # u[0,1]=1
185 U=eval(u,nodes)
186 gradu=algebraicGrad(u)
187 #
188 # test A:
189 #
190 for p in range(numEqu):
191 for q in range(numComp):
192 for i in range(dim):
193 for j in range(dim):
194
195 # check div( A grad(u) ) = div( X )
196 c_A=numarray.zeros([numEqu,dim,numComp,dim])
197 c_A[p,i,q,j]=1
198 if numEqu==1 and numComp==1:
199 c_A2=numarray.reshape(c_A,[dim,dim])
200 text="A[%d,%d]"%(i,j)
201 else:
202 c_A2=c_A
203 text="A[%d,%d,%d,%d]"%(p,i,q,j)
204 x=mult4(c_A,gradu)
205 mypde1=LinearPDE(domain=mydomain,A=c_A2,X=eval(x,elem))
206 mypde1.setReducedOrderForSolutionsTo(reduce)
207 checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
208 # check div( A grad(u) ) = Y
209 y=-algebraicDiv(x)
210 mypde2=LinearPDE(domain=mydomain,Y=eval(y,elem),y=matmult(eval(x,face_elem),nrml))
211 mypde2.setReducedOrderForSolutionsTo(reduce)
212 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())
213
214 # check div( B u ) = div( X )
215 c_B=numarray.zeros([numEqu,dim,numComp])
216 c_B[p,i,q]=1
217 if numEqu==1 and numComp==1:
218 c_B2=numarray.reshape(c_B,[dim])
219 text="B[%d]"%(i)
220 else:
221 c_B2=c_B
222 text="B[%d,%d,%d]"%(p,i,q)
223 x=mult3_2(c_B,u)
224 mypde1=LinearPDE(domain=mydomain,B=c_B2,X=eval(x,elem))
225 mypde1.setReducedOrderForSolutionsTo(reduce)
226 checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
227 # check div( B u ) = Y
228 y=-algebraicDiv(x)
229 mypde2=LinearPDE(domain=mydomain,Y=eval(y,elem),y=matmult(eval(x,face_elem),nrml))
230 mypde2.setReducedOrderForSolutionsTo(reduce)
231 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())
232
233 # check C grad(u) = Y
234 c_C=numarray.zeros([numEqu,numComp,dim])
235 c_C[p,q,i]=1
236 if numEqu==1 and numComp==1:
237 c_C2=numarray.reshape(c_C,[dim])
238 text="C[%d]"%(i)
239 else:
240 c_C2=c_C
241 text="C[%d,%d,%d]"%(p,q,i)
242 y=mult3_1(c_C,gradu)
243 mypde1=LinearPDE(domain=mydomain,C=c_C2,Y=eval(y,elem))
244 mypde1.setReducedOrderForSolutionsTo(reduce)
245 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())
246
247
248 # check D u= Y
249 c_D=numarray.zeros([numEqu,numComp])
250 c_D[p,q]=1
251 if numEqu==1 and numComp==1:
252 c_D2=numarray.reshape(c_D,[1])
253 text="D"
254 else:
255 c_D2=c_D
256 text="D[%d,%d]"%(p,q)
257 y=mult2(c_D,u)
258 mypde1=LinearPDE(domain=mydomain,D=c_D2,Y=eval(y,elem))
259 mypde1.setReducedOrderForSolutionsTo(reduce)
260 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())
261
262 # we start:
263 for dim in [2,3]:
264 len_Basis=dim*order_u+1
265 #
266 # the differential operators:
267 #
268 # D_j(matmul(sol,B))=matmul(sol,D_jB)=matmul(matmul(sol,D[j]),B)
269 #
270 D=()
271 for i in range(dim):
272 D=D+(numarray.zeros([len_Basis,len_Basis],numarray.Float64),)
273 for j in range(order_u):
274 if j==0:
275 D[i][0,i+j+1]=1
276 else:
277 D[i][(j-1)*dim+i+1,j*dim+i+1]=j+1
278 #
279 # generate mydomain:
280 #
281 for order in [1,2]:
282 for onElements in [False,True]:
283 if onElements==True:
284 onElmtext=", with elements on faces"
285 else:
286 onElmtext=""
287 if dim==3:
288 mydomain=finley.Brick(num_elem,num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
289 else:
290 mydomain=finley.Rectangle(num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
291 for reduce in [False,True]:
292 if reduce==True:
293 redtext=",reduced"
294 else:
295 redtext=""
296 #
297 # and start the test process:
298 #
299 for numEqu in range(1,num_equations+1):
300 print "@@@ Start testing assembling with dim=%d and %d equations (order=%d%s%s)"%(dim,numEqu,order,redtext,onElmtext)
301 TestSystem(numEqu,numEqu,mydomain,reduce)
302 print "@@@ end testing assembling (order=%d%s%s) with %d equations in %dD with maximum error %e for %s"%(order,redtext,onElmtext,numEqu,dim,maxError,coeff)
303 if maxError>=total_maxError:
304 total_maxError=maxError
305 total_coeff=coeff+" with %d equations in %dD (order=%d%s%s)"%(numEqu,dim,order,redtext,onElmtext)
306
307 print "@@@@ end testing assemblage with maximal error %e for %s"%(total_maxError,total_coeff)
308
309

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26