parallel

Здесь – коэффициент, определяющий вид разностной схемы (0 1), если 0, то схема явная, если 1, то схема неявная. Это два крайних случая, при 0 1 получаем схемы смешанного типа. Особого внимания заслуживает схема Кранка–Николсона, для которой 0.5. Эта схема имеет второй порядок аппроксимации.

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

Tn 1

2 Tn 1

Tn 1

i 1, j

i, j

i 1, j

2

hx

Tin1, j 2 Ti,nj Tin1, j

(1 )

Tn 1

Tn

2

hx

i, j

i, j

;

n 1

n 1

n 1

Ti, j 1

2 Ti, j

Ti, j 1

2

hy

Tn

2 Tn

Tn

(1 )

i, j 1

i, j

i, j 1

2

hy

i 1,Nx;

j 1,Ny; n 0,1,2,...

T (x , y

);

Tn 1

T (x

, y

);

j 0,Ny 1; n 0,1,2,...

Tn 1

j

Nx 1

j

0, j

b

0

Nx 1, j

b

Tb (xi , y0 );

Ti,nNy1 1

Tb (xi , yNy 1);

Ti,0n 1

i 0,Nx 1; n 0,1,2,...

T0

100;

i 0,Nx 1;

j 0,Ny 1.

i, j

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

171

Какую разностную схему получим после аппроксимации – явную, неявную или смешанную (явно-неявную), зависит от того, чему равен коэффициент при аппроксимации производных по пространству.

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

9.1Явная схема

Вэтом случае ( 0) для аппроксимации производных по пространству используются значения сеточной функции Ti, j с n -го

временного слоя. Тогда получим:

n

n

n

Tn 1

Tn

Ti 1, j 2 Ti, j

Ti 1, j

i, j

i, j

2

hx

Ti,nj 1 2 Ti,nj Ti,nj 1

;

2

hy

i 1,

Nx;

j 1,Ny;

(9.3)

;

T0,nj1 Tb x0, yj

TNxn 11, j Tb xNx 1, yj ;

j 0,...,Ny 1;

n 0,1,2,...;

n 1

Tb xi , y0

;

n 1

Tb xi , yNy 1 ;

Ti,0

Ti,Ny 1

n 0,1,2,...;

i 0,...,Nx 1;

j 0,...,Ny 1.

T0 100; i 0,...,Nx 1;

i, j

Таким образом, в (9.3) представлена явная формула для вычислений значений сеточной функции на новом временном слое. Перепишем полученную формулу в следующем виде:

Tn 1	(1 ap)Tn	ae Tn	aw Tn	an Tn	as Tn	.	(9.4)
i, j	i, j	i 1, j	i 1, j	i, j 1	i, j 1

Коэффициенты ap, ae, aw, an, as легко определяются из форму-

лы (9.3):

172

	ae		; aw		;	an		;	as		;
		hx2		hx2			hy2			hy2

ap ae aw an as .

Формула (9.4) легко программируется и очень проста для применения:

Do i=1,Nx Do j=1,Ny

Tnew(i,j)=(1+ap)*T(i,j)+ae*T(i+1,j)+aw*T(i-1,j) +an*T(i,j+1)+as*T(i,j-1)

End do End do

Далее обсудим возможности распараллеливания программы для нахождения приближенного решения поставленной задачи. В рассматриваемом примере возможны два различных способа разделе-

ния данных – одномерная или ленточная схема или двухмерное или

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

9.1.1 Одномерная декомпозиция

Для определенности декомпозицию расчетной области проведем по индексу j (рис. 9.1), хотя на практике выбор координаты для

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

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

173

каждого процессорного элемента требуется получить векторы сеточных значений приближенного решения с предшествующей и последующей подобластей декомпозиции. Например, для области R 1 – с областей R 0 и R 2.

R=0

(I,j+1)

(I,j)

R=1

(I-1,j) (I+1,j)

(I,j-1)

j R=2

i

Рис. 9.1 Одномерная декомпозиция расчетной области (открытые кружки представляют собой граничные узлы сетки)

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

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

174

Рис. 9.2 Схема обменов между вычислительными узлами

дают свои верхние граничные строки своим предыдущим соседям и принимают переданные строки отследующих процессорных элементов.

9.1.2 Двумерная декомпозиция

При большом числе используемых процессоров перспективнее становится двумерная декомпозиция. При этом по-прежнему необходимо в адресном пространстве каждого процессорного элемента создавать дополнительные фиктивные ячейки для обеспечения однородности вычислений. Операции обменов схематически изображены на рис. 9.3. На первый взгляд может показаться, что обменов стало больше. Проверим, что происходит на самом деле. Посчитав число одновременно пересылаемых элементов при каждом способе

175

декомпозиции, заметим, что в 1d-декомпозиции количество пересылаемых данных есть величина порядка 2n (n – размерность задачи, т.е. количество узлов сетки по каждому координатному направлению), а в 2d-декомпозиции количество пересылаемых данных есть

величина порядка 4n , где p – число используемых процессоров. p

Оценки приведены для решения двумерной задачи.

Таблица 9.1 Количество пересылаемых данных при различных способах декомпозиции и различном числе процессоров

Число процессоров	2	4	8	16
1d-декомпозиция	2 n	2 n	2 n	2 n

2d-декомпозиция

4 n

2 n

2 n

n

2

