fortranでプロセス並列計算

mpidiffusioneq1d.f90

program mpidiffusioneq1d

    implicit none

    include 'mpif.h'

    integer :: nprocs
    integer :: myrank
    integer :: ierr
    integer :: tag
    integer :: request1, request2
    integer :: status
    integer, allocatable :: count(:)
    integer, allocatable :: bc(:)
    integer :: division_x
    integer :: m
    integer :: n
    integer :: n_end, n_end_backup
    integer :: n_interval, n_interval_backup
    integer :: i
    integer :: i1, i2
    integer :: j
    integer :: j1, j2
    integer :: i_start, i_end
    integer :: cell_start, cell_end
    integer :: cell
    integer :: proc

    real(8) :: sbuffer1, sbuffer2
    real(8) :: rbuffer1, rbuffer2
    real(8) :: time_start, time_end
    real(8) :: totaltime_start, totaltime_end
    real(8), allocatable :: x(:)
    real(8) :: delta_x
    real(8), allocatable :: u(:)
    real(8), allocatable :: u_backup(:)
    real(8), allocatable :: du(:)
    real(8), allocatable :: u_exact(:)
    real(8) :: alpha
    real(8) :: t, delta_t
    real(8) :: d
    real(8) :: pi
    real(8) :: error, sum_error

    character(100) :: char_n
    character(1) :: dataname

    !------------------------------------------------------------------
    call mpi_init(ierr)
    call mpi_comm_size(mpi_comm_world, nprocs, ierr)
    call mpi_comm_rank(mpi_comm_world, myrank, ierr)
    !------------------------------------------------------------------
    pi = 4.0D0 * datan(1.0D0)
    !------------------------------------------------------------------

    if (myrank == 0) then
        open(10, file = 'input_diffusioneq1d.dat')

        read(10, *) dataname
        read(10, *) division_x

        read(10, *) dataname
        read(10, *) alpha

        read(10, *) dataname
        read(10, *) n_end

        read(10, *) dataname
        read(10, *) n_interval

        read(10, *) dataname
        read(10, *) d

        close(10)
    end if

    call mpi_bcast(division_x, 1, mpi_integer, 0, mpi_comm_world, ierr)
    call mpi_bcast(alpha, 1, mpi_double_precision, 0, mpi_comm_world, ierr)
    call mpi_bcast(n_end, 1, mpi_integer, 0, mpi_comm_world, ierr)
    call mpi_bcast(n_interval, 1, mpi_integer, 0, mpi_comm_world, ierr)
    call mpi_bcast(d, 1, mpi_double_precision, 0, mpi_comm_world, ierr)

    n_interval_backup = n_interval
    n_end_backup = n_end

    !------------------------------------------------------------------
    delta_x = (1.0D0 - 0.0D0) / dfloat(division_x)
    delta_t = d * delta_x * delta_x / alpha

    if (myrank == 0) then
        write(6, '(A, E15.8, A, A, E15.8)') 'delta_x = ', delta_x, ',', &
                'delta_t = ', delta_t
    end if
    !------------------------------------------------------------------
    m = int(dfloat(division_x) / dfloat(nprocs))
    !------------------------------------------------------------------

    if (myrank /= nprocs - 1) then
        allocate(x(m + 1))
        allocate(u(m + 1))
        allocate(u_backup(m + 1))
        allocate(du(m + 1))
        allocate(u_exact(m + 1))
        allocate(bc(m + 1))
        allocate(count(m + 1))
    else
        allocate(x(division_x - myrank * m + 1))
        allocate(u(division_x - myrank * m + 1))
        allocate(u_backup(division_x - myrank * m + 1))
        allocate(du(division_x - myrank * m + 1))
        allocate(u_exact(division_x - myrank * m + 1))
        allocate(bc(division_x - myrank * m + 1))
        allocate(count(division_x - myrank * m + 1))
    end if

    x = 0.0D0
    u = 0.0D0
    u_backup = 0.0D0
    du = 0.0D0
    u_exact = 0.0D0
    bc = 0
    count = 1

    if (myrank /= nprocs - 1) then
        i_start = myrank * m + 1
        i_end = (myrank + 1) * m

        cell_start = myrank * m + 1
        cell_end = (myrank + 1) * m
    else
        i_start = myrank * m + 1
        i_end = division_x + 1

        cell_start = myrank * m + 1
        cell_end = division_x
    end if

    count(1) = 2
    count(i_end - i_start + 1) = 2

    !------------------------------------------------------------------
    ! x-coordinate
    do i = i_start, i_end
        j = i - i_start + 1
        x(j) = 0.0D0 + delta_x * dfloat(i - 1)
    end do
    !------------------------------------------------------------------

    ! initial codition
    do i = i_start, i_end
        j = i - i_start + 1
        u(j) = dsin(pi * x(j)) ** 3
    end do

    ! boundary condition

    if (myrank == 0) then
        u(1) = 0.0D0
        bc(1) = 1
        count(1) = 1
    end if

    if (myrank == nprocs - 1) then
        u(i_end - i_start + 1) = 0.0D0
        bc(i_end - i_start + 1) = 1
        count(i_end - i_start + 1) = 1
    end if

    !------------------------------------------------------------------

    if (myrank == 0) then
        totaltime_start = mpi_wtime()
    end if

    !------------------------------------------------------------------
    n = 0
    t = 0.0D0

    do while(n < n_end)
        u_backup = u
        du = 0.0D0

        if (myrank == 0) then
            time_start = mpi_wtime()
        end if

        do cell = cell_start, cell_end
            i1 = cell
            i2 = cell + 1

            j1 = i1 - i_start + 1
            j2 = i2 - i_start + 1

            du(j1) = du(j1) + (-u(j1) + u(j2)) / (delta_x * delta_x)
            du(j2) = du(j2) + (u(j1) - u(j2)) / (delta_x * delta_x)
        end do

        tag = 1

        if (myrank /= 0) then
            sbuffer1 = du(1)
            call mpi_isend(sbuffer1, 1, mpi_double_precision, myrank - 1, tag, mpi_comm_world, request1, ierr)
            call mpi_irecv(rbuffer1, 1, mpi_double_precision, myrank - 1, tag, mpi_comm_world, request1, ierr)
            call mpi_wait(request1, status, ierr)

            du(1) = du(1) + rbuffer1
        end if

        if (myrank /= nprocs - 1) then
            sbuffer2 = du(i_end - i_start + 1)
            call mpi_isend(sbuffer2, 1, mpi_double_precision, myrank + 1, tag, mpi_comm_world, request2, ierr)
            call mpi_irecv(rbuffer2, 1, mpi_double_precision, myrank + 1, tag, mpi_comm_world, request2, ierr)
            call mpi_wait(request2, status, ierr)

            du(i_end - i_start + 1) = du(i_end - i_start + 1) + rbuffer2
        end if

        do i = i_start, i_end
            j = i - i_start + 1
            u(j) = u_backup(j) + dfloat(1 - bc(j)) * delta_t * alpha * du(j)
        end do

        if (myrank == 0) then
            time_end = mpi_wtime()
        end if

        n = n + 1
        t = t + delta_t

        n_end = n_end_backup
        n_interval = n_interval_backup

        if(mod(n, n_interval) == 0) then
            error = 0.0D0

            do i = i_start, i_end
                j = i - i_start + 1
                u_exact(j) = (3.0D0 / 4.0D0) * dsin(pi * x(j)) * dexp(-alpha * pi ** 2 * t) &
                        - (1.0D0 / 4.0D0) * dsin(3.0D0 * pi * x(j)) * dexp(-alpha * (3.0D0 * pi) ** 2 * t)

                error = error + (u(j) - u_exact(j)) ** 2 / dfloat(count(j))

            end do

            call mpi_allreduce(error, sum_error, 1, mpi_double_precision, mpi_sum, mpi_comm_world, ierr)
            error = dsqrt(sum_error) / dfloat(division_x + 1)

            do proc = 1, nprocs
                if(myrank == proc - 1) then
                    write(char_n, *) n
                    char_n = adjustl(char_n)

                    if (myrank == 0) then
                        open(90, file = 'output_diffusioneq1d_n=' // trim(char_n) // '.dat')
                        write(90, '(A, I8, A, A, E15.8)') '! n= ', n, ', ', 't=', t
                        write(90, '(A, E15.8, A, A, E15.8)') '! delta_x= ', delta_x, ', ', 'delta_t=', delta_t
                        write(90, '(7X, A, 11X, A, 11X, A, 8X, A)') 'i,', 'x(i), ', 'u(i), ', 'u_exact(i)'
                    else
                        open(90, file = 'output_diffusioneq1d_n=' // trim(char_n) // '.dat', access = 'append')
                    end if

                    do i = i_start, i_end
                        j = i - i_start + 1
                        write(90, '(I8, 3(A, E15.8))') j, ', ', x(j), ', ', u(j), ', ', u_exact(j)
                    end do
                    close(90)
                end if
                call mpi_barrier(mpi_comm_world, ierr)
            end do

            if (myrank == 0) then
                open(90, file = 'output_diffusioneq1d_n=' // trim(char_n) // '.dat', access = 'append')
                write(90, '(A, E15.8)') 'Error=', error
                close(90)
            end if
        end if
    end do

    if (myrank == 0) then
        totaltime_end = mpi_wtime()
    end if

    if (myrank == 0) then
        open(90, file = 'output_diffusioneq1d_exetime.dat')
        write(90, '(A, E15.8, 1X, A)') 'Execution time: ', time_end - time_start, 'sec'
        write(90, '(A, E15.8, 1X, A)') 'Total execution time: ', totaltime_end - totaltime_start, 'sec'
        close(90)
    end if

    deallocate(x)
    deallocate(u)
    deallocate(u_backup)
    deallocate(u_exact)
    deallocate(du)
    deallocate(bc)
    deallocate(count)

    call mpi_finalize(ierr)

    stop

end program mpidiffusioneq1d