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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (hide annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years ago) by jgs
File MIME type: text/x-python
File size: 4200 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

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_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     from esys.escript import *
29     from esys.linearPDEs import AdvectivePDE,LinearPDE
30     from esys import finley
31     from random import random
32    
33     def printError(u,u_ex):
34     if u.getRank()==0:
35     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))
36     else:
37     out="\n"
38     for i in range(u.getShape()[0]):
39     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]))
40     return out
41    
42    
43     def makeRandomFloats(n,val_low=0.,val_up=1.):
44     out=[]
45     for i in range(n):
46     out.append((val_up-val_low)*random()+val_low)
47     return out
48    
49     def makeRandomFloatMatrix(m,n,val_low=0.,val_up=1.):
50     out=[]
51     for i in range(m):
52     out.append(makeRandomFloats(n,val_low,val_up))
53     return out
54    
55     def makeRandomFloatTensor(l,k,m,n,val_low=0.,val_up=1.):
56     out=[]
57     for j in range(l):
58     out2=[]
59     for i in range(k): out2.append(makeRandomFloatMatrix(m,n,val_low,val_up))
60     out.append(out2)
61     return out
62    
63     ne=20
64 jgs 110 # for d in [2,3]:
65     for d in [3]:
66 jgs 108 # create domain:
67     if d==2:
68     mydomain=finley.Rectangle(ne,ne,1)
69     x=mydomain.getX()
70     msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero()
71     else:
72     mydomain=finley.Brick(ne,ne,ne,1)
73     x=mydomain.getX()
74     msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero()+x[2].whereZero()+(x[2]-1.).whereZero()
75     print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne)
76     # for ncomp in [1,2]:
77 jgs 110 for ncomp in [2]:
78 jgs 108 if ncomp==1:
79     maskf=1.
80     Z=makeRandomFloats(d,-1.,0.)
81     K_sup=makeRandomFloatMatrix(d,d,-1.,1.)
82     K=numarray.identity(d)*1.
83     else:
84     maskf=numarray.ones(ncomp)
85     Z=makeRandomFloatMatrix(ncomp,d,-1.,0.)
86     K_sup=makeRandomFloatTensor(ncomp,d,ncomp,d,-1.,1.)
87     K=numarray.zeros([ncomp,d,ncomp,d])*0.
88     for i in range(ncomp):
89     K[i,:,i,:]=numarray.identity(d)*1.
90 jgs 110 K_sup=numarray.array(K_sup)
91 jgs 108 Z=numarray.array(Z)
92     Z/=length(Z)
93     if ncomp==1:
94     Zx=Z[0]*x[0]
95     for j in range(1,d):
96     Zx+=Z[j]*x[j]
97     else:
98     Zx=x[0]*Z[:,0]
99     for j in range(1,d):
100     Zx+=x[j]*Z[:,j]
101 jgs 110 K+=0.02*K_sup/length(K_sup)
102 jgs 108 K/=length(K)
103     if ncomp==1:
104     U=numarray.matrixmultiply(numarray.transpose(K),Z)
105     else:
106 jgs 110 U=numarray.zeros([ncomp,d,ncomp])*1.
107 jgs 108 for m in range(ncomp):
108     for l in range(ncomp):
109     for j in range(d):
110     for i in range(d):
111     U[m,j,l]+=K[m,i,l,j]*Z[l,i]
112    
113     # create domain:
114     mypde=AdvectivePDE(mydomain)
115     # mypde.setSolverMethod(mypde.DIRECT)
116     mypde.setValue(q=msk*maskf,A=K)
117 jgs 110 mypde.checkSymmetry()
118 jgs 108 # run Peclet
119 jgs 110 for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.,10000000.]:
120     peclet=Pe*length(U)/2./length(K)/ne
121     print "@@@ Peclet Number :",peclet
122 jgs 108 u_ex=exp(Pe*Zx)
123     mypde.setValue(r=u_ex)
124 jgs 110 # mypde.setValue(B=Data(),C=Pe*U)
125     # u=mypde.getSolution()
126     # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
127 jgs 108 mypde.setValue(C=Data(),B=-Pe*U)
128     u=mypde.getSolution()
129 jgs 110 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