test/python/AssembleTest.py

# assemble test: 
#
# $Id$

import sys
import os
import unittest
                                                                                                                                                          
esys_root=os.getenv('ESYS_ROOT')
sys.path.append(esys_root+'/finley/lib')
sys.path.append(esys_root+'/escript/lib')
sys.path.append(esys_root+'/escript/py_src')
                                                                                                                                                          
from escript import *
from util import *
from linearPDE import *

import finley
from math import *
global seed
num_elem=2 # number of elements in each spatial direction
num_equations=3 # number of equations
integration_order=-1 # order of the integration scheme
order_u=1 # solution order


seed=1
rand_factor=sqrt(2.)
#
#   The test solution is represented in the basis 
#
#   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
#                 or 
#   B=[1,x_1,x_2,x_1**2,x_2**2,....,x_1**order_u,x_2**order_u]   for dim=2
#
#   length of B is len_Basis.
#
#   any test solution u with numComp components is represented by the numComp x len_Basis matrix sol.
#   the actual solution is then given by u=matmul(sol,B). 
#

global maxError,coeff,total_maxError,total_coeff
total_maxError=0
total_coeff=""
maxError=0
coeff=""
D=()


def randomNum():
     global seed
     seed+=1
     s=rand_factor*seed
     return s-int(s)
    
def generateRandom(arg):
    if len(arg)==0:
       return randomNum()
    out=numarray.ones(arg,numarray.Float64) # *randomNum()
    return out

def algebraicGrad(u):
   out=numarray.zeros([u.shape[0],dim,len_Basis],numarray.Float64)
   for s in range(u.shape[0]):
     for i in range(dim):
        for k in range(len_Basis):
          h=0
          for j in range(len_Basis):
             h+=u[s,j]*D[i][k,j]
          out[s,i,k]=h
   return out

def algebraicDiv(u):
   out=numarray.zeros([u.shape[0],len_Basis],numarray.Float64)
   for s in range(u.shape[0]):
     for k in range(len_Basis):
        h=0
        for i in range(dim):
          for j in range(len_Basis):
             h+=u[s,i,j]*D[i][k,j]
        out[s,k]=h
   return out

def mult4(A,u):
      out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
        for i in range(A.shape[1]):
           for k in range(len_Basis):
             h=0
             for t in range(A.shape[2]):
               for j in range(A.shape[3]):
                 h+=A[s,i,t,j]*u[t,j,k]
             out[s,i,k]=h
      return out

def mult3_2(A,u):
      out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
        for i in range(A.shape[1]):
           for k in range(len_Basis):
             h=0
             for t in range(A.shape[2]):
                 h+=A[s,i,t]*u[t,k]
             out[s,i,k]=h
      return out

def mult3_1(A,u):
      out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
           for k in range(len_Basis):
             h=0
             for t in range(A.shape[1]):
               for j in range(A.shape[2]):
                 h+=A[s,t,j]*u[t,j,k]
             out[s,k]=h
      return out

def mult2(A,u):
      out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
         for k in range(len_Basis):
           h=0
           for t in range(A.shape[1]):
               h+=A[s,t]*u[t,k]
           out[s,k]=h
      return out


def eval(u,this):
    x=this.getX()
    if u.rank==2:
       out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
       for i0 in range(u.shape[0]):
          out[i0]=u[i0,0]
          for p in range(order_u):
             for j in range(dim):
               out[i0]+=u[i0,p*dim+j+1]*x[j]**(p+1)
    else:
       if u.shape[0]==1:
          out=Data(value=0,shape=(u.shape[1],),what=this,expand=True)
          for i1 in range(u.shape[1]):
             out[i1]=u[0,i1,0]
             for p in range(order_u):
                for j in range(dim):
                  out[i1]+=u[0,i1,p*dim+j+1]*x[j]**(p+1)

       elif u.shape[1]==1:
          out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
          for i0 in range(u.shape[0]):
             out[i0]=u[i0,0,0]
             for p in range(order_u):
                for j in range(dim):
                  out[i0]+=u[i0,0,p*dim+j+1]*x[j]**(p+1)
       else:
          out=Data(value=0,shape=(u.shape[0],u.shape[1]),what=this,expand=True)
          for i0 in range(u.shape[0]):
           for i1 in range(u.shape[1]):
             out[i0,i1]=u[i0,i1,0]
             for p in range(order_u):
                for j in range(dim):
                  out[i0,i1]+=u[i0,i1,p*dim+j+1]*x[j]**(p+1)
    return out

def checkSystem(text,operator,u,rhs):
    global maxError,coeff
    error=Lsup(operator*u-rhs)
    print "@@ "+text+" error: ",error
    if error>=maxError:
       maxError=error
       coeff=text
