Add files via upload

calculateGradient.m

@@ -0,0 +1,56 @@
+%% 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
+

calculateStepLength.m

@@ -0,0 +1,32 @@
+%% This script is used to perform calculation on the steplength
+function [stepLength,newErr] = calculateStepLength(dims,model,gradient,oldErr,source,trueRec,unext,u,uprev)
+%% Depending on how you solve the problem you most probably need more input variables                      
+    recording = zeros(dims.nt,length(dims.recPos),'single');
+    stepLength = 256;
+    newErr = inf;
+    while (newErr > oldErr)
+        newErr = 0;
+        stepLength = stepLength/2;
+        %% Test model update calculation based on steplength and gradient
+        modelnew = model+stepLength*gradient;
+        %% Solve wave equation using test model update
+        for s = 1:dims.ds:length(dims.srcPos)
+            uprev(:)=0; u(:)=0; unext(:)=0;
+            for t = 1:dims.nt
+                unext = solveWaveEqn(modelnew,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);
+            end
+            %% Calculate new error and check against old
+            chi    = recording-trueRec(:,:,s);
+            newErr = newErr+norm(chi);
+        end
+    end
+end

runFWIGoldenRatio.m

@@ -0,0 +1,85 @@
+%% This script is used for running the Full Waveform Inversion method.
+% We're plotting the model & gradient in  this script.
+%% Setting up dimensions
+dims.dy =     10; % [m]
+dims.dx =     10; % [m]
+dims.dt = 1.0e-3; % [s]
+
+dims.ny = 201; % Cells in y-direction
+dims.nx = 301; % Cells in x-direction
+dims.nt = 801; % Amount of time steps
+
+%% Model dimensions
+dims.modely = 100:150;
+dims.modelx = 100:200;
+dims.my = length(dims.modely);
+dims.mx = length(dims.modelx);
+
+%% Source locations
+sx = min(dims.modelx):max(dims.modelx);
+sy = min(dims.modely)*ones(1,length(sx));
+dims.srcPos = sy + dims.ny*sx;
+
+%% Receiver locations
+rx = min(dims.modelx):max(dims.modelx);
+ry = min(dims.modely)*ones(1,length(rx));
+dims.recPos = ry+dims.ny*rx;
+
+%% Creating background model
+bg = zeros(dims.ny,dims.nx,'single');
+bg(:) = 2.0e3;         % [m/s] - Background
+bg(115:end,:) = 2.3e3; % [m/s] - Layer
+
+%% Begin iteration
+model = bg;     % Starting model
+dims.ds = 5;   % Grid point distance between sources
+maxIter = 10;   % Maximum number of iterations per frequency
+freqs = [4,6,8,10,12];  % Frequencies to use in inversion
+
+%% Initial value of u(x,t)
+unext = zeros(dims.ny,dims.nx,'single');
+u = zeros(dims.ny,dims.nx,'single');
+uprev = zeros(dims.ny,dims.nx,'single');
+
+errVec = zeros(1,maxIter*length(freqs));
+it = 1; 
+start=tic;
+fprintf('Iteration \t Error \t\t Runtime(s) \n');
+%% FWI loop
+for f = freqs
+    %% Generating ricker source signature wavelet 
+    source = rickerWave(f,dims);
+    %% Load true recording
+    load (['trueRec_',num2str(f),'Hz.mat']);   
+    for i = 1:maxIter
+        %% Calculate gradient
+        [gradient,err] = calculateGradient(dims,source,trueRec,model,unext,u,uprev);
+        %% Taper gradient
+        taperg = taperGradient(gradient);
+        %% Calculate error & step length
+        [stepLength,err] = calculateStepLength(dims,model,taperg,source,trueRec);
+        %% Update model
+        model = model+stepLength*taperg;
+        %% Gradient & Model plotting
+        figure(2)
+        subplot(2,1,1);
+        imagesc(taperg(dims.modely,dims.modelx));
+        colormap(jet);
+        title(['Gradient ',num2str(f),' Hz Plot']);
+        axis('image');
+        colorbar();
+        drawnow();
+        subplot(2,1,2);
+        imagesc(model(dims.modely,dims.modelx));
+        colormap(jet);
+        title(['Model ',num2str(f),' Hz Plot']);
+        axis('image');
+        colorbar();
+        drawnow();
+
+        errVec(it) = err;
+        it = it + 1;
+        runtime=toc(start);
+        fprintf('%i \t\t\t %6.4f \t %6.4f \n',i,err,runtime);
+    end
+end