% This file can be downloaded at

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

% choose "Synthetic Inlet Boundary Conditions" on the left side
%
%
% This Matlab script generates synthetic fluctuations used in 

% L. Davidson, "Two-equation hybrid RANS-LES Models: A novel way to treat k and omega
%  at the inlet", Turbulence, Heat and Mass Transfer, THMT-15, Sarajevo, Bosnia and Herzegovina,
% 2015.
%
% L. Davidson, "Zonal PANS: evaluation of different treatments of the {RANS-LES} interface (to appear)",
% Journal of Turbulence, 2015.
%
%
% Litterature
%
% L. Davidson, "Using Isotropic Synthetic Fluctuations as Inlet Boundary Conditions for Unsteady
% Simulations", "Advances and Applications in Fluid Mechanics",
% Vol 1, No  =1, pp. 1-35, 2007
%
% L. Davidson, "HYBRID LES-RANS:  Inlet Boundary Conditions for Flows With Recirculation",
% "Advances in Hybrid RANS-LES Modelling",
% Notes on Numerical Fluid Mechanics and Multidisciplinary Design,
% Springer Verlag, pp. 55-66, Vol. 97, 2008.
%
% L. Davidson, "Hybrid LES-RANS: Inlet Boundary Conditions for Flows Including Recirculation",
% 5th International Symposium on Turbulence and Shear Flow Phenomena,
% Vol. 2, pp. 689-694, Munich, Germany, 2007.
%
% M. Billson and L.-E. Eriksson and L. Davidson, "Modeling of Synthetic Anisotropic 
% Turbulence and its Sound Emission", The 10th AIAA/CEAS Aeroacoustics Conference, 
% AIAA 2004-2857, Manchester, United Kindom, 2004.
%
% M. Billson, "Computational Techniques for Turbulence Generated Noise",
% Dept. of Thermo and Fluid Dynamics, Chalmers University of Technology",
% Göteborg, Sweden, 2004.
%
% L. Davidson and S.-H. Peng, "Embedded Large-Eddy Simulation Using the Partially
% Averaged Navier-Stokes Model", AIAA J, vol. 51, number=5, pp=1066-1079,
%  doi="doi: 10.2514/1.J051864", 2013.
%
% L. Davidson, "Two-equation hybrid RANS-LES models: A novel way to treat k and omega 
% at inlets and at embedded interfaces", J. of Turbulence, Volume 18, Issue 4, pp. 291-315, 2017.
%

% In this program we generate fluctuation along one z line only


close all
clear all
amp=1.452762113;

%=========================== chapter 1 ============================================

%!!!  number of modes                        = nmodes
%!!!  smallest wavenumber                    = dxmin
%!!!  ratio  of ke and kmin (in wavenumber)  = wew1fct
%!!!  turb. velocity scale                   = up
%!!!  turb. kin. energy                      = qm
%!!!  diss. rate.                            = epsm
%!!!  kinetic viscosity                      = visc
%!!!  length scale                           = sli
%!!!  time step                              =dt


nmodes =150;
wew1fct=2;
visc=1/8000;
sli=0.3;    % make the length scaler 0.3 of the half channel width
up=1;
epsm=up^3/sli;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     number of  grid points in y, z

nj=97; % number of grid point in y direction
nk=33;  % number of grid point in z direction

% generate x-grid
xp(1)=0;
xp(2)=1;

% generate z-grid
dz=1.6/(nk-1);
zp(1)=0;
for k=2:nk
   zp(k)=zp(k-1)+dz;
end

% load RANS data created by rans.m (which can be downloaded)
load yp_u_k_om_vist_uv_yc_PDH_8000.dat

% viscosity
nu=visc;

y=yp_u_k_om_vist_uv_yc_PDH_8000(:,1);
u=yp_u_k_om_vist_uv_yc_PDH_8000(:,2);
k=yp_u_k_om_vist_uv_yc_PDH_8000(:,3);
om=yp_u_k_om_vist_uv_yc_PDH_8000(:,4);
uv_rans=yp_u_k_om_vist_uv_yc_PDH_8000(:,6);
yc=yp_u_k_om_vist_uv_yc_PDH_8000(:,7);
ed=0.09*k.*om;

