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

Revision 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 7656 byte(s)
```Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.

```
 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 ############################################################FILE HEADER 23 # example08b.py 24 # Antony Hallam 25 # Seismic Wave Equation Simulation using acceleration solution. 26 # Extend the solution in example 08a to use absorbing boundary 27 # conditions. 28 29 #######################################################EXTERNAL MODULES 30 from esys.escript import * 31 from esys.finley import Rectangle 32 from esys.weipa import saveVTK 33 import os 34 # smoothing operator 35 from esys.escript.pdetools import Projector, Locator 36 from esys.escript.unitsSI import * 37 import numpy as np 38 import matplotlib 39 matplotlib.use('agg') #It's just here for automated testing 40 41 import pylab as pl 42 import matplotlib.cm as cm 43 from esys.escript.linearPDEs import LinearPDE 44 45 ########################################################MPI WORLD CHECK 46 if getMPISizeWorld() > 1: 47 import sys 48 print("This example will not run in an MPI world.") 49 sys.exit(0) 50 51 #################################################ESTABLISHING VARIABLES 52 # where to save output data 53 savepath = "data/example08b" 54 mkDir(savepath) 55 #Geometric and material property related variables. 56 mx = 1000. # model lenght 57 my = 1000. # model width 58 ndx = 300 # steps in x direction 59 ndy = 300 # steps in y direction 60 xstep=mx/ndx # calculate the size of delta x 61 ystep=abs(my/ndy) # calculate the size of delta y 62 lam=3.462e9 #lames constant 63 mu=3.462e9 #bulk modulus 64 rho=1154. #density 65 # Time related variables. 66 testing=True 67 if testing: 68 print('The testing end time is currently selected. This severely limits the number of time iterations.') 69 print("Try changing testing to False for more iterations.") 70 tend=0.001 71 else: 72 tend=0.5 # end time 73 74 h=0.0001 # time step 75 # data recording times 76 rtime=0.0 # first time to record 77 rtime_inc=tend/50.0 # time increment to record 78 #Check to make sure number of time steps is not too large. 79 print("Time step size= ",h, "Expected number of outputs= ",tend/h) 80 81 U0=0.1 # amplitude of point source 82 ls=500 # length of the source 83 source=np.zeros(ls,'float') # source array 84 decay1=np.zeros(ls,'float') # decay curve one 85 decay2=np.zeros(ls,'float') # decay curve two 86 time=np.zeros(ls,'float') # time values 87 g=np.log(0.01)/ls 88 89 dfeq=50 #Dominant Frequency 90 a = 2.0 * (np.pi * dfeq)**2.0 91 t0 = 5.0 / (2.0 * np.pi * dfeq) 92 srclength = 5. * t0 93 ls = int(srclength/h) 94 print('source length',ls) 95 source=np.zeros(ls,'float') # source array 96 ampmax=0 97 for it in range(0,ls): 98 t = it*h 99 tt = t-t0 100 dum1 = np.exp(-a * tt * tt) 101 source[it] = -2. * a * tt * dum1 102 # source[it] = exp(-a * tt * tt) !gaussian 103 if (abs(source[it]) > ampmax): 104 ampmax = abs(source[it]) 105 #source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls))*(np.exp(-.1*g*t)-1) 106 #decay1[t]=np.exp(g*t) 107 #decay2[t]=(np.exp(-.1*g*t)-1) 108 time[t]=t*h 109 #tdecay=decay1*decay2*U0 110 #decay1=decay1*U0; decay2=decay2*U0 111 pl.clf(); 112 pl.plot(source) 113 #pl.plot(time,decay1);pl.plot(time,decay2); 114 #pl.plot(time,tdecay) 115 pl.savefig(os.path.join(savepath,'source.png')) 116 117 # will introduce a spherical source at middle left of bottom face 118 xc=[mx/2,0] 119 120 ####################################################DOMAIN CONSTRUCTION 121 domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain 122 x=domain.getX() # get the locations of the nodes in the domani 123 124 ##########################################################ESTABLISH PDE 125 mypde=LinearPDE(domain) # create pde 126 mypde.setSymmetryOn() # turn symmetry on 127 # turn lumping on for more efficient solving 128 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING) 129 kmat = kronecker(domain) # create the kronecker delta function of the domain 130 mypde.setValue(D=kmat*rho) #set the general form value D 131 132 ##########################################################ESTABLISH ABC 133 # Define where the boundary decay will be applied. 134 bn=50. 135 bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn) 136 # btop=ystep*bn # don't apply to force boundary!!! 137 138 # locate these points in the domain 139 left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot 140 141 tgamma=0.85 # decay value for exponential function 142 def calc_gamma(G,npts): 143 func=np.sqrt(abs(-1.*np.log(G)/(npts**2.))) 144 return func 145 146 gleft = calc_gamma(tgamma,bleft) 147 gright = calc_gamma(tgamma,bleft) 148 gbottom= calc_gamma(tgamma,ystep*bn) 149 150 print('gamma', gleft,gright,gbottom) 151 152 # calculate decay functions 153 def abc_bfunc(gamma,loc,x,G): 154 func=exp(-1.*(gamma*abs(loc-x))**2.) 155 return func 156 157 fleft=abc_bfunc(gleft,bleft,x[0],tgamma) 158 fright=abc_bfunc(gright,bright,x[0],tgamma) 159 fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma) 160 # apply these functions only where relevant 161 abcleft=fleft*whereNegative(left) 162 abcright=fright*wherePositive(right) 163 abcbottom=fbottom*wherePositive(bottom) 164 # make sure the inside of the abc is value 1 165 abcleft=abcleft+whereZero(abcleft) 166 abcright=abcright+whereZero(abcright) 167 abcbottom=abcbottom+whereZero(abcbottom) 168 # multiply the conditions together to get a smooth result 169 abc=abcleft*abcright*abcbottom 170 171 #visualise the boundary function 172 #abcT=abc.toListOfTuples() 173 #abcT=np.reshape(abcT,(ndx+1,ndy+1)) 174 #pl.clf(); pl.imshow(abcT); pl.colorbar(); 175 #pl.savefig(os.path.join(savepath,"abc.png")) 176 177 178 ############################################FIRST TIME STEPS AND SOURCE 179 # define small radius around point xc 180 src_length = 40; print("src_length = ",src_length) 181 # set initial values for first two time steps with source terms 182 y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length) 183 src_dir=numpy.array([0.,1.]) # defines direction of point source as down 184 y=y*src_dir 185 mypde.setValue(y=y) #set the source as a function on the boundary 186 # turn lumping on for more efficient solving 187 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING) 188 # for first two time steps 189 u=[0.0,0.0]*wherePositive(x) 190 u_m1=u 191 192 ####################################################ITERATION VARIABLES 193 n=0 # iteration counter 194 t=0 # time counter 195 ##############################################################ITERATION 196 while t= rtime): 206 saveVTK(os.path.join(savepath,"ex08b.%05d.vtu"%n),displacement=length(u),\ 207 acceleration=length(accel),tensor=stress) 208 rtime=rtime+rtime_inc #increment data save time 209 # increment loop values 210 t=t+h; n=n+1 211 if (n < ls): 212 y=source[n]*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length) 213 y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the boundary 214 print(n,"-th time step t ",t)