I'm not certains what this code is supposed to do or the result it's supposed to produce.
The problems were that the diff calls were shortening the matrices they operated on by one row or column, and that was throwing the dimension errors. This 'sort of' fixes that by appending a row or column of zeros to the respective matrices.
A better solution could be to use the gradient function instead of diff, since the results of those operations are the same sizes as the arguments. I leave that experiment to you. Note that if you use gradient, use the original statements with it (that I have left in and commented out), not the ones I latered to make the diff calls work.
Try this --
% solve_2d_wave_pml.m
% This code simulates the 2D acoustic wave equation with a Perfectly Matched Layer (PML)
% using the Finite-Difference Time-Domain (FDTD) method.
%
% `200`: The number of grid points in the x-direction (`nx`).
% `200`: The number of grid points in the y-direction (`ny`).
% `500`: The number of time steps for the simulation (`nt`).
% `343`: The speed of sound in the medium, in m/s (`c0`). This is approximately the speed of sound in air.
% `20`: The thickness of the Perfectly Matched Layer in grid points (`dpml`).
% `100`: The x-coordinate of the wave source (`source_x`).
% `100`: The y-coordinate of the wave source (`source_y`).
clear
clc
%Describe the experimental parameters
nx = 200;
ny = 200;
nt = 500;
c0 = 343;
dpml = 20;
source_x = 100;
source_y = 100;
% Physical parameters
dx = 1.0; % Spatial grid step (meters)
dy = 1.0; % Spatial grid step (meters)
dt = 0.5 * dx / c0; % Time step (stability condition: dt <= 1 / (c0 * sqrt(1/dx^2 + 1/dy^2)))
% Wave field variables
p = zeros(nx, ny); % Total pressure field
vx = zeros(nx, ny); % Particle velocity in x-direction
vy = zeros(nx, ny); % Particle velocity in y-direction
% PML parameters
m = 3; % Power for the damping profile (higher value -> faster ramp-up)
R0 = 1e-6; % Reflection coefficient
sigmax = zeros(nx, 1); % Damping profile for x-direction
sigmay = zeros(ny, 1); % Damping profile for y-direction
% The split-field PML requires auxiliary variables
p_x = zeros(nx, ny); % Split pressure field for x-direction
p_y = zeros(nx, ny); % Split pressure field for y-direction
% Split velocity fields
vx_x = zeros(nx, ny); % Split velocity in x-dir, x-part
vx_y = zeros(nx, ny); % Split velocity in x-dir, y-part
vy_x = zeros(nx, ny); % Split velocity in y-dir, x-part
vy_y = zeros(nx, ny); % Split velocity in y-dir, y-part
% PML damping profile calculation
x_range = 0:dx:dx*(dpml-1);
sigma_profile = ((x_range / (dpml*dx)).^m) * ((m+1) * c0 * log(1/R0) / (2 * (dpml*dx)));
% Apply PML damping profiles to the grid
sigmax(1:dpml) = fliplr(sigma_profile);
sigmax(nx-dpml+1:nx) = sigma_profile;
sigmay(1:dpml) = fliplr(sigma_profile);
sigmay(ny-dpml+1:ny) = sigma_profile;
% Create source signal (e.g., a simple Gaussian pulse)
f_center = 10; % Center frequency of the source (Hz)
t_source = (0:nt-1) * dt;
source_signal = exp(-((t_source - 1/f_center) / (0.1 * 1/f_center)).^2) * 1e-4;
% Visualization setup
figure;
h = imagesc(p, [-1e-4, 1e-4]);
axis equal tight;
title('2D Wave Propagation with PML');
xlabel('x');
ylabel('y');
caxis([-1e-4 1e-4]);
colormap(jet);
colorbar;
% Main FDTD time-stepping loop
for t = 1:nt
% Update velocity fields (x-direction)
% This is the core FDTD update equation, split for the PML
% q1 = vx_x
% q2 = dt * sigmax(1:nx)'
% q3 = dt/dx
% q4 = diff([p_x(:,1:end-1), zeros(nx,1)], 1, 2)
% vx_x = vx_x - dt * sigmax(1:nx)' .* vx_x - dt/dx * diff([p_x(:,1:end-1), zeros(nx,1)], 1, 2);
% vx_y = vx_y - dt * sigmay(1:ny) .* vx_y - dt/dy * diff([p_y(1:end-1,:); zeros(1,ny)], 1, 1);
vx_x = vx_x - dt * sigmax(1:nx)' .* vx_x - dt/dx * diff([p_x(:,1:end), zeros(nx,1)], 1, 2);
vx_y = vx_y - dt * sigmay(1:ny) .* vx_y - dt/dy * diff([p_y(1:end,:); zeros(1,ny)], 1, 1);
vx = vx_x + vx_y;
% Update velocity fields (y-direction)
% vy_x = vy_x - dt * sigmax(1:nx)' .* vy_x - dt/dx * diff([p_x(:,1:end-1), zeros(nx,1)], 1, 2);
% vy_y = vy_y - dt * sigmay(1:ny) .* vy_y - dt/dy * diff([p_y(1:end-1,:); zeros(1,ny)], 1, 1);
vy_x = vy_x - dt * sigmax(1:nx)' .* vy_x - dt/dx * diff([p_x(:,1:end), zeros(nx,1)], 1, 2);
vy_y = vy_y - dt * sigmay(1:ny) .* vy_y - dt/dy * diff([p_y(1:end,:); zeros(1,ny)], 1, 1);
vy = vy_x + vy_y;
% Update pressure fields (split)
% This is the core FDTD update equation, split for the PML
% q5 = p_x
% q6 = dt * sigmax(1:nx)'
% q7 = c0^2
% vx
% q8 = diff(vx, 1, 2)
% p_x = p_x - dt * sigmax(1:nx)' .* p_x + dt * c0^2 * (diff(vx, 1, 2) / dx);
% p_y = p_y - dt * sigmay(1:ny) .* p_y + dt * c0^2 * (diff(vy, 1, 1) / dy);
p_x = p_x - dt * sigmax(1:nx)' .* p_x + dt * c0^2 * (diff([vx zeros(size(vx,1),1)], 1, 2) / dx);
p_y = p_y - dt * sigmay(1:ny) .* p_y + dt * c0^2 * (diff([vy; zeros(1,size(vy,2))], 1, 1) / dy);
% Add source at specified location
p_x(source_y, source_x) = p_x(source_y, source_x) + source_signal(t);
% Combine split pressure fields for plotting
p = p_x + p_y;
% Visualize the wave field
set(h, 'CData', p);
drawnow;
end
.
