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

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

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

trunk/esys2/finley/test/python/AdvectivePDETest.py revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC temp/finley/test/python/AdvectivePDETest.py revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3    #
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    
17  # Test for the AdvectivePDE class  # Test for the AdvectivePDE class
18  #  #
# Line 25  Line 39 
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)  #   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    __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  from esys.escript import *  from esys.escript import *
48  from esys.linearPDEs import AdvectivePDE,LinearPDE  from esys.escript.linearPDEs import AdvectivePDE,LinearPDE
49  from esys import finley  from esys import finley
50  from random import random  from random import random
51    
# Line 61  def makeRandomFloatTensor(l,k,m,n,val_lo Line 80  def makeRandomFloatTensor(l,k,m,n,val_lo
80      return out      return out
81    
82  ne=20  ne=20
83  for d in [2,3]:  # for d in [2,3]:
84    for d in [3]:
85     # create domain:     # create domain:
86     if d==2:     if d==2:
87       mydomain=finley.Rectangle(ne,ne,1)       mydomain=finley.Rectangle(ne,ne,1)
88       x=mydomain.getX()       x=mydomain.getX()
89       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.)
90     else:     else:
91       mydomain=finley.Brick(ne,ne,ne,1)       mydomain=finley.Brick(ne,ne,ne,1)
92       x=mydomain.getX()       x=mydomain.getX()
93       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.)
94     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)
95     # for ncomp in [1,2]:     # for ncomp in [1,2]:
96     for ncomp in [1,2]:     for ncomp in [2]:
97          if ncomp==1:          if ncomp==1:
98             maskf=1.             maskf=1.
99             Z=makeRandomFloats(d,-1.,0.)             Z=makeRandomFloats(d,-1.,0.)
# Line 86  for d in [2,3]: Line 106  for d in [2,3]:
106             K=numarray.zeros([ncomp,d,ncomp,d])*0.             K=numarray.zeros([ncomp,d,ncomp,d])*0.
107             for i in range(ncomp):             for i in range(ncomp):
108               K[i,:,i,:]=numarray.identity(d)*1.               K[i,:,i,:]=numarray.identity(d)*1.
109            K_sup=numarray.array(K_sup)
110          Z=numarray.array(Z)          Z=numarray.array(Z)
111          Z/=length(Z)          Z/=length(Z)
112          if ncomp==1:          if ncomp==1:
# Line 96  for d in [2,3]: Line 117  for d in [2,3]:
117             Zx=x[0]*Z[:,0]             Zx=x[0]*Z[:,0]
118             for j in range(1,d):             for j in range(1,d):
119                Zx+=x[j]*Z[:,j]                Zx+=x[j]*Z[:,j]
120          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)  
121          K/=length(K)          K/=length(K)
122          if ncomp==1:          if ncomp==1:
123             U=numarray.matrixmultiply(numarray.transpose(K),Z)             U=numarray.matrixmultiply(numarray.transpose(K),Z)
124          else:          else:
125             U=numarray.zeros([ncomp,d,ncomp])*0.             U=numarray.zeros([ncomp,d,ncomp])*1.
126             for m in range(ncomp):             for m in range(ncomp):
127                for l in range(ncomp):                for l in range(ncomp):
128                   for j in range(d):                   for j in range(d):
# Line 112  for d in [2,3]: Line 132  for d in [2,3]:
132          # create domain:          # create domain:
133          mypde=AdvectivePDE(mydomain)          mypde=AdvectivePDE(mydomain)
134          # mypde.setSolverMethod(mypde.DIRECT)          # mypde.setSolverMethod(mypde.DIRECT)
135          mypde.setValue(q=msk*maskf,A=K)          print K
136            mypde.setValue(q=msk*maskf)
137            mypde.setValue(A=K)
138            mypde.setValue(A=K,q=msk*maskf)
139            mypde.checkSymmetry()
140          # run Peclet          # run Peclet
141          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.]:
142             print "@@@  Peclet Number :",Pe*length(U)/length(K)             peclet=Pe*length(U)/2./length(K)/ne
143               print "@@@  Peclet Number :",peclet
144             u_ex=exp(Pe*Zx)             u_ex=exp(Pe*Zx)
145             mypde.setValue(r=u_ex)             mypde.setValue(r=u_ex)
146             mypde.setValue(B=Data(),C=Pe*U)             # mypde.setValue(B=Data(),C=Pe*U)
147             u=mypde.getSolution()             # u=mypde.getSolution()
148             print "@@@@ C=U: Pe = ",Pe,printError(u,u_ex)             # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex)
149             mypde.setValue(C=Data(),B=-Pe*U)             mypde.setValue(C=Data(),B=-Pe*U)
150             u=mypde.getSolution()             u=mypde.getSolution()
151             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.1387

  ViewVC Help
Powered by ViewVC 1.1.26