program diff
      integer xmax,ymax
c     Размер сетки по оси X
      parameter (xmax=4000)
c     Размер сетки по оси Y
      parameter (ymax=6000)
      real*8 f0(xmax,ymax), f1(xmax,ymax)
      real*8 df
      integer x,y,n
c     Задаем начальные значения
      call initial(f0,xmax,ymax)
      call initial(f1,xmax,ymax)
c     Задаем граничные значения
      call boundary(f0,xmax,ymax)
      call boundary(f1,xmax,ymax)
c     Начинаем основной итерационный цикл
      df=0.0
      n=0
    1 continue
c     При смене итераций меняем местами массивы
      if ( ((n/2)*2) .EQ. n ) then
c       Обрабатываем f0 для каждой четной итерации
        call iter(f0,f1,xmax,ymax,df)
      else
c       Обрабатываем f1 для каждой нечетной итерации
        call iter(f1,f0,xmax,ymax,df)
      endif
c     Увеличиваем номер цикла на единицу
      n=n+1
c     Выводим на экран разницу между итерациями
      write(*,*) 'Diff:',df
c     Если разница больше заданной, то решение не найдено
c     и мы идем на следующий цикл
      if (df > 1e-2) goto 1
c     Печатаем количество итераций, потребовавшихся
c     для нахождения решения
      write(*,*) 'Iteration count:', n
      stop
      end
 
c     Подпрограмма вычисления искомой функции
      subroutine iter(f0,f1,xmax,ymax,df)
      integer xmax,ymax
      real*8 f0(xmax,ymax), f1(xmax,ymax)
      real*8 dt,dx,dy
      real*8 df,dff,df1
      integer x,y,n
      dt=0.01
      dx=0.5
      dy=0.5
      df1=0.0
c     Вычисляем функцию в ячейках сетки
      do y=2,ymax-1
        do x=2,xmax-1
        dff=dt*(
     x    (f0(x+1,y)-2*f0(x,y)+f0(x-1,y))/(dx*dx)
     x    +
     x    (f0(x,y+1)-2*f0(x,y)+f0(x,y-1))/(dy*dy)
     x    )
        f1(x,y)=f0(x,y)+dff
c     Находм максимальную дельту
        if( df1 < abs(dff) ) df1=abs(dff)
        end do
      end do
      df=df1
      return
      end
 
c     Подпрограмма задания начальных значений
      subroutine initial(f0,xmax,ymax)
      integer xmax,ymax,x,y
      real*8 f0(xmax,ymax)
c     Задаем начальные значения массива
      do x=1,xmax
        do y=1,ymax
          f0(x,y)=0.0
        end do
      end do
      return
      end
 
c     Подпрограмма задания граничных условий
      subroutine boundary(f0,xmax,ymax)
      integer rank,xmax,ymax,x,y
      real*8 f0(xmax,ymax)
      common myid,np
c     Задаем граничные значения массива
      do y=1,ymax
        f0(1,y)=0.0
        f0(xmax,y)=0.0
      end do
      do x=1,xmax
        f0(x,1)=sin((x*1.0)/(xmax/2))
	f0(x,ymax)=0.0
      end do
      f0(xmax/2,ymax)=-5.0
      return
      end
 
Загрузить исходный код программы

Как видите, программа достаточно проста. Здесь я ограничусь всего несколькими пояснениями. В качестве разностной сетки мы берем прямоугольник со сторонами 4000 ячеек по оси X и 6000 ячеек по оси Y. Размерность сетки задается в строка 4 и 6.

Далее, в строках 11 - 15 мы задаем начальные и граничные условия посредством вызова подпрограмм initial (строка 69) и boundary (строка 82) соответственно, заполняющих по некоторому правилу ячейки сетки.

Цикл с 19 по 34 строчку является основным итерационным процессом. В строке 34 принимается решение о продолжении итераций на основе максимальной разницы между значениями искомой функции на текущем и предыдущем слое по времени. Если разница больше некоторого заранее заданного значения, то итерации продолжаются. Само вычисление функции в ячейках оформлено в виде подпрограммы iter, код которой начинается со строки 42. Эта же подпрограмма вычисляет максимальное значение Δf.

В самом начале в качестве исходных данных берется массив f0, а вычисленные значения функции на следующем шаге по времени заносятся в массив f1 - строка 23. На следующей итерации за основу берется уже массив f1, то есть массивы меняются местами - строка 26. И так далее.

Теперь о вычислении самой функции. Выполнение этого действия выделено в подпрограмму iter (строка 42). Функция на k+1 слое по времени вычисляется на основе значений этой функции на k-ом слое. Делается это в два этапа. Сначала вычисляются значения внутренних ячеек на основе разностной формулы (строки 55-60), а затем находим максимальную разницу значений функции между слоями (строка 61).

В следующем разделе мы с вами посмотрим, как можно эту программу превратить в параллельную, что для этого надо и какие изменения в исходный код должны быть внесены.