2D acoustic wave equation with a Perfectly Matched Layer

Hi,

I got this code generated from Gemini AI:

% 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

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 = 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 = vy_x + vy_y;

% Update pressure fields (split)

% This is the core FDTD update equation, split for the PML

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);

% 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

And it came up with this error:

Arrays have incompatible sizes for this operation.

Error in solve_2d_wave_pml (line 80)

vx_x = vx_x - dt * sigmax(1:nx)' .* vx_x - dt/dx * diff([p_x(:,1:end-1), zeros(nx,1)], 1, 2);

^^^^

Can someone explain what is going on there. I would like to get my head around this.

If I can see the end result that would be brilliant

Thanks

Answer 1

[Star Strider]

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

●

.

[Andrew]

It is not what I expected. I was looking for a source at the centre that radiated out with no reflections.

This AI is very interesting though.

Thanks

[Star Strider]

As always, my pleasure!

I'm not certain what it's supposed to do. Feel free to experiment with it, now that it no lionger throws the dimension errors.

Also, in addition to using imagesc, consider using surf or surfc to draw the complete matrix. That might reveal a bit more detail. And of course, experiment with gradient, since that will likely produce a better result that diff, expecially since gradient uses a central difference approximation for the spatial derivative, as well as producing a result that's the same size as the argument (vector or matrix).