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

Revision 3911 - (show annotations)
Thu Jun 14 01:01:03 2012 UTC (7 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 6451 byte(s)
```Copyright changes
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2009-2012 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-2012 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 # example09.py 24 # Antony Hallam 25 # Seismic Wave Equation Simulation using acceleration solution. 26 # 3D model with multiple layers. 27 28 #######################################################EXTERNAL MODULES 29 from esys.escript import * 30 from esys.finley import Rectangle 31 import os 32 # smoothing operator 33 from esys.escript.pdetools import Projector, Locator 34 from esys.escript.unitsSI import * 35 import numpy as np 36 import matplotlib 37 matplotlib.use('agg') #It's just here for automated testing 38 39 import pylab as pl 40 import matplotlib.cm as cm 41 from esys.escript.linearPDEs import LinearPDE 42 from esys.finley import ReadMesh 43 from esys.weipa import saveVTK 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/example09a" 54 meshpath = "data/example09m" 55 mkDir(savepath) 56 #Geometric and material property related variables. 57 step=10.0 # the element size 58 59 vel2=1800.; vel1=3000. 60 rho2=2300.; rho1=3100. #density 61 mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus 62 lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant 63 64 ####################################################TESTING SWITCH 65 testing=True 66 if testing: 67 print('The testing end time is currently selected. This severely limits the number of time iterations.') 68 print("Try changing testing to False for more iterations.") 69 tend=0.001 70 #Model Parameters 71 mx=40. 72 my=40. 73 mz=20. 74 outputs=5 75 else: 76 tend=0.1 # end time 77 #Model Parameters 78 mx=100.0 #x width of model 79 my=100.0 #y width of model 80 mz=50.0 #depth of model 81 outputs=200 82 83 ####################################################TIME RELATED VARIABLES 84 h=0.00005 # time step 85 # data recording times 86 rtime=0.0 # first time to record 87 rtime_inc=tend/outputs # time increment to record 88 #Check to make sure number of time steps is not too large. 89 print("Time step size= ",h, "Expected number of outputs= ",tend/h) 90 91 ####################################################CREATING THE SOURCE FUNCTION 92 U0=0.1 # amplitude of point source 93 ls=500 # length of the source 94 source=np.zeros(ls,'float') # source array 95 decay1=np.zeros(ls,'float') # decay curve one 96 decay2=np.zeros(ls,'float') # decay curve two 97 time=np.zeros(ls,'float') # time values 98 g=np.log(0.01)/ls 99 100 dfeq=50 #Dominant Frequency 101 a = 2.0 * (np.pi * dfeq)**2.0 102 t0 = 5.0 / (2.0 * np.pi * dfeq) 103 srclength = 5. * t0 104 ls = int(srclength/h) 105 print('source length',ls) 106 source=np.zeros(ls,'float') # source array 107 ampmax=0 108 for it in range(0,ls): 109 t = it*h 110 tt = t-t0 111 dum1 = np.exp(-a * tt * tt) 112 source[it] = -2. * a * tt * dum1 113 if (abs(source[it]) > ampmax): 114 ampmax = abs(source[it]) 115 time[t]=t*h 116 117 # will introduce a spherical source at middle left of bottom face 118 xc=[mx/2,my/2,0] 119 120 ####################################################DOMAIN CONSTRUCTION 121 domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain 122 x=domain.getX() # get the locations of the nodes in the domain 123 124 lam=Scalar(0,Function(domain)) 125 lam.setTaggedValue("vintfa",lam1) 126 lam.setTaggedValue("vintfb",lam2) 127 mu=Scalar(0,Function(domain)) 128 mu.setTaggedValue("vintfa",mu1) 129 mu.setTaggedValue("vintfb",mu2) 130 rho=Scalar(0,Function(domain)) 131 rho.setTaggedValue("vintfa",rho1) 132 rho.setTaggedValue("vintfb",rho2) 133 134 ##########################################################ESTABLISH PDE 135 mypde=LinearPDE(domain) # create pde 136 mypde.setSymmetryOn() # turn symmetry on 137 # turn lumping on for more efficient solving 138 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING) 139 kmat = kronecker(domain) # create the kronecker delta function of the domain 140 mypde.setValue(D=rho*kmat) #set the general form value D 141 142 ############################################FIRST TIME STEPS AND SOURCE 143 # define small radius around point xc 144 src_rad = 20; print("sourc radius= ",src_rad) 145 # set initial values for first two time steps with source terms 146 xb=FunctionOnBoundary(domain).getX() 147 yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad) 148 stop=Scalar(0.0,FunctionOnBoundary(domain)) 149 stop.setTaggedValue("stop",1.0) 150 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down 151 152 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary 153 # initial value of displacement at point source is constant (U0=0.01) 154 # for first two time steps 155 u=[0.0,0.0,0.0]*wherePositive(x) 156 u_m1=u 157 158 ####################################################ITERATION VARIABLES 159 n=0 # iteration counter 160 t=0 # time counter 161 ##############################################################ITERATION 162 while t= rtime): 172 saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\ 173 acceleration=length(accel),tensor=stress) 174 rtime=rtime+rtime_inc #increment data save time 175 # increment loop values 176 t=t+h; n=n+1 177 if (n < ls): 178 mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary 179 print(n,"-th time step t ",t)