%  Homework on steady heat flow with source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Steady diffusion:
%   m-file solves 
%
%       d/dx[ k(x) d_phi/dx ] + S(x) = 0         (1)
%
%  on domain [x_0, x_1] with Boundary Conditions
%    d_phi/dx(x_0) = d_phi_dx_0
%    phi(x_1)      = phi_1
%
%  When S(x) and 1/k(x) are polynomials, there is an
%  analytical solution.
%  This m-file first develops this analytical solution, then
%  generates numerical solutions using 2 finite difference
%  formulations and 1 finite volume formulation.
%
%  User is prompted for the number of nodes Nx to use in the
%  numerical solution.
%
%
%                               E. Waddington  2008-04-25
%  ------------------------------------------------------
%
%  Define the domain
%  -----------------
    x_0 = 0;
    x_1 = 1;
%
%  set number of nodes Nx, and find node spacing dx  
    Nx = input('How many nodes?   ')
    dx = (x_1 - x_0)/(Nx-1);
%
%  set up nodal vector for plotting stuff at nodes
%  -----------------------------------------------
    x_vec = x_0: dx: x_1;
%
%  Set up detailed array x_fine_vec for plotting 
%  exact solutions
%  -----------------
    Nx_fine = 501;
    dx_fine = (x_1 - x_0)/(Nx_fine - 1 );
    x_vec_fine = x_0: dx_fine: x_1;
%
%
%  set Boundary conditions
%  -----------------------
%   specified gradient at x = x_0
      d_phi_dx_0 = 2;
%
%   specified value at x = x_1
      phi_1 = 5;
%

%  ================================
%  First get an analytical solution
%  ================================
%
%  Set the polynomials defining source and conductivity
%  ----------------------------------------------------
%   1/k(x) conductivity
      k_inv =[ 2 -2 1];
%
%   S(x) source
      S = 100*[ 2 -3 1 0];   
