program adiprc

c General information:

c     Program adiprc solves the two-dimensional convection equation:
c
c         du/dt + c_x du/dx + c_y du/dy =0
c     for
c         0 <= x <= x_max  0 <= y <= y_max  0 <= t <= t_max.

c     For simplicity it is assumed that the convection velocity
c     components c_x and c_y are positive, so that the boundary
c     conditions apply at x = 0 and at y = 0.  Reverse x or y to
c     handle negative c_x or c_y.

c     Program adiprc reads its parameters from an input file
c     "adiprcon.in".  See the comments in that file for the required
c     input parameters.

c     Program adiprc produces the output file "adiprcon.out". Any
c     existing file of that name will be deleted.

c     In the present implementation, program adiprc requires an exact
c     solution;
c        double precision function u_exa(x,y,t,c_x,c_y),
c     to give it its initial and boundary conditions. The numerical
c     errors are then determined by comparing the computed solution
c     to this exact solution.

c     Program adiprc also requires subroutines:
c        in.f: for reading the input file
c        out.f: for doing output
c        tridia.f: for solving tridiagonal systems of equations
c        openw.f: for opening files for writing

c     Compile as:
c        g77 -o adiprcon adiprcon.f u_exa.f openw.f tridia.f out.f in.f

c Numerical scheme:

c     Program adiprc uses the two-stage Peaceman-Rachford ADI method,
c     which advances the solution over a time step by means of two half
c     time steps:
c
c      (1/(dt/2) + c_x d/dx) u-tilde = (1/(dt/2) - c_y d/dy) u^n
c      (1/(dt/2) + c_y d/dy) u^n+1   = (1/(dt/2) - c_x d/dx) u-tilde
c
c     where u^n is the dependent variable u at the last known time,
c     u^n+1 is the dependent variable at the next time to be computed;
c     dt is the numerical time step; u-tilde a purely computational
c     "unknown at the middle of the time step", and the derivatives d/dx
c     and d/dy are evaluated using standard second order central
c     differences.  This method is second order accurate in space and
c     time and unconditionally stable.

c     A boundary condition for u-tilde at x=0 can be obtained by
c     substracting the two Peaceman-Rachford formulae from each other,
c     to give:
c
c        2 u-tilde = (1 - 0.5 dt c_y d/dy) u^n + (1 + 0.5 dt c_y d/dy) u^n+1
c
c     where u^n and u^n+1 should be known at x=0 from the given boundary
c     condition there.  (In this program, they are obtained from u_exa.)

c     Because of the central differencing, the scheme requires
c     artificial boundary conditions at x = x_max and at y = y_max.
c     At these boundaries, linear extrapolation from the interior
c     is used.

c Program variables:

c     Avoid typos:
      implicit none

c     Components of the convection velocity in the x and y directions:
      double precision c_x,c_y

c     Index in the x-direction with its maximum and dimensioned value:
      integer j,j_max,j_dim
      parameter (j_dim=128)

c     Index in the y-direction with its maximum and dimensioned value:
      integer l,l_max,l_dim
      parameter (l_dim=128)

c     Coordinates x and y and their limits:
      double precision x(0:j_dim),y(0:l_dim),x_max,y_max

c     Time t:
      double precision t

c     Time step number and its maximum value:
      integer n,n_max

c     Spatial mesh size in the x direction:
      double precision dx

c     Spatial mesh size in the y direction:
      double precision dy

c     The time step:
      double precision dt

c     The largest of the two Courant numbers c_x dt / dx and c_y dt / dy:
      double precision courmx

c     The values of the dependent variable for the last computed time:
      double precision u(0:j_dim,0:l_dim)

c     The intermediate unknowns u-tilde in the ADI method:
      double precision utilde(0:j_dim,0:l_dim)

c     The left (x=0) boundary values of u at the next time t^(n+1):
      double precision uleft(0:l_dim)

c     Storage for the tridiagonal equations to be solved:

c     Maximum tridiagonal system size (must be the maximum of j_dim and l_dim):
      integer td_dim
      parameter (td_dim=j_dim+l_dim)
c     The three nonzero matrix diagonals (the coefficients of the unknowns):
      double precision diags(3,0:td_dim)
c     The right hand sides of the equations:
      double precision b(0:td_dim)
c     The values of the unknowns as found by tridia:
      double precision sol(0:td_dim)

c     Number of nonzero output times:
      integer no_out
c     First nonzero output time:
      double precision t_out1
c     Corresponding time step index:
      integer n_out1

c     Whether an old output file exists:
      logical exists

c     An exact solution of the equation, used to test the method:
      double precision u_exa
      external u_exa

c     Constants, defined to simplify changing precision:
      double precision zero,one,two,o2
      parameter (zero=0.d0,one=1.d0,two=2.d0,o2=0.5d0)

