test/python/AdvectivePDETest.py

# $Id$

# Test for the AdvectivePDE class
#
#  for a single equation the test problem is
#
#   -(K_{ij}u_{,j})_{,i} - (w_i u)_{,i} + v_j u_{,j} =0 
#
#   + constraints on the surface
#
#  for system of two equation the test problem is
#
#   -(K_{milj}u_{l,j})_{,i} - (w_{mil} u_l)_{,i} + v_{mlj} u_{l,j} =0 
#
#   + constraints on the surface
#
#   K,w and v are constant (we will set v=0 or w=0)
#
#   the test solution is  u(x)=e^{z_i*x_i}  and u_l(x)=e^{z_{li}*x_i}
#
#   an easy caculation shows that    
#    
#     z_i*K_{ij}*z_j=(v_i-w_i)*z_i and  z_{li}*K_{milj}*z_{lj}=(v_{mjl}-w_{mlj})*z_{lj}
#
#   obviously one can choose: v_i-w_i=K_{ji}z_j and v_{mjl}-w_{mlj}=z_{li}*K_{milj} (no summation over l)
#

__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
from esys.escript import *
from esys.escript.linearPDEs import AdvectivePDE,LinearPDE
from esys import finley
from random import random

def printError(u,u_ex):
    if u.getRank()==0:
       out=" error = %e range = [%e:%e] [%e:%e]"%(Lsup(u-u_ex)/Lsup(u_ex),sup(u),inf(u),sup(u_ex),inf(u_ex))
    else:
       out="\n"
       for i in range(u.getShape()[0]):
          out+="     %d error = %e range = [%e:%e] [%e:%e]\n"%(i,Lsup(u[i]-u_ex[i])/Lsup(u_ex[i]),sup(u[i]),inf(u[i]),sup(u_ex[i]),inf(u_ex[i]))
    return out
    
    
def makeRandomFloats(n,val_low=0.,val_up=1.):
    out=[]
    for i in range(n):
         out.append((val_up-val_low)*random()+val_low)
    return out

def makeRandomFloatMatrix(m,n,val_low=0.,val_up=1.):
    out=[]
    for i in range(m):
         out.append(makeRandomFloats(n,val_low,val_up))
    return out

def makeRandomFloatTensor(l,k,m,n,val_low=0.,val_up=1.):
    out=[]
    for j in range(l):
         out2=[]
         for i in range(k): out2.append(makeRandomFloatMatrix(m,n,val_low,val_up))
         out.append(out2)
    return out

ne=20
# for d in [2,3]:
for d in [3]:
   # create domain:
   if d==2:
     mydomain=finley.Rectangle(ne,ne,1)
     x=mydomain.getX()
     msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)
   else:
     mydomain=finley.Brick(ne,ne,ne,1)
     x=mydomain.getX()
     msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)+whereZero(x[2])+whereZero(x[2]-1.)
   print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne)
   # for ncomp in [1,2]:
   for ncomp in [2]:
        if ncomp==1:
           maskf=1.
           Z=makeRandomFloats(d,-1.,0.)
           K_sup=makeRandomFloatMatrix(d,d,-1.,1.)
           K=numarray.identity(d)*1.
        else:
           maskf=numarray.ones(ncomp)
           Z=makeRandomFloatMatrix(ncomp,d,-1.,0.)
           K_sup=makeRandomFloatTensor(ncomp,d,ncomp,d,-1.,1.)
           K=numarray.zeros([ncomp,d,ncomp,d])*0.
           for i in range(ncomp):
             K[i,:,i,:]=numarray.identity(d)*1.
        K_sup=numarray.array(K_sup)
        Z=numarray.array(Z)
        Z/=length(Z)
        if ncomp==1:
           Zx=Z[0]*x[0]
           for j in range(1,d):
              Zx+=Z[j]*x[j]
        else:
           Zx=x[0]*Z[:,0]
           for j in range(1,d):
              Zx+=x[j]*Z[:,j]
        K+=0.02*K_sup/length(K_sup)
        K/=length(K)
        if ncomp==1:
           U=numarray.matrixmultiply(numarray.transpose(K),Z)
        else:
           U=numarray.zeros([ncomp,d,ncomp])*1.
           for m in range(ncomp):
              for l in range(ncomp):
                 for j in range(d):
                    for i in range(d):
                       U[m,j,l]+=K[m,i,l,j]*Z[l,i]

        # create domain:
        mypde=AdvectivePDE(mydomain)
        # mypde.setSolverMethod(mypde.DIRECT)
        print K
        mypde.setValue(q=msk*maskf)
        mypde.setValue(A=K)
        mypde.setValue(A=K,q=msk*maskf)
        mypde.checkSymmetry()
        # run Peclet
        for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.,10000000.]:
           peclet=Pe*length(U)/2./length(K)/ne
           print "@@@  Peclet Number :",peclet
           u_ex=exp(Pe*Zx)
           mypde.setValue(r=u_ex)
           # mypde.setValue(B=Data(),C=Pe*U)
           # u=mypde.getSolution()
           # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
           mypde.setValue(C=Data(),B=-Pe*U)
           u=mypde.getSolution()
           print "@@@@ B=-U: Pe = ",peclet,printError(u,u_ex)
