-
Notifications
You must be signed in to change notification settings - Fork 12
/
calculateGradient.m
56 lines (55 loc) · 2.26 KB
/
calculateGradient.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
%% This script is used to perform calculation on the gradient
function [gradient,err] = calculateGradient(dims,source,trueRec,model,unext,u,uprev)
recording = zeros(dims.nt,length(dims.recPos),'single');
gradient = zeros(dims.ny,dims.nx,'single');
forwardField = zeros(dims.my,dims.mx,dims.nt,'single');
adjointField = zeros(dims.my,dims.mx,dims.nt,'single');
err = 0;
for s = 1:dims.ds:length(dims.srcPos)
%% Run forward simulation on background model
uprev(:)=0; u(:)=0; unext(:)=0;
for t = 1:dims.nt
% Solve wave equation
unext = solveWaveEqn(model,dims,dims.srcPos(s),source,t,unext,u,uprev);
% Update u(x,t)
uprev = u; u = unext;
% Check wave equation stability
r = model*dims.dt/dims.dx + model*dims.dt/dims.dy;
if r>1
break
end
% Record traces
recording(t,:) = u(dims.recPos);
% Save forward field for use in correlation
forwardField(:,:,t) = u(dims.modely,dims.modelx);
end
%% Calculate difference and error
chi = recording-trueRec(:,:,s);
% Time reversal
chi = flipud(chi);
% Error calculation
err = err + norm(chi);
%% Run adjoint simulation
uprev(:)=0; u(:)=0; unext(:)=0;
for t = 1:dims.nt
% Solve wave equation using the difference (chi) as sources
unext = solveWaveEqn(model,dims,dims.recPos,chi,t,unext,u,uprev);
% Update u(x,t)
uprev = u; u = unext;
% Check wave equation stability
r = model*dims.dt/dims.dx + model*dims.dt/dims.dy;
if r>1
break
end
% Save adjoint field for use in correlation
adjointField(:,:,dims.nt-t+1) = u(dims.modely,dims.modelx);
end
%% Correlate
for t = 2:dims.nt-1
% Calculate the time derivative of the displacement to gradient.
dadt=(adjointField(:,:,t+1) - adjointField(:,:,t))/dims.dt;
dfdt=(forwardField(:,:,t+1) - forwardField(:,:,t-1))/dims.dt;
gradient(dims.modely,dims.modelx)=gradient(dims.modely,dims.modelx)+dadt.*dfdt;
end
end
end