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

Diff of /trunk/esys2/finley/test/python/AdvectivePDETest.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 109 by jgs, Thu Jan 27 06:21:59 2005 UTC revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC
# Line 61  def makeRandomFloatTensor(l,k,m,n,val_lo Line 61  def makeRandomFloatTensor(l,k,m,n,val_lo
61      return out      return out
62    
63  ne=20  ne=20
64  for d in [2,3]:  # for d in [2,3]:
65    for d in [3]:
66     # create domain:     # create domain:
67     if d==2:     if d==2:
68       mydomain=finley.Rectangle(ne,ne,1)       mydomain=finley.Rectangle(ne,ne,1)
# Line 73  for d in [2,3]: Line 74  for d in [2,3]:
74       msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero()+x[2].whereZero()+(x[2]-1.).whereZero()       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)     print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne)
76     # for ncomp in [1,2]:     # for ncomp in [1,2]:
77     for ncomp in [1,2]:     for ncomp in [2]:
78          if ncomp==1:          if ncomp==1:
79             maskf=1.             maskf=1.
80             Z=makeRandomFloats(d,-1.,0.)             Z=makeRandomFloats(d,-1.,0.)
# Line 86  for d in [2,3]: Line 87  for d in [2,3]:
87             K=numarray.zeros([ncomp,d,ncomp,d])*0.             K=numarray.zeros([ncomp,d,ncomp,d])*0.
88             for i in range(ncomp):             for i in range(ncomp):
89               K[i,:,i,:]=numarray.identity(d)*1.               K[i,:,i,:]=numarray.identity(d)*1.
90            K_sup=numarray.array(K_sup)
91          Z=numarray.array(Z)          Z=numarray.array(Z)
92          Z/=length(Z)          Z/=length(Z)
93          if ncomp==1:          if ncomp==1:
# Line 96  for d in [2,3]: Line 98  for d in [2,3]:
98             Zx=x[0]*Z[:,0]             Zx=x[0]*Z[:,0]
99             for j in range(1,d):             for j in range(1,d):
100                Zx+=x[j]*Z[:,j]                Zx+=x[j]*Z[:,j]
101          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)  
102          K/=length(K)          K/=length(K)
103          if ncomp==1:          if ncomp==1:
104             U=numarray.matrixmultiply(numarray.transpose(K),Z)             U=numarray.matrixmultiply(numarray.transpose(K),Z)
105          else:          else:
106             U=numarray.zeros([ncomp,d,ncomp])*0.             U=numarray.zeros([ncomp,d,ncomp])*1.
107             for m in range(ncomp):             for m in range(ncomp):
108                for l in range(ncomp):                for l in range(ncomp):
109                   for j in range(d):                   for j in range(d):
# Line 113  for d in [2,3]: Line 114  for d in [2,3]:
114          mypde=AdvectivePDE(mydomain)          mypde=AdvectivePDE(mydomain)
115          # mypde.setSolverMethod(mypde.DIRECT)          # mypde.setSolverMethod(mypde.DIRECT)
116          mypde.setValue(q=msk*maskf,A=K)          mypde.setValue(q=msk*maskf,A=K)
117            mypde.checkSymmetry()
118          # run Peclet          # run Peclet
119          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.]:
120             print "@@@  Peclet Number :",Pe*length(U)/length(K)             peclet=Pe*length(U)/2./length(K)/ne
121               print "@@@  Peclet Number :",peclet
122             u_ex=exp(Pe*Zx)             u_ex=exp(Pe*Zx)
123             mypde.setValue(r=u_ex)             mypde.setValue(r=u_ex)
124             mypde.setValue(B=Data(),C=Pe*U)             # mypde.setValue(B=Data(),C=Pe*U)
125             u=mypde.getSolution()             # u=mypde.getSolution()
126             print "@@@@ C=U: Pe = ",Pe,printError(u,u_ex)             # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
127             mypde.setValue(C=Data(),B=-Pe*U)             mypde.setValue(C=Data(),B=-Pe*U)
128             u=mypde.getSolution()             u=mypde.getSolution()
129             print "@@@@ B=-U: Pe = ",Pe,printError(u,u_ex)             print "@@@@ B=-U: Pe = ",peclet,printError(u,u_ex)

Legend:
Removed from v.109  
changed lines
  Added in v.110

  ViewVC Help
Powered by ViewVC 1.1.26