-
Notifications
You must be signed in to change notification settings - Fork 4
/
example_FDFD.m
156 lines (128 loc) · 5.59 KB
/
example_FDFD.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
149
150
151
152
153
154
155
156
clc;clear
% This example shows how to use FDFD_ver02 package to solve 2D
% visco-acoustic wave equation. First, we need to define model properties
% such as Vp, rho and Q factor models, then we need to set numerical modeling
% parameters.
% I tried to explain this example with appropriate comments. All input and output
% variables of the functions are also explained briefly inside each function.
%
% For more information about theory of this work please see the following paper:
%
% Amini, N. and Javaherian, A., “A MATLAB-based frequency-domain
% finite-difference package for solving 2D visco-acoustic wave equation”,
% Waves in Random and Complex Media, vol. 21, no. 1, pp. 161–183, 2011.
% doi:10.1080/17455030.2010.537708.
%
% Please cite the above paper when reporting, reproducing or extending the results.
%% model properties
% Homogeneous model
nx_org = 200; % horizontal dim of model
nz_org = 200; % vertical dim of model
vp = 2000*ones(nz_org,nx_org); % velocity model
rho = 2000*ones(nz_org,nx_org); % density model
q = 100*ones(nz_org,nx_org); % Q-factor model. For visco-acoustic modeling "atten_opt" should be set to 'KF'.
% Marmousi model
% load('data/Marmousi_vp.mat')
% vp = MARM_vp(1:2:end,1:2:end);
% rho = 310*vp.^0.25; % density model (using Gardner's equation)
% q = vp/15; % Q-factor model (just an assumption!). For visco-acoustic modeling "atten_opt" should be set to 'KF'.
%% modeling parameters
%%% boundary conditions
% hybrid PML+ABC absorbing boundary is used for attenuation of waves reflected from
% boundaries. Absorbing boundary is applied to the left, right and bottom boundaries automatically.
L = 30; % width of PML
alpha = 180; % amplitude of PML damping cosine function
% For top boundary we have 3 options:
% 'PML' : hybrid PML+ABC absorbing boundary
% 'Dirichlet' : fixed boundary
% 'Neumann' : free boundary
top_bc = 'Neumann'; % boundary condition at top of model
%%% source parameters
sname = 'ricker'; % type of source wavelet ('ricker' or 'gaussian')
f0 = 20; % dominant frequency of the source wavelet
%%% attenuation parameters
atten_opt = 'KF'; % attenuation option ('no_atten' for non-attenuation or 'KF' for Kolsky-Futterman)
wref= 5; % reference frequency for Kolsky-Futterman attenuation mechanism
%%% other parameters
fmax = 3*f0; % maximum frequency of the simulation (usually fmax = 3*fdom for 'ricker' and 5*fdom for 'gaussian')
freq_zpad = 0; % number of zeros to pad in the frequency domain to have a smooth waveform in time domain
twrap = 5; % damping to suppress time aliasing
%%% discretization parameters
G = 8; % number of samples per minimum wavelength
tmax = 1; % time of simulation
lambda = vp/fmax;
lambda_min = min(lambda(:));
dx = lambda_min/G;
%%% parallelize over frequencies
use_parfor = 0; % if 0 sequential "for" loop over frequencies and if 1 parallelize over frequencies using "parfor" loop
%% extend grids for boundary condition
vp = ext_pml(vp,L,top_bc);
rho = ext_pml(rho,L,top_bc);
q = ext_pml(q, L,top_bc);
[nz,nx] = size(vp); % extended model sizes
%% Sources and Receivers positions
%%% Sources positions
%Sx = 10:10:nx; % muti sources
%Sz = 35*ones(size(Sx)); % muti sources
Sx = fix(nx/2)+1; % single source
Sz = 70; % single source
%%% Receivers positions
Rx = 1:2:nx; % receivers on some of grids
Rz = 25*ones(size(Rx)); % receivers on some of grids
[Rx,Rz] = meshgrid(1:nx,1:nz); % receivers on all of grids. This case is
% good for generating wave propagation snapshots and visualizing as
% animation. this case requires a large amount of memory for many sources
%% FDFD modeling
[pf_r,w] = fdfd(vp,rho,q,atten_opt,wref,dx,tmax,twrap,...
sname,f0,fmax,L,alpha,top_bc,Rx(:),Rz(:),Sx(:),Sz(:),use_parfor);
%% convert frequency-domain to time-domain
[pt,t] = four2time(pf_r,tmax,twrap,freq_zpad);
dt = t(2) - t(1);
%% display options
%%% display wave propagation animation
shot_id = 1; % index number of source
if size(pt,1) == nz*nx % check if receivers are on all of grids
figure
for i = 1:length(t)
tmp = reshape(pt(:,i,shot_id),nz,nx);
imagesc(0:dx:(nx-1)*dx,0:dx:(nz-1)*dx,tmp)
colorbar
axis equal
axis tight
%caxis([-5 5]*2000) % for 'gaussian' source
caxis([-5 5]) % for 'ricker' source
title(['shot number = ' , num2str(shot_id) ...
' , time (s) = ' , num2str((i-1)*dt)])
xlabel('x (m)'), ylabel('z (m)')
pause(0.1)
end
else
warning('you need put to receivers at all of grids for animation!')
end
%%% display shot records
shot_id = 1; % index number of source (can be 1 to length(sx))
figure
if size(pt,1) == nz*nx % check if receivers are on all of grids
% in this case to display shot record we need to define locations
% of the receivers
Rx = 1:1:nx; % horizontal index of receivers
Rz = (5)*ones(size(Rx)); % vertical index of receivers
rec_shot_ind = sub2ind([nz,nx],Rz,Rx);
shot = squeeze(pt(rec_shot_ind,:,shot_id))';
imagesc((Rx-1)*dx,t,-shot(:,:))
xlabel('x (m)'), ylabel('t (s)')
colorbar
%caxis([-5 5]*2000) % for 'gaussian' source
caxis([-5 5]) % for 'ricker' source
grid on
title(['shot number = ' , num2str(shot_id)])
else % we already have defined locations of the receivers
shot = squeeze(pt(:,:,shot_id))';
imagesc((Rx-1)*dx,t,-shot(:,:))
xlabel('x (m)'), ylabel('t (s)')
colorbar
%caxis([-5 5]*2000) % for 'gaussian' source
caxis([-5 5]) % for 'ricker' source
grid on
title(['shot number = ' , num2str(shot_id)])
end