function exercise5 ( )

%*****************************************************************************80
%
%% exercise5 solves a fluid flow problem using stream function and vorticity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise5:\n' );
  fprintf ( 1, '  Solve a fluid flow using stream function and vorticity.\n' );
%
%  Set problem sizes.
%
  nc = 41;
  nr = 41;
  nt = 81;
%
%  Set values of discrete variables.
%
  xmin = -10.0;
  xmax = +10.0;
  x = linspace ( xmin, xmax, nc );
  dx = ( xmax - xmin ) / ( nc - 1 );

  ymin = -10.0;
  ymax = +10.0;
  y = linspace ( ymin, ymax, nr );
  dy = ( ymax - ymin ) / ( nr - 1 );

  tmin = 0.0;
  tmax = 20.0;
  t = linspace ( tmin, tmax, nt );
  dt = ( tmax - tmin ) / ( nt - 1 );
%
%  Deal with periodic boundary condition by suppressing last node.
%
  x = x(1:nc-1);
  y = y(1:nr-1);

  [ X, Y ] = ndgrid ( x, y );
  X = rot90 ( X );
  Y = rot90 ( Y );

  nc = nc - 1;
  nr = nr - 1;

  n = nr * nc;
%
%  Create space for vorticity and stream function.
%
  O = zeros ( nr, nc );
  P = zeros ( nr, nc );
%
%  Build the system matrix once.
%  The linear system A*PSI=OMEGA represents the mathematical equation
%    Del^2 Psi = Omega
%
%  C,N,S,E,W are the indices of variables in each equation.
%  We are assuming periodic boundary conditions, so the indices wrap around.
%  I, J, A are the row, column and value of each nonzero of the matrix A.
%  
  [ C, N, S, E, W ] = cnsew ( nr, nc );

  I = [ C; C; C; C; C ];
  J = [ C; N; S; E; W ];
  A = [-2.0/dx^2 - 2.0/dy^2; 1.0/dy^2; 1.0/dy^2; 1.0/dx^2; 1.0/dx^2 ];
  A = repmat ( A, 1, n );
  A(1,1) = 1.0;
  A(2:5,1) = 0.0;

  I = reshape ( I, 5 * n, 1 );
  J = reshape ( J, 5 * n, 1 );
  A = reshape ( A, 5 * n, 1 );

  A = sparse ( I, J, A, n, n, 5 * n );

  figure ( 1 );
  spy ( A );
  title ( 'Sparsity plot for matrix A:' );
%
%  Set viscosity.
%
  nu = 0.001;
%
%  Carry out the time loop.
%
  for k = 1 : nt
%
%  On the first step, OMEGA is set by an initial condition.
%  On later steps, we update OMEGA using a forward Euler step:
%    dO/dt = dP/dY * dO/dX - dPdX * dOdY + nu (d^2 O/dx^2 + d^2 O/dy^2)
%
    if ( k == 1 )
      O = exp ( - 2.0 * X.^2 - Y.^2 / 20.0 ); 
    else
      O(C) = O(C) + dt * ( ...
         - ( 1.0 / dx / dy / 4.0 ) * ( P(E) - P(W) ) .* ( O(N) - O(S) ) ...
         + ( 1.0 / dy / dx / 4.0 ) * ( P(N) - P(S) ) .* ( O(E) - O(W) ) ...
         + ( nu / dx^2 ) * ( O(E) - 2.0 * O(C) + O(W) ) ...
         + ( nu / dy^2 ) * ( O(N) - 2.0 * O(C) + O(S) ) ...
      );
    end

    figure ( 2 )
    contourf ( X, Y, O, 20, 'edgecolor', 'none' );
    label = sprintf ( '\\omega, k = %d', k );
    title ( label );
%
%  Get Psi: Solve A * P = O
%
    o = reshape ( O, n, 1 );
    o(1) = 0.0;
    p = A \ o;
    P = reshape ( p, nr, nc );
    figure ( 3 )
    contourf ( X, Y, P, 20 );
    label = sprintf ( '\\psi, k = %d', k );
    title ( label );

    pause ( );

  end

  [ U, V ] = uv_from_psi ( X, Y, P );

  return
end
function [ U, V ] = uv_from_psi ( X, Y, PSI )

%*****************************************************************************80
%
%% uv_from_psi computes a velocity field (u,v) from a stream function.
%
%  Discussion:
%
%    We assume the x and y grid spacings are each uniform.
%
%    U =   d PSI/dy
%    V = - d PSI/dx
%
%  Modified:
%
%    12 June 2020
%
%  Input:
%
%    real X(nr,nc), Y(nc,nc): the node coordinates.
%
%    real PSI(nr,nc): the stream function.
%
%  Output:
%
%    real U(nr,nc), V(nr,nc): the velocity.
%
  [ nr, nc ] = size ( PSI );

  dx = X(1,2) - X(1,1);
  dy = Y(1,1) - Y(2,1);
  
  U = zeros ( nr, nc );
  V = zeros ( nr, nc );
%
%  For U, use centered differences except in first and last rows.
%  For Y derivatives, the larger Y value is at the SMALLER r value.
%
  for r = 1 : nr
    if ( r == 1 )
      U(r,1:nc) = ( PSI(r,1:nc) - PSI(r+1,1:nc) ) / dy;
    elseif ( r < nr )
      U(r,1:nc) = ( PSI(r-1,1:nc) - PSI(r+1,1:nc) ) / dy / 2.0;
    else
      U(r,1:nc) = ( PSI(r-1,1:nc) - PSI(r,1:nc) ) / dy;
    end
  end
%
%  For V, use centered differences except in first and last columns.
%  For X derivatives, the larger X value is at the larger c value.
%
  for c = 1 : nc
    if ( c == 1 )
      V(1:nr,c) = - ( PSI(1:nr,c+1) - PSI(1:nr,c) ) / dx;
    elseif ( c < nc )
      V(1:nr,c) = - ( PSI(1:nr,c+1) - PSI(1:nr,c-1) ) / dx / 2.0;
    else
      V(1:nr,c) = - ( PSI(1:nr,c) - PSI(1:nr,c-1) ) / dx;
    end
  end
%
%  Display PSI and (U,V) together.
%
  figure ( );
  subplot ( 1, 2, 1 );
  hold ( 'on' );
  contour ( X, Y, PSI, 20 );
  quiver ( X, Y, U, V );
  hold ( 'off' );
  xlabel ( '<-- x -->' );
  ylabel ( '<-- y -->' );
  title ( 'Velocity field from stream function' );

  UN = U ./ sqrt ( U.^2 + V.^2 );
  VN = V ./ sqrt ( U.^2 + V.^2 );

  I = find ( ( U.^2 + V.^2 ) == 0 );
  UN(I) = 0.0;
  VN(I) = 0.0;

  subplot ( 1, 2, 2 );
  hold ( 'on' );
  contour ( X, Y, PSI, 20 );
  quiver ( X, Y, UN, VN );
  hold ( 'off' );
  xlabel ( '<-- x -->' );
  ylabel ( '<-- y -->' );
  title ( 'Velocity directions from stream function' );

  return
end

