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

Revision 4005 - (show annotations)
Fri Sep 28 06:09:03 2012 UTC (6 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 6471 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 # example09.py 25 # Antony Hallam 26 # Seismic Wave Equation Simulation using acceleration solution. 27 # 3D model with multiple layers. 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 from esys.finley import ReadMesh 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/example09c" 55 meshpath = "data/example09n" 56 mkDir(savepath) 57 #Geometric and material property related variables. 58 domain=ReadMesh(os.path.join(savepath,'example09n.fly')) # create the domain 59 x=Solution(domain).getX() 60 #parameters layers 1,2,3,4 and fault 61 prho=np.array([2200.,2500.,3200.,4500.,5500.]) #density 62 pvel=np.array([1500.,2200.,3000.,3200.,5000.]) #velocity 63 pmu=pvel**2.*prho/4. #bulk modulus 64 plam=pvel**2.*prho/2. #lames constant 65 nlayers=4 66 width=300.0 67 rho=Scalar(0,Function(domain)) 68 vel=Scalar(0,Function(domain)) 69 mu=Scalar(0,Function(domain)) 70 lam=Scalar(0,Function(domain)) 71 72 print(0.5*np.sqrt(prho/(plam+2*pmu))*0.5) 73 74 for i in range(0,nlayers): 75 rho.setTaggedValue('lblock%d'%i,prho[i]) 76 rho.setTaggedValue('rblock%d'%i,prho[i]) 77 vel.setTaggedValue('lblock%d'%i,pvel[i]) 78 vel.setTaggedValue('rblock%d'%i,pvel[i]) 79 mu.setTaggedValue('lblock%d'%i,pmu[i]) 80 mu.setTaggedValue('rblock%d'%i,pmu[i]) 81 lam.setTaggedValue('lblock%d'%i,plam[i]) 82 lam.setTaggedValue('rblock%d'%i,plam[i]) 83 i=nlayers 84 rho.setTaggedValue('fault',prho[i]) 85 vel.setTaggedValue('fault',pvel[i]) 86 mu.setTaggedValue('fault',pmu[i]) 87 lam.setTaggedValue('fault',plam[i]) 88 89 90 # Time related variables. 91 testing=False 92 if testing: 93 print('The testing end time is currently selected. This severely limits the number of time iterations.') 94 print("Try changing testing to False for more iterations.") 95 tend=0.1 96 else: 97 tend=0.1 # end time 98 99 h=0.00001 # time step 100 # data recording times 101 rtime=0.0 # first time to record 102 rtime_inc=tend/750.0 # time increment to record 103 #Check to make sure number of time steps is not too large. 104 print("Time step size= ",h, "Expected number of outputs= ",tend/h) 105 106 U0=0.1 # amplitude of point source 107 ls=500 # length of the source 108 source=np.zeros(ls,'float') # source array 109 decay1=np.zeros(ls,'float') # decay curve one 110 decay2=np.zeros(ls,'float') # decay curve two 111 time=np.zeros(ls,'float') # time values 112 g=np.log(0.01)/ls 113 114 dfeq=50 #Dominant Frequency 115 a = 2.0 * (np.pi * dfeq)**2.0 116 t0 = 5.0 / (2.0 * np.pi * dfeq) 117 srclength = 5. * t0 118 ls = int(srclength/h) 119 print('source length',ls) 120 source=np.zeros(ls,'float') # source array 121 ampmax=0 122 for it in range(0,ls): 123 t = it*h 124 tt = t-t0 125 dum1 = np.exp(-a * tt * tt) 126 source[it] = -2. * a * tt * dum1 127 if (abs(source[it]) > ampmax): 128 ampmax = abs(source[it]) 129 time[t]=t*h 130 131 # will introduce a spherical source at middle left of bottom face 132 xc=[150,0] 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_length = 10; print("src_length = ",src_length) 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_length)+1)*whereNegative(length(xb-xc)-src_length) 148 stop=Scalar(0.0,FunctionOnBoundary(domain)) 149 stop.setTaggedValue("top",1.0) 150 src_dir=numpy.array([0.,-1.]) # 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 154 # initial value of displacement at point source is constant (U0=0.01) 155 # for first two time steps 156 u=[0.0,0.0]*x 157 u_m1=u 158 159 ####################################################ITERATION VARIABLES 160 n=0 # iteration counter 161 t=0 # time counter 162 ##############################################################ITERATION 163 while t= rtime): 173 saveVTK(os.path.join(savepath,"ex09c.%05d.vtu"%n),displacement=length(u),\ 174 acceleration=length(accel),tensor=stress) 175 rtime=rtime+rtime_inc #increment data save time 176 # increment loop values 177 t=t+h; n=n+1 178 if (n < ls): 179 mypde.setValue(y=source[n]*yx*src_dir*stop) #set the source as a function on the boundary 180 print("time step %d, t=%s"%(n,t))

 ViewVC Help Powered by ViewVC 1.1.26