test/python/AssembleTest.py

# assemble test: 
#
# $Id$

from esys.escript import *
from esys.linearPDEs import *
import esys.finley as 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(mydomain)
               mypde1.setValue(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(mydomain)
               mypde2.setValue(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(mydomain)
            mypde1.setValue(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(mydomain)
            mypde2.setValue(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(mydomain)
            mypde1.setValue(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(mydomain)
         mypde1.setValue(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	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_12,x_22,x_32,....,x_1order_u,x_2order_u,x_3order_u] for dim=3
23	# or
24	# B=[1,x_1,x_2,x_12,x_22,....,x_1order_u,x_2order_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,pdim+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,pdim+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,pdim+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,pdim+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.setReducedOrderForSolutionsTo(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,jdim+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
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision