% Initialize_vortex
% Program to generate initial velocity field for vortex case

clear
clc
close all

Rc   = 0.005;               % Core radius of vortex (m)
Umax = 0.5;                 % Vortex maximum velocity (m/s)
Ucof = 0.1;                 % Coflow velocity (m/s)
L    = 0.1;                 % Length of domain (m)

Res  = [2 5 10 20 40];      % Number of grid cells across vortex diameter

for n = 1:numel(Res)
    
    dx  = (2*Rc)/Res(n);    % Grid spacing x (m)
    xv = 0;                 % Center of vortex, x (m)
    zv = 0;                 % Center of vortex, z (m)
    
    % Staggered grid format
    %
    %  *--NF--*
    %  |      |        z
    %  |  CC  EF       ^
    %  |      |        |
    %  *------*        *--> x
    %
    %  CC: cell center (xc(i),zc(k))
    %  EF: center of east face (x(i),zc(k))
    %  NF: center of north face (xc(i),z(k))
    %
    %  Ux-velocity calculated at EF
    %  Uz-velocity calculated at NF
    
    % Staggered grid
        
    Nx = L/dx;
    for i=1:Nx+1
        x(i)  = -L/2 + dx*(i-1);
        xc(i) = -L/2 + dx*(i-1) - dx/2 ;
    end
    dz = dx;
    Nz = L/dz;
    for k=1:Nz+1
        z(k)  = -L/2 + dz*(k-1);
        zc(k) = -L/2 + dz*(k-1) - dz/2;
    end
    Ny    = 1;
    y(1)  = -0.001;
    y(Ny) =  0.001;
    
    % Velocity field
    
    Ux = NaN(Nx+1,Ny+1,Nz+1);
    Uy = NaN(Nx+1,Ny+1,Nz+1);
    Uz = NaN(Nx+1,Ny+1,Nz+1);
    
    for i = 1:Nx+1
        for j = 1:Ny+1
            for k = 1:Nz+1
                Ux(i,j,k) = (Umax*((zc(k)-zv)/Rc)*exp((1/2)-...
                    ((((x(i)-xv)^2)+((zc(k)-zv)^2))/(2*(Rc^2)))))+Ucof;
                Uy(i,j,k) = 0;
                Uz(i,j,k) = -Umax*((xc(i)-xv)/Rc)*exp((1/2)-...
                    ((((xc(i)-xv)^2)+((z(k)-zv)^2))/(2*(Rc^2))));
            end
        end
    end
    
    % Store velocity field in CSV file
    clear uvw
    
    uvw(1,1) = 0;
    uvw(1,2) = Nx;
    uvw(1,3) = 0;
    uvw(1,4) = Ny;
    uvw(1,5) = 0;
    uvw(1,6) = Nz;
    
    uvw((2:(1+(Nx+1)*(Ny+1)*(Nz+1))),1) = Ux(:);
    uvw((2:(1+(Nx+1)*(Ny+1)*(Nz+1))),2) = Uy(:);
    uvw((2:(1+(Nx+1)*(Ny+1)*(Nz+1))),3) = Uz(:);
    
    % Write FDS uvw.csv input file
    FID = sprintf('init_vortex_%s_uvw.csv',num2str(Res(n)));
    csvwrite(FID,uvw);
    
end
