ACSI Update: 2D Diffusion Solver

In an effort learn the Fortran programming language and implement numerical solutions to PDEs, I have written a program to solve 2D diffusion problems. The user must define a computational domain for the solution and specify boundary and initial conditions. This is a major step towards implementing a fully functional CFD scheme.

      PROGRAM DiffusionSolver2

      !Numerically solves the PDE: Ut = alpha^2 * (Uxx + Uyy)

      IMPLICIT NONE

      DOUBLE PRECISION U(101, 101, 1001), X(101),

     .Y(101), T(1001)

      DOUBLE PRECISION alpha, Tmax, Xmax, Ymax, dx,

     .dy, dt, Uxx, Uyy, Rx, Ry,

     .top, bottom, left, right, initial

      INTEGER i, j, k, Imax, Jmax, Kmax

    

      alpha    = .1        !Choose an alpha

      Xmax    = 1.0        !Define the x domain of the solution

      Ymax    = 1.0        !Define the y domain of the solution

      Tmax    = 1.0        !Define the t domain of the solution

      Imax    = 20        !Specify the number of x finite elements < 101

      Jmax    = 20        !Specify the number of y finite elements < 101

      Kmax    = 50        !Specify the number of t finite elements < 1001

     

      dx = Xmax / (Imax-1)

      dy = Ymax / (Jmax-1)

      dt = Tmax / (Kmax-1)

      Rx = dt * alpha**2 / dx**2

      Ry = dt * alpha**2 / dy**2

     

      !setting up the grids

      DO i=1, Imax

          X(i) = (i - 1) * dx

      END DO

      DO j=1, Jmax

          Y(j) = (j - 1) * dy

      END DO

      DO k=1, Kmax

          T(k) = (k - 1) * dt

      END DO

     

      !initial conditions

      DO i=1, Imax

          DO j=1, Jmax

              U(i, j, 1) = initial(X(i), Y(j))

          END DO

      END DO

     

      !main time loop

      DO k = 1, Kmax

      !Boundary Conditions:

          !bottom

          DO i=1, Imax-1

              U(i, 1, k+1)    = bottom(X(i), T(k+1))

          END DO

          !top

          DO i=2, Imax

              U(i, Jmax, k+1)    = top(X(i), T(k+1))

          END DO

          !right

          DO j=1, Jmax-1

              U(Imax, j, k+1)    = right(Y(j), T(k+1))

          END DO

          !left

          DO j=2, Jmax

              U(1, j, k+1)    = left(Y(j), T(k+1))

          END DO

      !Middle Diffusion Scheme:

          DO i=2, Imax-1

              DO j=2, Jmax-1

              Uxx = U(i+1, j, k) - 2.0 * U(i, j, k) + U(i-1, j, k)

              Uyy = U(i, j+1, k) - 2.0 * U(i, j, k) + U(i, j-1, k)

              U(i, j, k+1) = U(i, j, k) + Rx * Uxx + Ry * Uyy

              END DO

          END DO

      END DO

     

      !Print to file

      OPEN(1, file = 'DiffusionSolver2D_Output.dat')

      WRITE(1, *) 't, x, y, u'

      DO k=1, Kmax

          DO i=1, Imax

              DO j=1, Jmax

                  WRITE(1, *) T(k), ', ', X(i), ', ', Y(j), ', ',

     .            U(i, j, k)

              END DO

          END DO

      END DO

      CLOSE(1)

     

      STOP

    END !PROGRAM DiffusionSolver2

   

    FUNCTION top(x, t)

          IMPLICIT NONE

          DOUBLE PRECISION top, x, t

          top = 0.0 !Write any function of x and t

          RETURN

      END

      FUNCTION bottom(x, t)

          IMPLICIT NONE

          DOUBLE PRECISION bottom, x, t

          bottom = 0.0 !Write any function of x and t

          RETURN

      END

      FUNCTION left(y, t)

          IMPLICIT NONE

          DOUBLE PRECISION left, y, t

          left = 0.0 !Write any function of y and t

          RETURN

      END

      FUNCTION right(y, t)

          IMPLICIT NONE

          DOUBLE PRECISION right, y, t

          right = 0.0 !Write any function of y and t

          RETURN

      END

      FUNCTION initial(x, y)

          IMPLICIT NONE

          DOUBLE PRECISION initial, x, y

          initial = 5 !Write any function of x and y

          RETURN

      END

The language is Fortran 95 and I am using GNU's gfortran compiler and Jedit to write the code