% This file can be downloaded at

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

% choose "Synthetic Inlet Boundary Conditions" on the left side
%
%
% 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.


% 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


nmodes =150;
wew1fct=2;
visc=15.2e-6;
qm=4;
delta=0.5;        % we assume that we have a boundary layer thickness of 0.5
sli=0.1*delta;    % make the length scaler one tenth of delta
up=sqrt(2*qm/3);
epsm=qm^1.5/sli;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     number of  grid points in y, z

nj=2; % number of grid point in y direction
nk=63;  % number of grid point in z direction

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

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

% generate y-grid
yp(1)=0;
yp(2)=1;

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

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


% generate nstep realizations
nstep=1;

% the Reynols stress tensor it taken DNS
% % Authors: Moser, Kim & Mansour
% Reference: DNS of Turbulent Channel Flow up to Re_tau=590, 1999,
%            Physics of Fluids, vol 11, 943-945.
%
%   y           y+        R_uu      R_vv       R_ww       R_uv        R_uw       R_vw
%
%2.7060e-02 1.5890e+01 7.6684e+00 3.1974e-01 1.4997e+00 -6.6206e-01 1.7267e-03 8.7423e-04

stress=[7.6684 -6.6206e-01 0; ...
-6.6206e-01 3.1974e-01 0; ...
0 0 1.4997];



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);

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

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

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

% 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);


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;


%a11=2.0781;
%a22=0.32418;
%a33=0.59775;
%
%r11=-0.99552;
%r12=-0.094591;
%r13=0;
%r21=0.094591;
%r22=-0.99552;
%r23=0;
%r31=0;
%r32=0;
%r33=1;

%a11=1;
%a22=1;
%a33=1;
%
%r11=1;
%r12=0;
%r13=0;
%
%r21=0;
%r22=1;
%r23=0;
%
%r31=0;
%r32=0;
%r33=1;
%

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)])

%
% 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);


% 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;

% store the fluctuations for each time step it (we have only one cell in y-dirrection; j=1)
% these will be used for plotting and evaluation below
      ut(k,it)=u(j,k);
      vt(k,it)=v(j,k);
      wt(k,it)=w(j,k);
   end % end j=1:nj-1 
   end % end k=1:nk-1


end % end it=1:nstep


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

%
%=========================== chapter 4 ============================================
%

umean=mean(mean(ut));
vmean=mean(mean(vt));
wmean=mean(mean(wt));

disp(['umean=',num2str(umean)])
disp(['vmean=',num2str(vmean)])
disp(['wmean=',num2str(wmean)])

disp('program stop');

break


%**************** time history  *************************************
%% plot u of one realization (time step 5)
%
figure(1)
hold off
plot(zp(1:nk-1),ut(:,5),'linewidth',2)
xlabel('z','fontsize',[20])
ylabel('u','fontsize',[20])
hold
%axis([0 0.2 -0.2 1])
%
handle=gca;
set(handle,'fontsize',[20])
print u_inst.ps -depsc

