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

2.5.2.4. ADDA

Программа ADDA подробно описана в подразделе 2.4.2. Важным отличием от других программ является использование метода интегрирования Ромберга для усреднения по ориентации, что имеет свои преимущества и недостатки (см. подраздел 2.4.2). В этом разделе использовалась версия 0.7a, итерационный метод КМН(КС), пакет FFTW 3 для вычисления БПФ и ДСР – для поляризуемости. Усреднение по ориентации проводилось по общей схеме, не используя симметрию рассеивателей, что, например, для осесимметричных рассеивателей могло бы сильно сократить время вычислений.

2.5.3. Сравнение программ

Мы сравниваем четыре реализации МДД друг с другом и с точными решениями, полученными на основе теории Ми или вычислении Т-матрицы, для пяти различных форм частиц и двух показателей преломления. Для четырёх форм из пяти проводилось усреднение по ориентации. Мы приводим ошибки вычисленных характеристик рассеяния, а также время вычислений и требуемую память.

Существует несколько факторов, которые нужно учитывать при создании хорошей программы на основе МДД. Они включают теоретические и практические аспекты, такие как:

1)Интерфейс пользователя программы (задача программирования).

2)Выбор выражения для поляризуемости диполей (физическая задача).

3)Выбор итерационного метода (задача прикладной математике).

4)Реализация итерационного метода (задача программирования).

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

2.5.3.1. Формы объектов и параметры

Использовались следующие пять форм частиц: шар, сфероид, цилиндр и кластеры из 4 и 50 одинаковых шаров – все они имеют одинаковый объёмный размерный параметр x = 5.1. Рассмотрены два показателя преломления: msi = 1.6 + 0.001i (силикат)

и mic = 1.313 (лёд). Размер диполя соответствовал kd = 0.3, но немного варьировался, чтобы объём дипольных дискретизаций всех форм был одинаков. При таком размере диполя ожидается хорошая точность (см. параграф 2.1.3.2).

116

Рис. 23. Формы частиц, используемых в сравнении. Верхний ряд слева: шар и сфероид, нижний ряд слева: цилиндр и кластеры из 4 и 50 шаров. Отношение диаметра к высоте для сфероида и цилиндра равно 0.5. Рисунки построены по сглаженной дипольной дискретизации частиц.

При указанных размере диполя и размерном параметре, каждая частица требует для своего описания примерно 21 000 непустых диполей. Требуемая компьютерная память для МДД определяется в основном размером прямоугольной решётки описанной вокруг частицы: для сферы это кубическая решётка с ребром в 35 диполей, а

для самой пористой формы, кластера из 50 шаров, требуется решётка 57×59×62 диполей. Дискретизации всех пяти форм показаны на рис. 23.

Шар обладает сферической симметрией, поэтому для него усреднение МДД результатов по ориентации кажется бессмысленным, однако его кубическая дискретизация не полностью сферически симметрична, и считается, что усреднение по ориентации улучшает точность, особенно для рассеяния в заднюю полусферу.* Тем не менее шар моделировался в одной ориентации. Для частиц других форм результаты усреднялись по ориентации и сравнивались с результатами аналитического усреднения ориентации по Т-матрице. Мы пробовали различное количество ориентаций для усреднения, но приводятся только результаты полученные с использованием 1024 ориентаций (при этом, по техническим причинам DDSCAT использовала 1089, а ADDA

– 1152 ориентаций). Критерий сходимости итерационного метода для всех программ составлял εit = 10−5.

Разные программы по разному усредняют по ориентации, перебирая три угла Эйлера: α, β и γ систематически или случайно. При этом поворот частицы по α можно заменить на поворот плоскости рассеяния, что намного быстрее (см. подраздел 2.4.2).

* А. Пентила (A. Penttila), частное сообщение (2007).

117

SIRRI абсолютно случайно перебирает все три угла. ZDD поступает так же, но для каждой ориентации частицы вычисляет рассеяние в четырёх плоскостях – таким образом 1024 суммарных ориентаций достигается за счёт 256 поворотов частиц. DDSCAT и ADDA систематически перебирают углы Эйлера. Для DDSCAT

использовалось 11 плоскостей рассеяния и 9 (β ) на 11 (γ ) ориентаций частицы (всего 1089). Для ADDA выбор количества значений углов ограничен, поскольку они должны быть степенями двойки (плюс один в некоторых случаях) из-за метода интегрирования Ромберга. Поэтому мы использовали 16 плоскостей рассеяния и 9 (β ) на 8 (γ ) ориентаций частицы (всего 1152).

2.5.3.2. Точные методы

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

1)Программа Ломпадо (Lompado)* на основе теории Ми для шаров (написана в

Mathematica).

