К вопросу преобразования и параллельного ввода граничных условий при решении краевых задач в едином вычислительном потоке Текст научной статьи по специальности «Математика»

УДК 681.3

Ю.Я. ЛЕДЯНКИН

К ВОПРОСУ ПРЕОБРАЗОВАНИЯ И ПАРАЛЛЕЛЬНОГО ВВОДА ГРАНИЧНЫХ УСЛОВИЙ ПРИ РЕШЕНИИ КРАЕВЫХ ЗАДАЧ В ЕДИНОМ ВЫЧИСЛИТЕЛЬНОМ ПОТОКЕ

Анотація. Розвивається ідея побудови єдиного обчислювального потоку при вирішенні задач математичної фізики. Запропонований метод передбачає паралельне введення у вигляді сплайнів крайових умов, обчислення коефіцієнтів вузлових функцій, рішення системи Ax = b з обробкою даних на рівні складних структур у спецпроцесорі зі скалярним помножувачем у складі процесорного елемента.

Ключові слова: сплайн, скалярний помножувач, задачі математичної фізики, крайові умови, єдиний обчислювальний потік, паралельне представлення та обробка складних структур даних.

Аннотация. Развивается идея построения единого вычислительного потока при решении задач математической физики. Предложенный метод предполагает параллельный ввод в виде сплайнов краевых условий, вычисление коэффициентов узловых функций, решение системы Ax = b с обработкой данных на уровне сложных структур в спецпроцессоре со скалярным умножителем в составе процессорного элемента.

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

Abstract. The idea of creating a unified computational flow solving the problems of mathematical physics is developed. A proposed method involves parallel input using the form of splines of boundary conditions, the coefficients calculation of the nodal functions, the solution for Ax = b with data processing at the level of complex structures in the specially designed processor with scalar multiplier in the processor element.

Keywords: spline, scalar multiplier, problems of mathematical physics, boundary conditions, single computational flow, parallel processing and presentation of complex data structures.

1. Введение

Идея параллельных вычислений и построения спецпроцессоров (СП) не нова [1], но актуальна и сейчас. Не является новостью и построение спецпроцессора для решения краевых задач, описываемых дифференциальными уравнениями в частных производных, что свойственно задачам математической физики (МФ). Разнообразие задач и методов их решения, способов аппроксимации краевых (начальных (НУ) и граничных (ГУ)) условий резко усложняет создание спецпроцессора, который должен их решать. Тем более решать на параллельной структуре и в едином вычислительном потоке (ЕВП): ввод и преобразование ГУ; расчет коэффициентов в узлах сетки; создание и расчет глобальной системы линейных алгебраических уравнений (СЛАУ) или коэффициентов узловых функций конечных элементов (КЭ) и глобальной матрицы жесткости с вектором нагрузок применительно к методу конечных элементов (МКЭ); решение СЛАУ (нахождение узловых значений) и обработка результатов счета. Задача сложная, важная, но необходимая.

Единая математическая линия и соответствующая ей структура, которая реализует такой вычислительный процесс в виде [1] единого технологического потока (ЕТП), предполагает высокую производительность СП, а массовое тиражирование элементов и узлов средств вычислительной техники (ВТ) - снижение стоимости изделия в целом. Также важно чтобы все базировалось на реализации простейших операций типа “+”, “-“, “ * ”, выполнение которых осуществляется с помощью элементов, освоенных в практике проекти-

© Ледянкин Ю.Я., 2012

ISSN 1028-9763. Математичні машини і системи, 2012, № 1

рования и производства и позволяющих выполнять требуемые операции в параллельном режиме.

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

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

2. Постановка задачи

В работе [1] предложен и описан подобный подход с примерами решения задач МФ, описываемых уравнениями вида

Эu / Эt = div(k * gradu), (1)

p

Lu = ^ Э / Э х( к Э u / Э х) (2)

а=1

