[m10_2]Implicit form solution of 2Dtemperature-equation.(variable k,den*Cp)

%This program complished the process of 2D Eulerian advection with method
%of allocating memory temporarily.
clear all;clf;
% Set parameters
Nt = 1;
output_interval = 1;
Lx = 1000000; 
Ly = 1500000;
x_midd = Lx/2;
y_midd = Ly/2;
Nx = 51; 
Ny = 31; 
kw = 10;%W/(m*K)
denw = 3300;%Kg/m^3
Cpw = 1100;%J/(Kg*K)
ks = 3;%W/(m*K)
dens = 3200;%Kg/m^3
Cps = 1000;%J/(Kg*K)
bt = 1000;%K,background temperature
wt = 1300;%K,wave temperature
Nx1 = Nx+1;
Ny1 = Ny+1;
dx = Lx / (Nx-1);
dy = Ly / (Ny-1);
% Basic E-nodes.
x = 0:dx:Lx;
y = 0:dy:Ly;
[X,Y] = meshgrid(0:dx:Lx,0:dy:Ly);
% Initialize k,den*Cp and temperature matrix at the temperature points. 
kh = zeros(Ny,Nx);
hc = zeros(Ny,Nx);
td = zeros(Ny,Nx);
T = zeros(Ny,Nx);
for j=1:1:Nx
    for i=1:1:Ny
        T(i,j) = bt;
        kh(i,j) = ks;
        hc(i,j) = dens*Cps;%unit:(Pa*s)
        if ( x(j)<=x_midd+100000 && x(j)>=x_midd-100000 && y(i)<=y_midd+150000 && y(i)>=y_midd-150000)
            kh(i,j) = kw;
            hc(i,j) = denw*Cpw;%unit:(Pa*s)
            T(i,j) = wt;
        end
        td(i,j) = kh(i,j)/hc(i,j);
    end
end
% 2.7548e-06(suround),9.3750e-07(wave)
dt = min(dx^2,dy^2)/(4*td(2,2));

figure(1);
colormap('Jet');
mesh(X,Y,T);
xlabel('Horizontal(m)')
ylabel('Vertical(m)')
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading interp;
colorbar
title(['T after 0 deltat'])
pause(0.01);
%grids has N2 rows and N1 columns
[x,y] = meshgrid(0:Lx/(Nx-1):Lx,0:Ly/(Ny-1):Ly);
b = zeros(Nx*Ny, 1);
coe = zeros(Nx*Ny, Nx*Ny);
b1 = zeros(Nx*Ny, 1);
coe1 = zeros(Nx*Ny, Nx*Ny);

% Constant temperature boundary condition.
% Tbc = constant = 1000
% up and down
for j = 1:Nx
    i = 1;
    k = (j-1)*Ny + i;
    coe(k, k) = 1;
    b(k, 1) = 1000;
    i = Ny;
    k = (j-1)*Ny + i;
    coe(k, k) = 1;  
    b(k, 1) = 1000;
end
% left and right
for i = 1:Ny
    j = 1;
    k = (j-1)*Ny + i;
    coe(k, k) = 1;
    b(k, 1) = 1000;
    j = Nx;
    k = (j-1)*Ny + i;
    coe(k, k) = 1;
    b(k, 1) = 1000;
end 
% coefficients of interval points
for j = 2:Nx-1
    for i = 2:Ny-1
        k = (j-1)*Ny + i; % k is the index
        coe(k, k) = (kh(i,j)+kh(i,j+1))/(2.*dx.^2)+(kh(i,j-1)+kh(i,j))/(2.*dx.^2)...
                   +(kh(i,j)+kh(i+1,j))/(2.*dy.^2)+(kh(i-1,j)+kh(i,j))/(2.*dy.^2)...
                   +hc(i,j)/dt; % middle T3
        coe(k, k-1) = -(kh(i-1,j)+kh(i,j))/(2.*dy.^2); % up T2
        coe(k, k+1) = -(kh(i,j)+kh(i+1,j))/(2.*dy.^2); % down T4
        coe(k, k-Ny) = -(kh(i,j-1)+kh(i,j))/(2.*dx.^2); % left T1
        coe(k, k+Ny) = -(kh(i,j)+kh(i,j+1))/(2.*dx.^2); % right T5
        b(k, 1) = T(i,j)*hc(i,j)/dt;
    end
end
%direct method:right martix is divided by coe martix on the left
u = coe \ b;
%transform vector to grids 
%grids has N2 rows and N1 columns
U = reshape(u, Ny, Nx);
figure(1);
mesh(x,y,U);
xlabel('Horizontal(Km)');
ylabel('Vertical(Km)');
zlabel('Gravitational potential');
title('(CTbc) Inplicit formsolution of T')
% Insulating boundary condition.
for j = 2:Nx-1
    i = 1;
    k = (j-1)*Ny + i;
    % T(1,j) - T(2,j) = 0;
    coe1(k, k) = 1;
    coe1(k, k+1) = -1;
    b1(k, 1) = 0;
    i = Ny;
    k = (j-1)*Ny + i;
    % T(Ny,j) - T(Ny-1,j) = 0;
    coe1(k, k) = 1;  
    coe1(k, k-1) = -1;
    b1(k, 1) = 0;
end
% left and right
for i = 1:Ny
    j = 1;
    k = (j-1)*Ny + i;
    coe1(k, k) = 1;
    coe1(k, k+Ny) = -1;
    b1(k, 1) = 0;
    j = Nx;
    k = (j-1)*Ny + i;
    coe1(k, k) = 1;
    coe1(k, k-Ny) = -1;
    b1(k, 1) = 0;
end 
% coefficients of interval points
for j = 2:Nx-1
    for i = 2:Ny-1
        k = (j-1)*Ny + i; % k is the index
        coe1(k, k) = (kh(i,j)+kh(i,j+1))/(2.*dx.^2)+(kh(i,j-1)+kh(i,j))/(2.*dx.^2)...
                   +(kh(i,j)+kh(i+1,j))/(2.*dy.^2)+(kh(i-1,j)+kh(i,j))/(2.*dy.^2)...
                   +hc(i,j)/dt; % middle T3
        coe1(k, k-1) = -(kh(i-1,j)+kh(i,j))/(2.*dy.^2); % up T2
        coe1(k, k+1) = -(kh(i,j)+kh(i+1,j))/(2.*dy.^2); % down T4
        coe1(k, k-Ny) = -(kh(i,j-1)+kh(i,j))/(2.*dx.^2); % left T1
        coe1(k, k+Ny) = -(kh(i,j)+kh(i,j+1))/(2.*dx.^2); % right T5
        b1(k, 1) = T(i,j)*hc(i,j)/dt;
    end
end
u1 = coe1 \ b1;
%transform vector to grids 
%grids has N2 rows and N1 columns
U1 = reshape(u1, Ny, Nx);
figure(2);
mesh(x,y,U1);
xlabel('Horizontal(Km)');
ylabel('Vertical(Km)');
zlabel('Gravitational potential');
title('(Ibc) Inplicit form solution of T')