trunk/esys2/finley/test/python/AdvectivePDETest.py revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC trunk/finley/test/python/AdvectivePDETest.py revision 617 by elspeth, Wed Mar 22 02:58:17 2006 UTC
# Line 25  Line 25
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)  #   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
29                        http://www.access.edu.au
33  from esys.escript import *  from esys.escript import *
35  from esys import finley  from esys import finley
36  from random import random  from random import random
37
# Line 61  def makeRandomFloatTensor(l,k,m,n,val_lo Line 66  def makeRandomFloatTensor(l,k,m,n,val_lo
66      return out      return out
67
68  ne=20  ne=20
69  for d in [2,3]:  # for d in [2,3]:
70    for d in [3]:
71     # create domain:     # create domain:
72     if d==2:     if d==2:
73       mydomain=finley.Rectangle(ne,ne,1)       mydomain=finley.Rectangle(ne,ne,1)
74       x=mydomain.getX()       x=mydomain.getX()
75       msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero()       msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)
76     else:     else:
77       mydomain=finley.Brick(ne,ne,ne,1)       mydomain=finley.Brick(ne,ne,ne,1)
78       x=mydomain.getX()       x=mydomain.getX()
79       msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero()+x[2].whereZero()+(x[2]-1.).whereZero()       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)     print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne)
81     # for ncomp in [1,2]:     # for ncomp in [1,2]:
82     for ncomp in [1,2]:     for ncomp in [2]:
83          if ncomp==1:          if ncomp==1:
85             Z=makeRandomFloats(d,-1.,0.)             Z=makeRandomFloats(d,-1.,0.)
# Line 86  for d in [2,3]: Line 92  for d in [2,3]:
92             K=numarray.zeros([ncomp,d,ncomp,d])*0.             K=numarray.zeros([ncomp,d,ncomp,d])*0.
93             for i in range(ncomp):             for i in range(ncomp):
94               K[i,:,i,:]=numarray.identity(d)*1.               K[i,:,i,:]=numarray.identity(d)*1.
95            K_sup=numarray.array(K_sup)
96          Z=numarray.array(Z)          Z=numarray.array(Z)
97          Z/=length(Z)          Z/=length(Z)
98          if ncomp==1:          if ncomp==1:
# Line 96  for d in [2,3]: Line 103  for d in [2,3]:
103             Zx=x[0]*Z[:,0]             Zx=x[0]*Z[:,0]
104             for j in range(1,d):             for j in range(1,d):
105                Zx+=x[j]*Z[:,j]                Zx+=x[j]*Z[:,j]
106          K_sup=numarray.array(makeRandomFloatMatrix(d,d,-1.,1.))          K+=0.02*K_sup/length(K_sup)
K+=0.05*K_sup/length(K_sup)
107          K/=length(K)          K/=length(K)
108          if ncomp==1:          if ncomp==1:
109             U=numarray.matrixmultiply(numarray.transpose(K),Z)             U=numarray.matrixmultiply(numarray.transpose(K),Z)
110          else:          else:
111             U=numarray.zeros([ncomp,d,ncomp])*0.             U=numarray.zeros([ncomp,d,ncomp])*1.
112             for m in range(ncomp):             for m in range(ncomp):
113                for l in range(ncomp):                for l in range(ncomp):
114                   for j in range(d):                   for j in range(d):
# Line 112  for d in [2,3]: Line 118  for d in [2,3]:
118          # create domain:          # create domain:
120          # mypde.setSolverMethod(mypde.DIRECT)          # mypde.setSolverMethod(mypde.DIRECT)
123            mypde.setValue(A=K)
125            mypde.checkSymmetry()
126          # run Peclet          # run Peclet
127          for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.]:          for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.,10000000.]:
128             print "@@@  Peclet Number :",Pe*length(U)/length(K)             peclet=Pe*length(U)/2./length(K)/ne
129               print "@@@  Peclet Number :",peclet
130             u_ex=exp(Pe*Zx)             u_ex=exp(Pe*Zx)
131             mypde.setValue(r=u_ex)             mypde.setValue(r=u_ex)
132             mypde.setValue(B=Data(),C=Pe*U)             # mypde.setValue(B=Data(),C=Pe*U)
133             u=mypde.getSolution()             # u=mypde.getSolution()
134             print "@@@@ C=U: Pe = ",Pe,printError(u,u_ex)             # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
135             mypde.setValue(C=Data(),B=-Pe*U)             mypde.setValue(C=Data(),B=-Pe*U)
136             u=mypde.getSolution()             u=mypde.getSolution()
137             print "@@@@ B=-U: Pe = ",Pe,printError(u,u_ex)             print "@@@@ B=-U: Pe = ",peclet,printError(u,u_ex)

Legend:
 Removed from v.108 changed lines Added in v.617