yp=yc; % here yp means face coordinates


% search min grid step
dminy=min(diff(yp(1:nj)));
dminz=min(diff(zp));
dxmin=min(dminy,dminz);

% don't let the minimum distance be too small. This will create unnecessay high wavenumbers
dxmin=max(dxmin,1.7e-3);

% create  seed from time (must be negative)
v=clock;
iseed=-(v(4)+v(5)+v(6));

% set iseed to a constant value; this gice the same synt. fluct every time you run the code
iseed=-5;

% generate nstep realizations
nstep=1;

for i=1:nstep
   dt(i)=1e-3;
end

ustar=(nu*u(2)/y(2))^0.5;

dudy=gradient(u,y);
nj_rans=length(k);

dudy=gradient(u,y);
nj_rans=length(k);

% set synth. fluct at old time step to zero
uprim_old=zeros(nj,nk);
vprim_old=zeros(nj,nk);
wprim_old=zeros(nj,nk);

% lowe Re EARSM by
% S. Wallin and A. V. Johansson,
%"A New Explicit Algebraic Reynolds Stress Model for Incompressible
%and Compressible Turbulent Flows", JFM,
%Vol. 403, pp. 89-132, 2000

for j=1:nj_rans
   ttau(j)=max(k(j)/ed(j),6*(nu/ed(j))^0.5);
   dist=min(y(j),y(nj_rans)-y(j));
   yplus(j)=dist*ustar/nu;
   om11=0;
   om12=ttau(j)*0.5*dudy(j);
   om21=-om12;
   om22=0;
   s11=0;
   s12=ttau(j)*0.5*dudy(j);
   s21=s12;
   s22=0;
   s33=0;
   vor=(-2*om12^2);
   str=(s11.^2+s12.^2+s21.^2+s22.^2);

   cpr1=9/4*(1.8-1.);
   p1=(1/27*cpr1^2+9/20*str-2/3*vor)*cpr1;
   p2=p1^2-(1/9*cpr1^2+9/10*str+2/3*vor)^.3;

   if p2 > 0 
      if p1-p2^0.5 >= 0
         sigg=1;
      else
         sigg=-1;
      end 
      un=cpr1/3+(p1+p2^0.5)^(1/3)+sigg*(abs(p1-p2^0.5))^(1/3);
   else
      un=cpr1/3+2*(p1^2-p2)^(1/6)*cos(1/3*acos(p1/(sqrt(p1^2-p2))));
   end 

   beta1=-6./5.*un/(un^2-2.*vor);
   beta4=beta1/un;


% low re
   aplus=26;
   f1(j)=1-exp(-yplus(j)/aplus);

   beta1_high(j)=beta1;

   beta1=f1(j)*beta1;
   b2=1.8;
   sigma2=str/2; % see Wallin&Johansson 2000, ekv. 3.2
   beta2=(3*b2-4)/2/sigma2*(1-f1(j)^2);
   beta4=f1(j)^2*beta4-(1-f1(j)^2)*b2/4/sigma2;

   term_lowre=beta2*(s11^2+s12^2-str/3);
   uu(j)=2/3*k(j)+k(j)*(beta1*s11+beta4*(s12*om21-om12*s21)+term_lowre);

   vv(j)=2/3*k(j)+k(j)*(beta1*s22+beta4*(s21*om12-om21*s12)+term_lowre);

   ww(j)=2/3*k(j);

   uv(j)=k(j)*(beta1*s21+beta4*(s11*om12-om12*s22));


end



% max(uv) at j=19
j=19;

stress=[uu(j) uv(j) 0; uv(j) vv(j) 0; 0 0 ww(j)];
diag_sum=trace(stress)/3
stress=stress/diag_sum

[V,lambda] = eig(stress)


v1=[V(1,1);V(2,1); V(3,1)];
v2=[V(1,2);V(2,2); V(3,2)];
v3=[V(1,3);V(2,3); V(3,3)];

lambda_1=lambda(1,1);
lambda_2=lambda(2,2);
lambda_3=lambda(3,3);

lam(1)=lambda_1;
lam(2)=lambda_2;
lam(3)=lambda_3;


% lam_1 < lam_2 < lam_3

% switch so that largest eigenvalue in x1*, smallest in x2*

