%%  ESS 590D: HEAT AND MASS FLOW MODELLING  SPRING 2006
%%
%%  Finite Volume solution for 1-D or 2_D steady phi(x,z)
%%  on a rectangle
%%  compared with Finite Difference solution
%%
%%       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
%%
%%  THIS CODE COMPARES FDM AND FVM SOLUTIONS TO ANALYTICAL SOLUTION
%%  BASED ON POLYNOMIALS FOR CONDUCTIVITY AND SOURCE TERM
%%
%% 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 is independent of phi, but S_C 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 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Choose Finite Volumes or Finite Differences
%  --------------------------------------------
     fvm = 1;    %   use Control Volumes
  %%   fvm = 0;    %  use Finite Differences
%

%  set up control volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This code sets up uniform-sized volumes 
%   (could be generalized)
%
     Lx = 1;                  %  x-extent of domain
%
     nx = 2;                  %  number of control volumes in x
     dx = Lx/max((nx-1),1 );  %  x-widths of volumes
%
     Lz = 10;                  %  z-extent of domain
%
  %   nz =  6;                  %  number of control volumes in z
  %   nz = 11;                  %  number of control volumes in z
     nz = 21;                  %  number of control volumes in z
  %   nz = 41;                  %  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);
%
%

%  use Finite-Volume method
%
      if(fvm == 1)
        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,:)' ]';
%
     
      else
          
%  use Finite-Difference method
%     (don't ever really need X_edges, Z_edges)
         X_edges = X_P;
         Z_edges = Z_P;
      end         
%

%%  Boundary Conditions 
%%%%%%%%%%%%%%%%%%%%%%
%%  flux through east, west and down boundaries
%   q_d is negative because heat flows toward decreasing depths
     q_d_0 = -0.10;            %   W m^-2
  %%   q_d_0 = 0.0;            %   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 );
     q_w = repmat( q_w_0, nz, 1 );
     q_e = repmat( q_e_0, nz, 1 );
%
%
%%  temperature phi_u on upper boundary
  %%    phi_u_0 = -20;                   %  deg C
      phi_u_0 = 10;                   %  deg C
      phi_u = repmat( phi_u_0, 1, nx );
%
%
%%  material parameters
%%%%%%%%%%%%%%%%%%%%%%%
% set up thermal-property matrix - uniform conductivity
  %%   K_0 = 2.0;                     %   W m^-1 deg^-1    
     K_0 = 1;                     %   W m^-1 deg^-1    
     K_P   = repmat( K_0, size(X_P) );
%

     if(nx <= 2)   %   (2 half-volumes)