1	jgs	108	# $Id$
2
3			# Test for the AdvectivePDE class
4			#
5			# for a single equation the test problem is
6			#
7			# -(K_{ij}u_{,j})_{,i} - (w_i u)_{,i} + v_j u_{,j} =0
8			#
9			# + constraints on the surface
10			#
11			# for system of two equation the test problem is
12			#
13			# -(K_{milj}u_{l,j})_{,i} - (w_{mil} u_l)_{,i} + v_{mlj} u_{l,j} =0
14			#
15			# + constraints on the surface
16			#
17			# K,w and v are constant (we will set v=0 or w=0)
18			#
19			# the test solution is u(x)=e^{z_ix_i} and u_l(x)=e^{z_{li}x_i}
20			#
21			# an easy caculation shows that
22			#
23			# z_iK_{ij}z_j=(v_i-w_i)z_i and z_{li}K_{milj}z_{lj}=(v_{mjl}-w_{mlj})z_{lj}
24			#
25			# obviously one can choose: v_i-w_i=K_{ji}z_j and v_{mjl}-w_{mlj}=z_{li}*K_{milj} (no summation over l)
26			#
27
28	elspeth	617	__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
29			http://www.access.edu.au
30			Primary Business: Queensland, Australia"""
31			__license__="""Licensed under the Open Software License version 3.0
32			http://www.opensource.org/licenses/osl-3.0.php"""
33	jgs	108	from esys.escript import *
34	gross	423	from esys.escript.linearPDEs import AdvectivePDE,LinearPDE
35	jgs	108	from esys import finley
36			from random import random
37
38			def printError(u,u_ex):
39			if u.getRank()==0:
40			out=" error = %e range = [%e:%e] [%e:%e]"%(Lsup(u-u_ex)/Lsup(u_ex),sup(u),inf(u),sup(u_ex),inf(u_ex))
41			else:
42			out="\n"
43			for i in range(u.getShape()[0]):
44			out+=" %d error = %e range = [%e:%e] [%e:%e]\n"%(i,Lsup(u[i]-u_ex[i])/Lsup(u_ex[i]),sup(u[i]),inf(u[i]),sup(u_ex[i]),inf(u_ex[i]))
45			return out
46
47
48			def makeRandomFloats(n,val_low=0.,val_up=1.):
49			out=[]
50			for i in range(n):
51			out.append((val_up-val_low)*random()+val_low)
52			return out
53
54			def makeRandomFloatMatrix(m,n,val_low=0.,val_up=1.):
55			out=[]
56			for i in range(m):
57			out.append(makeRandomFloats(n,val_low,val_up))
58			return out
59
60			def makeRandomFloatTensor(l,k,m,n,val_low=0.,val_up=1.):
61			out=[]
62			for j in range(l):
63			out2=[]
64			for i in range(k): out2.append(makeRandomFloatMatrix(m,n,val_low,val_up))
65			out.append(out2)
66			return out
67
68			ne=20
69	jgs	110	# for d in [2,3]:
70			for d in [3]:
71	jgs	108	# create domain:
72			if d==2:
73			mydomain=finley.Rectangle(ne,ne,1)
74			x=mydomain.getX()
75	gross	423	msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)
76	jgs	108	else:
77			mydomain=finley.Brick(ne,ne,ne,1)
78			x=mydomain.getX()
79	gross	423	msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)+whereZero(x[2])+whereZero(x[2]-1.)
80	jgs	108	print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne)
81			# for ncomp in [1,2]:
82	jgs	110	for ncomp in [2]:
83	jgs	108	if ncomp==1:
84			maskf=1.
85			Z=makeRandomFloats(d,-1.,0.)
86			K_sup=makeRandomFloatMatrix(d,d,-1.,1.)
87			K=numarray.identity(d)*1.
88			else:
89			maskf=numarray.ones(ncomp)
90			Z=makeRandomFloatMatrix(ncomp,d,-1.,0.)
91			K_sup=makeRandomFloatTensor(ncomp,d,ncomp,d,-1.,1.)
92			K=numarray.zeros([ncomp,d,ncomp,d])*0.
93			for i in range(ncomp):
94			K[i,:,i,:]=numarray.identity(d)*1.
95	jgs	110	K_sup=numarray.array(K_sup)
96	jgs	108	Z=numarray.array(Z)
97			Z/=length(Z)
98			if ncomp==1:
99			Zx=Z[0]*x[0]
100			for j in range(1,d):
101			Zx+=Z[j]*x[j]
102			else:
103			Zx=x[0]*Z[:,0]
104			for j in range(1,d):
105			Zx+=x[j]*Z[:,j]
106	jgs	110	K+=0.02*K_sup/length(K_sup)
107	jgs	108	K/=length(K)
108			if ncomp==1:
109			U=numarray.matrixmultiply(numarray.transpose(K),Z)
110			else:
111	jgs	110	U=numarray.zeros([ncomp,d,ncomp])*1.
112	jgs	108	for m in range(ncomp):
113			for l in range(ncomp):
114			for j in range(d):
115			for i in range(d):
116			U[m,j,l]+=K[m,i,l,j]*Z[l,i]
117
118			# create domain:
119			mypde=AdvectivePDE(mydomain)
120			# mypde.setSolverMethod(mypde.DIRECT)
121	gross	423	print K
122			mypde.setValue(q=msk*maskf)
123			mypde.setValue(A=K)
124			mypde.setValue(A=K,q=msk*maskf)
125	jgs	110	mypde.checkSymmetry()
126	jgs	108	# run Peclet
127	jgs	110	for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.,10000000.]:
128			peclet=Pe*length(U)/2./length(K)/ne
129			print "@@@ Peclet Number :",peclet
130	jgs	108	u_ex=exp(Pe*Zx)
131			mypde.setValue(r=u_ex)
132	jgs	110	# mypde.setValue(B=Data(),C=Pe*U)
133			# u=mypde.getSolution()
134			# print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
135	jgs	108	mypde.setValue(C=Data(),B=-Pe*U)
136			u=mypde.getSolution()
137	jgs	110	print "@@@@ B=-U: Pe = ",peclet,printError(u,u_ex)
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision