2D wave equation- tsunami simulation

I need to solve the 2D wave equation, which is written as follows:

Using Matlab, I have to simulate a 2D tsunami generated by an earthquake, given Neumann condition, a given function for lambda and a gaussian as initial displacement. The lambda function provided by the exercise is the one I've used in the code.

The code I've written is (I'm just learning how to program):

% Define simulation parameters
Lx = 1000; % domain size in x-direction (m)
Ly = 1000; % domain size in y-direction (m)
Nx = 100; % number of grid points in x-direction
Ny = 100; % number of grid points in y-direction
dx = Lx/Nx; % grid spacing in x-direction
dy = Ly/Ny; % grid spacing in y-direction
g = 9.81; % (m/s^2)
H0 = 5; % mean water depth (m)
c = sqrt(g*H0); % wave speed
theta = 1; % inclination angle (degree)
T = 150; % total simulation time (s)
dt = 0.150; % time step

% Define the variable diffusion coefficient
lambda = zeros(Ny, Nx); % variable diffusion coefficient
X_0 = Lx/8; % point where the sea bed changes
G = g*H0; % maximum value of lambda
for i = 1:Nx
    for j = 1:Ny
        if (i-1)*dx < X_0
            lambda(j,i) = G;
        else
            lambda(j,i) = g*(H0 - ((i-1)*dx-X_0)*0.0174550649282176);
        end
    end
end

% Initialize the solution
u = zeros(Ny, Nx); % current solution
u_old = zeros(Ny, Nx); % previous solution
u_xx = zeros(Ny, Nx); % second derivative in x-direction
u_yy = zeros(Ny, Nx); % second derivative in y-direction
D_xx = zeros(Ny, Nx); % diffusion term in x-direction
D_yy = zeros(Ny, Nx); % diffusion term in y-direction

% Define the initial conditions
sigmax = 10; % width of the Gaussian function (m)
sigmay = 30; % width of the Gaussian function (m)
x0 = Lx/4 ; % x-coordinate of the center of the Gaussian function
y0 = Ly/2; % y-coordinate of the center of the Gaussian function
u0 = 2*exp(-(dx*(0:Nx-1)-x0).^2/(sigmax^2))'*exp(-(dy*(0:Ny-1)-y0).^2/(sigmay^2)); % Gaussian function
u = u0; % set the initial condition


% Visualization
[X, Y] = meshgrid(dx*(0:Nx-1), dy*(0:Ny-1)); % create grid for plotting
surf(X, Y, u); % plot the initial condition
xlabel('x');
ylabel('y');
zlabel('Surface Elevation');
title(sprintf('Time = %.3f', 0));
axis([0 Lx 0 Ly -max(max(u0)) max(max(u0))]);
view(3);
colorbar;
drawnow;

% Initialize maximum wave elevation at coastline
max_elevation = max(u(end, :));

% Time-stepping loop
for n = 1:(T/dt)
    % Calculate the second derivatives using finite differences
    u_xx = (circshift(u, [0 -1]) - 2*u + circshift(u, [0 1]))/(dx^2);
    u_yy = (circshift(u, [-1 0]) - 2*u + circshift(u, [1 0]))/(dy^2);

    % Calculate the diffusion terms using the variable diffusion coefficient
    D_xx = lambda.*u_xx;
    D_yy = lambda.*u_yy;
     % Calculate the CFL condition
    max_lambda = max(max(lambda));
    CFL = dt <= 1/(c*sqrt((1/dx^2) + (1/dy^2)));
    if CFL > 1
        dt = 0.9/(c*sqrt((1/dx^2) + (1/dy^2)));
    end
    
    % Calculate the boundary terms using the Neumann boundary condition
    u(:, 1) = u(:, 2);
    u(:, end) = u(:, end-1);
    u(1, :) = u(2, :);
    u(end, :) = u(end-1, :);
    
    % Update the solution using the wave equation with variable diffusion coefficient
    u_new = 2*u - u_old + c^2*dt^2*(u_xx + u_yy) + dt^2*D_xx + dt^2*D_yy;
    
    % Update the solution for the next time step
    u_old = u; % store the current solution
    u = u_new; % update the solution
    
    % Update maximum wave elevation at coastline
    max_elevation = max(max_elevation, max(u(end, :)));
    
    % Visualization
    surf(X, Y, u); % plot the surface elevation
    xlabel('x');
    ylabel('y');
    zlabel('Surface Elevation');
    title(sprintf('Time = %.3f', n*dt));
    axis([0 Lx 0 Ly -50 50]);
    view(3);
    colorbar;
    drawnow;
