3
\$\begingroup\$

I'm a newbie playing around with Fortran 90 and openmp and wrote the code below (a simple 2D heat transfer simulation) for testing purposes.

So far I don't see any speedup by using openmp / parallel processing. In fact, it gets slower when I use more procs. I understand that the parallelization requires some extra overhead, but my T array should be big enough (2500 x 1500) to make a parallelization "worth it". What am I doing wrong? What can be improved?

I am compiling with: gfortran -fopenmp heat.f90 -o heat

and setting the number of threads in my shell with: export OMP_NUM_THREADS=4

heat.f90 is below:

program heat

    use omp_lib
    implicit none

    integer, parameter :: nx = 2500
    integer, parameter :: ny = 1500
    real, parameter :: H = 1.
    real, parameter :: alpha = 1.e-4
    real, parameter :: sigma = 0.2
    real, parameter :: dt = 0.001
    real, parameter :: physical_time = 50.

    real x, y
    integer i, j
    real aspect, L, dx, dy

    real, dimension (:,:), allocatable :: T, Tn
    real d2Tdx2, d2Tdy2

    ! -----

    integer ts, total_ts, ts_left, frames_total, output_period_ts, pic_counter
    real t1, t2, time_elapsed, run_rate, progress_pct
    real time_left_s, time_left_h, time_left_d, time_total_s, time_total_h, time_total_d
    character(len=16) :: pic_counter_str

    ! -----

    aspect = real(nx)/real(ny)
    L = H*aspect
    dx = L/real(nx)
    dy = H/real(ny)

    ! ----- initialize T

    allocate (T(nx,ny))
    allocate (Tn(nx,ny))

    do i=1, nx
        do j=1, ny
            x = real(i)/real(nx)*L
            y = real(j)/real(ny)*H
            T(i,j) = 10. + 2. * ( &
                                  1./(2.*3.14*sigma**2.) * &
                                  exp(-1. * ( (x-L/2.)**2. + (y-H/2.)**2. ) / (2.*sigma**2.)) &
                                )
        enddo
    enddo

    ! ----- report

    total_ts = ceiling(physical_time/dt)
    frames_total = 30*20
    output_period_ts = int(total_ts/frames_total)

    print '(A, ES10.4)', 'dt ................... ', dt
    print '(A, ES10.4)', 'physical time ........ ', physical_time
    print '(A, I2)',     'num procs ............ ', omp_get_num_procs()
    print '(A, I2)',     'num threads .......... ', omp_get_max_threads()
    print '(A, I6)',     'timesteps ............ ', total_ts
    print '(A, I4)',     'output period [ts] ... ', output_period_ts
    print '(A, I6)',     'num frames ........... ', frames_total
    print *, ' '
    print *, 'ts, %, time left [h], time total [h]'
    print *, '------------------------------------'

    ! ----- evolve T field

    call cpu_time (t1)
    pic_counter = 0

    do ts = 0, total_ts

        if (mod(ts,output_period_ts) .eq. 0) then

            ! ----- plot (optional)

            !open(3, file="T.dat", access="stream")
            !write(3) T(:,:)
            !close(3)

            !write (pic_counter_str,'(I5.5)') pic_counter
            !call system('gnuplot -c plot3d.gnu '//trim(pic_counter_str))
            !pic_counter = pic_counter + 1

            ! ----- quick check if everythings alright

            if (isnan(T(100,100))) then
                print *, "divergence! decrease ts."
            end if

            if (isnan(T(100,100))) call EXIT(1)

            ! ----- report run speed stats

            if ( ts .gt. 0) then
                call cpu_time (t2)

                time_elapsed = t2-t1
                ts_left = total_ts-ts
                run_rate = real(ts)/time_elapsed

                time_left_s = real(ts_left)/run_rate
                time_left_h = time_left_s/3600.
                time_left_d = time_left_h/24.

                time_total_s = real(total_ts)/run_rate
                time_total_h = time_total_s/3600.
                time_total_d = time_total_h/24.

                progress_pct = 100.*real(ts)/real(total_ts)

                print '(I8, F7.2, A, F8.2, F8.2)', ts, progress_pct, "%", time_left_h, time_total_h
           end if

        end if

        ! ----- energy eqn

        Tn(:,:) = T(:,:)

        !$omp parallel shared(T) private(i,j,d2Tdx2,d2Tdy2)
        !$omp do
        do i=2, nx-1
            do j=2, ny-1            
                d2Tdx2 = ( Tn(i+1,j) - 2.*Tn(i,j) + Tn(i-1,j) ) / dx**2.
                d2Tdy2 = ( Tn(i,j+1) - 2.*Tn(i,j) + Tn(i,j-1) ) / dy**2.
                T(i,j) =  Tn(i,j) + dt*(alpha*(d2Tdx2 + d2Tdy2))
            end do
        end do
        !$omp end do                     
        !$omp end parallel

        !! adiabatic BCs
        T(:, 1)  = T(:, 2)
        T(:, ny) = T(:, ny-1)
        T(1, :)  = T(2, :)
        T(nx, :) = T(nx-1, :)

        ! -----

    end do

end program heat

... the (optional) gnuplot (5.0) plot script plot3d.gnu is below, uncomment the 'plot' portion of the code above for image output if you'd like

#!/usr/bin/gnuplot --persist

pic_counter=ARG1
png_filename = sprintf("%s%s%s", "T_", pic_counter, ".png")

set terminal pngcairo background rgb "white" enhanced font "arial,16" fontscale 1.0 size 1280, 720

set key top
set border 1+2+4+8+16+32+64+256+512

set xrange [0:2500]
set yrange [0:1500]
set zrange [10:20]
set cbrange [10:18]
set view equal xy

set style data pm3d
set pm3d implicit at s hidden3d depth noborder
set pm3d lighting primary 0.5 specular 0.4

set view 57, 20, 2.05, 0.55
#show view

set colorbox horizontal user origin 0.05, 0.05 size 0.35, 0.05

set output png_filename
splot "T.dat" binary array=2500x1500 format="%f" notitle
```
\$\endgroup\$
1
\$\begingroup\$

One problem is that you're using d2Tdx2 and d2Tdy2 inside your omp loop but they're not listed as private in the parallel directive. This will cause all the threads to use the same variables with a big performance hit, and possibly errors in the calculations.

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ Great catch! This came from some last minute edits before posting on SE. Unfortunately fixing this does not help with the parallel performance. I've updated the question accordingly. \$\endgroup\$ – HotDogCannon Nov 24 '19 at 7:32
  • \$\begingroup\$ Assuming "cpu_time()" measures what the name says it does, then you should expect it to get worse when you add parallelism (due to some overhead in creating threads and so on). You need to measure elapsed (wall-clock) time to see a speedup, for which omp_get_wtime() is a simple and useful function. If you are measuring performance, you should also compile with optimization enabled (e.g. -O3). \$\endgroup\$ – Jim Cownie Nov 25 '19 at 9:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.