%  ESS 524: HEAT AND MASS FLOW MODELLING  SPRING 2010
%
%  Finite Volume solution for 2_D steady phi(x,z)
%  on a rectangle
%
%       div( K grad phi) + S = 0
% or
%       del^2( phi ) + (grad K) dot (grad phi ) + S = 0 
%
% For example,
%  phi   = temperature                         (deg)
%  K     = thermal conductivity                (W m^-1 deg^-1)
%  S     = sources                             (W m^-3)
%        = S_C + S_P * phi
%
%
% Boundary Conditions follow Patankar's approach B 
%  - First define volume edges, then put nodes in centers.
%  - boundary nodes at centers of boundary faces.
%
%  - 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(x,z) varies with position
%
% Source 
%  S_bar = S_C + S_P * phi
%  S_C(x,z) varies with position
%  S_P(x,z) varies 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      April 2010
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%  set up finite volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This code sets up uniform-sized volumes 
%   (could be generalized)
%
     Lx = 10;                  %  x-extent of domain
%
     nx = 26;                  %  number of control volume edges in x
     dx = Lx/max((nx-1),1 );   %  x-widths of volumes
%
     Lz = 10;                  %  z-extent of domain
%
     nz = 26;                  %  number of control volume edges in z
     dz = Lz/max((nz-1),1 );   %  z-heights of volumes


  
%
%  Use Practice B (Patankar)
% ----------------------------
%  Set up matrices X_edges, Z_edges - defining (x,z) coordinates 
%  of finite-volume edges
     x_edge_vec = dx * ( 0:(nx-1) ); 
     z_edge_vec = dz * ( 0:(nz-1) )'; 
%
% duplicate values at domain boundaries for volumes of zero thickness 
     x_edge_vec = [x_edge_vec(1) x_edge_vec x_edge_vec(end) ];
     z_edge_vec = [z_edge_vec(1) z_edge_vec' z_edge_vec(end) ]';
%
%  Now form the volume edge matrices
     X_edges = repmat(x_edge_vec,nz+1,1);
     Z_edges = repmat(z_edge_vec,1,nx+1);
%
% Then put nodes at volume centers
%   (including zero-thickness volumes on domain boundaries)
     x_P_vec = (x_edge_vec(1:end-1) + diff(x_edge_vec)/2);
     z_P_vec = ( (z_edge_vec(1:end-1)' + diff(z_edge_vec)'/2) )';
%
     X_P = repmat(x_P_vec, (nz+1), 1);
     Z_P = repmat(z_P_vec, 1, (nx+1) );

%

%%%  Boundary Conditions 
%%%%%%%%%%%%%%%%%%%%%%%%%
%  flux through east, west and down boundaries.
% Add 1 to account for the "volumes" of zero thickness at domain edges
%
%   q_d is negative because heat flows toward decreasing depths
     q_d_0 = -0.10;             %   W m^-2
     q_e_0 = 0.0;               %   W m^-2
     q_w_0 = 0.0;               %   W m^-2
%
     q_d = repmat( q_d_0, 1, nx+1 );
     q_w = repmat( q_w_0, nz+1, 1 );
     q_e = repmat( q_e_0, nz+1, 1 );
%
%
%%%  Specified temperature phi_u on upper boundary
      phi_u_0 = 10;                   %  deg C
      phi_u   = repmat( phi_u_0, 1, nx+1 );
%
%
%
%%% material parameters
%%%%%%%%%%%%%%%%%%%%%%%
%
%   set up uniform background conductivity
    K_ref = 0.1;                     %   W m^-1 deg^-1 
    K_P   = repmat( K_ref, nz+1, nx+1 );
%
%   set up central block with different conductivity
%   block spans 40%-60% of depth nodes and width nodes
    K_block = 100 * K_ref;
  %%  K_block = 0.01 * K_ref;
%
    n_top   = fix(0.4*nz+2);
    n_bot   = fix(0.6*nz+2);
    n_left  = fix(0.4*nx+2);
    n_right = fix(0.6*nx+2);
%
    K_P(n_top:n_bot, n_left:n_right )= K_block;
%
%
%%%  source term   
%%%%%%%%%%%%%%%%%%%%%
%   S = S_C + S_P * phi_P 

  % spatially uniform source term, no dependence on solution phi
      S_C_0  = 0.0;   %  W m^-3 s^-1
      S_P_0  = 0.0;   %  W m^-3 s^-1 deg^-1
%
      S_C  = repmat( S_C_0, nz+1, nx+1 );  %  constant source term
      S_P  = repmat( S_P_0, nz+1, nx+1 );  %  phi_P-dependent source
%

%
%%% This is it! solve for phi using Finite Volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    phi = fvm_soln_SS_2D_Practice_B( X_edges, Z_edges, X_P, Z_P, K_P, ...
                              S_C, S_P, phi_u, q_d, q_e, q_w );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%  Plot results
%%%%%%%%%%%%%%%%%
%
%  Plot phi(x,z) solution as perspective surface
%  --------------------------------------------------
   figure(1)
   surf(X_P, Z_P, phi )
   colorbar
   hold on; grid on; box on
   xlabel('Distance  (m)')
   ylabel('Depth  (m)')
   zlabel('Temperature  (deg C)')
 %%   shading interp
                

%  plot conductivity (K_z)
%  -----------------------
   figure(2)
   surf(X_P, Z_P, K_P )
   colorbar
   hold on; grid on; box on
   xlabel('Distance  (m)')
   ylabel('Depth  (m)')
   zlabel('Conductivity  (W m^{-1} deg^{-1})')
                


%  plot source (S_C)
%  -----------------
   figure(3)
   surf(X_P, Z_P, (S_C + S_P.*phi ) )
   colorbar
   hold on; grid on; box on

   xlabel('Distance  (m)')
   ylabel('Depth  (m)')
   zlabel('Source  (W m^{-3})')
%
%
%
%
%