c Executable statements:

c     Open the output file:
      call openw(1,'adiprcon.out')

c Initialize the variables:

c     Read the input file:
      call in(j_dim,l_dim,td_dim,
     & j_max,l_max,c_x,c_y,x_max,y_max,courmx,no_out,t_out1)

c     The mesh-size in the x-direction:
      dx=x_max/j_max

c     The mesh-size in the y-direction:
      dy=y_max/l_max

c     The time step:
      dt=min(courmx*dx/c_x,courmx*dy/c_y)

c     The time step index of the first nonzero output time:
      n_out1=nint(t_out1/dt)
      write(1,300)n_out1
 300  format('The time step of the first nonzero output time:',i4)
      if(n_out1.lt.1)then
         print310
         write(1,310)
 310     format('*** WARNING: The output times are spaced too closely;',
     &    ' output will be done '/'             every step instead.')
         n_out1=1
      endif

c     The maximum time step index:
      n_max=no_out*n_out1

c     Initialize the time index:
      n=0

c     Initialize the time:
      t=zero

c     Set the initial data from the exact solution:
      do 500 j=0,j_max
         x(j)=(j+zero)/j_max*x_max
         do 490 l=0,l_max
            y(l)=(l+zero)/l_max*y_max
            u(j,l)=u_exa(x(j),y(l),t,c_x,c_y)
 490     continue
 500  continue

c     Output the initial data:
      call out(j_dim,l_dim,j_max,l_max,c_x,c_y,x,y,t,n,u)

c Compute the solution at times greater than zero:

      do 5000 n=0,n_max-1

c        The time t^(n+1) being computed:
         t=(n+1)*dt

c        Store the left boundary condition at n+1, taken from u_exa, first:
         do 2100 l=0,l_max
            uleft(l)=u_exa(x(0),y(l),t,c_x,c_y)
 2100    continue

c        Compute the intermediate unknowns u-tilde in the PR method
c        for each horizontal line l=constant in turn:

         do 3000 l=1,l_max-1

c           Put the coefficients of the boundary condition at x=0 (j=0),
c              2 u-tilde
c                 = (1 - 0.5 dt c_y d/dy) u^n + (1 + 0.5 dt c_y d/dy) u^n+1
c           in the array diags, and the right hand side in array b:

c           Zero the coefficient of the nonexisisting unknown at j=-1:
            diags(1,0)=zero
c           The coefficient of unknown u-tilde_0 itself, as above:
            diags(2,0)=two
c           The coefficient of u-tilde_1, which is not there:
            diags(3,0)=zero
c           The right hand side:
            b(0)=u(0,l)-o2*dt*c_y/(two*dy)*(u(0,l+1)-u(0,l-1))
     &       +uleft(l)+o2*dt*c_y/(two*dy)*(uleft(l+1)-uleft(l-1))

c           Put the coefficients of the finite difference equation,
c              (1/(dt/2) + c_x d/dx) u-tilde = (1/(dt/2) - c_y d/dy) u^n
c           in the array diags, and the right hand side in array b:

            do 2500 j=1,j_max-1
c              The coefficient of u-tilde_(j-1):
               diags(1,j)=-c_x/(two*dx)
c              The coefficient of u-tilde_j:
               diags(2,j)=two/dt
c              The coefficient of u-tilde_(j+1):
               diags(3,j)=c_x/(two*dx)
c              The right hand side:
               b(j)=two/dt*u(j,l)-c_y/(two*dy)*(u(j,l+1)-u(j,l-1))
 2500       continue

c           We do not have a boundary condition at x_max, so we
c           make one up: we approximate u-tilde_(j_max) by linear
c           extrapolation from u-tilde_(j_max-1) and u-tilde_(j_max-2):
c
c              u-tilde_(j_max) = 2 u-tilde_(j_max-1) - u-tilde_(j_max-2)
c
c           This produces the equation:
c
c             1 u-tilde_(j_max-2) -2 u-tilde_(j_max-1) + 1 u-tilde_(j_max) = 0
c
c           Unfortunately, there is no way to accomodate the coefficient
c           1 of u-tilde_(j_max-2) in diags.  To fix this up, we
c           can substract a multiple of the second-last equation,
c
c              diags(1,j_max-1) u-tilde_(j_max-2)
c              + diags(2,j_max-1) u-tilde_(j_max-1)
c              + diags(3,j_max-1) u-tilde_(j_max)
c              = b(j_max-1)
c
c           Clearly, the multiple to substract to eliminate the 1 is
c           1/diags(1,j_max-1).  That gives the combined equation:
c
c              0 u-tilde_(j_max-2)
c              + (-2 - diags(2,j_max-1)/diags(1,j_max-1)) u-tilde_(j_max-1)
c              + ( 1 - diags(3,j_max-1)/diags(1,j_max-1)) u-tilde_(j_max)
c              = - b(j_max-1)/diags(1,j_max-1)
c
c           so we get:
            diags(1,j_max)=-two-diags(2,j_max-1)/diags(1,j_max-1)
            diags(2,j_max)=one-diags(3,j_max-1)/diags(1,j_max-1)
            diags(3,j_max)=zero
            b(j_max)=-b(j_max-1)/diags(1,j_max-1)

