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

Revision 108 - (show annotations)
Thu Jan 27 06:21:59 2005 UTC (16 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/AdvectivePDETest.py
File MIME type: text/x-python
File size: 4129 byte(s)
```*** empty log message ***

```
 1 # \$Id\$ 2 3 # Test for the AdvectivePDE class 4 # 5 # for a single equation the test problem is 6 # 7 # -(K_{ij}u_{,j})_{,i} - (w_i u)_{,i} + v_j u_{,j} =0 8 # 9 # + constraints on the surface 10 # 11 # for system of two equation the test problem is 12 # 13 # -(K_{milj}u_{l,j})_{,i} - (w_{mil} u_l)_{,i} + v_{mlj} u_{l,j} =0 14 # 15 # + constraints on the surface 16 # 17 # K,w and v are constant (we will set v=0 or w=0) 18 # 19 # the test solution is u(x)=e^{z_i*x_i} and u_l(x)=e^{z_{li}*x_i} 20 # 21 # an easy caculation shows that 22 # 23 # 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} 24 # 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) 26 # 27 28 from esys.escript import * 29 from esys.linearPDEs import AdvectivePDE,LinearPDE 30 from esys import finley 31 from random import random 32 33 def printError(u,u_ex): 34 if u.getRank()==0: 35 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)) 36 else: 37 out="\n" 38 for i in range(u.getShape()[0]): 39 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])) 40 return out 41 42 43 def makeRandomFloats(n,val_low=0.,val_up=1.): 44 out=[] 45 for i in range(n): 46 out.append((val_up-val_low)*random()+val_low) 47 return out 48 49 def makeRandomFloatMatrix(m,n,val_low=0.,val_up=1.): 50 out=[] 51 for i in range(m): 52 out.append(makeRandomFloats(n,val_low,val_up)) 53 return out 54 55 def makeRandomFloatTensor(l,k,m,n,val_low=0.,val_up=1.): 56 out=[] 57 for j in range(l): 58 out2=[] 59 for i in range(k): out2.append(makeRandomFloatMatrix(m,n,val_low,val_up)) 60 out.append(out2) 61 return out 62 63 ne=20 64 for d in [2,3]: 65 # create domain: 66 if d==2: 67 mydomain=finley.Rectangle(ne,ne,1) 68 x=mydomain.getX() 69 msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero() 70 else: 71 mydomain=finley.Brick(ne,ne,ne,1) 72 x=mydomain.getX() 73 msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero()+x[2].whereZero()+(x[2]-1.).whereZero() 74 print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne) 75 # for ncomp in [1,2]: 76 for ncomp in [1,2]: 77 if ncomp==1: 78 maskf=1. 79 Z=makeRandomFloats(d,-1.,0.) 80 K_sup=makeRandomFloatMatrix(d,d,-1.,1.) 81 K=numarray.identity(d)*1. 82 else: 83 maskf=numarray.ones(ncomp) 84 Z=makeRandomFloatMatrix(ncomp,d,-1.,0.) 85 K_sup=makeRandomFloatTensor(ncomp,d,ncomp,d,-1.,1.) 86 K=numarray.zeros([ncomp,d,ncomp,d])*0. 87 for i in range(ncomp): 88 K[i,:,i,:]=numarray.identity(d)*1. 89 Z=numarray.array(Z) 90 Z/=length(Z) 91 if ncomp==1: 92 Zx=Z[0]*x[0] 93 for j in range(1,d): 94 Zx+=Z[j]*x[j] 95 else: 96 Zx=x[0]*Z[:,0] 97 for j in range(1,d): 98 Zx+=x[j]*Z[:,j] 99 K_sup=numarray.array(makeRandomFloatMatrix(d,d,-1.,1.)) 100 K+=0.05*K_sup/length(K_sup) 101 K/=length(K) 102 if ncomp==1: 103 U=numarray.matrixmultiply(numarray.transpose(K),Z) 104 else: 105 U=numarray.zeros([ncomp,d,ncomp])*0. 106 for m in range(ncomp): 107 for l in range(ncomp): 108 for j in range(d): 109 for i in range(d): 110 U[m,j,l]+=K[m,i,l,j]*Z[l,i] 111 112 # create domain: 113 mypde=AdvectivePDE(mydomain) 114 # mypde.setSolverMethod(mypde.DIRECT) 115 mypde.setValue(q=msk*maskf,A=K) 116 # run Peclet 117 for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.]: 118 print "@@@ Peclet Number :",Pe*length(U)/length(K) 119 u_ex=exp(Pe*Zx) 120 mypde.setValue(r=u_ex) 121 mypde.setValue(B=Data(),C=Pe*U) 122 u=mypde.getSolution() 123 print "@@@@ C=U: Pe = ",Pe,printError(u,u_ex) 124 mypde.setValue(C=Data(),B=-Pe*U) 125 u=mypde.getSolution() 126 print "@@@@ B=-U: Pe = ",Pe,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