Annotation of /trunk/doc/examples/cookbook/example09a.py

Revision 3148 - (hide annotations)
Fri Sep 3 02:09:47 2010 UTC (11 years ago) by jfenwick
File MIME type: text/x-python
File size: 6228 byte(s)
```Another attempt to patch the X issue

```
 1 ahallam 3055 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 # 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 jfenwick 3148 import matplotlib 37 matplotlib.use('agg') #It's just here for automated testing 38 39 ahallam 3055 import pylab as pl 40 import matplotlib.cm as cm 41 from esys.escript.linearPDEs import LinearPDE 42 ahallam 3069 from esys.finley import ReadMesh 43 ahallam 3055 44 ########################################################MPI WORLD CHECK 45 if getMPISizeWorld() > 1: 46 import sys 47 print "This example will not run in an MPI world." 48 sys.exit(0) 49 50 #################################################ESTABLISHING VARIABLES 51 # where to save output data 52 ahallam 3069 savepath = "data/example09" 53 ahallam 3055 mkDir(savepath) 54 #Geometric and material property related variables. 55 ahallam 3089 mx = 200. # model lenght 56 my = 200. # model width 57 mz=100.0 58 step=10.0 # the element size 59 ahallam 3069 ndx = int(mx/step) # steps in x direction 60 ndy = int(my/step) # steps in y direction 61 ndz = int(mz/step) 62 63 vel2=1800.; vel1=3000. 64 rho2=2300.; rho1=3100. #density 65 ahallam 3089 mu2=vel2**2.*rho2/4.; mu1=vel1**2.*rho1/4. #bulk modulus 66 lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2. #lames constant 67 ahallam 3069 68 ahallam 3055 # Time related variables. 69 ahallam 3089 tend=0.1 # end time 70 h=0.0001 # time step 71 ahallam 3055 # data recording times 72 rtime=0.0 # first time to record 73 rtime_inc=tend/50.0 # time increment to record 74 #Check to make sure number of time steps is not too large. 75 print "Time step size= ",h, "Expected number of outputs= ",tend/h 76 77 U0=0.1 # amplitude of point source 78 ls=500 # length of the source 79 source=np.zeros(ls,'float') # source array 80 decay1=np.zeros(ls,'float') # decay curve one 81 decay2=np.zeros(ls,'float') # decay curve two 82 time=np.zeros(ls,'float') # time values 83 g=np.log(0.01)/ls 84 85 dfeq=50 #Dominant Frequency 86 a = 2.0 * (np.pi * dfeq)**2.0 87 t0 = 5.0 / (2.0 * np.pi * dfeq) 88 srclength = 5. * t0 89 ls = int(srclength/h) 90 print 'source length',ls 91 source=np.zeros(ls,'float') # source array 92 ampmax=0 93 for it in range(0,ls): 94 t = it*h 95 tt = t-t0 96 dum1 = np.exp(-a * tt * tt) 97 source[it] = -2. * a * tt * dum1 98 if (abs(source[it]) > ampmax): 99 ampmax = abs(source[it]) 100 time[t]=t*h 101 #tdecay=decay1*decay2*U0 102 #decay1=decay1*U0; decay2=decay2*U0 103 ahallam 3089 #pl.clf(); 104 #pl.plot(source) 105 ahallam 3055 #pl.plot(time,decay1);pl.plot(time,decay2); 106 #pl.plot(time,tdecay) 107 ahallam 3089 #pl.savefig(os.path.join(savepath,'source.png')) 108 ahallam 3055 109 # will introduce a spherical source at middle left of bottom face 110 ahallam 3089 xc=[50,50,0] 111 ahallam 3055 112 ####################################################DOMAIN CONSTRUCTION 113 ahallam 3069 domain=ReadMesh(os.path.join(savepath,'example09m.fly')) # create the domain 114 x=domain.getX() # get the locations of the nodes in the domain 115 ahallam 3055 116 ahallam 3069 lam=Scalar(0,Function(domain)) 117 lam.setTaggedValue("vintfa",lam1) 118 lam.setTaggedValue("vintfb",lam2) 119 mu=Scalar(0,Function(domain)) 120 mu.setTaggedValue("vintfa",mu1) 121 mu.setTaggedValue("vintfb",mu2) 122 rho=Scalar(0,Function(domain)) 123 rho.setTaggedValue("vintfa",rho1) 124 rho.setTaggedValue("vintfb",rho2) 125 126 ahallam 3055 ##########################################################ESTABLISH PDE 127 mypde=LinearPDE(domain) # create pde 128 mypde.setSymmetryOn() # turn symmetry on 129 # turn lumping on for more efficient solving 130 ahallam 3075 #mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) 131 ahallam 3055 kmat = kronecker(domain) # create the kronecker delta function of the domain 132 ahallam 3069 mypde.setValue(D=rho*kmat) #set the general form value D 133 ahallam 3055 134 ############################################FIRST TIME STEPS AND SOURCE 135 # define small radius around point xc 136 ahallam 3089 src_length = 20; print "src_length = ",src_length 137 ahallam 3055 # set initial values for first two time steps with source terms 138 ahallam 3069 xb=FunctionOnBoundary(domain).getX() 139 ahallam 3089 #sy=source[0]*(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-src_length)) 140 y=Vector(0.0,FunctionOnBoundary(domain)) 141 142 ahallam 3069 src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down 143 ahallam 3089 144 #sy=sy*src_dir 145 #sy.setTaggedValue("stop") 146 y.setTaggedValue("stop",src_dir*source[0])#)*(cos(length(xb-xc)*3.1415/src_length)+1)) 147 ahallam 3055 mypde.setValue(y=y) #set the source as a function on the boundary 148 # initial value of displacement at point source is constant (U0=0.01) 149 # for first two time steps 150 ahallam 3089 u=[0.0,0.0,0.0]*wherePositive(x) 151 ahallam 3055 u_m1=u 152 153 ####################################################ITERATION VARIABLES 154 n=0 # iteration counter 155 t=0 # time counter 156 ##############################################################ITERATION 157 while t= rtime): 167 ahallam 3069 saveVTK(os.path.join(savepath,"ex09.%05d.vtu"%n),displacement=length(u),\ 168 ahallam 3055 acceleration=length(accel),tensor=stress) 169 rtime=rtime+rtime_inc #increment data save time 170 # increment loop values 171 t=t+h; n=n+1 172 if (n < ls): 173 ahallam 3089 y.setTaggedValue("stop",source[n]*src_dir) 174 mypde.setValue(y=y) #set the source as a function on the boundary 175 ahallam 3055 print n,"-th time step t ",t