2

Рис. 9.3 Схема обменов при двумерной декомпозиции

176

Таким образом, из табл. 9.1 видно, что при p 4 двумерная де-

композиция является более перспективной для распараллеливания. В то же время реализация 2d-декомпозиции является более трудоемким (с точки зрения программирования) процессом по сравнению с реализацией 1d-декомпозиции. Основная возникающая сложность состоит в организации межпроцессорных обменов. Есть различные способы решения этой проблемы. Ниже представлен вариант с использованием декартовой топологии и новых типов для организации обменов, а также вариант с перенумерацией неизвестных значений сеточной функции для сбора решения на одном процессорном элементе. Эти два способа организации обменов являются взаимозаменяемыми.

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

9.2 Неявная схема

В этом случае ( 1) для аппроксимации производных по пространству используются значения сеточной функции Ti, j с n +1-го временного слоя. Тогда получаем:

							n 1	n 1	n 1
Tn 1			Tn			Ti 1, j 2 Ti, j			Ti 1, j
	i, j		i, j					2
							hx
				Ti,nj 11 2 Ti,nj 1 Ti,nj 11
										;
								2
							hy
				i 1,Nx;			j 1,Ny;			(9.5)
		1 Tb x0,yj ;					TNxn 11, j Tb xNx 1, yj ;
T0,nj
n 0,1,2,...; j 0,...,Ny 1;								Tb xi , yNy 1 ;
	n 1		Tb xi , y0 ;				n 1
Ti,0							Ti,Ny 1
n 0,1,2,...; i 0,...,Nx 1;
		100;		i 0,...,Nx 1;				j 0,...,Ny 1.
T0
	i, j

177

Для большей наглядности перепишем полученные формулы в следующем виде:

(1 ap)Tn 1	ae Tn 1	aw Tn 1	an Tn 1	as Tn 1	b	. (9.6)
i, j	i 1, j	i 1, j	i, j 1	i, j 1	i, j

Коэффициенты ap,ae,aw,an,as,bi, j легко определяются из фор-

мулы (9.5):

ae

; aw

; an

; as

;b

Tn

h2

h2

h2

h2

i, j

i, j

x

x

y

y

ap ae aw an as ;

i 1,Nx; j 1,Ny;

i 0;

ap 1;

bi, j Tb (x0 , yj );

ae aw as an 0;

i Nx 1;

ap 1;

bi, j Tb (xNx 1, yj );

ae aw as an 0;

j 0;

ap 1;

bi, j Tb (xi , y0 );

ae aw as an 0;

j Ny 1;

ap 1;

bi, j Tb (xi , yNy 1);

ae aw as an 0.

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

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

178

9.2.1 Метод сопряженных градиентов

Одними из перспективных методов итерационного решения систем вида Ax b являются алгоритмы, построенные на основе выбора итерационных параметров из условия минимизации функционалов, определяющих точность текущих последовательных при-

ближений.
	В дальнейшем будут использованы следующие			обозначения:
rk	b Axk – невязка; k	B 1rk – поправка и	zk	xk x – по-
грешность (ошибка). Здесь		x – точное решение системы Ax b ;
xk	– k -е приближение к точному решению; B		– предобуславли-

вающая матрица.

К открытию метода сопряженных градиентов (CG) независимо пришли М. Хестенес и Э. Штифель. Он является наиболее предпочтительным по быстродействию для симметричных положительно определенных систем. Формулы классического метода сопряженных градиентов имеют следующий вид:

p x;
v Ap;
p b v;
r p;
r	2 2	;

for ( (i i_max).and.( ) )

v Ap;

v, p 2 ; x x p; r r v;

new r2 2 ;

p r new p;

new;

end do.

179

xk ! xk anew a pk

Достаточным условием сходимости метода сопряженных градиентов являются симметричность и положительная определенность матрицы A , при этом если для спектра матрицы выполняется условие 0 m A M , то скорость сходимости можно определить

по формуле

	2	n
xn				x0	;
	1	2n

1 m M

. 1 mM

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

Метод CG имеет то особенно привлекательное свойство, что в его реализации предусмотрено одновременное хранение в памяти лишь четырех векторов. Кроме того, в его внутреннем цикле помимо матрично-векторного произведения вычисляются только два скалярных произведения, три операции типа «saxpy» (сложение вектора, умноженного на число, с другим вектором) и небольшое количество скалярных операций. Таким образом, и необходимая оперативная память, и объем вычислительной работы в методе не очень велики. В CG алгоритме xk 1 вычисляется с использованием рекуррентных соотношений для трех групп векторов. Одновременно в памяти требуется хранить лишь самые последние векторы из каждой группы, которые переписываются поверх ранее рассчитанных значений. Первая группа векторов – это приближенные решения xk . Вторая группа – это невязки rk b Axk . Третья группа – это сопряженные градиенты pk . Эти векторы называют градиентами по следующей причине: шаг метода CG можно рассматривать как выбор числа anew a из условия, чтобы новое приближенное решение

минимизировало норму невязки

1/2

rk A 1 rkT A 1rk . Иными словами, pk используется как направ-

ление градиентного поиска. Они называются сопряженными или, точнее, А-сопряженными, потому что pkT Apj 0 при j k .

Параллельную реализацию метода сопряженных градиентов рассмотрим на примере приближенного решения уравнения теплопро-

180