СИНТЕЗ БЛОЧНОГО АЛГОРИТМА РАЗНОСТНОГО РЕШЕНИЯ УРАВНЕНИЯ ДАЛАМБЕРА

АННОТАЦИЯ

Работа посвящена блочному «алмазному» алгоритму для решения двумерного уравнения Даламбера. Целью эксперимента является демонстрация ускорения не блочного алгоритма по сравнению с «алмазным» алгоритмом решения уравнения Даламбера в двумерном случае для прямоугольной и квадратной областей.  Благодаря учету структуры памяти используемых аппаратных средств удалось получить практически двух кратное ускорение.

Ключевые слова: уравнение Даламбера; блочные алгоритмы; ускорение вычислений.

Рассмотри двумерное волновое уравнение следующего вида [1]:

,                                        (1)

Положим вид правой части, начальные и краевые условия в прямоугольной области  следующего вида:

     (2)

Нетрудно проверить, что задача (1), (2) в области  имеет аналитическое решение:

.

Определим сетку по пространству , , где шаги по направлениям  и  определяются согласно формулам:

, .

Запишем разностную схему для уравнения (1) [2], тогда

.                                             (3) 

Рассмотрим «алмазный» алгоритм прохода по узлам сеточной области, схема которого представлены на рисунке 1.

Рисунок 1. Схема «алмазного» алгоритма

Вычисления производятся от блока под номером 1, к блоку под номером 2 и так далее к блоку под номером m слева направо.

Каждый блок данного алгоритма пересекает сразу несколько временных слоев, что позволяет уменьшить объем коммуникаций между кэш памятью и ОЗУ.

Код программы решения задачи (1) и (2) с помощью разностной схемы (3) «алмазным» алгоритмом с комментариями приведен ниже:

do t=1,Nt,n ! переход по слоям по времени    

! РАССЧЕТ ЛЕВОГО ТРЕУГОЛЬНИКА

wl=2;

do wr=n+1,3,-2 ! перемещение вверх по слоям по времени

E1(2:Ny-1,wl:wr)=2.0*E2(2:Ny-1,wl:wr)-E1(2:Ny-1,wl:wr)+ & ! рассчет нечетного слоя по времени

(E2(1:Ny-2,wl:wr)+E2(3:Ny,wl:wr)+E2(2:Ny-1,wl-1:wr-1)+E2(2:Ny-1,wl+1:wr+1)-&

4.0*E2(2:Ny-1,wl:wr))*Eps(2:Ny-1,wl:wr);

E1(Izl,2)=sin(c5*float(t+n+1-wr));  ! "жесткое" излучающее условие

wr2=wr-1; ! рассчет четного слоя по времени

 E2(2:Ny-1,wl:wr2)=2.0*E1(2:Ny-1,wl:wr2)-E2(2:Ny-1,wl:wr2)+ &

 (E1(1:Ny-2,wl:wr2)+E1(3:Ny,wl:wr2)+E1(2:Ny-1,wl-1:wr2-1)+E1(2:Ny-1,wl+1:wr2+1)- &

4.0*E1(2:Ny-1,wl:wr2))*Eps(2:Ny-1,wl:wr2);

E2(Izl,2)=sin(c5*float(t+n+1-wr2));    ! "жесткое" излучающее условие

end do

! РАССЧЕТ ЦЕНТРАЛЬНЫХ ПОДОБЛАСТЕЙ

do g=2,p-1 ! перемещение по центральным подобластям

Nl=(g-1)*n+2; Nr=g*n+1; ! определение границ нижней грани алмаза

do k=0,n-2,2 ! перемещение вверх по временным слоям

wl=Nl-k; wr=Nr-k; ! рассчет четного слоя по времени

E1(2:Ny-1,wl:wr)=2.0*E2(2:Ny-1,wl:wr)-E1(2:Ny-1,wl:wr)+ &

(E2(1:Ny-2,wl:wr)+E2(3:Ny,wl:wr)+E2(2:Ny-1,wl-1:wr-1)+E2(2:Ny-1,wl+1:wr+1)- &

4.0*E2(2:Ny-1,wl:wr))*Eps(2:Ny-1,wl:wr);

wl=wl-1; wr=wr-1; ! рассчет нечетного слоя по времени