2)Программа МРГУ Мищенко и Тревиса (Travis) для цилиндра и сфероида [186] (написана на Фортране 77).

3)Программа ОММЧ Маковски и Мищенко для кластеров шаров [187] (написана на Фортране 77).

2.5.3.3. Точность

Результаты программ МДД сравнивались с точными методами. Все программы выдают полную матрицу Мюллера для заданного набора углов рассеяния θ (мы использовали диапазон от 0° до 180° с шагом в один градус), однако мы приводим лишь интенсивность рассеяния S11 и степень линейной поляризации P = −S21/S11. Поскольку S11 имеет большой динамический диапазон, обычно её отображают в логарифмическом масштабе. Поэтому мы приводим относительные ошибки S11 (со знаком), в то время как для P более подходящими являются абсолютные отклонения. Эти ошибки вместе с точными значениями приведены на рис. 24–28 для всех форм и обоих показателей преломления.

* http://diogenes.iwt.uni-bremen.de/vt/laser/codes/Mie%20Code%20-%20Lompado-new.zip

118

	103								100
			(а)			msi			75	(б)
	102					mic		P, %	50
									25
)								Поляризация
(θ
11									0
S
	101								-25
									-50
	100								-75
,%	15			SIRRI		DDSCAT				(г) msi
11				ZDD		ADDA
S	10								5
Относительная ошибка
	5
	0								0
	-5	(в) msi
	-10								-5
	-15							Ошибка%P,
, %	15	(д) mic								(е) mic
11
S	10								5
Относительная ошибка
	5
	0								0
	-5
	-10								-5
	-15							Ошибка%P,
	-200	30	60	90	120	150	180		-100	30	60	90	120	150	180
			Угол рассеяния θ, градусы								Угол рассеяния θ, градусы
Рис. 24. Сравнение результатов для шара. Точные значения S11 (в логарифмическом масштабе)
и P показаны в частях (а) и (б). Ошибки интенсивности и поляризации МДД программ
приведены в (в) и (г) для msi и в (д) и (е) для mic соответственно. Результаты DDSCAT, ZDD и
ADDA практически совпадают для частиц с фиксированной ориентацией. Заметим, что
масштабы обеих осей одинаковы на рис. 24–28.

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

В целом видно, что типичные относительные ошибки S11 в МДД составляют от 2%

до 6% для типичных показателей преломления силикатов и льда и размеров x ≈ 5, а для P типичные ошибки – от 1% до 3%. Обе ошибки в среднем больше для большего показателя преломления. Для частиц в фиксированной ориентации DDSCAT, ZDD и

ADDA приводят к практически одинаковым результатам, что и ожидается, так как они

119

	103								100
			(а)			msi			75	(б)
	102					mic		P, %	50
									25
)								Поляризация
(θ
11									0
S
	101								-25
									-50
	100								-75
, %	15	(в) m								(г) m
11		si								si
S	10								5
Относительная ошибка
	5
	0								0
	-5
	-10			SIRRI	DDSCAT				-5
	-15			ZDD	ADDA			Ошибка%P,
, %	15 (д) mic									(е) mic
11
S	10								5
Относительная ошибка
	5
	0								0
	-5
	-10								-5
	-15							Ошибка%P,
	-200	30	60	90	120	150	180		-100	30	60	90	120	150	180
			Угол рассеяния θ, градусы								Угол рассеяния θ, градусы
Рис. 25. То же, что и рис. 24, но для сфероида в случайной ориентации.

основаны на одинаковых физических принципах. Их точность немного лучше чем SIRRI, что согласуется с результатами других исследователей, которые сравнивали ДСР и ЛАК (см. параграф 2.1.3.2).

Точность при усреднении по ориентации разная для всех четырёх программ, что в основном связано с различными методами усреднения. DDSCAT и ZDD в целом сравнимы по точности, а SIRRI систематически немного менее точна, возможно из-за использования поляризуемости ЛАК вместо ДСР. Точность ADDA сравнима или немного хуже чем ZDD и DDSCAT для кластеров из 4 и 50 шаров, но намного хуже для сфероида и цилиндра. Возможное количество значений углов Эйлера для усреднения не может принимать любые значения для ADDA, поэтому приходится выбирать компромисс между скоростью (больше плоскостей рассеяния) и точностью (больше значений углов Эйлера). Наш выбор 16 плоскостей и 9×8 углов для ADDA приводит к

120

	103								100
			(а)			msi			75	(б)
	102					mic		P, %	50
									25
S								Поляризация
)
(θ
11									0
	101								-25
									-50
	100								-75
, %	15	(в) msi								(г) msi
11	10								5
S
Относительнаяошибка	-15			ZDD		ADDA		ОшибкаP, %
	5
	0								0
	-5
	-10			SIRRI		DDSCAT			-5
, %	15	(д) m								(е) m
11		ic								ic
S	10								5
ошибкаОтносительная
	5							,PОшибка%
	0								0
	-5
	-10								-5
	-15
	-200	30	60	90	120	150	180		-100	30	60	90	120	150	180
			Угол рассеяния θ, градусы								Угол рассеяния θ, градусы

Рис. 26. То же, что и рис. 24, но для цилиндра в случайной ориентации.

хорошей скорости, быстрее других программ от 2.3 до 16 раз, но немного худшей точности для несимметричных частиц. Удвоение числа реальных ориентаций в ADDA значительно улучшает точность и при этом она остаётся самой быстрой программой (данные не приведены).

В дальнейшем мы объясним особенно плохую точность ADDA для осесимметричных рассеивателей. Для таких форм дискретизация симметрична относительно поворота на 90° вокруг оси, а диапазон γ от 0° to 90° симметричен относительно середины (угла 45°). Поэтому наиболее эффективным способом усреднения является равномерное расположения значений γ в интервале от 0° до 45°,

используя при этом меньше значений γ и больше – β. Случайный выбор ориентации удовлетворяет обоим этим условиям. DDSCAT использовала 11 значений γ равномерно расположенных каждые 32.7°, которые при отображении в интервал от 0° до 45° 121

	103								100
			(а)			msi			75	(б)
	102					mic		P, %	50
									25
)								Поляризация
(θ
11									0
S
	101								-25
									-50
	100								-75
,%	15	(в) msi		SIRRI		DDSCAT				(г) msi
11				ZDD		ADDA
S	10								5
Относительная ошибка
	5
	0								0
	-5
	-10								-5
	-15							Ошибка%P,
, %	15	(д) mic								(е) mic
11
S	10								5
Относительная ошибка
	5
	0								0
	-5
	-10								-5
	-15							Ошибка%P,
	-200	30	60	90	120	150	180		-100	30	60	90	120	150	180
			Угол рассеяния θ, градусы								Угол рассеяния θ, градусы
Рис. 27. То же, что и рис. 24, но для кластера из 4 шаров в случайной ориентации.

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

была совершенно другая: использовалось 8 значений γ, расположенных через 45°, т.е.

они соответствовали только двум разным значениям (0° и 45°) в уменьшенном интервале. Это значит, что для осесимметричных частиц ADDA повторяла вычисления для четырёх эквивалентных наборов ориентаций. Более того метод интегрирования Ромберга, который является методом высокого порядка, приводит к неудовлетворительным результатам при применении к периодической функции на интервале, который содержит несколько её периодов и только несколько точек в каждом периоде. Аналогичные аргументы применимы и к β, так как и сфероид, и цилиндр имеют плоскость симметрии перпендикулярную их оси.

122

	103
				(а)
	102
(θ )
11
S
	101
	100
, %	15	(в)	msi
11	10
S
ошибка	5
Относительная	0
	-15					ZDD
	-5
	-10					SIRRI
, %	15	(д)	mic
11	10
S
ошибка	5
Относительная	0
	-15
	-5
	-10
	-200		30	60	90	120
				Угол рассеяния θ, градусы

	100
msi			75	(б)
mic	%		50
	P,		25
	Поляризация
			0
		-25
		-50
		-75
				(г)	msi
			5
	P, %
	Ошибка		0
DDSCAT			-5
ADDA
				(е)	mic
			5
	P, %
	Ошибка		0
			-5
150	180	-100			30	60	90	120	150	180
						Угол рассеяния θ, градусы

Рис. 28. То же, что и рис. 24, но для кластера из 50 шаров в случайной ориентации.

Важно отметить, что ADDA полностью подходит для эффективного усреднения по ориентации осесимметричных частиц, но для этого пользователь должен вручную изменить интервалы для углов Эйлера во входном файле. В частности, для рассматриваемых рассеивателей более подходящий выбор это 17 значений β от 0° до

90° и 5 значений γ от 0° до 45°. Для ледяного сфероида это приводит к СМОО S11 равным 0.82% (сравните с 4.82% в таблице 6) используя лишь на 20% больше времени. Однако, подобная модификация входных файлов рассматривалась в данной работе как «тонкая настройка» и поэтому не использовалась. Уменьшение интервалов может помочь также и DDSCAT, но в намного меньшей степени (для этих конкретных задач).

При внимательном рассмотрении структуры ошибок на рис. 24–28 видно, что ошибки в основном распределены случайно. В частности, средние ошибки (с учётом знаков) бывают как положительные, так и отрицательные примерно в равных

123