Revision 617 - (hide annotations)
Wed Mar 22 02:58:17 2006 UTC (13 years, 7 months ago) by elspeth
File MIME type: text/x-python
File size: 4556 byte(s)
```More copyright.

```
 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 elspeth 617 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF 29 30 Primary Business: Queensland, Australia""" 31 __license__="""Licensed under the Open Software License version 3.0 32 33 jgs 108 from esys.escript import * 34 gross 423 from esys.escript.linearPDEs import AdvectivePDE,LinearPDE 35 jgs 108 from esys import finley 36 from random import random 37 38 def printError(u,u_ex): 39 if u.getRank()==0: 40 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)) 41 else: 42 out="\n" 43 for i in range(u.getShape()[0]): 44 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])) 45 return out 46 47 48 def makeRandomFloats(n,val_low=0.,val_up=1.): 49 out=[] 50 for i in range(n): 51 out.append((val_up-val_low)*random()+val_low) 52 return out 53 54 def makeRandomFloatMatrix(m,n,val_low=0.,val_up=1.): 55 out=[] 56 for i in range(m): 57 out.append(makeRandomFloats(n,val_low,val_up)) 58 return out 59 60 def makeRandomFloatTensor(l,k,m,n,val_low=0.,val_up=1.): 61 out=[] 62 for j in range(l): 63 out2=[] 64 for i in range(k): out2.append(makeRandomFloatMatrix(m,n,val_low,val_up)) 65 out.append(out2) 66 return out 67 68 ne=20 69 jgs 110 # for d in [2,3]: 70 for d in [3]: 71 jgs 108 # create domain: 72 if d==2: 73 mydomain=finley.Rectangle(ne,ne,1) 74 x=mydomain.getX() 75 gross 423 msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.) 76 jgs 108 else: 77 mydomain=finley.Brick(ne,ne,ne,1) 78 x=mydomain.getX() 79 gross 423 msk=whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1])+whereZero(x[1]-1.)+whereZero(x[2])+whereZero(x[2]-1.) 80 jgs 108 print "@ generated %d-dimension mesh with %d elements in each direction"%(d,ne) 81 # for ncomp in [1,2]: 82 jgs 110 for ncomp in [2]: 83 jgs 108 if ncomp==1: 84 maskf=1. 85 Z=makeRandomFloats(d,-1.,0.) 86 K_sup=makeRandomFloatMatrix(d,d,-1.,1.) 87 K=numarray.identity(d)*1. 88 else: 89 maskf=numarray.ones(ncomp) 90 Z=makeRandomFloatMatrix(ncomp,d,-1.,0.) 91 K_sup=makeRandomFloatTensor(ncomp,d,ncomp,d,-1.,1.) 92 K=numarray.zeros([ncomp,d,ncomp,d])*0. 93 for i in range(ncomp): 94 K[i,:,i,:]=numarray.identity(d)*1. 95 jgs 110 K_sup=numarray.array(K_sup) 96 jgs 108 Z=numarray.array(Z) 97 Z/=length(Z) 98 if ncomp==1: 99 Zx=Z[0]*x[0] 100 for j in range(1,d): 101 Zx+=Z[j]*x[j] 102 else: 103 Zx=x[0]*Z[:,0] 104 for j in range(1,d): 105 Zx+=x[j]*Z[:,j] 106 jgs 110 K+=0.02*K_sup/length(K_sup) 107 jgs 108 K/=length(K) 108 if ncomp==1: 109 U=numarray.matrixmultiply(numarray.transpose(K),Z) 110 else: 111 jgs 110 U=numarray.zeros([ncomp,d,ncomp])*1. 112 jgs 108 for m in range(ncomp): 113 for l in range(ncomp): 114 for j in range(d): 115 for i in range(d): 116 U[m,j,l]+=K[m,i,l,j]*Z[l,i] 117 118 # create domain: 119 mypde=AdvectivePDE(mydomain) 120 # mypde.setSolverMethod(mypde.DIRECT) 121 gross 423 print K 122 mypde.setValue(q=msk*maskf) 123 mypde.setValue(A=K) 124 mypde.setValue(A=K,q=msk*maskf) 125 jgs 110 mypde.checkSymmetry() 126 jgs 108 # run Peclet 127 jgs 110 for Pe in [0.001,1.,1.,10.,100,1000.,10000.,100000.,1000000.,10000000.]: 128 peclet=Pe*length(U)/2./length(K)/ne 129 print "@@@ Peclet Number :",peclet 130 jgs 108 u_ex=exp(Pe*Zx) 131 mypde.setValue(r=u_ex) 132 jgs 110 # mypde.setValue(B=Data(),C=Pe*U) 133 # u=mypde.getSolution() 134 # print "@@@@ C=U: Pe = ",peclet,printError(u,u_ex) 135 jgs 108 mypde.setValue(C=Data(),B=-Pe*U) 136 u=mypde.getSolution() 137 jgs 110 print "@@@@ B=-U: Pe = ",peclet,printError(u,u_ex)

## Properties

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