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	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			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