%%  ESS 590D: HEAT AND MASS FLOW MODELLING  SPRING 2006
%%
%%  Finite Volume solution for 1-D or 2_D transient phi(x,z)
%%  on a rectangle
%%  compared with Finite Difference solution
%%
%%       div( Gamma grad phi) + S = d(rho phi)/dt
%%
%% For example,
%%  phi   = internal energy                     (J kg^-1)
%%        = int_0^T c(T') dT'
%%  T     = 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
%%  rho   = density
%%
%%  THIS CODE (what's special?)
%%
%% Boundary Conditions follow Patankar's approach - half-volumes on
%% boundaries
%%  - Specified quantity phi_u(x) on upper surface
%%  - specified fluxes q_e, q_w, and q_d through other 3 sides of domain
%%
%% Material Properties
%%  K varies with position (x,z)
%%
%% Source 
%%  Source S = S_C + S_P * T
%%
%%  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.
%%
%%
%%  This main program sets and delivers:
%%    - Mesh geometry (X_P, Z_P, X_edges, Z_edges)
%%    - time-marching info (nt, dt)
%%    - constitutive properties (K_P)
%%    - sources (S_C, S_P)
%%    - initial conditions (phi_0) 
%%    - boundary conditions (phi_u, q_d, q_e, q_w)
%%
%%
%%                                   Ed Waddington      April 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%
%%  set up finite volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This code sets up uniform-sized volumes in both space and time 
%   (could be generalized)
%
%  Space dimensions
%  -----------------
     Lx = 20;                  %  x-extent of domain
%
     nx = 21;                  %  number of control volumes in x
     dx = Lx/max((nx-1),1 );  %  x-widths of volumes
%                                 max function prevents division by zero 
%                                 in 1-D (i.e. when nx=1)
%
     Lz = 20;                 %  z-extent of domain
%
     nz = 21;                 %  number of control volumes in z
     dz = Lz/max((nz-1),1 );  %  z-heights of volumes


%%  set up matrices X_P, Z_P - defining (x,z) coordinates 
%%    of finite-volume centers
     x_P_vec = dx * (0:(nx-1) );
     z_P_vec = dz * (0:(nz-1) )'; 
     X_P = repmat(x_P_vec,nz,1);
     Z_P = repmat(z_P_vec,1,nx);
%
%
%%  set up matrices X_edges, Z_edges - defining (x,z) coordinates 
%%    of finite-volume interfaces
        X_edges = [ X_P(:,1) (X_P(:, 1:end-1) + diff(X_P')'/2) ...
                             X_P(:,end) ];
%
        Z_edges = [ Z_P(1,:)' (Z_P(1:end-1, :) + diff(Z_P)/2)' ...
                             Z_P(end,:)' ]';
%
%
%  Time dimension
%  ---------------
     sec_per_year = 365.25*24*3600;
  %%   Lt = 5;                     % Total duration of time integration
  %%   nt = 101;                     % number of future points
     Lt = 1;                     % Total duration of time integration
     nt = 51;                     % number of future points
     dt_yr = Lt/max((nt-1),1);   % time step in years
     dt = dt_yr * sec_per_year;
%

%
%%  Boundary Conditions 
%%%%%%%%%%%%%%%%%%%%%%%%%
%%  flux through east, west and down boundaries
%   q_d is negative because heat flows toward decreasing depths
     q_d_0 = -0.50;            %   W m^-2
     q_e_0 = 0.0;               %   W m^-2
   %%  q_w_0 = -0.50;               %   W m^-2
     q_w_0 = 0.0;               %   W m^-2
%
%  This code sets of simple steady B.C.s
%  put "interesting" transient B.C.s into these matrices
     q_d = repmat( q_d_0, nt, nx );
     q_w = repmat( q_w_0, nz, nt );
     q_e = repmat( q_e_0, nz, nt );
%
%
%%  temperature phi_u on upper boundary
%%  put your favorite pattern here
%   -------------------------------------
      phi_u_0 = -20;                   %  deg C
%
  %  constant and uniform surface temperature
     %% phi_u = repmat( phi_u_0, nt, nx );
%
%    seasonally varying but spatially uniform
      phi_range = 5;
      time_vec = dt_yr * (0:(nt-1));
      phi_seasonal = phi_u_0 + phi_range * sin( 2*pi * time_vec );
      phi_u = repmat( phi_seasonal', 1, nx );
%

%
%%  material parameters
%%%%%%%%%%%%%%%%%%%%%%%
% set up thermal-property matrix - e.g. uniform conductivity
     K_0 = 2.0;                     %   W m^-1 deg^-1    
     K_P = repmat( K_0, size(X_P) );
%
%
% set up density matrix - e.g. uniform on domain 
     rho_0 = 900;                   %  kg m^-3   (e.g. ice)
     rho_P = repmat( rho_0, size(X_P) );
%
%
% set up specific heat capacity - uniform on domain 
     c_0 = 2000;                    %  J kg^-1 deg^-1
     c_P = repmat( c_0, size(X_P) );
%
%
% set up generic Gamma matrix - for heat flow, Gamma = K/c 
  %%   Gamma_P = K_P./c_P;
     rho_c_P = rho_P.* c_P;



%
%%  source term   
%%%%%%%%%%%%%%%%%%%%%
%   S = S_C + S_P * phi_P 
%  If you want transient source terms, you need to build that information
%  into a transient vector S(t), or (storage expensive) 3-D arrays
%   S_C(x,z,t) and S_P(x,z,t)

  % spatially uniform source term, no dependence on solution phi
      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
%
%
%
%  initial condition
%  -----------------
  %%   phi_0 = zeros( size(X_P) );
  %%   phi_0 = repmat( phi_u, nz,1 )...
  %%       - [ Z_P(1,:)' cumsum(diff(Z_P))']'.* repmat(q_d, nz,1)./ K_P;
%
%  For example, set constant initial value, equal to BC on upper boundary
      phi_0 = repmat( phi_u(1,:), nz, 1 );

%

%
% This is it! solve for phi(x,z,t) using Finite Volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    phi = fvm_soln_transient_bdry_nodes( X_edges, Z_edges, X_P, Z_P, ...
                                        nt, dt, rho_c_P, K_P, S_C, S_P, ...
                                        phi_0, phi_u, q_d, q_e, q_w );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
dt
rho_c_dz_2_over_2K = rho_c_P(1,1)*(Z_P(2,1)-Z_P(1,1)).^2/(2 * K_P(1,1))

%
%%  Plot results
%%%%%%%%%%%%%%%%%
%
   if(nx <= 2)
%
%% 1-D line graphs
%%%%%%%%%%%%%%%%%%
%
%  plot conductivity (K_z)
%  -----------------------
      figure(1)
      hold on; grid on; box on
      plot( Z_P(:,1), K_P(:,1), 'ro' )
      ylabel('Conductivity  (W m^{-1} deg^{-1})')
      xlabel('Depth  (m)')
      title('Conductivity  K_P')


%  plot source (S_C)
%  -----------------
      figure(2)
      hold on; grid on; box on
      plot( Z_P(:,1), S_C(:,1), 'rx' )
      ylabel('Source S_C  (W m^{-3})')
      xlabel('Depth  (m)')
      title('Source Term S_C')
%

%  plot source (S_C)
%  -----------------
      figure(3)
      hold on; grid on; box on
      plot( Z_P(:,1), S_P(:,1), 'bx' )
      ylabel('Source S_P  (W m^{-3} deg^{-1})')
      xlabel('Depth  (m)')
      title('Source Term S_P')
%
%  1-D, plot phi(z) solution as line graph
%  ---------------------------------------
      figure(4)
      plot(Z_P(:,1), phi(:,1), 'bo-' )
      hold on; grid on; box on              
      xlabel('Depth (m)')
      ylabel('Temperature (deg C)') 
      title('Finite Volume Method')
%
%
%
   else
%
%%  plot in 2-D
%%%%%%%%%%%%%%%%%%%
%
%  plot conductivity (K_z)
%  -----------------------
      figure(1)
      hold on; grid on; box on
        surf(X_P, Z_P, K_P )
        colorbar
      zlabel('Conductivity  (W m^{-1} deg^{-1})')
      xlabel('Depth  (m)')
      title('Conductivity K_P')


%  plot source (S_C)
%  -----------------
      figure(2)
      hold on; grid on; box on
        surf(X_P, Z_P, S_C )
        colorbar
      zlabel('Source S_C  (W m^{-3})')
      xlabel('Distance  (m)')
      ylabel('Depth  (m)')
      title('Source Term S_C')
%

%  plot source (S_P)
%  -----------------
      figure(3)
      hold on;
      surf(X_P, Z_P, S_P )
      colorbar
      zlabel('Source S_P  (W m^{-3} deg^{-1})')
      xlabel('Distance  (m)')
      ylabel('Depth  (m)')
      title('Source Term S_P')
%
       
%  2-D, plot phi(x,z) solution as perspective surface
%  --------------------------------------------------
      figure(4)
      hold on
      surf(X_P, Z_P, phi )
      colorbar
      xlabel('Distance  (m)')
      ylabel('Depth  (m)')
      title('Control Volume Method')

   end   %   if( nx <= 2 ) ...

%
%
%
%

