% This file can be downloaded at

% http://www.tfd.chalmers.se/~lada/projects/proind.html

% choose "Synthetic Inlet Boundary Conditions" on the left side
%
%

% exemple of 1d Channel flow with a k-omega model [1]. Re=u_tau*h/nu=30000 (h=half 
% channel height).
%
% Discretization described in detail in
% http://www.tfd.chalmers.se/~lada/comp_fluid_dynamics/
%
%
% [1} S-H Peng and L. Davidson and S. Holmberg",
% "A Modified Low-{R}eynolds-Number k-omega Model for Recirculating Flows",
% Journal of Fluids Engng, vol.  119, pp. 867-875, 1997 




close all
clear

% max number of iterations
niter=60000;

% friction velocity u_*=1
% half channel width=1
%

% create the grid

nj=82;
njm1=nj-1;

% equi-distant grid
for j=2:njm1
    yc(j)=2*(j-1)/(nj-2);
end

yc(nj)=0;

% cell centres
for j=2:nj-1
   yp(j)=0.5*(yc(j)+yc(j-1));
end
yp(nj)=yc(njm1);


% viscosity
viscos=1/30000;
%viscos=1/300;

% initialize ustar
ustar=1;
yplus=yp*ustar/viscos;
for j=2:nj
   u_reichard(j)=1/0.4*log(1+0.4*yplus(j))+7.8*(1-exp(-yplus(j)/11))-yplus(j)*exp(-yplus(j)/3);
end

% under-relaxation
urf=0.9;

% plot k for each iteration at node jmon
jmon=10; 

% turbulent constants 
sigma_k=0.8;
sigma_om=1.35;
betas=0.09;
kappa=0.42

gama0 = 0.42;
beta = 0.075;

small=1e-10;
great=1e10;

% init u, k & omega
k(2:nj)=1;
u(2:nj)=1;
om(1:nj)=10;
vist(1:nj)=100*viscos;

k(nj)=0;


% do a loop over all nodes (except the boundary nodes)
for j=2:nj-1

% compute dy_s
   dy_s=yp(j)-yp(j-1);

% compute dy_n
   dy_n=yp(j+1)-yp(j);

% compute deltay
   delta_y(j)=yc(j)-yc(j-1);
 
   dn(j)=1/dy_n;
   ds(j)=1/dy_s;

% interpolation factor
   del1=yc(j)-yp(j);
   del2=yp(j+1)-yc(j);
   fy(j)=del1/(del1+del2);

end


u(1)=0;


vist(1)=0;
vist(nj)=0;


% do max. niter iterations
for n=1:niter

   for j=2:nj-1

% compute turbulent viscosity
      vist_old=vist(j);
      vist(j)=urf*k(j)/om(j)+(1-urf)*vist_old;

    end

% boundary conditions for vist from wall functions
   j=1;
   yplus_2=yp(2)*ustar/viscos;
   vist(j)=ustar*kappa*yp(2)/log(9*yplus_2)-viscos;  %wall viscosity. Subtract 'viscos' since it's added in as below

   j=nj;
   vist(j)=vist(1);

% solve u
    for j=2:nj-1

% driving pressure gradient
      su(j)=delta_y(j);

      sp(j)=0;

% interpolate turbulent viscosity to faces
      vist_n=fy(j)*vist(j+1)+(1-fy(j))*vist(j);
      vist_s=fy(j-1)*vist(j)+(1-fy(j-1))*vist(j-1);

% compute an & as
      an(j)=(vist_n+viscos)*dn(j);
      as(j)=(vist_s+viscos)*ds(j);
    end

% boundary conditions for u
   u(1)=0;
   u(nj)=0;


   for j=2:nj-1
% compute ap
      ap(j)=an(j)+as(j)-sp(j);

% under-relaxate
      ap(j)= ap(j)/urf;
      su(j)= su(j)+(1-urf)*ap(j)*u(j);

% use Gauss-Seidel
      u(j)=(an(j)*u(j+1)+as(j)*u(j-1)+su(j))/ap(j);
   end

% use TDMA
%   u=tdma(ap,an,as,su,u,nj);


% monitor the development of tau_w
    tau_w(n)=ustar^2;
    vist_iter(n)=vist(jmon);


% print iteration info
    disp(['Iteration number: ',num2str(n)])

% check for convergence (when converged, the difference should be less then 1e-4)
    if abs(tau_w(n)-1) < 1.e-4
% do at least 1000 iter 
        if n > 1000
           disp(' ')
           disp('Converged!')
           break
        end
    end

% solve k
    dudy=gradient(u,yp);
    dudy2=dudy.^2;
    for j=2:nj-1

