fortranでスレッド並列計算

ompdiffusioneq1d.f90

program ompdiffusioneq1d

    !$ use omp_lib

    implicit none

    !$ integer :: nthreads
    !$ integer :: mythreadno
    integer, allocatable :: bc(:)
    integer :: division_x
    integer :: n
    integer :: n_end
    integer :: n_interval
    integer :: i
    integer :: i1, i2
    integer :: cell

    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
    !$ real(8) :: sum_error

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

    !------------------------------------------------------------------
    pi = 4.0D0 * datan(1.0D0)
    !------------------------------------------------------------------

    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)

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

    write(6, '(A, E15.8, A, A, E15.8)') 'delta_x = ', delta_x, ',', &
            'delta_t = ', delta_t
    !------------------------------------------------------------------

    allocate(x(division_x + 1))
    allocate(u(division_x + 1))
    allocate(u_backup(division_x + 1))
    allocate(du(division_x + 1))
    allocate(u_exact(division_x + 1))
    allocate(bc(division_x + 1))

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

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

    ! initial codition
    do i = 1, division_x + 1
        u(i) = dsin(pi * x(i)) ** 3
    end do

    ! boundary condition
    u(1) = 0.0D0
    bc(1) = 1

    u(division_x + 1) = 0.0D0
    bc(division_x + 1) = 1

    !$omp parallel &
    !$omp private(nthreads, mythreadno)
    !$ nthreads = omp_get_num_threads()
    !$ mythreadno = omp_get_thread_num()
    !$ write(6, '(A, I6, A, A, I6)') 'mythreadno: ', mythreadno, ',' , 'nthreads: ', nthreads
    !$omp end parallel

    !$ if(nthreads .eq. 1) then

    call cpu_time(totaltime_start)

    !$ else
    !$ totaltime_start = omp_get_wtime()
    !$ end if

    n = 0
    t = 0.0D0

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

        !$ if(nthreads .eq. 1) then
        call cpu_time(time_start)
        !$ else
        !$ time_start = omp_get_wtime()
        !$ end if

        !$omp parallel do &
        !$omp private(cell, i1, i2)

        do cell = 1, division_x
            i1 = cell
            i2 = cell + 1
            du(i1) = du(i1) + (-u(i1) + u(i2)) / (delta_x * delta_x)
            du(i2) = du(i2) + (u(i1) - u(i2)) / (delta_x * delta_x)
        end do

        !$omp end parallel do

        !$omp parallel do &
        !$omp private(i)

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

        !$omp end parallel do

        !$ if (nthreads .eq. 1) then
        call cpu_time(time_end)
        !$ else
        !$ time_end = omp_get_wtime()
        !$ end if

        n = n + 1
        t = t + delta_t

        if(mod(n, n_interval) .eq. 0) then
            write(90, '(7X, A, 11X, A, 11X, A, 8X, A)') 'i,', 'x(i), ', 'u(i), ', 'u_exact'
            error = 0.0D0

            !$omp parallel do &
            !$omp reduction(+: sum_error) &
            !$omp private(i, error)

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

                error = error + (u(i) - u_exact(i)) ** 2
                !$ sum_error = error
            end do

            !$omp end parallel do
            !$ error = sum_error
            error = dsqrt(error) / dfloat(division_x + 1)

            write(char_n, *) n
            char_n = adjustl(char_n)
            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

            do i = 1, division_x + 1
                write(90, '(I8, 3(A, E15.8))') i, ', ', x(i), ', ', u(i), ', ', u_exact(i)
            end do
            write(90, '(A, E15.8)') 'Error=', error
            close(90)

        end if
    end do

    !$ if  (nthreads .eq. 1) then
    call cpu_time(totaltime_end)
    !$ else
    !$ totaltime_end = omp_get_wtime()
    !$ end if

    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)

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

    stop

end program ompdiffusioneq1d