%% Evolving the variable diffusivity Diffusion equation
% A tutorial of the variable diffusivity diffusion equation
% based on Script 36 from Trefethen's Book. This version steps using LU
% decomposition. Timing is done separately in the accompanying script. This
% is the continuation of the diff_bc_step script.
% Please note that this script defines functions at the end, which is only
% supported by MATLAB 2016b or later.
%%
clear all,close all
%% Set up grid and 2D Laplacian, boundary points included:
% Notice the low number of points.
N = 64;
%%
% Use the cheb function to get the Chebyshev points.
[D,x] = cheb(N);
y = x;
%%
% Create the 2D grid.
[xx,yy] = meshgrid(x,y);
Lx=1;
Ly=1;
[xxx,yyy] = meshgrid(linspace(-Lx,Lx,500),linspace(-Ly,Ly,500));
%% Create the differentiation matrices
% Make them sparse to save space.
D2 = sparse(D^2); 
I = speye(N+1);
% Use the outer product to create the matrices.
Dy=kron(I,D); 
Dx=kron(D,I);
Dyy=kron(I,D2); 
Dxx=kron(D2,I);
%% Create variable diffusion constant
r=sqrt((xx-0.5).^2+(yy+0.75).^2);
kappa0=1e-2;
kappa=kappa0*(1+0.5*sech((r/0.5).^2));
kappav=kappa(:);
mykappax=Dx*kappa(:);
mykappay=Dy*kappa(:);
kappaxv=mykappax(:);
kappayv=mykappay(:);
kappavmat=sparse(diag(kappav));
kappaxmat=sparse(diag(kappaxv));
kappaymat=sparse(diag(kappayv));
%% Create the Laplacian operator
% Combine the last two parts and build the Laplacian operator.
Llap = kappavmat*(Dxx+Dyy)+kappaxmat*Dx+kappaymat*Dy;
%% Create the time evolving Laplace operator
dt=0.1; 
L=speye(size(Llap))-dt*Llap;
%% Boundary Conditions
% Impose boundary conditions by replacing appropriate rows of L:
bxp1 = find(xx==1); % Rows of the big matrix at which x=1.
bxn1 = find(xx==-1); % Rows of the big matrix at which x=-1.
byp1 = find(yy==1); % Rows of the big matrix at which y=1.
byn1 = find(yy==-1); % Rows of the big matrix at which y=-1.
u = zeros((N+1)^2,1); 
rhs=u;
L(bxp1,:) = zeros(N+1,(N+1)^2);
L(bxp1,bxp1) = eye(N+1);
rhs(bxp1)=-yy(bxp1);
% Dirichlet BC at x=-1 notice how RHS is modified.
for ii=1:length(byp1)
    L(byp1(ii),:) = Dy(byp1(ii),:); 
    rhs(byp1)=2.5; % Set L=Dy at y=1, notice RHS is set to zero above.
    L(bxn1(ii),:) = Dx(bxn1(ii),:); 
    rhs(bxn1)=sin(pi*(1+yy(bxn1))); % Set L=Dy at x=-1, notice RHS modifed.
end
L(byn1,:) = zeros(N+1,(N+1)^2); 
L(byn1,byn1) = eye(N+1); 
rhs(byn1)=-cos(3*pi*xx(byn1));
%% Solve the equation
% Factor the Laplace equation.
[ll,uu,pp,qq]= lu(L);
numsteps=10;
numouts=100;
t=0;
for ii=1:numouts
    % Plotting here shows the ICs.
    uuu = interp2(xx,yy,reshape(u,size(xx)),xxx,yyy,'spline');
    figure(1),clf
    pcolor(xxx,yyy,uuu),  shading flat,caxis([-1 1]*1),colorbar;
    title(t)
    drawnow
    % Solve the equation for the new u.
    for jj=1:numsteps
        t=t+dt;
        u = qq*(uu\(ll\(pp*rhs)));
        rhs=u;
    end
    
end
% Plot one last time.
uuu = interp2(xx,yy,reshape(u,size(xx)),xxx,yyy,'spline');
figure(1),clf
pcolor(xxx,yyy,uuu),  shading flat,caxis([-1 1]*1),colorbar;
title(t)
drawnow
%% Chebyshev function
% CHEB  compute D = differentiation matrix, x = Chebyshev grid
function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
D  = D - diag(sum(D'));                 % diagonal entries
end
