Skip to content
Snippets Groups Projects
Commit be315f3b authored by Jan Thorbecke's avatar Jan Thorbecke
Browse files

added cylinder comparison example

parent 5807e882
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fdelmodc/applySource.c
+ 1
− 1
View file @ be315f3b
......@@ -114,6 +114,7 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
if (verbose>=4 && itime==0) {
vmess("Source %d positioned at grid ix=%d iz=%d",isrc, ix, iz);
}
/* cosine squared windowing to reduce edge effects on shot arrays */
if ( (src.n>1) && src.window) {
scl = 1.0;
......@@ -139,7 +140,6 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
vmess("Source %d at grid [ix=%d,iz=%d] at itime %d has value %e",isrc, ix,iz, itime, src_ampl);
}
/* Force source */
if (src.type == 6) {
......
fdelmodc/demo/matlab/FD_mod_grid.m 0 → 100644
+ 138
− 0
View file @ be315f3b
function [ Ptsct, Pfinc, Pfsct, f_out] = FD_mod_grid( xS, xR, ntrcv, dtrcv, dx, cgrid, rhogrid, f )
%Summary of this function goes here
% Detailed explanation goes here
% save Velocity and density grid
dims = size(cgrid);
fileID =fopen('mod_cp.bin','w+','l');
fwrite(fileID,cgrid,'float32');
fclose(fileID);
fileID =fopen('mod_ro.bin','w+','l');
fwrite(fileID,rhogrid,'float32');
fclose(fileID);
% Compute sizes for makemod
sizez=(dims(1)-1)*dx;
sizex=(dims(2)-1)*dx;
origz=-(dims(1)-1)*dx/2
origx=-(dims(2)-1)*dx/2
zrcv1=xR(1)-origz;
zrcv2=xR(1)-origz;
xrcv1=0.0;
xrcv2=xrcv1+(dims(2)-1)*dx;
zsrc=xS(1)-origz;
xsrc=xS(2)-origx;
tmod=(ntrcv-1)*dtrcv;
%compute dt for modeling dt < 0.606*dx/Cmax
Cmax=max(cgrid(:));
dxmod=dx;
dtmax=0.606*dxmod/Cmax;
dtmod=dtrcv/(ceil(dtrcv/dtmax));
ntwave=16384;
ntfft=4*ntrcv;
fileID = fopen('run.scr','w+');
fprintf(fileID,'#!/bin/bash\n');
fprintf(fileID,'export PATH=$HOME/src/OpenSource/bin:/opt/CWP/bin/:.:$PATH\n');
fprintf(fileID,'which fdelmodc\n');
fprintf(fileID,'which surange\n');
fprintf(fileID,'set -x\n');
fprintf(fileID,'df=0.25\n');
fprintf(fileID,'dt=%e\n',dtmod);
fprintf(fileID,'suaddhead < mod_ro.bin ntrpr=%d ns=%d | \\\n',dims(2), dims(1));
fprintf(fileID,'sushw key=d1,d2,gx,scalco a=%f,%f,%d,-1000 b=0,0,%d,0 > mod_ro.su\n',dx, dx, int32(xrcv1*1000), int32(dx*1000));
fprintf(fileID,'suaddhead < mod_cp.bin ntrpr=%d ns=%d | \\\n',dims(2), dims(1));
fprintf(fileID,'sushw key=d1,d2,gx,scalco a=%f,%f,%d,-1000 b=0,0,%d,0 > mod_cp.su\n',dx, dx, int32(xrcv1*1000), int32(dx*1000));
fprintf(fileID,'makewave w=fw fmin=0 flef=6 frig=94 fmax=100 dt=$dt file_out=wavefw.su nt=%d t0=0.4 scale=0 scfft=1 verbose=1\n', ntwave);
fprintf(fileID,'makewave w=fw fmin=0 flef=6 frig=94 fmax=100 dt=%f file_out=wavefwdt.su nt=%d t0=0.4 scale=0 scfft=1 verbose=1\n', dtrcv, ntfft);
fprintf(fileID,'export filecp=mod_cp.su\n');
fprintf(fileID,'export filero=mod_ro.su\n');
fprintf(fileID,'export OMP_NUM_THREADS=2\n');
fprintf(fileID,'fdelmodc \\\n');
fprintf(fileID,'file_cp=$filecp file_den=$filero \\\n');
fprintf(fileID,'ischeme=1 \\\n');
fprintf(fileID,'file_src=wavefw.su verbose=1 \\\n');
fprintf(fileID,'file_rcv=recvgrid.su \\\n');
fprintf(fileID,'rec_type_vz=0 rec_type_p=1 rec_int_vz=2 \\\n');
fprintf(fileID,'xrcv1=%d xrcv2=%d zrcv1=%d zrcv2=%d \\\n',xrcv1,xrcv2,zrcv1,zrcv2);
fprintf(fileID,'dxrcv=%d \\\n', dx);
fprintf(fileID,'dtrcv=%e \\\n', dtrcv);
fprintf(fileID,'xsrc=%d zsrc=%d\\\n', xsrc, zsrc);
fprintf(fileID,'src_type=1 tmod=%e \\\n', tmod);
fprintf(fileID,'ntaper=100 \\\n');
fprintf(fileID,'left=2 right=2 bottom=2 top=2\n\n');
fprintf(fileID,'\n');
fprintf(fileID,'makemod file_base=hom.su \\\n');
fprintf(fileID,'cp0=1500 ro0=3000 sizex=%d sizez=%d dx=%.1f dz=%.1f orig=0,0 verbose=1\n', sizex,sizez, dxmod,dxmod);
fprintf(fileID,'export filecp=hom_cp.su\n');
fprintf(fileID,'export filero=hom_ro.su\n');
fprintf(fileID,'fdelmodc \\\n');
fprintf(fileID,'file_cp=$filecp file_den=$filero \\\n');
fprintf(fileID,'ischeme=1 \\\n');
fprintf(fileID,'file_src=wavefw.su verbose=1 \\\n');
fprintf(fileID,'file_rcv=recvhom.su \\\n');
fprintf(fileID,'rec_type_vz=0 rec_type_p=1 \\\n');
fprintf(fileID,'xrcv1=%d xrcv2=%d zrcv1=%d zrcv2=%d \\\n',xrcv1,xrcv2,zrcv1,zrcv2);
fprintf(fileID,'dxrcv=%d \\\n', dx);
fprintf(fileID,'dtrcv=%e \\\n', dtrcv);
fprintf(fileID,'xsrc=%d zsrc=%d\\\n', xsrc, zsrc);
fprintf(fileID,'src_type=1 tmod=%e \\\n', tmod);
fprintf(fileID,'ntaper=100 \\\n');
fprintf(fileID,'left=2 right=2 bottom=2 top=2\n\n');
fprintf(fileID,'suop2 recvgrid_rp.su recvhom_rp.su op=diff > recv_rp.su\n');
fprintf(fileID,'sustrip < recv_rp.su > recv_rp.bin\n');
fprintf(fileID,'sustrip < recvhom_rp.su > recvhom_rp.bin\n');
fprintf(fileID,'sustrip < wavefwdt.su > wavefwdt.bin\n');
fclose(fileID);
!chmod +x run.scr
system('./run.scr');
path = getenv('PATH');
path = [path ':$HOME/src/OpenSource/bin:/opt/CWP/bin/:.:'];
setenv('PATH', path);
% Scattered field
ns=ntrcv;
ntr=dims(2);
file='recv_rp.bin';
fid=fopen(file,'r');
Ptsct=fread(fid,[ns,ntr],'float32');
fclose(fid);
% Direct field
file='recvhom_rp.bin';
fid=fopen(file,'r');
Ptinc=fread(fid,[ns,ntr],'float32');
fclose(fid);
df=double(1.0/(ntfft*dtrcv));
a=round(f/df)+1;
f_out=(a-1)*df
fprintf('selected discrete frequency P %f\n', (a-1)*df)
% f2=a*df %these are the selected discrete frequencies
Pfft=fft(Ptsct,ntfft);
Pfsct=Pfft(a,:);
%
Pfft=fft(Ptinc,ntfft);
Pfinc=Pfft(a,:);
ns=ntfft;
ntr=1;
file='wavefwdt.bin';
fid=fopen(file,'r');
Wt=fread(fid,[ns,ntr],'float32');
fclose(fid);
df=double(1.0/(ntfft*dtrcv));
c=round(f/df)+1; % select frequencies as close as the one's in Pf
fprintf('selected discrete frequency W %f\n', (c-1)*df)
Wfft=fft(Wt,ntfft);
Wf=Wfft(c,1);
Pfsct = (Pfsct)./(1.0*Wf); % deconvolve for the wavelet
Pfinc = (Pfinc)./(1.0*Wf); % deconvolve for the wavelet
end
fdelmodc/demo/matlab/ForwardCircle.m 0 → 100644
+ 152
− 0
View file @ be315f3b
function [P_inc,P_sct,xR] = ForwardCircle( xS, xR, gsize, dxR, c_0, c_s, rho_0, rho_s, f )
% xS = source position
% dxR = receiver grid size
% c_0 = wave speed in embedding
% c_s = wave speed of scattering object
% rho_0 = mass density in embedding
% rho_s = mass density of scattering object
% f = temporal frequency;
wavelength = c_0 / f; % wavelength
s = 1e-16 + 1i*2*pi*f; % LaPlace parameter
gam_0 = s/c_0; % propagation coefficient 0
gam_s = s/c_s; % propagation coefficient s
a = 40; % radius circle cylinder
depth = 20; % distance betwen circle and origin
disp(['wavelength = ' num2str(wavelength)]);
% Thorbecke coordinates
xSource1 = xS(1); % source position
xSource2 = xS(2);
xReceiver1 = xR(1); % receiver positions
nr=(gsize(2)-1)/2;
xReceiver2 = (-nr:nr) * dxR;
NR = size(xReceiver2,2);
% Van den Berg coordinates
xS(1) = xSource1;
xS(2) = xSource2;
xR(1,1:NR) = xReceiver1;
xR(2,1:NR) = xReceiver2(1:NR);
% Compute 2D incident field ---------------------------------------------
DIS = sqrt( (xR(1,:)-xS(1)).^2 + (xR(2,:)-xS(2)).^2 );
p_inc = 1/(2*pi).* besselk(0,gam_0*DIS);
% multiply by rho_0 and plot
P_inc = rho_0 * p_inc ; % Take care of factor \rho
% Make grid
N1 = gsize(1); % number of samples in x_1
N2 = gsize(2); % number of samples in x_2
dx = dxR; % with meshsize dx
x1 = -(N1+1)*dx/2 + (1:N1)*dx;
x2 = -(N2+1)*dx/2 + (1:N2)*dx;
[X1,X2] = ndgrid(x1,x2);
% Now array subscripts are equivalent with Cartesian coordinates
% x1 axis points downwards and x2 axis is in horizontal direction
% x1 = X1(:,1) is a column vector in vertical direction
% x2 = X2(1,:) is a row vector in horizontal direction
R = sqrt(X1.^2 + X2.^2);
cgrid = c_s * (R < a) + c_0 * (R >= a);
rhogrid = rho_s * (R < a) + rho_0 * (R >= a);
x1 = X1(:,1); x2 = X2(1,:);
set(figure(8),'Units','centimeters','Position',[5 5 18 12]);
subplot(1,2,1)
imagesc(x2,x1,cgrid);
title('\fontsize{13} c-grid');
xlabel('x_1 \rightarrow');
ylabel('\leftarrow x_3');
axis('equal','tight');
colorbar('hor'); colormap jet;
subplot(1,2,2)
imagesc(x2,x1,rhogrid);
title('\fontsize{13} \rho-grid');
xlabel('x_1 \rightarrow');
ylabel('\leftarrow x_3');
axis('equal','tight');
colorbar('hor'); colormap jet;
%-------------------------------------------------------------------------
% CHECK EXACT INCIDENT FIELD AND BESSEL FUNCTION REPRESENTATION for r = a
%-------------------------------------------------------------------------
% Transform Cartesian coordinates to polar ccordinates
rS = sqrt(xS(1)^2+xS(2)^2); phiS = atan2(xS(2),xS(1));
r = a; phi = 0:.01:2*pi;
% (1) Compute incident wave in closed form --------------------------------
DIS = sqrt(rS^2 + r.^2 - 2*rS*r.*cos(phiS-phi));
p_inc_exact = 1/(2*pi) .* besselk(0,gam_0*DIS);
dp_inc_exact = - gam_0 *(r-rS*cos(phiS-phi))./DIS ...
.* 1/(2*pi) .* besselk(1,gam_0*DIS);
% (2) Compute incident wave as Bessel series with M+1terms --------------
M = 100; % increase M for more accuracy
p_inc = besselk(0,gam_0*rS) .* besseli(0,gam_0*r);
dp_inc = gam_0 * besselk(0,gam_0*rS) .* besseli(1,gam_0*r);
for m = 1 : M;
Ib0 = besseli(m,gam_0*r);
dIb0 = gam_0 * (besseli(m+1,gam_0*r) + m/(gam_0*r) * Ib0);
p_inc = p_inc + 2 * besselk(m,gam_0*rS) * Ib0 .* cos(m*(phiS-phi));
dp_inc = dp_inc + 2 * besselk(m,gam_0*rS) .* dIb0 .* cos(m*(phiS-phi));
end % m_loop
p_inc = 1/(2*pi) * p_inc;
dp_inc = 1/(2*pi) * dp_inc;
% (3) Determine mean error and plot error in domain -----------------------
Error_p = p_inc - p_inc_exact;
disp(['normalized norm of error = ' ...
num2str(norm(Error_p(:),1)/norm(p_inc_exact(:),1))]);
Error_dp = dp_inc - dp_inc_exact;
disp(['normalized norm of error = ' ...
num2str(norm(Error_dp(:),1)/norm(dp_inc_exact(:),1))]);
set(figure(9),'Units','centimeters','Position',[5 5 18 14]);
subplot(2,1,1)
angle = phi * 180 / pi;
semilogy(angle,abs(Error_p)./abs(p_inc_exact)); axis tight;
xlabel('observation angle in degrees \rightarrow');
ylabel('abs(p_{inc}-p_{inc}^{exact}) / abs(p_{inc}^{exact}) \rightarrow');
subplot(2,1,2)
semilogy(angle,abs(Error_dp)./abs(dp_inc_exact)); axis tight;
title('\fontsize{12} relative error on circle boundary');
xlabel('observation angle in degrees \rightarrow');
ylabel('abs(dp_{inc}-dp_{inc}^{exact}) / abs(dp_{inc}^{exact}) \rightarrow');
%--------------------------------------------------------------------------
% COMPUTE SCATTERED FIELD WITH BESSEL FUNCTION REPRESENTATION for r > a
%--------------------------------------------------------------------------
% (4) Compute coefficients of series expansion ----------------------------
Z_0 = c_0 * rho_0; Z_s = c_s * rho_s;
arg0 = gam_0 * a; args = gam_s *a;
A = zeros(1,M+1);
for m = 0 : M;
Ib0 = besseli(m,arg0); dIb0 = besseli(m+1,arg0) + m/arg0 * Ib0;
Ibs = besseli(m,args); dIbs = besseli(m+1,args) + m/args * Ibs;
Kb0 = besselk(m,arg0); dKb0 = -besselk(m+1,arg0) + m/arg0 * Kb0;
A(m+1) = - ((1/Z_s) * dIbs*Ib0 - (1/Z_0) * dIb0*Ibs) ...
/((1/Z_s) * dIbs*Kb0 - (1/Z_0) * dKb0*Ibs);
end
% (5) Compute scattered field at receivers (data) -------------------------
rR = sqrt(xR(1,:).^2 + xR(2,:).^2); phiR = atan2(xR(2,:),xR(1,:));
rS = sqrt(xS(1)^2 + xS(2)^2); phiS = atan2(xS(2),xS(1));
p_sct = A(1) * besselk(0,gam_0*rS).* besselk(0,gam_0*rR);
for m = 1 : M;
factor = 2 * besselk(m,gam_0*rS) .* cos(m*(phiS-phiR));
p_sct = p_sct + A(m+1) * factor .* besselk(m,gam_0*rR);
end % m_loop
p_sct = 1/(2*pi) * p_sct;
P_sct = rho_0 * p_sct; % Take care of factor \rho
end
\ No newline at end of file
fdelmodc/demo/matlab/comparison.m 0 → 100644
+ 83
− 0
View file @ be315f3b
clear all; close all; clc;
display('Running test');
xS = [-100,0]; % source position: 100 m above center
xR = [-60,0]; % central point of receiver array (-50*dxR:50*dxR)
c_0 = 1500; % wave speed in embedding
c_s = 3000; % wave speed in scattering object
rho_0 = 3000; % mass density of enbedding
rho_s = 1500; % mass density of scattering object
dxR = 0.5; % receiver grid size
gsize = [1+240/dxR,1+200/dxR]; % gridsize in [z,x], center of model at coordinate (0,0)
f_in = 50; % selected frequency to compare, returned in Pf
ntrcv = 256; % number of time samples in FD modeling
dtrcv = 0.004; % dt in receivers
% Make grid
a = 40; % radius circle cylinder
N1 = gsize(1); % number of samples in x_1
N2 = gsize(2); % number of samples in x_2
dx = dxR; % with meshsize dx
x1 = -(N1+1)*dx/2 + (1:N1)*dx;
x2 = -(N2+1)*dx/2 + (1:N2)*dx;
[X1,X2] = ndgrid(x1,x2);
% Now array subscripts are equivalent with Cartesian coordinates
% x1 axis points downwards and x2 axis is in horizontal direction
% x1 = X1(:,1) is a column vector in vertical direction
% x2 = X2(1,:) is a row vector in horizontal direction
R = sqrt(X1.^2 + X2.^2);
cgrid = c_s * (R < a) + c_0 * (R >= a);
rhogrid = rho_s * (R < a) + rho_0 * (R >= a);
% DATA from Thorbecke's finite difference code
[Ptsct, Pfinc, Pfsct, f_out]=FD_mod_grid( xS, xR, ntrcv, dtrcv, dxR, cgrid, rhogrid, f_in );
f=f_out % nearest computed discrete frequency
% Compare with analytical solution ---------------------------------------
[P_inc,P_sct,xR] = ForwardCircle( xS, xR, gsize, dxR, c_0, c_s, rho_0, rho_s, f );
set(figure(1),'Units','centimeters','Position',[1 1 30 10]);
subplot(1,3,1)
plot(xR(2,:),real(Pfinc),'LineWidth',1.2);
title('Real part of P^{inc}'); axis tight; hold on;
plot(xR(2,:),real(P_inc),'--r','LineWidth',1.2);
axis tight; hold off;
subplot(1,3,2)
plot(xR(2,:),imag(Pfinc),'LineWidth',1.2);
title('Imaginary part of P^{inc}'); axis tight; hold on;
plot(xR(2,:),imag(P_inc),'--r','LineWidth',1.2);
axis tight; hold off;
subplot(1,3,3)
plot(xR(2,:), abs(Pfinc),'LineWidth',1.2);
title('Absolute value of P^{inc}'); axis tight; hold on;
plot(xR(2,:), abs(P_inc),'--r','LineWidth',1.2);
axis tight; hold off
legendtitle1 = sprintf('f=%.2fHz FiniteDiff', f);
legendtitle2 = sprintf('f=%.2fHz Analytic ', f);
legend(legendtitle1,legendtitle2,'Location','Best');
set(figure(2),'Units','centimeters','Position',[9 12 10 10]);
error = abs(P_sct(:)-Pfsct(:))./abs(Pfsct(:));
plot(error,'LineWidth',1.2); title('Relative error'); axis tight;
set(figure(3),'Units','centimeters','Position',[1 1 30 10]);
subplot(1,3,1)
plot(xR(2,:),real(Pfsct),'LineWidth',1.2);
title('Real part of P_{sct}'); axis tight; hold on;
plot(xR(2,:),real(P_sct),'LineWidth',1.2);
subplot(1,3,2)
plot(xR(2,:),imag(Pfsct),'LineWidth',1.2);
title('Imaginary part of P^{sct}'); axis tight; hold on;
plot(xR(2,:),imag(P_sct),'LineWidth',1.2);
subplot(1,3,3)
plot(xR(2,:), abs(Pfsct),'LineWidth',1.2);
title('Absolute value of P^{sct}'); axis tight; hold on;
plot(xR(2,:), abs(P_sct),'LineWidth',1.2);
axis tight; hold off
legendtitle1 = sprintf('f=%.2fHz FiniteDiff', f);
legendtitle2 = sprintf('f=%.2fHz Analytic ', f);
legend(legendtitle1,legendtitle2,'Location','Best');
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment