%FINAL PROJECT #9--------------COMPUTATION OF LINEARIZED EULER EQUATIONS
%Lance Simms 12/2005
%
%u=velocity in x direction v=velocity in y direction
%Flow is flowing in from left. Assumed to be uniform rho=1, u=M_inf, v=0
%Linearization about p=1, u=M_inf, v=0
%
%Thin airfoil present having shape given by
% ywall=Tx(1-x)/2 0<x<1
% ywall=0 x<0 x>1
%
%Euler Equations are written as
%
% DQ DQ DQ
% -- + A -- + B -- = 0
% Dt Dx Dy
%
% where: | 0 1 0 | | rho |
% A = | -M_inf^2 2 M_inf 0 |, Q = | rho u |
% | 0 0 M_inf | | rho v |
%
% | 0 0 1 |
% B = | 0 0 M_inf |
% | 1 0 0 |
% Inputs:
% M_inf : Free Stream Mach number
% nx : Number of points on airfoil (x=0,1)
% cfl : dt = min(dx,dy)*cfl/(M_inf+1 + 1)
% where Spectral radius of A is (M_inf+1) and B is 1
% nmax : Number of time steps to go
% nout : Output frequency of plot of cp (0 for no plots)
% meth : 1: Explicit Euler / 2: RK4 / 3: Implicit Euler
% order : 0 : central 2nd / 1 : 1st order upwind / 2 : 2nd order upwind
% ichar : 0 : simple bc / 1 : char bc / -1 simpler bc
%
%
clear;

%MACH NUMBER

M_inf=.8;
M_inf=input('Enter M_inf');

%Grid Generation
xmin=-1.0; xmax=3.0;
ymin=0.0; ymax=2.0;

nx=10;
nx=input('Enter number of points on body:');

dx=1.0/(nx-1); dy=dx;
x=[xmin:dx:xmax]'; y=[ymin:dy:ymax]';
jmax=size(x,1); kmax=size(y,1);

%Generate meshgrid having x values as rows and y values as columns
[X,Y]=meshgrid(x,y);

jle=nx; jte=2*nx-1;
fprintf('Leading edge point= %g trailing edge point = %g \n', jle,jte);
fprintf('X of Leading edge point =%g X of trailing edge point = %g \n', x(jle),x(jte));

