It appears to be fortuitous... When I use the velocity as source term in 2D, it does not match...
I think the following observations are reliable (in 1D and 2D):
- using only p as source term in the second simulation results in too low power.
- using only u results in higher power (too much in 2D, the good one in 1D, probably by chance)
- using both results in even higher power...
Note : I always use Dirichlet source terms.
However, I think I have a workaround: record the signal and set it as source term not only in one plane but in several consecutive planes (with pressure for example, I don't think it matters...).
--> It works nicely for around 10 planes (source_width = 7 for example). I think it has a more "coercive effect" on the field in the second domain. I don't know if it makes sense considering the spatial colocation method, but it seems to work. I will try it in 3D and keep you informed.
Best regards,
Anthony
clear all
close all
clc
% Parameters
CFL = 0.3;
velocity_source = false;
pressure_source = true;
source_width = 7;
% create the computational grid
PMLsize = 20;
Nx = 88;
Ny = 88;
dx = 0.1e-3;
dy = 0.1e-3;
kgrid = makeGrid(Nx, dx, Ny, dy);
% define the properties of the propagation medium
medium.sound_speed = 1500;
medium.density = 1000;
% create time array
t_end = 1.1*sqrt(kgrid.x_size.^2 + kgrid.y_size.^2)/min(medium.sound_speed(:));
pixel_dim = min([kgrid.dx, kgrid.dy]);
dt = CFL*pixel_dim/max(medium.sound_speed(:));
source_freq = 3e6;
periode = 1/source_freq;
NtimeStepsPerPeriod = ceil(periode/dt);
dt = periode/NtimeStepsPerPeriod;
kgrid.t_array = 0:dt:t_end;
% create initial pressure distribution using makeDisc
R = Ny*0.45;
Rap = 28/38*R;
source.p_mask = zeros(Nx, Ny);
xc = round(Nx/2);
yc = round(Ny/2);
xgrid = (1:Nx)-xc;
ygrid = ((1:Ny)-yc)';
rgrid = sqrt(bsxfun(@plus,xgrid.^2,ygrid.^2));
source.p_mask(logical(bsxfun(@times,(abs(rgrid-R)<1/2), bsxfun(@times,ygrid<0,abs(xgrid)<Rap)))) = 1;
source.p_mode = 'dirichlet';
% define source
source_mag = 0.1e6;
source_cycles = 300;
source_padding = 10;
source.p = source_mag*sin(2*pi*kgrid.t_array*source_freq);
% source.p = [source_mag * toneBurst(1/kgrid.dt, source_freq, source_cycles, 'Envelope', 'Rectangular'), zeros(1, source_padding)];
% define sensor
sensor.mask = ones(Nx, Ny);
sensor.record = {'p', 'u_non_staggered','u'};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,...
'DisplayMask', 'off', 'PlotScale', [-0.5, 0.5]*source_mag,'PMLInside',false,'PlotPML',false,'PMLSize',[PMLsize PMLsize]);
% Get the simulated pressure field over a plane
planeIdx = 25; % arbitrary
p_input = reshape(sensor_data.p,Nx,Ny,size(sensor_data.p,2));
p_input = p_input(planeIdx:planeIdx+source_width-1,:,:);
p_input = reshape(p_input,[source_width*Ny,size(p_input,3)]);
% Simulate the rest of the domain starting from the "intercepted plane"
sourceShift = 10;
Nx_new = Nx - planeIdx + 1 + sourceShift;
kgrid2 = makeGrid(Nx_new, dx, Ny, dy);
kgrid2.t_array = kgrid.t_array;
% Set input
if pressure_source
source2.p_mask = zeros(Nx_new, Ny);
source2.p_mask(1+sourceShift:1+sourceShift+source_width-1,:) = true;
source2.p = p_input;
source2.p_mode = 'dirichlet';
end
if velocity_source
ux_input = reshape(sensor_data.ux,Nx,Ny,size(sensor_data.ux,2));
ux_input = ux_input(planeIdx:planeIdx+source_width-1,:,:);
ux_input = reshape(ux_input,[source_width*Ny,size(ux_input,3)]);
uy_input = reshape(sensor_data.uy,Nx,Ny,size(sensor_data.uy,2));
uy_input = uy_input(planeIdx:planeIdx+source_width-1,:,:);
uy_input = reshape(uy_input,[source_width*Ny,size(uy_input,3)]);
if pressure_source
source2.u_mask = source2.p_mask;
else
source2.u_mask = zeros(Nx_new, Ny);
source2.u_mask(1+sourceShift:1+sourceShift+source_width-1,:) = true;
end
source2.ux = ux_input;
source2.uy = uy_input;
source2.u_mode = 'dirichlet';
end
% set sensor2
sensor2.mask = ones(Nx_new, Ny);
sensor2.record = {'p', 'u_non_staggered'};
% run the simulation
sensor_data2 = kspaceFirstOrder2D(kgrid2, medium, source2, sensor2,...
'DisplayMask', 'off', 'PlotScale', [-0.5, 0.5]*source_mag,'PMLInside',false,'PlotPML',false,'PMLSize',[PMLsize PMLsize]);
% time shift
ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt);
uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt);
ux2 = timeShift(sensor_data2.ux_non_staggered, kgrid2.dt);
uy2 = timeShift(sensor_data2.uy_non_staggered, kgrid2.dt);
% Intensity
NpointsToKeep = NtimeStepsPerPeriod;
Ix = reshape(sensor_data.p.*ux,[Nx,Ny,size(ux,2)]);
Iy = reshape(sensor_data.p.*uy,[Nx,Ny,size(uy,2)]);
Ix = mean(Ix(:,:,end-NpointsToKeep+1:end),3);
Iy = mean(Iy(:,:,end-NpointsToKeep+1:end),3);
Ix2 = reshape(sensor_data2.p.*ux2,[Nx_new,Ny,size(ux2,2)]);
Iy2 = reshape(sensor_data2.p.*uy2,[Nx_new,Ny,size(uy2,2)]);
Ix2 = mean(Ix2(:,:,end-NpointsToKeep+1:end),3);
Iy2 = mean(Iy2(:,:,end-NpointsToKeep+1:end),3);
%% Display
maskIdx = find(source.p_mask);
[maskI,~] = ind2sub(size(source.p_mask),maskIdx);
clear maskIdx
endTrdI = max(maskI);
figure(2)
clf
sidePowerLoss = ( cumsum(-Iy(:,1),1)+cumsum(Iy(:,end),1) ) * dx;
compensatedPower = squeeze(sum(Ix,2)) * dy + sidePowerLoss;
sidePowerLoss2 = ( cumsum(-Iy2(:,1),1) + cumsum(Iy2(:,end),1) ) * dx;
compensatedPower2 = squeeze(sum(Ix2,2)) * dy + sidePowerLoss2 + sidePowerLoss(planeIdx);
plot(compensatedPower, '-', 'LineWidth',2)
hold on
plot(planeIdx-sourceShift:planeIdx-sourceShift+numel(compensatedPower2)-1,compensatedPower2, '-', 'LineWidth',2)
ylim([0 1.1*max([compensatedPower;compensatedPower2])])
ylabel('Integrated intensity [W]')
xlabel('Main propagation axis [pixels]')
set(gca,'FontSize',18)
title(['Blue and red should be superimposed'])
){kind=link}