c           Now let tridia solve this system:
            call tridia(j_max+1,diags,b,sol)

c           Put the solution found by tridia into array u-tilde:
            do 2700 j=0,j_max
               utilde(j,l)=sol(j)
 2700       continue

 3000    continue

c        Compute the u^(n+1) values in the PR method for each vertical
c        line j=constant in turn:

         do 4000 j=1,j_max-1

c           Put the coefficients of the boundary condition at y=0 (l=0),
c              u_0^(n+1) = given (taken here from u_exa)
c           in the array diags, and the right hand side in array b:

c           Zero the coefficient of the nonexisisting unknown at l=-1:
            diags(1,0)=zero
c           The coefficient of unknown u_0 itself, as above:
            diags(2,0)=one
c           The coefficient of u_1, which is not there:
            diags(3,0)=zero
c           The right hand side:
            b(0)=u_exa(x(j),y(0),t,c_x,c_y)

c           Put the coefficients of the finite difference equation,
c              (1/(dt/2) + c_y d/dy) u^n+1   = (1/(dt/2) - c_x d/dx) u-tilde
c           in the array diags, and the right hand side in array b:

            do 3500 l=1,l_max-1
c              The coefficient of u_(l-1):
               diags(1,l)=-c_y/(two*dy)
c              The coefficient of u-tilde_l:
               diags(2,l)=two/dt
c              The coefficient of u-tilde_(l+1):
               diags(3,l)=c_y/(two*dy)
c              The right hand side:
               b(l)=two/dt*utilde(j,l)-c_x/(two*dx)
     &          *(utilde(j+1,l)-utilde(j-1,l))
 3500       continue

c           We do not have a boundary condition at y_max, so we
c           make one up: we approximate u-tilde_(l_max) by linear
c           extrapolation from u-tilde_(l_max-1) and u-tilde_(l_max-2):
c
c              u-tilde_(l_max) = 2 u-tilde_(l_max-1) - u-tilde_(l_max-2)
c
c           This produces the equation:
c
c             1 u-tilde_(l_max-2) -2 u-tilde_(l_max-1) + 1 u-tilde_(l_max) = 0
c
c           Unfortunately, there is no way to accomodate the coefficient
c           1 of u-tilde_(l_max-2) in diags.  To fix this up, we
c           can substract a multiple of the second-last equation,
c
c              diags(1,l_max-1) u-tilde_(l_max-2)
c              + diags(2,l_max-1) u-tilde_(l_max-1)
c              + diags(3,l_max-1) u-tilde_(l_max)
c              = b(l_max-1)
c
c           Clearly, the multiple to substract to eliminate the 1 is
c           1/diags(1,l_max-1).  That gives the combined equation:
c
c              0 u-tilde_(l_max-2)
c              + (-2 - diags(2,l_max-1)/diags(1,l_max-1)) u-tilde_(l_max-1)
c              + ( 1 - diags(3,l_max-1)/diags(1,l_max-1)) u-tilde_(l_max)
c              = - b(l_max-1)/diags(1,l_max-1)
c
c           so we get:
            diags(1,l_max)=-two-diags(2,l_max-1)/diags(1,l_max-1)
            diags(2,l_max)=one-diags(3,l_max-1)/diags(1,l_max-1)
            diags(3,l_max)=zero
            b(l_max)=-b(l_max-1)/diags(1,l_max-1)

c           Now let tridia solve this system:
            call tridia(l_max+1,diags,b,sol)

c           Put the solution found by tridia into array u:
            do 3700 l=0,l_max
               u(j,l)=sol(l)
 3700       continue

 4000    continue

c        Find u^(n+1) also on the left and right boundaries j=0 and j=j_max:
         do 4100 l=0,l_max
            u(0,l)=uleft(l)
            u(j_max,l)=two*u(j_max-1,l)-u(j_max-2,l)
 4100    continue

c        Perform output if n+1 is an integer multiple of n_out1:
         if(mod(n+1,n_out1).eq.0)
     &    call out(j_dim,l_dim,j_max,l_max,c_x,c_y,x,y,t,n+1,u)

c        Go compute the next time plane, if any:

 5000 continue

      stop
      end