[m10_4](Ibc) Inplicit form solution of T(advection) with method of marker-in-cell(recycle marker).

%This program complished the process of 2D Eulerian advection with method
%of allocating memory temporarily.
clear all;clf;

Nt = 100;
output_interval = 1;
% BASIC NODES
Lx = 1000000; 
Ly = 1500000;
x_midd = Lx/2;
y_midd = Ly/2;
Nx = 51; 
Ny = 31; 
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;
g_y = 10;
% PARAMETERS
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
% L-MARKERS
num_per_cell = 4;
Nx_marker = (Nx-1)*num_per_cell;
Ny_marker = (Ny-1)*num_per_cell;
marknum=Nx_marker*Ny_marker;
x_m_step = Lx/Nx_marker;
y_m_step = Ly/Ny_marker;
% Marker points.
xm=zeros(1,marknum); % Horizontal coordinates, m
ym=zeros(1,marknum);
rhom=zeros(1,marknum); % K, kg/m^3  kt
etam=zeros(1,marknum); % den*Cp     hct
tdm = zeros(1,marknum); % K/den/Cp  td
Tm = zeros(1,marknum); % temperature
% INITIALIZE K,Den*Cp in the L-markers.
m = 1;
for jm=1:1:Nx_marker
    for im=1:1:Ny_marker
        % Define marker coordinates
        xm(m)=x_m_step/2+(jm-1)*x_m_step+(rand-0.5)*x_m_step;
        ym(m)=y_m_step/2+(im-1)*y_m_step+(rand-0.5)*y_m_step;
        % Marker properties
        if ( xm(m)<=x_midd+100000 && xm(m)>=x_midd-100000 && ym(m)<=y_midd+150000 && ym(m)>=y_midd-150000)
            rhom(m)=kw; % K wave
            etam(m)=denw*Cpw; % den*Cp wave
            Tm(m) = wt;
        else
            rhom(m)=ks; % K surround
            etam(m)=dens*Cps; % den*Cp surround
            Tm(m) = bt;
        end
        % Update marker counter
        m=m+1;
    end
end
% Define matrix
hct = zeros(Ny,Nx);
kht = zeros(Ny,Nx);
T = zeros(Ny,Nx);
td = zeros(Ny,Nx);
b = zeros(Nx*Ny, 1);
coe = zeros(Nx*Ny, Nx*Ny);
vx = zeros(Ny,Nx);
vy = zeros(Ny,Nx);

for j = 1:Nx
    for i = 1:Ny
        % Constant velocity
        vx(i,j) = 1e-9;
        vy(i,j) = 1e-9;
    end
end

% DEFINE dt
% 2.7548e-06(suround),9.3750e-07(wave)
dtexp = min(dx,dy)^2/(3*td(2,2));
if(vx(1,1)~=0)
    dtexp=min(dtexp,abs(dx/vx(1,1))); % Limitation for horizontal advection timestep
end
if(vy(1,1)~=0)
    dtexp=min(dtexp,abs(dy/vy(1,1))); % Limitation for vertical advection timestep
end
dt = 0.7*dtexp;

