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

Annotation of /trunk/finley/test/python/AssembleTest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (hide annotations)
Wed Dec 15 07:08:39 2004 UTC (16 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/AssembleTest.py
File MIME type: text/x-python
File size: 10384 byte(s)
*** empty log message ***

1 jgs 82 # 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 jgs 102 from linearPDEs import *
17 jgs 82
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 jgs 102 mypde1=LinearPDE(domain=mydomain,A=c_A2,X=eval(x,elem))
206 jgs 82 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 jgs 102 mypde2=LinearPDE(domain=mydomain,Y=eval(y,elem),y=matmult(eval(x,face_elem),nrml))
211 jgs 82 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 jgs 102 mypde1=LinearPDE(domain=mydomain,B=c_B2,X=eval(x,elem))
225 jgs 82 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 jgs 102 mypde2=LinearPDE(domain=mydomain,Y=eval(y,elem),y=matmult(eval(x,face_elem),nrml))
230 jgs 82 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 jgs 102 mypde1=LinearPDE(domain=mydomain,C=c_C2,Y=eval(y,elem))
244 jgs 82 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 jgs 102 mypde1=LinearPDE(domain=mydomain,D=c_D2,Y=eval(y,elem))
259 jgs 82 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