Hi @peyton,
To increase the speed of the for loop, you can use vectorization technique.
Here are some of the key steps you can take to vectorize the code:
- 2D Laplacian with conv2:Instead of manually computing the second derivative with multiple matrix slices, I used conv2, which is a built-in MATLAB function to perform a 2D convolution.
- Vectorized Update for the Concentration Field:The concentration update can be done for the entire grid in one step as shown in the code below.
- Boundary Condition Handling: Boundary conditions are applied directly using matrix indexing without looping over individual grid points. After each time step, we simply reset the boundary values for the grid.
Here is the updated code using vectorization which helps to increase the performance of the model:
tmax = 60*30;
c0 = 0;
D0 = 2.4e-3;
Q = 3.87;
T = 1173;
kb = 8.6173303e-5;
xmax = 1000;
ymax = 300;
h = 1;
ct = 650;
cb = 0;
D = D0 * exp(-Q / (kb * T)) * (1e9)^2;
dt = h^2 / (2 * D) / 2;
nsteps = ceil(tmax / dt);
x = (0:h:xmax);
y = (0:h:ymax);
x = [-h, x, xmax + h];
[X, Y] = meshgrid(x, y);
c = zeros(size(X));
c(:) = c0;
x1 = xmax * (0.5 - open_width / 2);
x2 = xmax * (0.5 + open_width / 2);
top_vls = (x < x2) & (x > x1) * ct;
c(1,:) = top_vls;
c(end,:) = cb;
c(:,1) = c(:,2);
c(:,end) = c(:,end-1);
% Precompute matrices for the 2D Laplacian operator (second derivative)
laplacian_filter = [1, -4, 1]; % 1D Laplacian kernel for finite difference
for i = 1:nsteps
if i == nsteps
dt = tmax - t; % Adjust last time step
end
t = t + dt; % Update the current time
% Calculate the second spatial derivative using the 2D Laplacian
D2c = conv2(c, laplacian_filter, 'same') / h^2;
% Update concentration for the interior points
c(2:end-1, 2:end-1) = c(2:end-1, 2:end-1) + dt * D * D2c(2:end-1, 2:end-1);
c(1,:) = top_vls;
c(end,:) = cb;
c(:,1) = c(:,2);
c(:,end) = c(:,end-1);
end
You can learn more about the conv2 function here:
