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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 5 months ago) by jgs
File MIME type: text/x-python
File size: 10065 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

1 jgs 82 # assemble test:
2     #
3     # $Id$
4    
5 jgs 149 from esys.escript import *
6     from esys.escript.linearPDEs import *
7     from esys import finley
8 jgs 104 from math import *
9 jgs 82
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 jgs 147 def evaluate(u,this):
120 jgs 82 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 jgs 147 U=evaluate(u,nodes)
176 jgs 82 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 jgs 104 mypde1=LinearPDE(mydomain)
196 jgs 147 mypde1.setValue(A=c_A2,X=evaluate(x,elem))
197 jgs 108 mypde1.setReducedOrderForSolutionTo(reduce)
198 jgs 82 checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
199     # check div( A grad(u) ) = Y
200     y=-algebraicDiv(x)
201 jgs 104 mypde2=LinearPDE(mydomain)
202 jgs 147 mypde2.setValue(Y=evaluate(y,elem),y=matrixmult(evaluate(x,face_elem),nrml))
203 jgs 82 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 jgs 104 mypde1=LinearPDE(mydomain)
217 jgs 147 mypde1.setValue(B=c_B2,X=evaluate(x,elem))
218 jgs 82 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 jgs 104 mypde2=LinearPDE(mydomain)
223 jgs 147 mypde2.setValue(Y=evaluate(y,elem),y=matrixmult(evaluate(x,face_elem),nrml))
224 jgs 82 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 jgs 104 mypde1=LinearPDE(mydomain)
238 jgs 147 mypde1.setValue(C=c_C2,Y=evaluate(y,elem))
239 jgs 82 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 jgs 104 mypde1=LinearPDE(mydomain)
254 jgs 147 mypde1.setValue(D=c_D2,Y=evaluate(y,elem))
255 jgs 82 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