%%  ESS 590D: HEAT AND MASS FLOW MODELLING  SPRING 2006
%%
%%  Finite-Volume solution for steady phi(x,z)
%%  diffusion with advection
%%
%%       div( Gamma grad phi) + S = div( rho u phi )
%%
%% 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 INCLUDES 2 OPTIONS FOR ADVECTIVE FLUX 
%%   (1) (upwind = 0) CENTERED DIFFERENCE (NO UPWINDING SCHEME)
%%       *** DON'T USE THIS AT HOME FOR SERIOUS WORK!!  ***
%%
%%   (2) (upwind = 1) PATANKAR POWER-LAW SCHEME
%%
%%
%%    EXTERNALLY SPECIFIED VELOCITY MUST SATISFY CONTINUITY,
%%     i.e. u IS A CONSTANT IN 1-D
%%
%%
%%  Boundary Conditions:
%%  vector phi_u       is fixed value of phi on upper boundary
%%  vector phi_d       is fixed value of phi on lower boundary
%%  vector dphi_dx_w   is fixed gradient on west boundary
%%  vector dphi_dx_e   is fixed gradient on east boundary
%%
%% In this treatment, gradient boundary conditions (e-w) are applied on 
%% edges of finite volumes. This means that I don't have to use half-volumes.
%% Function-value boundary conditions (u-d) are applied on nodes, i.e. with
%% half-volumes at Up and Down boundaries.
%%
%%  THIS CODE APPLIES GRADIENT BOUNDARY CONDITIONS ON VOLUME INTERFACES
%%  on East and West
%%  APPLIES FIXED-VALUE BOUNDARY CONDITIONS ON NODES
%%  in Up and Down directions
%%
%% Material Properties
%%  K can vary with position - K(x,z)
%%
%%
%%  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 6 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    global upwind   Peclet

% decide whether to use centered difference   (upwind = 0)
%                  or power-law upwind scheme (upwind = 1)
     upwind = 1
        
%  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 ' ]; 


%  set up finite volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This code sets up uniform-sized volumes 
%   (could be generalized)
%
   Lx = 10;                   %  x-extent of domain
%
   nx = 1;                   %  number of control volumes in x
   dx = Lx/nx;               %  x-widths of volumes
%
   Lz = 10;                  %  z-extent of domain
%
%
% how many different grids to solve with
   jz_n = [3 6 21 51 ];
%
   n_models = length(jz_n);
%

%  Loop over models - various grid spacings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   for jz = 1:n_models
%
      pause
%
       nz = jz_n(jz);              %  number of control volumes in z
%
     
       dz = Lz/max((nz-1),1 );   %  z-heights of volumes

%  in vertical direction, use half volumes and put nodes on boundaries 
       z_P_vec = dz * (0:(nz-1) )'; 
       Z_P = repmat(z_P_vec,1,nx);
%
       Z_edges = [ Z_P(1,:)' (Z_P(1:end-1, :) + diff(Z_P)/2)' ...
                          Z_P(end,:)' ]';
%

%  set up matrices X_edges, Z_edges - defining (x,z) coordinates 
%  of control-volume edges
       x_edges_vec = dx * (0:nx);
       X_edges = repmat(x_edges_vec,nz+1,1);
       X_P = X_edges(:, 1:end-1) + diff(X_edges')'/2;
       X_P = X_P(1:end-1,:);
%
%

%%  Boundary Conditions 
%%%%%%%%%%%%%%%%%%%%%%%%%
%%  phi gradients at east, west  boundaries
     dphi_dx_e_0 = 0.0;               %   deg m^-1
     dphi_dx_w_0 = 0.0;               %   deg m^-1
%
     dphi_dx_w = repmat( dphi_dx_w_0, nz, 1 );
     dphi_dx_e = repmat( dphi_dx_e_0, nz, 1 );
%
%
%%  temperature phi_u on upper boundary
      phi_u_0 = 20;                   %  deg C
      phi_u   = repmat( phi_u_0, 1, nx );
%
%
%%  temperature phi_d on lower boundary
      phi_d_0 = 40;                   %  deg C
      phi_d   = repmat( phi_d_0, 1, nx );
%
%
%%  material parameters
%%%%%%%%%%%%%%%%%%%%%%%
% set up thermal-property matrix - uniform conductivity
     Gamma_0 = 1;                     %   W m^-1 deg^-1    
     Gamma_P = repmat( Gamma_0, size(X_P) );
%
%
% set up density matrix
     rho_0 = 1;    %   kg m^-3
     rho = repmat( rho_0, size(X_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, size(X_P) );  %  constant source term
     S_P  = repmat( S_P_0, size(X_P) );  %  phi_P-dependent source
%


%
% This is it! solve for phi using Finite Volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    phi = fvm_soln_SS_adv( X_edges, Z_edges, X_P, Z_P, Gamma_P, rho, ...
                    S_C, S_P, phi_u, phi_d, dphi_dx_e, dphi_dx_w );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jz
dz
Peclet

%%  Plot results
%%%%%%%%%%%%%%%%%
%
    if(upwind == 0)
       figure(1)
    else
        figure(2)
    end
%
   if( nx <= 2 )
%
%  1-D, plot phi(z) solution as line graph
%  ---------------------------------------
      plot(Z_P(:,1), phi(:,1), 'Color', color(:,jz), 'Marker', marker(jz) )
      hold on; grid on; box on

      xlabel('Depth (m)')
      ylabel('Temperature (deg C)')
      set(gca,'xlim',[min(min(Z_edges)), max(max(Z_edges))] )
      set(gca, 'ylim', [0, max(phi_u_0,phi_d_0)])

   else
       
%  2-D, plot phi(x,z) solution as perspective surface
%  --------------------------------------------------
        surf(X_P, Z_P, phi )
        colorbar
        hold on; grid on; box on
        xlabel('Distance  (m)')
        ylabel('Depth  (m)')
        zlabel('Temperature  (deg C)')
                
   end   %   if( nx <= 2 ) ...
%
   end    %  for jz = 1:n_models
%
%
%
%