end

% Print the maximum wave elevation at the coastline
fprintf('Maximum wave elevation at coastline: %.3f\n', max_elevation);

And it seems to work fine. The problem is that the exercise asks me to consider sigmax=80000 m, sigmay=200000 m and H0= 5000m . When I change the code with these parameters, the algorithm is no more stable. I've tried to change the number of grid points and time steps too, but doesn't work.

If someone could help me, I would really appreciate it. Thank you.

[Florian Bidaud]

Hi

Are you also increasing the domain size when you change sigmax, sigmay and H0 ?

[Torsten]

"lambda" depends on x and y. So you can't take it out of the divergence operator:

D_xx = lambda.*u_xx;
D_yy = lambda.*u_yy;

[Davide]

Thank you for the answers, I really appreciete it!

Actually, I have considered an incorrect discretization.

The exercise suggests this discretization for the second member of the wave equation:

The lambda functions is defined costant before and after it just depends on the variable "x".

So I change the code in this way:

% Calculate the second derivatives using finite differences
    u_xx = (circshift(u, [0 -1]) - 2*u + circshift(u, [0 1]))/(dx^2);
    u_yy = (circshift(u, [-1 0]) - 2*u + circshift(u, [1 0]))/(dy^2);

    % Calculate the diffusion terms using the variable diffusion coefficient
    D_xx = (1/dx) * (lambda(:, [2:Nx, Nx]) .* ((u(:, [2:Nx, Nx]) - u) ./ dx) - lambda(:, [1, 1:Nx-1]) .* ((u - u(:, [1, 1:Nx-1])) ./ dx));
    D_yy = (1/dy^2) * (lambda .* (circshift(u, [-1, 0]) - 2*u + circshift(u, [1, 0])));

But the problem still persists when I try to enter the values of the parameters required by the exercise I showed in the question. Clearly, I have also considered larger domains, but that has not been enough.

[Torsten]

You can't change dt within the loop because the loop index depends on dt:

for n = 1:(T/dt)
    ...
  if CFL > 1
    dt = 0.9/(c*sqrt((1/dx^2) + (1/dy^2)));
  end
   ...
end

Use a while-loop instead:

t = 0;
while t < T
  ...
  t = t + dt;
end

[Davide]

I'm sorry if I continue to bother you, but I have also changed this following your suggestion, but it still doesn't work.

[Torsten]

Did you try to use the PDE toolbox for a reference solution ?

Are you sure that an initial condition for u is sufficient ? Since you solve a 2nd order PDE in time, I would have expected that you also have to prescribe du/dt.

[Davide]

You are right! I had not considered the condition concerning du/dt... As I mentioned, I have recently started using MATLAB, so I'm constantly searching for things online through various forums. After some research, I rewrote the code in the following way, but it only plots the initial condition and not its evolution, as the previous code did.

% Parameters
Lx = 1000000; % Length of the square region in x-direction (m)
Ly = 1000000; % Length of the square region in y-direction (m)
Nx = 1000; % Number of grid points in x-direction
Ny = 1000; % Number of grid points in y-direction
dx = Lx / Nx; % Grid spacing in x-direction
dy = Ly / Ny; % Grid spacing in y-direction
G = 9.81; % Acceleration due to gravity (m/s^2)
H0 = 5000; % Still water depth (m)
X0 = 30; % Point where the sea bed changes (m)
A0 = 1; % Elevation of the surface (m)
sigmax = 80000; % Extension of the surface elevation in x-direction (m)
sigmay = 200000; % Extension of the surface elevation in y-direction (m)
tFinal = 300; % Final simulation time (s)
dt = 0.1; % Time step (s)

