1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
"""
simple seismic solver in frequency domain
with reflectors and PML.
"""
from esys.escript import *
from esys.finley import Rectangle
from esys.downunder.apps import SonicWaveInFrequencyDomain, PMLCondition
from esys.weipa import saveSilo
from esys.escript.pdetools import Locator
import numpy as np
# horizontal and vertical grid spacing in
dx=10.
dz=10.
# grid size
NEx=400
NEz=180
# element order (1 or 2):
Order=2
# extend of domain:
Width=NEx*dx
Depth=NEz*dz
# position of geophones:
StationOffset = 60*dx # position of first phone from left boundary
StationSpacing = dx # phone spacing
NStations = int((Width-2*StationOffset)//StationSpacing+0.5)+1 # number of stations
# velocity:
Layers=[ 30*dx , 40*dx, 90*dx ] # thickness of layer from the top down
Vs =[ 1500., 2000., 3000. ] # corresponding velocities
Vbase=3500 # base velocity
Frequency = 8. # frequency
Amplitude= 1. # and amplitude
# thickness of PML zone:
L_PML=40*dx
print("Domain = %g x %g m"%(Width, Depth))
print("Grid Cells= %d x %d"%(NEx, NEz))
print("Grid Nodes= %d x %d"%(Order*NEx+1, Order*NEz+1))
print("Grid Spacing = %g x %g m"%(dx, dz))
print("PML zone = %g m"%(L_PML))
print("Number of Stations = %d"%(NStations,))
print("Station spacing = %g m"%(StationSpacing,))
print("Station offset = %g m"%(StationOffset,))
print("Reflector Layers [m]: ", Layers)
print("Velocities [m/s]: ", Vs)
print("Base velocity : [m/s]", Vbase)
print("Frequency = %g Hz "%Frequency)
print("Amplitude = %g"%Amplitude)
# generate domain:
stationPositions = [ (StationOffset+s*StationSpacing, Depth ) for s in range(NStations)]
stationTags= [ "P%4.4d"%s for s in range(NStations)]
domain=Rectangle(NEx,NEz,l0=Width,l1=Depth, diracPoints=stationPositions, diracTags=stationTags, order=Order, fullOrder=True)
print("Domain has been generated.")
Dim=domain.getDim()
# setup velocity field:
vp=Scalar(Vbase, Function(domain))
X=vp.getX()
Top=sup(domain.getX()[Dim-1])
for l in range(len(Vs)-1, -1, -1 ):
mask=wherePositive(X[Dim-1]-(Top-Layers[l]))
vp=vp*(1-mask)+mask*Vs[l]
print("Velocity field generated:",vp)
# check spacing vs, wavelength:
wavelength=inf(vp/(2*np.pi*Frequency))
print("minimum wavelength = %g"%wavelength)
assert wavelength > max(dx,dz)/Order
# locator to grap amplitudes at geophones:
geophones=Locator(Solution(domain), stationPositions)
# create PML condition:
pml_condition=PMLCondition(sigma0=vp*100/Frequency, Lleft=[L_PML,L_PML], Lright=[L_PML, None], m=4)
print("PMLConditon generated ...")
# sonic wave model
wave=SonicWaveInFrequencyDomain(domain, vp=vp, pml_condition=pml_condition)
wave.setFrequency(Frequency)
print("model created...")
# set source and get wave response:
source=Scalar(0j, DiracDeltaFunctions(domain))
source.setTaggedValue(stationTags[NStations//2], Amplitude)
u=wave.getWave(source)
print("model solved...")
# get amplitude and phase
amps=abs(u)
phase=phase(u)/np.pi*180.
# save to file:
saveSilo("solution", velocity=vp, amplitude=amps, phase=phase, pmlzone=pml_condition.getPMLMask(domain))
print("solution saved to file ...")
# get the amplitude and phase at geophones:
import matplotlib.pyplot as plt
fig=plt.gcf()
ax1=plt.gca()
geophoneX=[ x[0] for x in geophones.getX()]
color = 'tab:red'
ax1.set_xlabel('offset')
ax1.set_ylabel('amplitude', color=color)
ax1.plot(geophoneX , geophones(amps), color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('phase [deg]', color=color) # we already handled the x-label with ax1
ax2.plot(geophoneX , geophones(phase), color=color, marker=".")
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
|
