раз, а разлагает его на одномерные преобразования, перемежающиеся перестановкой данных, – это позволяет использовать наличие нулей во входных данных для прямого БПФ и избыточность выходных данных для обратного, чтобы оптимизировать вычисления [35]. Во-вторых, наличие четырёх различных итерационных методов позволяет выбрать наиболее подходящий для конкретной задачи, поскольку, как известно из литературы (см. параграф 2.1.4.1) и будет показано в подразделе 2.4.3, нельзя однозначно выбрать лучший итерационный метод в МДД – в зависимости от конкретной задачи любой метод может оказаться лучше других. В-третьих, динамическое выделение памяти и оптимальная структура данных приводят к тому, что все вычисления кроме БПФ выполняются только для реальных (непустых) диполей, а не для всей прямоугольной решётки – это также уменьшает требуемую ADDA память. В-четвёртых, симметрия матрицы взаимодействия позволяет использовать меньше памяти для хранения её преобразования Фурье. И последнее, все переменные с плавающей точкой представлены в ADDA в двойной точности – это ускоряет сходимость в тех случаях, когда ограничением является машинная точность. Более того основные операции с числами двойной точности выполняются быстрее чем с числами одиночной точности на современных процессорах. Хотя это ускорение возникает за счёт увеличения используемой памяти, она тем не менее ниже чем для аналогичных программ (см. раздел 2.5).

Дополнительная информация об ADDA представлена в обширной инструкции пользователя, входящей в распространяемый комплект [181]. Текущая версия ADDA 0.77 (на момент написания – лето 2007 г.) была выпущена в июне 2007 г.

2.4.3. Численное моделирование

2.4.3.1. Параметры моделирования

Для тестовых задач использовалась ADDA 0.75, скомпилированная с помощью компилятора Intel C 9.0 с максимальными возможными оптимизациями (настройки по умолчанию в файле Makefile). Все вычисления в данной диссертации были проведены на Голландском вычислительном кластере LISA,* используя 32 узла (на каждом два процессора Intel Xeon 3.4 ГГц и 4 ГБ памяти), если не указано обратное. Мы использовали ДСР, как наиболее часто используемое выражение для поляризуемости, и три разных итерационных метода: КМН, Би-СГ и Би-СГСтаб. Для всех трёх использовался критерий остановки по умолчанию – εit = 10−5.

* http://www.sara.nl/userinfo/lisa/description/

102

Показатель преломления m

2.0 2 ГБ

70 ГБ

1.8

1.6	εit (10−5,10−3)

1.4

εit =10−5

1.2

1.0

0	20	40	60	80	100	120	140	160
			Размерный параметр x

Рис. 16. Текущие возможности ADDA для шаров с различными x и m. Заштрихованная область соответствует полной сходимости, а перекрёстно заштрихованная – неполной. Пунктирные линии отображают два уровня постоянного объёма памяти в соответствии с эмпирическим правилом Дрейна (описание в тексте).

В качестве тестовых задач использовались шары с x от 20 до 160 и m от 1.05 до 2, мы ограничились вещественными m. Текущие возможности ADDA показаны в виде области на плоскости (x,m) на рис. 16: заштрихованная область соответствует полной сходимости, а перекрёстно заштрихованная – случаям, когда ADDA не смогла уменьшить норму невязки до требуемого уровня, а только до εit (10−5,10−3). Хотя неполная сходимость возможно незначительно влияет на итоговую точность характеристик рассеяния, мы не рассматриваем в дальнейшем эти результаты, так как требуется отдельное исследование, чтобы численно описать эффекты неполной сходимости (см. подраздел 2.4.4). Для полностью сошедшихся результатов вклад ненулевой невязки в суммарные ошибки характеристик рассеяния пренебрежимо мал (данные не приведены).

Полный набор пар значений (x,m), для которых ADDA сошлась, приведён в таблице 5 – она также отображает количество диполей на длину волны в частице

(λ/md). Мы старались выбрать его равным 10 в соответствии с эмпирическим правилом Дрейна, но реально оно немного отличалось, так как мы варьировали размер решётки для оптимизации параллельной эффективности ADDA.* Полное количество диполей в прямоугольной решётке, приведённое в таблице 5, изменялось от 2.6×105 до 1.3×108, его можно оценить как (3.18xm)3. Этому количеству пропорциональны и требуемая

