%% ============================================================
%  Unsteady Couette Flow solved by Crank-Nicolson method
%
%  控制方程：
%      du/dt = (1/Re_D) d2u/dy2
%
%  边界条件：
%      u(0,t) = 0
%      u(1,t) = 1
%
%  初始条件：
%      u(y,0) = 0,  0 <= y < 1
%      u(1,0) = 1
%
%  稳态解析解：
%      u_exact = y
%% ============================================================

clear; clc; close all;

%% ---------------- 参数设置 ----------------
ReD = 5000;          % 基于板间距 D 的雷诺数
N   = 80;            % 区间划分数，网格点数为 N+1
dy  = 1 / N;         % 无量纲网格间距

E   = 1;             % E = dt / (ReD*dy^2)，控制时间步长
dt  = E * ReD * dy^2;

tEnd = 4e6;         % 无量纲终止时间
Nt   = ceil(tEnd/dt);

plotEvery = 1;       % 每隔多少步刷新一次图像

%% ---------------- 网格 ----------------
y = linspace(0,1,N+1)';      % y/D
% u = zeros(N+1,1);            % u/u_e
u = 2*rand(N+1,1); 

% 边界条件
u(1)   = 0;                  % 下壁面
u(end) = 1;                  % 上壁面

u_exact = y;                 % 稳态解析解 u/u_e = y/D

%% ---------------- CN格式矩阵 ----------------
% 控制方程：
%   du/dt = alpha d2u/dy2, alpha = 1/ReD
%
% Crank-Nicolson:
%   (u^{n+1}-u^n)/dt
% = 1/(2 ReD) [Dyy u^{n+1} + Dyy u^n]


% 只求解内部点 j = 2,3,...,N
M = N - 1;

mainL = (1 + E) * ones(M,1);
offL  = (-E / 2) * ones(M,1);

A = spdiags([offL mainL offL], [-1 0 1], M, M);

mainR = (1 - E) * ones(M,1);
offR  = (E/2) * ones(M,1);

B = spdiags([offR mainR offR], [-1 0 1], M, M);

%% ---------------- 作图初始化 ----------------
figure('Color','w','Position',[100 100 1200 500]);

subplot(1,2,1);
h_num = plot(u, y, 'LineWidth', 2); hold on;
h_exa = plot(u_exact, y, 'k--', 'LineWidth', 1.8);
h_pts = plot(u, y, 'o', 'MarkerSize', 4);

xlabel('u/u_e');
ylabel('y/D');
title('Crank-Nicolson solution of Couette flow');
legend('Numerical solution','Steady exact solution','Grid values', ...
       'Location','southeast');
grid on;
axis([0 1.05 0 1]);

subplot(1,2,2);
h_err = semilogy(0, norm(u-u_exact,inf), 'o-', 'LineWidth', 1.8);
xlabel('time step n');
ylabel('max error to steady solution');
title('Convergence to steady state');
grid on;
hold on;

err_history = [];
time_history = [];

%% ---------------- 时间推进 ----------------
for n = 1:Nt

    u_old = u;

    % 右端项：B * u^n
    rhs = B * u_old(2:N);

    % 边界条件贡献
    % 下边界 u(1)=0，对第一行无贡献
    rhs(1) = rhs(1) + (E/2) * u_old(1) + (E/2) * u(1);

    % 上边界 u(N+1)=1，对最后一行有贡献
    rhs(end) = rhs(end) + (E/2) * u_old(end) + (E/2) * u(end);

    % 求解内部点
    u(2:N) = A \ rhs;

    % 强制边界条件
    u(1)   = 0;
    u(end) = 1;

    % 误差记录
    err = norm(u - u_exact, inf);
    err_history(end+1) = err;
    time_history(end+1) = n;

    % ---------------- 动态显示 ----------------
    if mod(n, plotEvery) == 0 || n == 1

        subplot(1,2,1);
        set(h_num, 'XData', u, 'YData', y);
        set(h_pts, 'XData', u, 'YData', y);

        title(sprintf(['Crank-Nicolson Couette Flow\n', ...
                       'Re_D = %g, E = %.2f, dt = %.2f, step = %d'], ...
                       ReD, E, dt, n));

        subplot(1,2,2);
        set(h_err, 'XData', time_history, 'YData', err_history);
        xlim([1 max(10,n)]);
        ylim([max(min(err_history)*0.5,1e-8), 1]);

        drawnow;
    end

    % 若已经非常接近稳态，提前停止
    if err < 1e-6
        fprintf('Converged at step %d, time = %.4f, max error = %.3e\n', ...
                n, n*dt, err);
        break;
    end
end

% ---------------- 最终结果 ----------------
fprintf('\n===== Final Result =====\n');
fprintf('Re_D      = %.2f\n', ReD);
fprintf('N         = %d\n', N);
fprintf('dy        = %.5f\n', dy);
fprintf('E         = %.2f\n', E);
fprintf('dt        = %.5f\n', dt);
fprintf('steps     = %d\n', n);
fprintf('final err = %.3e\n', err);