E2(2:Ny-1,wl:wr)=2.0*E1(2:Ny-1,wl:wr)-E2(2:Ny-1,wl:wr)+ &

 (E1(1:Ny-2,wl:wr)+E1(3:Ny,wl:wr)+E1(2:Ny-1,wl-1:wr-1)+E1(2:Ny-1,wl+1:wr+1)- &

4.0*E1(2:Ny-1,wl:wr))*Eps(2:Ny-1,wl:wr);

end do

end do

! РАССЧЕТ ПРАВОГО ТРЕУГОЛЬНИКА

Nl=(p-1)*n+2; wr=Nz-1;

do wl=Nl,Nl-n+2,-2 ! перемещение вверх по слоям по времени

E1(2:Ny-1,wl:wr)=2.0*E2(2:Ny-1,wl:wr)-E1(2:Ny-1,wl:wr)+ &! рассчет нечетного слоя по времени

 (E2(1:Ny-2,wl:wr)+E2(3:Ny,wl:wr)+E2(2:Ny-1,wl-1:wr-1)+E2(2:Ny-1,wl+1:wr+1)- &

4.0*E2(2:Ny-1,wl:wr))*Eps(2:Ny-1,wl:wr);

wl2=wl-1; ! рассчет четного слоя по времени

E2(2:Ny-1,wl2:wr)=2.0*E1(2:Ny-1,wl2:wr)-E2(2:Ny-1,wl2:wr)+ &

 (E1(1:Ny-2,wl2:wr)+E1(3:Ny,wl2:wr)+E1(2:Ny-1,wl2-1:wr-1)+E1(2:Ny-1,wl2+1:wr+1)- &

4.0*E1(2:Ny-1,wl2:wr))*Eps(2:Ny-1,wl2:wr);

end do

end do

Обход внутри областей  задается внутренним циклом по k. Перемещение между блоками сеточной области задается внешним циклом по g.

Рассмотрим решение уравнения Даламбера (1) в двумерном случае для прямоугольной области при помощи «алмазного» алгоритма.

Целью эксперимента является демонстрация ускорения «алмазного» алгоритма по сравнению с не блочным решением уравнения Даламбера в двумерном случае для прямоугольной области. 

Для проведения эксперимента в качестве аппаратного обеспечения был выбран процессор, данные которого представлены в таблице 1.

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

Таблица 1.

 Данные процессора для проведения эксперимента с алмазным алгоритмом

В качестве операционной системы была выбрана ОС Ubuntu 16.04.

В качестве компилятора был выбран GFortran 5.3 - компилятора языка программирования Фортран, входящего в коллекцию компиляторов GNU.

Представленный выбор обусловлен широким распространением  линейки процессоров Intel Core i7, хорошей воспроизводимостью результатов экспериментов в операционной системе Linux и доступностью GCC.

Были выбраны следующие параметры для проведения эксперимента (Таблица 2):

Таблица 2.

 Параметры «алмазного» алгоритма для решения уравнения Даламбера в двумерном случае для прямоугольной области

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

Теоретически, в первой серии экспериментов для прямоугольной области, значение длины грани  должно быть оптимальным и иметь наименьшее время вычисления. Действительно, при указанном , размер блока в памяти компьютера равен 5,7 Мб, что не превышает объем кэш памяти выбранного процессора, указанного в Таблица 1 и в тоже время  достаточно велико для уменьшения коммуникаций между кэш памятью и ОЗУ.

При проведении численного эксперимента, время решения уравнения Даламбера в двумерном случае для прямоугольной области не блочным алгоритмом  секунд.

Отметим, что столь  малое время выбрано из соображения быстроты и надежности нахождения ускорения искомого алгоритма.

В практике оптического вычислительного эксперимента нередки длительности вычислений равные суткам и более. Так при моделировании распространения электромагнитной волны в кварцевом волноводе с указанными параметрами длительность вычислений составит почти пять часов. Время решения уравнения Даламбера в двумерном случае для прямоугольной области «алмазным» алгоритмом при изменяемом параметре n представлено в таблице 3.

Таблица 3.

Время решения уравнения Даламбера в двумерном случае для прямоугольной области при помощи «алмазного» алгоритма

График зависимости времени вычисления от длины грани блока представлен на рисунке 2.

Рисунок 2. График зависимости времени вычисления от длины грани блока

Найдем наилучшее ускорение «алмазного» алгоритма по сравнению с не блочным алгоритмом:

,

при .

Мы видим, что оптимальным значением длины грани блока, при котором время решения уравнения Даламбера в двумерном случае для прямоугольной области минимально, является значение .