for t = 1:Nt
    % Indispensable
    % Eq.8.18
    w_m_x_node = zeros(Ny,Nx);
    den_w_m = zeros(Ny,Nx);
    vis_w_m = zeros(Ny,Nx);
    % Eq.10.12
    T_w_m = zeros(Ny,Nx);
    w_m_t = zeros(Ny,Nx);
    % Searching the coordinate of upper-left node of every L-marker.
    for m=1:1:marknum
        % Define i,j indexes for the upper left BASIC NODE(Nx*Ny)
        j=fix((xm(m)-x(1))/dx)+1;
        i=fix((ym(m)-y(1))/dy)+1;
        if(j<1)
            j=1;
        elseif(j>Nx-1)
            j=Nx-1;
        end
        if(i<1)
            i=1;
        elseif(i>Ny-1)
            i=Ny-1;
        end
        % Compute distances between marker and upper-left node.
        dis_x = xm(m)-x(j);
        dis_y = ym(m)-y(i);
        wtmij = (1-dis_x/dx)*(1-dis_y/dy);
        wtmi1j = (1-dis_x/dx)*(dis_y/dy);
        wtmij1 = (dis_x/dx)*(1-dis_y/dy);
        wtmi1j1 = (dis_x/dx)*(dis_y/dy);
        
        % Sum of (K/Den*Cp)*(weight) and sum of weight. 
        w_m_x_node(i,j) = w_m_x_node(i,j) + wtmij;   % denominator% numerator
        den_w_m(i,j) = den_w_m(i,j) + rhom(m).*wtmij;% K
        vis_w_m(i,j) = vis_w_m(i,j) + etam(m).*wtmij;% Den*Cp
        T_w_m(i,j) = T_w_m(i,j) + Tm(m).*etam(m)*wtmij;% T
        w_m_t(i,j) = w_m_t(i,j) + etam(m)*wtmij;
        
        w_m_x_node(i+1,j) = w_m_x_node(i+1,j) + wtmi1j;
        den_w_m(i+1,j) = den_w_m(i+1,j) + rhom(m).*wtmi1j;
        vis_w_m(i+1,j) = vis_w_m(i+1,j) + etam(m).*wtmi1j;
        T_w_m(i+1,j)   = T_w_m(i+1,j) + Tm(m).*wtmi1j*etam(m);% T
        w_m_t(i+1,j) = w_m_t(i+1,j) + etam(m)*wtmi1j;
        
        w_m_x_node(i,j+1) = w_m_x_node(i,j+1) + wtmij1;
        den_w_m(i,j+1) = den_w_m(i,j+1) + rhom(m).*wtmij1;
        vis_w_m(i,j+1) = vis_w_m(i,j+1) + etam(m).*wtmij1;
        T_w_m(i,j+1)   = T_w_m(i,j+1) + Tm(m).*wtmij1*etam(m);% T
        w_m_t(i,j+1) = w_m_t(i,j+1) + etam(m)*wtmij1;
        
        w_m_x_node(i+1,j+1) = w_m_x_node(i+1,j+1) + wtmi1j1;
        den_w_m(i+1,j+1) = den_w_m(i+1,j+1) + rhom(m).*wtmi1j1;
        vis_w_m(i+1,j+1) = vis_w_m(i+1,j+1) + etam(m).*wtmi1j1;
        T_w_m(i+1,j+1)   = T_w_m(i+1,j+1) + Tm(m).*wtmi1j1*etam(m);% T
        w_m_t(i+1,j+1) = w_m_t(i+1,j+1) + etam(m)*wtmi1j1;
    end

    % Computing the value of K and Den*Cp and T of E-nodes. 
    for j = 1:Nx
        for i = 1:Ny
            if(w_m_x_node(i,j)>0)
                kht(i,j) = den_w_m(i,j)/w_m_x_node(i,j);%K
                hct(i,j) = vis_w_m(i,j)/w_m_x_node(i,j);%Den*Cp
            end
            td(i,j) = kht(i,j)/hct(i,j);%K/Den/Cp
            if(w_m_t(i,j)>0)
                T(i,j)  = T_w_m(i,j)/w_m_t(i,j);%T
            end
        end
    end
    % Apply the boundary conditions to T at every t+dt.
    % Constant temperature boundary
    % T = constant
    T(1,:) = 1000;
    T(Ny,:) = 1000;
    T(:,1) = 1000;
    T(:,Nx) = 1000;
    % Insulting boundary condition
    % T1(boudnary) = T2(internal)
