/[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 108 - (show annotations)
Thu Jan 27 06:21:59 2005 UTC (14 years, 10 months ago) by jgs
File MIME type: text/x-python
File size: 10016 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26