При этом длительность вычислений с ростом  квадратично уменьшается, (при )  затем стабилизируется у отметки приблизительно 108-105 (при ), после чего резко линейно возрастает (при ).

График зависимости времени от блочного параметра имеет u-образный характер. Правая часть графика, представленного на рисунке 2, имеет линейную зависимость.

При длине грани  блок занимает приблизительно 2,28 Мб памяти и умещается в L3 кэш память. Объем коммуникаций минимальный и, следовательно, время вычислений минимально.

При длине грани  блок не умещается в кэш память L3, что приводит к резкому увеличению объема коммуникаций и как следствие, времени вычислений.

При  кэш память используется не оптимальным образом, что приводит к большим коммуникационным издержкам

Теоретические ожидания почти подтвердились.

Незначительное уменьшение оптимального значения блочного параметра приблизительно в 2,5 раза  связано с тем, что кроме хранения блока в кэш памяти операционная система занята другими системными процессами, часть из которых имеет доступ к кэшу и занимает там место.

Рассмотрим решение уравнения Даламбера в двумерном случае для квадратной области при помощи «алмазного» алгоритма.

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

Так, целью следующего эксперимента является демонстрация ускорения «алмазного» алгоритма по сравнению с не блочным решением уравнения Даламбера в двумерном случае для квадратной области.  

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

Были выбраны следующие параметры для проведения эксперимента (Таблица 4):

Таблица 4.

 Параметры «алмазного» алгоритма для решения уравнения Даламбера в двумерном случае для квадратной области

Заданная сеточная область отличается от сеточной области предыдущего эксперимента по геометрическим параметрам. Однако совпадает с ней по количеству узлов сеточной области.

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

Теоретически, в первой серии экспериментов для квадратной области, значение длины грани  должно быть оптимальным и иметь наименьшее время вычисления.

При проведении численного эксперимента, время решения уравнения Даламбера в двумерном случае для квадратной области блочным алгоритмом  секунд.                                                                                               

Время решения уравнения Даламбера в двумерном случае для квадратной области «алмазным» алгоритмом при изменяемом параметре n представлено в таблице 5.

Таблица 5.

Время решения уравнения Даламбера в двумерном случае для квадратной области при помощи «алмазного» алгоритма

График зависимости времени вычисления от длины грани блока представлен на рисунке 3.

Рисунок 3. График зависимости времени вычисления от длины грани блока

Найдем наилучшее ускорение «алмазного» алгоритма по сравнению с не блочным алгоритмом:

,

при .

Мы видим, что оптимальным значением длины грани блока, при котором время решения уравнения Даламбера в двумерном случае для прямоугольной области минимально, является значение .

Аналогично предыдущему случаю, график зависимости времени от блочного параметра имеет u-образный характер. Однако правые части графиков, представленных на рисунке 2 и 3, отличаются.  Если для предыдущего эксперимента правая часть имеет линейную зависимость, то здесь она по виду соответствует функции  . Кроме того, «дно у ковша» менее пологое.

Так как при длине грани  блок занимает приблизительно 5,86 Мб памяти и умещается в L3 кэш память.

Отмеченные выше особенности объясняются резко возросшей шириной сеточной области.

Если ранее увеличение блочного параметра допустим на  приводило к увеличению объема блока на , то теперь те же модификации вызывают увеличение объема на величину . Именно поэтому график оказался менее пологим.

По сравнению с прямоугольной областью, время не блочного алгоритма и время «алмазного» алгоритма при оптимальном значении  приблизительно одинаковое. Разница составляет примерно 10%.

Теоретические ожидания почти подтвердились.

Незначительное уменьшение оптимального значения блочного параметра приблизительно на 20% связано с тем, что кроме хранения блока в кэш памяти, операционная система занята другими системными процессами, часть из которых имеет доступ к кэшу и занимает там место.

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

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

Следовательно, классическое понятие вычислительной сложности как количество арифметических операций не может в полной мере характеризовать численный метод.

Список литературы:

1. Плохотников К.Э. Вычислительные методы. Теория и практика в среде MATLAB: курс лекций / К.Э. Плохотников. – 2-е изд. – М.: МГУ, 2015. – 496 с.
2. Неганов В.А. Линейная макроскопическая электродинамика / В.А. Неганов, С.Б. Раевский, Г.П. Яровой // Т. 1, – М.: Радио и связь, 2000. – 509 с.