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

Contents of /trunk/finley/test/python/AdvectivePDETest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 617 - (show annotations)
Wed Mar 22 02:58:17 2006 UTC (13 years, 8 months ago) by elspeth
File MIME type: text/x-python
File size: 4556 byte(s)
More copyright.

1 # $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_i*x_i} and u_l(x)=e^{z_{li}*x_i}
20 #
21 # an easy caculation shows that
22 #
23 # 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}
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 __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 from esys.escript import *
34 from esys.escript.linearPDEs import AdvectivePDE,LinearPDE
35 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 # for d in [2,3]:
70 for d in [3]:
71 # create domain:
72 if d==2:
73 mydomain=finley.Rectangle(ne,ne,1)
74 x=mydomain.getX()
75 msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)
76 else:
77 mydomain=finley.Brick(ne,ne,ne,1)
78 x=mydomain.getX()
79 msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)+whereZero(x[2])+whereZero(x[2]-1.)
80 print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne)
81 # for ncomp in [1,2]:
82 for ncomp in [2]:
83 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 K_sup=numarray.array(K_sup)
96 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 K+=0.02*K_sup/length(K_sup)
107 K/=length(K)
108 if ncomp==1:
109 U=numarray.matrixmultiply(numarray.transpose(K),Z)
110 else:
111 U=numarray.zeros([ncomp,d,ncomp])*1.
112 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 print K
122 mypde.setValue(q=msk*maskf)
123 mypde.setValue(A=K)
124 mypde.setValue(A=K,q=msk*maskf)
125 mypde.checkSymmetry()
126 # run Peclet
127 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 u_ex=exp(Pe*Zx)
131 mypde.setValue(r=u_ex)
132 # mypde.setValue(B=Data(),C=Pe*U)
133 # u=mypde.getSolution()
134 # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
135 mypde.setValue(C=Data(),B=-Pe*U)
136 u=mypde.getSolution()
137 print "@@@@ B=-U: Pe = ",peclet,printError(u,u_ex)

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26