examples/cookbook/example09a.py


########################################################
#
# Copyright (c) 2009-2010 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2009-2010 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

############################################################FILE HEADER
# example09.py
# Antony Hallam
# Seismic Wave Equation Simulation using acceleration solution.
# 3D model with multiple layers.

#######################################################EXTERNAL MODULES
from esys.escript import *
from esys.finley import Rectangle
import os
# smoothing operator 
from esys.escript.pdetools import Projector, Locator
from esys.escript.unitsSI import *
import numpy as np
import matplotlib
matplotlib.use('agg') #It's just here for automated testing

import pylab as pl
import matplotlib.cm as cm
from esys.escript.linearPDEs import LinearPDE
from esys.finley import ReadMesh
from esys.weipa import saveVTK

########################################################MPI WORLD CHECK
if getMPISizeWorld() > 1:
    import sys
    print("This example will not run in an MPI world.")
    sys.exit(0)

#################################################ESTABLISHING VARIABLES
# where to save output data
savepath = "data/example09a"
meshpath = "data/example09m"
mkDir(savepath)
#Geometric and material property related variables.
step=10.0 # the element size

vel2=1800.;   vel1=3000.
rho2=2300.;   rho1=3100. #density
mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus
lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant

####################################################TESTING SWITCH
testing=True
if testing:
    print('The testing end time is currently selected. This severely limits the number of time iterations.')
    print("Try changing testing to False for more iterations.")
    tend=0.001
    #Model Parameters    
    mx=40.
    my=40.
    mz=20.
    outputs=5
else:
    tend=0.1    # end time
    #Model Parameters
    mx=100.0   #x width of model
    my=100.0   #y width of model
    mz=50.0   #depth of model
    outputs=200

####################################################TIME RELATED VARIABLES 
h=0.00005    # time step
# data recording times
rtime=0.0 # first time to record
rtime_inc=tend/outputs # time increment to record
#Check to make sure number of time steps is not too large.
print("Time step size= ",h, "Expected number of outputs= ",tend/h)

####################################################CREATING THE SOURCE FUNCTION
U0=0.1 # amplitude of point source
ls=500   # length of the source
source=np.zeros(ls,'float') # source array
decay1=np.zeros(ls,'float') # decay curve one
decay2=np.zeros(ls,'float') # decay curve two
time=np.zeros(ls,'float')   # time values
g=np.log(0.01)/ls

dfeq=50 #Dominant Frequency
a = 2.0 * (np.pi * dfeq)**2.0
t0 = 5.0 / (2.0 * np.pi * dfeq)
srclength = 5. * t0
ls = int(srclength/h)
print('source length',ls)
source=np.zeros(ls,'float') # source array
ampmax=0
for it in range(0,ls):
    t = it*h
    tt = t-t0
    dum1 = np.exp(-a * tt * tt)
    source[it] = -2. * a * tt * dum1
    if (abs(source[it]) > ampmax):
        ampmax = abs(source[it])
    time[t]=t*h

# will introduce a spherical source at middle left of bottom face
xc=[mx/2,my/2,0]

####################################################DOMAIN CONSTRUCTION
domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
x=domain.getX() # get the locations of the nodes in the domain

lam=Scalar(0,Function(domain))
lam.setTaggedValue("vintfa",lam1)
lam.setTaggedValue("vintfb",lam2)
mu=Scalar(0,Function(domain))
mu.setTaggedValue("vintfa",mu1)
mu.setTaggedValue("vintfb",mu2)
rho=Scalar(0,Function(domain))
rho.setTaggedValue("vintfa",rho1)
rho.setTaggedValue("vintfb",rho2)

##########################################################ESTABLISH PDE
mypde=LinearPDE(domain) # create pde
mypde.setSymmetryOn() # turn symmetry on
# turn lumping on for more efficient solving
#mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().HRZ_LUMPING)
kmat = kronecker(domain) # create the kronecker delta function of the domain
mypde.setValue(D=rho*kmat) #set the general form value D

############################################FIRST TIME STEPS AND SOURCE
# define small radius around point xc
src_rad = 20; print("sourc radius= ",src_rad)
# set initial values for first two time steps with source terms
xb=FunctionOnBoundary(domain).getX()
yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)
stop=Scalar(0.0,FunctionOnBoundary(domain))
stop.setTaggedValue("stop",1.0)
src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down

mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
# initial value of displacement at point source is constant (U0=0.01)
# for first two time steps
u=[0.0,0.0,0.0]*wherePositive(x)
u_m1=u

####################################################ITERATION VARIABLES
n=0 # iteration counter
t=0 # time counter
##############################################################ITERATION
while t<tend:
    # get current stress
    g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
    mypde.setValue(X=-stress) # set PDE values
    accel = mypde.getSolution() #get PDE solution for accelleration
    u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
    u_p1=u_p1#*abc          # apply boundary conditions
    u_m1=u; u=u_p1 # shift values by 1
    # save current displacement, acceleration and pressure
    if (t >= rtime):
        saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
                                    acceleration=length(accel),tensor=stress)
        rtime=rtime+rtime_inc #increment data save time
    # increment loop values
    t=t+h; n=n+1
    if (n < ls):
        mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
    print(n,"-th time step t ",t)
1
2	########################################################
3	#
4	# Copyright (c) 2009-2010 by University of Queensland
5	# Earth Systems Science Computational Center (ESSCC)
6	# http://www.uq.edu.au/esscc
7	#
8	# Primary Business: Queensland, Australia
9	# Licensed under the Open Software License version 3.0
10	# http://www.opensource.org/licenses/osl-3.0.php
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	http://www.opensource.org/licenses/osl-3.0.php"""
20	__url__="https://launchpad.net/escript-finley"
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]yxsrc_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<tend:
163	# get current stress
164	g=grad(u); stress=lamtrace(g)kmat+mu(g+transpose(g))#abc
165	mypde.setValue(X=-stress) # set PDE values
166	accel = mypde.getSolution() #get PDE solution for accelleration
167	u_p1=(2.u-u_m1)+hh*accel #calculate displacement
168	u_p1=u_p1#*abc # apply boundary conditions
169	u_m1=u; u=u_p1 # shift values by 1
170	# save current displacement, acceleration and pressure
171	if (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]yxsrc_dir*stop) #set the source as a function on the boundary
179	print(n,"-th time step t ",t)