#
#
def TestSystem(numEqu,numComp,mydomain,reduce):
   elem=Function(mydomain)
   face_elem=FunctionOnBoundary(mydomain)
   nodes=ContinuousFunction(mydomain)
   nrml=face_elem.getNormal()
   #
   #   test solution:
   #
   u=generateRandom([numComp,len_Basis])
   # u=numarray.zeros([numComp,len_Basis])
   # u[0,0]=1
   # u[0,1]=1
   U=eval(u,nodes)
   gradu=algebraicGrad(u)
   #
   #   test A:
   #
   for p in range(numEqu):
      for q in range(numComp):
         for i in range(dim):
            for j in range(dim):

               # check div( A grad(u) ) = div( X )
               c_A=numarray.zeros([numEqu,dim,numComp,dim])
               c_A[p,i,q,j]=1
               if numEqu==1 and numComp==1:
                  c_A2=numarray.reshape(c_A,[dim,dim])
                  text="A[%d,%d]"%(i,j)
               else:
                  c_A2=c_A
                  text="A[%d,%d,%d,%d]"%(p,i,q,j)
               x=mult4(c_A,gradu)
               mypde1=linearPDE(domain=mydomain,A=c_A2,X=eval(x,elem))
               mypde1.setReducedOrderForSolutionsTo(reduce)
               checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
               # check div( A grad(u) ) = Y
               y=-algebraicDiv(x)
               mypde2=linearPDE(domain=mydomain,Y=eval(y,elem),y=matmult(eval(x,face_elem),nrml))
               mypde2.setReducedOrderForSolutionsTo(reduce)
               checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())

            # check div( B u ) = div( X )
            c_B=numarray.zeros([numEqu,dim,numComp])
            c_B[p,i,q]=1
            if numEqu==1 and numComp==1:
               c_B2=numarray.reshape(c_B,[dim])
               text="B[%d]"%(i)
            else:
               c_B2=c_B
               text="B[%d,%d,%d]"%(p,i,q)
            x=mult3_2(c_B,u)
            mypde1=linearPDE(domain=mydomain,B=c_B2,X=eval(x,elem))
            mypde1.setReducedOrderForSolutionsTo(reduce)
            checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
            # check div( B u ) = Y
            y=-algebraicDiv(x)
            mypde2=linearPDE(domain=mydomain,Y=eval(y,elem),y=matmult(eval(x,face_elem),nrml))
            mypde2.setReducedOrderForSolutionsTo(reduce)
            checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())

            # check C grad(u) = Y
            c_C=numarray.zeros([numEqu,numComp,dim])
            c_C[p,q,i]=1
            if numEqu==1 and numComp==1:
               c_C2=numarray.reshape(c_C,[dim])
               text="C[%d]"%(i)
            else:
               c_C2=c_C
               text="C[%d,%d,%d]"%(p,q,i)
            y=mult3_1(c_C,gradu)
            mypde1=linearPDE(domain=mydomain,C=c_C2,Y=eval(y,elem))
            mypde1.setReducedOrderForSolutionsTo(reduce)
            checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())


         # check D u= Y
         c_D=numarray.zeros([numEqu,numComp])
         c_D[p,q]=1
         if numEqu==1 and numComp==1:
            c_D2=numarray.reshape(c_D,[1])
            text="D"
         else:
            c_D2=c_D
            text="D[%d,%d]"%(p,q)
         y=mult2(c_D,u)
         mypde1=linearPDE(domain=mydomain,D=c_D2,Y=eval(y,elem))
         mypde1.setReducedOrderForSolutionsTo(reduce)
         checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())

#  we start:
for dim in [2,3]:
    len_Basis=dim*order_u+1
    #
    #  the differential operators:
    #   
    #  D_j(matmul(sol,B))=matmul(sol,D_jB)=matmul(matmul(sol,D[j]),B)
    #
    D=()
    for i in range(dim):
      D=D+(numarray.zeros([len_Basis,len_Basis],numarray.Float64),)
      for j in range(order_u):
         if j==0:
           D[i][0,i+j+1]=1
         else:
           D[i][(j-1)*dim+i+1,j*dim+i+1]=j+1
    #
    #   generate mydomain:
    #
    for order in [1,2]:
     for onElements in [False,True]:
        if onElements==True:
            onElmtext=", with elements on faces"
        else:
            onElmtext=""
        if dim==3:
           mydomain=finley.Brick(num_elem,num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
        else:
           mydomain=finley.Rectangle(num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
        for reduce in [False,True]:
           if reduce==True:
              redtext=",reduced"
           else:
              redtext=""
           #
           #  and start the test process:
           #
           for numEqu in range(1,num_equations+1):
               print "@@@ Start testing assembling with dim=%d and %d equations (order=%d%s%s)"%(dim,numEqu,order,redtext,onElmtext)
               TestSystem(numEqu,numEqu,mydomain,reduce)
               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)
               if maxError>=total_maxError:
                  total_maxError=maxError
                  total_coeff=coeff+" with %d equations in %dD (order=%d%s%s)"%(numEqu,dim,order,redtext,onElmtext)
    
print "@@@@ end testing assemblage with maximal error %e for %s"%(total_maxError,total_coeff)


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 linearPDE 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_12,x_22,x_32,....,x_1order_u,x_2order_u,x_3order_u] for dim=3
33	# or
34	# B=[1,x_1,x_2,x_12,x_22,....,x_1order_u,x_2order_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,pdim+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,pdim+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,pdim+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,pdim+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,jdim+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
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision