forked from hunterakins/pynm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmesh_routines.py
303 lines (265 loc) · 10.5 KB
/
mesh_routines.py
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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
"""
Description:
Functions to get the discretized wave equation on a grid
Date:
4/18/2023 (copied from sturm_seq)
Author: Hunter Akins
Institution: Scripps Institution of Oceanography, UC San Diego
Copyright (C) 2023 F. Hunter Akins
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import rc
rc('text', usetex=True)
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'STIXGeneral'
import numba as nb
from numba import njit
@njit(nb.f8[:](nb.f8, nb.f8, nb.f8[:], nb.f8[:], nb.f8))
def get_a(omega, h, c, rho, h0):
"""
Compute the diagonal elements a used for depths within a layer
For the pseudo-Sturm Liouville problem
I use the layer mesh scaling so that multiple layers can be used with different
meshes and operate on the same eigenvalue kr^2 h0^2
which results in
d_{i} = -2 h_{0}^{2}/h_{i}^{2} + h_{0}^{2}(\omega^{2} / c^{2}(z_{i})) - \lambda
\implies a_{i} = -2 h_{0}^{2}/h_{i}^{2} + h_{0}^{2}(\omega^{2}
The first and last entries of c and rho are the values
at the layers. The interface value is handled separately,
so the first depth returned here is z[1] and the last is z[-2]
To accomplish that I compute for every detph and discard first and last pt
I compute the value for the interface depth using a_last or a_bdry
Input
omega - float
source freq (rad/s)
c - np 1d array
ssp vals at depths z0, z1, ..., zN-1
as mentioned, depths z0 and zN-1 are interface depths
rho - np 1d array
density at depths z0, z1, ..., zN-1
h0 - float
spacing of first mesh grid
Output -
avec - np 1d array
diagonal values of matrix
"""
avec = (-2.0*h0*h0/h/h + h0*h0*omega*omega / c/c)
return avec[1:-1]
@njit(nb.f8(nb.f8, nb.f8, nb.f8, nb.f8, nb.f8, nb.f8, nb.f8, nb.f8))
def get_a_last(omega, h, cu, rhou, cb, rhob, h0, lam):
"""
Enforce halfspace boundary condition
Eqn. 5.110 but multiplied by 2*h*rho * h0^2 / h^2
fb(kr^2) = 1
gb(kr^2) = rhob / sqrt(kr^2- (w/cb)^2) 5.62
a_{N} = $ -h_{0}^{2} / h^{2} + \frac{1}{2}h_{0}^{2} \omega^{2} / c_{u}^{2}
-\frac{h_{0}^{2} \rho_{u} h}{h^{2} \rho_{b}} \gamma_{b}^{2} $
"""
kr_sq = lam /(h0*h0)
if kr_sq < np.square(omega / cb):
raise ValueError('cmax not set properly (mode is not exponentially decaying in halfspace). kr_sq < (omega / cb)^2 ({0} < {1}')
gamma_b = np.sqrt(kr_sq - np.square(omega/cb))
term1 = -2*h0*h0 / (h*h)
term2 = h0*h0*omega*omega / cu /cu
term3 = -gamma_b * 2* h0*h0 * rhou / ( h * rhob )
alast = term1 + term2 + term3
return alast
@njit(nb.f8(nb.f8, nb.f8, nb.f8, nb.f8, nb.f8, nb.f8, nb.f8, nb.f8))
def get_a_bdry(omega, hu, hb, cu, rhou, cb, rhob, h0):
"""
Get center term for finite difference approximation to layer condition
Equation 5.131 multiplied through by alpha = h0^2 (h_u / (2 \rho_u) + h_b / 2 \rho_b)^{-1}
omega - float
hu - float
mesh grid stepsize in above layer
hb - float
mesh grid stepsize in below layer
cu - float
sound speed above
rhou
density above
cb - float
speed blow
rhob - float
density below
h0 - float
reference meshgrid step size for A
"""
om_sq = np.square(omega)
cu_sq = np.square(cu)
cb_sq = np.square(cb)
alpha = h0*h0 / ((hu/2/rhou + hb/2/rhob))
term1 = -1/(hu*rhou)
term2 = -1/(hb*rhob)
term3 = .5*hu*om_sq / cu_sq / rhou
term4 = .5*hb*om_sq / cb_sq / rhob
a = alpha*(term1 + term2 + term3 + term4)
return a
def get_A_size(z_list):
""" Each layer in z_list shares an interface point
Since I don't use the first depth (pressure release) in z_list[0],
just counting z_list[i][1:] will include all depths needed (including
last boundary condition)
"""
a_size = 0
for x in z_list:
a_size += x.size -1
return a_size
@njit
def get_A_size_numba(z_arr, ind_arr):
a_size = z_arr.size - ind_arr.size # subtract
return a_size
def get_A(omega, h_list, z_list, c_list, rho_list, c_hs, rho_hs, lam):
"""
Compute diagonal and off diagonal terms for the matrix
required in the Sturm sequence solution method
The way the indexing works:
a[i] is the ith diagonal element
d[i] is the element of the matrix directly beneath a[i]
e[i] is the element of the matrix directly to the right of a[i]
omega - float
source frequency (rad/s)
h_list - list of floats
mesh width for each layer
z_list - list of 1d numpy ndarrays
each element is the depths of the c and rho vals
c_list - list of 1d numpy ndarrays
each element is the discretized SSP array for each layer
rho_list - list of 1d numpy ndarrays
elements are discretized density for each layer
c_hs - float
halfspace speed
rho_hs - float
halfspace density
"""
num_layers = len(h_list)
h0 = h_list[0]
a_size = get_A_size(z_list)
a_diag = np.zeros((a_size))
e1 = np.zeros((a_size)-1) # upper off-diagonal
d1 = np.zeros((a_size)-1) # lower off-diagonal
upper_layer_ind = 0 # index of upper layer (entry in diag a)
for i in range(num_layers):
"""
First compute the diagonal terms
"""
z = z_list[i]
c = c_list[i]
h = h_list[i]
rho = rho_list[i]
"""
Fill z.size - 2 entries corresponding to depths starting below the upper layer
interface and going down to the grid point above the bottom of the layer interface
Then add in line for boundary value
"""
a_layer = get_a(omega, h, c, rho, h0) # remember this contains no interface points
a_inds = (upper_layer_ind, upper_layer_ind + z.size-2) # z includes the interface points, so exclude those...
a_diag[a_inds[0]:a_inds[1]] = a_layer[:]
""" Now add the final entry for the interface beneath the layer """
hu = h
rhou = rho[-1]
cu = c[-1]
""" If it's not the bottom halfspace """
if i < num_layers - 1:
hb = h_list[i+1]
rhob = rho_list[i+1][0]
cb = c_list[i+1][0]
a_bdry = get_a_bdry(omega, hu, hb, cu, rhou, cb, rhob, h0)
else: # last layer
rhob =rho_hs
cb = c_hs
a_bdry = get_a_last(omega, hu, cu, rhou, cb, rhob, h0, lam)
a_diag[a_inds[1]] = a_bdry
""" Now compute off diagonal terms """
e1[a_inds[0]:a_inds[1]] = h0*h0 / (h*h) # above diag
d1[a_inds[0]:a_inds[1]-1] = h0*h0 / (h*h) # below diag
if i < num_layers - 1: # if not on the last layer, e1
alpha = h0*h0 / (.5 * (hu / rhou + hb / rhob))
f = alpha / (hu*rhou)
d1[a_inds[1]-1] = f
d1[a_inds[1]] = h0*h0 / (hb*hb) # this point extends into the next layer
g = alpha / (hb*rhob)
e1[a_inds[1]] = g
else: # for bottom boundary condition
d1[a_inds[1]-1] = 2*h0*h0 / (h*h)
upper_layer_ind += z.size - 1
return a_diag, e1, d1
@njit
def get_A_numba(omega, h_arr, ind_arr, z_arr, c_arr, rho_arr, c_hs, rho_hs, lam):
"""
the arrays are the concatenated list elements in the equivalent function above
ind_arr gives the index of the ith layer
so the first element is always zero
z_arr[ind_arr[i]] is the first value in the
"""
num_layers = h_arr.size
h0 = h_arr[0]
#a_size = get_A_size(z_list)
a_size = z_arr.size - ind_arr.size # subtract
a_diag = np.zeros((a_size))
e1 = np.zeros((a_size)-1) # upper off-diagonal
d1 = np.zeros((a_size)-1) # lower off-diagonal
upper_layer_ind = 0 # index of upper layer (entry in diag a)
for i in range(num_layers):
"""
First compute the diagonal terms
"""
if i < num_layers-1:
z = z_arr[ind_arr[i]:ind_arr[i+1]]
c = c_arr[ind_arr[i]:ind_arr[i+1]]
h = h_arr[i]
rho = rho_arr[ind_arr[i]:ind_arr[i+1]]
else:
z = z_arr[ind_arr[i]:]
c = c_arr[ind_arr[i]:]
h = h_arr[i]
rho = rho_arr[ind_arr[i]:]
"""
Fill z.size - 2 entries corresponding to depths starting below the upper layer
interface and going down to the grid point above the bottom of the layer interface
Then add in line for boundary value
"""
a_layer = get_a(omega, h, c, rho, h0) # remember this contains no interface points
a_inds = (upper_layer_ind, upper_layer_ind + z.size-2) # z includes the interface points, so exclude those...
a_diag[a_inds[0]:a_inds[1]] = a_layer[:]
""" Now add the final entry for the interface beneath the layer """
hu = h
rhou = rho[-1]
cu = c[-1]
""" If it's not the bottom halfspace """
if i < num_layers - 1:
hb = h_arr[i+1]
rhob = rho_arr[ind_arr[i+1]]
cb = c_arr[ind_arr[i+1]]
a_bdry = get_a_bdry(omega, hu, hb, cu, rhou, cb, rhob, h0)
else: # last layer
rhob =rho_hs
cb = c_hs
a_bdry = get_a_last(omega, hu, cu, rhou, cb, rhob, h0, lam)
a_diag[a_inds[1]] = a_bdry
""" Now compute off diagonal terms """
e1[a_inds[0]:a_inds[1]] = h0*h0 / (h*h) # above diag
d1[a_inds[0]:a_inds[1]-1] = h0*h0 / (h*h) # below diag
if i < num_layers - 1: # if not on the last layer, e1
alpha = h0*h0 / (.5 * (hu / rhou + hb / rhob))
f = alpha / (hu*rhou)
d1[a_inds[1]-1] = f
d1[a_inds[1]] = h0*h0 / (hb*hb) # this point extends into the next layer
g = alpha / (hb*rhob)
e1[a_inds[1]] = g
else: # for bottom boundary condition
d1[a_inds[1]-1] = 2*h0*h0 / (h*h)
upper_layer_ind += z.size - 1
return a_diag, e1, d1