% 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 PDH k-omegaa model [1]. Re=u_tau*h/nu=8000 (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=30000;

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

% create the grid

nj=98;
njm1=nj-1;
yfac=1.15;
dy=0.1;
yc(1)=0;
for j=2:nj/2
    yc(j)=yc(j-1)+dy;
    dy=yfac*dy;
end

ymax=yc(nj/2);

% cell faces
for j=2:nj/2
   yc(j)=yc(j)/ymax;
   yc(nj+1-j)=2-yc(j-1);
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/8000;

% under-relaxation
urf=0.5;

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

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

gama0 = 0.42;
beta = 0.075;

cr1=0.75 

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

% init u, k & omega
k(2:nj)=1;
u(2:njm1)=1;
om(2:nj)=10;
vist(1: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);
      ret=k(j)/(viscos*om(j));
      ret=max(ret,1.e-5);
      fmu=0.025+(1-exp(-(ret/10.0)^0.75))*(0.975+(1.e-3/ret)*exp(-(ret/200.)^2));
      vist(j)=urf*fmu*k(j)/om(j)+(1-urf)*vist_old;

    end

% 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 TDrMmA
%   u=tdma(ap,an,as,su,u,nj);


% monitor the development of u_tau in node jmon
    tau_w(n)=viscos*u(2)/yp(2);


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

% check for convergence (when converged, the wall shear stress must be one)
    if abs(tau_w(n)-1) < 0.001
% 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);

% dissipation term
      ret=k(j)/(viscos*om(j));
      ret=max(ret,1.e-5);

      fk=1-0.722*exp(-(ret/10.)^4);
      sp(j)=-fk*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
   k(1)=0;
   k(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)*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.

   dkdy=gradient(k,yp);
   domdy=gradient(om,yp);
   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);

      ret=k(j)/(viscos*om(j));
      ret=max(ret,1.e-5);

      fmu=0.025+(1-exp(-(ret/10.0)^0.75))*(0.975+(1.e-3/ret)*exp(-(ret/200.)^2));
      fomega=1+4.3*exp(-(ret/1.5)^0.5);
      gama = gama0*fomega;

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

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

% cross diffusion term
      crosv=dkdy(j)*domdy(j);

      crost=cr1*vist(j)*crosv/k(j)*delta_y(j);

      su(j)=su(j)+max(crost,0);
      sp(j)=sp(j)+min(crost,0)/om(j);

    end


% b.c. south wall
      dy=yp(2);
      omega=6.*viscos/0.075/dy^2;
      sp(2)=-great;
      su(2)=great*omega;

% b.c. north wall
      dy=yp(nj)-yp(njm1);
      omega=6.*viscos/0.075/dy^2;
      sp(njm1)=-great;
      su(njm1)=great*omega;

   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;


% plot u
%load ~/calc_channel/PDH-re8000/u_8000_PDH.dat
figure(1)
tau_turb=viscos*u(2)/yp(2);
plot(u,yp,'linewidth',2)
hold
plot(u_8000_PDH(:,2),yp(2:njm1),'o')
xlabel('u')
ylabel('y')
legend('PDH','Matlab')
print u_8000.ps -deps

% plot k
%load ~/calc_channel/PDH-re8000/k_8000_PDH.dat
figure(2)
plot(k,yp,'linewidth',2)
hold
plot(k_8000_PDH,yp(2:njm1),'o')
xlabel('k')
ylabel('y')
legend('PDH','Matlab')
print k_8000.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

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

%save yp_u_k_om_vist_uv_yc_PDH_8000.dat ff -ASCII
save yp_u_k_om_vist_uv_yc_PDH_8000.dat ff -ascii