/[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 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years ago) by ksteube
File MIME type: text/x-python
File size: 5313 byte(s)
Copyright updated in all python files

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