%Surface Definition: tau is the thickness
tau=0.1;
%tau=input ('tau);
s=[jle:jte];
ys=zeros(jmax,1); ys(s)=tau*0.5*x(s).*(1.0-x(s));

%Try to figure out indices of airfoil
foil=0;
for xind=1:length(x)
newfoil=find(y<ys(xind));
foil=[foil,xind*length(y)+newfoil'];
end
foil=foil(2:length(foil));

%Initialze elements of Q
Q(:,:,1)=ones(jmax,kmax,1); %rho
Q(:,:,2)=M_inf*ones(jmax,kmax,1); %rho u
Q(:,:,3)=zeros(jmax,kmax,1); %rho v

%Initialize storage for Runge-Kutta Methods
Q1=zeros(jmax,kmax,3);
Q2=zeros(jmax,kmax,3);
Q3=zeros(jmax,kmax,3);

%Jacobian Matrices
A=zeros(3,3); B=zeros(3,3);

A(1,:)=[0, 1, 0 ]; B(1,:)=[0, 0, 1 ];
A(2,:)=[-M_inf^2+1 2.0*M_inf, 0 ]; B(2,:)=[0, 0, M_inf];
A(3,:)=[0, 0, M_inf]; B(3,:)=[1, 0, 0 ];

%FLUX SPLIT
[evecA,evalA]=eig(A); [evecB,evalB]=eig(B);
evalAp=(evalA+abs(evalA))/2; evalBp=(evalB+abs(evalB))/2;
evalAm=(evalA-abs(evalA))/2; evalBm=(evalB-abs(evalB))/2;

Ap=evecA*evalAp*inv(evecA); Bp=evecB*evalBp*inv(evecB);
Am=evecA*evalAm*inv(evecA); Bm=evecB*evalBm*inv(evecB);

%Thin airfoil BC v=M_inf D(ys)/Dx
Q(s,1,3)=M_inf*(ys(s+1)-ys(s-1))*0.5/dx;

%Indices
J=[1:jmax]; K=[1:kmax];
JM=[2:jmax-1]; KM=[2:kmax-1];

%INPUTS
nmax=300;
nmax=input('Enter nmax ');
nmax_tot=nmax;

%cfl=0.5;
cfl=input('Enter cfl');

meth=0;
meth=input('Enter 0 for RK4 \n 1 for Euler Explicit\n 2 for Euler implicit: ');

order=1;
epse=0.0;
%epse=Enter('2nd order diss coef epse ',epse);

Eig_t=2.0+M_inf;

dt=min(dx,dy)*cfl/Eig_t;

ncount=0;
cpu0=cputime;

% Initial Residual Computation
nstep=1;
R=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
resid_0 = (norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax);

figure;
clf;
set (gcf,'doublebuffer','on');

%Density Plot
subplot(2,1,1);
im=imagesc(x,y,Q(:,:,1));
maxp=text(0,1.08,'p_{max}=','units','normalized');
maxpu=text(.22,1.08,'pu_{max}=','units','normalized');
maxpv=text(.46,1.08,'pv_{max}=','units','normalized');
CFLtext=text(.7,1.10,'CFL=','units','normalized');
set(CFLtext,'string',sprintf('CFL=%2.2f',cfl));
numpoints=text(.88,1.10,'npoints=','units','normalized');
set(numpoints,'string',sprintf('npoints=%2.2f',nx));
xlabel('x');ylabel('y');
colorbar;
colorbartext=text(1.1,-0.1,'\rho','units','normalized');

%Residual Plot
subplot(2,1,2);
residp=semilogy([1:nstep],resid_0);
xlabel('n');
ylabel('resid');
residtext=text(.6,.5,'Residual=','units','normalized');
set(residtext,'string',sprintf('Residual=%2.8f',resid_0));
residslopetext=text(.6,.62,'Slope=','units','normalized');

%set (gcf,'Units','normalized');
%set (gcf,'Position',[-1,0,.8,.8]);
%subplot(2,1,2);

%%RUNGE-KUTTA 4
if (meth==0)
for nstep=1:nmax_tot

resid(nstep) = ...
(norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax);
resid(nstep) = resid(nstep)/resid_0;
ncount=ncount+1;

Q1=zeros(jmax,kmax,3);
Q2=zeros(jmax,kmax,3);
Q3=zeros(jmax,kmax,3);

%First RK4 Step: Q'n+1/2=Qn+.5hR(Q)n
R1=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
Q1=Q+0.5*dt*R1;

%Second RK4 Step: Q-n+1/2=Qn+.5hR(Q')n
R2=rhs_flux_split_so(Q1,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
Q2=Q+0.5*dt*R2;

%Third RK4 Step: Q/n+1=Qn+hR(Q-)n+1/2
R3=rhs_flux_split_so(Q2,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
Q3=Q+dt*R3;

%Last RK4 Step: Qn=1/6*h*[R(Q)n+2(R(Q')n+1/2+R(Q-)n+1/2)+R(Q)/n+1]
R4=rhs_flux_split_so(Q3,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
Q(JM,KM,:) = Q(JM,KM,:) + (1/6)*dt*(R1(JM,KM,:)+2*(R2(JM,KM,:)+R3(JM,KM,:))+R4(JM,KM,:));

R=(1/6)*dt*(R1(JM,KM,:)+2*(R2(JM,KM,:)+R3(JM,KM,:))+R4(JM,KM,:));
%Apply Boundary conditions
Q=bc8(Q,JM,KM,jmax,kmax,M_inf);

PlotQ=Q(:,:,1)';
PlotQ(foil)=1;
%subplot(2,1,2);

set(gcf,'CurrentObject',im);
%imagesc(x,y,flipud(PlotQ));
%imagesc(x,y(kmax:-1:1),Q(:,:,1)');
set(gcf,'units','normalized');
set(im,'cdata',flipud(PlotQ),'xdata',x,'ydata',y);
set(maxp,'string',sprintf('p_{max}=%2.3f',max(max(Q(:,:,1)))));
set(maxpu,'string',sprintf('pu_{max}=%2.3f',max(max(Q(:,:,2)))));
set(maxpv,'string',sprintf('pv_{max}=%2.3f',max(max(Q(:,:,3)))));
set(residtext,'string',sprintf('Residual=%2.12f',resid(nstep)));

% set(gcf,'CurrentObject',resid);
set(residp,'XData',[1:nstep],'YData',resid);
if nstep>2
residslope=resid(nstep)-resid(nstep-1);
set(residslopetext,'string',sprintf('Slope=%2.12f',residslope));
end

refresh;
drawnow;
%[flag]=figflag(gcf,0);
pause(0.1);
end
end

%%EULER EXPLICIT
if (meth==1)

for nstep=1:nmax_tot

resid(nstep) = ...
(norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax);
resid(nstep) = resid(nstep)/resid_0;

R=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
Q(JM,KM,:) = Q(JM,KM,:) + dt*R(JM,KM,:);

%Boundary Conditions
Q=bc8(Q,JM,KM,jmax,kmax,M_inf);

%Position Airfoil in Plot
PlotQ=Q(:,:,1)';
PlotQ(foil)=1;

set(gcf,'CurrentObject',im);
set(im,'cdata',flipud(PlotQ),'xdata',x,'ydata',y);
set(maxp,'string',sprintf('p_{max}=%2.3f',max(max(Q(:,:,1)))));
set(maxpu,'string',sprintf('pu_{max}=%2.3f',max(max(Q(:,:,2)))));
set(maxpv,'string',sprintf('pv_{max}=%2.3f',max(max(Q(:,:,3)))));
% set(gca,'title',text('string',sprintf('b=%2.2f',max(max(max(Q(:,:,1)))),...
% 'color','r','fontsize',14)));
% set(gcf,'CurrentObject',resid);
set(residp,'XData',[1:nstep],'YData',resid);
set(residtext,'string',sprintf('Residual=%2.9f',resid(nstep)));
if nstep>2
residslope=resid(nstep)-resid(nstep-1);
set(residslopetext,'string',sprintf('Slope=%2.9f',residslope));
end
colorbar;
refresh;
drawnow;
pause(0.1);
end
end

%%EULER IMPLICIT
if meth == 2
% Delta Form Implicit Approximate Factorization
[a_x,b_x,c_x,d_x,e_x] = fill_x(JM,KM,A,dx,dt,jmax,kmax,epse);
[a_y,b_y,c_y,d_y,e_y] = fill_y(JM,KM,B,dy,dt,jmax,kmax,epse);

for nstep=1:nmax_tot

resid(nstep) = ...
(norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax);
resid(nstep) = resid(nstep)/resid_0;

% [I + dt A Del_cen_x ]^(-1) R
% First Step of Implicit Factored Scheme : x block tridiagonal step
R=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
R = dt*R;
R = btrix_3x3(2,jmax-1,KM,b_x,c_x,d_x,R);

% [I + dt B Del_cen_y ]^(-1) R
% Second Step of Implicit Factored Scheme : y block tridiagonal step
DelQ = btriy_3x3(2,kmax-1,JM,b_y,c_y,d_y,R);
% Q^(n+1) = Q^n + Delta Q^n
Q(JM,KM,:) = Q(JM,KM,:) + DelQ(JM,KM,:);

%Boundary Conditions
Q=bc8(Q,JM,KM,jmax,kmax,M_inf);

%subplot(2,1,1);
%quiver(X,Y,Q(:,:,2)',Q(:,:,3)');
%contourf(X,Y,Q(:,:,2)',50);
%xlabel('rho');
PlotQ=Q(:,:,1)';
PlotQ(foil)=1;
%subplot(2,1,2);
%imagesc(x,y,flipud(PlotQ));
%imagesc(x,y(kmax:-1:1),Q(:,:,1)');
%set(gcf,'CurrentObject',im);
%imagesc(x,y,flipud(PlotQ));
%imagesc(x,y(kmax:-1:1),Q(:,:,1)');
set(im,'cdata',flipud(PlotQ),'xdata',x,'ydata',y);
set(maxp,'string',sprintf('p_{max}=%2.3f',max(max(Q(:,:,1)))));
set(maxpu,'string',sprintf('pu_{max}=%2.3f',max(max(Q(:,:,2)))));
set(maxpv,'string',sprintf('pv_{max}=%2.3f',max(max(Q(:,:,3)))));
set(residtext,'string',sprintf('Residual=%2.9',resid(nstep)));
%text(.7,1.1,sprintf('max_p=%2.2f',max(max(Q(:,:,1)))));
streamslice(X,Y,Q(:,:,2)',Q(:,:,3)')
% set(gcf,'CurrentObject',resid);
set(residp,'XData',[1:nstep],'YData',resid);
if nstep>2
residslope=resid(nstep)-resid(nstep-1);
set(residslopetext,'string',sprintf('Slope=%2.9f',residslope));
end
refresh;
drawnow;
pause(0.1);
end
end