-
Notifications
You must be signed in to change notification settings - Fork 4
/
getA.m
148 lines (121 loc) · 3.86 KB
/
getA.m
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
function A = getA(f, m, h, n, fs)
% Create 2D Helmholtz matrix with 2nd order Clayton-Enquist absorbing boundary
%
% use:
% A = getA(f,m,h,n,fs)
%
% input:
% f - frequency [Hz]
% m - squared-slownes [s^2/km^2]
% h - gridspacing in each direction d = [d1, d2];
% n - number of gridpoints in each direction n = [n1, n2]
% fs - free surface flag, bool, true or false
%
% output:
% A - sparse matrix
%
% Based on https://github.com/TristanvanLeeuwen/SimpleFWI
%
% Vladimir Kazei, Oleg Ovcharenko, 2019
%%
defval('fs',false);
% scale frequency to handle velocity in km/s
omega = 1e-3*2*pi*f;
k = omega.*sqrt(m);
% z derivative operator
D1 = spdiags(ones(n(1),1)*[1 -1]/h(1),[0:1],n(1)-1,n(1));
% x derivative operator
D2 = spdiags(ones(n(2),1)*[1 -1]/h(2),[0:1],n(2)-1,n(2));
% internal points of the domain
spy2 = speye(n(2));
% zeroing out the boundary
spy2([1 end],[1 end])=0;
% 2nd derivative in z direction
L11 = -kron(spy2,D1'*D1);
spy1 = speye(n(1));
spy1([1 end],[1 end])=0;
% 2nd derivative in x direction
L22 = -kron(D2'*D2,spy1);
% Laplacian inside the domain
L = L11 + L22;
% boundary
L1 = -kron(speye(n(2)),D1'*D1);
L2 = -kron(D2'*D2,speye(n(1)));
% 2nd order derivative along the boundary
L_BOUND = L2 + L1 - L11 - L22;
% this is a trick to have half Helmholtz in the second variable
% Matrix a is not equal to 1 at the boundary
a = ones(n);
a(:,[1 end]) = .5; a([1 end],:) = .5; a = a(:);
% boundary points of the domain
BND = find(a~=1);
% derivative with respect to the normal of the boundary
L_n = 0*L;
% Matrix a_corners is non-zero at the corners
a_corners = zeros(n);
a_corners(1,1) = 1; % top left
a_corners(1,end) = 1; % top right
a_corners(end,1) = 1; % bottom left
a_corners(end,end) = 1; % bottom right
% find linear indices of corner points
CORNERS = find(a_corners==1);
% adjusting 2nd derivative in the corners (distance is sqrt(2) shorter)
L_BOUND(CORNERS, :) = 2 * L_BOUND(CORNERS, :);
% Normal 1st derivatives at the corners
L_corner = 0*L;
% upper left corner
% compute the first derivative
L_corner(1,1) = 1;
L_corner(1,n(1)+2) = -1;
% distribute the second derivative
L_BOUND(1,1) = L_BOUND(1,1)/2;
L_BOUND(1,n(1)+2) = L_BOUND(1,1);
% lower left corner
% compute the first derivative
L_corner(n(1),n(1)) = 1;
L_corner(n(1),2*n(1) -1) = -1;
% distribute the second derivative
L_BOUND(n(1),n(1)) = L_BOUND(n(1),n(1))/2;
L_BOUND(n(1),2*n(1) -1) = L_BOUND(n(1),n(1));
% lower right corner
% compute the first derivative
L_corner(prod(n),prod(n)) = 1;
L_corner(prod(n),prod(n)-n(1)-1) = -1;
% distribute the second derivative
L_BOUND(prod(n),prod(n)) = L_BOUND(prod(n),prod(n))/2;
L_BOUND(prod(n),prod(n)-n(1)-1) = L_BOUND(prod(n),prod(n));
% upper right corner
% compute the first derivative
L_corner(prod(n)-n(1)+1, prod(n)-n(1)+1) = 1;
L_corner(prod(n)-n(1)+1, prod(n)-2*n(1)+2) = -1;
% distribute the second derivative
L_BOUND(prod(n)-n(1)+1, prod(n)-n(1)+1) = L_BOUND(prod(n)-n(1)+1, prod(n)-n(1)+1)/2;
L_BOUND(prod(n)-n(1)+1, prod(n)-2*n(1)+2) = L_BOUND(prod(n)-n(1)+1, prod(n)-n(1)+1);
%
L_corner = L_corner/(sqrt(2)*h(1));
% for the boundary points L is working as normal derivative
L_n(BND,:) = L(BND,:);
%% ASSEMBLE HELMHOLTZ MATRIX WITH 2ND ORDER ABC
% inside domain
A = diags(k.^2) + L-L_n;
% Clayton-Enquist 1977
% P_tt + (1/v) P_xt - (v/2) P_zz = 0
% we translate with the rules _t -> -i* omega;
% k^2 P + 1i * k P_x + 1/2 P_zz = 0
% where k = omega/v.
% Then we use prepared derivative operators
% _x -> L_n (normal derivative), _zz -> L_BOUND
A = A - 1i*(k(1))*(L_n * h(1) - L_corner) + 0.5*L_BOUND;
%A = prod(h)*A;
% Free surface
if fs
b = zeros(n); b(1,:) = 1; b=b(:);
FSP = find(b);
A(FSP,:) = 0;
A(:,FSP) = 0;
A = omega^2*diags(b.*m) + A;
end
% check the number of grid points per wavelength
if (f > min(1e3*1./sqrt(m))/(5*h(1)))
warning('Dispersion! f>min(1e3*1./sqrt(m))/(5*h(1))=%e',min(1e3*1./sqrt(m))/(5*h(1)));
end