v1_new=v3;
v2_new=v1;
v3_new=v2;

% 1st eigenvector in 1st or 3rd quadrant
% 2nd eigenvector in 2nd or 4th quadrant

% switch sign on 12 and 21 to fix the above requirements
v1_new(2)=-v1_new(2);
v2_new(1)=-v2_new(1);


lambda_1_new=lambda_3;
lambda_2_new=lambda_1;
lambda_3_new=lambda_2;

r11=v1_new(1);
r12=v1_new(2);
r13=v1_new(3);

r21=v2_new(1);
r22=v2_new(2);
r23=v2_new(3);

r31=v3_new(1);
r32=v3_new(2);
r33=v3_new(3);

a11=lambda_1_new;
a22=lambda_2_new;
a33=lambda_3_new;


%

disp(['eigenvector 1= ',num2str(r11),', ', num2str(r21),', ',num2str(r31)])
disp(['eigenvector 2= ',num2str(r12),', ', num2str(r22),', ',num2str(r32)])
disp(['eigenvector 3= ',num2str(r13),', ', num2str(r23),', ',num2str(r33)])

disp(' ')

disp(['eigenvalue 1= ',num2str(a11)])
disp(['eigenvalue 2= ',num2str(a22)])
disp(['eigenvalue 3= ',num2str(a33)])



umean=u;
% clear variables
clear u v 

%
% zero all arrays to zero

wnr=zeros(nmodes,1);
fi=zeros(nmodes,1);
teta=zeros(nmodes,1);
psi=zeros(nmodes,1);

iy=0;
iv=zeros(32,1);
%
%=========================== chapter 2 ============================================
%
disp(['No of time steps: ',num2str(nstep)])

for it=1:nstep

   disp(['time step no: ',num2str(it)])

% compute random angles
   [fi,psi,alfa,teta,iy,iv] = angles(nmodes,iseed,iy,iv);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%     highest wave number
   wnrn=2*pi/dxmin;
%
%     k_e (related to peak energy wave number)
   wnre=9*pi*amp/(55*sli);
%
% wavenumber used in the viscous expression (high wavenumbers) in the von Karman spectrum
   wnreta=(epsm/visc^3)^0.25;

%     smallest wavenumber 
   wnr1=wnre/wew1fct;

% wavenumber step
   dkl=(wnrn-wnr1)/nmodes;

% wavenumber at faces
   for m=1:nmodes+1
       wnrf(m)=wnr1+dkl*(m-1);
   end

% wavenumber at cell centers
   for m=1:nmodes
       wnr(m)=.5*(wnrf(m+1)+wnrf(m));
       dkn(m)=wnrf(m+1)-wnrf(m);
   end

% invert the eigenvalue matrix (anisotropic)
   a11i=1/a11;
   a22i=1/a22;
   a33i=1/a33;

%   wavenumber vector from random angles
   for m=1:nmodes
       kxio(m)=sin(teta(m))*cos(fi(m));
       kyio(m)=sin(teta(m))*sin(fi(m));
       kzio(m)=cos(teta(m));
%
% sigma (s=sigma) from random angles. sigma is the unit direction which gives the direction
% of the synthetic velocity vector (u, v, w)
       sxio(m)=cos(fi(m))*cos(teta(m))*cos(alfa(m))-sin(fi(m))*sin(alfa(m));
       syio(m)=sin(fi(m))*cos(teta(m))*cos(alfa(m))+cos(fi(m))*sin(alfa(m));
       szio(m)=-sin(teta(m))*cos(alfa(m));
   
   end 


%=========================== chapter 3 ============================================
%
%     loop through grid 
   for k=1:nk-1
   for j=1:nj-1

%     cell center coordinates
       xc=0.5;
       yc=0.5*(yp(j)+yp(j+1));
       zc=0.5*(zp(k)+zp(k+1));

       xcp=r11*xc+r21*yc+r31*zc;
       ycp=r12*xc+r22*yc+r32*zc;
       zcp=r13*xc+r23*yc+r33*zc;


%     initiate turbulent velocities to zero.
       utrp=0;
       vtrp=0;
       wtrp=0;
   
% loop over all wavenumbers
       for m=1:nmodes
   