% Stability condition
lambda_max = G * H0;
CFL = dt / sqrt(lambda_max) / sqrt(1/dx^2 + 1/dy^2);

if CFL > 1
    error('Stability condition is not satisfied! Reduce time step (dt) or increase grid resolution.');
end

% Initialize the grid
x = linspace(0, Lx, Nx+1);
y = linspace(0, Ly, Ny+1);
[X, Y] = meshgrid(x, y);

% Initialize the solution matrix (including ghost cells)
u = zeros(Nx+3, Ny+3, 2);

% Apply the initial conditions
u(2:Nx+2, 2:Ny+2, 1) = A0 * exp(-((X-Lx/2)/sigmax).^2 - ((Y-Ly/2)/sigmay).^2);
u(:,:,2) = u(:,:,1);

% Plot the initial condition
figure;
surf(X, Y, u(2:Nx+2, 2:Ny+2, 1), 'EdgeColor', 'none');
colormap jet;
colorbar;
title('Initial Condition');
xlabel('x (m)');
ylabel('y (m)');
zlabel('Surface Elevation (m)');
axis([0 Lx 0 Ly -A0 A0]);
view(3);
drawnow;

% Time-stepping loop
numSteps = round(tFinal / dt);
for step = 1:numSteps
    for i = 2:Nx+2
        for j = 2:Ny+2
            lambda_x1 = G * H0;
            lambda_x2 = G * (H0 - max(0, (x(i-1) - X0)) * tand(1));
            lambda_y1 = G * H0;
            lambda_y2 = G * (H0 - max(0, (y(j-1) - X0)) * tand(1));
            
            uxx = (lambda_x1 * (u(i+1,j,1) - u(i,j,1)) - lambda_x2 * (u(i,j,1) - u(i-1,j,1))) / dx^2;
            uyy = (lambda_y1 * (u(i,j+1,1) - u(i,j,1)) - lambda_y2 * (u(i,j,1) - u(i,j-1,1))) / dy^2;
            u(i,j,2) = 2 * u(i,j,1) - u(i,j,2) + dt^2 * (uxx + uyy);
        end
    end
    
    % Apply the initial condition du/dt = 0
    u(:,:,2) = u(:,:,1);
    
    % Update u for the next time step
    u(:,:,1) = u(:,:,2);
    
    % Apply the Neumann boundary conditions
    u(1,:,:) = u(2,:,:); % x = 0
    u(Nx+3,:,:) = u(Nx+2,:,:); % x = Lx
    u(:,1,:) = u(:,2,:); % y = 0
    u(:,Ny+3,:) = u(:,Ny+2,:); % y = Ly

    % Plot the current solution
    if mod(step, 10) == 0
        figure;
        surf(X, Y, u(2:Nx+2, 2:Ny+2, 1), 'EdgeColor', 'none');
        colormap jet;
        colorbar;
        title(sprintf('Time = %.1f s', step * dt));
        xlabel('x (m)');
        ylabel('y (m)');
        zlabel('Surface Elevation (m)');
        axis([0 Lx 0 Ly -A0 A0]);
        view(3);
        drawnow;
    end
end

% Find the maximum wave elevation that reaches the coastline
max_elevation = max(max(u(2:Nx+2, 2:Ny+2, 1)));

fprintf('Maximum wave elevation that reaches the coastline: %.2f m\n', max_elevation);

Answer 1

[Yatharth]

Hi Davide, I understand that you are able to solve the initial problem of using larger parameters in the comments itself, however you are not getting the desired output of plotting the evolution in a single graph.

The reason you are getting individual plots for each time step is that you have included “figure” inside your "for loop" where you are plotting the current solution.

    % Plot the current solution
    if mod(step, 10) == 0
        % figure; % comment this figure and include it outside the for loop
        % where you have defined numSteps
        
        surf(X, Y, u(2:Nx+2, 2:Ny+2, 1), 'EdgeColor', 'none');
        colormap jet;
        colorbar;
        title(sprintf('Time = %.1f s', step * dt));

I hope this helps.