* Наибольшая параллельная эффективность достигается, когда размер решётки делится на количество процессоров. Однако ADDA может работать с любым размером решётки.

103

Таблица 5. Параметры численного моделирования.

	m	x	λ/md	N a	Итерационный метод	Niter
	1.05	20	9.6	2.6×105	Би-СГСтаб	6
		30	9.6	8.8×105	Би-СГСтаб	7
		40	9.6	2.1×106	Би-СГСтаб	9
		60	9.6	7.1×106	Би-СГСтаб	14
		80	9.6	1.7×107	Би-СГСтаб	20
		100	9.6	3.3×107	Би-СГСтаб	27
		130	10.3	9.0×107	Би-СГСтаб	40
		160	9.6	1.3×108	Би-СГСтаб	65
1.2		20	10.5	5.1×105	КМН	86
		30	11.2	2.1×106	КМН	223
		40	10.5	4.1×106	КМН	598
		60	9.8	1.1×107	КМН	2120
		80	10.5	3.3×107	Би-СГСтаб	21748
		100	10.1	5.7×107	Би-СГСтаб	6169
		130	10.3	1.3×108	Би-СГСтаб	29200
1.4		20	10.8	8.8×105	КМН	1344
		30	10.8	3.0×106	КМН	16930
		40	10.8	7.1×106	КМН	8164
		60	9.6	1.7×107	Би-СГ	127588
1.6		20	11.0	1.4×106	КМН	8496
		30	10.5	4.1×106	Би-СГ	69748
1.8		20	11.2	2.1×106	КМН	28171
		30	10.2	5.5×106	Би-СГ	118383
2		20	10.1	2.1×106	КМН	58546

a Это количество диполей во всей прямоугольной решётке, которое в основном определяет время одной итерации. Для шаров количество диполей, реально занятых рассеивателем, примерно в два раза меньше.

память, и время одной итерации. Две пунктирные линии на рис. 16 отображают требуемую память для разных x и m, они соответствуют типичному объёму памяти современного ПК (2 ГБ) и максимальной памяти, использованной в наших вычислениях, (70 ГБ) соответственно.

Для каждого шара вычислялись: эффективность экстинкции, параметр асимметрии, и матрица Мюллера в одной плоскости рассеяния, которая является плоскостью симметрии для кубической дискретизации шара. Точные результаты для этих же шаров были получены с помощью теории Ми [14]. Мы использовали сферическую симметрию задачи, чтобы получить все результаты из МДД моделирования одной поляризации падающей волны, поэтому время вычисления примерно в два раза меньше чем для несимметричных частиц с такими же x и m. Мы использовали корректировку объёма, чтобы гарантировать равенство объёмов шара и его дипольного представления [47], однако для очень больших шаров эта поправка чрезвычайно мала.

2.4.3.2. Результаты

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

104

невязки	100
	10-1
норма	10-2
Относительная	10-4
	10-3
	10-5	0	5000	10000	15000	20000	25000
				Количество итераций

Рис. 17. Сходимость метода КМН для шаров с x = 20 и m = 1.8 – показана невязка как функция количества итераций. Система состоит из 3×106 линейных уравнений.

	106						1 неделя
	105						1 сутки
t, с	104
вычислений							1 час
	103							m =
								1.05
Время	102					1 мин
								1.2
								1.4
	101							1.6
								1.8
	100							2
		40	60	80	100	120	140	160
	20
			Размерный параметр x

Рис. 18. Суммарное (астрономическое) время вычислений (на 64 процессорах) в логарифмическом масштабе для шаров с разными x и m. Горизонтальные пунктирные линии, приведённые для удобства, соответствуют минуте, часу, суткам и неделе.

конкретный пример сходимости итерационного метода – применение КМН к системе из 3×106 линейных уравнений для шара с x = 20 и m = 1.8. Суммарное время вычислений t для всех частиц приведено на рис. 18, а относительные ошибки для эффективности экстинкции Qext и параметра асимметрии <cosθ > – на рис. 19. Максимальные и среднеквадратичные (СК) относительные ошибки S11 по всему диапазону углов рассеяния показаны на рис. 20. Ошибки других нетривиальных элементов матрицы Мюллера ведут себя схожим образом (данные не приведены).

Результаты МДД для S11(θ ) для шара с x = 160 и m = 1.05 сравниваются с теорией Ми на рис. 21, врезка показывает в увеличении область углов вблизи направления

105