при х е G : U(х,0) = u0(х), ха е Г; Эи/Эп + а(ы - т(х) = 0, к = k(х, t). (3)

Он предлагает применение МКЭ и реализацию всех вычислений на спецпроцессоре, состоящем из процессорных элементов (ПЭ), а каждый из ПЭ включает скалярный умножитель (СУ).

В [1] на примере уравнения Пуассона

Э2 u/Э х2+ Э2 u/Э у2 -ф = 0 (4)

<р = 0 при х = 0, (р = 1 при х = 1 применительно к методу Галеркина предложен и описан способ получения матриц жесткости КЭ, а также глобальной матрицы жесткости и вектора нагрузок путем аппроксимации сеточных функций с помощью МКЭ. Предложенная процедура, кроме расчета коэффициентов узловых функций и решения полученных систем алгебраических уравнений, включает и численное интегрирование в ПЭ, дифференцирование и интегрирование сплайнов, которыми аппроксимируют сеточные функции. Все процедуры выполняются на параллельной структуре с обработкой информации на уровне ССД.

Каждое из этих направлений - параллельная реализация, обработка ССД и ЕТП, включающее и расчет, и обработку коэффициентов (решение СЛАУ), обеспечивает в отдельности выигрыш в производительности СП, а совмещение их в одном устройстве - это дополнительные преимущества. Очень важно, что все этапы по обработке информации выполняют, а полученные результаты хранят в памяти КЭ (то есть в каждом ПЭ), сокращая до минимума пересылки данных. А пересылки, кроме временных затрат, это и дополнительные затраты по организации линий связи в пределах площади ЧИПа, это и задержки электрических сигналов в линиях связи, и возможные ошибки при передаче информации, а

линии связи, подключенные к наружным выводам, могут иметь и ограничения на количество наружных выводов в ЧИПе.

3. Краевые условия задач МФ

Ниже предлагается рассмотреть вопросы, связанные с вводом краевых условий в СП параллельного типа при решении задач МФ.

Триангуляция области решения G и границ Г - это процедуры, с которых начинается процесс решения задачи МФ. В [2, 3] исследованы вопросы, связанные с различными методами решения задач МФ и ГУ.

На этапе исследований с аппроксимацией уравнений (1), (2) методами конечноразностной аппроксимации (КРА) предполагалось решение задачи разбить на два этапа. На первом - подготовка, расчет ГУ и ввод их в СП, на втором - расчет системы линейных алгебраических уравнений. При этом каждому этапу ведения работы соответствовала своя техническая реализация.

Выполнение работ по расчету и представлению ГУ с единых позиций в рамках единого вычислительного потока [1] обеспечит дополнительные преимущества СП, ориентированному на решение задач МФ.

Задача Дирихле [3] (первая краевая задача) связана с представлением значений некоторой функции (или постоянной величины), которыми описываются всевозможные физические процессы. Для описания этих процессов используют аналитические уравнения. Применение ЭВМ дало мощный толчок в цифровом описании, расчете и трансформации их. Но ЭВМ, последовательно обрабатывающие число за числом, плохо вписываются в идеологию параллельных структур.

Применение ЭВМ в качестве HOST-машин для расчета и ввода ГУ в спецпроцессор задерживает и удорожает весь процесс решения задачи. Ранее в случаях, когда при моделировании или решении не требовалась высокая скорость обработки информации и отсутствовали ограничения на время решения задачи, их применение можно было считать оправданным. Однако переход пользователей к решению трехмерных задач, а машин - к работе в режиме реального времени создает определенные трудности при организации вычислительного процесса.

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

Смешанная (третья краевая) задача так же, как и сопряженная задача (ГУ IV-го рода), усложняет задачи 1 и 2-го рода, создает трудности при расчете, преобразовании и вводе ГУ в спецпроцессор параллельного типа, увеличивает общее время решения. Возникает потребность в согласовании исходных данных, представленных числом из HOST -машины, со сплайнами, поскольку применение МКЭ предполагает использование полиномов в качестве значений узловых функций, а их использование в СП предполагает обработку в виде сложных структур данных.

Поэтому процедура расчета, преобразования исходных данных в СП параллельного типа, получение ГУ в виде ССД - важная проблема. Ее рассмотрению посвящена данная работа.

4. Описание полиномов с помощью регулярного матричного представления (РМП)

Предлагается базу данных для функций, которые могут быть использованы в СП в качестве ГУ, формировать и задавать в ПЭ в виде полиномов (степенных рядов, многочленов и т.д.). Они хорошо известны и описаны в [4, 5]. С помощью полиномов Чебышева, ряда

Тейлора можно аппроксимировать тригонометрические, логарифмические, показательные, гиперболические и пр. функции. А поскольку они могут быть записаны в виде полиномов, то и в СП их можно вводить в виде полиномов с использованием регулярного матричного представления [1].

Полиномы вида

/ 1(X) = а 0 Х° + а1 Х + а2 Х 2+---, І2(Х)= Ь0 Х° + Ь1 X ' + Ь2 Х 2+ ••••,

(5)

записанные в виде циркулянта

І1 (X )= [І і(Х)] *

а.

а

а

а

а1...

а

X0

Xі, І 2 ( X )=[ І 2^)]* X2

Ь0 Ь1 Ь2 •••

Ь0 Ь1..

*0...

X0

Xі

X2

(6)

произведение которых у

у=І (X)" І l(X) * І2 (^)=Уo+УlX,+У2 X ’-+•••=

= а0 Ь0 +(а1 Ь0 + а0 Ь1) X '+(а2 Ь0 + а1 *1 + а0 Ь2) X 2 +•

(7)

представлены в матричной форме в виде РМП:

У=[І l(X)]* [ І 2(£)]*

а0 а1

а

а0...

Ь0 Ь1 Ь2...

Ь0 Ь1..

Ь0...

X0

Xі"

X2

(а0Ь0) (а1 Ь0 + а0Ь1) (а0Ь1 + а1Ь1 + а1Ь0)...

(а0Ь0)

(а1Ь0 + а0 Ь1). (а0Ь0)...

X0

Xі

X2

(8)

Замечание. Поскольку запись в виде циркулянта избыточна и при умножении матриц Ат*т и Вт*т необходимо выполнить т 3 операций умножения и иметь т2 скалярных умножителей, в [1] предложена иная запись циркулянтов и, соответственно, иные правила их умножения. С этой целью полином (6), который в (8) записан в виде матрицы, является множимым (и записывается слева), заменяем на матрицу (6) /1 (X), транспонированную к ней

[/1X :

а

0

а1

а

а

(9)

а полином І2 (X) - множитель, определенный первой строкой матрицы [ І2 (X) ] в (6), записывается в виде вектора столбца [ І-^^) ] (запись справа):

[ І2* (X)]*

Т огда произведение функций І1(X) и І2 ( X )

у=І1^) * Л^)* /Т (X) * І2* (X)