% r21, r31, r12, r32, r13, r23 = 0 in isotropic turbulence
          kxi=r11*kxio(m)+r21*kyio(m)+r31*kzio(m);
          kyi=r12*kxio(m)+r22*kyio(m)+r32*kzio(m);
          kzi=r13*kxio(m)+r23*kyio(m)+r33*kzio(m);

          sxi=r11*sxio(m)+r21*syio(m)+r31*szio(m);
          syi=r12*sxio(m)+r22*syio(m)+r32*szio(m);
          szi=r13*sxio(m)+r23*syio(m)+r33*szio(m);

% a11, a22, a33=1 in isotropic turbulence
          sx=sqrt(a11)*sxi;
          sy=sqrt(a22)*syi;
          sz=sqrt(a33)*szi;


          kx=sqrt(a11i)*kxi*wnr(m);
          ky=sqrt(a22i)*kyi*wnr(m);
          kz=sqrt(a33i)*kzi*wnr(m);
          rk=sqrt(kx^2+ky^2+kz^2);

          sxm(m)=sx;
          sxim(m)=sxi;
          a11m(m)=a11;
          a11im(m)=a11i;

          kzm(m)=kz;
          kzim(m)=kzi;
          a33im(m)=a33i;
    
          rkm(m)=rk;

% if the wavenumber, rk, is smaller than the largest wavenumber, then create fluctuations
          if rk < wnrn

             arg=kx*xcp+ky*ycp+kz*zcp+psi(m);

             tfunk=cos(arg);


% von Karman spectrum
             e=amp/wnre*(wnr(m)/wnre)^4/((1+(wnr(m)/wnre)^2)^(17/6))*exp(-2*(wnr(m)/wnreta)^2);

             utn=sqrt(e*up^2*dkn(m));

% synthetic velocity field 
             utrp=utrp+2*utn*tfunk*sx;
             vtrp=vtrp+2*utn*tfunk*sy;
             wtrp=wtrp+2*utn*tfunk*sz;

          end  % end if rk < wnrn

      end  % end for m=1:nmodes


       

      utrpn=r11*utrp+r12*vtrp+r13*wtrp;
      vtrpn=r21*utrp+r22*vtrp+r23*wtrp;
      wtrpn=r31*utrp+r32*vtrp+r33*wtrp;

% store the fluctuations u, v, w in 2D arrays 
      u(j,k)=utrpn;
      v(j,k)=vtrpn;
      w(j,k)=wtrpn;

   end % end j=1:nj-1 
   end % end k=1:nk-1




   disp(' ');
   disp('synthetic fluctuations created');
   disp(' ');

% Now scale the fluctuation using the RANS uv profile

% the inlet fluct should be scaled so that the max shear stress of the synthetic
% fluctuatons agree with the RANS value. Averaging the synthetic fluct over many
% realizations gives
   uv_synt_max=0.375;

% max uv_rans
   uv_rans_max=max(uv_rans);

% scaling factor
   ss_u=(uv_rans_max/uv_synt_max)^0.5;

% turbulent length scale
   rl=sli;

% turbulent time scale from Taylor's hypothesis using bulk
% velocity ubulk

   ubulk=trapz(yp,umean)/y(end);
   tturb=rl/ubulk;

   a=exp(-dt(it)/tturb);
   b=(1-a^2)^0.5;

%     loop through grid 
   for k=1:nk-1
   for j=1:nj-1

       ruv=uv_rans(j);
       ss_k=ss_u*(ruv/uv_rans_max)^0.5;

       uprim=u(j,k)*ss_k;
       vprim=v(j,k)*ss_k;
       wprim=w(j,k)*ss_k;

% assume that we have the synthetic fluct from the previous time step (uprim_old, ...)
       unew=a*uprim_old(j,k)+b*uprim;
       vnew=a*vprim_old(j,k)+b*vprim;
       wnew=a*wprim_old(j,k)+b*wprim;

% add unew, vnew, wnew to mean inlet profiles


% store the fluct to be used in the above expressionm at the next time step
       uprim_old(j,k)=unew;
       vprim_old(j,k)=vnew;
       wprim_old(j,k)=wnew;

   end 
   end

% end of time step

end % end it=1:nstep
% 

