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 423 by gross, Fri Jan 6 08:09:43 2006 UTC
# Line 26  Line 26
26  #  #
27
28  from esys.escript import *  from esys.escript import *
30  from esys import finley  from esys import finley
31  from random import random  from random import random
32
# 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)
69       x=mydomain.getX()       x=mydomain.getX()
70       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.)
71     else:     else:
72       mydomain=finley.Brick(ne,ne,ne,1)       mydomain=finley.Brick(ne,ne,ne,1)
73       x=mydomain.getX()       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()       msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)+whereZero(x[2])+whereZero(x[2]-1.)
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:
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 112  for d in [2,3]: Line 113  for d in [2,3]:
113          # create domain:          # create domain:
115          # mypde.setSolverMethod(mypde.DIRECT)          # mypde.setSolverMethod(mypde.DIRECT)
118            mypde.setValue(A=K)
120            mypde.checkSymmetry()
121          # run Peclet          # run Peclet
122          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.]:
123             print "@@@  Peclet Number :",Pe*length(U)/length(K)             peclet=Pe*length(U)/2./length(K)/ne
124               print "@@@  Peclet Number :",peclet
125             u_ex=exp(Pe*Zx)             u_ex=exp(Pe*Zx)
126             mypde.setValue(r=u_ex)             mypde.setValue(r=u_ex)
127             mypde.setValue(B=Data(),C=Pe*U)             # mypde.setValue(B=Data(),C=Pe*U)
128             u=mypde.getSolution()             # u=mypde.getSolution()
129             print "@@@@ C=U: Pe = ",Pe,printError(u,u_ex)             # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
130             mypde.setValue(C=Data(),B=-Pe*U)             mypde.setValue(C=Data(),B=-Pe*U)
131             u=mypde.getSolution()             u=mypde.getSolution()
132             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.423