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

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

69

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

Афанасьев В. П., Будак В. П., Ефременко Д. С., Лубенченко А. В.

Московский энергетический институт (технический университет), 111250 Москва, Россия Поступила в редакцию 17.05.2011, после доработки 20.05.2011.

Представлена чл.-корр. РАЕНЯА. Илюшиным 23.05.2011

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

Ключевые слова: уравнение переноса излучения, плоская геометрия, анизотропное рассеяние

PACS 42.68.Ay____________________________

содержание

1. Введение (69).

2. Формулировка ВУПИ (70).

3. Дискретизация ВУПИ (71).

4. Решение для гладкой части (72).

5. Аналитическое представление анизотропной части (72).

6. Анализ влияния программных и аппаратных способов ускорения (73).

7. Результаты и обсуждение (74).

8. Заключение (76).

Литература (76).

1. ВВЕДЕНИЕ

Решение задач восстановления послойного состава атмосферы на основе данных спутникового зондирования является обратной задачей теории переноса излучения. Как правило, обратная задача несопоставимо сложнее прямой задачи. Из-за плохой обусловленности стандартные методы, основанные на линеаризации модели, могут привести к решениям, лишенным физического смысла [1, 2]. Поэтому часто используются методы подбора, основанные на многократном решении прямой задачи. Кроме того, с развитием экспериментальной техники возрастает массив данных, получаемых со спутников. В силу этого к

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

С другой стороны, возрастают требования к точности вычислений. Так в проекте GOSAT (The Greenhouse gases Observing SATellite) требуется рассчитывать интенсивность рассеянного излучения с погрешностью не более 1% [3]. Для обратных задач важно отсутствие априорных ограничений на параметры среды. При малых оптических толщах хорошо зарекомендовал себя метод TMS [4]. Для больших толщ метод оказывается не столь эффективным. Для аккуратного решения обратных задач требуется точное решение уравнения переноса, которое может оказаться трудоемким. Ускорить алгоритмы точного решения можно за счет программно-аппаратных средств. Поскольку высота атмосферы существенно меньше радиуса Земли, то в настоящее время в задачах дистанционного зондирования принята модель плоской геометрии среды. Соответственно мы сосредоточим наше внимание на методах, алгоритмах решения векторного уравнения переноса излучения (ВУПИ) и особенностях их реализации для плоского слоя мутной среды.

РЕНСИТ | 2011 | ТОМ 3 | НОМЕР 1

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

70 АФАНАСЬЕВ В.П., БУДАК В.П.,

ЕФРЕМЕНКО Д.С., ЛУБЕНЧЕНКО А.В.

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

2. ФОРМУЛИРОВКА ВУПИ

Имеем краевую задачу ВУПИ для слоя мутной среды толщиной т , на который падает плоская волна произвольно поляризованного света в направлении 1° [5]:

д Л

Ц — L(t, i) + L(t, i) = -j R(x) £(U')R(x')L(T,

L(T,i)|

lx=0, (Zi)>0

LoSd - io),

L(t, i)

T=V

(Zj)<0

= 0;

(1)

где L (t,I) - вектор параметров Стокса светового поля в среде на оптической глубине т по направлению I. Введена декартова система координат OXYZ, ось OZ направлена вниз перпендикулярно границе слоя, z - единичный орт вдоль OZ, I = {yj 1 -ц2 cos9, -y/l -ц2 sin9, д}, I° = ц/1 -ц2 °, д(1,z), z). х ДДО - матри-

ца рассеяния, Л — альбедо однократного рассеяния. R (х)- матрица изменения (ротатор) параметров Стокса при повороте плоскости референции. Угол у - двухгранный между плоскостями (z XI) и (IXI) а х’ — между плоскостями (IxI') и (I'хz). Единичные вектора отмечены символом крышечка сверху, столбцы стрелкой вправо, строки — влево, матрицы — двойной стрелкой над символом.

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

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

L (тД) = L а(тД) + L г(тД). (2)

Предполагается, что анизотропную часть можно найти аналитически, что позволит выразить интеграл от нее также аналитически - метод выделения сингулярности [7].

Впервые такое представление предложили Eddington [8] и Milne [9], вычитая из решения прямое нерассеянное излучение. В случае поляризованного излучения выражение для прямого излучения имеет вид

