%% Boundary conditions with diffusion operator
% A tutorial version of Script 36 from Spectral Methods in Matlab by Trefethen. This solves the
% variables coefficient diffusion equation on a box with non-trivial
% boundary conditions. We also have the "\" matlab operator explain its
% logic when solving the problem. There is no time dependence.
% Please note that this script defines functions at the end, which is only
% supported by MATLAB 2016b or later. NOTE: IF YOU WOULD LIKE THIS SCRIPT
% TO WORK WITH OLDER VERSIONS OF MATLAB REMOVE THE FUNCTION FROM THE END
% AND GIVE IT ITS OWN .m FILE.
%%
clear,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);
[xxx,yyy] = meshgrid(-1:.04:1,-1:.04:1);
%% 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.
L = kappavmat*(Dxx+Dyy)+kappaxmat*Dx+kappaymat*Dy;
%% 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
rhs = zeros((N+1)^2,1);
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 modified
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
% Use spparms to show the logic used by Matlab to solve the linear system
% with the "\" command.
spparms('spumoni',2)
%%
% Just solve the equations with nothing extra.
tic
u1 = L\rhs; 
uu1 = reshape(u1,N+1,N+1);
justsolve=toc
%%
% Solve the equations with factoring.
tic
[ll,uu,pp,qq]= lu(L);
factortime=toc
tic
u2 = qq*(uu\(ll\(pp*rhs))); 
uu2 = reshape(u2,N+1,N+1);
factorsolve=toc
%%
% The difference between the solutions.
figure(1)
clf
pcolor(xx,yy,uu1-uu2),shading interp,colorbar
title('Difference between solution without and with factoring')
%%
% The individual solutions.
figure(2)
clf
subplot(2,1,1)
pcolor(xx,yy,uu1),shading interp,caxis([-1 1])
subplot(2,1,2)
pcolor(xx,yy,uu2),shading interp,caxis([-1 1])
%%
% The diffusivity.
figure(3)
clf
pcolor(xx,yy,kappa),shading flat,colorbar
ylabel('y')
%%
% Create the derivatives.
u1y=Dy*u1;
u1x=Dx*u1;
uy=reshape(u1y,N+1,N+1);
ux=reshape(u1x,N+1,N+1);
%%
figure(4)
clf
% Plot the derivatives.
subplot(2,2,1)
pcolor(xx,yy,ux),shading interp,caxis([-10 10])
subplot(2,2,2)
pcolor(xx,yy,uy),shading interp,caxis([-10 10])
subplot(2,2,3)
% Plot the derivatives at the boundaries.
plot(yy(bxn1),ux(bxn1),'bo-',xx(byp1),uy(byp1),'rs-')
subplot(2,2,4)
plot(xx(byn1),u1(byn1),'bo-',yy(bxp1),u1(bxp1),'rs-')

%% 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