%%  ESS 590D: HEAT AND MASS FLOW MODELLING  SPRING 2006
%%
%%  Control Volume solution for 1-D or 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
%%
%%  THIS CODE SETS UP POSITIVE SOURCE FEEDBACK S_P
%%
%% 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 can vary with position (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      April 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%  set up Finite 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
     Lz = 4 * pi;                  %  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
  %   nz = 81;                  %  number of control volumes in z
  %   nz = 161;                  %  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 with fdm)
         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) );
%

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

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

S_P_0
two_K_over_dz_squared = 2*K_0/dz^2
%
   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 ) ...
%


%%  Plot results
%%%%%%%%%%%%%%%%%
%
   if(fvm == 1)
       figure(1)   %  Control-Volume solution
   else
       figure(2)   %  Finite-Difference solution (under construction)
   end
%

 if( nx <= 2 )

%  1-D, plot phi(z) solution as line graph
%  ---------------------------------------
   %%     plot(Z_P(:,1), phi(:,1), 'bo-' )
        plot(Z_P(:,1), phi(:,1), 'ro-' )
        hold on; grid on; box on

   %%     plot analytical solution
        zz = linspace(0, Lz);
%
        if( S_P_0 >= 0)
        beta = sqrt(S_P_0/K_0);            
   %%      (assumes q_d=0, S_P_0=1, and Lz is multiple of 2 pi))
           phi_zz = phi_u(1,1) * cos(beta * zz);
        else
        beta = sqrt(-S_P_0/K_0);
           b = phi_u(1,1)/(1 + exp(2*Lz) );
           a = phi_u(1,1) - b;
           phi_zz = a*exp(-beta.*zz) + b*exp(beta.*zz);
        end
%
           plot(zz,phi_zz, 'g')

%
        xlabel('Depth (m)')
        ylabel('Temperature (deg C)')
        
        if(fvm == 1)
           title('Finite Volume Method')
        else 
           title('Finite Difference Method')
        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 <= 1 ) ...


%
%
%
%
%

