%%  ESS 590: HEAT AND MASS FLOW MODELLING  SPRING 2008
%%
%%  Finite-Volume solution for Steady-State phi(x,z)
%%  diffusion with advection
%%
%%   EXPLORE FALSE DIFFUSION
%%
%%     div( rho u phi ) = div( Gamma grad phi) + S 
%%
%% For example,
%%  phi   = temperature                         (deg)
%%  Gamma = k/c
%%             k = thermal conductivity         (W m^-1 deg^-1)
%%             c = specific heat capacity       (J kg^-1 deg^-1)
%%  S     = sources                             (W m^-3)
%%        + S_C + S_P phi
%%
%%  This code uses PATANKAR POWER-LAW SCHEME for advective fluxes
%%
%%  Externally specified velocity field (u,w) must satisfy 
%%   mass conservation
%%
%%    div( rho u) = 0
%%
%%  Boundary Conditions:
%%   All boundaries are finite-volume interfaces
%%   (Patankar Practice B, page 69)
%%   However, there is an additional "point" X_B located on boundary.
%%   Its value and coefficients are determined by phi value or gradient
%%   boundary conditions there.
%%   Adjacent interior point is X_I
%%   (X_I - X_B) = dX_B
%%
%%    phi_B = phi_0
%%        a_B = 1,  all other a_J = 0
%%        b   = phi_0
%%
%%    d phi/dx = dphi_dx_0
%%        use gradient between X_B and X_I (interior point)
%%        d phi/dx|B = Gamma_I (phi_I - phi_B)/dx_B
%%    or
%%        phi_B = dphi_dx_0 * dx_B + phi_I
%%
%%  For each boundary there is a matrix bc_x
%%     with either 2 rows (Up,Down) or 2 columns (East,West).
%%  e.g. East  
%%   bc_e is [nz_P x 2] where nz_P is number of nodes on boundary 
%%     - first column contains value of BC at that node on boundary
%%     - second column specifies type of BC
%%            - 1 -> specified phi
%%            - 2 -> specified phi gradient across boundary
%%
%%   Similarly bc_w, bc_u, bc_d on other boundaries
%%
%%   fvm_SS_adv_diff_2D.m will correctly apply any combination of BC
%%        that are passed to it in matrices bc_e etc
%%
%%
%% Material Properties
%%  Gamma(x,z) and rho((x,z) can vary with position 
%%
%%
%%  In the spirit of matlab's vector capabilities, the program
%%   avoids loops over spatial indices.
%%  Instead, matrices are formed for all spatially varying parameters
%%   and variables.
%%
%%                                   Ed Waddington      May  7, 2006
%%                                      updated         May 22, 2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
   global c_x  c_z  c_prime    %  velocity parameterization
%
%
% False Diffusion Test
% bring in 2 streams at different phi on boundary
% When Gamma is very small, the 2 streams should stay distinct
% ------------------------------------------------------------
  % Gamma_j = [ 0.1 0.1 0.01 0.001*[1 1 1 1] ];
   Gamma_j = [ 0.1 0.1 0.01 0.0001*[1 1 1 1] ];
%
% set up numbers of volumes in each direction
% for each of several runs
%   decreasing volume size should decrease false diffusion
% --------------------------------------------------------
  % N_vol = [ 9 9 9 9 9 25 100 ];
   N_vol = [ 10 10 10 10 10 50 1000 ];
%
% set up velocity patterns
%   1,2 velocity parallel to grid
%           same behavior with u=0 or w=0
%   3,4 velocity parallel to grid (w=0)
%           decreasing Gamma
%   5,6,7 velocity at 45 degrees to grid, decreasing dx,dz
% -------------------------------------------------------
   c_x_vec = [ 0 sqrt(2)*[1 1 1] 1 1 1 ]; 
   c_z_vec = [ sqrt(2) 0*[1 1 1] 1 1 1 ]; 
%


% set values of phi in 2 incoming streams
     bc_phi_0 = 50;                  %  deg C
     bc_phi_1 = 100;                 %  deg C


   for i_count = 1:length(N_vol)  %  loop over different models
      pause
%
% set parameters for velocity calculations
%   u = c_x + c_prime * x
%   w = c_z - c_prime * z
% --------------------------
      c_x = c_x_vec(i_count);
      c_z = c_z_vec(i_count);
      c_prime = 0;
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%  set up finite volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This code sets up uniform-sized volumes (could be generalized)
%   There are points on all boundaries, which are also
%   volume interfaces  (Patankar Practice B, page 69)
%
      Lx    = 10;               %  x-extent of domain
      nx_fv = N_vol(i_count);   %  number of finite volumes in x
      nx_P  = nx_fv + 2;        %  number of nodes in x
%
%
      Lz    = 10;               %  z-extent of domain
      nz_fv = N_vol(i_count);   %  number of finite volumes in z
      nz_P  = nz_fv + 2;        %  number of nodes in z
%
%
%  set up vectors x_edges_vec, z_edges_vec - 
%  defining (x,z) coordinates of finite-volume edges
%
%    set x-widths of volumes 
%    (this bit of code sets up uniform edge spacing)
%   -------------------------------------------------
      dx = Lx/nx_fv;              %  z-widths of volumes
      x_edges_vec = dx * (0:nx_fv);
%
  %%   X_edges = repmat(x_edges_vec, nz_fv+1, 1);
%
%
%    set x-widths of volumes 
%    (this bit of code sets up uniform edge spacing)
%   -------------------------------------------------
      dz = Lz/nz_fv;               %  z-widths of volumes
      z_edges_vec = dz * (0:nz_fv)';
%
  %%   Z_edges = repmat(z_edges_vec, 1, nx_fv+1);

%
% set nodes at finite-volume centers and on all boundaries
      x_P_vec = [ x_edges_vec(1) ...
                 (x_edges_vec(1:end-1) + diff(x_edges_vec)/2) ...
                  x_edges_vec(end) ]; 
%
%
      z_P_vec = [ z_edges_vec(1) ...
                 (z_edges_vec(1:end-1) + diff(z_edges_vec)/2)' ...
                  z_edges_vec(end) ]'; 
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Boundary Conditions 
%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%  Boundary conditions along each boundary are carried in an array
%  labelled bc_x.  For each boundary point, 
%      - first entry is numerical value of BC,
%      - second entry indicates type of BC
%               1 for specified value
%               2 for specified gradient
%  ---------------------------------------------------
%

   if(i_count == 1)
%

%  Flow in z only
%  ----------------
%%  phi gradients at east boundary
     bc_e_0  = 0;                %   deg m^-1
     bc_type = 2;                %  gradient BC
     bc_e = repmat( [ bc_e_0  bc_type ], nz_P, 1 );
%
%
%%  phi gradients at west boundary
     bc_w_0  = 0;                %   deg m^-1
     bc_type = 2;                %  gradient BC
     bc_w = repmat( [ bc_w_0  bc_type ], nz_P, 1 );
%
%
%%  temperature phi_u on upper boundary
%   west half cold, east half hot
     n_cold = fix(nx_P/2) + 1;
     n_hot = nx_P - n_cold;
%
     bc_type = 1;                  %  specified-value BC
  %%   bc_u   = repmat( [ bc_phi_0  bc_type ]', 1, nx_P );
%
       bc_u   = [ repmat( [ bc_phi_0  bc_type ]', 1, n_cold ) ...
                  repmat( [ bc_phi_1  bc_type ]', 1, n_hot ) ];
%
%
%%  phi gradient on lower boundary
       bc_d_0 = 0;                  %  deg C
       bc_type = 2;                  %  specified-value BC
       bc_d   = repmat( [ bc_d_0  bc_type ]', 1, nx_P );
%
%
   elseif( (i_count > 1) & (i_count < 5) )
%
%
%  Flow in x only
%  ----------------
%  phi gradient at east boundary
       bc_e_0  = 0;                  %   deg m^-1
       bc_type = 2;                 %  gradient BC
%
       bc_e = repmat( [ bc_e_0  bc_type ], nz_P, 1 );
%
%%  phi at west boundary
       n_cold = fix(nx_P/2) + 1;
       n_hot = nx_P - n_cold;
%
       bc_type   = 1;                 %  fixed phi BC
%
       bc_w   = [ repmat( [ bc_phi_0  bc_type ]', 1, n_cold ) ...
                  repmat( [ bc_phi_1  bc_type ]', 1, n_hot ) ]';
%
%
%%  phi gradient on upper boundary
       bc_u_0 = 0;                   %  deg C m^-1
       bc_type = 2;                  %  specified-flux BC
       bc_u   = repmat( [ bc_u_0  bc_type ]', 1, nx_P );
%
%
%%  phi gradient on lower boundary
     bc_d_0 = 0;                   %  deg C m^-1
     bc_type = 2;                  %  specified-flux BC
     bc_d   = repmat( [ bc_d_0  bc_type ]', 1, nx_P );
%
%
   else
%
%
%  flow is at 45 degres on 2-D domain
% ------------------------------------
%  phi gradient at east boundary
     bc_e_0 = 0;                  %   deg m^-1
     bc_type   = 2;                 %  gradient BC
     bc_e = repmat( [ bc_e_0  bc_type ], nz_P, 1 );
%
%
%%  phi at west boundary
     bc_w_0 = bc_phi_1;                  %   deg m^-1
     bc_type   = 1;                 %  gradient BC
     bc_w = repmat( [ bc_w_0  bc_type ], nz_P, 1 );
%
%
%%  phi on upper boundary
     bc_u_0 = bc_phi_0;                   %  deg C
     bc_type = 1;                  %  specified-value BC
     bc_u   = repmat( [ bc_u_0  bc_type ]', 1, nx_P );
     
  %%   bc_u_1 = 100;
  %%   bc_u(1, 1:fix(nx_P/2)) =  bc_u_1;
     
%
%
%%  phi gradient on lower boundary
     bc_d_0 = 0;                   %  deg C
     bc_type = 2;                  %  specified-value BC
     bc_d   = repmat( [ bc_d_0  bc_type ]', 1, nx_P );
%
   end   %     if(i_count == 1)  ...
%
%
%%  material parameters
%%%%%%%%%%%%%%%%%%%%%%%
% set up thermal-property matrix - uniform conductivity
     Gamma_0 = Gamma_j(i_count);
     Gamma_P = repmat( Gamma_0, nz_P, nx_P );
%
%
% set up density matrix
     rho_0 = 1;                       %   kg m^-3
     rho = repmat( rho_0, nz_P, nx_P );
%

%
%%  source term   
%%%%%%%%%%%%%%%%%%%%%
%   S = S_C + S_P * phi_P 
%   ----------------------
     S_C_0  = 0;                   %  W m^-3 s^-1
     S_P_0  = 0;                   %  W m^-3 s^-1 deg^-1
%
     S_C  = repmat( S_C_0, nz_P, nx_P );  %  constant source term
     S_P  = repmat( S_P_0, nz_P, nx_P );  %  phi_P-dependent source
%

%
% This is it! solve for phi using Finite Volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    phi = fvm_SS_adv_diff_2D( x_edges_vec, z_edges_vec, ...
                       x_P_vec, z_P_vec, Gamma_P, rho, S_C, S_P, ...
                       bc_u, bc_d, bc_e, bc_w );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%
%%  Plot results
%%%%%%%%%%%%%%%%%
%
     figure(i_count)
%
%  set up colors for plotting
    red     = [ 1 0 0 ];
    green   = [ 0 1 0 ];
    blue    = [ 0 0 1 ];
    magenta = [ 1 0 1 ];
    black   = [ 0 0 0 ];
    color = [ red' blue' magenta' green' black'];
%
%
%  set up marker symbols for plotting
    marker = [ 'o' 'x' '*' '+' 'none ' ]; 
%
%

     X_P = repmat(x_P_vec, nz_P, 1);
     Z_P = repmat(z_P_vec, 1, nx_P);

%
       
%  2-D, plot phi(x,z) solution as perspective surface
%  --------------------------------------------------
      %%  surf(X_P, Z_P, phi )
         i_skip = 1;
         if(i_count == length(N_vol) )
            surf(X_P(1:i_skip:end, 1:i_skip:end), ...
                Z_P(1:i_skip:end, 1:i_skip:end), ...
                phi(1:i_skip:end, 1:i_skip:end), ...
                'linestyle', 'none' )
         else
         surf(X_P(1:i_skip:end, 1:i_skip:end), ...
             Z_P(1:i_skip:end, 1:i_skip:end), ...
             phi(1:i_skip:end, 1:i_skip:end) )
         end    
     %%    shading interp
         colorbar
         hold on; grid on; box on
         set(gca,'zlim', [ bc_phi_0, bc_phi_1 ] )
         xlabel('Distance  (m)', 'FontSize', 12)
         ylabel('Depth  (m)', 'FontSize', 12)
         zlabel('Temperature  ( ^oC)', 'FontSize', 12)
%
%   show parameters of current model
   display( '  ' );
   display( ['Gamma = ' num2str(Gamma_0) '   dx = '  num2str(dx) ] );
   display( [ 'Velocity (v_x, v_z) = (' num2str(c_x) ', ' num2str(c_z) ')' ] );
%
   P_x = rho_0*c_x*dx/Gamma_0;
   P_z = rho_0*c_z*dz/Gamma_0; 
   display( [ 'Grid Peclet numbers (P_x, P_z) = (' num2str(P_x) ', ' num2str(P_z) ')' ]);
%
      end   %   for_i_count = 1:length(N_vol)  ...
%
%
%
%
%