%
   % set up inverse-polynomial conductivity
   % -------------------------------------- 
    %%    z_K = z_edges_vec;
    %%    K_z = repmat(K_0, size(z_K) );
     
        z_K = [ 0 0.35 0.7 1 2 3 4 5 6.3 7  8  9   9.5  10 ];
        K_z = [ 2 1.7  1.4 1 3 1 1 2  2  3  1  2.5 1.5 0.5 ];
    %
    % get polynomial through K_z(z_K) (check it doesn't go negative!)
        P_K_inv = polyfit( z_K, 1./K_z, (length(z_K)-1) );
        
    %  get polynomial values at volume centers    
        K_P = 1./polyval( P_K_inv, Z_P );       

     end   %  if(nx <= 2)
     


%
%%  source term   
%%%%%%%%%%%%%%%%%%%%%
%   S = S_C + S_P * phi_P 

  % spatially uniform source term, no dependence on solution phi
      S_C_0  = 0;   %  W m^-3 s^-1
   %%   S_P_0  = 10;  %  W m^-3 s^-1 deg^-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
%

  %%   if(nx == 1)
     if(nx <= 2)
  %%   if(nx <= 0)
   % set up polynomial source
    %%    z_S = z_edges_vec;
    %%    S_z = repmat(S_P_0, size(z_S), 1 );
     
        z_S =       [ 0 0.5 1 2 3 4 5  6   7  8.5 9.7 10 ];
        S_z = 0.1 * [ 2 1.5 1 3 1 1 2  2.1 1.5  3 3.6 0 ];
  
    %  get polynomial fitting S_z(z_S)
   %%     P_S = polyfit( z_S, S_z, 2 ); 
        P_S = polyfit( z_S, S_z, (length(z_S)-1) );

    %
    %  get polynomial source term S_C at control volume centers
        S_C = polyval( P_S, Z_P );
     
     end   %  if(nx == 1)
%


%
   if (fvm == 1)
%
% This is it! solve for phi using Finite Volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    phi = fvm_soln_SS_bdry_nodes( X_edges, Z_edges, X_P, Z_P, K_P, ...
                              S_C, S_P, phi_u, q_d, q_e, q_w );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
   else
%
%
% This is it! solve for phi using Finite Differences
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      phi = fdm_soln_SS( X_edges, Z_edges, X_P, Z_P, K_P, S_C, S_P, ...
                            phi_u, q_d, q_e, q_w );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
   end    % if (fvm == 1 ) ...
%



    if(nx <= 2)

%%    find analytical solution in 1-D
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      zz = linspace(z_P_vec(1), z_P_vec(end), 200 );
      
  %  integrate polynomial for source term S_C  
       P_S_int = polyint( P_S );
%
  %  multiply S_int and 1/K to get int_x1^x [S(x')] dx' / K(x)
       S_int_over_K = polyval(P_S_int, zz).* polyval(P_K_inv, zz); 
%     
  %  find polynomial representing S_int_over_K
       n_poly = length(P_S_int) + length(P_K_inv) - 1;
       P_S_int_over_K = polyfit( zz, S_int_over_K, n_poly);
%
  %  now integrate it again
       P_S_int_over_K_int = polyint( P_S_int_over_K );
%
  %  integrate polynomial for inverse conductivity
       P_K_inv_int = polyint( P_K_inv ); 
%  
 
      phi_analytic = phi_u_0 - ( polyval(P_S_int_over_K_int, zz ) ...
                             - polyval(P_S_int_over_K_int, zz(1) ) ) ...
                           + ( polyval(P_K_inv_int, zz ) ...
                             - polyval(P_K_inv_int, zz(1) ) ) ...
                          .* ( polyval(P_S_int, zz(end) ) - q_d_0 );
      
    end   %    if(nx <= 2) ...



%%  Plot results
%%%%%%%%%%%%%%%%%
%

   if(fvm == 1)
       figure(1)   %  Finite-Volume solution
   else
       figure(2)   %  Finite-Difference solution
   end
%

   if( nx <= 2 )

%  1-D, plot phi(z) solution as line graph
%  ---------------------------------------
        plot(Z_P(:,1), phi(:,1), 'bo-' )
        hold on; grid on; box on
        
  %  analytical solution
        plot(zz, phi_analytic, 'r-')
        
        xlabel('Depth (m)')
        ylabel('Temperature (deg C)')
        
        if(fvm == 1)
           title('Finite Volume Method')
          legend( 'blue - fvm', 'Red - analytical', 'location', 'northwest')
        else 
           title('Finite Difference Method')
           legend( 'blue - fdm', 'Red - analytical', 'location', 'northwest')
        end
%

   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)')
                
        if(fvm == 1)
           title('Control Volume Method')
        else 
           title('Finite Difference Method')
        end

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

%
   if(nx <= 2)

%  plot conductivity (K_z)
%  -----------------------
      figure(3)
      hold on; grid on; box on
      plot( Z_P(:,1), 1./polyval( P_K_inv, Z_P(:,1) ), 'ro' )
      plot( zz, 1./polyval( P_K_inv, zz ), 'b-' )
      ylabel('Conductivity  (W m^{-1} deg^{-1})')
      xlabel('Depth  (m)')
      title('Polynomial Conductivity K_P(z)')


%  plot source (S_C)
%  -----------------
      figure(4)
      hold on; grid on; box on
      plot( Z_P(:,1), S_C(:,1), 'gx' )
      plot( Z_P, polyval( P_S, Z_P(:,1) ), 'ro' )
      
      plot( zz, polyval( P_S, zz ), 'b-' )

      ylabel('Source S_C  (W m^{-3})')
      xlabel('Depth  (m)')
      title('Polynomial Source Term S_C(z)')
   end  %    if(nx <= 2) ...
%
%
%
%
%