L = L °(w) = eT/ll°R (Ф) L °S(I -

1 " p cos ф0 p sin Фо

q _

L0 = E

где L0 - вектор параметров Стокса падающего излучения, Eq- нормальная к пучку облученность, р, q — степень линейной и циркулярной поляризации пучка, ф0 — азимут поляризации.

Подстановка (2) в (1) изменяет краевую задачу на эквивалентную

+ L, (т,1) =

дт

= 4- ( R (х) х(1, l') R (X') L г (x,l') &' + Д (т, 1);

4п

L (т, 1) 0 0 = 0; L (т, 1)

т=0, ц>0

т=т0, ц<0

= 0,

с функцией источников в правой части

Д(т, 1) = — (| R(x) Х(1, 1') R(x 'Ж (x,l') d\ ' -

-Ц-

4п

д La (Т,1)

дт

- La (Т,1),

(3)

(4)

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

Л(т, 1) = —e-^R(9)Jx(l,l0).

4п (5)

1 НОМЕР | ТОМ 3 | 2011 | РЕНСИТ

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

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

3. ДИСКРЕТИЗАЦИЯ ВУПИ

Kuscer-Ribaric [10] предложили использовать для определения поляризации циркулярный базис

" L+2" Q - iU

о + -1 I -V

CP L-o -2 I + V

_ L-2 _ Q + iU

"0 1 -i 0 " " I"

1 1 0 0 -1 Q

2 1 0 0 1 U

0 1 i 0 V

- T Т

1CS ^SP’

T = T-1

He _ Hs’ (6)

что изменяет вид ротатора и позволяет использовать для представления матрицы рассеяния обобщенные сферические функции (cos 0), для ко-

торых выполняется специальный вариант теоремы сложения

e-"" Pi (l • =

= z (->)’ pm, (' • *)Pi (z • i v*,

q=-k

Аф = ф-ф'. (7)

Однако в этом случае все коэффициенты в ВУПИ становятся комплексными, что затрудняет использование эффективных численных методов решения систем уравнений, поэтому после применения теоремы сложения для обобщенных сферических функций лучше вернуться в стоксово (SP — Stokes presentation) представление. Представим решение задачи (3) в форме [11-13]

M Гф1(тф)Ь'т (т, ц) +

_+Ф2(тФ)ц;(т, ц)

что приведет ВУПИ к виду

IVI

L (т,1) = £ (2 -8„m )

m=0

(8)

ВЛИЯНИЕ АППАРАТНО-ПРОГРАММНЫХ 71 СРЕДСТВ НА СКОРОСТЬ ВЫЧИСЛЕНИЯ

д Lm (т, ц) , L ( ,

Р—-— + Lc (Т Р) =

дт

Л к V

= Т£ (2к + 1)Пm (ц)х, [йm (ц') L (т, ц')ф' +

2 к=0 -1

+А(т, ц), c = 1,2, (9)

где M — число членов ряда Фурье по азимуту, m е 0, M; элементы матрицы рассеяния представляются в виде разложения по обобщенным сферическим функциям

Хк = TSC Хк TCS , Хк = [Xrs ] , [XCP (c0S У)]га =

= Z (2к +1xk Ks(cos YX

к =max(r ,s) (10)

r и s принимают только значения ±0 и ±2, K — количество гармоник в разложении матрицы рассеяния по обобщенным сферическим функциям;

"q: (р) 0 0 0 "

о r: (р) - т; (р) о

о - т; (р) r; (р) о

_ 0 0 о q: (р)_

ф1(ф) = diag (cos ф, cos ф, sin ф, sin ф), ф2(ф) = diag (-sin ф, - sin ф, cos ф, cos ф),

п; (р)=

(11)

Q" (р) =

(k - т)!

К (р),

(k + т)!

R" (р) = 0.5im (?m,2(v) + P,k,.,(v)),

T," (р) = 0.5im (Pm,2 (v) - P,k,-2 (v)).

Для решения полученного уравнения с помощью ДОМ представим интегралы, вхо -дящие в уравнение (9) , с учетом плоской симметрии задачи в виде двойной гауссовой ква -дратуры [14]:

|п m (ц')Е (т, цу

-1

* 0.б£ WjЙ:(ц+)С (т,ц*) +

j-1

+0.5£ Wj П ’ (ц- )Ц (т, ц-),

(12)

где р± = 0.5(^ ± 1), Zj, w — нули и веса гауссовой квадратуры порядка N/2. Отметим, что р± соответствуют нисходящим «+» и нисходящим «-» потокам.

РЕНСИТ | 2011 | ТОМ 3 | НОМЕР 1

72 АФАНАСЬЕВ В.П., БУДАК В.П.,

ЕФРЕМЕНКО Д.С., ЛУБЕНЧЕНКО А.В.

Соответственно для фиксированных m и c (поэтому их опускаем в записи уравнений) имеем систему обыкновенных дифференциальных уравнений d

L± (т) + L± (т) = d т

Л N2 K

= Т X W X (2k + 1)П Г (Е7 )х k

4 j=1 k=0

ГПГ (м+ )L+(t) й V+nГ (E-)L-(т),

+А±(т),

где L± (т) - L” (т, ), А± (т) = А~ (т, ).

(13)

Введем следующие матрицы и векторы

L =

L+ (т)" , А = а+со"

L_(x)_ _ А_ (т)_

"цГ 0 " - Л Щ 0"

М = 0 -ц! , w = — 4 1 0 w,._

X (2k+i)n m (ц+)Xk п m (ц+) х (2k+i)nm (ц+)Xk п ^)

_ k=0 k=0

X (2k+i)n m (e- )x k п m (e+) X (2k+i)n m (e- )x k п m (e- )

k=0 k=0

где

L±(t) =

I(t, e± ), Q(t, e± ), U(t, е± ), V(t, e± ), ...,

1 (T Ел/2 X Q(т, Ел/2 X U(т, Ел/2 X V(т, Ел/2)

= - BL(t) + M-1 A(t),

d т

Ё = M-1(I - AW).

(15)

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Для устранения этого используется масштабное преобразование [15] — умножим уравнение (16) слева на матрицу

S = diag(l, e-r+To), (17)

что приведет (16) к следующему матричному выражению

- S U-1

(18)

"l + (0)" + й U-1 L+СО " J_"

L - (0)_ _L- (0_ _ J+_

где Й = Sefx° = diag(ef-x°, 1),

J т0

A = S fefxU-1M-1 А(т,p0)dx. J 0

J =

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

(19)

" L- (0)" " F_" R _ T_' " L+ (0)"

(14) _L + (to)_ A_ + _ T+ R+_ _L- CO_

и аналогично упорядочены остальные массивы.

Поскольку регулярная часть является гладкой функцией угла, то все матрицы конечны, что позволяет записать систему (13) в матрич-но м виде d L (т)

- T где

F- " J-"

= A

_ F+_ _ J+_

R- f"

= A

_ t+ R+ _

uii

-Г +т0

- еГ-Тои

12

U

21

- U

22

U-1 ^

A =

4. РЕШЕНИЕ ДЛЯ ГЛАДКОЙ ЧАСТИ

Решение полученной системы уравнений имеет анал итический в л д в матричной форме [15]

- U-1L (0) + efx°lJ-1L (т0) =

т0 а

= I" еГт U-1M-1 А(т, |a0)d т,

^ (16) где U - матрица собственных векторов матрицы g, Г = diag(T_, Г+) - диагональная матрица собственных значений упорядоченных по возрастанию так, что Г— = —Г+ .

Проблема решения (16) связана с наличием как отрицательных Г— , так положительных Г+ экспонент в решении, что приводит к быстрому ухудшению обусловленности системной матрицы.

U11 U12

U12

e-r+X0U

е

Г-т0

-]-1

U

22

U

21

Аналогичные выражения для скалярного случая приведены в работе [5].

5. АНАЛИТИЧЕСКОЕ ПРЕДСТАВЛЕНИЕ АНИЗОТРОПНОЙ ЧАСТИ

Аналитическую часть представляли в виде прямой нерассеянной компоненты излучения [8, 9], учитывали одну кратность [16], и две кратности в малоугловом приближении [4]. Логичным шагом будет учесть в L а(т,1) все кратности рассеяния в малоугловом приближении [12, 17, 18]. Такой подход подробно описан в работах [5, 11, 12] на основе малоугловой модификации метода сферических гармоник (МСГ). В этом случае анизотропная часть имеет вид

+

е

U21 U22

1 НОМЕР | ТОМ 3 | 2011 | РЕНСИТ

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

ВЛИЯНИЕ АППАРАТНО-ПРОГРАММНЫХ 73 СРЕДСТВ НА СКОРОСТЬ ВЫЧИСЛЕНИЯ

La (T,l, lo) =

= £ £ (2 - A,o Ж («ф)£ Пm (p) Zk (т)Пm (p) Dc L0,

c=1,2 m=0 к=0

здесь

Zk (x) = exp Г-(1 - лхk) VЖ

(20)

Dj=diag(1 1 0 0), D2=diag(0 0 1 1).

Подставляя последнее выражение в (4) и произведя необходимые преобразования [6, 14], можно п олучить выражение для функции источников

р а m (т, р;)=

= Z (2k +№ - Ро )Пm (Рг) dk Zk (т)П m (Ро )Dc Lо +

k=m

+ A

m

K+1

nm (p, )d кa+i(x)fi m+i(Pc)

-n m+i(p, )d к z к (т)п m (Po)

DL,

0 ’

(21)

где dк = 1 - ЛХk > [A* ]rs = 2 - m2)(k2 - s2)brs.

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

отметим, что все способы выделения анизотропной части не меняют вид матричного уравнения (19). Меняется способ вычисления функции источников. Поскольку (19) относится только к регулярной части, то источники связаны с невязкой анизотропной части. Заметим, что (19) является строгим аналитическим решением краевой задачи для дискретизованного ВУПИ по методу дискретных ординат.

Полученная система (19) является искомой. Из ее решения можно определить распределение яркости отраженного и прошедшего излучений. Перейти от УПИ (1) к системе (19) позволило устранений особенностей из решения (2) и замена интеграла рассеяния (дискретизация) гауссовой квадратурой (12).

При реализации полученного решения (19) алгоритм включает:

• вычисление весов корней и весов квадратурной формулы для дискретизации ВУПИ (12);

• вычисление функции источников (4);

• определение собственных значений и векторов матрицы системы (16);

• перемножение матриц (19).

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

(19) записано для регулярной части, поэтому ее размер определяется не столько оптическими параметрами среды и граничными условиями, сколько тем, насколько точно мы можем выделить анизотропную часть решения. В общем случае M—N—K, но при удачном выделении анизотропной части решения, когда гладкая регулярная часть близка к изотропному распределению по углу, возможно, что M<<N<<K, что существенно повысит эффективность алгоритма. С этой точки зрения очевидно, что наибольшие трудности представляет реализация алгоритма для сильно анизотропного рассеяния.

6. АНАЛИЗ ВЛИЯНИЯ ПРОГРАММНЫХ И АППАРАТНЫХ СПОСОБОВ УСКОРЕНИЯ

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

Все решения этих задач имеют программную реализацию на языках Си и ФОРТРАН [19]. Существует несколько библиотек, в которых собраны требуемые подпрограммы (LAPACK, IMSL, NAG и др.). Для разных платформ существуют свои реализации библиотек. Для платформы Intel существует Math Kernel Library (MKL) [http:// software.intel.com/en-us/intel-mkl], оптимизирован-

ная специально под процессоры Intel. К очевидным достоинствам библиотеки можно отнести тот факт, что она поддерживает многопоточность и автоматически распараллеливает задачу на доступное число ядер, что существенно увеличивает скорость. Также в состав MKL входит библиотека векторных операций VML (например, расчет корня квадратного от элементов массива) в параллельном режиме. Эффективность MKL на процессорах Intel столь высока, что популярные математические пакеты, такие как MATLAB и MAPLE, используют MKL, выполняя функцию оболочки библиотеки.

Нельзя не отметить SCALAPACK — версию LAPACK, поддерживающую параллельные вычисления. Для её компиляции необходима версия MPICH. SCALAPACK приводит к примерно кратному увеличению производительности при выполнении стандартных функций библиотеки пропорционально количеству ядер.

РЕНСИТ | 2011 | ТОМ 3 | НОМЕР 1

74 АФАНАСЬЕВ В.П., БУДАК В.П.,

ЕФРЕМЕНКО Д.С., ЛУБЕНЧЕНКО А.В.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

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

• определение номера текущего процесса (MPI_ COMM_RANK)

• определение общего количества процессов (MPI_COMM_SIZE)

• Блокирующие и неблокирующие операции рассылки и сбора данных (MPI_(I)SEND, MPI_(I)Recv)

• Барьер (MPI_Barrier)

Технология MPI идеально подходит для распараллеливания вычислений по азимутальным гармоникам.

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

Существует несколько форматов хранения разреженных матриц: compressed sparse row (CSR) format, compressed sparse column format (CSC), Coordinate Format. Тот или иной формат хранения предпочтительней в зависимости от вида операции. Например, для умножения на разреженную матрицу более всего удобен формат CSC, так как элементы в этом случае упорядочены по каждому столбцу. Например, матрица

B =

( 1 0 0 0

0 0 7 0

-3 0 0 4

^0 5 0 2

в формате CSC представима в виде четырех массивов: Values = [1, -3, 5, 7, 4, 2], Rows = [1, 3, 4, 2, 3, 4], Pointer B = [1, 3, 4, 5], Pointer E = [3, 4, 5, 7].

Тот или иной формат хранения предпочтительней в зависимости от вида операции. Например, для умножения на разреженную матрицу более 1

всего удобен формат CSC, так как элементы в этом случае упорядочены по каждому столбцу.

В последнее время активно используется технология вычислений на графических процессорах. Необходимо отметить получившую распространение технологию CUDA. MATLAB 2011a поддерживает работу с CUDA на видеокартах с совместимостью выше 2.0.

7. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Алгоритм, описанный в п. 5.3, был запрограммирован в виде расчетного кода MVDOM на языке FORTRAN и в среде Matlab. В таблице 1 представлено сравнение времен счета для двух тестов (тест 1 — N = 101, K=500, M=32; тест 2 — N = 101, K=1000, M=32). Тесты производились на системе Ubuntu 10.04, Intel Core 2 Duo 3GHz, 2 Gb RAM, Intel Fortran Compiler 11.1 with MKL 10.2. Использовались два компилятора (gfortran и ifort), различные режимы оптимизации и аппарат разреженных матриц.

Заметим, что MKL использует в расчетах все доступные ядра процессора (в нашем случае 2), что приводит к двукратному ускорению. Существенно позволяет сократить время счета аппарат разреженных матриц, используемый для вычисления анизотропной части (20). На рис. 1 приведена зависимость времени счета от K. Разреженные матрицы представляют двухмерные массивы в одномерные, поэтому время счета расчет линейно, а не квадратично.

Сравнение компиляторов ifort и gfortran показывает, что их эффективность примерно одинакова с небольшим преимуществом ifort.

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

Таблица 1

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

Features Runtime I, sec Runtime II, sec

gfortran + LAPACK 240 530

gfortran + LAPACK + optimization 230 505

ifort + LAPACK 210 490

ifort + MKL 115 250

ifort + MKL + optimization 105 230

ifort + MKL + optimization + sparse matrix 33 44

MATLAB 2010a 27 45

1 НОМЕР | ТОМ 3 | 2011 | РЕНСИТ

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

ВЛИЯНИЕ АППАРАТНО-ПРОГРАММНЫХ 75 СРЕДСТВ НА СКОРОСТЬ ВЫЧИСЛЕНИЯ

Рис. 1. Зависимость времени счета от K. программами. Но определение точности расчета является тяжелой задачей, так как определяется рядом факторов, например, вычислением сферических функций. Минимальное расхождение между программами MDOM и DISORT составляет 10-4. При больших M расхождение только увеличивается. Это можно объяснить малыми значениями сферических функций при больших M, и как следствие большой относительной ошибкой, возникающей при их рекуррентном вычислении.

Код DISORT доступен по адресу [ftp://dimate1. gsfc.nasa.gov/wiscombe/Multiple Scatt]. Интересно отметить, что время счета может быть снижено на порядок с помощью небольшой модификации процедуры ZEROIT. В случае достаточно гладких индикатрис (g<0.9) время счета MDOM и DISORT практически одинаковое. В DISORT2.0 реализованы дельта-M метод и TMS метод. Для индикатрисы [20] оказалось, что метод выделения анизотропной части в малоугловом приближении оказался существенно более эффективным. Результаты сравнения программ представлены на рис. 2. Скорость счета MDOM (10 секунд) при этом на порядок выше, чем у DISORT (100 секунд).

т=10.0, ©O=45.0o, Л=1.00

a

b

Рис. 2. Сравнение MDOM и DISORT: а) отраженное и b) прошедшее излучение.

Таблица 2

Сравнение компиляторов

Компилятор и библиотека Время счета, сек.

ifort + mkl 1.00

ifort+ lapack 1.86

qfortran + lapack 2.12

f95 + lapack 2.34

Результаты сравнения компиляторов представлены в таблице 2.

Схожие результаты получены и в векторном случае — сравнивались коды PSTAR и MVDOM для случая фазовая матрица Water Cloud C1, в0=85°, т0=10.0. Результаты расчета представлены на рис. 3. Интересно отметить, что с ростом толщины слоя растет расхождение в прошедшем излучении. Времена счета для MVDOM составляют 50 секунд, для PSTAR —150 секунд.

^ X 10 3 „X10 6

b

Рис. 3. Сравнение MVDOM и Pstar: a) отраженное и b) прошедшее излучение

РЕНСИТ | 2011 | ТОМ 3 | НОМЕР 1

76 АФАНАСЬЕВ В.П., БУДАК В.П.,

ЕФРЕМЕНКО Д.С., ЛУБЕНЧЕНКО А.В.

С учетом аналитического выделения анизотропной части на основе МСГ занимаемая память программой составляет около 100 Мб, то есть объем заведомо больше кэша процессора. В настоящее время наиболее доступными средствами распараллеливания вычислений являются MPI, OpenMP и CUDA.

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

При использовании стандарта OpenMP главный процесс создает несколько (в данном случае два) процессов. Такой подход применительно к рассматриваемому алгоритму не имеет преимуществ, так как эти процессы работают с одними и теми же массивами, а очевидные места кода (например, поэлементное присваивание A(1:N) = B(1:N)+C) компилятор ifort может распараллелить самостоятельно.

В статье [21] обсуждается проблема оптимального распределения данных внутри памяти видеокарты с тем, чтобы скорость передачи данных между процессами была наибольшей. В случае реализации алгоритма решения ВУПИ для слоя , основанного на операциях линейной алгебры, исполь-

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

реализации на ней библиотек BLAS и LAPACK. В этом случае оптимальное распределение данных берет на себя библиотека CUDA. ограничение на использование “своих” подпрограмм и действий над данными в среде MATLAB заставили нас совершать процедуру “пересылка в GPU - расчет -пересылка в CPU”. Использование CUDA в среде MATLAB привело к ускорению на 20 %. Столь незначительное ускорение по сравнению с кратным выигрышем по скорости вычислений по методу Монте-Карло объясняется двумя причинами. Наиболее трудоемкая процедура — вычисление собственных векторов — не поддерживается в среде MATLAB на графическом профессоре. Применение же метода выделения анизотропной части привело к существенному уменьшению размера матриц, участвующих в решении. Поэтому выигрыш от использования CUDA получился не столь большим.

8. ЗАКЛЮЧЕНИЕ

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

ЛИТЕРАТУРА

1. Doicu A, Schreier F, Hilgers S, Bargen A, Slijkhuis S, Hess M. An efficient inversion algorithm for atmospheric remote sensing with application to UV limb observations. J Quant Spectrosc Radiat Transf, 2007, 103:193-208.

2. Spurr RJD, Kurosu TP, Chance KV. A linearized discrete ordinate radiative transfer model for atmospheric remote-sensing retrieval. J Quant Spectrosc Radiat Transf 2001, 68:689-735.

3. Yokota T, Oguma HM, Morino I, Inoue G. A nadirlooking ‘SWIR’ sensor to monitor CO2 column density for Japanese GOSAT project. Proceedings of the XXIV international symposium on space tech science, Miyazaki, Japan, 2004, p. 887.

4. Nakajima T, Tanaka M. Algorithms for radiative intensity calculations in moderately thick atmos using a trun-cation approximation. J Quant Spectrosc Radiat Transf, 1988, 40:51-69.

1 НОМЕР | ТОМ 3 | 2011 | РЕНСИТ

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

5. Budak VP, Klyuykov DA, Korkin SV. Complete matrix solution of radiative transfer equation for pile of horizontally homogeneous slabs. JQuant Spectrosc Radiat Transf 2011, 112:1141-1148.

6. Wang MC, Guth E. On the theory of multiple scattering, particularly of charged particles. Phys.Rev, 1951, 84:1092-111.

7. Krylov VI. Approximate calculation of integrals. New York, Macmillan, 1962, 276 p.

8. Eddington AS. On the radiative equilibrium of the stars. Month Not R Astroph Soc, 1916, 77:16-35.

9. Milne EA. The reflection effect in eclipse binaries. Mon Not R Astroph Soc, 1926, 87:43-55.

10. Kuscer I, Ribaric M. Matrix formalism in the theory of diffusion of light. Opt Acta, 1959, 6:42-51.

11 .Budak VP, Korkin SV On the solution of a vectorial radiative transfer equation in an arbitrary three-dimensional turbid medium with anisotropic scattering. J Quant Spectrosc Radiat Transf, 2008, 109:220-34.

12.Budak VP, Klyuykov DA, Korkin SV Convergence acceleration of radiative transfer equation solution at strongly anisotropic scattering. In: Light Scattering Reviews 5: Single Light Scattering and Radiative Transfer. Ed. A.A. Kokhanovsky. Berlin, Springer Praxis Books, 2010, pp. 147-204.

13.Siewert CE. A discrete-ordinates solution for radiative-transfer models that include polarization effects. J Quant Spectrosc Radiat Transf, 2000, 64:227-254.

14.Sykes JB. Approximate integration of the equation of transfer. Month Not RAstroph Soc, 1951 (111):378—386.

15. Karp AH, Greenstadt J, Fillmore JA. Radiative transfer through an arbitrary thick scattering atmosphere. J Quant Spectrosc Radiat Transf, 1980, 24:391-406.

16. Sobolev VV A Treatise on Radiative Transfer. Princeton, NJ, Van Nostrand, 1963, 171 p.

17. Kisselev V Peaked phase function approximation in the solution of radiative transfer equation. Proc of SPIE, 2005, 5829:63-73.

18. Budak VP, Kozelskii AV, Savitskii EN. Improvement of the spherical harmonics method convergence at strongly anisotropic scattering. Atm Ocean Opt J, 2004, 17:28-33.

19. Press W Teukolsky S, Vetterling W Flannery B. Numerical Recipes in Fortran 77: The Art of Scientific Computing. Cambridge: University Press, 2007, 818 p.

20. Mobley CD, Sundman LK, Boss E. Phase function effects on oceanic light fields. Appl. Opt., 2002, 41:10351050.

21. Peng K, Gao X, Qu X, Ren N, Chen X, He X, Wang X, Liang J, Tian J. Graphics processing unit parallel accelerated solution of the discrete ordinates for photon transport in biological tissues. Appl. Opt, 2011, 50:3808-3823.

ВЛИЯНИЕ АППАРАТНО-ПРОГРАММНЫХ 77 СРЕДСТВ НА СКОРОСТЬ ВЫЧИСЛЕНИЯ

Афанасьев Виктор Петрович,

МЭИ (ТУ), кафедра общей физики и ядерного синтеза, 111250 Москва, Красноказарменная, 14, тел. +7 903 222 0899, [email protected]

Будак Владимир Павлович,

д.т.н, профессор,

МЭИ (ТУ), ИРЭ, кафедра светотехники,

111250 Москва, Красноказарменная, 13, к.Е, тел. +7 985 763 5239, [email protected]

Ефременко Дмитрий Сергеевич,

МЭИ (ТУ), кафедра общей физики и ядерного синтеза, 111250 Москва, Красноказарменная, 14, тел. +7 905 521 6815, [email protected]

Лубенченко Александр Владимирович,

МЭИ (ТУ), кафедра общей физики и ядерного синтеза,

111250 Москва, Красноказарменная, 14, тел. +7 905 541 0495, [email protected]

РЕНСИТ | 2011 | ТОМ 3 | НОМЕР 1

78

information technologies

HARD-SOFTWARE INFLUENCE ON CALCULATING SPEED OF RADIATIVE TRANSFER EQUATION SOLUTION ALGORITHM IN MEDIUM SLAB

Afanasiev V. P., Budak V. P., Efremenko D. S., Lubenchenko A. V.

Moscow Power Engineering Institute (Technical University), Krasnokazarmennaya, 13, 111250 Moscow, Russian Federation [email protected]

It is shown that in slab geometry the discretized radiative transfer equation (RTE) has the unique analytical solution. Since the solution represents as the linear matrix equation, where there are optimized libraries, the only possible algorithm solution exists. The RTE discretization is probable only after the solution anisotropic part selection, including all the singularities. The various implementations of the RTE solution algorithms differ by methods of the anisotropic part selection, among which one the most effective is the small angle modification of spherical harmonics method. The effect of hard-software on the implementation efficiency of RTE solution algorithms in the slab is analyzed.

Keywords: radiative transfer equation, slab geometry, anisotropic scattering

PACS 42.68.Ay

Bibliography — 21 references

RENSIT, 2011, 3(1):69-78______________________________

REFERENCES

1. Doicu A, Schreier F, Hilgers S, Bargen A, Slijkhuis S, Hess M. J Quant Spectrosc Radiat Transf, 2007, 103(2):193-208.

2. Spurr RJD, Kurosu TP, Chance KV. J Quant Spectrosc Radiat Transf 2001, 68:689-735.

3. Yokota T, Oguma HM, Morino I, Inoue G. Proc. 24th Intern. Symp. on Space Tech. Science. Miyazaki, Japan, 2004, p. 887.

4. Nakajima T, Tanaka M. J Quant Spectrosc Radiat Transf 1988, 40:51-69.

5. Budak VP, Klyuykov DA, Korkin SV. J Quant Spectrosc Radiat Transf 2011, 112:1141-1148.

6. Wang MC, Guth E. Phys.Rev, 1951, 84:1092-1110.

7. Krylov VI. Approximate calculation of integrals. New York, Macmillan, 1962, 276 p.

8. Eddington AS. Month Not R Astroph Soc, 1916, 77:1635.

9. Milne EA. Mon Not R Astroph Soc, 1926, 87:43-55.

10. Kuscer I, Ribaric M. Opt Acta, 1959, 6:42-51.

11. Budak VP, Korkin SV J Quant Spectrosc Radiat Transf, 2008, 109:220-34.

12. Budak VP, Klyuykov DA, Korkin SV In: Light Scattering Reviews 5: Single Light Scattering and Radiative Transfer. Ed. A.A. Kokhanovsky. Berlin, Springer Praxis Books, 2010, pp. 147204.

13.Siewert CE. J Quant Spectrosc Radiat Transf, 2000, 64:227254. 1

Received 17.05.2011, revised 20.05.2011

14.Sykes JB. Month Not RAstroph Soc, 1951, 111:378-386.

15. Karp AH, Greenstadt J, Fillmore JA. J Quant Spectrosc Radiat Transf 1980, 24:391-406.

16. Sobolev VV. A Treatise on Radiative Tranfer. Princeton, NJ, Van Nostrand, 1963, 171 p.

17. Kisselev V Proc of SPIE, 2005, (5829):63-73.

18. Budak VP, Kozelskii AV, Savitskii EN. Atm Ocean Opt J, 2004, 17:28-33.

19. Press W Teukolsky S, Vetterling W Flannery B.

Numerical Recipes in Fortran 77: The Art of Scientific Computing. Cambridge, University Press, 2007, 818 p.

20. Mobley CD, Sundman LK, Boss E. Appl Opt, 2002, 41:1035-1050.

21. Peng K, Gao X, Qu X, Ren N, Chen X, He X, Wang X, Liang J, Tian J. Appl Opt, 2011, 50:38083823.

1 НОМЕР | ТОМ 3 | 2011 | РЕНСИТ