Hi, Thank you very much for this code, its really helpfull!
I have concerns about using it in a specific calculation. The idea is to simulate the propagation of ultrasound at very low frequency (8kHz) inside a tube filled with water. There are two opposite sources on the side of the tube. Everything is supposed to be symetric here but the two sources behave very differently. How do you exlain the asymmetry in the calculation? Is there something wrong in my code?
'Nx = 101; % number of grid points in the x (row) direction
Ny = 101; % number of grid points in the y (column) direction
dx = 0.0005; % grid point spacing in the x direction [m]
dy = 0.0005; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);
% define the properties of the propagation medium
ray_int=18e-3; % [m]
medium.sound_speed_compression = 330*ones(Nx, Ny); % [m/s]
medium.sound_speed_compression(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx+3)==1) = 5650;
medium.sound_speed_compression(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx)==1) = 1480;
medium.sound_speed_shear = 0*ones(Nx, Ny); % [m/s]
medium.sound_speed_shear(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx+3)==1) = 3100;
medium.sound_speed_shear(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx)==1) = 0;
medium.density = 2*ones(Nx, Ny); % [kg/m^3]
medium.density(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx+3)==1) = 8010;
medium.density(makeDisc(Nx,Ny,(Nx+1)/2,(Ny+1)/2,ray_int/dx)==1) = 1000;
% define a curved transducer element
source.u_mask = zeros(Nx,Ny);
source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2-ray_int/dx-2)=1;% create the time array
source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2+ray_int/dx+2)=1;% create the time array
% source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2-ray_int/dx-2)=1;% create the time array
% source.u_mask((Nx+1)/2-20:(Nx+1)/2+20,(Nx+1)/2+ray_int/dx+1)=1;% create the time array
% define a time varying sinusoidal source
source_freq1 = 8e3; % [Hz]
source_freq2 = 8e3; % [Hz]
source_mag = 0.03
t_end=1/source_freq2;
% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, 16000, [], t_end);
signal_uy_plus = source_mag*sin(2*pi*source_freq1*kgrid.t_array);
signal_uy_moins = -source_mag*sin(2*pi*source_freq2*kgrid.t_array);
% filter the source to remove any high frequencies not supported by the grid
source.uy = [repmat(filterTimeSeries(kgrid, medium, signal_uy_plus),41,1);repmat(filterTimeSeries(kgrid, medium, signal_uy_moins),41,1)];
display_mask = source.u_mask;
% define a centered circular sensor
% create a sensor mask covering the entire computational domain using the
% opposing corners of a rectangle
Pos_x=(Nx+1)/2;
Pos_y=(Ny+1)/2;
sensor.mask = zeros(Nx,Ny);
sensor.mask (Pos_x, :)=1;
sensor.mask (:,Pos_y)=1;
sensor.record={'p'};
input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false,'PlotScale', [-100e3,100e3,-100e3,100e3],};
sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});'
I supposed the definition of the sources is not correct. By changing the location of the sources (commented lines), the results are slightly more symetric..
Subsidiary question: Would it be better with a coupling between sources and tube?
Thank you in advance.
Damien