(10)

(11)

а0 а1

Ь

0

Ь

Ь

2

а0 Л" (а0Ь0) >0 "

/ * (X)* а1 а0 * ь = (аА + а0Л1) = >

а2 а1 а0 ъ2 (а2Л0 + а1Ъ + а0Л2)1 _ >2 _

=/* (X),

(12)

[/1*(X)]Т =[«0 «1 а2...]. (13)

В итоге имеем полную адекватность с результатом, получаемым при выполнении операции умножения верхних треугольных матриц с РМП. Введение замены / (X) на

/ т (X) и /*(X) переводит процедуру умножения двух матриц размерности т к процедуре умножения матрицы [ / (X)]т на вектор [ /* (X)], сокращая в т раз число СУ, количество операций умножения и память для хранения соответствующих операндов.

Описание рядов с помощью РМП. Примеры

Таким образом, переходя к аппроксимирующим функциям, описанным, например, в [4], для синуса первых пяти членов разложение в ряд Тейлора имеет вид

Ф ( х )= біп р х /2=р х /2—1/3 !(р х /2)3+1/5!(р х /2)5-1/7!(р х /2)7+1/9!(р х /2)9= (14)

13 5 7 9

= а 0 + а 1 х - а 3 х + а 5 х - а 7 х + а 9 х =

= а п х 0+1,5707963 х 1-0,64596410 х 3+0,079692626 х

(14а)

-0,0046817541 х 7+0,00016044118 х 9 *

(14б)

а 0 а^ 0 а3 0 5а 0 а7 0 а9

а о а 1 0 а3 0 а5 0 а7 0

а0 а0 1 а3 0 а5 0 а7

а0 а 1

а0 _

(14в)

Запись функции синуса с помощью разложения чебышевским рядом [4] до 9-го порядка приведения может выглядеть так:

т

Р п (х) = Ъ ак х Л ,

к=0

біп(р х/2)=1,5707963 х '-0,64596336 х 3+0,079688475 х 5-0,0046722203 х 7 +0,00015081716 х 9 =

х

х

= а 0 х 0 + а1 х1 -а 3 х 3 +а 5 х 5-а 7 х 7 + а 9 х 9 *.

0 а1 0 - а3 0 а5 0 - а7 0 а9

0 а1 0 - а3 0 а5 0 - а7 0

0 а1 0 - а3 0 а5 0 -а

(15)

(16) (16а)

а,

хк-1

х

5

0

х

х

х

8

х

9

х

х

7

0

0

=[Ф][Р(£)]г. (1бб)

Когда ГУ представлены с помощью рациональных функций [5] вида

т п

/(х)= Ё акх'/(Ё Ьрхр)’ (17)

к=0 р=0

в ЭВМ с последовательным вычислением функций типа

/п ( х )= Р0( х УО.( х ) (17а)

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

В записи с РМП /п (х) п (вида (17а)), когда приходится делить на полином О(х), следует переходить к записи умножения на обратную матрицу:

/п (х)= Р 0(х)* О "Чх), (18)

вычислив для этого обратную матрицу О"'(х). Такая процедура сохраняет технологию ввода и обработки значений ГУ по структуре с остальными данными и с заданием их в параллельный спецпроцессор в виде сложных структур данных.

5. Рекомендации по вводу ГУ для задач МФ

Если остановиться на предложенной выше структуре ГУ (не описывая весь набор функций, представленных в [4, 5]), предложим несколько рекомендаций, которые могут быть полезны при таком способе задания ГУ.

1. Из уравнений (14, 16) видно, что при задании функций с точностью до 9 степени при неизвестном, на практике мы часто имеем половину коэффициентов, значения которых отличны от нуля. Тем не менее, матрицы для их представления в виде РМП должны учитывать нулевые коэффициенты и иметь порядок 9, то есть представление избыточно. А полиномы, которые используют для аппроксимации узловых функций с помощью МКЭ, имеют второй или третий порядок. Поэтому необходимо контролировать потребность в задании высокой точности полинома на границе области; использовать такие описания с помощью многочленов (разложения функции в ряд), у которых коэффициенты отличны от нуля, а размерность матрицы циркулянта ниже и согласована с размерностью функций, аппроксимирующих узловые функции по всей области решения.

Кроме этого, с целью согласования размерности матриц циркулянта на границе области с размерностью СУ (с разрядной сеткой его) можно использовать свойства матриц:

[а А ]=а[ А ], [ А + В ]=[ А ]+[В ], [ А + В ][С ]=[ А С ]+[В С ]. (19)

2. Коэффициенты при неизвестных следует заранее рассчитывать и нормировать. Это согласуется с тем, что числа, представленные в ПЭ в форме фиксированной запятой (ФЗ), должны быть нормированы. Надо помнить, что числа матрицы СЛАУ, которую следует решить, нормированы, матрица - положительно определена, на главной диагонали стоят коэффициенты с наибольшим значением и работает принцип максимума (для эллиптических и параболических уравнений).

3. Произведение матриц размерности тх и т2 имеет размерность т = т1 + т2+1.

4. Полиномы, записанные в матричном виде с помощью РМП (особенно, если матрицы от них подлежат обращению) с коэффициентами на главной диагонали, равными 1 (ьо = 1), имеют преимущество. Оно определено тем, что при обращении матрицы по схеме Гаусса используется только обратный ход (прямой ход для получения верхней диагональ-

ной матрицы отсутствует). Поэтому не требуется операция деления (взятия обратной величины от коэффициента Ь0, стоящего на главной диагонали). Всё сводится к операциям умножения и сложения, на выполнение которых и ориентирован скалярный умножитель ПЭ.

6. Применение для решения СЛАУ по СПИ и МП

Часто для решения (1) используют схему простой итерации (СПИ) с конечно-разностной аппроксимацией уравнений [2, 3, 6]:

(5+1) (5) (0) (0)

У = 5 У +@, 5 = Е-т А, (р = т /, ^ =0. 1. ... w, У = У (х 0), х . е ^, (21)

где т - итерационный параметр.

К (21) сводится большинство сеточных аналогов задач МФ с оператором 5, который представляют в виде

(о (.)

5у (х )= 2 В(х х )у(х ^ (22)

к к

хке Ш (хг), где В ( х 1. х^ ) обычно удовлетворяют условия

0< В ( хг . хкг))= В (х . ) < 1. В ( х 0)<1. В ( хг . х кг)) ^ 0. х г , хг 0 е ™ .

Для аппаратной реализации имеется [6] метод приращений (МП), который отрабатывает оператор (21, 22). Метод итерационный, близок к схеме простой итерации (СПИ) по Гауссу. Он имеет большую скорость сходимости и выше по точности решения, но требует некоторой специфики в структуре узлового процессорного элемента. При его реализации в каждом ПЭ выделяют старший значащий разряд узловой функции на данной итерации. После этого в каждом ПЭ складывают с учетом веса выделенного разряда значения именно тех коэффициентов, в которых он выделен. Практически это сложение адекватно операции умножения на полноразрядное число, равное значению узловой функции (в которой выделен старший значащий разряд). Таким образом, ускоряется вычислительный процесс, ибо взамен умножения полноразрядных значений функций на соответствующий им коэффициент они одноразово складываются с учетом веса разряда, выделенного в текущей итерации. А число разрядов может быть от 128 и более.

7. Обсуждение

Решение краевой задачи начинается с ввода значений граничных условий в контурные узлы. В ходе решения задачи (с помощью ЭВМ или СП) в процессе вычислений сеточных функций в узлах области сначала «отрабатывается» низкочастотная часть численных значений функций, задаваемых с помощью ГУ, поскольку им соответствуют старшие значащие разряды граничных условий. Далее, по мере «проникновения» численных значений ГУ внутрь области решения, в узловых ПЭ появляется высокочастотная составляющая ГУ. Поэтому в течение первых итераций, практически до получения «грубого» решения, во всех узлах, лежащих в области вблизи контура, часто многие ПЭ «перемалывают» вхолостую старшие разряды, несмотря на то, что их значение уже давно установилось. Разница только в младших значащих разрядах узловых функций, окончательное значение в которых получают в конце решения.

Представление ГУ в виде степенного ряда (полинома) или в записи РМП, где численные значения а0 х0, а1 х1, а2 х2 . начинаются с главной части функции и далее идут

по убывающей величине, хорошо вписывается в алгоритм решения задач МФ с помощью МП.

Полиномиальное представление функции обладает определенным свойством. В этом представлении каждая следующая компонента численного значения ai при xг сплайна, вводимая в узловой процессор, дискретна и по абсолютной величине меньше предыдущего значения аг-1 при хг-1. Это свойство обеспечивает задание ГУ в виде РМП и для параллельных структур, ориентированных на решение задач МФ методом приращений. Объективно оно хорошо сочетается с МП, в котором обработка последующего значения функции в узле сетки начинается со значения ее на данной итерации или, как минимум, вдвое меньше.

Данное свойство хорошо сочетается с вычислением численных значений сеточных функций, размерность которых превышает размерность регистровых структур в ПЭ и СУ в их составе. Кроме того, оно является дополнительным преимуществом при обработке «безразмерных» значений сеточных функций на структурах конечной разрядности.

8. Выводы

Предложенный и описанный способ задания граничных условий хорошо вписывается в структуру параллельного спецпроцессора.

Он органически вписывается в единый технологический поток решения задач математической физики.

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Ледянкин Ю.Я. Единый технологический поток в организации вычислений - способ повышения производительности параллельных структур на процессорных элементах транспьютерного типа / Ледянкин Ю.Я. - Киев, 1989, 1989. -20 с. - (Препринт / АН УССР, Ин-т кибернетики им. В. М. Глушкова; 89-57).

2. Ледянкин Ю.Я. Вопросы создания спецпроцессоров, проблемно-ориентированных на решение задач математической физики. - Киев: ИК АН УССР, 1980. - 54 с. - (Препринт / АН УССР, Ин-т кибернетики; 80-15).

3. Малиновский Б.Н Анализ и синтез цифрового сеточного микропроцессора для решения уравнений теории поля / Малиновский Б.Н., Ледянкин Ю.Я. и др. - Киев: ИК АН УССР, 1978. -55 с. -(Препринт; 78-75).

4. Численные методы и программирование на ФОРТРАНе / пер. с англ. - М.:Мир, 1977. - 584 с.

5. Благовещенский Ю.В. Вычисление элементарных функций на ЭВМ / Ю.В. Благовещенский, Г.С. Теслер. - К.: Техника, 1977. - 208 с.

6. Ледянкин Ю.Я. Метод приращений для решения уравнений Ау = / / Ю.Я. Ледянкин // ЖВМ и МФ, 1979. - Т. 19, № 4. - С. 1043 - 1047.

Стаття надійшла до редакції 05.10.2011