%
%
%  Integrate (1) from x_0 to x to find d_phi_dx,
%   then divide by k(x).
%  Boundary Condition d_phi_dx_0 appears naturally.
%
%   d_phi/dx =
%      [k(x_0) d_phi_dx_0 - int_(x_0)^x S(x') dx'] /k(x) (2)
%
%   The first term ( k(x_0) d_phi_dx_0 ) is flux q_0 at x_0  
%  --------------------------------------------------------
    q_0 = - d_phi_dx_0 / polyval(k_inv, x_0);
%
%  Integrate the source polynomial S(x)
%   Result S_int is also a polynomial 
%  ------------------------------------
    S_int = polyint(S);
%
%  Put q_0 and S_int (evaluated at the limits [x_0, x]) into (2)
%
%   d_phi/dx = [- q_0 - S_int(x) + S_int(x_0)] /k(x)   (3)
%
%   RHS of (3) is also a polynomial in x. 
%   We need to incorporate the constant terms into the polynomial
%   const_1 is a constant "polynomial"
%   (its only nonzero coefficient is coefficient of x^0)
%  --------------------------------------------------------
     const_1 =zeros( size(S_int) );
     const_1(end) = - polyval(S_int, x_0 ) + q_0;
%
%   Add const_1 to S_int(x) and divide by k(x) to get 
%   polynomial representation of d_phi/dx in (3)
%   ----------------------------------------------
     d_phi_dx = - conv( (S_int + const_1), k_inv );
%
%
%  Our other BC is at x_1, so now integrate d_phi_dx
%  from x to x_1 using phi_1 as BC. 
%   Result is phi(x) (also a polynomial) evaluated at 
%   the limits x_1 and x.
%
%   phi(x_1) - phi(x) = polyint(d_phi_dx)
%                       - polyval(polyint(d_phi_dx), x_1)
%
%  where the second BC phi(x_1) = phi_1 has appeared naturally
%  phi_const is constant polynomial used to add the constant 
%   term to phi polynomial
%     phi_const = desired value at x_1 - natural value of
%                 phi at x_1 
%  Natural polynomial phi(x) goes through origin, by construction
%  -------------------------------------------------------------
    phi_natural    = polyint(d_phi_dx);
    phi_const      = zeros( size(phi_natural) );
    phi_const(end) = phi_1 - polyval(phi_natural, x_1);
    phi = phi_natural + phi_const;
%

%  ===========
%  plot stuff
%  ===========
%
%  k(x)
%  ----
    k_vec      = 1./ polyval(k_inv, x_vec);
    k_vec_fine = 1./ polyval(k_inv, x_vec_fine);
%
%  Source S(x)
%  -----------
   S_vec      = polyval(S, x_vec);
   S_vec_fine = polyval(S, x_vec_fine);
%
%  solution phi(x)
%  ---------------
    phi_vec = polyval(phi, x_vec);
    phi_vec_fine = polyval(phi, x_vec_fine);
%
%
%
%  Conductivity k(x)
%  ------------------
    figure(1)
    hold on; grid on;
    plot(x_vec_fine, k_vec_fine, 'b-' )
    plot(x_vec, k_vec, 'rx' )
    xlabel('distance x')
    ylabel('conductivity   k(x)')
    set(gca,'ylim', [0,2])
%
%  Source  S(x)
%  ------------
    figure(2)
    hold on; grid on;
    plot(x_vec_fine, S_vec_fine, 'b-' )
    plot(x_vec, S_vec, 'rx' )
    xlabel('distance x')
    ylabel('source   S(x)')
%
%  Solution phi(x)
%  ---------------
    figure(3)
    hold on; grid on;
    plot(x_vec_fine, phi_vec_fine, 'k')
  %  plot(x_vec, phi_vec, 'ko')
    xlabel('distance  x')
    ylabel('solution   \phi(x)')
    title(' k- exact solution; r*-- fd solution #1; g+-- fd solution #2; bx-- fv solution')
%
%
%  ====================================================   
%    Now solve Equation (1) numerically with 3 schemes
%  ====================================================
%
%  Set up first 1-D finite-difference solution
%  ===========================================
%
%   Use product rule on derivatives term
%   -------------------------------------------
%  d/dx[k(x) d/dx(phi(x)) ] = dk/dx d_phi/dx + k(x) d^2phi/dx^2
%
%  Each derivative is then represented in standard fd way
%    [df/dx]_i = ( f_{i+1} - f_{i-1} )/ (2 dx_i)
%  
%    d^2_f/dx^2 = ( f_{i+1} - 2 f_i + f_{i-1} )/ dx_i^2
%  --------------------------------------------------
%
%  Initialize A matrix and RHS vector with all elements zero
%  ---------------------------------------------------------
    A_fd1 = zeros( length(x_vec) );
    RHS_fd1 = zeros(size(x_vec) )';
%
%  Put in the nonzero coefficients for interior points. 
%  This operation could be vectorized to avoid the loops,
%    but loops are easier to program quickly
%  ------------------------------------------
    for i = 2:length(x_vec)-1
        A_fd1(i,i-1) = ...
            ( k_vec(i) - ( k_vec(i+1) - k_vec(i-1) )/4 )/dx^2;
        A_fd1(i,i)   = -2*k_vec(i)/dx^2;
        A_fd1(i,i+1) = ...
            ( k_vec(i) + ( k_vec(i+1) - k_vec(i-1) )/4 )/dx^2;
  %      
        RHS_fd1(i) = - S_vec(i);
    end
%
%  Include boundary conditions
%  ---------------------------
%   At x_0 boundary point:  [ phi(2) - phi(1) ]/dx = d_phi_dx_0
       A_fd1(1,1) = -1/dx;
       A_fd1(1,2) =  1/dx;
       RHS_fd1(1) = d_phi_dx_0;
%
%   At x_1 boundary point:  phi(x_1) = phi_1
      A_fd1(end,end) = 1;
      RHS_fd1(end) = phi_1;                       ;
%
%
%   Now solve the finite-difference equations
%  ------------------------------------------
      phi_fd1_vec = A_fd1 \ RHS_fd1;

%
%   Add phi_fd1_vec to figure(3)
%  -----------------------------
     figure(3)
     plot(x_vec, phi_fd1_vec, 'r*--')
%
%
%
% Second fd solution
% ===================
%  Use first-order approximation to derivative term
%   [df/dx]_i = ( f_{i+1} - f_{i-1} ) /(2 dx)
%
%  So
%    d/dx[ k(x) d/dx(phi(x)) ] = 
%       k_{i+1} [d_phi/dx]_{i+1} - k_{i-1} [d_phi/dx]_{i-1} 
%
%  Then use first-order approximation on [d_phi_dx]_{i-1} etc
%    =  phi_{i-2} ( k_{i-1} /(4 dx^2) )
%     + phi_i     ( k_{i-1} - k_{i+1} ) /(4 dx^2)
%     + phi_{i+2} ( k_{i+1} /(4 dx^2) )  
%  --------------------------------------------------------------
%
%  Initialize matrices
%  -------------------
    A_fd2 = zeros( length(x_vec) );
    RHS_fd2 = zeros(size(x_vec) )';
%
%  Now, put in coefficients for interior points
%  --------------------------------------------
    for i = 3:length(x_vec)-2
        A_fd2(i,i-2) = k_vec(i-1)/(4*dx^2);
        A_fd2(i,i)   = - (k_vec(i+1) + k_vec(i-1)) / (4*dx^2);
        A_fd2(i,i+2) = k_vec(i+1)/(4*dx^2);
        
        RHS_fd2(i) = - S_vec(i);
    end
%
%  Put in coefficients for flux boundary at x_0.
%  With this fd approximation, which stretches over 5 grid points,
%  we need 2 equations for each boundary.
%   --------------------------------------
%   At x_0 boundary point:  [ phi(2) - phi(1) ]/dx = d_phi_dx_0
%   Same as in method #1
       A_fd2(1,1) = -1/dx;
       A_fd2(1,2) =  1/dx;
       
       RHS_fd2(1) = d_phi_dx_0;
%
%     At adjacent interior point:  
%     *** explain this ***
       A_fd2(2,1) = -k_vec(3)/(4*dx^2);
       A_fd2(2,3) =  k_vec(3)/(4*dx^2);
       
       RHS_fd2(2) = - S_vec(2) + d_phi_dx_0 * ( k_vec(1)/(2*dx) );
       
%
%   put in coefficients for fixed phi_1 at x_1 boundary
%   ---------------------------------------------------
 %  
%   x_{end-1} equation     
      A_fd2(end-1,end)   =  k_vec(end) /(1.5*dx^2);
      A_fd2(end-1,end-1) = - ( 2*k_vec(end) + k_vec(end-2) ) /(3*dx^2); 
      A_fd2(end-1,end-3) =  k_vec(end-2) /(3*dx^2); 
%         
      RHS_fd2(end-1) = - S_vec(end-1);        
%
%    x_end equation
      A_fd2(end,end) = 1;
      RHS_fd2(end) = phi_1;
%
%
%
%   Now solve the finite-difference equations
%  ------------------------------------------
      phi_fd2_vec = A_fd2 \ RHS_fd2;

%
%   Add phi_fd2_vec to figure(3)
%  -----------------------------
     figure(3)
     plot(x_vec, phi_fd2_vec, 'g+--')
%
%
%
%  Now set up finite-volume solution
%  ==================================
%
%  get interface conductivities k_e and k_w
%   using gemetric mean (Patankar, p. 44)
%  -----------------------------------------
      k_E = [ k_vec(2:end) k_vec(end) ];
      k_P = k_vec;
      k_W = [ k_vec(1) k_vec(1:(end-1) ) ];
      k_e = 2 * k_P .* k_E ./ (k_E + k_P); 
      k_w = 2 * k_P .* k_W ./ (k_W + k_P); 
%
%  get source vector
%  -----------------
      S_C = S_vec;
%
%  get coefficients for interior points
%  ------------------------------------
      a_E = k_e ./ dx;
      a_W = k_w ./ dx;
      a_P = a_E + a_W; 
      b   = S_C .* dx;
%
%   stuff coefficients into matrix A_fd2
%  -------------------------------------
      A_fvm   = zeros(length(x_vec) );
      RHS_fvm = zeros(size(x_vec))';
%
      for i=2:length(x_vec)-1
          A_fvm(i,i)   = -a_P(i);
          A_fvm(i,i-1) =  a_W(i);
          A_fvm(i,i+1) =  a_E(i);
          
          RHS_fvm(i)  = -b(i);
      end
%
%  now set up flux BC at x_0
   %  q_B - k_W(2) * ( phi(2) - phi(1) ) / dx + S_C(1) * dx = 0
% -------------------------------------------------------------
     A_fvm(1,1) = -k_w(2)/dx;
     A_fvm(1,2) =  k_w(2)/dx;
     RHS_fvm(1) = -(S_C(1)*dx - k_P(1) * d_phi_dx_0 );
%
%  now set up set phi at x_1
%  -------------------------
     A_fvm(end,end) = 1;
     RHS_fvm(end)   = phi_1;
%
%  Solve it!
%  ----------
     phi_fvm_vec = A_fvm \ RHS_fvm;
%

%   Add plot to figure(3)
%  ----------------------
    figure(3)
    plot( x_vec, phi_fvm_vec, 'bx--')
%
%
%  Now set up finite element solution
%  ----------------------------------
%   ha ha ha !!
%   Later, maybe...
%
%
%