%     T(1,2:Nx-1) = T(2,2:Nx-1);
%     T(Ny,2:Nx-1) = T(Ny,2:Nx-1);
%     T(:,1) = T(:,2);
%     T(:,Nx) = T(:,Nx-1);
    
    khtdt = kht;
    hctdt = hct;
    % Change k and den*Cp at internal nodes with upwind difference
    for j = 2:Nx-1
        for i = 2:Ny-1
            khtdt(i,j) = dt*( -vx(i,j)*( kht(i,j)-kht(i,j-1) )/dx -vy(i,j)*( kht(i,j)-kht(i-1,j) )/dy ) +kht(i,j) ;
            hctdt(i,j) = dt*( -vx(i,j)*( hct(i,j)-hct(i,j-1) )/dx -vy(i,j)*( hct(i,j)-hct(i-1,j) )/dy ) +hct(i,j) ;
        end
    end
    kht = khtdt;
    hct = hctdt;
    % 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
            if vx(i,j)>0
                dtdx = (T(i,j)-T(i,j-1))/dx;
            elseif vx(i,j)<0
                dtdx = (T(i,j+1)-T(i,j))/dx;
            end
            if vy(i,j)>0
                dtdy = (T(i,j)-T(i-1,j))/dy;
            elseif vy(i,j)<0
                dtdy = (T(i+1,j)-T(i,j))/dy;
            end
            coe(k, k) = (khtdt(i,j)+khtdt(i,j+1))/(2.*dx.^2)+(khtdt(i,j-1)+khtdt(i,j))/(2.*dx.^2)...
                +(khtdt(i,j)+khtdt(i+1,j))/(2.*dy.^2)+(khtdt(i-1,j)+khtdt(i,j))/(2.*dy.^2)...
                +hctdt(i,j)/dt; % middle T3
            coe(k, k-1) = -(khtdt(i-1,j)+khtdt(i,j))/(2.*dy.^2); % up T2
            coe(k, k+1) = -(khtdt(i,j)+khtdt(i+1,j))/(2.*dy.^2); % down T4
            coe(k, k-Ny) = -(khtdt(i,j-1)+khtdt(i,j))/(2.*dx.^2); % left T1
            coe(k, k+Ny) = -(khtdt(i,j)+khtdt(i,j+1))/(2.*dx.^2); % right T5
            b(k, 1) = T(i,j)*hctdt(i,j)/dt - hctdt(i,j)*(vx(i,j)*dtdx+vy(i,j)*dtdy);
        end
    end
    %direct method:right martix is divided by coe martix on the left
    u = coe \ b;
    %transform vector to grids
    %grids has Ny rows and Nx columns
    U = reshape(u, Ny, Nx);
    T = U;
    figure(1);
    pcolor(x,y,T);colorbar;
    xlabel('Horizontal(Km)');
    ylabel('Vertical(Km)');
    zlabel('Gravitational potential');
    title(['(CTbc) Inplicit form solution of T(advection) after',num2str(t),' deltat'])
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    pause(0.01);
    % Marker transport
    for m = 1:marknum
        vxm = 1e-9;
        vym = 1e-9;
        xm(m) = xm(m)+vxm*dt;
        ym(m) = ym(m)+vym*dt;
        if (xm(m) > Lx)
            xm(m) = xm(m) - Lx;
        elseif(xm(m) < 0)
            xm(m) = xm(m) + Lx;
        end
        if (ym(m) > Ly)
            ym(m) = ym(m) - Ly;
        elseif(ym(m) < 0)
            ym(m) = ym(m) + Ly;
        end
    end
    if(t==1)
        % Interpolate T from E-nodes to L-markers.
        for m = 1:marknum
            % Interpolate T to BASIC NODES(Nx*Ny)
            % Define i,j indexes for the upper left node
            j=fix((xm(m)-x(1))/dx)+1;
            i=fix((ym(m)-y(1))/dy)+1;
            if(j<1)
                j=1;
            elseif(j>Nx-1)
                j=Nx-1;
            end
            if(i<1)
                i=1;
            elseif(i>Ny-1)
                i=Ny-1;
            end
            % Compute distances
            dis_x=xm(m)-x(j);
            dis_y=ym(m)-y(i);
            % Compute weights
            wtmij=(1-dis_x/dx)*(1-dis_y/dy);
            wtmi1j=(1-dis_x/dx)*(dis_y/dy);
            wtmij1=(dis_x/dx)*(1-dis_y/dy);
            wtmi1j1=(dis_x/dx)*(dis_y/dy);
            % Compute Tm of next t+dt
            Tm(m)=T(i,j)*wtmij + T(i+1,j)*wtmi1j + T(i,j+1)*wtmij1 + T(i+1,j+1)*wtmi1j1;
        end
    end

end