/[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 617 - (hide annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years, 11 months ago) by elspeth
File MIME type: text/x-python
File size: 10353 byte(s)
More copyright.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26