# Contents of /trunk/doc/examples/cookbook/example07.py

Revision 3001 - (show annotations)
Wed Mar 31 04:29:10 2010 UTC (10 years, 5 months ago) by ahallam
File MIME type: text/x-python
File size: 4593 byte(s)
```Pressue wave problem, looking at sampling theorem, and Crank-Nicolson
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2009-2010 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2009-2010 by University of Queensland 15 Earth Systems Science Computational Center (ESSCC) 16 http://www.uq.edu.au/esscc 17 Primary Business: Queensland, Australia""" 18 __license__="""Licensed under the Open Software License version 3.0 19 20 __url__= 21 22 # Antony Hallam 23 # Acoustic Wave Equation Simulation 24 25 # Importing all the necessary modules required. 26 from esys.escript import * 27 from esys.finley import Rectangle 28 import sys 29 import os 30 # smoothing operator 31 from esys.escript.pdetools import Projector 32 import numpy as np 33 import pylab as pl 34 import matplotlib.cm as cm 35 from esys.escript.linearPDEs import LinearPDE 36 37 # Establish a save path. 38 savepath = "data/example07" 39 mkDir(savepath) 40 41 #Geometric and material property related variables. 42 mx = 1000. # model lenght 43 my = 1000. # model width 44 ndx = 400 # steps in x direction 45 ndy = 400 # steps in y direction 46 47 xstep=mx/ndx 48 ystep=my/ndy 49 50 c=380.0 51 csq=c*c 52 # Time related variables. 53 tend=1.5 #end time 54 #calculating )the timestep 55 h=tend/1000. 56 #Check to make sure number of time steps is not too large. 57 print "Time step size= ",h, "Expected number of outputs= ",tend/h 58 59 #uncomment the following lines to give the user a chance to stop 60 #proceeder = raw_input("Is this ok?(y/n)") 61 #Exit if user thinks too many outputs. 62 #if proceeder == "n": 63 # sys.exit() 64 65 U0=0.01 # amplitude of point source 66 # spherical source at middle of bottom face 67 68 xc=[500,500] 69 70 mydomain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 71 #wavesolver2d(mydomain,h,tend,lam,mu,rho,U0,xc,savepath,output="mpl") 72 x=mydomain.getX() 73 74 # ... open new PDE ... 75 mypde=LinearPDE(mydomain) 76 print mypde.isUsingLumping() 77 print mypde.getSolverOptions() 78 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) 79 mypde.setSymmetryOn() 80 #kmat = kronecker(mydomain) 81 mypde.setValue(D=1.)#kmat) 82 83 # define small radius around point xc 84 # Lsup(x) returns the maximum value of the argument x 85 src_radius = 30#2*Lsup(domain.getSize()) 86 print "src_radius = ",src_radius 87 88 # ... set initial values .... 89 n=0 90 # for first two time steps 91 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius) 92 u_m1=u 93 t=0 94 95 #plot source shape 96 uT=np.array(u.toListOfTuples()) 97 uT=np.reshape(uT,(ndx+1,ndy+1)) 98 source_line=uT[ndx/2,:] 99 pl.plot(source_line) 100 pl.plot(source_line,'ro') 101 pl.axis([70,130,0,0.2]) 102 pl.savefig(os.path.join(savepath,"source_line.png")) 103 #~ u_pc_x1 = u_pot[0,0] 104 #~ u_pc_y1 = u_pot[0,1] 105 #~ u_pc_x2 = u_pot[1,0] 106 #~ u_pc_y2 = u_pot[1,1] 107 #~ u_pc_x3 = u_pot[2,0] 108 #~ u_pc_y3 = u_pot[2,1] 109 #~ 110 #~ # open file to save displacement at point source 111 #~ u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w') 112 #~ u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3)) 113 114 while t