Equation FKPP difference finie 2D inplicite Crank-Nicolson avec approximation ADI: alternate direction implicit method

FKPP_2D_fd_implicite_ADI.m

% Equation FKPP difference finie 2D inplicite Crank-Nicolson avec
% approximation ADI approximation: alternate direction implicit method
% FKPP equation finite-difference explicite 2D

% parametre des equations/equation parameters
r = 0.5; % taux de croissance
D = 0.1; % coefficient de diffusion

% parametres de simulation, espace/space simulation parameters
% x in Omega [0,Sx]x[0,Sy]
h = 0.025;  % intervalle de discretisation spatiale/space step
Sx = 10.0;
Sy = 6.0;
x0 = 0;   % bord gauche de l'intervalle/left end of the interval
x1 = Sx;  % bord droit de l'intervalle/right end of the interval
y0 = 0;
y1 = Sy;
x = x0:h:x1; % discretisation/discretization
y = y0:h:y1; 
[X,Y] = meshgrid(x,y);
J1 = length(x); 
J2 = length(y);
J  = J1*J2;  % nombre de points de discretisation/number of points

% variable dynamiques/dynamical variables
u = zeros(J,1); % stocke seulement l'etat au temps t/store only current sol
u1 = zeros(J,1);

% discretisation du Laplacien avec conditions aux bord de Neumann
% no flux, i.e. du/dx = 0 en x0 et en x1
% Neumann, no flux boundary conditions

% Le Laplacien discretise est une matrice L de taille JxJ
% Discretized Laplacian

cornertopleft = 1;
cornerbottomleft = J2;
cornertopright = (J1-1)*J2+1;
cornerbottomright = J;
sideleft = 2:J2-1;
sidetop = J2+1:J2:J2*(J1-2)+1;
sidebottom = 2*J2:J2:J2*(J1-1);
sideright = J2*(J1-1)+2:J-1;
side = [cornertopleft, cornertopright, cornerbottomleft, cornerbottomright, ...
    sideleft, sidetop, sidebottom, sideright];
interiorY = setdiff(1:J, side);

[ix,iy] = ind2sub([J2,J1],interiorY);
interiorX = sub2ind([J1,J2],iy,ix);

%% test indices Y
% MY = zeros(J2,J1);
% MY(interiorY) = 1;
% spy(MY)
%% test indices X
% MX = zeros(J1,J2);
% MX(interiorX) = 1;
% spy(MX)

%% interior
% On construit deux matrices du Laplacien discretise: LX et LY. LX est le Laplacien
% discretise en prenant les indices en X d'abord:
% We build two matrices for the discretized Laplacian: LX and LY. LX
% is the Laplacian taking indices along the x-axis first.
% 
%   1    2 3 4 5 ...   J1
%   J1+1 ...         2*J1
%   ...
%   J1*(J2-1)+1 ... J1*J2
%
% LY est le Laplacien discretise en prenant les indices le long de Y d'abord:
% LY is the disctretized Laplacian taking indices along the y-axis first 
%  
%   1 J2+1 ...(J1-1)*J2+1
%   2 J2+2            ...
%   3 
%   ...               ...
%   J2 2*J2 3*J ... J1*J2
%
% LX et LY sont des matrices tridiagonales
% LX and LY are tri-diagonal matrices 

LX = sparse(interiorX,interiorX,-2,J,J); % matrice creuse, compacte en memoire
LX = LX + sparse(interiorX,interiorX+1,1,J,J);
LX = LX + sparse(interiorX,interiorX-1,1,J,J);

LY = sparse(interiorY,interiorY,-2,J,J); % matrice creuse, compacte en memoire
LY = LY + sparse(interiorY,interiorY+1,1,J,J);
LY = LY + sparse(interiorY,interiorY-1,1,J,J);

% condition initiales/initial condtions
% u( X(:).^2 + Y(:).^2 < 3 ) = 0.8; % condition initiale constante par morceaux
                                    % piecewise constant initial condition
u = exp( - X(:).^2/0.5 - Y(:).^2/2 ); % condition initiale lisse/smooth IC


% parametres de simulation, temps/time simulation parameters
t0 = 0;
tfinal = 20; 
t = t0;
tp = t0;
k = .1;

% Schema implicite Crank-Nicolson implicit scheme
% ADI: alternate direction implicit method
AX = (speye(J) - k/h^2*D/2*LX);
AY = (speye(J) - k/h^2*D/2*LY);

% figures
figure(2); clf;
image(x,y,255*reshape(u,J2,J1));

% uncomment the next two lines to pause the simulation at the initial condition
% disp('appuyer sur une touche pour continuer');
% pause

% BOUCLE PRINCIPALE/MAIN LOOP
while t < tfinal
    drawnow;
    % Methode ADI:
    % On resout en 2 etapes: d'abord on fait par rapport a X:
    % u12 = u + k/2*r*u*(1-u) + k/2/h^2*D*(LX*u12 + LY*u)
    % ensuite par rapport a Y:
    % u = u12 + k/2*r*u12*(1-u12) + k/2/h^2*D*(LX*u12 + LY*u)
    % L'evaluation se fait a k/2 a chaque etape

    % ADI method: 
    % Two step solver: first we apply diffusion only along the x-axis
    % u12 = u + k/2*r*u*(1-u) + k/2/h^2*D*(LX*u12 + LY*u)
    % then we make a change of index and apply diffusion along the y-axis
    % u = u12 + k/2*r*u12*(1-u12) + k/2/h^2*D*(LX*u12 + LY*u)
    % Each step advances by k/2 time step 

    b = ( u + k/2*r*u.*(1-u) + k/2/h^2*D*LY*u ); % En X: on calcule les donnee du systeme implicite AX*u12 = b
                                                   % Along X: evaluate data b for the implicite system AX*u12 = b            
    b = reshape(reshape(b,J2,J1)',J,1);            % On convertit les donnees en indices X-dominant
                                                   % Convert data to X-dominant index             
    u12 = AX\b;                                    % On resout
                                                   % solve
    LXu12 = LX*u12;                                % On stocke LX*u12
                                                   % Store LX*u12 for later
    u12 = reshape(reshape(u12,J1,J2)',J,1);        % On reconvertit u12 en Y-dominant
                                                   % Convert back u12 to Y-dominant index
    LXu12 = reshape(reshape(LXu12,J1,J2)',J,1);    % on reconvertit LXu12 en Y-dominant
                                                   % Convert LXu12 to Y-dominant
    u =   AY\( u12 + k/2*r*u12.*(1-u12) + k/2/h^2*D*LXu12 ); % En Y: on resout AY*u = b 
                                                               % Along Y: solve AY*u = b
    % end of ADI 

   
    u(sideleft) = u(sideleft+J2);              % condition Neumann no flux condition
    u(sideright) = u(sideright-J2);
    u(sidetop) = u(sidetop+1);
    u(sidebottom) = u(sidebottom-1);
    if ( fix(10*t) > fix(10*tp) )                  % affichage toutes les 0.1 unite de temps/display every 0.1 time unit
        image(x,y,255*reshape(u,J2,J1));
    end
    tp = t;
    t = t + k;
    fprintf("t = %.5f\n",t);
end