arma-thesis-ru.org (347312B)
1 #+TITLE: Имитационное моделирование нерегулярного волнения для программ динамики морских объектов 2 #+AUTHOR: Ганкевич Иван Геннадьевич 3 #+DATE: Санкт-Петербург, 2017 4 #+LANGUAGE: ru 5 #+SETUPFILE: setup.org 6 #+LATEX_HEADER_EXTRA: \organization{Санкт-Петербургский государственный университет} 7 #+LATEX_HEADER_EXTRA: \manuscript{на правах рукописи} 8 #+LATEX_HEADER_EXTRA: \degree{Диссертация на соискание ученой степени\\кандидата физико-математических наук} 9 #+LATEX_HEADER_EXTRA: \speciality{Специальность 05.13.18\\Математическое моделирование, численные методы и комплексы программ} 10 #+LATEX_HEADER_EXTRA: \supervisor{Научный руководитель\\д.т.н. Дегтярев Александр Борисович} 11 #+LATEX_HEADER_EXTRA: \newcites{published}{Список опубликованных по теме диссертации работ} 12 13 #+latex: \maketitlepage 14 15 * Введение 16 **** Актуальность темы. 17 Программы, моделирующие поведение судна на морских волнах, широко применяются 18 для расчета качки судна, оценки величины воздействия внешних сил на плавучую 19 платформу или другой морской объект, а также для оценки вероятности 20 опрокидывания судна при заданных погодных условиях; однако, большинство из них 21 используют линейную теорию для моделирования морского 22 волнения\nbsp{}cite:shin2003nonlinear,van2007forensic,kat2001prediction,van2002development, 23 в рамках которой сложно воспроизвести определенные особенности ветро-волнового 24 климата. Среди них можно выделить переход от нормальных погодных условий к 25 шторму и волнение, вызванное наложением множества систем ветровых волн и волн 26 зыби, распространяющихся в нескольких направлениях. Другой недостаток линейной 27 теории волн заключается в предположении, что высота волн много меньше их длины. 28 Это делает расчеты грубыми при моделировании качки судна в условиях 29 нерегулярного волнения, когда такое предположение несправедливо. Разработка 30 новых и более совершенных моделей и методов, используемых в программах расчета 31 динамики судна, увеличит количество сценариев применения таких программ и, в 32 частности, поспособствует исследованию поведения судна в экстремальных условиях. 33 34 **** Степень разработанности. 35 Модель авторегрессии скользящего среднего (АРСС) является ответом на сложности, 36 с которыми на практике сталкиваются ученые, использующие в своей работе модели 37 морского волнения, разработанные в рамках линейной теории волн. Проблемы, с 38 которыми они сталкиваются при использовании модели Лонге---Хиггинса (которая 39 полностью основана на линейной теории волн) перечислены ниже. 40 1. /Периодичность/. В рамках линейной теории волны аппроксимируются суммой 41 гармоник, из-за чего период реализации взволнованной поверхности зависит от 42 их количества. Чем больше размер реализации, тем больше коэффициентов 43 требуется для исключения периодичности, поэтому с увеличением размера 44 реализации время ее генерации растет нелинейно. Это приводит к тому, что 45 любая модель, основанная на линейной теории, неэффективна при генерации 46 больших реализаций взволнованной поверхности, независимо от того, насколько 47 оптимизирован исходный код программы. 48 2. /Линейность/. В рамках линейной теории волн дается математическое определение 49 морским волнам в предположении малости их амплитуд по сравнению с длинами. 50 Такие волны, в основном, характерны для открытого океана, а волны в 51 прибрежных районах и штормовые волны, для которых это предположение 52 несправедливо, грубо описываются в рамках линейной теории. 53 3. /Вероятностная сходимость/. Фаза волны, значение которой получается с помощью 54 генератора псевдослучайных чисел (ГПСЧ), имеет равномерное распределение, что 55 приводит к медленной сходимости интегральных характеристик взволнованной 56 поверхности (распределение высот волн, их периодов, длин и т.п.). 57 58 Эти сложности стали отправной точкой в поиске модели, не основанной на линейной 59 теории волн, и в исследованиях процесса АРСС был найден необходимый 60 математический аппарат. 61 1. Входным параметром процесса АРСС является автоковариационная функция (АКФ), 62 которая может быть напрямую получена из энергетического или 63 частотно-направленного спектра морского волнения (который, в свою очередь 64 является входным параметром для модели Лонге---Хиггинса). Так что входные 65 параметры одной модели могут быть легко преобразованы во входные параметры 66 другой. 67 2. Процесс АРСС не имеет ограничение на амплитуду генерируемых волн: их крутизна 68 может быть увеличена на столько, на сколько это позволяет АКФ реальных 69 морских волн. 70 3. Период реализации равен периоду ГПСЧ, поэтому время генерации растет линейно 71 с увеличением размера реализации. 72 4. Белый шум, который является единственным вероятностным членом формулы 73 процесса АРСС, имеет нормальное распределение, из-за чего скорость сходимости 74 не носит вероятностный характер. 75 76 **** Цели и задачи. 77 Процесс АРСС является основой модели морского волнения АРСС, но, ввиду своей 78 нефизической природы, модель нуждается в доработке перед тем, как ее можно 79 использовать для генерации взволнованной поверхности. 80 1. Исследовать, как различные формы АКФ влияют на выбор параметров процесса АРСС 81 (количество коэффициентов процесса скользящего среднего и процесса 82 авторегрессии). 83 2. Исследовать возможность генерации волн с произвольными профилями, а не только 84 с профилем синусоиды (учесть асимметричность распределения волновых аппликат 85 взволнованной поверхности). 86 3. Разработать метод для определения поля давлений под дискретно заданной 87 взволнованной поверхностью. Такие формулы обычно выводятся для конкретной 88 модели путем подстановки формулы профиля волны в систему уравнений для 89 определения давлений\nbsp{}eqref:eq-problem, однако процесс АРСС не содержит 90 в себе формулу профиля волны в явном виде, поэтому для него необходимо 91 получить решение для взволнованной поверхности общего вида (для которой не 92 существует математического выражения) без линеаризации граничных условий (ГУ) 93 и предположении о малости амплитуд волн. 94 4. Верифицировать соответствие интегральных характеристик взволнованной 95 поверхности реальным морским волнам. 96 5. Разработать комплекс программ, реализующий созданную модель и метод расчета 97 давлений и позволяющий проводить расчеты на вычислительной системе как с 98 общей, так и распределенной памятью. 99 100 **** Научная новизна. 101 Модель АРСС в отличие от других моделей ветрового волнения не основана на 102 линейной теории волн, что позволяет 103 - генерировать волны произвольной амплитуды, регулируя крутизну посредством АКФ; 104 - генерировать волны произвольной формы, регулируя асимметричность распределения 105 волновых аппликат посредством нелинейного безынерционного преобразования 106 (НБП). 107 В то же время математический аппарат процесса АРСС хорошо изучен в других 108 научных областях, что позволяет его обобщить для моделирования развития морского 109 волнения в условиях шторма с учетом климатических спектров и данных ассимиляции 110 определенных районов мирового океана, что невозможно сделать с помощью модели, 111 основанной на линейной теории волн. 112 113 Отличительными особенностями данной работы являются 114 - использование во всех экспериментах /трехмерных/ моделей АР и СС, 115 - использование метода вычисления поля потенциала скорости, работающего с 116 дискретно заданной взволнованной поверхностью, и 117 - разработка комплекса программ, реализующего генерацию взволнованной морской 118 поверхности и вычисления поля давлений на системах с общей (SMP) и 119 распределенной (MPP) памятью, а также гибридных системах (GPGPU), использующих 120 графические сопроцессоры для ускорения вычислений. 121 122 **** Теоретическая и практическая значимость работы. 123 Применение модели АРСС и формулы поля давлений, не использующей предположения 124 линейной теории волн, качественно повысит работу комплексов программ для расчета 125 воздействия океанских волн на морские объекты. 126 127 1. Поскольку метод для поля давлений выводится для дискретно заданной 128 взволнованной поверхности и без каких-либо предположений об амплитудах волн, 129 то он применим для любой взволнованной поверхности невязкой несжимаемой 130 жидкости (в частности он применим для поверхности, генерируемой моделью 131 Лонге---Хиггинса). Это позволяет использовать этот метод без привязки к 132 модели АРСС. 133 2. С вычислительной точки зрения этот метод более эффективен, чем 134 соответствующий метод для модели ЛХ, поскольку интегралы в формуле сводятся к 135 преобразованиям Фурье, для которых существует семейство алгоритмов быстрого 136 преобразования Фурье (БПФ), оптимизированных под разные архитектуры 137 процессоров. 138 3. Поскольку используемая в методе формула явная, то обмена данными между 139 параллельными процессами можно избежать, что позволяет масштабировать 140 производительность на большое количество процессорных ядер. 141 4. Наконец, сама модель АРСС более эффективна, чем модель ЛХ, ввиду отсутствия 142 тригонометрических функций в ее формуле. Взволнованная поверхность 143 вычисляется как сумма большого числа многочленов, для которых существует 144 низкоуровневая ассемблерная инструкция FMA (Fused Multiply-Add), а шаблон 145 доступа к памяти позволяет эффективно использовать кэш центрального 146 процессора. 147 148 **** Методология и методы исследования. 149 Программная реализация модели АРСС и метода вычисления поля давлений создавалась 150 поэтапно: прототип, написанный на высокоуровневом инженерных языках 151 (Mathematica\nbsp{}cite:mathematica10 и Octave\nbsp{}cite:octave2015), был 152 преобразован в программу на языке более низкого уровня (C++). Ввиду 153 использования различных абстракций и языковых примитивов реализация одних и тех 154 же алгоритмов и методов на языках разного уровня позволяет выявить и исправить 155 ошибки, которые остались бы незамеченными в случае использования лишь одного 156 языка. Генерируемая моделью АРСС взволнованная поверхность, а также все входные 157 параметры (АКФ, формула распределения волновых аппликат и т.п.) были проверены с 158 помощью встроенных в язык программирования графических средств. 159 160 **** Положения, выносимые на защиту. 161 - Имитационная модель морского волнения, способная генерировать реализации 162 взволнованной морской поверхности, имеющие большой период и состоящие из волн 163 произвольной амплитуды; 164 - Метод вычисления поля давлений, разработанный для этой модели без 165 предположений линейной теории волн; 166 - Программная реализация имитационной модели и метода для вычислительных систем 167 с общей и распределенной памятью. 168 **** Степень достоверности и апробация результатов. 169 Верификация модели АРСС проводится путем сравнения интегральных характеристик 170 (распределений волновых аппликат, высот и длин волн и т.п.) генерируемой 171 взволнованной поверхности с характеристиками реальных морских волн. Метод 172 вычисления поля давлений был разработан с помощью языка Mathematica, в котором 173 полученные выражения проверяются с помощью встроенных в язык графических 174 средств. 175 176 Модель АРСС и метод вычислений поля давлений были реализованы в Large Amplitude 177 Motion Programme (LAMP), программе для моделирования качки судна, и сопоставлены 178 с используемой ранее моделью ЛХ. Численные эксперименты показали более высокую 179 вычислительную эффективность модели АРСС. 180 181 * Постановка задачи 182 Задача состоит в применении математического аппарата процесса АРСС для 183 моделирования морских волн и в разработке метода вычисления поля давлений под 184 дискретно заданной взволнованной морской поверхностью для случая идеальной 185 несжимаемой жидкости без предположений линейной теории волн. 186 - Для случая волн малых амплитуд полученное по новому методу поле давлений 187 должно быть сопоставимо с полем, полученным по формуле из линейной теории 188 волн; для остальных случаев значения поля не должны стремиться к 189 бесконечности. 190 - Интегральные характеристики генерируемой взволнованной поверхности должны 191 совпадать с характеристиками реальных морских волн. 192 - Программная реализация модели АРСС и метода вычисления давлений должна 193 работать на системах с общей и распределенной памятью. 194 195 **** Формула для поля давлений. 196 Задача определения поля давлений под взволнованной морской поверхностью 197 представляет собой обратную задачу гидродинамики для несжимаемой невязкой 198 жидкости. Система уравнений для нее в общем виде записывается 199 как\nbsp{}cite:kochin1966theoretical 200 \begin{align} 201 & \nabla^2\phi = 0,\nonumber\\ 202 & \phi_t+\frac{1}{2} |\vec{\upsilon}|^2 + g\zeta=-\frac{p}{\rho}, & \text{на }z=\zeta(x,y,t),\label{eq-problem}\\ 203 & D\zeta = \nabla \phi \cdot \vec{n}, & \text{на }z=\zeta(x,y,t),\nonumber 204 \end{align} 205 где \(\phi\)\nbsp{}--- потенциал скорости, \(\zeta\)\nbsp{}--- подъем 206 (аппликата) взволнованной поверхности, \(p\)\nbsp{}--- давление жидкости, 207 \(\rho\)\nbsp{}--- плотность жидкости, 208 \(\vec{\upsilon}=(\phi_x,\phi_y,\phi_z)\)\nbsp{}--- вектор скорости, 209 \(g\)\nbsp{}--- ускорение свободного падения и \(D\)\nbsp{}--- субстанциональная 210 производная (производная Лагранжа). Первое уравнение является уравнением 211 неразрывности (уравнение Лапласа), второе\nbsp{}--- законом сохранения импульса 212 (которое иногда называют динамическим граничным условием); третье 213 уравнение\nbsp{}--- кинематическое граничное условие, которое сводится к 214 равенству скорости перемещения этой поверхности (\(D\zeta\)) нормальной 215 составляющей скорости жидкости (\(\nabla\phi\cdot\vec{n}\), 216 см.\nbsp{}разд.\nbsp{}[[#directional-derivative]]). 217 218 Обратная задача гидродинамики заключается в решении этой системы уравнений 219 относительно \(\phi\). В такой постановке динамическое ГУ становится явной 220 формулой для определения поля давлений по значениям производных потенциалов 221 скорости, получаемых из оставшихся уравнений. С математической точки зрения 222 обратная задача гидродинамики сводится к решению уравнения Лапласа со смешанным 223 ГУ\nbsp{}--- задаче Робена. 224 225 Обратная задача возможна, поскольку модель АРСС генерирует гидродинамически 226 адекватную взволнованную морскую поверхность: распределения интегральных 227 характеристик и дисперсионное соотношение соответствуют реальным 228 волнам\nbsp{}cite:boukhanovsky1997thesis,degtyarev2011modelling. 229 230 * Модель АРСС в задаче имитационного моделирования морского волнения 231 ** Анализ моделей морского волнения 232 Вычисление давлений возможно только при условии знания формы взволнованной 233 поверхности, которая задается либо дискретно в каждой точке пространственной 234 сетки, либо непрерывно с помощью аналитической формулы. Как будет показано в 235 разделе\nbsp{}[[#linearisation]], знание такой формулы может упростить вычисление 236 давлений, фактически сведя задачу к генерации поля давлений, а не самой 237 взволнованной поверхности. 238 239 **** Модель Лонге---Хиггинса. 240 Наиболее простой моделью, формула которой выводится в рамках линейной теории 241 волн (см.\nbsp{}разд.\nbsp{}[[#longuet-higgins-derivation]]), является модель 242 Лонге---Хиггинса (ЛХ)\nbsp{}cite:longuet1957statistical. Подробный сравнительный 243 анализ этой модели и модели АРСС проведен в 244 работах\nbsp{}cite:degtyarev2011modelling,boukhanovsky1997thesis. 245 246 Модель ЛХ представляет взволнованную морскую поверхность в виде суперпозиции 247 элементарных гармонических волн случайных амплитуд \(c_n\) и фаз \(\epsilon_n\), 248 равномерно распределенных на интервале \([0,2\pi]\). Подъем (координата \(z\)) 249 поверхности определяется формулой 250 \begin{equation} 251 \zeta(x,y,t) = \sum\limits_n c_n \cos(u_n x + v_n y - \omega_n t + \epsilon_n). 252 \label{eq-longuet-higgins} 253 \end{equation} 254 Здесь волновые числа \((u_n,v_n)\) непрерывно распределены на плоскости 255 \((u,v)\), т.е.\nbsp{}площадка \(du\times{}dv\) содержит бесконечно большое 256 количество волновых чисел. Частота связана с волновыми числами дисперсионным 257 соотношением \(\omega_n=\omega(u_n,v_n)\). Функция \(\zeta(x,y,t)\) является 258 трехмерным эргодическим стационарным однородным гауссовым процессом, 259 определяемым соотношением 260 \begin{equation*} 261 2E_\zeta(u,v)\, du\, dv = \sum\limits_n c_n^2, 262 \end{equation*} 263 где \(E_\zeta(u,v)\)\nbsp{}--- двухмерная спектральная плотность энергии волн. 264 Коэффициенты \(c_n\) определяются из энергетического спектра волнения 265 \(S(\omega)\) по формуле 266 \begin{equation*} 267 c_n = \sqrt{ \textstyle\int\limits_{\omega_n}^{\omega_{n+1}} S(\omega) d\omega}. 268 \end{equation*} 269 270 **** Недостатки модели Лонге---Хиггинса. 271 Несмотря на то что модель ЛХ отличается простотой численного алгоритма, на 272 практике она обладает рядом недостатков. 273 1. Модель рассчитана на представление стационарного гауссова поля. Это является 274 следствием центральной предельной теоремы (ЦПТ): сумма большого числа 275 гармоник со случайными амплитудами и фазами имеет нормальное распределение в 276 независимости от спектра, подаваемого на вход модели. Использование меньшего 277 количества коэффициентов может решить проблему, но также уменьшит период 278 реализации. Использование модели ЛХ для генерации волн с негауссовым 279 распределением аппликат (которое имеют реальные морские 280 волны\nbsp{}cite:huang1980experimental,rozhkov1996theory) не реализуемо на 281 практике. 282 2. С вычислительной точки зрения, недостатком модели является нелинейный рост 283 времени генерации поверхности с увеличением размера реализации. Чем больше 284 размер реализации, тем больше коэффициентов (дискретных точек 285 частотно-направленного спектра) требуется для исключения периодичности. Это 286 делает модель неэффективной для проведения длительных численных 287 экспериментов. 288 3. Наконец, с инженерной точки зрения, модель обладает рядом особенностей, 289 которые не позволяют использовать ее в качестве фундамента для построения 290 более совершенных моделей. 291 - В программной реализации скорость сходимости 292 выражения\nbsp{}[[eqref:eq-longuet-higgins]] низка, т.к.\nbsp{}фазы 293 \(\epsilon_n\) имеют распределены равномерно. 294 - Обобщение модели для негауссовых и нелинейных процессов затруднено, 295 поскольку требует включения нелинейных членов 296 в\nbsp{}[[eqref:eq-longuet-higgins]], для которых не известна формула 297 вычисления коэффициентов\nbsp{}cite:rozhkov1990stochastic. 298 299 Таким образом, модель ЛХ применима для решения задачи генерации взволнованной 300 морской поверхности только в рамках линейной теории волн, неэффективна для 301 длительных экспериментов и имеет ряд недостатков, не позволяющих использовать ее 302 в качестве основы для построения более совершенных моделей. 303 304 **** Модель АРСС. 305 В\nbsp{}cite:spanos1982arma модель АРСС используется для генерации временного 306 ряда, спектр которого совпадает со спектром Пирсона---Московица (ПМ). Авторы 307 проводят эксперименты для одномерных моделей АР, СС и АРСС. Они отмечают 308 превосходное совпадение полученного и исходного спектров и более высокую 309 вычислительную эффективность модели АРСС по сравнению с моделями, основанными на 310 суммировании большого числа гармоник со случайными фазами. Также отмечается, что 311 для того чтобы спектр полученного временного ряда совпадал с заданным, модели СС 312 требуется меньшее количество коэффициентов, чем модели АР. 313 В\nbsp{}cite:spanos1996efficient автор обобщает формулы для нахождения 314 коэффициентов модели АРСС для случая нескольких переменных (векторов). 315 316 Отличие данной работы от вышеперечисленных состоит в исследовании трехмерной 317 модели АРСС (два пространственных и одно временное измерение), что во многом 318 является другой задачей. 319 1. Система уравнений Юла---Уокера, используемая для определения коэффициентов 320 АР, имеет более сложную блочно-блочную структуру. 321 2. Оптимальный (для совпадения заданного и исходного спектров) порядок модели 322 определяется вручную. 323 3. Вместо аппроксимации ПМ в качестве входа модели используются аналитические 324 выражения для АКФ стоячих и прогрессивных волн. 325 4. Трехмерная взволнованная поверхность должна быть сопоставима с реальной 326 морской поверхностью не только по спектральным характеристикам, но и по форме 327 волновых профилей, поэтому верификация модели производится и для 328 распределений различных параметров генерируемых волн (длин, высот, периодов и 329 др.). 330 Многомерность исследуемой модели не только усложняет задачу, но и позволяет 331 провести визуальную проверку генерируемой взволнованной поверхности. Именно 332 возможность визуализировать результат работы программы позволяет удостовериться, 333 что генерируемая поверхность действительно похожа на реальное морское волнение, 334 а не является абстрактным многомерным случайным процессом, совпадающим с 335 реальным лишь статистически. 336 337 В\nbsp{}cite:fusco2010short модель АР используется для прогнозирования волн зыби 338 для управления преобразователем энергии волн (ПЭВ) в реальном времени. Для 339 эффективной работы ПЭВ необходимо, чтобы частота встроенного осциллятора 340 совпадала с частотой морских волн. Авторы статьи представляют подъем волны как 341 временной ряд и сравнивают эффективность модели АР, нейронных сетей и 342 циклических моделей в прогнозировании будущих значений ряда. Модель АР дает 343 наиболее точный прогноз для низкочастотных волн зыби вплоть до двух типовых 344 периодов волн. Это пример успешного применения процесса АР для моделирования 345 морских волн. 346 347 ** Основные формулы трехмерного процесса АРСС 348 Модель АРСС для морского волнения определяет взволнованную морскую поверхность 349 как трехмерный (два пространственных и одно временное измерение) процесс 350 авторегрессии скользящего среднего: каждая точка взволнованной поверхности 351 представляется в виде взвешенной суммы предыдущих по времени и пространству 352 точек и взвешенной суммы предыдущих по времени и пространству нормально 353 распределенных случайных импульсов. Основным уравнением для трехмерного процесса 354 АРСС является 355 \begin{equation} 356 \zeta_{\vec i} 357 = 358 \sum\limits_{\vec j = \vec 0}^{\vec N} 359 \Phi_{\vec j} \zeta_{\vec i - \vec j} 360 + 361 \sum\limits_{\vec j = \vec 0}^{\vec M} 362 \Theta_{\vec j} \epsilon_{\vec i - \vec j} 363 , 364 \label{eq-arma-process} 365 \end{equation} 366 где \(\zeta\)\nbsp{}--- подъем (аппликата) взволнованной поверхности, 367 \(\Phi\)\nbsp{}--- коэффициенты процесса АР, \(\Theta\)\nbsp{}--- коэффициенты 368 процесса СС, \(\epsilon\)\nbsp{}--- белый шум, имеющий Гауссово распределение, 369 \(\vec{N}\)\nbsp{}--- порядок процесса АР, \(\vec{M}\)\nbsp{}--- порядок 370 процесса СС, причем \(\Phi_{\vec{0}}\equiv0\), \(\Theta_{\vec{0}}\equiv0\). 371 Здесь стрелки обозначают многокомпонентные индексы, содержащие отдельную 372 компоненту для каждого измерения. В общем случае в качестве компонент могут 373 выступать любые скалярные величины (температура, соленость, концентрация 374 какого-либо раствора в воде и т.п.). Параметрами уравнения служат коэффициенты и 375 порядки процессов АР и СС. 376 377 Имитационная модель морского волнения предназначена для воспроизведения 378 реалистичных реализаций ветро-волнового поля и пригодна для использования в 379 расчетах динамики судна. Свойства стационарности и обратимости являются 380 основными критериями выбора того или иного процесса для моделирования волн 381 разных профилей, которые обсуждаются в разделе\nbsp{}[[#sec-process-selection]]. 382 383 **** Процесс авторегрессии (АР). 384 Процесс АР\nbsp{}--- это процесс АРСС только лишь с одним случайным импульсом 385 вместо их взвешенной суммы: 386 \begin{equation} 387 \zeta_{\vec i} 388 = 389 \sum\limits_{\vec j = \vec 0}^{\vec N} 390 \Phi_{\vec j} \zeta_{\vec i - \vec j} 391 + 392 \epsilon_{i,j,k} 393 . 394 \label{eq-ar-process} 395 \end{equation} 396 Коэффициенты авторегрессии \(\Phi\) определяются из многомерных уравнений 397 Юла---Уокера (ЮУ), получаемых после домножения на \(\zeta_{\vec{i}-\vec{k}}\) 398 обеих частей уравнения и взятия математического ожидания. В общем виде уравнения 399 ЮУ записываются как 400 \begin{equation} 401 \label{eq-yule-walker} 402 \gamma_{\vec k} 403 = 404 \sum\limits_{\vec j = \vec 0}^{\vec N} 405 \Phi_{\vec j} 406 \text{ }\gamma_{\vec{k}-\vec{j}} 407 + 408 \Var{\epsilon} \delta_{\vec{k}}, 409 \qquad 410 \delta_{\vec{k}} = 411 \begin{cases} 412 1, \quad \text{if } \vec{k}=0 \\ 413 0, \quad \text{if } \vec{k}\neq0, 414 \end{cases} 415 \end{equation} 416 где \(\gamma\)\nbsp{}--- АКФ процесса \(\zeta\), \(\Var{\epsilon}\)\nbsp{}--- 417 дисперсия белого шума. Матричная форма трехмерной системы уравнений ЮУ, 418 используемой в данной работе, имеет следующий вид. 419 \begin{equation*} 420 \Gamma 421 \left[ 422 \begin{array}{l} 423 \Phi_{\vec 0}\\ 424 \Phi_{0,0,1}\\ 425 \vdotswithin{\Phi_{\vec 0}}\\ 426 \Phi_{\vec N} 427 \end{array} 428 \right] 429 = 430 \left[ 431 \begin{array}{l} 432 \gamma_{0,0,0}-\Var{\epsilon}\\ 433 \gamma_{0,0,1}\\ 434 \vdotswithin{\gamma_{\vec 0}}\\ 435 \gamma_{\vec N} 436 \end{array} 437 \right], 438 \qquad 439 \Gamma= 440 \left[ 441 \begin{array}{llll} 442 \Gamma_0 & \Gamma_1 & \cdots & \Gamma_{N_1} \\ 443 \Gamma_1 & \Gamma_0 & \ddots & \vdotswithin{\Gamma_0} \\ 444 \vdotswithin{\Gamma_0} & \ddots & \ddots & \Gamma_1 \\ 445 \Gamma_{N_1} & \cdots & \Gamma_1 & \Gamma_0 446 \end{array} 447 \right], 448 \end{equation*} 449 где \(\vec N = \left( N_1, N_2, N_3 \right)\) и 450 \begin{equation*} 451 \Gamma_i = 452 \left[ 453 \begin{array}{llll} 454 \Gamma^0_i & \Gamma^1_i & \cdots & \Gamma^{N_2}_i \\ 455 \Gamma^1_i & \Gamma^0_i & \ddots & \vdotswithin{\Gamma^0_i} \\ 456 \vdotswithin{\Gamma^0_i} & \ddots & \ddots & \Gamma^1_i \\ 457 \Gamma^{N_2}_i & \cdots & \Gamma^1_i & \Gamma^0_i 458 \end{array} 459 \right] 460 \qquad 461 \Gamma_i^j= 462 \left[ 463 \begin{array}{llll} 464 \gamma_{i,j,0} & \gamma_{i,j,1} & \cdots & \gamma_{i,j,N_3} \\ 465 \gamma_{i,j,1} & \gamma_{i,j,0} & \ddots & \vdotswithin{\gamma_{i,j,0}} \\ 466 \vdotswithin{\gamma_{i,j,0}} & \ddots & \ddots & \gamma_{i,j,1} \\ 467 \gamma_{i,j,N_3} & \cdots & \gamma_{i,j,1} & \gamma_{i,j,0} 468 \end{array} 469 \right], 470 \end{equation*} 471 Поскольку по определению \(\Phi_{\vec 0}\equiv0\), то первую строку и столбец 472 матрицы \(\Gamma\) можно отбросить. Матрица \(\Gamma\), как и оставшаяся от нее 473 матрица, будут блочно-блочно-теплицевы, положительно определены и симметричны, 474 поэтому систему уравнений можно эффективно решить методом Холецкого, специально 475 предназначенного для таких матриц. 476 477 После нахождения решения системы уравнений дисперсия белого шума определяется из 478 уравнения\nbsp{}eqref:eq-yule-walker при \(\vec k = \vec 0\) как 479 \begin{equation*} 480 \Var{\epsilon} = 481 \Var{\zeta} 482 - 483 \sum\limits_{\vec j = \vec 0}^{\vec N} 484 \Phi_{\vec j} 485 \text{ }\gamma_{\vec{j}}. 486 \end{equation*} 487 488 **** Процесс скользящего среднего (СС). 489 Процесс СС\nbsp{}--- это процесс АРСС, в котором \(\Phi\equiv0\): 490 \begin{equation} 491 \zeta_{\vec i} 492 = 493 \sum\limits_{\vec j = \vec 0}^{\vec M} 494 \Theta_{\vec j} \epsilon_{\vec i - \vec j} 495 . 496 \label{eq-ma-process} 497 \end{equation} 498 Коэффициенты СС \(\Theta\) определяются неявно из системы нелинейных уравнений 499 \begin{equation*} 500 \gamma_{\vec i} = 501 \left[ 502 \displaystyle 503 \sum\limits_{\vec j = \vec i}^{\vec M} 504 \Theta_{\vec j}\Theta_{\vec j - \vec i} 505 \right] 506 \Var{\epsilon}. 507 \end{equation*} 508 Система решается численно с помощью метода простой итерации по формуле 509 \begin{equation*} 510 \Theta_{\vec i} = 511 -\frac{\gamma_{\vec 0}}{\Var{\epsilon}} 512 + 513 \sum\limits_{\vec j = \vec i}^{\vec M} 514 \Theta_{\vec j} \Theta_{\vec j - \vec i}. 515 \end{equation*} 516 Здесь новые значения коэффициентов \(\Theta\) вычисляются, начиная с последнего: 517 от \(\vec{i}=\vec{M}\) до \(\vec{i}=\vec{0}\). Дисперсия белого шума вычисляется 518 из 519 \begin{equation*} 520 \Var{\epsilon} = \frac{\gamma_{\vec 0}}{ 521 1 522 + 523 \sum\limits_{\vec j = \vec 0}^{\vec M} 524 \Theta_{\vec j}^2 525 }. 526 \end{equation*} 527 Авторы\nbsp{}cite:box1976time предлагают использовать метод Ньютона---Рафсона 528 для решения этого уравнения с большей точностью, однако, в данном случае этот 529 метод сложно обобщить для трех измерений. Использование более медленного метода 530 не оказывает большого эффекта на общую производительность программы, потому что 531 количество коэффициентов мало, и большую часть времени программа тратит на 532 генерацию взволнованной поверхности. 533 534 **** Стационарность и обратимость процессов АР и СС. 535 Для того чтобы моделируемая взволнованная поверхность представляла собой 536 физическое явление, соответствующий процесс должен быть стационарным и 537 обратимым. Если процесс /обратим/, то существует разумная связь текущих событий 538 с событиями в прошлом, и, если процесс /стационарен/, то амплитуда моделируемого 539 физического сигнала не увеличивается бесконечно в пространстве и времени. 540 541 Процесс АР всегда обратим, а для стационарности необходимо, чтобы корни 542 характеристического уравнения 543 \begin{equation*} 544 1 - \Phi_{0,0,1} z - \Phi_{0,0,2} z^2 545 - \cdots 546 - \Phi_{\vec N} z^{N_0 N_1 N_2} = 0, 547 \end{equation*} 548 лежали \emph{вне} единичного круга. Здесь \(\vec{N}\)\nbsp{}--- порядок процесса 549 АР, а \(\Phi\)\nbsp{}--- коэффициенты. 550 551 Процесс СС всегда стационарен, а для обратимости необходимо, чтобы корни 552 характеристического уравнения 553 \begin{equation*} 554 1 - \Theta_{0,0,1} z - \Theta_{0,0,2} z^2 555 - \cdots 556 - \Theta_{\vec M} z^{M_0 M_1 M_2} = 0, 557 \end{equation*} 558 лежали \emph{вне} единичного круга. Здесь \(\vec{M}\)\nbsp{}--- порядок процесса 559 СС, а \(\Theta\)\nbsp{}--- коэффициенты. 560 561 Свойства стационарности и обратимости являются основными критериями при выборе 562 процесса для моделирования разных профилей волн, которые обсуждаются в 563 разделе\nbsp{}[[#sec-process-selection]]. 564 565 **** Смешанный процесс авторегрессии скользящего среднего (АРСС). 566 :PROPERTIES: 567 :CUSTOM_ID: sec-how-to-mix-arma 568 :END: 569 570 В общем и целом, процесс АРСС получается путем подстановки сгенерированной 571 процессом СС взволнованной поверхности в качестве случайного импульса процесса 572 АР, однако, для того чтобы АКФ результирующего процесса соответствовала 573 заданной, необходимо предварительно скорректировать значения коэффициентов АР. 574 Существует несколько способов "смешивания" процессов АР и СС. 575 - Подход, предложенный авторами\nbsp{}cite:box1976time, который включает в себя 576 разделение АКФ на часть для процесса АР и часть для процесса СС по каждому из 577 измерений, не подходит в данной ситуации, поскольку в трех измерениях 578 невозможно таким образом разделить АКФ: всегда останутся части, которые не 579 будут учтены ни в процессе АР, ни в процессе СС. 580 - Альтернативный подход состоит в использовании одной и той же (неразделенной) 581 АКФ для процессов АР и СС процессов разных порядков, однако, тогда 582 характеристики реализации (математической ожидание, дисперсия и др.) будут 583 смещены: они станут характеристиками двух наложенных друг на друга процессов. 584 Для первого подхода авторами\nbsp{}cite:box1976time предложена формула 585 корректировки коэффициентов процесса АР, для второго же подхода такой формулы 586 нет. Таким образом, единственным проверенным решением на данный момент является 587 использование процессов АР и СС по отдельности. 588 589 **** Критерии выбора процесса. 590 :PROPERTIES: 591 :CUSTOM_ID: sec-process-selection 592 :END: 593 594 Одной из проблем в применении модели АРСС для генерации взволнованной морской 595 поверхности является то, что для разных профилей волн /необходимо/ использовать 596 разные процессы: стоячие волны моделируются только процессом АР, а прогрессивные 597 волны\nbsp{}--- только процессом СС. Это утверждение пришло из практики: если 598 попытаться использовать процессы наоборот, результирующая реализация либо 599 расходится, либо не представляет собой реальные морские волны (такое происходит 600 в случае необратимого процесса СС, который всегда стационарен). Таким образом, 601 процесс АР может быть использован только для моделирования стоячих волн, а 602 процесс СС\nbsp{}--- для прогрессивных волн. 603 604 Другой проблемой является сложность определения оптимального количества 605 коэффициентов для трехмерных процессов АР и СС. Для одномерных процессов 606 существуют простые в реализации итеративные методы\nbsp{}cite:box1976time, а для 607 трехмерного случая существуют аналогичные 608 методы\nbsp{}cite:byoungseon1999arma3d,liew2005arma3d, которые сложны в 609 реализации, из-за чего не были использованы в данной работе. Ручной выбор 610 порядка модели тривиален: задание заведомо высокого порядка сказывается на 611 производительности, но и не на качестве реализации, поскольку периодичность 612 процессов не зависит от количества коэффициентов. Программная реализация 613 итеративного метода поиска оптимального порядка модели повысила бы уровень 614 автоматизации генератора взволнованной поверхности, но не качество проводимых 615 экспериментов. 616 617 Практика показывает, что некоторые утверждения авторов\nbsp{}cite:box1976time не 618 выполняются для трехмерной модели АРСС. Например, авторы утверждают, что АКФ 619 процесса СС обрывается на отсчете \(q\), а АКФ процесса АР затухает на 620 бесконечности, однако, на практике при использовании слабо затухающей и 621 обрывающейся на отсчете \(q\) АКФ для трехмерного процесса СС получается 622 необратимый процесс СС и реализация, не соответствующая реальными морским 623 волнам, в то время как при использовании той же самой АКФ для трехмерного 624 процесса АР получается стационарный обратимый процесс и адекватная реализация. 625 Также, авторы утверждают, что первые \(q\) точек АКФ смешанного процесса 626 необходимо выделить процессу СС (поскольку он обычно используется для описания 627 пиков АКФ) и отдать остальные точки процессу АР, однако, на практике в случае 628 АКФ прогрессивной волны процесс АР стационарен только для начального временного 629 среза АКФ, а остальные точки отдаются процессу СС. 630 631 Суммируя вышесказанное, наиболее разработанным сценарием применения модели АРСС 632 для генерации взволнованной морской поверхности является использование процесса 633 АР для стоячих волн и процесса СС для прогрессивных волн. Использование 634 смешанного процесса АРСС для обоих типов волн может сделать модель более точной 635 при условии наличия соответствующих формул корректировки коэффициентов, что 636 является целью дальнейших исследований. 637 638 **** Верификация интегральных характеристик взволнованной поверхности. 639 Для модели АР в 640 работах\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis 641 экспериментальным путем были верифицированы 642 - распределения различных характеристик волн (высоты волн, длины волн, длины 643 гребней, период волн, уклон волн, показатель трехмерности), 644 - дисперсионное соотношение, 645 - сохранение интегральных характеристик для случая смешанного волнения. 646 В данной работе верифицируются как модель АР, так и СС путем сравнения 647 распределений различных характеристик волн. 648 649 В\nbsp{}cite:rozhkov1990stochastic авторы показывают, что некоторые 650 характеристики морских волн (перечисленные в табл.\nbsp{}[[tab-weibull-shape]]) 651 имеют распределение Вейбулла, а подъем взволнованной поверхности\nbsp{}--- 652 нормальное распределение. Для верификации генерируемых моделями АР и СС 653 реализаций используются спрямленные диаграммы (графики, в которых по оси \(OX\) 654 откладываются квантили функции распределения, вычисленные аналитически, а по оси 655 \(OY\), вычисленные экспериментально). Если экспериментально полученное 656 распределение соответствует аналитическому, то график представляет собой прямую 657 линию. Концы графика могут отклоняться от прямой линии, поскольку не могут быть 658 надежно получены из реализации конечной длины. 659 660 #+name: tab-weibull-shape 661 #+caption[Значение коэффициента формы распределения Вейбулла]: 662 #+caption: Значение коэффициента формы распределения Вейбулла для различных 663 #+caption: характеристик волн. 664 #+attr_latex: :booktabs t 665 | Характеристика | Коэффициент формы \(k\) | 666 |-------------------------+-------------------------| 667 | | <l> | 668 | Высота волны | 2 | 669 | Длина волны | 2,3 | 670 | Длина гребня волны | 2,3 | 671 | Период волны | 3 | 672 | Уклон волны | 2,5 | 673 | Показатель трехмерности | 2,5 | 674 675 Верификация производится для стоячих и прогрессивных волн. Соответствующие АКФ и 676 спрямленные диаграммы распределений характеристик волн представлены на 677 рис.\nbsp{}[[acf-slices]],\nbsp{}[[standing-wave-distributions]],\nbsp{}[[propagating-wave-distributions]]. 678 679 #+name: propagating-wave-distributions 680 #+begin_src R :file build/propagating-wave-qqplots-ru.pdf 681 source(file.path("R", "common.R")) 682 par(pty="s", mfrow=c(2, 2), family="serif") 683 arma.qqplot_grid( 684 file.path("build", "arma-benchmarks", "verification-orig", "propagating_wave"), 685 c("elevation", "heights_y", "lengths_y", "periods"), 686 c("подъем", "высота по Y", "длина по Y", "период"), 687 xlab="x", 688 ylab="y" 689 ) 690 #+end_src 691 692 #+caption[Спрямленные диаграммы для прогрессивных волн]: 693 #+caption: Спрямленные диаграммы для прогрессивных волн. 694 #+name: propagating-wave-distributions 695 #+RESULTS: propagating-wave-distributions 696 [[file:build/propagating-wave-qqplots-ru.pdf]] 697 698 #+name: standing-wave-distributions 699 #+begin_src R :file build/standing-wave-qqplots-ru.pdf 700 source(file.path("R", "common.R")) 701 par(pty="s", mfrow=c(2, 2), family="serif") 702 arma.qqplot_grid( 703 file.path("build", "arma-benchmarks", "verification-orig", "standing_wave"), 704 c("elevation", "heights_y", "lengths_y", "periods"), 705 c("подъем", "высота по Y ", "длина по Y", "период"), 706 xlab="x", 707 ylab="y" 708 ) 709 #+end_src 710 711 #+caption[Спрямленные диаграммы для стоячих волн]: 712 #+caption: Спрямленные диаграммы для стоячих волн. 713 #+name: standing-wave-distributions 714 #+RESULTS: standing-wave-distributions 715 [[file:build/standing-wave-qqplots-ru.pdf]] 716 717 #+name: acf-slices 718 #+header: :width 6 :height 9 719 #+begin_src R :file build/acf-slices-ru.pdf 720 source(file.path("R", "common.R")) 721 propagating_acf <- read.csv(file.path("build", "arma-benchmarks", "verification-orig", "propagating_wave", "acf.csv")) 722 standing_acf <- read.csv(file.path("build", "arma-benchmarks", "verification-orig", "standing_wave", "acf.csv")) 723 par(mfrow=c(5, 2), mar=c(0,0,0,0), family="serif") 724 for (i in seq(0, 4)) { 725 arma.wavy_plot(standing_acf, i, zlim=c(-5,5)) 726 arma.wavy_plot(propagating_acf, i, zlim=c(-5,5)) 727 } 728 #+end_src 729 730 #+caption[Временные срезы АКФ для стоячих и прогрессивных волн]: 731 #+caption: Временные срезы АКФ для стоячих (слева) и прогрессивных (справа) волн. 732 #+name: acf-slices 733 #+RESULTS: acf-slices 734 [[file:build/acf-slices-ru.pdf]] 735 736 Хвосты распределений на рис.\nbsp{}[[propagating-wave-distributions]] отклоняются от 737 оригинального распределения для характеристик отдельных волн, поскольку каждую 738 волну необходимо извлечь из полученной взволнованной поверхности, чтобы измерить 739 ее длину, период и высоту. Алгоритм, который бы гарантировал безошибочное 740 извлечение всех волн, не известен, поскольку волны могут и часто накладываются 741 друг на друга. Правый хвост распределения Вейбулла отклоняется больше, поскольку 742 он представляет редко возникающие волны. 743 744 Степень соответствия для стоячих волн (рис.\nbsp{}[[standing-wave-distributions]]) 745 ниже для высот и длин, примерно одинакова для подъема поверхности и выше для 746 периодов волн. Более низкая степень соответствия длин и высот может быть 747 результатом того, что распределения были получены эмпирически для морских волн, 748 которые, в основном, являются прогрессивными, и аналогичные распределения для 749 стоячих волн могут отличаться. Более высокая степень соответствия периодов волн 750 является следствием того, что периоды стоячих волн извлекаются более точно, 751 поскольку волн не перемещаются вне моделируемой области взволнованной 752 поверхности. Одинаковая степень соответствия для подъема поверхности получается 753 из-за того, что это характеристика поверхности (и соответствующего процесса АР 754 или СС), и она не зависит от типа волн. 755 756 Таким образом, программная реализация моделей АР и СС (которая подробно описана 757 в разд.\nbsp{}[[#sec-hpc]]) генерирует взволнованную поверхность, распределения 758 характеристик отдельных волн которой соответствуют натурным данным. 759 760 ** Моделирование нелинейности морских волн 761 Модель АРСС позволяет учесть асимметричность распределения волновых аппликат, 762 т.е.\nbsp{}генерировать морские волны, закон распределения аппликат которых 763 имеет ненулевой эксцесс и асимметрию. Такой закон распределения характерен для 764 реальных морских волн\nbsp{}cite:longuet1963nonlinear и задается либо 765 полиномиальной аппроксимацией натурных данных, либо аналитически. 766 Асимметричность волн моделируется с помощью нелинейного безынерционного 767 преобразования (НБП) случайного процесса, однако, любое нелинейное 768 преобразование случайного процесса приводит к преобразованию его АКФ. Для того 769 чтобы подавить этот эффект, необходимо предварительно преобразовать АКФ, как 770 показано в\nbsp{}cite:boukhanovsky1997thesis. 771 772 **** Преобразование взволнованной поверхности. 773 Формула \(z=f(y)\) преобразования взволнованной поверхности к необходимому 774 одномерному закону распределения \(F(z)\) получается путем решения нелинейного 775 трансцендентного уравнения \(F(z) = \Phi(y)\), где \(\Phi(y)\)\nbsp{}--- функция 776 одномерного нормального закона распределения. Поскольку функция распределения 777 аппликат морских волн часто задается некоторой аппроксимацией, основанной на 778 натурных данных, то это уравнение целесообразно решать численно в каждой точке 779 \(y_k|_{k=0}^N\) сетки сгенерированной поверхности относительно \(z_k\). Тогда 780 уравнение запишется в виде 781 \begin{equation} 782 \label{eq-distribution-transformation} 783 F(z_k) 784 = 785 \frac{1}{\sqrt{2\pi}} 786 \int\limits_0^{y_k} \exp\left[ -\frac{t^2}{2} \right] dt 787 . 788 \end{equation} 789 Поскольку функции распределения монотонны, для решения этого уравнения 790 используется простейший численный метод половинного деления (метод бисекции). 791 792 **** Предварительное преобразование АКФ. 793 Для преобразования АКФ \(\gamma_z\) процесса ее необходимо разложить в ряд по 794 полиномам Эрмита (ряд Грама---Шарлье) 795 \begin{equation*} 796 \gamma_z \left( \vec u \right) 797 = 798 \sum\limits_{m=0}^{\infty} 799 C_m^2 \frac{\gamma_y^m \left( \vec u \right)}{m!}, 800 \end{equation*} 801 где 802 \begin{equation*} 803 C_m = \frac{1}{\sqrt{2\pi}} 804 \int\limits_{0}^\infty 805 f(y) H_m(y) \exp\left[ -\frac{y^2}{2} \right], 806 \end{equation*} 807 \(H_m\)\nbsp{}--- полином Эрмита, а \(f(y)\)\nbsp{}--- решение 808 уравнения\nbsp{}eqref:eq-distribution-transformation. Воспользовавшись 809 полиномиальной аппроксимацией \(f(y) \approx \sum\limits_i d_i y^i\) и 810 аналитическими выражениями для полиномов Эрмита, формулу определения 811 коэффициентов можно упростить, используя следующее равенство: 812 \begin{equation*} 813 \frac{1}{\sqrt{2\pi}} 814 \int\limits_\infty^\infty 815 y^k \exp\left[ -\frac{y^2}{2} \right] 816 = 817 \begin{cases} 818 (k-1)!! & \text{для четных }k,\\ 819 0 & \text{для нечетных }k. 820 \end{cases} 821 \end{equation*} 822 Оптимальное количество коэффициентов \(C_m\) определяется путем вычисления их 823 последовательно и критерий прекращения счета определяется совпадением дисперсий 824 обоих полей с требуемой точностью \(\epsilon\): 825 \begin{equation} 826 \label{eq-nit-error} 827 \left| \Var{z} - \sum\limits_{k=0}^m 828 \frac{C_k^2}{k!} \right| \leq \epsilon. 829 \end{equation} 830 831 В\nbsp{}cite:boukhanovsky1997thesis автор предлагает использовать полиномиальную 832 аппроксимацию для \(f(y)\) также для преобразования поверхности, однако на 833 практике в реализации взволнованной поверхности часто находятся точки, 834 выпадающие за промежуток на котором построена аппроксимация, что приводит к 835 резкому уменьшению ее точности. В этих точках 836 уравнение\nbsp{}eqref:eq-distribution-transformation эффективнее решать методом 837 бисекции. Использование полиномиальной аппроксимации в формулах для 838 коэффициентов ряда Грама---Шарлье не приводит к аналогичным ошибкам. 839 840 **** Разложение в ряд Грама---Шарлье. 841 В\nbsp{}cite:huang1980experimental было экспериментально показано, что 842 распределение аппликат морской поверхности отличается от нормального ненулевым 843 эксцессом и асимметрией. В\nbsp{}cite:rozhkov1996theory показано, что такое 844 распределение раскладывается в ряд Грама---Шарлье (РГШ): 845 \begin{align} 846 \label{eq-skew-normal-1} 847 & F(z; \mu=0, \sigma=1, \gamma_1, \gamma_2) \approx 848 \Phi(z; \mu, \sigma) \frac{2 + \gamma_2}{2} 849 - \frac{2}{3} \phi(z; \mu, \sigma) 850 \left(\gamma_2 z^3+\gamma_1 851 \left(2 z^2+1\right)\right) \nonumber 852 \\ 853 & f(z; \gamma_1, \gamma_2) \approx 854 \frac{1}{\sigma\sqrt{2 \pi}}e^{-\frac{(z-\mu)^2}{2\sigma^2}} 855 \left[ 856 1+ 857 \frac{1}{6} \gamma_1 z \left(z^2-3\right) 858 + \frac{1}{24} \gamma_2 \left(z^4-6z^2+3\right) 859 \right], 860 \end{align} 861 где \(\Phi(z)\)\nbsp{}--- функция распределения (ФР) нормального закона, 862 \(\phi\)\nbsp{}--- ФПР нормального закона, \(\gamma_1\)\nbsp{}--- асимметрия, 863 \(\gamma_2\)\nbsp{}--- эксцесс, \(f\)\nbsp{}--- ФПР, \(F\)\nbsp{}--- функция 864 распределения (ФР). Согласно\nbsp{}cite:rozhkov1990stochastic для аппликат 865 морских волн значение асимметрии выбирается на интервале 866 \(0,1\leq\gamma_1\leq{0,52}]\), а значение эксцесса на интервале 867 \(0,1\leq\gamma_2\leq{0,7}\). Семейство плотностей распределения при различных 868 параметрах показано на рис.\nbsp{}[[fig-skew-normal-1]]. 869 870 #+name: fig-skew-normal-1 871 #+begin_src R :file build/skew-normal-1-ru.pdf 872 source(file.path("R", "common.R")) 873 par(family="serif") 874 x <- seq(-3, 3, length.out=100) 875 params <- data.frame( 876 skewness = c(0.00, 0.52, 0.00, 0.52), 877 kurtosis = c(0.00, 0.00, 0.70, 0.70), 878 linetypes = c("solid", "dashed", "dotdash", "dotted") 879 ) 880 arma.skew_normal_1_plot(x, params) 881 legend( 882 "topleft", 883 mapply( 884 function (s, k) { 885 as.expression(bquote(list( 886 gamma[1] == .(arma.fmt(s, 2)), 887 gamma[2] == .(arma.fmt(k, 2)) 888 ))) 889 }, 890 params$skewness, 891 params$kurtosis 892 ), 893 lty = paste(params$linetypes) 894 ) 895 #+end_src 896 897 #+caption[График плотности распределения на основе РГШ]: 898 #+caption: График плотности распределения закона, основанного на РГШ, 899 #+caption: при различных значениях асимметрии \(\gamma_1\) и эксцесса 900 #+caption: \(\gamma_2\). 901 #+name: fig-skew-normal-1 902 #+RESULTS: fig-skew-normal-1 903 [[file:build/skew-normal-1-ru.pdf]] 904 905 **** Асимметричное нормальное распределение. 906 Альтернативной аппроксимацией распределения волновых аппликат служит формула 907 асимметричного нормального распределения (АНР): 908 \begin{align} 909 \label{eq-skew-normal-2} 910 F(z; \alpha) & = \frac{1}{2} 911 \mathrm{erfc}\left[-\frac{z}{\sqrt{2}}\right]-2 T(z,\alpha ), \nonumber \\ 912 f(z; \alpha) & = \frac{e^{-\frac{z^2}{2}}}{\sqrt{2 \pi }} 913 \mathrm{erfc}\left[-\frac{\alpha z}{\sqrt{2}}\right], 914 \end{align} 915 где \(T\)\nbsp{}--- функция Оуэна\nbsp{}cite:owen1956tables. Эта формула не 916 позволяет задать значения асимметрии и эксцесса по отдельности\nbsp{}--- оба 917 значения регулируются параметром \(\alpha\). Преимущество данной формулы лишь в 918 относительной простоте вычисления: эта функция встроена в некоторые программы и 919 библиотеки математических функций. График функции для разных значений \(\alpha\) 920 представлен на рис.\nbsp{}[[fig-skew-normal-2]]. 921 922 #+name: fig-skew-normal-2 923 #+begin_src R :file build/skew-normal-2-ru.pdf 924 source(file.path("R", "common.R")) 925 par(family="serif") 926 x <- seq(-3, 3, length.out=100) 927 alpha <- c(0.00, 0.87, 2.25, 4.90) 928 params <- data.frame( 929 alpha = alpha, 930 skewness = arma.bits.skewness_2(alpha), 931 kurtosis = arma.bits.kurtosis_2(alpha), 932 linetypes = c("solid", "dashed", "dotdash", "dotted") 933 ) 934 arma.skew_normal_2_plot(x, params) 935 legend( 936 "topleft", 937 mapply( 938 function (a, s, k) { 939 as.expression(bquote(list( 940 alpha == .(arma.fmt(a, 2)), 941 gamma[1] == .(arma.fmt(s, 2)), 942 gamma[2] == .(arma.fmt(k, 2)) 943 ))) 944 }, 945 params$alpha, 946 params$skewness, 947 params$kurtosis 948 ), 949 lty = paste(params$linetypes) 950 ) 951 #+end_src 952 953 #+caption[График плотности АНР]: 954 #+caption: График плотности АНР при различных значениях коэффициента 955 #+caption: асимметрии \(\alpha\). 956 #+name: fig-skew-normal-2 957 #+RESULTS: fig-skew-normal-2 958 [[file:build/skew-normal-2-ru.pdf]] 959 960 **** Тестирование. 961 Для того чтобы измерить влияние НБП на форму результирующей взволнованной 962 поверхности, было сгенерировано три реализации: 963 - реализация с Гауссовым распределением (без НБП), 964 - реализация с распределением на основе ряда Грама---Шарлье (РГШ), и 965 - реализация с асимметричным нормальным распределением (АНР). 966 Начальные состояния ГПСЧ были заданы одинаковыми для всех запусков программы, 967 чтобы модель АРСС выдавала одни и те же значения для каждой реализации. Было 968 проведено два эксперимента: для стоячих и прогрессивных волн с АКФ, заданными 969 формулами из раздела\nbsp{}[[#sec-wave-acfs]]. 970 971 В то время как эксперимент показал, что применение НБП с распределением РГШ 972 увеличивает крутизну волн, то же самое нельзя сказать об асимметричном 973 нормальном распределении (рис.\nbsp{}[[fig-nit]]). Использование этого распределения 974 приводит к взволнованной поверхности, в которой аппликаты всегда больше или 975 равны нулю. Таким образом, асимметричное нормальное распределение не подходит 976 для НБП. НБП увеличивает высоту и крутизну как прогрессивных, так и стоячих 977 волн. Увеличение асимметрии или эксцесса РГШ приводит в увеличению как высоты, 978 так и крутизны волн. Ошибка аппроксимации АКФ (ур.\nbsp{}eqref:eq-nit-error) 979 принимает значения от 0,20 для РГШ до 0,70 для АНР (табл.\nbsp{}[[tab-nit-error]]). 980 981 #+name: fig-nit 982 #+header: :width 7 :height 7 983 #+begin_src R :file build/nit-ru.pdf 984 source(file.path("R", "nonlinear.R")) 985 par(mfrow=c(2, 1), mar=c(4,4,4,0.5), family='serif') 986 args <- list( 987 graphs=c('Гауссово', 'РГШ', 'АНР'), 988 linetypes=c('solid', 'dashed', 'dotted') 989 ) 990 args$title <- 'Прогрессивные волны' 991 arma.plot_nonlinear(file.path("build", "arma-benchmarks", "verification-orig", "nit-propagating"), args) 992 args$title <- 'Стоячие волны' 993 arma.plot_nonlinear(file.path("build", "arma-benchmarks", "verification-orig", "nit-standing"), args) 994 #+end_src 995 996 #+name: fig-nit 997 #+caption[Срезы поверхности с различными распределениями аппликат]: 998 #+caption: Срезы взволнованной поверхности с различными распределениями 999 #+caption: волновых аппликат (Гауссово, РГШ и асимметричное нормальное). 1000 #+RESULTS: fig-nit 1001 [[file:build/nit-ru.pdf]] 1002 1003 #+name: tab-nit-error 1004 #+caption[Ошибки аппроксимации АКФ]: 1005 #+caption: Ошибки аппроксимации АКФ (разность дисперсий) для различных 1006 #+caption: распределений волновых аппликат. \(N\)\nbsp{}--- количество 1007 #+caption: коэффициентов аппроксимации АКФ. 1008 #+attr_latex: :booktabs t 1009 | | | <r> | <r> | <r> | <r> | | <r> | 1010 | Тип волн | Распред. | \(\gamma_1\) | \(\gamma_2\) | \(\alpha\) | Ошибка | \(N\) | Высота волн | 1011 |---------------+----------+--------------+--------------+------------+--------+-------+-------------| 1012 | прогрессивные | Гауссово | | | | | | 2,41 | 1013 | прогрессивные | РГШ | 2,25 | 0,4 | | 0,20 | 2 | 2,75 | 1014 | прогрессивные | АНР | | | 1 | 0,70 | 3 | 1,37 | 1015 | стоячие | Гауссово | | | | | | 1,73 | 1016 | стоячие | РГШ | 2,25 | 0,4 | | 0,26 | 2 | 1,96 | 1017 | стоячие | АНР | | | 1 | 0,70 | 3 | 0,94 | 1018 1019 Таким образом, единственный тестовый сценарий, который показал приемлемые 1020 результаты\nbsp{}--- это реализации с распределением на основе РГШ для 1021 прогрессивных и стоячих волн. АНР искажает взволнованную поверхность для обоих 1022 типов волн. Реализации с распределением на основе РГШ характеризуются большой 1023 ошибкой аппроксимации АКФ, что приводит к увеличению высоты волн. Причина 1024 большой ошибки заключается в неточности аппроксимации РГШ, которая не сходится 1025 для всевозможных функций\nbsp{}cite:wallace1958asymptotic. Несмотря на большую 1026 ошибку, изменение высоты волн невелико (табл.\nbsp{}[[tab-nit-error]]). 1027 1028 ** Форма АКФ для разных волновых профилей 1029 :PROPERTIES: 1030 :CUSTOM_ID: sec-wave-acfs 1031 :END: 1032 1033 **** Аналитический метод. 1034 Прямой способ нахождения АКФ, соответствующей заданному профилю морской волны, 1035 заключается в применении теоремы Винера---Хинчина. Согласно этой теореме 1036 автокорреляционная функция \(K\) функции \(\zeta\) равна преобразованию Фурье от 1037 квадрата модуля этой функции: 1038 \begin{equation} 1039 K(t) = \Fourier{\left| \zeta(t) \right|^2}. 1040 \label{eq-wiener-khinchin} 1041 \end{equation} 1042 Если заменить \(\zeta\) на формулу для волнового профиля, то это выражение даст 1043 аналитическую формулу для соответствующей АКФ. 1044 1045 Для трехмерного волнового профиля (два пространственных и одно временное 1046 измерение) аналитическая формула представляет собой многочлен высокой степени, и 1047 ее лучше всего вычислять с помощью программы для символьных вычислений. Затем, 1048 для практического применения она может быть аппроксимирована суперпозицией 1049 экспоненциально затухающих косинусов (именно так выглядит АКФ стационарного 1050 процесса АРСС\nbsp{}cite:box1976time). 1051 1052 **** Эмпирический метод. 1053 Впрочем, для трехмерного случая существует более простой эмпирический метод 1054 нахождения формы АКФ, не требующий использования сложного программного 1055 обеспечения. Известно, что АКФ, представляющая собой суперпозицию 1056 экспоненциально затухающих косинусов, является решением уравнения Стокса для 1057 гравитационных волн\nbsp{}cite:boccotti1983wind. Значит, если в моделируемом 1058 морском волнении важна только форма волны, а не точные ее характеристики, то 1059 заданный волновой профиль можно просто домножить на затухающую экспоненту, чтобы 1060 получить подходящую АКФ. Эта АКФ не отражает параметры волн, такие как высота и 1061 период, зато это открывает возможность моделировать волны определенных форм, 1062 определяя профиль волны дискретно заданной функцией, домножая его на экспоненту 1063 и используя результирующую функцию в качестве АКФ. Таким образом, эмпирический 1064 метод неточен, но более простой по сравнению с применением теоремы 1065 Винера---Хинчина; он, в основном, полезен для тестирования модели АРСС. 1066 1067 **** АКФ стоячей волны. 1068 Профиль трехмерной плоской стоячей волны задается как 1069 \begin{equation} 1070 \zeta(t, x, y) = A \sin (k_x x + k_y y) \sin (\sigma t). 1071 \label{eq-standing-wave} 1072 \end{equation} 1073 Найдем АКФ с помощью аналитического метода. Домножив формулу на затухающую 1074 экспоненту (поскольку преобразование Фурье определено для функции \(f\), для 1075 которой справедливо \(f\underset{x\rightarrow\pm\infty}{\longrightarrow}0\)), 1076 получим 1077 \begin{equation} 1078 \zeta(t, x, y) = 1079 A 1080 \exp\left[-\alpha (|t|+|x|+|y|) \right] 1081 \sin (k_x x + k_y y) \sin (\sigma t). 1082 \label{eq-decaying-standing-wave} 1083 \end{equation} 1084 Затем, применяя трехмерное преобразование Фурье к обоим частям уравнения с 1085 помощью программы для символьных вычислений, получим многочлен высокой степени, 1086 который аппроксимируем выражением 1087 \begin{equation} 1088 K(t,x,y) = 1089 \gamma 1090 \exp\left[-\alpha (|t|+|x|+|y|) \right] 1091 \cos \beta t 1092 \cos \left[ \beta x + \beta y \right]. 1093 \label{eq-standing-wave-acf} 1094 \end{equation} 1095 Таким образом, после применения теоремы Винера---Хинчина получаем исходную 1096 формулу, но с косинусами вместо синусов. Это различие важно, поскольку значение 1097 АКФ в точке \((0,0,0)\) равно дисперсии процесса АРСС, которое при использовании 1098 синусов было бы неверным. 1099 1100 Если попытаться получить ту же самую формулу с помощью эмпирического метода, то 1101 выражение\nbsp{}eqref:eq-decaying-standing-wave необходимо адаптировать для 1102 соответствия\nbsp{}eqref:eq-standing-wave-acf. Это можно осуществить либо, 1103 изменяя фазу синуса, либо заменой синуса на косинус, чтобы сдвинуть максимум 1104 функции в начало координат. 1105 1106 **** АКФ прогрессивной волны. 1107 Профиль трехмерной плоской прогрессивной волны задается как 1108 \begin{equation} 1109 \zeta(t, x, y) = A \cos (\sigma t + k_x x + k_y y). 1110 \label{eq-propagating-wave} 1111 \end{equation} 1112 Для аналитического метода повторение шагов из предыдущих двух параграфов дает 1113 \begin{equation} 1114 K(t,x,y) = 1115 \gamma 1116 \exp\left[-\alpha (|t|+|x|+|y|) \right] 1117 \cos\left[\beta (t+x+y) \right]. 1118 \label{eq-propagating-wave-acf} 1119 \end{equation} 1120 Для эмпирического метода профиль волны можно просто домножить на затухающую 1121 экспоненту, не изменяя положение максимума АКФ (как это требуется для стоячей 1122 волны). 1123 1124 **** Сравнение изученных методов. 1125 Итого, аналитический метод нахождения АКФ морских волн сводится к следующим 1126 шагам. 1127 - Обеспечить затухание выражения для профиля волны на \(\pm\infty\), домножив 1128 его на затухающую экспоненту. 1129 - Взять преобразование Фурье от квадрата модуля получившегося профиля, 1130 воспользовавшись программой для символьных вычислений. 1131 - Аппроксимировать получившийся многочлен подходящим выражением для АКФ. 1132 1133 Два примера этого раздела показывают, что затухающие профили стоячих и 1134 прогрессивных волн схожи по форме с соответствующими АКФ с тем лишь различием, 1135 что максимум АКФ должен быть перенесен в начало координат, чтобы сохранить 1136 дисперсию моделируемого процесса. Применение эмпирического метода нахождения АКФ 1137 сводится к следующим шагам. 1138 - Обеспечить затухание выражения для профиля волны на \(\pm\infty\), домножив его 1139 на затухающую экспоненту. 1140 - Перенести максимум получившейся функции в начало координат, используя свойства 1141 тригонометрических функций для сдвига фазы. 1142 1143 ** Выводы 1144 В силу своей нефизической природы модель АРСС не включает в себя понятие морской 1145 волны; вместо этого она моделирует взволнованную поверхность как единое целое. 1146 Движения отдельных волн и их форма часто получаются грубыми, а точное количество 1147 генерируемых волн неизвестно. Несмотря на это, интегральные характеристики 1148 взволнованной поверхности соответствуют реальным морским волнам. 1149 1150 Теоретически, профили самих морских волн могут быть использованы в качестве АКФ, 1151 если предварительно обеспечить их экспоненциальное затухание. Это может 1152 позволить генерировать волны произвольных профилей, не прибегая к НБП, и 1153 является одной из тем дальнейших исследований. 1154 1155 * Поле давлений под дискретно заданной взволнованной поверхностью 1156 ** Известные формулы определения поля давлений 1157 **** Теория волн малых амплитуд. 1158 В\nbsp{}cite:stab2012,degtyarev1998modelling,degtyarev1997analysis дается 1159 решение обратной задачи гидродинамики для случая идеальной несжимаемой жидкости 1160 в рамках теории волн малых амплитуд (в предположении, что длина волны много 1161 больше ее высоты: \(\lambda \gg h\)). В этом случае обратная задача линейна и 1162 сводится к уравнению Лапласа со смешанным граничным условием, а уравнение 1163 движения используется только для нахождения давлений по известным значениям 1164 производных потенциала скорости. Предположение о малости амплитуд волн означает 1165 слабое изменение локального волнового числа во времени и пространстве по 1166 сравнению с подъемом (аппликатой) взволнованной поверхности. Это позволяет 1167 вычислить производную подъема поверхности по \(z\) как \(\zeta_z=k\zeta\), где 1168 \(k\)\nbsp{}--- волновое число. В двухмерном случае решение записывается явной 1169 формулой 1170 \begin{align} 1171 \left.\frac{\partial\phi}{\partial x}\right|_{x,t}= & 1172 -\frac{1}{\sqrt{1+\alpha^{2}}}e^{-I(x)} 1173 \int\limits_{0}^x\frac{\partial\dot{\zeta}/\partial 1174 z+\alpha\dot{\alpha}}{\sqrt{1+\alpha^{2}}}e^{I(x)}dx,\label{eq-old-sol-2d}\\ 1175 I(x)= & \int\limits_{0}^x\frac{\partial\alpha/\partial z}{1+\alpha^{2}}dx,\nonumber 1176 \end{align} 1177 где \(\alpha\)\nbsp{}--- уклоны волн. В трехмерном случае решение записывается в 1178 виде эллиптического дифференциального уравнения в частных производных (ДУЧП): 1179 \begin{align*} 1180 & \frac{\partial^2 \phi}{\partial x^2} \left( 1 + \alpha_x^2 \right) + 1181 \frac{\partial^2 \phi}{\partial y^2} \left( 1 + \alpha_y^2 \right) + 1182 2\alpha_x\alpha_y \frac{\partial^2 \phi}{\partial x \partial y} + \\ 1183 & \left( 1184 \frac{\partial \alpha_x}{\partial z} + 1185 \alpha_x \frac{\partial \alpha_x}{\partial x} + 1186 \alpha_y \frac{\partial \alpha_x}{\partial y} 1187 \right) \frac{\partial \phi}{\partial x} + \\ 1188 & \left( 1189 \frac{\partial \alpha_y}{\partial z} + 1190 \alpha_x \frac{\partial \alpha_y}{\partial x} + 1191 \alpha_y \frac{\partial \alpha_y}{\partial y} 1192 \right) \frac{\partial \phi}{\partial y} + \\ 1193 & \frac{\partial \dot{\zeta}}{\partial z} + 1194 \alpha_x \dot{\alpha_x} + \alpha_y \dot{\alpha_y} = 0. 1195 \end{align*} 1196 Уравнение предполагается решать численно путем сведения к разностному. 1197 1198 Как будет показано в разделе\nbsp{}[[#sec-compare-formulae]] 1199 формула\nbsp{}eqref:eq-old-sol-2d расходится при попытке вычислить поле 1200 скоростей для волн больших амплитуд, а значит не может быть использована 1201 совместно с моделью морского волнения, генерирующей волны произвольных амплитуд. 1202 1203 **** Линеаризация граничного условия. 1204 :PROPERTIES: 1205 :CUSTOM_ID: linearisation 1206 :END: 1207 1208 Модель ЛХ позволяет вывести явную формулу для поля скоростей путем линеаризации 1209 кинематического граничного условия. Формула для потенциала скорости запишется 1210 как 1211 \begin{equation*} 1212 \phi(x,y,z,t) = \sum_n \frac{c_n g}{\omega_n} 1213 e^{\sqrt{u_n^2+v_n^2} z} 1214 \sin(u_n x + v_n y - \omega_n t + \epsilon_n). 1215 \end{equation*} 1216 Формула дифференцируется для получения производных потенциала, которые 1217 подставляются в динамическое граничное условие для определения поля давлений. 1218 1219 ** Определение поля давлений под дискретно заданной взволнованной поверхностью 1220 Аналитические решения граничных задач для классических уравнений часто 1221 используются для исследования различных свойств уравнений, и для таких 1222 исследований запись формулы общего решения неудобна ввиду своей сложности и 1223 наличия интегралов от неизвестных функций. Одним из методов нахождения 1224 аналитических решений ДУЧП является метод Фурье. Основой метода служит 1225 преобразование Фурье, применение которого к некоторым ДУЧП позволяет свести их к 1226 алгебраическому, а его решение записывается как обратное преобразование Фурье от 1227 некоторой функции (которая может содержать преобразования Фурье от других 1228 функций). Поскольку эти преобразования не всегда можно записать аналитически, то 1229 вместо этого ищутся частные решения задачи и анализируется их поведение в 1230 различных областях. В то же время, численный расчет дискретных преобразований 1231 Фурье возможен для любой дискретно заданной функции, используя алгоритмы БПФ. 1232 Эти алгоритмы используют симметрию комплексных экспонент для понижения 1233 асимптотической сложности с \(\mathcal{O}(n^2)\) до \(\mathcal{O}(n\log_{2}n)\). 1234 Таким образом, даже если общее решение содержит преобразования Фурье от 1235 неизвестных функций, они все равно могут быть взяты численно, а использование 1236 алгоритмов БПФ делает этот подход эффективным. 1237 1238 Альтернативным подходом к решению ДУЧП является их сведение к разностным 1239 уравнениям, решаемым с помощью построения различных численных схем. При этом 1240 решение получается приближенным, а асимптотическая сложность соответствующих 1241 алгоритмов сопоставима со сложностью алгоритма БПФ. Например, для стационарного 1242 эллиптического уравнения в частных производных строится неявная численная схема; 1243 эта схема расчитывается итерационным методом, на каждом шаге которого ищется 1244 решение трехдиагональной или пятидиагональной СЛАУ методом прогонки. 1245 Асимптотическая сложность алгоритма составляет \(\mathcal{O}({n}{m})\), где 1246 \(n\)\nbsp{}--- количество точек на сетке взволнованной поверхности, а 1247 \(m\)\nbsp{}--- число итераций. 1248 1249 Несмотря на широкое распространение, итеративные алгоритмы неэффективно 1250 отображаются на архитектуру параллельных машин ввиду неизбежной синхронизации 1251 процессов после каждой итерации; в частности, отображение на сопроцессоры может 1252 включать в себя копирование данных на сопроцессор и обратно на каждой итерации, 1253 что отрицательно сказывается на их производительности. В то же время, наличие 1254 большого количества преобразований Фурье в решении является скорее 1255 преимуществом, чем недостатком. Во-первых, решения, полученные с помощью метода 1256 Фурье, явные, а значит хорошо масштабируются на большое количество параллельно 1257 работающих вычислительных ядер с использованием простейших приемов параллельного 1258 программирования. Во-вторых, для алгоритмов БПФ существуют готовые 1259 оптимизированные реализации для различных архитектур процессоров и сопроцессоров 1260 (GPU, MIC). Эти преимущества обусловливают использование метода Фурье для 1261 получения явного решения задачи определения давлений под взволнованной морской 1262 поверхностью. 1263 1264 #+LATEX: \pagebreak 1265 *** Двухмерное поле потенциала скорости 1266 :PROPERTIES: 1267 :CUSTOM_ID: sec-pressure-2d 1268 :END: 1269 1270 **** Формула для жидкости бесконечной глубины. 1271 Задача Робена для уравнения Лапласа в двух измерениях записывается как 1272 \begin{align} 1273 \label{eq-problem-2d} 1274 & \phi_{xx}+\phi_{zz}=0,\\ 1275 & \zeta_t + \zeta_x\phi_x 1276 = \FracSqrtZetaX{\zeta_x} \phi_x 1277 - \FracSqrtZetaX{1} \phi_z, 1278 & \text{на }z=\zeta(x,t).\nonumber 1279 \end{align} 1280 Для ее решения воспользуемся методом Фурье. Возьмем преобразование Фурье от 1281 обоих частей уравнений Лапласа и получим 1282 \begin{equation*} 1283 -4 \pi^2 \left( u^2 + v^2 \right) 1284 \FourierY{\phi(x,z)}{u,v} = 0, 1285 \end{equation*} 1286 откуда имеем \(v = \pm i u\). Здесь и далее будет использоваться следующая 1287 симметричная форма преобразования Фурье: 1288 \begin{equation*} 1289 \FourierY{f(x,y)}{u,v} = 1290 \iint\limits_{-\infty}^{\phantom{--}\infty} 1291 f(x,y) 1292 e^{-2\pi i (x u + y v)} 1293 dx dy. 1294 \end{equation*} 1295 Ищем решение уравнения в виде обратного преобразования Фурье 1296 \(\phi(x,z)=\InverseFourierY{E(u,v)}{x,z}\). Подставляя[fn::Выражение 1297 \(v={-i}{u}\) не подходит в данной задаче, поскольку потенциал скорости должен 1298 стремиться к нулю с увеличением глубины до бесконечности.] \(v={i}{u}\) в 1299 формулу, решение перепишется как 1300 \begin{equation} 1301 \label{eq-guessed-sol-2d} 1302 \phi(x,z) = \InverseFourierY{e^{2\pi u z}E(u)}{x}. 1303 \end{equation} 1304 Для того чтобы подстановка \(z=\zeta(x,t)\) не помешала использованию 1305 преобразований Фурье в решении, перепишем\nbsp{}eqref:eq-guessed-sol-2d в виде 1306 свертки: 1307 \begin{equation*} 1308 \phi(x,z) 1309 = 1310 \Fun{z} 1311 \ast 1312 \InverseFourierY{E(u)}{x}, 1313 \end{equation*} 1314 где \(\Fun{z}\)\nbsp{}--- некоторая функция, вид которой будет определен в 1315 разделе\nbsp{}[[#sec-compute-delta]] и для которой выполняется соотношение 1316 \(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Подставляя выражение для \(\phi\) в 1317 граничное условие, получим 1318 \begin{equation*} 1319 \zeta_t 1320 = 1321 \left( i f(x) - 1/\SqrtZetaX \right) 1322 \left[ 1323 \Fun{z} 1324 \ast 1325 \InverseFourierY{2\pi u E(u)}{x} 1326 \right], 1327 \end{equation*} 1328 где \(f(x) = {\zeta_x}/{\sqrt{1 + \zeta_x^2}} - \zeta_x\). Применяя преобразование 1329 Фурье к обеим частям, получаем выражение для коэффициентов \(E\): 1330 \begin{equation*} 1331 E(u) = 1332 \frac{1}{2\pi u} 1333 \frac{ 1334 \FourierY{\zeta_t / \left(i f(x) - 1/\SqrtZetaX\right)}{u} 1335 }{ 1336 \FourierY{\Fun{z}}{u} 1337 } 1338 \end{equation*} 1339 Выполняя подстановку \(z=\zeta(x,t)\) и подставляя полученное выражение 1340 в\nbsp{}eqref:eq-guessed-sol-2d, получаем окончательное выражение для 1341 \(\phi(x,z)\): 1342 \begin{equation} 1343 \label{eq-solution-2d} 1344 \boxed{ 1345 \phi(x,z) 1346 = 1347 \InverseFourierY{ 1348 \frac{e^{2\pi u z}}{2\pi u} 1349 \frac{ 1350 \FourierY{ \zeta_t / \left(i f(x) - 1/\SqrtZetaX\right) }{u} 1351 }{ 1352 \FourierY{ \Fun{\zeta(x,t)} }{u} 1353 } 1354 }{x}. 1355 } 1356 \end{equation} 1357 1358 Множитель \(e^{2\pi u z}/(2\pi u)\) делает график функции, от которой берется 1359 обратное преобразования Фурье, несимметричным относительно оси \(OY\). Это 1360 затрудняет применение БПФ, поскольку оно требует периодичную функцию, 1361 принимающую нулевые значения на концах промежутка. В то же время, использование 1362 численного интегрирования вместо БПФ не позволит получить преимущество над 1363 решением всей системы уравнений с помощью разностных схем. Эту проблему можно 1364 обойти, используя формулу\nbsp{}eqref:eq-solution-2d-full для жидкости конечной 1365 глубины с заведомо большим значением глубины водоема \(h\). Вывод формулы дан в 1366 следующем разделе. 1367 1368 **** Формула для жидкости конечной глубины. 1369 На дне водоема вертикальная составляющая скорости перемещения жидкости должна 1370 равняться нулю, т.е.\nbsp{}\(\phi_z=0\) на \(z=-h\), где \(h\)\nbsp{}--- глубина 1371 водоема. В этом случае пренебречь равенством \(v = -i u\), полученным из 1372 уравнения Лапласа, нельзя, и решение ищется в виде 1373 \begin{equation} 1374 \phi(x,z) 1375 = 1376 \InverseFourierY{ 1377 \left( C_1 e^{2\pi u z} + C_2 e^{-2\pi u z} \right) 1378 E(u) 1379 }{x}. 1380 \label{eq-guessed-sol-2d-full} 1381 \end{equation} 1382 Подставляя \(\phi\) в условие на дне водоема, получим 1383 \begin{equation*} 1384 C_1 e^{-2\pi u h} - C_2 e^{2\pi u h} = 0, 1385 \end{equation*} 1386 откуда имеем \(C_1=\frac{1}{2}C{e}^{2\pi{u}{h}}\) и 1387 \(C_2=-\frac{1}{2}C{e}^{-2\pi{u}{h}}\). Константа \(C\) здесь произвольна, 1388 поскольку при подстановке станет частью неизвестных коэффициентов \(E(u)\). 1389 Подставляя полученные выражения для \(C_1\) и \(C_2\) 1390 в\nbsp{}eqref:eq-guessed-sol-2d-full, получаем выражение 1391 \begin{equation*} 1392 \phi(x,z) = \InverseFourierY{ \Sinh{2\pi u (z+h)} E(u) }{x}. 1393 \end{equation*} 1394 Подставляя \(\phi\) в граничное условие на свободной поверхности, получаем 1395 \begin{align*} 1396 \zeta_t & = f(x) \InverseFourierY{ 2\pi i u \Sinh{2\pi u (z+h)} E(u) }{x} \\ 1397 & - \frac{1}{\SqrtZetaX} 1398 \InverseFourierY{ 2\pi u \SinhX{2\pi u (z+h)} E(u) }{x}. 1399 \end{align*} 1400 Здесь \(\sinh\) и \(\cosh\) дают схожие результаты вблизи свободной поверхности, и, 1401 поскольку эта область является наиболее интересной с точки зрения практического 1402 применения, положим \(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\). Выполняя 1403 аналогичные предыдущему разделу операции, получаем окончательное выражение для 1404 \(\phi(x,z)\): 1405 \begin{equation} 1406 \boxed{ 1407 \phi(x,z,t) 1408 = 1409 \InverseFourierY{ 1410 \frac{\Sinh{2\pi u (z+h)}}{2\pi u} 1411 \frac{ 1412 \FourierY{ \zeta_t / \left(i f(x) - 1/\SqrtZetaX\right) }{u} 1413 }{ 1414 \FourierY{ \FunSecond{\zeta(x,t)} }{u} 1415 } 1416 }{x}, 1417 } 1418 \label{eq-solution-2d-full} 1419 \end{equation} 1420 где \(\FunSecond{z}\)\nbsp{}--- некоторая функция, вид которой будет определен 1421 в\nbsp{}[[#sec-compute-delta]] и для которой выполняется соотношение 1422 \(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\). 1423 1424 **** Сведение к формулам линейной теории волн. 1425 Справедливость полученных формул проверим, подставив в качестве \(\zeta(x,t)\) 1426 известные аналитические выражения для плоских волн. Символьные вычисления 1427 преобразований Фурье в этом разделе производились с помощью пакета 1428 Mathematica\nbsp{}cite:mathematica10. В линейной теории широко используется 1429 предположение о малости амплитуд волн, что позволяет упростить исходную систему 1430 уравнений\nbsp{}eqref:eq-problem-2d до 1431 \begin{align*} 1432 & \phi_{xx}+\phi_{zz}=0,\\ 1433 & \zeta_t = -\phi_z & \text{на }z=\zeta(x,t), 1434 \end{align*} 1435 решение которой запишется как 1436 \begin{equation*} 1437 \phi(x,z,t) 1438 = 1439 -\InverseFourierY{ 1440 \frac{e^{2\pi u z}}{2\pi u} 1441 \FourierY{\zeta_t}{u} 1442 }{x} 1443 . 1444 \end{equation*} 1445 Профиль прогрессивной волны задается формулой 1446 \(\zeta(x,t)=A\cos(2\pi(kx-t))\). Подстановка этого выражения 1447 в\nbsp{}eqref:eq-solution-2d дает равенство 1448 \(\phi(x,z,t)=-\frac{A}{k}\sin(2\pi(kx-t))\Sinh{2\pi{k}{z}}\). Чтобы свести его 1449 к формуле линейной теории волн, представим гиперболический косинус в 1450 экспоненциальной форме и отбросим член, содержащий \(e^{-2\pi{k}{z}}\), как 1451 противоречащий условию \(\phi\underset{z\rightarrow-\infty}{\longrightarrow}0\). 1452 После взятия действительной части выражения получится известная формула линейной 1453 теории \(\phi(x,z,t)=\frac{A}{k}e^{2\pi{k}{z}}\sin(2\pi(kx-t))\). Аналогично, 1454 предположение о малости амплитуд волн позволяет упростить 1455 формулу\nbsp{}eqref:eq-solution-2d-full до 1456 \begin{equation*} 1457 \phi(x,z,t) 1458 = 1459 -\InverseFourierY{ 1460 \frac{\Sinh{2\pi u (z+h)}}{2\pi u \Sinh{2\pi u h}} 1461 \FourierY{\zeta_t}{u} 1462 }{x}. 1463 \end{equation*} 1464 Подстановка формулы для прогрессивной плоской волны вместо \(\zeta(x,t)\) дает 1465 равенство 1466 \begin{equation} 1467 \label{eq-solution-2d-linear} 1468 \phi(x,z,t)=\frac{A}{k} 1469 \frac{\Sinh{2 \pi k (z+h)}}{ \Sinh{2 \pi k h} } 1470 \sin(2 \pi (k x-t)), 1471 \end{equation} 1472 что соответствует формуле линейной теории для конечной глубины. 1473 1474 Различные записи решения уравнения Лапласа, в которых затухающая экспонента 1475 может встречаться как со знаком "+", так и со знаком "-", могут стать причиной 1476 разницы между формулами линейно теории и формулами, выведенными в данной работе, 1477 где вместо \(\sinh\) используется \(\cosh\). Выражение 1478 \(\frac{\Sinh{2\pi{k}(z+h)}}{\Sinh{2\pi{k}{h}}}\approx\frac{\sinh(2\pi{k}(z+h))}{\sinh(2\pi{k}{h})}\) 1479 превращается в строгое равенство на поверхности, и разница между правой левой 1480 частью увеличивается при приближении к дну водоема (для достаточно большой 1481 глубины ошибка вблизи поверхности жидкости незначительна). Поэтому для 1482 достаточно большой глубины можно использовать любую из функций (\(\cosh\) или 1483 \(\sinh\)) для вычисления потенциала скорости вблизи взволнованной поверхности. 1484 1485 Сведение формул\nbsp{}eqref:eq-solution-2d и\nbsp{}eqref:eq-solution-2d-full к 1486 формулам линейной теории волн показывает, что формула\nbsp{}eqref:eq-solution-2d 1487 для жидкости бесконечной глубины не подходит для вычисления потенциала скорости 1488 с использованием метода Фурье, т.к.\nbsp{}не обладает необходимой для 1489 преобразования Фурье симметрией. Однако, для такого случая можно использовать 1490 формулу для конечной глубины, полагая \(h\) равным некоторому характерному 1491 значению глубины. Для стоячих волн сведение к формулам линейной теории 1492 происходит с аналогичными предположениями. 1493 1494 *** Трехмерное поле потенциала скорости 1495 В трех измерениях исходная система уравнений\nbsp{}eqref:eq-problem 1496 переписывается как 1497 \begin{align} 1498 \label{eq-problem-3d} 1499 & \phi_{xx} + \phi_{yy} + \phi_{zz} = 0,\\ 1500 & \zeta_t + \zeta_x\phi_x + \zeta_y\phi_y 1501 = \FracSqrtZetaY{\zeta_x} \phi_x 1502 + \FracSqrtZetaY{\zeta_y} \phi_y 1503 - \FracSqrtZetaY{1} \phi_z, \nonumber\\ 1504 & \text{на }z=\zeta(x,y,t).\nonumber 1505 \end{align} 1506 Для ее решения также воспользуемся методом Фурье. Возьмем преобразование Фурье 1507 от обоих частей уравнений Лапласа и получим 1508 \begin{equation*} 1509 -4 \pi^2 \left( u^2 + v^2 + w^2 \right) 1510 \FourierY{\phi(x,y,z)}{u,v,w} = 0, 1511 \end{equation*} 1512 откуда имеем \(w=\pm{i}\sqrt{u^2+v^2}\). Решение уравнения будем искать в виде 1513 обратного преобразования Фурье 1514 \(\phi(x,y,z)=\InverseFourierY{E(u,v,w)}{x,y,z}\). Подставляя 1515 \(w=i\sqrt{u^2+v^2}=i\Kveclen\) в исходную формулу, получаем 1516 \begin{equation*} 1517 \phi(x,y,z) = \InverseFourierY{ 1518 \left( 1519 C_1 e^{2\pi \Kveclen z} 1520 -C_2 e^{-2\pi \Kveclen z} 1521 \right) 1522 E(u,v) 1523 }{x,y}. 1524 \end{equation*} 1525 Подставляя \(\phi\) в условие на дне водоема аналогично двухмерному случаю, 1526 получаем 1527 \begin{equation} 1528 \label{eq-guessed-sol-3d} 1529 \phi(x,y,z) = \InverseFourierY{ 1530 \Sinh{2\pi \Kveclen (z+h)} E(u,v) 1531 }{x,y}. 1532 \end{equation} 1533 Подставляя выражение для \(\phi\) в граничное условие, получим 1534 \begin{equation*} 1535 \arraycolsep=1.4pt 1536 \begin{array}{rl} 1537 \zeta_t = & i f_1(x,y) \InverseFourierY{2 \pi u \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\ 1538 + & i f_2(x,y) \InverseFourierY{2 \pi v \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\ 1539 - & f_3(x,y) \InverseFourierY{2 \pi \sqrt{u^2+v^2} \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} 1540 \end{array} 1541 \end{equation*} 1542 где \(f_1(x,y)=\zeta_x/{\SqrtZetaY}-\zeta_x\), 1543 \(f_2(x,y)=\zeta_y/{\SqrtZetaY}-\zeta_y\) и \(f_3(x,y)=1/\SqrtZetaY\). 1544 1545 Также как и в разделе\nbsp{}[[#sec-pressure-2d]] мы предполагаем, что 1546 \(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\) вблизи свободной поверхности, 1547 однако в трехмерном случае этого недостаточно для решения задачи. Для того чтобы 1548 получить явную формулу для коэффициентов \(E\), мы должны предположить, что 1549 преобразования Фурье в равенстве имеют радиально симметричные ядра, 1550 т.е.\nbsp{}заменить \(u\) и \(v\) на \(\Kveclen\). Есть два момента, 1551 поддерживающих это предположение. Во-первых, в численной реализации 1552 интегрирование ведется по положительным волновым числам, так что знак \(u\) и 1553 \(v\) не влияет на решение. Во-вторых, скорость роста \(\cosh\) в ядре интеграла 1554 значительно выше, чем скорость роста \(u\) или \(\Kveclen\), так что замена 1555 слабо влияет на величину решения. Несмотря на эти два момента, использование 1556 более математически строго подхода было бы предпочтительнее. 1557 1558 Выполняя замену, применяя преобразование Фурье к обеим частям равенства, и, 1559 подставляя результат в\nbsp{}eqref:eq-guessed-sol-3d, получаем выражение для 1560 \(\phi\): 1561 \begin{equation*} 1562 \label{eq-phi-3d} 1563 \phi(x,y,z,t) = \InverseFourierY{ 1564 \frac{ \Sinh{\smash{2\pi \Kveclen (z+h)}} }{ 2\pi\Kveclen } 1565 \frac{ \FourierY{ \zeta_t / \left( i f_1(x,y) + i f_2(x,y) - f_3(x,y) \right)}{u,v} } 1566 { \FourierY{\mathcal{D}_3\left( x,y,\zeta\left(x,y\right) \right)}{u,v} } 1567 }{x,y}, 1568 \end{equation*} 1569 где 1570 \(\FourierY{\mathcal{D}_3\left(x,y,z\right)}{u,v}=\Sinh{\smash{2\pi\Kveclen{}z}}\). 1571 1572 *** Формулы нормировки для потенциалов скоростей 1573 :PROPERTIES: 1574 :CUSTOM_ID: sec-compute-delta 1575 :END: 1576 1577 В решениях\nbsp{}eqref:eq-solution-2d и\nbsp{}eqref:eq-solution-2d-full 1578 двухмерной задачи определения поля давлений присутствуют функции 1579 \(\Fun{z}=\InverseFourierY{e^{2\pi{u}{z}}}{x}\) и 1580 \(\FunSecond{z}=\InverseFourierY{\Sinh{2\pi{u}{z}}}{x}\), которые могут быть 1581 записаны несколькими различными аналитическими выражениями и представляют 1582 сложность при вычислении на компьютере. Каждая функция\nbsp{}--- это 1583 преобразование Фурье от линейной комбинации экспонент, которое сводится к плохо 1584 определенной дельта функции комплексного аргумента 1585 (см.\nbsp{}табл.\nbsp{}[[tab-delta-functions]]). Обычно такого типа функции 1586 записывают как произведение дельта функций от действительной и мнимой части, 1587 однако, в данном случае такой подход не работает, поскольку взятие обратного 1588 преобразования Фурье не даст экспоненту, что сильно исказит результирующее поле 1589 скоростей. Для получения однозначного аналитического выражения можно 1590 воспользоваться нормировкой \(1/\Sinh{2\pi{u}{h}}\) (которая также включается в 1591 выражение для коэффициентов \(E(u)\)). Численные эксперименты показывают, что 1592 нормировка хоть и позволяет получить адекватное поле скоростей, оно мало 1593 отличается от выражений из линейной теории волн, в которых члены с \(\zeta\) 1594 опускаются. Как следствие, формула для трехмерного случая не выводилась. 1595 1596 #+name: tab-delta-functions 1597 #+caption[Формулы для вычисления \(\Fun{z}\) и \(\FunSecond{z}\)]: 1598 #+caption: Формулы для вычисления \(\Fun{z}\) и \(\FunSecond{z}\) из 1599 #+caption: разд.\nbsp{}[[#sec-pressure-2d]], использующие нормировку для 1600 #+caption: исключения неоднозначности определения дельта функции комплексного 1601 #+caption: аргумента. 1602 #+attr_latex: :booktabs t 1603 | Функция | Без нормировки | С нормировкой | 1604 |-------------------+--------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------| 1605 | \(\Fun{z}\) | \(\delta (x+i z)\) | \(\frac{1}{2 h}\mathrm{sech}\left(\frac{\pi (x-i (h+z))}{2 h}\right)\) | 1606 | \(\FunSecond{z}\) | \(\frac{1}{2}\left[\delta (x-i z) + \delta (x+i z) \right]\) | \(\frac{1}{4 h}\left[\text{sech}\left(\frac{\pi (x-i (h+z))}{2 h}\right)+\text{sech}\left(\frac{\pi (x+i(h+z))}{2 h}\right)\right]\) | 1607 1608 ** Верификация полей потенциалов скоростей 1609 :PROPERTIES: 1610 :CUSTOM_ID: sec-compare-formulae 1611 :END: 1612 1613 Сравнение полученных общих формул\nbsp{}eqref:eq-solution-2d 1614 и\nbsp{}eqref:eq-solution-2d-full с известными формулами линейной теории волн 1615 позволяет оценить различие между полями скоростей для волн как больших, так и 1616 малых амплитуд. В общем случае аналитическое выражение для потенциала скорости 1617 неизвестно даже для плоских волн, поэтому сравнение производится численно. Имея 1618 ввиду выводы раздела\nbsp{}[[#sec-pressure-2d]], сравниваются только формулы для 1619 конечной глубины. 1620 1621 **** Отличие от формул линейной теории волн. 1622 Для того чтобы получить поля потенциалов скоростей, плоские волны различных 1623 амплитуд были сгенерированы. Волновые числа в преобразованиях Фурье выбирались 1624 на интервале от \(0\) до максимального волнового числа, определяемого численно 1625 из полученной взволнованной поверхности. Эксперименты проводились для волн малых 1626 и больших амплитуд. 1627 1628 Эксперимент показал, что поля потенциалов скоростей для волн больших амплитуд, 1629 полученные по формуле\nbsp{}eqref:eq-solution-2d-full для конечной глубины и по 1630 формуле\nbsp{}eqref:eq-solution-2d-linear линейной теории, качественно 1631 отличаются (см.\nbsp{}рис.\nbsp{}[[fig-potential-field-nonlinear]]). Во-первых, 1632 контуры потенциала скорости имеют вид затухающей синусоиды, что отличается от 1633 овальной формы, описываемой линейной теории волн. Во-вторых, по мере приближения 1634 к дну водоема потенциал скорости затухает гораздо быстрее, чем в линейной 1635 теории, а область, где сконцентрирована большая часть энергии волны, еще больше 1636 приближена к ее гребню. Аналогичный численный эксперимент, в котором из 1637 формулы\nbsp{}eqref:eq-solution-2d-full были исключены члены, которыми 1638 пренебрегают в рамках линейной теории волн, показал полное соответствие 1639 получившихся полей потенциалов скоростей (насколько это позволяет сделать 1640 машинная точность). 1641 1642 #+name: fig-potential-field-nonlinear 1643 #+header: :width 8 :height 11 1644 #+begin_src R :file build/plain-wave-velocity-field-comparison-ru.pdf 1645 source(file.path("R", "velocity-potentials.R")) 1646 par(pty="s", family="serif") 1647 nlevels <- 41 1648 levels <- pretty(c(-200,200), nlevels) 1649 palette <- colorRampPalette(c("blue", "lightyellow", "red")) 1650 col <- palette(nlevels-1) 1651 1652 # linear solver 1653 par(fig=c(0,0.95,0,0.5),new=TRUE) 1654 arma.plot_velocity_potential_field( 1655 file.path("build", "arma-benchmarks", "verification-orig", "plain_wave_linear_solver"), 1656 levels=levels, 1657 col=col 1658 ) 1659 1660 # high-amplitude solver 1661 par(fig=c(0,0.95,0.5,1),new=TRUE) 1662 arma.plot_velocity_potential_field( 1663 file.path("build", "arma-benchmarks", "verification-orig", "plain_wave_high_amplitude_solver"), 1664 levels=levels, 1665 col=col 1666 ) 1667 1668 # legend 1 1669 par(pty="m",fig=c(0.80,1,0.5,1), new=TRUE) 1670 arma.plot_velocity_potential_field_legend( 1671 levels=levels, 1672 col=col 1673 ) 1674 1675 # legend 2 1676 par(pty="m",fig=c(0.80,1,0,0.5), new=TRUE) 1677 arma.plot_velocity_potential_field_legend( 1678 levels=levels, 1679 col=col 1680 ) 1681 #+end_src 1682 1683 #+name: fig-potential-field-nonlinear 1684 #+caption[Поля потенциала скорости прогрессивной волны]: 1685 #+caption: Поля потенциала скорости прогрессивной волны 1686 #+caption: \(\zeta(x,y,t) = \cos(2\pi x - t/2)\). Поле, полученное по общей 1687 #+caption: формуле (сверху) и по формуле из линейной теории волн (снизу). 1688 #+attr_latex: :width \textwidth 1689 #+RESULTS: fig-potential-field-nonlinear 1690 [[file:build/plain-wave-velocity-field-comparison-ru.pdf]] 1691 1692 **** Отличие от формул теории волн малой амплитуды. 1693 Эксперимент, в котором сравнивались поля потенциалов скоростей, полученные 1694 численно различными формулами, показал, что поля скоростей, полученные по 1695 формуле\nbsp{}eqref:eq-solution-2d-full и формуле для волн малой 1696 амплитуды\nbsp{}eqref:eq-old-sol-2d, сопоставимы для волн малых амплитуд. В этом 1697 эксперименте использовались две реализации взволнованной морской поверхности, 1698 полученные по модели АР: одна содержала волны малой амплитуды, другая\nbsp{}--- 1699 большой. Интегрирование в формуле\nbsp{}eqref:eq-solution-2d-full велось по 1700 диапазону волновых чисел, полученному из морской поверхности. Для волн малой 1701 амплитуды обе формулы показали сопоставимые результаты (разница в значениях 1702 скорости приписывается стохастической природе модели АР), в то время как для 1703 волн больших амплитуд устойчивое поле скоростей дала только 1704 формула\nbsp{}eqref:eq-solution-2d-full (рис.\nbsp{}[[fig-velocity-field-2d]]). 1705 Таким образом, общая формула\nbsp{}eqref:eq-solution-2d-full показывает 1706 удовлетворительные результаты, не вводя ограничения на амплитуду волн. 1707 1708 #+name: fig-velocity-field-2d 1709 #+header: :width 8 :height 11 1710 #+begin_src R :file build/large-and-small-amplitude-velocity-field-comparison-ru.pdf 1711 source(file.path("R", "velocity-potentials.R")) 1712 linetypes = c("solid", "dashed") 1713 par(mfrow=c(2, 1), family="serif") 1714 arma.plot_velocity( 1715 file.path("data", "velocity", "low-amp"), 1716 file.path("data", "velocity", "low-amp-0"), 1717 linetypes=linetypes, 1718 ylim=c(-2,2) 1719 ) 1720 arma.plot_velocity( 1721 file.path("data", "velocity", "high-amp"), 1722 file.path("data", "velocity", "high-amp-0"), 1723 linetypes=linetypes, 1724 ylim=c(-2,2) 1725 ) 1726 #+end_src 1727 1728 #+name: fig-velocity-field-2d 1729 #+caption[Поля потенциала скорости волн малых и больших амплитуд]: 1730 #+caption: Сравнение полей потенциала скорости, полученных 1731 #+caption: по общей формуле (\(u_1\)) и формуле для волн малой амплитуды 1732 #+caption: (\(u_2\)): поле скоростей для поверхности волн малой (сверху) 1733 #+caption: и большой амплитуды (снизу). 1734 #+attr_latex: :width \textwidth 1735 #+RESULTS: fig-velocity-field-2d 1736 [[file:build/large-and-small-amplitude-velocity-field-comparison-ru.pdf]] 1737 1738 ** Выводы 1739 Полученные в данном разделе формулы позволяют произвести численный расчет поля 1740 скоростей (а значит и поля давлений) вблизи дискретно или математически заданной 1741 взволнованной морской поверхности, минуя предположения линейной теории и теории 1742 волн малых амплитуд. Для волн малых амплитуд новые формулы дают то же самое поле 1743 потенциала скорости, что и формулы из линейной теории волн. Для волн больших 1744 амплитуд использование новых формул приводит к сдвигу области, где 1745 сконцентрирована больше всего энергии волны, ближе к ее гребню. 1746 1747 * Высокопроизводительный программный комплекс для моделирования морского волнения 1748 :PROPERTIES: 1749 :CUSTOM_ID: sec-hpc 1750 :END: 1751 1752 ** Реализация для систем с общей памятью (SMP) 1753 *** Генерация взволнованной поверхности 1754 **** Параллельные алгоритмы для моделей АР, СС и ЛХ. 1755 :PROPERTIES: 1756 :CUSTOM_ID: sec-parallel 1757 :END: 1758 1759 Несмотря на то что модели АР и СС являются частью одной смешанной модели, они 1760 имеют разные параллельные алгоритмы, которые отличаются от тривиального 1761 алгоритма модели ЛХ. Алгоритм АР заключается в разбиении взволнованной 1762 поверхности на части одинакового размера вдоль каждой из координатных осей и их 1763 параллельном вычислении с учетом каузальных ограничений, накладываемых 1764 авторегрессионными зависимостями между точками поверхности. В модели СС такие 1765 зависимости отсутствуют, а ее формула представляет собой свертку белого шума с 1766 коэффициентами модели, которая сводится к вычислению трех преобразований Фурье 1767 посредством теоремы о свертке. Таким образом, алгоритм СС заключается в 1768 параллельном вычислении свертки, которая основана на вычислении БПФ. Наконец, 1769 алгоритм ЛХ делается параллельным простым вычислением каждой точки параллельно в 1770 нескольких потоках. Таким образом, параллельная реализация модели АРСС включает 1771 в себя два параллельных алгоритма, по одному для каждой составляющей модели, 1772 которые сложнее, чем алгоритм для модели ЛХ. 1773 1774 Основная особенность формулы модели АР\nbsp{}--- авторегрессионные зависимости 1775 между точками взволнованной поверхности по каждому из измерений,\nbsp{}--- 1776 который препятствуют параллельному вычислению каждой точки поверхности. Вместо 1777 этого поверхность разделяется на равные части по каждому из измерений, и для 1778 каждой части устанавливаются информационные зависимости, определяющие порядок 1779 вычисления. На рис.\nbsp{}[[fig-ar-cubes]] показаны эти зависимости. Стрелка 1780 обозначает зависимость одной части от наличия другой, т.е.\nbsp{}вычисление 1781 части может начаться, только если все части, от которых она зависит, уже 1782 вычислены. Здесь часть \(A\) не имеет зависимостей, части \(B\) и \(D\) зависят 1783 только от \(A\), а часть \(E\) зависит от \(A\), \(B\) и \(C\). В общем случае, 1784 каждая часть зависит от всех частей, имеющих предыдущий индекс хотя бы по одному 1785 из измерений (если такие части существуют). Первая часть не имеет никаких 1786 зависимостей; размер каждой части по каждому из измерений выбирается большим или 1787 равным количеству коэффициентов по соответствующему измерению, чтобы принимать 1788 во внимание только прилегающие части в разрешении зависимостей. 1789 1790 #+name: fig-ar-cubes 1791 #+header: :width 5 :height 5 1792 #+begin_src R :file build/ar-cubes-ru.pdf 1793 source(file.path("R", "common.R")) 1794 par(family="serif") 1795 arma.plot_ar_cubes_2d(3, 3, xlabel="Индекс части (X)", ylabel="Индекс части (Y)") 1796 #+end_src 1797 1798 #+name: fig-ar-cubes 1799 #+caption[Авторегрессионные зависимости между частями поверхности]: 1800 #+caption: Авторегрессионные зависимости между частями взволнованной поверхности. 1801 #+RESULTS: fig-ar-cubes 1802 [[file:build/ar-cubes-ru.pdf]] 1803 1804 Каждая часть имеет трехмерный индекс и состояние завершения. Алгоритм начинается 1805 с отправки всех объектов, содержащих эту информацию, в очередь. После этого 1806 параллельные потоки запускаются, каждый поток последовательно ищет первую часть, 1807 для которой все зависимости удовлетворены (путем проверки состояния каждой из 1808 частей), извлекает эту часть из очереди, генерирует взволнованную поверхность 1809 для этой части и устанавливает состояние завершения. Алгоритм заканчивается, 1810 когда очередь становится пустой. Доступ к очереди из разных потоков 1811 синхронизируется посредством блокировок. Алгоритм подходит для SMP машин; в 1812 случае MPP части, от которых зависит данная, должны быть предварительно 1813 скопированы на узел, на котором будут проводится вычисления. 1814 1815 Таким образом, параллельный алгоритм модели АР сводится к созданию 1816 минималистичного планировщика задач, в котором 1817 - каждая задача соответствует части взволнованной поверхности, 1818 - порядок выполнения задач определяется авторегрессионными зависимостями, и 1819 - очередь обрабатывается простым пулом потоков, в котором каждый поток в цикле 1820 извлекает из очереди первую задачу, для которой все зависимые задачи 1821 завершены, и выполняет ее. 1822 1823 В модели СС, в отличие от АР, отсутствуют авторегрессионные зависимости между 1824 точками; вместо этого, каждая точка поверхности зависит от предыдущих по времени 1825 и пространству значений белого шума. Формула модели СС может быть переписана как 1826 свертка белого шума с коэффициентами модели в качестве ядра. Используя теорему о 1827 свертке, свертка переписывается как обратное преобразование Фурье от 1828 произведения прямых преобразования Фурье от белого шума и коэффициентов. 1829 Поскольку количество коэффициентов СС много меньше, чем количество точек 1830 поверхности, то параллельное БПФ не подходит, поскольку требует дополнение 1831 массива коэффициентов нулями для того чтобы его размер совпадал с размером 1832 массива точек поверхности. Вместо этого, поверхность разбивается на части по 1833 каждому из измерений, которые дополняются нулями, чтобы получить размер равный 1834 количеству коэффициентов домноженному на два. Затем, преобразование Фурье 1835 вычисляется параллельно для каждой части, домножается на заранее вычисленное 1836 преобразование Фурье от коэффициентов и обратное преобразование Фурье 1837 вычисляется от результата. После этого, каждая часть записывается в выходной 1838 массив, а перекрывающие друг друга (из-за заполнения нулями) точки складываются 1839 друг с другом. Этот алгоритм известен в области обработки сигналов как 1840 "overlap-add"\nbsp{}cite:oppenheim1989discrete,svoboda2011efficient,pavel2013algorithms. 1841 Заполнение нулями необходимо для предотвращения маскированных ошибок: без него 1842 результатом вычислений была бы циклическая свертка. 1843 1844 Несмотря на то что алгоритм модели СС разбивает поверхность на те же самые части 1845 (но, возможно, другого размера), что и алгоритм модели АР, отсутствие 1846 авторегрессионных зависимостей между ними позволяет вычислять их параллельно без 1847 использования специализированного планировщика задач. Однако, этот алгоритм 1848 также требует дополнение каждой части нулями, чтобы результат вычислений 1849 совпадал с результатом, полученным по исходной формуле модели СС. Таким образом, 1850 алгоритм модели СС лучше масштабируется на большое количество узлов ввиду 1851 отсутствия информационных зависимостей между частями, но размер частей больше, 1852 чем в алгоритме модели АР. 1853 1854 Отличительной особенностью алгоритма модели ЛХ является его простота: чтобы 1855 сделать его параллельным, взволнованная поверхность разделяется на части равного 1856 размера, каждая из которых генерируется параллельно. Между частями отсутствуют 1857 информационные зависимости, что делает этот алгоритм подходящим для вычисления 1858 на видеокарте: каждый аппаратный поток вычисляет свою точку поверхности. При 1859 этом, функции синуса и косинуса в формуле модели медленно вычисляются на 1860 процессоре, что делает видеокарту еще более выгодным выбором. 1861 1862 Итого, несмотря на то что модели АР и СС являются составными частями одной 1863 смешанной модели, их параллельные алгоритмы фундаментально отличаются и являются 1864 более сложными, чем тривиальный параллельный алгоритм модели ЛХ. Эффективная 1865 реализация алгоритма АР требует специализированного планировщика задач для учета 1866 авторегрессионных зависимостей между частями взволнованной поверхности, в то 1867 время как алгоритм СС требует дополнения каждой части нулями, чтобы иметь 1868 возможность обработать их параллельно. В отличие от этих моделей, в модели ЛХ 1869 отсутствуют информационные зависимости между частями, но она также требует 1870 больше вычислительных ресурсов (операций с плавающей точкой в секунду) для 1871 трансцендентных функций в формуле. 1872 1873 **** Параллельный алгоритм генерации белого шума 1874 Чтобы исключить периодичность из генерируемой моделью морского волнения 1875 реализации взволнованной поверхности, для генерации белого шума необходимо 1876 использовать ГПСЧ с достаточно большим периодом. В качестве такого генератора в 1877 работе используется параллельная реализация вихря 1878 Мерсенна\nbsp{}cite:matsumoto1998mersenne с периодом \(2^{19937}-1\). Это 1879 позволяет создавать апериодические реализации взволнованной морской поверхности 1880 для любых сценариев применения, встречаемых на практике. 1881 1882 Запуск нескольких ГПСЧ с разными начальными состояниями в параллельных потоках 1883 не гарантирует некоррелированность генерируемых последовательностей 1884 псевдослучайных чисел, и для предоставления такой гарантии используется алгоритм 1885 динамического создания вихрей Мерсенна\nbsp{}cite:matsumoto1998dynamic. Суть 1886 алгоритма заключается в поиске таких матриц начальных состояний генераторов, 1887 которые бы дали максимально некоррелированные последовательности псевдослучайных 1888 чисел при параллельном запуске нескольких вихрей Мерсенна с этими начальными 1889 состояниями. Поскольку на поиск начальных состояний тратится значительное 1890 количество процессорного времени, то вектор состояний создается предварительно 1891 для заведомо большего количества параллельных потоков и сохраняется в файл, 1892 который впоследствии считывается основной программой перед началом генерации 1893 белого шума. 1894 1895 **** Устранение интервала разгона. 1896 В модели АР значение подъема взволнованной поверхности в каждой точке зависит от 1897 предыдущих по пространству и времени значений, из-за чего в начале реализации 1898 образуется так называемый /интервал разгона/ 1899 (см.\nbsp{}рис.\nbsp{}[[fig-ramp-up-interval]])\nbsp{}--- промежуток, на котором 1900 реализация не соответствует заданной АКФ. Способ решения этой проблемы зависит 1901 от контекста, в котором производится имитационное моделирование. Если реализация 1902 используется в контексте расчета остойчивости судна без учета маневрирования, то 1903 интервал никак не повлияет результаты эксперимента, поскольку находится на 1904 границе (далеко от исследуемого морского объекта). Альтернативным подходом 1905 является генерация взволнованной поверхности на интервале разгона моделью ЛХ и 1906 генерация остальной реализации с помощью модели АР. Если изучается остойчивость 1907 судна в условиях маневрирования, то интервал проще всего исключить из реализации 1908 (размер интервала примерно равен числу коэффициентов АР). Однако, это приводит к 1909 потере большого числа точек, поскольку исключение происходит по каждому из трех 1910 измерений. 1911 1912 #+name: fig-ramp-up-interval 1913 #+begin_src R :file build/ramp-up-interval-ru.pdf 1914 source(file.path("R", "common.R")) 1915 par(family="serif") 1916 arma.plot_ramp_up_interval(label="Интервал разгона") 1917 #+end_src 1918 1919 #+caption[Интервал разгона в начале реализации процесса АР]: 1920 #+caption: Интервал разгона в начале реализации процесса АР. 1921 #+name: fig-ramp-up-interval 1922 #+RESULTS: fig-ramp-up-interval 1923 [[file:build/ramp-up-interval-ru.pdf]] 1924 1925 **** Производительность реализаций на OpenMP и OpenCL. 1926 :PROPERTIES: 1927 :header-args:R: :results output raw :exports results 1928 :END: 1929 1930 Различия в параллельных алгоритмах моделей делает их эффективными на разных 1931 архитектурах процессоров, и для того чтобы найти наиболее эффективную из них, 1932 все модели были протестированы как на процессоре, так и на видеокарте. 1933 1934 Модели АР и СС не требуют высоко оптимизированных кодов, для того чтобы быть 1935 эффективными, их производительность высокая и без использования сопроцессоров; 1936 на это есть две причины. Во-первых, эти модели не используют трансцендентные 1937 функции (синусы, косинусы и экспоненты) в отличие от модели ЛХ. Все вычисления 1938 (исключая коэффициенты модели) производятся через полиномы, которые эффективно 1939 вычисляются на современных процессорах, используя инструкции FMA. Во-вторых, 1940 вычисление давлений происходит по явной формуле с использованием нескольких 1941 вложенных БПФ. Поскольку двухмерное БПФ одного и того же размера постоянно 1942 вычисляется на каждом временном срезе, его коэффициенты (комплексные экспоненты) 1943 вычисляются один раз для всех временных срезов, а дальнейшие вычисления включают 1944 в себя лишь несколько трансцендентных функций. В случае модели СС, 1945 производительность также повышается за счет вычисления свертки с помощью БПФ. 1946 Таким образом, высокая производительность моделей АР и СС обусловлена скудным 1947 использованием трансцендентных функций и интенсивным использованием БПФ, не 1948 говоря уже о том, что высокая скорость сходимости и отсутствие периодичности 1949 позволяют использовать гораздо меньше коэффициентов по сравнению с моделью ЛХ. 1950 1951 #+name: tab-gpulab 1952 #+caption[Конфигурация системы "Gpulab"]: 1953 #+caption: Конфигурация системы "Gpulab". 1954 #+attr_latex: :booktabs t 1955 | Процессор | AMD FX-8370 | 1956 | Память процессора | 16ГБ | 1957 | Видеокарта | GeForce GTX 1060 | 1958 | Память видеокарты | 6ГБ | 1959 | Жесткий диск | WDC WD40EZRZ, 5400об./мин. | 1960 | Количество процессорных ядер | 4 | 1961 | Количество потоков на ядро | 2 | 1962 1963 Программная реализация использует несколько библиотек математических функций, 1964 численных алгоритмов и примитивов визуализации (перечисленных в 1965 таб.\nbsp{}[[tab-arma-libs]]), и написана с использованием нескольких технологий 1966 параллельного программирования (OpenMP, OpenCL) для возможности выбора наиболее 1967 эффективной. Для каждой технологии и каждой модели была оптимизирована генерация 1968 взволнованной поверхности (за исключением модели СС, для которой была сделана 1969 только реализация OpenMP). Вычисление потенциала скорости было реализовано на 1970 OpenMP, реализация на OpenCL была сделана только для визуализации взволнованной 1971 поверхности в реальном времени. Для каждой технологии программа 1972 перекомпилировалась, запускалась несколько раз и измерялась производительность 1973 каждой высокоуровневой функции посредством системных часов. 1974 1975 Результаты тестирования представлены в таб.\nbsp{}[[tab-arma-performance]]. Все 1976 тесты запускались на машине с видеокартой, характеристики которой собраны в 1977 таб.\nbsp{}[[tab-gpulab]]. Все тесты запускались с одними и теми же входными 1978 параметрами для всех моделей: реализация длиной в 10000 сек. и размер выходной 1979 сетки \(40\times40\)м. Единственный параметр, который отличался, это порядок 1980 (количество коэффициентов): порядок моделей АР и СС был выбран равным \(7,7,7\), 1981 а порядок модели ЛХ \(40,40\). Это связано с б\oacute{}льшим количеством 1982 коэффициентов, необходимым для исключения периодичности модели ЛХ. 1983 1984 Во всех тестах генерация взволнованной поверхности и НБП занимают 1985 б\oacute{}льшую часть времени работы, а вычисление потенциала скорости вместе с 1986 остальными подпрограммами лишь малую его часть. 1987 1988 #+name: tab-arma-libs 1989 #+caption[Список библиотек, используемых в программной реализации]: 1990 #+caption: Список библиотек, используемых в программной реализации. 1991 #+attr_latex: :booktabs t 1992 | Библиотека | Назначение | 1993 |--------------------------------------------------------------+----------------------------------| 1994 | DCMT\nbsp{}cite:matsumoto1998dynamic | параллельный ГПСЧ | 1995 | Blitz\nbsp{}cite:veldhuizen1997will,veldhuizen2000techniques | многомерные массивы | 1996 | GSL\nbsp{}cite:gsl2008scientific | вычисление ФПР, ФР, БПФ | 1997 | | проверка стационарности процесса | 1998 | LAPACK, GotoBLAS\nbsp{}cite:goto2008high,goto2008anatomy | определение коэффициентов АР | 1999 | GL, GLUT\nbsp{}cite:kilgard1996opengl | трехмерная визуализация | 2000 | CGAL\nbsp{}cite:fabri2009cgal | триангуляция волновых чисел | 2001 2002 Модель АР показывает наибольшую производительность в реализации на OpenMP и 2003 наименьшую в реализации на OpenCL, что также является наибольшей и наименьшей 2004 производительностью среди всех комбинаций моделей и технологий. В самой 2005 оптимальной комбинации производительность АР в 4,5 раз выше, чем 2006 производительность СС, и в 20 раз выше, чем производительность ЛХ; в самой 2007 неоптимальной комбинации\nbsp{}--- в 137 раз медленней, чем СС и в два раза 2008 медленней, чем ЛХ. Отношение между наибольшей (OpenMP) и наименьшей (OpenCL) 2009 производительностью модели АР составляет несколько сотен. Это объясняется тем, 2010 что формула модели\nbsp{}eqref:eq-ar-process эффективно отображается на 2011 архитектуру центрального процессора, который отличается от видеокарты наличием 2012 нескольких кэшей, памятью с низкой пропускной способностью и небольшим 2013 количеством модулей для операций с плавающей точкой. 2014 - Эта формула не содержит трансцендентных функций (синусов, косинусов и 2015 экспонент), 2016 - все операции умножения и сложения в формуле реализуются посредством FMA 2017 инструкций процессора, и 2018 - эффективное использование (локальность) кэша достигается путем использования 2019 библиотеки Blitz, которая реализует оптимизированный обход элементов 2020 многомерного массива, основанный на заполняющей пространство кривой Гильберта. 2021 В отличие от центрального процессора, видеокарта имеет меньшее количество кэшей, 2022 память с высокой пропускной способностью и большое количество модулей для 2023 операций с плавающей точкой, что является наименее благоприятным сценарием для 2024 модели АР. 2025 - Формула модели АР не содержит трансцендентных функций, которые могли бы 2026 компенсировать высокие задержки памяти, 2027 - в видеокарте присутствуют инструкции FMA, но они не увеличивают 2028 производительность из-за высоких задержек памяти, и 2029 - оптимальный обход многомерного массива не использовался ввиду отсутствия 2030 библиотек, реализующих его для видеокарты. 2031 Наконец, архитектура видеокарты не содержит примитивы синхронизации, позволяющих 2032 эффективно реализовать авторегрессионные зависимости между отдельными частями 2033 взволнованной поверхности; вместо этого отдельная подпрограмма OpenCL 2034 запускается для каждой части, а управление информационными зависимостями между 2035 ними осуществляется на стороне центрального процессора. Таким образом, в случае 2036 модели АР архитектура центрального процессора превосходит архитектуру 2037 видеокарты, поскольку более эффективно обрабатывает сложные информационные 2038 зависимости, простые вычисления (сложения и умножения) и сложные шаблоны доступа 2039 к памяти. 2040 2041 #+header: :results output raw :exports results 2042 #+name: tab-arma-performance 2043 #+begin_src R :results output org :exports results 2044 source(file.path("R", "benchmarks.R")) 2045 options(arma.mark=",") 2046 model_names <- list( 2047 ar.x="АР", 2048 ma.x="СС", 2049 lh.x="ЛХ", 2050 ar.y="АР", 2051 ma.y="СС", 2052 lh.y="ЛХ", 2053 Row.names="\\orgcmidrule{2-4}{5-6}Подпрограмма" 2054 ) 2055 row_names <- list( 2056 determine_coefficients="Определение коэффициентов", 2057 validate="Проверка модели", 2058 generate_surface="Генерация поверхности", 2059 nit="НБП", 2060 write_all="Запись вывода в файл", 2061 copy_to_host="Копирование данных с GPU", 2062 velocity="Выч. потенциалов скорости" 2063 ) 2064 arma.print_openmp_vs_opencl(model_names, row_names) 2065 #+end_src 2066 2067 #+name: tab-arma-performance 2068 #+caption[Производительность реализаций OpenMP и OpenCL]: 2069 #+caption: Время работы (с.) реализаций OpenMP и OpenCL моделей АР, СС и ЛХ. 2070 #+attr_latex: :booktabs t 2071 #+RESULTS: tab-arma-performance 2072 | | | | OpenMP | | OpenCL | 2073 | | <r> | <r> | <r> | <r> | <r> | 2074 | \orgcmidrule{2-4}{5-6}Подпрограмма | АР | СС | ЛХ | АР | ЛХ | 2075 |------------------------------------+------+------+--------+--------+--------| 2076 | Определение коэффициентов | 0,02 | 0,01 | 0,19 | 0,01 | 1,19 | 2077 | Проверка модели | 0,08 | 0,10 | | 0,08 | | 2078 | Генерация поверхности | 1,26 | 5,57 | 350,98 | 769,38 | 0,02 | 2079 | НБП | 7,11 | 7,43 | | 0,02 | | 2080 | Копирование данных с GPU | | | | 5,22 | 25,06 | 2081 | Выч. потенциалов скорости | 0,05 | 0,05 | 0,06 | 0,03 | 0,03 | 2082 | Запись вывода в файл | 0,27 | 0,27 | 0,27 | 0,28 | 0,27 | 2083 2084 В отличие от модели АР, модель ЛХ показывает наибольшую производительность на 2085 видеокарте и наименьшую производительность на центральном процессоре. Причиной 2086 этому служат 2087 - большое количество трансцендентных функций в ее формуле, которые компенсируют 2088 высокие задержки памяти, 2089 - линейный шаблон доступа к памяти, который позволяет векторизовать вычисления и 2090 объединить операции доступа к памяти из разных потоков, и 2091 - отсутствие информационных зависимостей между точками взволнованной 2092 поверхности. 2093 Несмотря на то что видеокарта на тестовой системе более производительна, чем 2094 центральный процессор (с точки зрения операций с плавающей точкой в секунду), 2095 общая производительность модели ЛХ по сравнению с моделью АР меньше. Причиной 2096 этому служит медленная передача данных между памятью видеокарты и центрального 2097 процессора. 2098 2099 Модель СС быстрее, чем модель ЛХ, но медленнее, чем модель АР. Поскольку свертка 2100 в ее формуле реализована с помощью БПФ, ее производительность зависит от 2101 производительности реализации БПФ: библиотека GSL в случае центрального 2102 процессора и clFFT в случае видеокарты. В данной работе производительность 2103 модели СС на видеокарте не была протестирована ввиду недоступности трехмерного 2104 БПФ в библиотеке clFFT; если бы преобразование было доступно, оно могло бы 2105 сделать модель даже более быстрой, чем АР. 2106 2107 НБП занимает меньше времени на видеокарте, чем на центральном процессоре, 2108 однако, если принять во внимание время передачи данных между ними, время 2109 становится сопоставимым. Это объясняется большим количеством трансцендентных 2110 функций, которые необходимо вычислить для каждой точки взволнованной поверхности 2111 для преобразования ее координат \(z\). Для каждой точки нелинейное 2112 трансцендентное уравнение\nbsp{}eqref:eq-distribution-transformation решается 2113 методом бисекции. Видеокарта выполняет эту задачу в несколько сотен раз быстрее, 2114 чем центральный процессор, но тратит много времени на передачу результата 2115 обратно в память процессора. Таким образом, единственная возможность 2116 оптимизировать эту подпрограмму заключается в использовании метода поиска корня 2117 уравнения с квадратичной сходимостью, чтобы уменьшить количество трансцендентных 2118 функций, которые необходимо вычислить. 2119 2120 **** Производительность ввода-вывода. 2121 :PROPERTIES: 2122 :header-args:R: :results output raw :exports results 2123 :CUSTOM_ID: sec-io-perf 2124 :END: 2125 2126 Несмотря на то что в тестах из предыдущего раздела запись данных в файлы не 2127 занимает большое количество времени, использование монтируемых по сети файловых 2128 систем может замедлить этот этап программы. Для его оптимизации части 2129 взволнованной поверхности записываются в файл, как только генерация всего 2130 временного среза завершена (рис.\nbsp{}[[fig-arma-io-events]]): запускается 2131 отдельный поток, который начинает запись, как только первый временной срез 2132 становится доступен, и завершает, после того как основная группа потоков 2133 заканчивает вычисления и последний срез записывается в файл. Суммарное время, 2134 затрачиваемое на ввод/вывод, увеличивается, но суммарное время работы программы 2135 уменьшается, поскольку ввод/вывод производится параллельно вычислениям 2136 (таб.\nbsp{}[[tab-arma-io-performance]]). Использование этого подхода для локальных 2137 систем имеет тот же эффект, но выигрыш в производительности меньше. 2138 2139 #+header: :results output raw :exports results 2140 #+name: tab-arma-io-performance 2141 #+begin_src R 2142 source(file.path("R", "benchmarks.R")) 2143 options(arma.mark=",") 2144 suffix_names <- list( 2145 xfs.x="XFS", 2146 xfs.y="XFS", 2147 nfs.x="NFS", 2148 nfs.y="NFS", 2149 gfs.x="GlusterFS", 2150 gfs.y="GlusterFS", 2151 Row.names="\\orgcmidrule{2-4}{5-7}Подпрограмма" 2152 ) 2153 top_names <- c("I", "II") 2154 row_names <- list( 2155 generate_surface="Генерация поверхности", 2156 write_all="Запись вывода в файл" 2157 ) 2158 arma.print_sync_vs_async_io(suffix_names, row_names, top_names) 2159 #+end_src 2160 2161 #+name: tab-arma-io-performance 2162 #+caption[Производительность ввода/вывода для XFS, NFS и GlusterFS]: 2163 #+caption: Время работы подпрограмм (сек.) при использовании файловых систем 2164 #+caption: XFS, NFS и GlusterFS совместно с последовательным (I) и 2165 #+caption: параллельным (II) выводом в файл. 2166 #+attr_latex: :booktabs t 2167 #+RESULTS: tab-arma-io-performance 2168 | | | | I | | | II | 2169 | | <r> | <r> | <r> | <r> | <r> | <r> | 2170 | \orgcmidrule{2-4}{5-7}Подпрограмма | XFS | NFS | GlusterFS | XFS | NFS | GlusterFS | 2171 |------------------------------------+------+------+-----------+------+------+-----------| 2172 | Генерация поверхности | 1,26 | 1,26 | 1,33 | 1,33 | 3,30 | 11,06 | 2173 | Запись вывода в файл | 0,28 | 2,34 | 10,95 | 0,00 | 0,00 | 0,00 | 2174 2175 #+name: fig-arma-io-events 2176 #+header: :height 9 :results output graphics 2177 #+begin_src R :file build/arma-io-events-ru.pdf 2178 source(file.path("R", "benchmarks.R")) 2179 par(family="serif") 2180 names <- list(xlab="Время, сек.", ylab="Поток") 2181 par(mfrow=c(3, 1), mar=c(4,4,4,0.5), family='serif', cex=1) 2182 arma.plot_io_events(names) 2183 #+end_src 2184 2185 #+name: fig-arma-io-events 2186 #+caption[График событий для XFS, NFS и GlusterFS]: 2187 #+caption: График событий для XFS, NFS и GlusterFS, показывающий 2188 #+caption: временные интервалы для каждого из потоков, на протяжении которых 2189 #+caption: производится запись в файл (красный) и вычисления (черный). 2190 #+RESULTS: fig-arma-io-events 2191 [[file:build/arma-io-events-ru.pdf]] 2192 2193 *** Вычисление поля потенциала скорости 2194 **** Параллельное вычисление поля потенциала скорости. 2195 Тесты для моделей АР, СС, и ЛХ показали, что вычисление поля потенциала скорости 2196 занимает лишь малую долю суммарного времени работы программы, однако, абсолютное 2197 время вычисления на плотной сетке \(XY\) может быть большим. Одно из приложений, 2198 в котором используется плотная сетка,\nbsp{}--- это имитационное моделирование и 2199 визуализация в реальном времени взволнованной поверхности. Визуализация в 2200 реальном времени позволяет 2201 - настроить параметры модели и АКФ, мгновенно получая результат изменений, и 2202 - сравнить размер и форму областей, в которых сконцентрирована основная часть 2203 энергии волн. 2204 2205 Поскольку визуализация производится на видеокарте, вычисление потенциала 2206 скорости на центральном процессоре может сделать передачу данных между памятью 2207 двух устройств узким местом. Для того чтобы обойти это, программа использует 2208 программный интерфейс взаимодействия OpenCL/OpenGL, позволяющий создавать 2209 буферы, используемые совместно контекстами OpenCL и OpenGL, что исключает 2210 копирование данных между устройствами. Кроме того, формула трехмерного поля 2211 потенциала скорости\nbsp{}eqref:eq-phi-3d особенно подходит для вычисления на 2212 видеокарте: 2213 - она содержит трансцендентные функции (гиперболические косинусы и комплексные 2214 экспоненты); 2215 - она вычисляется на большой четырехмерной области \((t,x,y,z)\); 2216 - она явная и не имеет информационных зависимостей между отдельными точками в 2217 измерениях \(t\) и \(z\). 2218 2219 Для того чтобы выяснить, на сколько использование видеокарты может ускорить 2220 вычисления поля потенциала скорости, была протестирована упрощенная 2221 версия\nbsp{}eqref:eq-phi-3d: 2222 \begin{align} 2223 \label{eq:phi-linear} 2224 \phi(x,y,z,t) &= \InverseFourierY{ 2225 \frac{ \Sinh{\smash{2\pi \Kveclen (z+h)}} } 2226 { 2\pi\Kveclen \Sinh{\smash{2\pi \Kveclen h}} } 2227 \FourierY{ \zeta_t }{u,v} 2228 }{x,y}\nonumber \\ 2229 &= \InverseFourierY{ g_1(u,v) \FourierY{ g_2(x,y) }{u,v} }{x,y}. 2230 \end{align} 2231 Код, вычисляющий потенциал скорости, был переписан на языке OpenCL и его 2232 производительность сравнивалась с реализацией на OpenMP. 2233 2234 Для каждой реализации измерялось время работы соответствующих подпрограмм и 2235 время передачи данных между устройствами. Поле потенциала скорости вычислялось 2236 для одной точки по оси \(t\), 128 точек по оси \(z\), расположенных под 2237 взволнованной поверхностью, и для каждой точки по оси \(x\) и \(y\) 2238 четырехмерной сетки \((t,x,y,z)\). Между запусками программы изменялся размер 2239 сетки по оси \(x\). 2240 2241 В реализациях используются разные библиотеки БПФ: GNU Scientific Library 2242 (GSL)\nbsp{}cite:galassi2015gnu для OpenMP и clFFT\nbsp{}cite:clfft для OpenCL. 2243 Подпрограммы БПФ из этих библиотек имеют следующие особенности. 2244 - Порядок частот в БПФ у обоих библиотек разный. В случае clFFT элементы 2245 результирующего массива дополнительно сдвигаются, чтобы соответствовать 2246 корректному полю потенциала скорости. В случае GSL никакого сдвига не 2247 требуется. 2248 - Разрыв в точках сетки \((x,y)=(0,0)\) обрабатывается библиотекой clFFT 2249 автоматически, в то время как библиотека GSL выдает искаженные значения в этой 2250 точке, поэтому с случае GSL значения в этих точках интерполируются. 2251 Другие отличия подпрограмм БПФ, влияющие на производительность, не были 2252 выявлены. 2253 2254 **** Производительность кода, вычисляющего поле потенциала скорости. 2255 :PROPERTIES: 2256 :header-args:R: :results output raw :exports results 2257 :END: 2258 2259 Эксперименты показали, что реализация OpenCL превосходит OpenMP по 2260 производительности в 2--6 раз (рис.\nbsp{}[[fig-arma-realtime-graph]]), однако, 2261 распределение времени работы между подпрограммами отличается 2262 (таб.\nbsp{}[[tab-arma-realtime]]). В случае процессора больше всего времени 2263 тратится на вычисление \(g_1\), а в случае видеокарты время вычисления \(g_1\) 2264 сопоставимо с \(g_2\). Копирование результирующего поля потенциала скорости 2265 между процессором и видеокартой занимает \(\approx{}70\%\) общего времени 2266 вычисления этого поля. Вычисление \(g_2\) занимает больше всего времени в случае 2267 OpenCL и меньше всего времени в случае OpenMP. В обоих реализациях \(g_2\) 2268 вычисляется на центральном процессоре, поскольку готовая подпрограмма вычисления 2269 производной на OpenCL не была найдена. В случае OpenCL результат вычислений 2270 дублируется для каждой точки сетки по оси \(z\), для того чтобы перемножить все 2271 точки одного временного среза в одной подпрограмме OpenCL, а, затем, копируется 2272 в память видеокарты, что негативно сказывается на производительности. Все тесты 2273 запускались на машине с видеокартой, характеристики которой просуммированы в 2274 таб.\nbsp{}[[tab-storm]]. 2275 2276 #+name: tab-storm 2277 #+caption[Конфигурация системы "Storm"]: 2278 #+caption: Конфигурация системы "Storm". 2279 #+attr_latex: :booktabs t 2280 | Процессор | Intel Core 2 Quad Q9550 | 2281 | Память процессора | 8ГБ | 2282 | Видеокарта | AMD Radeon R7 360 | 2283 | Память видеокарты | 2ГБ | 2284 | Жесткий диск | Seagate Barracuda, 7200 об./мин. | 2285 | Количество процессорных ядер | 4 | 2286 2287 #+name: fig-arma-realtime-graph 2288 #+header: :results output graphics 2289 #+begin_src R :file build/realtime-performance-ru.pdf 2290 source(file.path("R", "benchmarks.R")) 2291 par(family="serif") 2292 data <- arma.load_realtime_data() 2293 arma.plot_realtime_data(data) 2294 title(xlab="Размер взволнованной поверхности по оси OX", ylab="Время, сек.") 2295 #+end_src 2296 2297 #+name: fig-arma-realtime-graph 2298 #+caption[Производительность кода для потенциала скорости]: 2299 #+caption: Сравнение производительности версий кода, вычисляющего поле 2300 #+caption: потенциала скорости, для центрального процессора (OpenMP) и 2301 #+caption: видеокарты (OpenCL). 2302 #+RESULTS: fig-arma-realtime-graph 2303 [[file:build/realtime-performance-ru.pdf]] 2304 2305 Причина разного распределения времени работы между подпрограммами OpenCL и 2306 OpenMP та же, что и в случае разной производительности модели АР на центральном 2307 процессоре и видеокарте: видеокарта выполняет больше операций с плавающей точкой 2308 в секунду и имеет больше модулей трансцендентных функций, чем процессор, что 2309 ускоряет вычисление \(g_1\), но в ней отсутствует кэш, который необходим для 2310 оптимизации нерегулярного шаблона доступа к памяти при вычислении \(g_2\). В 2311 отличие от модели АР, производительность вычисления многомерной производной на 2312 видеокарте легче увеличить ввиду отсутствия информационных зависимостей между 2313 точками: в данной работе оптимизация не была проведена ввиду отсутствия готовой 2314 реализации. Кроме того, такая реализация может позволить эффективно вычислить 2315 неупрощенную формулу полностью на видеокарте, поскольку опущенные в формуле 2316 функции также содержат производные. 2317 2318 #+name: tab-arma-realtime 2319 #+begin_src R 2320 source(file.path("R", "benchmarks.R")) 2321 options(arma.mark=",") 2322 routine_names <- list( 2323 harts_g1="Функция \\(g_1\\)", 2324 harts_g2="Функция \\(g_2\\)", 2325 harts_fft="БПФ", 2326 harts_copy_to_host="Копирование данных с видекарты" 2327 ) 2328 column_names <- c("Подпрограмма", "OpenMP", "OpenCL") 2329 data <- arma.load_realtime_data() 2330 arma.print_table_for_realtime_data(data, routine_names, column_names) 2331 #+end_src 2332 2333 #+name: tab-arma-realtime 2334 #+caption[Производительность вычисления поля потенциала скорости]: 2335 #+caption: Время работы (сек.) подпрограмм вычисления поля потенциала скорости 2336 #+caption: в реальном времени для взволнованной поверхности (размер по 2337 #+caption: \(OX\) равен 16384). 2338 #+attr_latex: :booktabs t 2339 #+RESULTS: tab-arma-realtime 2340 | | <r> | <r> | 2341 | Подпрограмма | OpenMP | OpenCL | 2342 |--------------------------------+--------+--------| 2343 | Функция \(g_1\) | 4,6730 | 0,0038 | 2344 | Функция \(g_2\) | 0,0002 | 0,8253 | 2345 | БПФ | 2,8560 | 0,3585 | 2346 | Копирование данных с видекарты | | 2,6357 | 2347 2348 Как и ожидалось, совместное использование одного буфера контекстами OpenCL и 2349 OpenGL увеличивает производительность путем исключения копирования данных между 2350 памятью центрального процессора и видеокарты, но также требует, чтобы данные 2351 были в формате вершин, с которым непосредственно работает OpenGL. Преобразование 2352 в этот формат выполняется быстро, однако результирующий массив занимает больше 2353 памяти, поскольку каждая точка записывается как вектор из трех компонент. Другим 2354 недостатком совместного использования OpenCL и OpenGL является требование ручной 2355 блокировки общего буфера: невыполнение этого требования может стать причиной 2356 появления артефактов изображения на экране, которые можно убрать, только 2357 перезагрузив компьютер. 2358 2359 *** Выводы 2360 Тесты показали, что видеокарта превосходит центральный процессор по 2361 производительности в задачах с большим количеством арифметических операций, 2362 требующих большое количество операций с плавающей точкой в секунду, однако, 2363 производительность падает, когда объем данных, которые нужно скопировать между 2364 памятью процессора и видеокарты, увеличивается или, когда шаблон доступа к 2365 памяти отличается от линейного. Первая проблема может быть решена путем 2366 использования сопроцессора, в котором память с высокой пропускной способностью 2367 расположена на одной плате вместе с процессором и основной памятью. Это 2368 исключает узкое место при пересылке данных, но может также увеличить время 2369 выполнения ввиду меньшего количества модулей для операций с плавающей точкой. 2370 Вторую проблему можно решить программно, при наличии библиотеки для вычисления 2371 многомерных производных на OpenCL. 2372 2373 Модели АР и СС превосходят модель ЛХ в тестах производительности и для этого не 2374 требуется наличие видеокарты. Их сильными сторонами с вычислительной точки 2375 зрения являются 2376 - отсутствие трансцендентных функций, и 2377 - простые алгоритмы, производительность которых зависит от производительности 2378 библиотеки для многомерных массивов и библиотеки для БПФ. 2379 Обеспечение основного функционала посредством низкоуровневых библиотек делает 2380 производительность программы переносимой: поддержка новых архитектур процессоров 2381 может быть добавлена путем замены библиотек. Наконец, использование явной 2382 формулы позволяет тратить лишь небольшую долю суммарного времени работы 2383 программы на вычисление поля потенциала скорости. Если бы такой формулы не было 2384 или она не содержала бы интегралы в виде преобразований Фурье, на вычисление 2385 поля потенциала скорости требовалось бы гораздо больше времени. 2386 2387 ** Отказоустойчивый планировщик пакетных задач 2388 *** Архитектура системы 2389 **** Физический уровень. 2390 Состоит из узлов и прямых/маршрутизируемых физических сетевых соединений. На 2391 этом уровне предполагается полная сетевая связность, т.е.\nbsp{}возможность 2392 отправлять сетевые пакеты между любой парой узлов кластера. 2393 2394 **** Уровень резидентных процессов. 2395 :PROPERTIES: 2396 :CUSTOM_ID: sec-daemon-layer 2397 :END: 2398 2399 Состоит из процессов, запущенных на узлах кластера и иерархических логических 2400 связей (главный/подчиненный) между ними. На каждом узле запущен ровно один 2401 процесс, поэтому эти понятия взаимозаменяемы в данной работе. Роль главного и 2402 подчиненного назначается резидентным процессам динамически, т.е.\nbsp{}любой 2403 физический узел может стать главным или подчиненным, или быть и тем, и другим 2404 одновременно. Для динамического назначения ролей используется алгоритм выбора 2405 лидера, не требующий периодической отправки широковещательных сообщений всем 2406 узлам кластера, вместо этого роль определяется IP-адресом узла. Подробное 2407 описание алгоритма приведено в разделе\nbsp{}[[#sec-node-discovery]]. Его сильными 2408 сторонами являются масштабируемость на большое количество узлов и низкие 2409 накладные расходы, что делает его пригодным для высокопроизводительных 2410 вычислений; его слабой стороной является искусственная зависимость места, 2411 занимаемого узлом в иерархии, от его IP-адреса, что делает его непригодным для 2412 виртуальных сред, где IP-адрес узла может динамически меняться. 2413 2414 Единственное предназначение древовидной иерархии состоит в балансировке нагрузки 2415 между узлами кластера. Нагрузка распределяется от текущего узла к его соседям в 2416 иерархии путем простой итерации по ним. Когда в иерархии появляются новые связи 2417 или меняются старые (из-за того что новый узел присоединяется к кластеру или 2418 работающий узел выходит из строя), резидентные процессы сообщают друг другу 2419 количество узлов, находящихся /за/ соответствующей связью в иерархии. Эта 2420 информация используется для равномерного распределения нагрузки, даже если 2421 распределенная программа запускается на подчиненном узле. Кроме того, иерархия 2422 уменьшает количество одновременных сетевых соединений между узлами: на каждую 2423 связь в иерархии устанавливается и поддерживается только одно 2424 соединение,\nbsp{}--- что уменьшает вероятность возникновения перегрузки сети 2425 при большом количестве узлов в кластере. 2426 2427 Балансировка нагрузки реализована следующим образом. Когда узел \(A\) пытается 2428 стать подчиненным узлу \(B\), он отправляет сообщение на соответствующий 2429 IP-адрес, передавая, сколько узлов уже связаны с ним в иерархии (включая себя 2430 самого). После того как все связи в иерархии установлены, каждый узел имеет 2431 достаточно информации, чтобы сказать, сколько узлов находится за каждой из его 2432 связей, а суммарное количество узлов за всеми связями равно общему количеству 2433 узлов в кластере. Если это связь между главным и подчиненным узлом, и 2434 подчиненный узел хочет узнать, сколько узлов находится за связью с главным, то 2435 из суммарного количества узлов за связями главного со всеми его подчиненным 2436 узлами он вычитает количество узлов за главным, чтобы получить правильную сумму. 2437 Для распределения управляющих объектов (см.\nbsp{}разд.\nbsp{}[[#sec-kernel-layer]]) 2438 между узлами используется циклический алгоритм (round-robin), 2439 т.е.\nbsp{}производится итерация по всем связям текущего узла, включая связь с 2440 главным узлом, и, принимая во внимание количество узлов, находящихся за каждой 2441 из связей. Таким образом, даже если программа запускается на подчиненном узле, 2442 находящимся внизу иерархии, ее управляющие объекты распределяются равномерно 2443 между всеми узлами кластера. 2444 2445 Предложенный подход может быть расширен путем включения сложной логики в 2446 алгоритм распределения нагрузки. Вместе с количеством узлов, находящихся за 2447 связью в иерархии, может быть отправлена любая метрика, которая необходима для 2448 реализации этой логики. Однако, если эта метрика обновляется чаще, чем 2449 изменяется иерархия, или периодически, то и сообщения о текущем состоянии будут 2450 отправляться чаще. Для того чтобы эти пересылки не повлияли на 2451 производительность программ, для них можно использовать отдельную физическую 2452 сеть. Также, когда происходит изменение иерархии, управляющие объекты, которые 2453 уже отправлены на узлы до этого изменения, не учитываются при распределении 2454 нагрузки, из-за чего частые изменения иерархии могут стать причиной 2455 неравномерной загрузке узлов кластера (которая, однако, станет равномерной со 2456 временем). 2457 2458 Динамическое назначение ролей главного и подчиненного узла вместе с 2459 распределением управляющих объектов делает архитектуру системы однородной в 2460 рамках одного кластера. На каждом узле запущен один и тот же резидентный 2461 процесс, и никакой предварительной конфигурации не требуется, чтобы объединить 2462 процессы, запущенные на разных узлах, в иерархию. 2463 2464 **** Уровень управляющих объектов. 2465 :PROPERTIES: 2466 :CUSTOM_ID: sec-kernel-layer 2467 :END: 2468 2469 Ключевая особенность, отсутствующая в текущих технологиях параллельного 2470 программирования,\nbsp{}--- это возможность задавать иерархические зависимости 2471 между параллельно выполняющимися задачами. Когда такая зависимость есть, 2472 руководящая задача становится ответственной за перезапуск подчиненной, 2473 завершившейся неудачно, на оставшихся узлах кластера. Для того чтобы иметь 2474 возможность перезапустить задачу, у которой нету главенствующей над ней задачи, 2475 создается ее копия, и отправляется на другой узел 2476 (см.\nbsp{}разд.\nbsp{}[[#sec-fail-over]]). Существует несколько систем, которые 2477 способны выполнять направленные ациклические графы задач 2478 параллельно\nbsp{}cite:acun2014charmpp,islam2012oozie, однако, графы не подходят 2479 для определения связей руководитель-подчиненный между задачами, поскольку 2480 вершина графа может иметь несколько входящих ребер (а значит несколько 2481 руководящих узлов). 2482 2483 Основное назначение предлагаемого подхода состоит в упрощении разработки 2484 распределенных приложений, работающих в пакетном режиме, и промежуточного 2485 программного обеспечения. Идея состоит в обеспечении отказоустойчивости на самом 2486 низком из возможных уровней. Реализация разделена на два слоя: нижний слой 2487 состоит из подпрограмм и классов для приложений, работающих на одном узле (без 2488 взаимодействия по сети), а верхний слой предназначен для приложений, работающих 2489 на произвольном количестве узлов. Имеется два типа сильно связанных 2490 сущностей\nbsp{}--- это /управляющие объекты/ (или /объекты/ для краткости) и 2491 /конвейеры/,\nbsp{}--- которые составляют основу программы. 2492 2493 Управляющие объекты реализуют логику (порядок выполнения) программы в методах 2494 ~act~ и ~react~ и хранят состояние текущей ветки исполнения. Как логика, так и 2495 состояние задаются программистом. В методе ~act~ какая-либо функция либо 2496 вычисляется непосредственно, либо разлагается на вызовы вложенных функций 2497 (представляемые подчиненными управляющими объектами), которые впоследствии 2498 отправляются на конвейер. В методе ~react~ подчиненные управляющие объекты, 2499 вернувшиеся с конвейера, обрабатываются их родительским объектом. Вызовы методов 2500 ~act~ и ~react~ производятся асинхронно внутри потоков, присоединенных к 2501 конвейеру. Для каждого управляющего объекта метод ~act~ вызывается только один 2502 раз, и для нескольких объектов вызовы происходят параллельно друг другу, в то 2503 время как метод ~react~ вызывается один раз для каждого подчиненного объекта, и 2504 все вызовы происходят в одном потоке для предотвращения одновременного изменения 2505 состояния несколькими потоками (для разных родительских объектов могут 2506 использоваться разные потоки). 2507 2508 Конвейеры осуществляют асинхронные вызовы методов ~act~ и ~react~, стараясь 2509 сделать как можно больше вызовов одновременно, учитывая предоставляемый 2510 платформой параллелизм (количество процессорных ядер на узле и количество узлов 2511 в кластере). Конвейер включает в себя пул управляющих объектов, содержащий все 2512 подчиненные объекты, отправленные в него родителями, и пул потоков, 2513 обрабатывающий эти объекты в соответствии с правилами, описанными в предыдущем 2514 параграфе. Для каждого устройства используется отдельный конвейер. Существуют 2515 конвейеры для параллельной обработки, обработки по расписанию (периодические и 2516 отложенные задачи) и промежуточный конвейер для обработки управляющих объектов 2517 на других узлах кластера (см.\nbsp{}рис.\nbsp{}[[fig-subord-ppl]]). 2518 2519 По принципу работы механизм управляющих объектов и конвейеров напоминает 2520 механизм работы процедур и стеков вызовов, с тем лишь преимуществом, что методы 2521 объектов вызываются асинхронно и параллельно друг другу (насколько это позволяет 2522 логика программы). Поля управляющего объекта\nbsp{}--- это локальные переменные 2523 стека, метод ~act~\nbsp{}--- это последовательность процессорных инструкций 2524 перед вложенным вызовом процедуры, а метод ~react~\nbsp{}--- это 2525 последовательность инструкций после вложенного вызова. Создание и отправка на 2526 конвейер подчиненного объекта\nbsp{}--- это вложенный вызов процедуры. Наличие 2527 двух методов обусловливается асинхронностью вложенных вызовов и помогает 2528 заменить активное ожидание завершения подчиненных объектов пассивным при помощи 2529 конвейеров. Конвейеры, в свою очередь, позволяют реализовать пассивное ожидание 2530 и вызывают правильные методы, анализируя внутреннее состояние объектов. 2531 2532 #+name: fig-subord-ppl 2533 #+begin_src dot :exports results :file build/subord-ppl-ru.pdf 2534 graph G { 2535 2536 node [fontname="Old Standard",fontsize=14,margin="0.055,0",shape=box] 2537 graph [nodesep="0.25",ranksep="0.25",rankdir="LR" margin=0] 2538 edge [arrowsize=0.66] 2539 2540 subgraph cluster_daemon { 2541 label="Родительский процесс" 2542 style=filled 2543 color=lightgrey 2544 graph [margin=8] 2545 2546 factory [label="Фабрика"] 2547 parallel_ppl [label="Параллельный\nконвейер"] 2548 io_ppl [label="Конвейер\nввода/вывода"] 2549 sched_ppl [label="Конвейер\nдля таймера"] 2550 net_ppl [label="Конвейер для\nсетевых устройств"] 2551 proc_ppl [label="Конвейер\nдля процессов"] 2552 2553 upstream [label="Пул потоков\nupstream"] 2554 downstream [label="Пул потоков\ndownstream"] 2555 } 2556 2557 factory--parallel_ppl 2558 factory--io_ppl 2559 factory--sched_ppl 2560 factory--net_ppl 2561 factory--proc_ppl 2562 2563 subgraph cluster_hardware { 2564 label="Вычислительные устройства" 2565 style=filled 2566 color=lightgrey 2567 graph [margin=8] 2568 2569 cpu [label="CPU"] 2570 core0 [label="Ядро 0"] 2571 core1 [label="Ядро 1"] 2572 core2 [label="Ядро 2"] 2573 core3 [label="Ядро 3"] 2574 2575 storage [label="Устройства\nхранения"] 2576 disk0 [label="Диск 0"] 2577 2578 network [label="Сетевые\nкарты"] 2579 nic0 [label="СК 0"] 2580 2581 timer [label="Таймер"] 2582 2583 } 2584 2585 core0--cpu 2586 core1--cpu 2587 core2--cpu 2588 core3--cpu 2589 2590 disk0--storage 2591 nic0--network 2592 2593 parallel_ppl--upstream 2594 parallel_ppl--downstream 2595 2596 upstream--{core0,core1,core2,core3} [style="dashed"] 2597 downstream--core0 [style="dashed"] 2598 2599 io_ppl--core0 [style="dashed"] 2600 io_ppl--disk0 [style="dashed"] 2601 sched_ppl--core0 [style="dashed"] 2602 sched_ppl--timer [style="dashed"] 2603 net_ppl--core0 [style="dashed"] 2604 net_ppl--nic0 [style="dashed"] 2605 proc_ppl--core0 [style="dashed"] 2606 2607 subgraph cluster_children { 2608 style=filled 2609 color=white 2610 2611 subgraph cluster_child0 { 2612 label="Дочерний процесс 0" 2613 style=filled 2614 color=lightgrey 2615 2616 app0_factory [label="Фабрика"] 2617 app0 [label="Конвейер\nдочернего\nпроцесса"] 2618 } 2619 2620 # subgraph cluster_child1 { 2621 # label="Дочерний процесс 1" 2622 # style=filled 2623 # color=lightgrey 2624 # 2625 # app1_factory [label="Фабрика"] 2626 # app1 [label="Конвейер\nдочернего процесса"] 2627 # } 2628 } 2629 2630 proc_ppl--app0 2631 # proc_ppl--app1 2632 2633 app0_factory--app0 [constraint=false] 2634 # app1_factory--app1 [constraint=false] 2635 2636 } 2637 #+end_src 2638 2639 #+caption[Отображение конвейеров на вычислительные устройства]: 2640 #+caption: Отображение конвейеров родительского и дочернего процессов 2641 #+caption: на вычислительные устройства. Сплошные линии обозначают агрегацию, 2642 #+caption: пунктирные линии\nbsp{}--- отображение между логическими и 2643 #+caption: физическими сущностями. 2644 #+attr_latex: :width \textwidth 2645 #+label: fig-subord-ppl 2646 #+RESULTS: fig-subord-ppl 2647 [[file:build/subord-ppl-ru.pdf]] 2648 2649 **** Программная реализация. 2650 Из соображений эффективности конвейер объектов и методы обеспечения 2651 отказоустойчивости (которые будут описаны далее) были реализованы во фреймворке 2652 на языке C++: с точки зрения автора язык C слишком низкоуровневый для написания 2653 распределенных программ, а использование языка Java влечет за собой накладные 2654 расходы, и не популярно в высокопроизводительных вычислениях. Фреймворк 2655 называется Bscheduler и находится на этапе проверки концепции. 2656 2657 **** Программный интерфейс. 2658 Каждый управляющий объект имеет четыре типа полей (перечисленных в 2659 табл.\nbsp{}[[tab-kernel-fields]]): 2660 - поля, относящиеся к управлению потоком выполнения, 2661 - поля, определяющие исходное местонахождение управляющего объекта, 2662 - поля, определяющие текущее местонахождение управляющего объекта и 2663 - поля, определяющие целевое местонахождение управляющего объекта. 2664 2665 #+name: tab-kernel-fields 2666 #+caption[Описание полей управляющего объекта]: 2667 #+caption: Описание полей управляющего объекта. 2668 #+attr_latex: :booktabs t :align lp{0.7\textwidth} 2669 | Поле | Описание | 2670 |-------------------+----------------------------------------------------------------------------------------------------------------| 2671 | ~process_id~ | Идентификатор процесса (приложения) в рамках кластера, которому принадлежит управляющий объект. | 2672 | ~id~ | Идентификатор управляющего объекта внутри процесса. | 2673 | ~pipeline~ | Идентификатор конвейера, который обрабатывает управляющий объект. | 2674 | | | 2675 | ~exit_code~ | Результат выполнения управляющего объекта. | 2676 | ~flags~ | Вспомогательный битовые флаги, используемые при планировании. | 2677 | ~time_point~ | Момент времени, в который запланировано выполнение управляющего объекта. | 2678 | | | 2679 | ~parent~ | Адрес/идентификатор родительского объекта. | 2680 | ~src_ip~ | IP-адрес исходного узла кластера. | 2681 | ~from_process_id~ | Идентификатор процесса в рамках кластера, из которого пришел управляющий объект. | 2682 | | | 2683 | ~principal~ | Адрес/идентификатор целевого управляющего объекта (объекта, к которому текущий направляется или возвращается). | 2684 | ~dst_ip~ | IP-адрес целевого узла кластера. | 2685 2686 При создании каждому управляющему объекту присваивается адрес родительского 2687 объекта и идентификатор конвейера. Если другие поля не установлены, то 2688 управляющий объект является /upstream/-объектом\nbsp{}--- объектом, который 2689 можно распределить на любой узел кластера и любое ядро процессора для извлечения 2690 параллелизма. Если адрес целевого управляющего объекта установлен, то объект 2691 является /downstream/-объектом\nbsp{}--- объектом, который можно направить 2692 только к целевому, а ядро процессора, на которое будет направлен этот объект, 2693 определяется адресом/идентификатором целевого объекта. Если downstream-объект 2694 предполагается направить на другой узел, то IP-адрес целевого узла должен быть 2695 установлен, иначе система не сможет найти целевой объект. 2696 2697 Когда выполнение управляющего объекта завершается (после завершения метода 2698 ~act~), этот объект явно направляется к какому-либо другому объекту, эта 2699 директива явно вызывается в методе ~act~. Обычно, по окончании выполнения 2700 объекта, он направляется к своему родительскому объекту путем установки поля 2701 ~principal~ в значение адреса/идентификатора родителя, целевого IP-адреса в 2702 значение исходного IP-адреса и идентификатора процесса в значение идентификатора 2703 исходного процесса. После этого управляющий объект становится 2704 downstream-объектом, направляется системой на узел, где находится его текущий 2705 руководитель, без использования алгоритма балансировки нагрузки. Для 2706 downstream-объектов метод ~react~ их родителя вызывается конвейером с текущим 2707 объектом в качестве аргумента, чтобы позволить родителю собрать результат 2708 выполнения объекта. 2709 2710 Не существует способа предоставить мелкозернистую отказоустойчивость к сбоям 2711 узлов кластера, если в программе присутствуют downstream-объекты кроме тех, что 2712 возвращаются к своим родителям. Вместо этого, проверяется код завершения 2713 управляющего объекта и выполняется пользовательская подпрограмма для 2714 восстановления. Если проверка на ошибку отсутствует, система перезапускает 2715 выполнение, начиная с первого родительского объекта, который не создавал 2716 downstream-объектов. Это означает, что если решаемая программой задача имеет 2717 информационные зависимости между частями вычисляемыми параллельно, и узел 2718 выходит из строя во время вычисления этих частей, то эти вычисления 2719 перезапускается с самого начала, отбрасывая части, вычисленные ранее. Такого не 2720 происходит в высоко параллельных программах, где параллельные части не имеют 2721 таких информационных зависимостей между друг другом: в этом случае только части 2722 с вышедших из строя узлов вычисляются заново, а все ранее вычисленные части 2723 сохраняются. 2724 2725 В отличие от функции ~main~ программ, основанных на библиотеке MPI, первый 2726 (главный) управляющий объект изначально запускается только на одном узле, а 2727 другие узлы кластера используются при необходимости. Это позволяет использовать 2728 больше узлов для участков кода с большой степенью параллелизма, и меньше для 2729 других участков. Похожий подход применяется во фреймворках для обработки больших 2730 объемов данных\nbsp{}cite:dean2008mapreduce,vavilapalli2013yarn --- узлы, на 2731 которых будет выполняться задача выбираются в зависимости от того, где физически 2732 находятся входные файлы. 2733 2734 **** Отображение имитационных моделей на архитектуру системы. 2735 Модели АР и СС реализованы в программном комплексе, работающем по принципу 2736 вычислительного конвейера, в котором каждое звено применяет некоторую функцию к 2737 выходным данным предыдущего звена. Звенья конвейера распределяются по узлам 2738 вычислительного кластера, чтобы сделать возможным параллелизм по операциям, а 2739 затем данные, перемещающиеся между звеньями конвейера распределяются между 2740 ядрами процессора, чтобы сделать возможным параллелизм по данным. На 2741 рис.\nbsp{}[[fig-pipeline]] представлена схема такого конвейера. Здесь 2742 прямоугольниками со скругленными углами обозначены звенья конвейера, обычными 2743 прямоугольниками\nbsp{}--- массивы объектов из предметной области задачи, 2744 передаваемые от одного звена к другому, а стрелками\nbsp{}--- направление 2745 передачи данных. Некоторые звенья разделены на /секции/, каждая из которых 2746 обрабатывает отдельную часть массива. Если звенья соединены без использования 2747 /барьера/ (горизонтальная или вертикальная полоса), то передача отдельных 2748 объектов между такими звеньями происходит параллельно с вычислениями, по мере их 2749 готовности. Секции работают параллельно на нескольких ядрах процессора 2750 (нескольких узлах кластера). Таким образом, между множеством ядер процессора, 2751 секций звеньев конвейера и объектами устанавливается сюръективное отображение, 2752 т.е.\nbsp{}на одном ядре процессора может работать несколько секций звеньев 2753 конвейера, каждая из которых может обрабатывать несколько объектов 2754 последовательно, но одна секция не может работать сразу на нескольких ядрах, а 2755 объект не может обрабатываться сразу несколькими секциями конвейера. Таким 2756 образом, программа представляет собой конвейер, через который проходят 2757 управляющие объекты. 2758 2759 #+name: fig-pipeline 2760 #+begin_src dot :exports results :file build/pipeline-ru.pdf 2761 digraph { 2762 2763 node [fontname="Old Standard"fontsize=14,margin="0.055,0"] 2764 graph [nodesep="0.25",ranksep="0.25",rankdir="TB",margin=0] 2765 edge [arrowsize=0.66] 2766 2767 # data 2768 subgraph xcluster_linear { 2769 label="Линейная модель" 2770 2771 start [label="",shape=circle,style=filled,fillcolor=black,width=0.23] 2772 spectrum [label="S(ω,θ)",shape=box] 2773 acf [label="K(i,j,k)",shape=box] 2774 phi [label="Φ(i,j,k)",shape=box] 2775 2776 # transformations 2777 fourier_transform [label="Преобразование Фурье",shape=box,style=rounded] 2778 solve_yule_walker [label="Решение уравнений\nЮла—Уокера",shape=box,style=rounded] 2779 2780 subgraph cluster_nonlinear_1 { 2781 label="Моделир. нелинейности\l" 2782 labeljust=left 2783 style=filled 2784 color=lightgrey 2785 graph [margin=8] 2786 acf2 [label="K*(i,j,k)",shape=box] 2787 transform_acf [label="Преобразование АКФ",shape=box,style=rounded] 2788 } 2789 } 2790 2791 subgraph xcluster_linear2 { 2792 2793 eps_parts [label="<e1> ε₁|<e2> ε₂|<e3> …|<e4> εₙ|<e> ε(t,x,y)",shape=record] 2794 end [label="",shape=doublecircle,style=filled,fillcolor=black,width=0.23] 2795 2796 generate_white_noise [label="<g1> g₁|<g2> g₂|<g3> …|<g4> gₙ|<gen> Генерация\lбелого шума",shape=record,style=rounded] 2797 generate_zeta [label="<g1> g₁|<g2> g₂|<g3> …|<g4> gₙ|<gen> Генерация частей\lвзволнованной мор-\lской поверхности\l",shape=record,style=rounded] 2798 2799 zeta_parts [label="<g1> ζ₁|<g2> ζ₂|<g3> …|<g4> ζₙ|<gen> Несшитые части реализации",shape=record] 2800 overlap_add [label="<g1> ζ₁|<g2> ζ₂|<g3> …|<g4> ζₙ|<gen> Сшивание час-\lтей реализации\l",shape=record,style=rounded] 2801 2802 zeta_parts:g1->overlap_add:g1 2803 zeta_parts:g2->overlap_add:g2 2804 zeta_parts:g3->overlap_add:g3 2805 zeta_parts:g4->overlap_add:g4 2806 2807 zeta_parts:g2->overlap_add:g1 [constraint=false] 2808 zeta_parts:g3->overlap_add:g2 [constraint=false] 2809 zeta_parts:g4->overlap_add:g3 [constraint=false] 2810 2811 overlap_add:g1->zeta2_parts:g1 2812 overlap_add:g2->zeta2_parts:g2 2813 overlap_add:g3->zeta2_parts:g3 2814 overlap_add:g4->zeta2_parts:g4 2815 2816 zeta2_parts:g1->transform_zeta:g1->zeta3_parts:g1->write_zeta:g1->eps_end 2817 zeta2_parts:g2->transform_zeta:g2->zeta3_parts:g2->write_zeta:g2->eps_end 2818 zeta2_parts:g3->transform_zeta:g3->zeta3_parts:g3->write_zeta:g3->eps_end 2819 zeta2_parts:g4->transform_zeta:g4->zeta3_parts:g4->write_zeta:g4->eps_end 2820 2821 } 2822 2823 subgraph part3 { 2824 2825 zeta2_parts [label="<g1> ζ₁|<g2> ζ₂|<g3> …|<g4> ζₙ|<gen> Поверхность с нормаль-\lным законом распреде-\lления\l",shape=record] 2826 2827 subgraph cluster_nonlinear_2 { 2828 label="Моделир. нелинейности\r" 2829 labeljust=right 2830 style=filled 2831 color=lightgrey 2832 graph [margin=8] 2833 zeta3_parts [label="<g1> ζ₁|<g2> ζ₂|<g3> …|<g4> ζₙ|<gen> ζ(t,x,y)",shape=record] 2834 transform_zeta [label="<g1> g₁|<g2> g₂|<g3> …|<g4> gₙ|<gen> Преобразование за-\lкона распределения\lвзволнованной мор-\lской поверхности\l",shape=record,style=rounded] 2835 } 2836 2837 # barriers 2838 eps_start [label="",shape=box,style=filled,fillcolor=black,height=0.05] 2839 eps_end [label="",shape=box,style=filled,fillcolor=black,height=0.05] 2840 2841 write_zeta [label="<g1> g₁|<g2> g₂|<g3> …|<g4> gₙ|<gen> Запись готовых\lчастей в файл\l",shape=record,style=rounded] 2842 } 2843 2844 # edges 2845 start->spectrum->fourier_transform->acf->transform_acf 2846 transform_acf->acf2 2847 acf2->solve_yule_walker 2848 solve_yule_walker->phi 2849 phi->eps_start [constraint=false] 2850 eps_start->generate_white_noise:g1 2851 eps_start->generate_white_noise:g2 2852 eps_start->generate_white_noise:g3 2853 eps_start->generate_white_noise:g4 2854 generate_white_noise:g1->eps_parts:e1->generate_zeta:g1->zeta_parts:g1 2855 generate_white_noise:g2->eps_parts:e2->generate_zeta:g2->zeta_parts:g2 2856 generate_white_noise:g3->eps_parts:e3->generate_zeta:g3->zeta_parts:g3 2857 generate_white_noise:g4->eps_parts:e4->generate_zeta:g4->zeta_parts:g4 2858 2859 eps_end->end 2860 } 2861 #+end_src 2862 2863 #+caption[Схема конвейера обработки данных]: 2864 #+caption: Схема конвейера обработки данных, реализующего генерацию 2865 #+caption: взволнованной морской поверхности с помощью АР модели. 2866 #+name: fig-pipeline 2867 #+RESULTS: fig-pipeline 2868 [[file:build/pipeline-ru.pdf]] 2869 2870 Конвейер объектов можно считать развитием модели Bulk Synchronous Parallel 2871 (BSP)\nbsp{}cite:valiant1990bridging, применяемой в системах обработки 2872 графов\nbsp{}cite:malewicz2010pregel,seo2010hama. Конвейер позволяет исключить 2873 глобальную синхронизацию (где это возможно) между последовательно идущим этапами 2874 вычислений путем передачи данных между звеньями параллельно с вычислениями, в то 2875 время как в модели BSP глобальная синхронизация происходит после каждого шага. 2876 2877 Конвейер объектов ускоряет программу путем параллельного выполнения блоков кода, 2878 работающих с разными вычислительными устройствами: в то время как текущая часть 2879 взволнованной поверхности генерируется на процессоре, предыдущая часть 2880 записывается на диск. Такой подход позволяет получить ускорение 2881 (см.\nbsp{}разд.\nbsp{}[[#sec-io-perf]]), потому что различные вычислительные 2882 устройства работают асинхронно, и их параллельное использование увеличивает 2883 производительность программы. 2884 2885 Поскольку передача данных между звеньями конвейера происходит параллельно с 2886 вычислениями, то на одном и том же конвейере можно запустить сразу несколько 2887 копий приложения с разными параметрами (генерировать сразу несколько 2888 взволнованных морских поверхностей с разными характеристиками). На практике 2889 оказывается, что высокопроизводительные приложения не всегда загружают 2890 процессор на 100%, тратя время на синхронизацию параллельных процессов и 2891 запись данных на диск. Использование конвейера в таком случае позволит на одном 2892 и том же множестве процессов запустить сразу несколько расчетов и максимально 2893 эффективно использовать все устройства компьютера. Например, во время записи в 2894 файл одной задачей может производиться расчет на процессоре другой задачей. Это 2895 минимизирует время простоя процессора и других устройств компьютера и повышает 2896 общую пропускную способность кластера. 2897 2898 Конвейеризация шагов программы, которые в противном случае последовательны, 2899 выгодно не только для кода, работающего с различными устройствами, но и для 2900 кода, различные ветки которого могут быть запущены на нескольких аппаратных 2901 потоках одного процессорного ядра, т.е.\nbsp{}ветки, осуществляющие доступ к 2902 различным блокам памяти или использующие смешанную арифметику (целочисленную и с 2903 плавающей точкой). Ветки кода, которые используют различные модули центрального 2904 процессора, являются подходящими кандидатами для параллельного запуска на 2905 процессорном ядре с несколькими аппаратными потоками. 2906 2907 Таким образом, вычислительную модель на основе конвейера можно рассматривать как 2908 /массивно асинхронную модель/ (bulk-asynchronous model) из-за параллельной 2909 природы шагов программы. Эта модель является основой модели отказоустойчивости, 2910 которая будет описана далее. 2911 2912 *** Обнаружение узлов кластера 2913 :PROPERTIES: 2914 :CUSTOM_ID: sec-node-discovery 2915 :END: 2916 2917 **** Алгоритмы выбора лидера. 2918 Многие системы пакетной обработки задач построены по принципу /субординации/: в 2919 каждом кластере выбирается главный узел, который управляет очередью задач, 2920 планирует их запуск на подчиненных узлах и следит за их состоянием. Роль 2921 главного узла назначается либо /статически/ системным администратором 2922 определенному физическому узлу, либо /динамически/, путем периодического 2923 избрания какого-либо из узлов кластера главным. В первом случае 2924 отказоустойчивость обеспечивается посредством резервирования дополнительного 2925 свободного узла, который выполнит роль главного в случае отказа текущего. Во 2926 втором случае отказоустойчивость обеспечивается выбором нового главного узла из 2927 оставшихся. Несмотря на то что динамическое задание ролей требует наличия 2928 алгоритма выбора лидера, этот подход становится все более и более популярным, 2929 поскольку не требует наличия простаивающих резервных узлов на случай отказа 2930 главного 2931 узла\nbsp{}cite:hunt2010zookeeper,lakshman2010cassandra,divya2013elasticsearch и 2932 в общем случае приводит к симметричной архитектуре системы, в которой один и тот 2933 же стек программного обеспечения с одними и теми же настройками установлен на 2934 каждом узле кластера\nbsp{}cite:boyer2012glusterfs,ostrovsky2015couchbase. 2935 2936 Алгоритмы выбора лидера (которые иногда называют алгоритмами /распределенного 2937 консенсуса/) являются частными случаями волновых алгоритмов. 2938 В\nbsp{}cite:tel2000introduction автор определяет их как алгоритмы, в которых 2939 событие завершения программы предваряется хотя бы одним каким-либо другим 2940 событием, происходящем в /каждом/ параллельном процессе. Волновые алгоритмы не 2941 определены для анонимных сетей, т.е.\nbsp{}они работают только с теми 2942 параллельными процессами, которые могут себя уникально идентифицировать. Однако, 2943 количество процессов, которые затрагивает "волна", может быть определено по мере 2944 выполнения алгоритма. В рамках распределенных систем это означает, что волновые 2945 алгоритмы подходят для вычислительных кластеров с динамически меняющимся 2946 количеством узлов, так что включение и выключение отдельных узлов не влияет на 2947 работу алгоритма. 2948 2949 Подход к поиску узлов кластера не использует волновые алгоритмы, а значит не 2950 требует опроса всех узлов кластера для выбора лидера. Вместо этого каждый узел 2951 кластера нумерует все узлы подсети, в которой он находится, и преобразует список 2952 в /древовидную иерархию/ с заданным значением /ветвления/ (максимальным 2953 количеством подчиненных вершин, которых может иметь узел). Затем узел определяет 2954 свой уровень иерархии и пытается соединиться с вышестоящими узлами, чтобы стать 2955 их подчиненным. Сначала он проверяет близко расположенные к нему узлы, а потом 2956 все остальные узлы вплоть до вершины иерархии. Если вышестоящих узлов нет или с 2957 ними невозможно соединиться, то узел сам становится главой всей иерархии. 2958 2959 Древовидная иерархия узлов подсети определяет отношение строгого порядка на 2960 множестве всех узлов кластера. Несмотря на то что технически любая функция может 2961 быть выбрана для присвоения узлу подсети номера в списке, на практике эта 2962 функция должна быть достаточно гладкой вдоль временной оси и иметь лишь редкие 2963 скачки: быстрые изменения в структуре иерархии узлов (которые часто являются 2964 следствием погрешности измерений) могут привести постоянной передаче роли 2965 главного узла от одного узла к другому, что сделает иерархию непригодной для 2966 распределения нагрузки. Простейшей такой функцией является позиция IP-адреса 2967 узла в диапазоне всех IP-адресов подсети. 2968 2969 **** Алгоритм создания древовидной иерархии. 2970 Отношение строго порядка на множестве \(\mathcal{N}\) узлов одной подсети 2971 определяется как 2972 \begin{equation*} 2973 \forall n_1 \forall n_2 \in \mathcal{N}, 2974 \forall f \colon \mathcal{N} \rightarrow \mathcal{R}^n 2975 \Rightarrow (f(n_1) < f(n_2) \Leftrightarrow \neg (f(n_1) \geq f(n_2))), 2976 \end{equation*} 2977 где \(f\)\nbsp{}--- отображение узла на его ранг, а \(<\)\nbsp{}--- оператор, 2978 определяющий отношение строго порядка на множестве \(\mathcal{R}^n\). Функция 2979 \(f\) присваивает узлу порядковый номер, а оператор \(<\) делает этот номер 2980 уникальным. 2981 2982 Простейшее отображение \(f\) ставит в соответствие каждому узлу подсети позицию 2983 его IP-адреса в диапазоне всех адресов подсети. Без использования древовидной 2984 иерархии (когда в подсети выбирается только один лидер) узел, адрес которого 2985 занимает наименьшую позицию в диапазоне, становится главным. Если адрес узла 2986 занимает первую позицию в диапазоне, то для него невозможно выбрать лидера, и он 2987 будет находится на вершине иерархии вплоть до выхода из строя. Несмотря на то 2988 что идентификацию узлов на основе их IP-адресов легко реализовать в программе, 2989 такой подход устанавливает искусственную зависимость роли главного узла от 2990 IP-адреса узла. Тем не менее, подход полезен для первичного объединения узлов в 2991 кластер, когда более сложные отображения неприменимы. 2992 2993 Для того чтобы алгоритм обнаружения масштабировался на большое количество узлов, 2994 диапазон IP адресов подсети отображается на древовидную иерархию. В такой 2995 иерархии каждый узел определяется уровнем \(l\) иерархии, на котором он 2996 находится, и отступом \(o\), который равен порядковому номеру узла на его 2997 уровне. Значения уровня и отступа определяются из следующей задачи оптимизации. 2998 \begin{equation*} 2999 n = \sum\limits_{i=0}^{l(n)} p^i + o(n), \quad 3000 l \rightarrow \min, \quad 3001 o \rightarrow \min, \quad 3002 l \geq 0, \quad 3003 o \geq 0 3004 \end{equation*} 3005 где \(n\)\nbsp{}--- позиция IP-адреса узла в диапазоне IP-адресов подсети и 3006 \(p\)\nbsp{}--- значение ветвления. Главный узел для узла на уровне \(l\) с 3007 отступом \(o\) имеет уровень \(l-1\) и отступ \(\lfloor{o/p}\rfloor\). 3008 Расстояние между любыми двумя узлами в иерархии, адреса которых занимают позиции 3009 \(i\) и \(j\) в диапазоне определяется как 3010 \begin{align*} 3011 & \langle 3012 \text{lsub}(l(j), l(i)), \quad 3013 \left| o(j) - o(i)/p \right| 3014 \rangle,\\ 3015 & \text{lsub}(l_1, l_2) = 3016 \begin{cases} 3017 \infty & \quad \text{if } l_1 \geq l_2, \\ 3018 l_1 - l_2 & \quad \text{if } l_1 < l_2. 3019 \end{cases} 3020 \end{align*} 3021 Расстояние является составным, чтобы уровень иерархии учитывался в первую 3022 очередь. 3023 3024 Для выбора главного каждый узел ранжирует все узлы подсети в соответствии с их 3025 позицией \(\langle{l(n),o(n)}\rangle\) и, используя формулу для определения 3026 расстояния, выбирает ближайший к потенциальному главному узлу узел, имеющий 3027 наименьший ранг. Это позволяет пропустить IP-адреса выключенных узлов, однако, 3028 для разреженных сетей (в которых узлы занимают непоследовательные IP-адреса) 3029 сбалансированность дерева не гарантируется. 3030 3031 Формулу для вычисления расстояния можно сделать сколь угодно сложной (например, 3032 чтобы учесть задержки и пропускную способность сети или географическое 3033 местоположение узла), однако, для ее простейшего вида более выгодно использовать 3034 /алгоритм обхода/ узлов кластера. Этот алгоритм требует меньшего количества 3035 памяти, поскольку не нужно хранить ранжированный список всех узлов, вместо этого 3036 он перебирает все IP-адреса сети в порядке, определяемом значением ветвления. 3037 Алгоритм работает следующим образом. Сначала базовый узел (узел, который ищет 3038 главного) вычисляет адрес своего потенциального главного узла и пытается 3039 установить соединение с ним. Если соединение не удается, базовый узел 3040 последовательно пытается соединиться с каждым узлом, находящимся на более 3041 высоком уровне иерархии, пока не достигнет вершины иерархии (корня дерева). Если 3042 ни одно из соединений не удается, базовый узел последовательно пытается 3043 соединиться с каждым узлом на своем уровне, имеющим более низкую позицию в 3044 диапазоне всех IP-адресов подсети. Если ни один из узлов не отвечает, базовый 3045 узел занимает вершину иерархии, а обход иерархии повторяется через заданный 3046 промежуток времени. Пример порядка обхода для кластера из 11 узлов и древовидной 3047 иерархии с ветвлением 2 показан на рис.\nbsp{}[[fig-tree-hierarchy-11]]. 3048 3049 #+name: fig-tree-hierarchy-11 3050 #+begin_src dot :exports results :file build/tree-hierarchy-11-ru.pdf 3051 digraph { 3052 3053 node [fontname="Old Standard",fontsize=14,margin="0.055,0.055",shape=box,style=rounded] 3054 graph [nodesep="0.30",ranksep="0.30",rankdir="BT"] 3055 edge [arrowsize=0.66] 3056 3057 m1 [label="10.0.0.1"] 3058 m2 [label="10.0.0.2"] 3059 m3 [label="10.0.0.3"] 3060 m4 [label="10.0.0.4"] 3061 m5 [label="10.0.0.5"] 3062 m6 [label="10.0.0.6"] 3063 m7 [label="10.0.0.7"] 3064 m8 [label="10.0.0.8"] 3065 m9 [label="10.0.0.9"] 3066 m10 [label="10.0.0.10",fillcolor="#c0c0c0",style="filled,rounded"] 3067 m11 [label="10.0.0.11",shape=Mrecord] 3068 3069 m2->m1 3070 m3->m1 3071 m4->m2 3072 m5->m2 3073 m6->m3 3074 m7->m3 3075 m8->m4 3076 m9->m4 3077 m10->m5 3078 m11->m5 3079 3080 m5->m4->m6 [style="dashed,bold",color="#c00000",constraint=false] 3081 {rank=same; m6->m7 [style="dashed,bold",color="#c00000"]} 3082 m7->m2->m3->m1->m8->m9 [style="dashed,bold",color="#c00000",constraint=false] 3083 3084 } 3085 #+end_src 3086 3087 #+caption[Древовидная иерархия для кластера из 11 узлов]: 3088 #+caption: Древовидная иерархия для кластера из 11 узлов со 3089 #+caption: значением ветвления 2. Красными стрелками обозначен порядок 3090 #+caption: обхода иерархии узлом с IP-адресом 10.0.0.10. 3091 #+name: fig-tree-hierarchy-11 3092 #+RESULTS: fig-tree-hierarchy-11 3093 [[file:build/tree-hierarchy-11-ru.pdf]] 3094 3095 **** Результаты тестирования. 3096 Для тестирования производительности алгоритма обхода на большом количестве 3097 узлов, на каждом физическом узле кластера запускалось несколько резидентных 3098 процессов, каждый из которых был привязан к отдельному IP-адресу. Количество 3099 процессов на одно физическое ядро варьировалось от 2 до 16. Каждый процесс был 3100 привязан к определенному физическому ядру, чтобы уменьшить накладные расходы на 3101 миграцию процессов между ядрами. Алгоритм имеет низкие требования к 3102 процессорному времени и пропускной способности сети, поэтому запуск нескольких 3103 процессов на одном физическом ядре целесообразен, в отличие от кодов 3104 высокопроизводительных приложений, в которых это часто снижает 3105 производительность. Конфигурация тестовой системы показана в 3106 табл.\nbsp{}[[tab-ant]]. 3107 3108 #+name: tab-ant 3109 #+caption[Конфигурация системы "Ant"]: 3110 #+caption: Конфигурация системы "Ant". 3111 #+attr_latex: :booktabs t 3112 | Процессор | Intel Xeon E5440, 2,83ГГц | 3113 | Память | 4ГБ | 3114 | Жесткий диск | ST3250310NS, 7200об./мин. | 3115 | Количество узлов | 12 | 3116 | Количество ядер на узел | 8 | 3117 3118 Похожий подход используется 3119 в\nbsp{}cite:lantz2010network,handigol2012reproducible,heller2013reproducible, 3120 где авторы воспроизводят разнообразные практические эксперименты на виртуальных 3121 кластерах, созданных на основе пространств имен Linux, и сопоставляют результаты 3122 с физическими. Преимущество данного подхода заключается в возможности проведения 3123 экспериментов на больших виртуальных кластерах, используя сравнительно небольшое 3124 количество физических узлов. Преимущество же подхода, используемого в данной 3125 работе (в котором не применяются пространства имен Linux) заключается в том, что 3126 он более легковесный и большее количество резидентных процессов можно 3127 протестировать на одном и том же физическом кластере. 3128 3129 Производительность алгоритма обхода была протестирована путем измерения времени, 3130 которое необходимо для того чтобы все узлы кластера нашли друг друга, 3131 т.е.\nbsp{}времени, которое необходимо для того чтобы древовидная иерархии узлов 3132 достигла устойчивого состояния. Каждое изменение иерархии, то, как его видит 3133 каждый узел, записывалось в файл журнала, и по прошествии заданного промежутка 3134 времени все резидентные процессы (каждый из которых моделирует узел кластера) 3135 принудительно завершались. Процессы запускались последовательно с задержкой в 3136 100мс., чтобы удостовериться, что главные узлы становятся доступны раньше 3137 подчиненных, а иерархия не меняется произвольным образом в результате разного 3138 времени запуска процессов. Эта искусственная задержка впоследствии была вычтена 3139 из результатов тестирования. Таким образом, результаты теста представляют собой 3140 время обнаружения узлов в "идеальном" кластере, в котором каждый резидентный 3141 процесс находит главного с первой попытки. 3142 3143 Тест запускался несколько раз, варьируя количество резидентных процессов на 3144 физический узел. Эксперимент показал, что обнаружение 512 узлами (8 физических 3145 узлов по 64 процесса на узел) друг друга занимает не более двух секунд 3146 (рис.\nbsp{}[[fig-discovery-benchmark]]). Это значение меняется незначительно с 3147 увеличением количества физических узлов. Использование более 8 узлов по 64 3148 процесса на узел приводит к большим колебаниям времени обнаружения, ввиду того 3149 что большое количество процессов одновременно устанавливает соединение с одним и 3150 тем же главным процессом (значение ветвления во всех тестах равнялось 10000), 3151 поэтому эти результаты были исключены из рассмотрения. 3152 3153 #+name: fig-discovery-benchmark 3154 #+header: :width 7 :height 5 3155 #+begin_src R :file build/discovery-benchmark-ru.pdf 3156 source(file.path("R", "discovery.R")) 3157 par(family="serif") 3158 bscheduler.plot_discovery( 3159 xlabel="Количество физических узлов", 3160 ylabel="Время, сек.", 3161 toplabel="ppn") 3162 #+end_src 3163 3164 #+caption[Время обнаружения всех резидентных процессов на кластере]: 3165 #+caption: Время обнаружения всех резидентных процессов, запущенных на кластере, 3166 #+caption: в зависимости от количества процессов. Пунктирная линия обозначает 3167 #+caption: минимальное и максимальное значение за все проведенные тесты, 3168 #+caption: "ppn"\nbsp{}--- количество резидентных процессов на узел. 3169 #+name: fig-discovery-benchmark 3170 #+RESULTS: fig-discovery-benchmark 3171 [[file:build/discovery-benchmark-ru.pdf]] 3172 3173 **** Обсуждение. 3174 Поскольку узлу для выбора главного нужно соединиться с узлом, адрес которого 3175 известен заранее, то алгоритм обхода масштабируется на большое количество узлов. 3176 Соединение с другими узлами происходит только в том случае, если текущий главный 3177 узел выходит из строя. Таким образом, если адреса узлов кластера расположены 3178 непрерывно в диапазоне адресов подсети, каждый узел устанавливает соединение 3179 только со своим главным узлом, и неэффективного сканирования всей сети каждым 3180 узлом не происходит. 3181 3182 Следующие ключевые особенности отличают предложенный подход от некоторых 3183 существующих\nbsp{}cite:brunekreef1996design,aguilera2001stable,romano2014design. 3184 - *Многоуровневая иерархия.* Количество главных узлов в сети зависит от значения 3185 ветвления. Если оно меньше количества IP-адресов в подсети, то в кластере 3186 будет несколько главных узлов. Если оно больше или равно количеству IP-адресов 3187 в подсети, то в кластере будет только один главный узел. Когда какой-либо узел 3188 выходит из строя, многоуровневая иерархия изменятся локально, только узлы, 3189 примыкающие к вышедшему из строя, взаимодействуют друг с другом. 3190 - *Отображение IP-адресов.* Поскольку структура иерархии зависит только от 3191 IP-адресов узлов, то в алгоритме отсутствует фаза выбора лидера. Чтобы сменить 3192 главного, каждый узел отправляет сообщение только прежнему и новому главному 3193 узлу. 3194 - *Полностью основан на событиях.* Сообщения отправляются только при введении 3195 нового узла в кластер или при выходе из строя существующего, поэтому 3196 постоянной нагрузки на сеть нету. Поскольку алгоритм допускает ошибку при 3197 отправке любого сообщения, то нет необходимости в heartbeat-пакетах, 3198 являющихся индикацией работоспособности узла в сети; вместо этого все 3199 сообщения выполняют роль heartbeat-пакетов и настраивается время ожидания 3200 отправки пакета\nbsp{}cite:rfc5482. 3201 - *Отсутствие ручной конфигурации.* Узлу не требуется никаких предварительных 3202 знаний, чтобы найти главного: он определяет сеть, узлом которой он является, 3203 вычисляет IP-адрес потенциального главного узла и отправляет ему сообщение. 3204 Если это не срабатывает, то процесс повторяется для следующего потенциального 3205 главного узла. Таким образом, алгоритм способен выполнить начальную загрузку 3206 кластера (сделать так, чтобы все узлы узнали друг о друге) без предварительной 3207 ручной настройки, для этого требуется только назначить IP-адрес каждому узлу и 3208 запустить резидентный процесс на нем. 3209 Суммируя вышесказанное, достоинством алгоритма является то, что он 3210 - масштабируется на большое количество узлов посредством иерархии с несколькими 3211 главными узлами, 3212 - не нагружает сеть отправкой сообщений с текущим состоянием узлов и 3213 heartbeat-пакетами, 3214 - не требует ручной настройки для первичной загрузки кластера. 3215 3216 Недостатком алгоритма является то, что он требует, чтобы IP-адрес изменялся 3217 редко, чтобы быть полезным для распределения нагрузки. Он не подходит для 3218 облачной среды, в которой только DNS имя узла сохраняется, а IP-адрес может 3219 меняться со временем. Когда IP-адрес меняется, текущие соединения могут 3220 закрыться, сигнализируя о "выходе из строя" узла и перестраивая иерархию узлов. 3221 Таким образом, окружения, в которых узлы не идентифицируются IP-адресами, не 3222 подходят для алгоритма. 3223 3224 Другим недостатком алгоритма является искусственная зависимость ранга узла от 3225 IP-адреса: замена отображения IP-адресов на что-то более сложное. Если 3226 отображение использует загрузку текущего узла и сети для ранжирования узлов, то 3227 погрешность измерений может стать причиной неустойчивой иерархии, а алгоритм 3228 перестанет быть полностью основан на событиях, поскольку уровни загрузки 3229 необходимо измерять периодически на каждом узле и распространять на все 3230 остальные узлы кластера. 3231 3232 Алгоритм обнаружения узлов спроектирован для балансировки нагрузки на кластер 3233 вычислительных узлов (см.\nbsp{}разд.\nbsp{}[[#sec-daemon-layer]]), и его применение 3234 в других приложениях не рассматривается в данной работе. Когда распределенная 3235 или параллельная программа запускается на одном из узлов кластера, ее подзадачи 3236 распределяются между всеми примыкающими узлами иерархии (включая главный узел, 3237 если есть). Для того чтобы равномерно распределить нагрузку при запуске 3238 программы на подчиненном узле, каждый узел хранит вес каждого из примыкающих 3239 узлов иерархии. Вес равен количеству узлов дерева, находящегося "за" примыкающим 3240 узлом. Например, если вес первого примыкающего узла равен 2, то циклический 3241 алгоритм балансировки нагрузки распределит две подзадачи на первый узел перед 3242 тем как перейти к следующему узлу. 3243 3244 Суммируя вышесказанное, алгоритм обхода 3245 - спроектирован для облегчения распределения нагрузки на кластер, 3246 - полностью отказоустойчивый, состояние каждого узла можно вычислить заново в 3247 любой момент времени, 3248 - полностью основан на событиях, а значит не нагружает сеть периодической 3249 отправкой сообщений. 3250 3251 *** Алгоритм восстановления после сбоев 3252 :PROPERTIES: 3253 :CUSTOM_ID: sec-fail-over 3254 :END: 3255 3256 **** Контрольные точки восстановления. 3257 Сбои узлов распределенной системы можно разделить на два типа: сбой подчиненного 3258 узла и сбой главного узла. Для того чтобы запущенная на кластере задача могла 3259 пережить сбой подчиненного узла, планировщик задач периодически создает для нее 3260 контрольные точки восстановления и записывает их в надежное (резервируемое) 3261 хранилище. Для того чтобы создать контрольную точку, планировщик временно 3262 останавливает все параллельные процессы задачи, копирует все страницы памяти и 3263 все структуры ядра операционной системы, выделенные для этих процессов, на диск, 3264 и продолжает выполнение задачи. Для того чтобы пережить сбой главного узла, 3265 резидентный процесс планировщика задач непрерывно копирует свое внутреннее 3266 состояние на резервный узел, который становится главным после сбоя. 3267 3268 Оптимизации работы контрольных точек восстановления посвящено большое количество 3269 работ\nbsp{}cite:egwutuoha2013survey, а альтернативным подходам уделяется меньше 3270 внимания. Обычно высокопроизводительные приложения используют передачу сообщений 3271 для обмена данными между параллельными процессами и хранят свое текущее 3272 состояние в глобальной памяти, поэтому не существует способа перезапустить 3273 завершившийся процесс, не записав образ всей выделенной для него памяти на диск. 3274 Обычно общее число процессов фиксировано и задается планировщиком, и в случае 3275 отказа перезапускаются сразу все процессы. Существуют некоторые обходные 3276 решения, которые позволяют перезапустить только часть 3277 процессов\nbsp{}cite:meyer2012radic, восстановив их из контрольной точки на 3278 выживших узлах, однако это может привести к перегрузке, если на этих узлах уже 3279 запущены другие задачи. Теоретически, перезапуск процесса необязателен если 3280 задача может быть продолжена на выживших узлах, но библиотека передачи сообщений 3281 не позволяет изменять количество параллельных процессов во время работы 3282 программы, и большинство программ все равно предполагают, что это значение 3283 является константой. Таким образом, не существует надежного способа обеспечения 3284 отказоустойчивости на уровне библиотеки передачи сообщений, кроме как путем 3285 перезапуска всех параллельных процессов из контрольной точки восстановления. 3286 3287 В то же время, существует возможность продолжить выполнение задачи на меньшем 3288 количестве узлов, чем было изначально выделено изначально, реализовав 3289 отказоустойчивость на уровне планировщика задач. В этом случае роли главного и 3290 подчиненного динамически распределяются между резидентными процессами 3291 планировщика, работающими на каждом узле кластера, образуя древовидную иерархию 3292 узлов кластера, а параллельная программа состоит из управляющих объектов, 3293 использующих иерархию узлов для динамического распределения нагрузки и свою 3294 собственную иерархию для перезапуска управляющих объектов в случае сбоя узла. 3295 3296 **** Динамическое распределение ролей. 3297 Отказоустойчивость параллельной программы\nbsp{}--- это одна из проблем, которая 3298 решается планировщиками задач обработки больших данных или 3299 высокопроизводительных вычислений, однако, большинство планировщиков 3300 обеспечивают только отказоустойчивость подчиненных узлов. Такого рода сбои 3301 обычно обрабатываются путем перезапуска затронутой задачи (из контрольной точки 3302 восстановления) или ее части на оставшихся узлах, а выход из строя главного узла 3303 считается либо маловероятным, либо слишком сложным для обработки и настройки на 3304 целевой платформе. Системные администраторы обычно находят альтернативы 3305 отказоустойчивости на уровне приложения: они изолируют главный процесс 3306 планировщика от остальных узлов кластера, размещая его на специально выделенной 3307 машине, или, вместо этого, используют технологии виртуализации. Все эти 3308 альтернативы усложняют конфигурацию и обслуживание, и, уменьшая вероятность 3309 выхода из строя машины, приводящей к выходу из строя всей системы, увеличивают 3310 вероятность ошибки оператора. 3311 3312 С этой точки зрения более практично реализовать отказоустойчивость главного узла 3313 на уровне приложения, однако, не существует универсального зарекомендовавшего 3314 себя решения. Большинство реализаций слишком привязаны к конкретному приложению, 3315 чтобы стать повсеместно применяемыми. Часто кластер представляется, как машина, 3316 в которой есть управляющий модуль (главный узел), отличный от всех остальных, и 3317 несколько одинаковых вычислительных модулей (подчиненных узлов). На практике же 3318 оказывается, что главный и подчиненные узлы физически одинаковы, и различаются 3319 лишь процессами планировщика пакетных задач, запущенными на них. Понимание того, 3320 что узлы кластера являются лишь вычислительными модулями, позволяет реализовать 3321 промежуточное программное обеспечение, которое автоматически распределяет роли 3322 главного и подчиненного узла и универсальным способом обрабатывает сбои узлов. 3323 Это программное обеспечение предоставляет программный интерфейс и распределяет 3324 управляющие объекты между доступными на данный момент узлами. Используя этот 3325 интерфейс, можно написать программу, которая запускается на кластере, не зная 3326 точного количества работающих узлов. Это промежуточное программное обеспечение 3327 работает как кластерная операционная система в пользовательском пространстве, 3328 позволяющая запускать распределенные приложения прозрачно на любом количестве 3329 узлов. 3330 3331 **** Симметричная архитектура. 3332 Многие распределенные хранилища типа ключ-значение и параллельные файловые 3333 системы имеют симметричную архитектуру, в которой роли главного и подчиненного 3334 распределяются динамически, так что любой узел может выступать в роли главного, 3335 если текущий главный узел выходит из 3336 строя\nbsp{}cite:ostrovsky2015couchbase,divya2013elasticsearch,boyer2012glusterfs,anderson2010couchdb,lakshman2010cassandra. 3337 Однако, такая архитектура до сих пор не используется в планировщиках задач 3338 обработки больших данных и высокопроизводительных вычислений. Например, в 3339 планировщике YARN\nbsp{}cite:vavilapalli2013yarn, роли главного и подчиненного 3340 являются статическими. Восстановление после сбоя подчиненного узла 3341 осуществляется путем перезапуска работавшей на нем части задачи на одном из 3342 выживших узлов, а восстановление после сбоя главного узла осуществляется путем 3343 установки резервного главного узла\nbsp{}cite:murthy2011architecture. Оба 3344 главных узла управляются сервисом Zookeeper, который использует динамическое 3345 распределение ролей для обеспечения своей 3346 отказоустойчивости\nbsp{}cite:okorafor2012zookeeper. Таким образом, отсутствие 3347 динамического распределения ролей у планировщика YARN усложняет конфигурацию 3348 всего кластера: если бы динамические роли были доступны, Zookeeper был бы лишним 3349 в данной конфигурации. 3350 3351 Такая же проблема возникает в планировщиках задач для высокопроизводительных 3352 вычислений, главный узел (на котором запущен главный процесс планировщика задач) 3353 является единой точкой сбоя. 3354 В\nbsp{}cite:uhlemann2006joshua,engelmann2006symmetric авторы копируют состояние 3355 планировщика задач на резервный узел, чтобы обеспечить высокую доступность 3356 главного узла, но роль резервного узла задается статически. Такое решение близко 3357 к симметричной архитектуре, поскольку не использует внешний сервис для 3358 обеспечения высокой доступности, но далеко от идеала, в котором резервный узел 3359 выбирается динамически. 3360 3361 Наконец, наиболее простой вариант высокой доступности главного узла 3362 реализован в протоколе VRRP (Virtual Router Redundancy 3363 Protocol)\nbsp{}cite:rfc2338,rfc3768,rfc5798. Несмотря на то что протокол VRRP 3364 предоставляет динамическое распределение ролей, он не может быть использован в 3365 планировщиках задач, поскольку спроектирован для маршрутизаторов, за которыми 3366 стоят инвертированные прокси-серверы. В таких серверах отсутствует состояние 3367 (очередь задач), которое необходимо восстановить после выхода из строя узла, 3368 поэтому их высокую доступность обеспечить проще. Это может быть реализовано даже 3369 без маршрутизаторов, используя вместо этого сервис 3370 Keepalived\nbsp{}cite:cassen2002keepalived. 3371 3372 Симметричная архитектура выгодна для планировщиков задач, поскольку позволяет 3373 - сделать физические узлы взаимозаменяемыми, 3374 - реализовать динамическое распределение ролей главного и подчиненного узла и 3375 - реализовать автоматическое восстановление после сбоя любого из узлов. 3376 В последующих разделах описаны компоненты необходимые для написания параллельной 3377 программы и планировщика, которые устойчивы к сбоям узлов кластера. 3378 3379 **** Определения иерархий. 3380 Для устранения неоднозначности иерархических связей между резидентными 3381 процессами и управляющими объектами и для того чтобы упростить изложение, в 3382 тексте используются следующие условные обозначения. Если связь установлена между 3383 двумя резидентными процессами, то отношения обозначаются /главный-подчиненный/. 3384 Если связь установлена между двумя управляющими объектами, то отношения 3385 обозначаются либо /руководитель-подчиненный/, либо /родитель-потомок/. Две 3386 иерархии ортогональны друг к другу в том смысле, что ни один управляющий объект 3387 не может иметь связь с резидентным процессом, и наоборот. Поскольку иерархия 3388 резидентных процессов используется для распределения нагрузки на узлы кластера, 3389 иерархия управляющих объектов отображается на нее, и это отображение может быть 3390 произвольным: обычна ситуация, когда руководящий управляющий объект находится на 3391 подчиненном узле, а его подчиненные управляющие объекты распределены равномерно 3392 между всеми узлами кластера (включая узел, где находится руководящий объект). 3393 Обе иерархии может быть сколь угодно глубокими, но "неглубокие" являются 3394 предпочтительными для программ с высокой степенью параллелизма, так как в них 3395 меньше количество промежуточных узлов, через которые должны пройти управляющие 3396 объекты при распределении между узлами кластера. Поскольку существует 3397 однозначное соответствие между резидентными процессами и узлами кластера, в 3398 тексте они используются как взаимозаменяемые термины. 3399 3400 **** Обработка выхода узлов из строя. 3401 Основным методом восстановления при выходе из строя подчиненного узла является 3402 перезапуск выполнявшихся на нем объектов на рабочих узлах (такой же метод 3403 использует язык Erlang при перезапуске подчиненных 3404 процессов\nbsp{}cite:armstrong2003thesis). Для того чтобы реализовать этот метод 3405 в рамках иерархии управляющих объектов, узел-отправитель сохраняет каждый 3406 объект, передаваемый на другие узлы кластера, и в случае отказа произвольного 3407 количества узлов, на которые были переданы объекты, их копии перераспределяются 3408 между оставшимися узлами без индивидуальной обработки программистом. Если больше 3409 не осталось узлов, на которые можно отправить объекты, то они выполняются 3410 локально. В отличие от "тяжеловесного" метода контрольных точек восстановления, 3411 используемого планировщиками задач для высокопроизводительных кластеров, 3412 древовидная иерархия узлов в паре с иерархией объектов позволяет автоматически 3413 продолжить выполнение программы при выходе из строя произвольного количества 3414 подчиненных узлов без перезапуска каких-либо процессов параллельной программы, а 3415 процесссы выполняют роль единиц выделения ресурсов. 3416 3417 Возможный подход к обработке выхода из строя главного узла (узла, на котором 3418 запускается главный управляющий объект) заключается в копировании этого главного 3419 объекта на резервный узел и синхронизации любых изменений между двумя копиями 3420 объекта посредством распределенных транзакций, однако, этот подход не 3421 соотносится с асинхронностью вычислительных ядер и слишком сложен в реализации. 3422 На практике оказывается, что главный управляющий объект обычно не выполняет 3423 операции параллельно, а последовательно переходит от вычисления одного шага 3424 программы к вычислению другого, и, значит, имеет не больше одного подчиненного в 3425 каждый момент времени. (Каждый подчиненный объект представляет собой 3426 последовательный шаг вычислений, который может быть, а может не быть 3427 параллельным внутри.) Имея это ввиду, можно упростить синхронизацию состояния 3428 главного объекта программы: отправить главный объект на подчиненный узел вместе 3429 с его подчиненным объектом. Тогда при выходе из строя главного узла, копия 3430 главного объекта принимает подчиненный объект (поскольку оба объекта находятся 3431 на одном и том же узле), и время на восстановление не тратится. Если же выходит 3432 из строя подчиненный узел, на который был отправлен подчиненный объект вместе с 3433 копией главного объекта, то подчиненный объект отправляется на один из 3434 оставшихся узлов, и в худшем случае текущий шаг вычислений выполняется заново. 3435 3436 Описанный выше подход предназначен для объектов, у которых нет родителя и 3437 которые имеют только один подчиненный объект в каждый момент времени, и 3438 повторяет механизм работы контрольных точек восстановления. Преимуществом 3439 данного подхода является то, что он 3440 - сохраняет состояние программы только между последовательными шагами вычислений 3441 (когда оно занимает минимальный объем памяти), 3442 - сохраняет только актуальное данные и 3443 - использует для сохранения состояния оперативную память другого узла кластера, 3444 а не дисковое хранилище. 3445 Этот подход позволяет выдержать выход из строя не более одного /любого/ узла 3446 кластера за один шаг вычислений или произвольного количества подчиненных узлов в 3447 любой момент работы программы. 3448 3449 Для кластера из четырех узлов и гипотетической параллельной программы алгоритм 3450 восстановления после сбоев работает следующим образом 3451 (рис.\nbsp{}[[fig-fail-over-example]]). 3452 1. /Исходное состояние./ На начальном этапе вычислительный кластер не требует 3453 никакой настройки за исключением настройки сети. Алгоритм предполагает полную 3454 связность узлов кластера и лучше всего работает с древовидными топологиями 3455 сети, в которых все узлы кластера соединены несколькими сетевыми 3456 коммутаторами. 3457 2. /Построение иерархии узлов./ При первичной загрузке на всех узлах кластера 3458 запускаются резидентные процессы, которые совместно строят свою иерархию 3459 поверх топологии сети кластера. Положение резидентного процесса в иерархии 3460 определяется позицией IP-адреса его узла в диапазоне IP-адресов сети. Для 3461 установления связи каждый из процессов соединяется только с предполагаемым 3462 главным процессом. В данном случае процесс на узле \(A\) становится главным 3463 процессом для всех остальных. Иерархия изменяется, только если новый узел 3464 присоединяется к кластеру или какой-либо из существующих узлов выходит из 3465 строя. 3466 3. /Запуск главного управляющего объекта./ Первый управляющий объект запускается 3467 на одном из подчиненных узлов (узел \(B\)). Главный объект может иметь только 3468 один подчиненный объект в любой момент времени, а резервная копия главного 3469 объекта посылается вместе с этим подчиненным объектом \(T_1\) на главный узел 3470 \(A\). \(T_1\) представляет собой последовательный шаг программы. В программе 3471 может быть произвольное количество последовательных шагов, и, когда узел 3472 \(A\) выходит из строя, текущий шаг перезапускается с начала. 3473 4. /Запуск подчиненных управляющих объектов./ Управляющие объекты \(S_1\), 3474 \(S_2\), \(S_3\) запускаются на подчиненных узлах кластера. Когда узел \(B\), 3475 \(C\) или \(D\) выходит из строя, соответствующий руководящий управляющий 3476 объект перезапускает завершившиеся некорректно подчиненные объекты (\(T_1\) 3477 перезапускает \(S_1\), главный объект перезапускает \(T_1\) и т.д.). Когда 3478 выходит из строя узел \(B\), главный объект восстанавливается из резервной 3479 копии. 3480 3481 #+name: fig-fail-over-example 3482 #+header: :headers '("\\input{preamble}\\setdefaultlanguage{russian}") 3483 #+begin_src latex :file build/fail-over-example-ru.pdf :exports results :results raw 3484 \input{tex/preamble} 3485 \newcommand*{\spbuInsertFigure}[1]{% 3486 \vspace{2\baselineskip}% 3487 \begin{minipage}[b]{0.5\linewidth}% 3488 \Large% 3489 \input{#1}% 3490 \end{minipage}% 3491 }% 3492 \noindent% 3493 \spbuInsertFigure{tex/cluster-0}~\spbuInsertFigure{tex/frame-0}\newline 3494 \spbuInsertFigure{tex/frame-3-ru}~\spbuInsertFigure{tex/frame-4-ru}\newline 3495 \spbuInsertFigure{tex/legend-ru} 3496 #+end_src 3497 3498 #+caption[Схема работы алгоритма восстановления после сбоев]: 3499 #+caption: Схема работы алгоритма восстановления после сбоев. 3500 #+name: fig-fail-over-example 3501 #+attr_latex: :width \textwidth 3502 #+RESULTS: fig-fail-over-example 3503 [[file:build/fail-over-example-ru.pdf]] 3504 3505 **** Результаты тестирования. 3506 Алгоритм обеспечения отказоустойчивости был протестирован на физическом кластере 3507 (см.\nbsp{}табл.\nbsp{}[[tab-ant]]) на примере распределенной программы для модели 3508 АР, подробно описанной в разделе\nbsp{}[[#sec-arma-mpp]]. Программа состоит из серии 3509 функций, каждая из которых применяется к результату работы предыдущей. Некоторые 3510 из функций вычисляются параллельно, так что вся программа состоит из 3511 последовательно выполняющихся шагов, некоторые из которых внутри реализованы 3512 параллельно, чтобы получить большую производительность. Только наиболее 3513 ресурсоемкий этап программы (генерация взволнованной морской поверхности) 3514 выполняется параллельно на всех узлах, другие этапы выполняются параллельно на 3515 всех процессорных ядрах главного узла. 3516 3517 Программа была переписана для распределенной версии фреймворка, что потребовало 3518 добавления методов чтения/записи для каждого управляющего объекта, который 3519 передается по сети и небольших изменений исходного кода для корректной обработки 3520 выхода из строя узла с главным объектом. Главный объект был помечен, чтобы 3521 фреймворк смог передать его на подчиненный узел вместе с подчиненным ему 3522 объектом. Другие изменения исходного кода были связаны с изменением программного 3523 интерфейса фреймворка. Таким образом, обеспечение отказоустойчивости посредством 3524 иерархии управляющих объектов, в основном, прозрачно для программиста: требует 3525 маркировки главного объекта для его репликации на резервный узел и добавления 3526 кода для чтения/записи объектов в байтовый буфер. 3527 3528 В ряде экспериментов была измерена производительность новой версии программы при 3529 выходе из строя различных типов узлов во время выполнения: 3530 - без выхода из строя узлов, 3531 - выход из строя главного узла (на котором запускается главный объект), 3532 - выход из строя подчиненного узла (на который копируется главный объект 3533 программы). 3534 Только два напрямую соединенных узла кластера были использованы в тесте. Выход 3535 из строя узла имитировался путем отправки сигнала ~SIGKILL~ резидентному 3536 процессу на соответствующем узле, сразу после того как копия главного объекта 3537 создана. Приложение сразу обнаруживало выход из строя узла, поскольку 3538 соответствующее соединение закрывалось; на практике, однако, выход узла из строя 3539 обнаруживается только по прошествии настраиваемого время ожидания протокола 3540 TCP\nbsp{}cite:rfc5482. Время выполнения этих запусков сравнивалось со временем 3541 выполнения без имитирования выхода из строя узлов. Результаты тестов 3542 представлены на рис.\nbsp{}[[fig-master-slave-failure]], а схема распределения 3543 управляющих объектов между двумя узлами на рис.\nbsp{}[[fig-master-slave-backup]]. 3544 3545 #+name: fig-master-slave-backup 3546 #+begin_src dot :exports results :file build/master-slave-backup-ru.pdf 3547 digraph { 3548 3549 node [fontname="Old Standard",fontsize=14,margin="0.055,0.055",shape=box,style=rounded] 3550 graph [nodesep="0.30",ranksep="0.30",rankdir="BT"] 3551 edge [arrowsize=0.66] 3552 3553 m1 [label="{{<master>Master | 10.0.0.1} | {<main>M | <slave1>S₁ | <slave3>S₃}}",shape=record] 3554 m2 [label="{{<slave>Slave | 10.0.0.2} | {M' | <step>N | <slave2>S₂ | <slave4>S₄}}",shape=record] 3555 3556 m2->m1 3557 3558 } 3559 #+end_src 3560 3561 #+name: fig-master-slave-backup 3562 #+caption[Отображение объектов на главный и подчиненный узлы]: 3563 #+caption: Главный и подчиненный узлы, а также отображение главного объекта 3564 #+caption: \(M\), его копии \(M'\), объекта, соответствующего текущему шагу 3565 #+caption: выполнения \(N\) и подчиненных объектов \(S_{1,2,3}\) на эти узлы. 3566 #+RESULTS: fig-master-slave-backup 3567 [[file:build/master-slave-backup-ru.pdf]] 3568 3569 Как и ожидалось, существует большая разница в производительности приложения при 3570 выходе из строя различных типов узлов. В случае отказа подчиненного узла главный 3571 управляющий объект вместе с некоторыми подчиненными (которые были распределены 3572 на подчиненный узел) теряются, но главный узел имеет копию главного объекта и 3573 использует ее, чтобы продолжить выполнение. Таким образом, при выходе из строя 3574 подчиненного узла ничего не теряется, за исключением потенциала 3575 производительности подчиненного узла. В случае выхода из строя главного узла 3576 копия главного объекта, а также подчиненный объект, который переносит эту копию, 3577 теряются, но подчиненный узел имеет оригинальный главный объект и использует его 3578 для перезапуска вычислений с текущего шага, т.е.\nbsp{}отправляет подчиненный 3579 объект на один из оставшихся узлов кластера (в случае двух напрямую соединенных 3580 узлов он отправляет объект сам себе). Таким образом, разница в 3581 производительности приложения объясняется разным количеством и разными ролями 3582 объектов, которые теряются при выходе из строя того или иного узла. 3583 3584 Обнаружение выхода из строя подчиненного узла требует некоторого времени: это 3585 происходит, только когда подчиненный объект, переносящий копию главного, 3586 заканчивает выполнение и пытается добраться до родительского объекта. 3587 Мгновенное обнаружение требует принудительной остановки подчиненного объекта, 3588 что может быть неприменимо для программ со сложной логикой. 3589 3590 #+name: fig-master-slave-failure 3591 #+begin_src R :file build/master-slave-failure-ru.pdf 3592 source(file.path("R", "benchmarks.R")) 3593 par(family="serif") 3594 data <- arma.load_master_slave_failure_data() 3595 arma.plot_master_slave_failure_data( 3596 data, 3597 list( 3598 master="Bscheduler (главный узел)", 3599 slave="Bscheduler (подчиненный узел)", 3600 nofailures="Bscheduler (без выхода из строя)" 3601 ) 3602 ) 3603 title(xlab="Размер взволнованной поверхности", ylab="Время, сек.") 3604 #+end_src 3605 3606 #+caption[Производительность программы при выходе из строя узлов]: 3607 #+caption: Производительность модели АР при выходе из строя узлов. 3608 #+name: fig-master-slave-failure 3609 #+RESULTS: fig-master-slave-failure 3610 [[file:build/master-slave-failure-ru.pdf]] 3611 3612 Итого, если выход из строя узла происходит сразу после того как копия главного 3613 объекта сделана, лишь малая часть производительности теряется в независимости от 3614 того, теряется главный объект или его копия. 3615 3616 **** Обсуждение результатов тестирования. 3617 Поскольку выход из строя имитируется, сразу после того как первый подчиненный 3618 объект достигает точки назначения (узла, на котором предполагается его 3619 выполнить), выход из строя подчиненного узла приводит к потере небольшой доли 3620 производительности; на практике, где выход из строя может произойти в середине 3621 генерации взволнованной поверхности, потери производительности при выходе из 3622 строя /резервного/ узла (узла, где находится копия главного объекта) были бы 3623 выше. Аналогично, на практике количество узлов в кластере больше, а значит 3624 меньшее количество подчиненных объектов теряется при выходе из строя главного 3625 узла, из-за чего потери производительности меньше. В тесте потери в случае 3626 выхода из строя подчиненного узла выше, что является результатом отсутствия 3627 параллелизма в начале генерации взволнованной поверхности моделью АР: первая 3628 часть вычисляется последовательно, а другие части вычисляются, только когда 3629 первая доступна. Таким образом, потеря первого подчиненного объекта замедляет 3630 выполнение всех зависимых объектов в программе. 3631 3632 Алгоритм восстановления после сбоев гарантирует обработку выхода из строя одного 3633 узла на один последовательный шаг программы; больше сбоев может быть выдержано, 3634 если они не затрагивают главный узел. Алгоритм обрабатывает одновременный выход 3635 из строя всех подчиненных узлов, однако, если главный и резервный узлы вместе 3636 выходят из строя, у программы нет ни единого шанса продолжить работу. В этом 3637 случае состояние текущего шага вычислений теряется полностью, и его можно 3638 восстановить только перезапуском программы с начала (что на данный момент не 3639 реализовано в Bscheduler). 3640 3641 Управляющие объекты являются абстракциями, отделяющими распределенное приложение 3642 от физических устройств: для непрерывной работы программы не важно, сколько 3643 узлов кластера в данный момент работают. Управляющие объекты позволяют 3644 отказаться от выделения физического резервного узла для обеспечения устойчивости 3645 к выходу из строя главного узла: в рамках иерархии управляющих объектов любой 3646 физический узел (кроме главного) может выполнять роль резервного. Наконец, 3647 иерархия управляющих объектов позволяет обрабатывать сбои прозрачно для 3648 программиста, определяя порядок действий из внутреннего состояния объекта. 3649 3650 Проведенные эксперименты показывают, что параллельной программе необходимо иметь 3651 несколько последовательных этапов выполнения, чтобы сделать ее устойчивой к 3652 сбоям узлов, иначе выход из строя резервного узла фактически вызывает 3653 восстановление исходного состояния программы. Несмотря на то что вероятность 3654 сбоя резервного узла меньше вероятности сбоя одного из подчиненных узлов, это не 3655 повод потерять все данные, когда выполнявшаяся продолжительное время программа 3656 почти завершилась. В общем случае, чем больше последовательных этапов вычислений 3657 содержит параллельная программа, тем меньше времени потеряется в случае сбоя 3658 резервного узла, и, аналогично, чем больше параллельных частей содержит каждый 3659 последовательный этап, тем меньше времени потеряется при сбое главного или 3660 подчиненного узла. Другими словами, чем больше количество узлов, на которое 3661 масштабируется программа, тем она становится более устойчива к их сбоям. 3662 3663 Хотя это не было показано в экспериментах, Bscheduler не только обеспечивает 3664 устойчивость к выходу из строя узлов кластера, но и позволяет автоматически 3665 вводить новые узлы в кластер и распределять на них часть управляющих объектов из 3666 уже запущенных программ. В контексте фреймворка этот процесс тривиален, 3667 поскольку не требует перезапуска незавершившихся управляющих объектов и 3668 копирования их состояния, и не изучался экспериментально в данной работе. 3669 3670 Теоретически, основанная на иерархии отказоустойчивость может быть реализована 3671 поверх библиотеки передачи сообщений без потери общности. Хотя использование 3672 незагруженных узлов вместо вышедших из строя в рамках такой библиотеки 3673 представляет определенную сложность, поскольку количество узлов, на которых 3674 запущена программа, в таких библиотеках фиксировано, выделение достаточно 3675 большого количества узлов для программы будет достаточно для обеспечения ее 3676 отказоустойчивости. В то же время, реализация основанной на иерархии 3677 отказоустойчивости внутри самой библиотеки передачи сообщений не практично, 3678 поскольку это потребует сохранения текущего состояния параллельной программы, 3679 объем которого эквивалентен всей занимаемой ей памятью на каждом узле кластера, 3680 что, в свою очередь, не позволит сделать такой подход эффективнее контрольных 3681 точек восстановления. 3682 3683 Слабым местом описанного метода является период времени, начиная с отказа 3684 главного узла и заканчивая обнаружением сбоя подчиненным узлом, восстановлением 3685 главного объекта из копии и получением нового подчиненного объекта вместе с 3686 копией его родителя подчиненным узлом. Если на протяжении этого промежутка 3687 резервный узел выходит из строя, то состояние выполнения программы полностью 3688 теряется без возможности его восстановить, кроме как перезапуском с самого 3689 начала. Протяженность этого опасного промежутка времени может быть 3690 минимизирована, но полностью исключить вероятность внезапного завершения 3691 программы невозможно. Этот результат согласуется с исследованиями /теории 3692 невыполнимости/ в рамках которой доказывается невозможность распределенного 3693 консенсуса с хотя бы одним процессом, дающим 3694 сбой\nbsp{}cite:fischer1985impossibility и невозможность надежной передачи 3695 данных в случае сбоя одного из узлов\nbsp{}cite:fekete1993impossibility. 3696 3697 *** Выводы 3698 Современный подход к разработке и запуску параллельных программ на кластере 3699 заключается в использовании библиотеки передачи сообщений и планировщика задач, 3700 и, несмотря на то что этот подход имеет высокую эффективность с точки зрения 3701 параллельных вычислений, он недостаточно гибок, чтобы вместить в себя 3702 динамическую балансировку нагрузки и автоматическое обеспечение 3703 отказоустойчивости. Программы, написанные с помощью библиотеки передачи 3704 сообщений обычно предполагают 3705 - равномерную загрузку каждого процессора, 3706 - бесперебойное и надежное выполнение пакетных задач, и 3707 - постоянное число параллельных процессов во время выполнения, равное общему 3708 количеству процессоров. 3709 Первое предположение несправедливо для программы моделирование морского 3710 волнения, поскольку модель АР требует динамической балансировки нагрузки между 3711 процессорами для генерации каждой части поверхности только когда генерация всех 3712 зависимых частей уже закончена. Последнее предположение также несправедливо, 3713 поскольку в угоду эффективности каждая часть записывается в файл отдельным 3714 потоком асинхронно. Оставшееся предположение относится не к самой программе, а к 3715 планировщику задач, и несправедливо для больших вычислительных кластеров, в 3716 которых узлы часто выходят из строя, а планировщик восстанавливает задачу из 3717 контрольной точки, сильно увеличивая время ее выполнения. Идея предлагаемого 3718 подхода\nbsp{}--- дать параллельным программам больше гибкости: 3719 - предоставить динамическую балансировку нагрузки путем выполнения 3720 последовательных, параллельных изнутри шагов программы в режиме конвейера, 3721 - перезапускать только затронутые выходом из строя узла процессы, и 3722 - выполнять программу на как можно большем количестве узлов, которое доступно в 3723 кластере. 3724 В данном разделе обсуждаются преимущества и недостатки этого подхода. 3725 3726 В сравнении с системами пакетной обработки заданий (PBS) для распределения 3727 нагрузки на узлы кластера предлагаемый подход использует легковесные управляющие 3728 объекты вместо тяжеловесных параллельных задач. Во-первых, это позволяет иметь 3729 очереди объектов на каждом узле, вместо того чтобы иметь несколько очередей 3730 задач на весь кластер. Зернистость управляющих объектов гораздо выше, чем у 3731 пакетных задач, и, несмотря на то что время их выполнения не может быть надежно 3732 спрогнозировано (также как и время выполнения пакетных 3733 задач\nbsp{}cite:zotkin1999job), объекты из нескольких параллельных программ 3734 могут быть динамически распределены между одним и тем же множеством узлов 3735 кластера, делая нагрузку более равномерной. Недостатком является необходимость в 3736 большем количестве оперативной памяти для выполнения нескольких задач на одних и 3737 тех же узлах, а также в том что выполнение каждой программы может занять больше 3738 времени из-за общих очередей управляющих объектов. Во-вторых, предлагаемый 3739 подход использует динамическое распределение ролей главного и подчиненного между 3740 узлами кластера вместо их статического присвоения конкретным физическим узлам. 3741 Это делает узлы взаимозаменяемыми, что необходимо для обеспечения 3742 отказоустойчивости. Одновременное выполнение нескольких параллельных программ на 3743 одном и том же множестве узлов увеличивает пропускную способность кластера, но 3744 также уменьшает их производительность, взятую по отдельности; динамическое 3745 распределение ролей является основанием, на котором строится устойчивость к 3746 сбоям. 3747 3748 В сравнении с MPI для разбиения программы на отдельные сущности предлагаемый 3749 подход использует легковесные управляющие объекты вместо тяжеловесных процессов. 3750 Во-первых, это позволяет определить число обрабатываемых параллельно сущностей, 3751 исходя из задачи, а не архитектуры компьютера или кластера. Это поощряет 3752 программиста создавать столько объектов, сколько необходимо, руководствуясь 3753 алгоритмом или ограничениями на размер структур данных из предметной области 3754 задачи. В программе моделирования морского волнения минимальный размер каждой 3755 части поверхности зависит от числа коэффициентов вдоль каждой из осей, и, в то 3756 же время, количество частей должно быть больше, чем количество процессоров, для 3757 того чтобы сделать нагрузку на процессоры более равномерной. Учитывая эти 3758 ограничения оптимальный размер части определяется во время выполнения, и, в 3759 общем случае, не совпадает с количеством параллельных процессов. Недостатком 3760 является то, что, чем больше управляющих объектов в программе, тем больше общих 3761 структур данных копируется на один и тот же узел вместе с подчиненными 3762 объектами; проблема решается введением промежуточного слоя объектов, что в свою 3763 очередь влечет увеличивает сложность программы. Во-вторых, иерархия управляющих 3764 объектов совместно с иерархией узлов позволяет автоматически пересчитывать 3765 завершившиеся некорректно подчиненные объекты на выживших узлах кластера в 3766 случае выхода из строя оборудования. Это возможно, поскольку ход выполнения 3767 программы сохраняется в объектах, а не в глобальных переменных, как это делается 3768 в программах MPI. Дублируя состояние на подчиненные узлы, система пересчитывает 3769 только объекты из поврежденных процессов, а не программу целиком. Переход от 3770 процессов к управляющим объектам увеличивает производительность параллельной 3771 программы посредством динамической балансировки нагрузки, но также может 3772 повлиять на ее масштабируемость на большое количество узлов из-за дублирования 3773 состояния хода выполнения. 3774 3775 Разбиение программы на отдельные сущности применяется во фреймворке 3776 Charm++\nbsp{}cite:kale2012charm и в модели 3777 актеров\nbsp{}cite:hewitt1973universal,agha1985actors, однако ни один из 3778 подходов не использует иерархические связи для перезапуска обработки сущностей 3779 после ошибки. Вместо использования древовидной иерархии управляющих объектов эти 3780 подходы позволяют обмениваться сообщениями любой паре объектов. Такая 3781 беспорядочная схема обмена сообщениями не позволяет решить, какой объект 3782 ответственен за перезапуск другого, завершившегося ошибкой, из-за чего 3783 неуниверсальные подходы обеспечения отказоустойчивости используются вместо 3784 этого\nbsp{}cite:lifflander2014scalable. 3785 3786 Три составляющих предлагаемого подхода\nbsp{}--- управляющие объекты, конвейеры 3787 и иерархии\nbsp{}--- дополняют друг друга. Если бы управляющие объекты не 3788 содержали в себе состояние хода выполнения программы, то было бы невозможно 3789 пересчитать завершившиеся некорректно подчиненные объекты и обеспечить 3790 отказоустойчивость. Если бы иерархии узлов не было, то было бы невозможно 3791 распределить нагрузку между узлами кластера, поскольку все узлы одинаковы без 3792 иерархии. Если бы для каждого устройства не было конвейера, то было бы 3793 невозможно обрабатывать управляющие объекты асинхронно и реализовать 3794 динамическую балансировку нагрузки. Эти три сущности образуют замкнутую систему, 3795 в которой логика программы реализуется либо в управляющих объектах, либо в 3796 конвейерах, логика восстановления после сбоев в иерархии объектов, а логика 3797 передачи данных в иерархии узлов кластера. 3798 3799 ** Реализация для систем с распределенной памятью (MPP) 3800 :PROPERTIES: 3801 :CUSTOM_ID: sec-arma-mpp 3802 :END: 3803 3804 **** Распределенный алгоритм для модели АР. 3805 :PROPERTIES: 3806 :CUSTOM_ID: sec-distributed 3807 :END: 3808 3809 Этот алгоритм в отличие от параллельной версии, использует копирование данных 3810 для выполнения вычислений на других узлах кластера, и, поскольку пропускная 3811 способность сети гораздо меньше, чем у памяти, размер передаваемых по сети 3812 данных должен быть оптимизирован для получения большей производительности, чем 3813 на системе с общей памятью. Один из способов добиться этого\nbsp{}--- это 3814 распределить части взволнованной поверхности между узлами кластера, копируя на 3815 узлы коэффициенты и необходимые точки на границах, и, копируя обратно 3816 сгенерированную часть взволнованной поверхности. Авторегрессионные зависимости 3817 не позволяют создать все части сразу и статически распределить их между узлами 3818 кластера, поэтому части создаются динамически на первом узле, когда точки, от 3819 которых они зависят, становятся доступны. Таким образом, распределенный алгоритм 3820 для модели АР является алгоритмом типа ведущий-ведомый, в котором ведущий 3821 динамически создает задачи для каждой части взволнованной поверхности, принимая 3822 во внимание авторегрессионные зависимости между точками, и отправляет их на 3823 подчиненные узлы, а ведомые вычисляют каждую часть взволнованной поверхности и 3824 отправляют ее обратно ведущему. 3825 3826 В реализации для систем с распределенной памятью каждая задача моделируется 3827 управляющим объектом: существует руководящий объект, создающий подчиненные при 3828 необходимости, и подчиненные объекты, генерирующие части взволнованной 3829 поверхности. В методе ~act~ главного объекта создается подчиненный объект для 3830 первой части\nbsp{}--- части, которая не зависит ни от каких других точек. Когда 3831 этот объект возвращается, руководящий объект в методе ~react~ определяет, какие 3832 части могут быть вычислены сейчас, создает подчиненный объект для каждой части и 3833 отправляет его на конвейер. В методе ~act~ подчиненного объекта генерируется 3834 часть взволнованной поверхности, а затем объект отправляет самого себя 3835 руководителю. Метод ~react~ подчиненного объекта пустой. 3836 3837 Реализация распределенного алгоритма АР имеет ряд преимуществ по сравнению с 3838 параллельной. 3839 - Конвейеры автоматически распределяют подчиненные объекты между доступными 3840 узлами кластера, и приложение не имеет дела с такими низкоуровневыми деталями. 3841 - Нет необходимости реализовать минималистичный планировщик задач, который 3842 определяет последовательность выполнения задач (объектов), учитывая 3843 авторегрессионные зависимости: порядок выполнения полностью определяется в 3844 методе ~react~ руководящего объекта. 3845 - Нет необходимости в отдельной версии программы для машины с общей памятью: 3846 реализация работает прозрачно на любой количестве узлов, даже если планировщик 3847 задач не запущен. 3848 3849 **** Производительность реализации распределенного алгоритма для АР модели. 3850 Реализация распределенной модели АР была протестирована на двух узлах системы 3851 Ant (таб.\nbsp{}[[tab-ant]]). Для увеличения пропускной способности сети, два узла 3852 были соединены друг с другом напрямую, а максимальный размер единицы 3853 передаваемых данных (MTU) установлен в 9200 байт. Рассматривались два случая: с 3854 одним резидентным процессом Bscheduler, запущенным на локальном узле, и с двумя 3855 резидентными процессами, запущенными на каждом узле. Производительность программы 3856 сравнивалась с производительностью версии на OpenMP, запускаемой на одном узле. 3857 3858 Bscheduler превосходит OpenMP в обоих случаях 3859 (рис.\nbsp{}[[fig-bscheduler-performance]]). В случае одного узла более высокая 3860 производительность объясняется тем, что Bscheduler не сканирует очередь задач в 3861 поисках частей, для которых готовы все зависимости (как это происходит в 3862 параллельной версии алгоритма), вместо этого он для каждой части обновляет 3863 счетчик готовых частей, от которых она зависит. Такой же подход может быть 3864 использован для версии OpenMP, однако он был обнаружен только в более новой 3865 версии программы для Bscheduler, поскольку сканирование очереди задач не может 3866 быть эффективно реализовано в рамках этого планировщика. В случае двух узлов 3867 более высокая производительность объясняется большим суммарным количеством 3868 процессорных ядер (16) и высокой пропускной способностью прямого сетевого 3869 соединения. Таким образом, реализация распределенного алгоритма модели АР на 3870 Bscheduler быстрее на системе с общей памятью ввиду более эффективной обработки 3871 авторегрессионных зависимостей, а на системе с распределенной памятью его 3872 производительность масштабируется на большее количество ядер ввиду низких 3873 накладных расходов на передачу данных по прямому сетевому соединению. 3874 3875 #+name: fig-bscheduler-performance 3876 #+begin_src R :file build/bscheduler-performance-ru.pdf 3877 source(file.path("R", "benchmarks.R")) 3878 par(family="serif") 3879 data <- arma.load_bscheduler_performance_data() 3880 arma.plot_bscheduler_performance_data( 3881 data, 3882 list( 3883 openmp="OpenMP", 3884 bsc1="Bscheduler (один узел)", 3885 bsc2="Bscheduler (два узла)" 3886 ) 3887 ) 3888 title(xlab="Размер взволнованной поверхности", ylab="Время, сек.") 3889 #+end_src 3890 3891 #+name: fig-bscheduler-performance 3892 #+caption[Производительность Bscheduler и OpenMP версий программы]: 3893 #+caption: Сравнение производительности Bscheduler и OpenMP версий программы 3894 #+caption: для модели АР. 3895 #+RESULTS: fig-bscheduler-performance 3896 [[file:build/bscheduler-performance-ru.pdf]] 3897 3898 * Заключение 3899 В изучении возможностей математического аппарата для имитационного моделирования 3900 морского волнения, выходящего за рамки линейной теории волн, были достигнуты 3901 следующие основные результаты. 3902 - Процессы АР и СС были использованы для моделирования морских волн произвольных 3903 амплитуд. Интегральные характеристики взволнованной поверхности были 3904 верифицированы путем сопоставления с характеристиками реальной морской 3905 поверхности. 3906 - Новый метод был использован для вычисления поля потенциала скорости под 3907 генерируемой поверхностью. Получившееся поле было верифицировано путем 3908 сравнения с полем, вычисляемым по формулам из линейной теории волн. Новый 3909 метод эффективен с вычислительной точки зрения, поскольку все интегралы в его 3910 формуле записываются как преобразования Фурье, для которого существуют 3911 высокопроизводительные реализации. 3912 - Модель и метод были реализованы для систем как с общей, так и с распределенной 3913 памятью, и в нескольких тестах показали масштабируемость на различном 3914 количество ядер, которая близка к линейной. Модель АР более эффективна с 3915 вычислительной точки зрения на центральном процессоре, нежели на видеокарте, и 3916 превосходит по производительности модель ЛХ. 3917 3918 Одной из тем дальнейших исследований является изучение возможности генерации 3919 волн произвольных профилей на базе смешанного процесса АРСС. Другим направлением 3920 является интеграция разработанной модели и метода расчета давлений в 3921 существующие пакеты прикладного программного обеспечения. 3922 3923 * Выводы 3924 Задача вычисления давлений под реальной морской поверхностью решается явно, 3925 минуя предположения линейной теории волн и теории волн малой амплитуды. Это 3926 решение в паре с моделями АР и СС, генерирующими морские волны произвольных 3927 амплитуд, может быть использовано для расчета влияния колебаний волн на 3928 поведение динамического объекта в открытом море, и в случае волн больших 3929 амплитуд дает более правдоподобное поле потенциала скорости, чем решения, 3930 полученные в рамках линейной теории волн и теории волн малой амплитуды. 3931 3932 Численные эксперименты показывают, что как генерация взволнованной поверхности 3933 так и расчет гидродинамических давлений эффективно реализуются как на системах с 3934 общей, так и с распределенной памятью, без проведения вычислений на видеокарте. 3935 Высокая производительность в случае генерации взволнованной поверхности 3936 обеспечивается использованием специализированного планировщика задач и 3937 библиотеки для работы с многомерными массивами, а в случае вычисления поля 3938 потенциала скорости\nbsp{}--- использованием алгоритмов БПФ для вычисления 3939 интегралов. 3940 3941 Разработанный в работе математический аппарат и его численная реализация могут 3942 стать основой виртуального полигона, предназначенного для расчетов динамики 3943 морских объектов. Использование новых моделей в виртуальном полигоне позволит 3944 - проводить более длительные сеансы моделирования без потери эффективности ввиду 3945 отсутствия периодичности в моделях, 3946 - получить более правдоподобные поля давлений благодаря использованию нового 3947 метода вычисления поля потенциала скорости, и 3948 - сделать программный комплекс более надежным благодаря быстрой сходимости 3949 моделей и использованию отказоустойчивого планировщика пакетных задач. 3950 3951 * Благодарности 3952 Графики в этой работе были подготовлены с помощью языка для статистических 3953 вычислений R\nbsp{}cite:rlang2016,sarkar2008lattice и программного обеспечения 3954 Graphviz\nbsp{}cite:gansner00anopen. Документ был подготовлен с использованием 3955 Org-mode\nbsp{}cite:schulte2011org2,schulte2011org1,dominik2010org для GNU 3956 Emacs, предоставляющего вычислительное окружение для воспроизводимых 3957 исследований. Это означает, что все графики можно воспроизвести и 3958 соответствующие утверждения проверить на других вычислительных системах, 3959 скопировав репозиторий диссертации[fn:repo], установив Emacs и экспортировав 3960 документ. 3961 3962 [fn:repo] [[https://github.com/igankevich/arma-thesis]] 3963 3964 * Список сокращений и условных обозначений 3965 - <<<MPP>>> :: Massively Parallel Processing, класс вычислительных систем с 3966 разделенной памятью. 3967 - <<<SMP>>> :: Symmetric Multi-Processing, класс вычислительных систем с общей 3968 памятью. 3969 - <<<GPGPU>>> :: General-purpose computing on graphics processing units, 3970 вычисления общего назначения на графических ускорителях. 3971 - АКФ :: автоковариационная функция. 3972 - <<<БПФ>>> :: быстрое преобразование Фурье. 3973 - <<<ГПСЧ>>> :: генератор псевдослучайных чисел. 3974 - <<<ГУ>>> :: граничное условие. 3975 - <<<ДУЧП>>> :: дифференциальное уравнение в частных производных. 3976 - <<<НБП>>> :: нелинейное безынерционное преобразование, позволяющее задать 3977 произвольный закон распределения аппликат взволнованной 3978 поверхности без изменения ее исходной автоковариационной функции. 3979 - АР :: процесс авторегрессии. 3980 - АРСС :: процесс авторегрессии скользящего среднего. 3981 - СС :: процесс скользящего среднего. 3982 - ЛХ :: модель Лонге---Хиггинса, формула которой выводится в рамках линейной 3983 теории волн. 3984 - <<<LAMP>>> :: Large Amplitude Motion Programme, программа для моделирования 3985 качки судна на морских волнах. 3986 - <<<ЦПТ>>> :: центральная предельная теорема. 3987 - <<<ПМ>>> :: аппроксимация Пирсона---Московица для спектра морского волнения. 3988 - <<<ЮУ>>> :: уравнения Юла---Уокера, используемые для определения коэффициентов 3989 авторегрессионной модели по заданной автоковариационной функции. 3990 - <<<МНК>>> :: метод наименьших квадратов. 3991 - <<<ФПР>>> :: функция плотности распределения. 3992 - <<<ФР>>> :: функция распределения. 3993 - <<<BSP>>> :: Bulk Synchronous Parallel. 3994 - OpenCL :: Open Computing Language, технология параллельного программирования 3995 для гибридных систем с видеокартами или другими сопроцессорами. 3996 - <<<OpenGL>>> :: Open Graphics Library. 3997 - OpenMP :: Open Multi-Processing, технология параллельного программирования 3998 для многопроцессорных систем. 3999 - <<<MPI>>> :: Message Passing Interface. 4000 - <<<FMA>>> :: Fused multiply-add. 4001 - <<<DCMT>>> :: Dynamic creation of Mersenne Twisters, алгоритм создания 4002 генераторов псевдослучайных чисел, которые выдают 4003 некоррелированные последовательности, будучи запущенными 4004 параллельно. 4005 - <<<GSL>>> :: GNU Scientific Library. 4006 - <<<BLAS>>> :: Basic Linear Algebra Sub-programmes. 4007 - <<<LAPACK>>> :: Linear Algebra Package. 4008 - <<<DNS>>> :: Dynamic name resolution. 4009 - РГШ :: Распределение на основе ряда Грама---Шарлье. 4010 - АНР :: Асимметричное нормальное распределение. 4011 - <<<PBS>>> :: Portable batch system, система выделения и распределения ресурсов 4012 кластера под параллельные программы. 4013 - Трансцендентные функции :: математические функции, не являющиеся 4014 алгебраическими (т.е.\nbsp{}логарифмические, тригонометрические, 4015 экспоненциальные и др.). 4016 4017 #+begin_export latex 4018 \input{postamble} 4019 #+end_export 4020 4021 bibliographystyle:ugost2008 4022 bibliography:bib/refs.bib 4023 4024 * Приложение 4025 ** Вывод формулы модели Лонге---Хиггинса 4026 :PROPERTIES: 4027 :CUSTOM_ID: longuet-higgins-derivation 4028 :END: 4029 4030 Двухмерная система уравнений\nbsp{}eqref:eq-problem в рамках линейной теории 4031 волн записывается как 4032 \begin{align*} 4033 & \phi_{xx} + \phi_{zz} = 0,\\ 4034 & \zeta(x,t) = -\frac{1}{g} \phi_t, & \text{на }z=\zeta(x,t), 4035 \end{align*} 4036 где \(\frac{p}{\rho}\) включено в \(\phi_t\). Решение уравнения Лапласа ищется в 4037 виде ряда Фурье cite:kochin1966theoretical: 4038 \begin{equation*} 4039 \phi(x,z,t) = \int\limits_{0}^{\infty} e^{k z} 4040 \left[ A(k, t) \cos(k x) + B(k, t) \sin(k x) \right] dk. 4041 \end{equation*} 4042 Подставляя его в граничное условие, получаем 4043 \begin{align*} 4044 \zeta(x,t) &= -\frac{1}{g} \int\limits_{0}^{\infty} 4045 \left[ A_t(k, t) \cos(k x) + B_t(k, t) \sin(k x) \right] dk \\ 4046 &= -\frac{1}{g} \int\limits_{0}^{\infty} C_t(k, t) \cos(kx + \epsilon(k, t)). 4047 \end{align*} 4048 Здесь \(\epsilon\)\nbsp{}--- белый шум, а \(C_t\) включает в себя значение 4049 \(dk\). Подставляя бесконечную сумму вместо интеграла, получаем двухмерную 4050 форму\nbsp{}[[eqref:eq-longuet-higgins]]. 4051 4052 ** Производная в направлении нормали к поверхности 4053 :PROPERTIES: 4054 :CUSTOM_ID: directional-derivative 4055 :END: 4056 4057 Производная от \(\phi\) в направлении вектора \(\vec{n}\) определяется как 4058 \(\nabla_n\phi=\nabla\phi\cdot\frac{\vec{n}}{|\vec{n}|}\). Вектор \(\vec{n}\), 4059 направленный по нормали к поверхности \(z=\zeta(x,y)\) в точке \((x_0,y_0)\) 4060 определяется как 4061 \begin{equation*} 4062 \vec{n} = \begin{bmatrix}\zeta_x(x_0,y_0)\\\zeta_y(x_0,y_0)\\-1\end{bmatrix}. 4063 \end{equation*} 4064 Отсюда производная в направлении нормали к поверхности определяется 4065 \begin{equation*} 4066 \nabla_n \phi = \phi_x \frac{\zeta_x}{\sqrt{\zeta_x^2+\zeta_y^2+1}} 4067 + \phi_y \frac{\zeta_y}{\sqrt{\zeta_x^2+\zeta_y^2+1}} 4068 + \phi_z \frac{-1}{\sqrt{\zeta_x^2+\zeta_y^2+1}}, 4069 \end{equation*} 4070 где производные \(\zeta\) вычисляются в \((x_0,y_0)\).