/[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 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 9 months ago) by phornby
Original Path: temp_trunk_copy/finley/test/python/AdvectivePDETest.py
File MIME type: text/x-python
File size: 4976 byte(s)
Make a temp copy of the trunk before checking in the windows changes


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