Revision 110 - (hide annotations)
Mon Feb 14 04:14:42 2005 UTC (16 years, 2 months ago) by jgs
File MIME type: text/x-python
File size: 4200 byte(s)
```*** empty log message ***

```
 1 jgs 108 # \$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 jgs 110 # for d in [2,3]: 65 for d in [3]: 66 jgs 108 # create domain: 67 if d==2: 68 mydomain=finley.Rectangle(ne,ne,1) 69 x=mydomain.getX() 70 msk=x[0].whereZero()+(x[0]-1.).whereZero()+x[1].whereZero()+(x[1]-1.).whereZero() 71 else: 72 mydomain=finley.Brick(ne,ne,ne,1) 73 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() 75 print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne) 76 # for ncomp in [1,2]: 77 jgs 110 for ncomp in [2]: 78 jgs 108 if ncomp==1: 79 maskf=1. 80 Z=makeRandomFloats(d,-1.,0.) 81 K_sup=makeRandomFloatMatrix(d,d,-1.,1.) 82 K=numarray.identity(d)*1. 83 else: 84 maskf=numarray.ones(ncomp) 85 Z=makeRandomFloatMatrix(ncomp,d,-1.,0.) 86 K_sup=makeRandomFloatTensor(ncomp,d,ncomp,d,-1.,1.) 87 K=numarray.zeros([ncomp,d,ncomp,d])*0. 88 for i in range(ncomp): 89 K[i,:,i,:]=numarray.identity(d)*1. 90 jgs 110 K_sup=numarray.array(K_sup) 91 jgs 108 Z=numarray.array(Z) 92 Z/=length(Z) 93 if ncomp==1: 94 Zx=Z[0]*x[0] 95 for j in range(1,d): 96 Zx+=Z[j]*x[j] 97 else: 98 Zx=x[0]*Z[:,0] 99 for j in range(1,d): 100 Zx+=x[j]*Z[:,j] 101 jgs 110 K+=0.02*K_sup/length(K_sup) 102 jgs 108 K/=length(K) 103 if ncomp==1: 104 U=numarray.matrixmultiply(numarray.transpose(K),Z) 105 else: 106 jgs 110 U=numarray.zeros([ncomp,d,ncomp])*1. 107 jgs 108 for m in range(ncomp): 108 for l in range(ncomp): 109 for j in range(d): 110 for i in range(d): 111 U[m,j,l]+=K[m,i,l,j]*Z[l,i] 112 113 # create domain: 114 mypde=AdvectivePDE(mydomain) 115 # mypde.setSolverMethod(mypde.DIRECT) 116 mypde.setValue(q=msk*maskf,A=K) 117 jgs 110 mypde.checkSymmetry() 118 jgs 108 # run Peclet 119 jgs 110 for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.,10000000.]: 120 peclet=Pe*length(U)/2./length(K)/ne 121 print "@@@ Peclet Number :",peclet 122 jgs 108 u_ex=exp(Pe*Zx) 123 mypde.setValue(r=u_ex) 124 jgs 110 # mypde.setValue(B=Data(),C=Pe*U) 125 # u=mypde.getSolution() 126 # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex) 127 jgs 108 mypde.setValue(C=Data(),B=-Pe*U) 128 u=mypde.getSolution() 129 jgs 110 print "@@@@ B=-U: Pe = ",peclet,printError(u,u_ex)

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision