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

Revision 4005 - (show annotations)
Fri Sep 28 06:09:03 2012 UTC (8 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 7722 byte(s)
```test fixes, doco updates, annoyance removals.

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