% production term
      su(j)=vist(j)*dudy2(j)*delta_y(j);

      sp(j)=-betas*om(j)*delta_y(j);

% compute an & as
      vist_n=fy(j)*vist(j+1)+(1-fy(j))*vist(j);
      an(j)=(vist_n/sigma_k+viscos)*dn(j);
      vist_s=fy(j-1)*vist(j)+(1-fy(j-1))*vist(j-1);
      as(j)=(vist_s/sigma_k+viscos)*ds(j);
    end

% boundary conditions for k from wall functions
% U+ = 1/kappa*log10(E*yplus) => ustar = kappa*U/log10(E*yplus)
   yplus_2=yp(2)*ustar/viscos;
   ustar=kappa*u(2)/log(9*yplus_2);
   j=2;
   k(j)=betas^(-0.5)*ustar^2;
   sp(j)=-great;
   su(j)=great*k(j);

   j=njm1;
   k(j)=betas^(-0.5)*ustar^2;
   sp(j)=-great;
   su(j)=great*k(j);


   for j=2:nj-1
% compute ap
      ap(j)=an(j)+as(j)-sp(j);

% under-relaxate
      ap(j)= ap(j)/urf;
      su(j)= su(j)+(1-urf)*ap(j)*k(j);

% use Gauss-Seidel
      k(j)=(an(j)*k(j+1)+as(j)*k(j-1)+su(j))/ap(j);
   end

% use TDmA
%  k=tdma(ap,an,as,su,k,nj);

% monitor the development of k in node jmon
   k_iter(n)=k(jmon);

%****** solve om-eq.

   for j=2:nj-1
% compute an & as
      vist_n=fy(j)*vist(j+1)+(1-fy(j))*vist(j);
      an(j)=(vist_n/sigma_om+viscos)*dn(j);
      vist_s=fy(j-1)*vist(j)+(1-fy(j-1))*vist(j-1);
      as(j)=(vist_s/sigma_om+viscos)*ds(j);

% production term
      su(j)=gama0*dudy2(j)*delta_y(j);

% dissipation term
      sp(j)=-beta*om(j)*delta_y(j);

    end

% boundary conditions for om from wall functions
   j=2;
   om(j)=ustar/kappa/yp(2)/betas^0.5;
   sp(j)=-great;
   su(j)=great*om(j);

   j=njm1;
   om(j)=ustar/kappa/yp(2)/betas^0.5; % same wall distant
   sp(j)=-great;
   su(j)=great*om(j);



   for j=2:nj-1
% compute ap
      ap(j)=an(j)+as(j)-sp(j);

% under-relaxate
      ap(j)= ap(j)/urf;
      su(j)= su(j)+(1-urf)*ap(j)*om(j);

% use Gauss-Seidel
      om(j)=(an(j)*om(j+1)+as(j)*om(j-1)+su(j))/ap(j);

   end

% use TDmA
%  om=tdma(ap,an,as,su,om,nj);

   om_iter(n)=om(jmon);

% flush the buffer (in octave)
% fflush(stdout)

end

% compute shear stress
uv=-vist.*dudy;
% set uv=1 at first cell
uv(2)=-1;


yplus=yp*ustar/viscos;
% plot u
figure(1)
tau_turb=viscos*u(2)/yp(2);
plot(u,yp,'linewidth',2)
hold
plot(u_reichard(2:nj/2),yp(2:nj/2),'r--','linewidth',2)
xlabel('u')
ylabel('y')
legend('PDH','Reichard')
print u_30000.ps -deps

% plot k
figure(2)
plot(k,yp,'linewidth',2)
hold
xlabel('k')
ylabel('y')
print k_30000.ps -deps

% plot tau_w versus iteration number
figure(3)
plot(tau_w,'linewidth',2)
title('wall shear stress')
xlabel('Iteration number')
ylabel('tauw')
print tauw.ps -deps

% plot k(jmon) versus iteration number
figure(4)
plot(k_iter,'linewidth',2)
title('k in node jmon')
xlabel('Iteration number')
ylabel('k')
print k_iter.ps -deps

% plot om(jmon) versus iteration number
figure(5)
plot(om_iter,'linewidth',2)
title('omega in node jmon')
xlabel('Iteration number')
ylabel('omega')
print om_iter.ps -deps

eps=0.09*k.*om;
% save data
ff(:,1)=yp;
ff(:,2)=u;
ff(:,3)=uv;
ff(:,4)=k;
ff(:,5)=eps;
ff(:,5)=vist;
ff(:,6)=yc;

save yp_u_uv_k_eps_vist_yc_wallf_30000.dat ff -ascii

