CC BY
17
УДК 639.3.06
А. А. Недоступ, А. О. Ражев
ОПТИМИЗАЦИЯ РАСЧЕТА СИЛ ГИДРОДИНАМИЧЕСКОГО СОПРОТИВЛЕНИЯ ПРИ КОМПЬЮТЕРНОМ МОДЕЛИРОВАНИИ САДКОВ АКВАКУЛЬТУРЫ1
Предложенный алгоритм расчета сил гидродинамического сопротивления садка аквакульту-ры, оптимизированный для работы на современных персональных компьютерах, основан на разложении зависимостей сил от скорости потока в ряд Тейлора 4-го порядка. Вычисление производилось с использованием SIMD-инструкций одновременно для двух гидродинамических коэффициентов: гидродинамического коэффициента для сети, расположенной параллельно потоку воды, и гидродинамического коэффициента для сети, расположенной перпендикулярно потоку воды. При помощи компьютерной программы был проведен численный эксперимент и определены скоростные характеристики оптимизированного алгоритма в сравнении с алгоритмом, использующим операцию возведения в степень. Замер времени выполнения производился для одного миллиарда узлов. Скорость работы оптимизированного алгоритма оказалась в 35 раз выше. Аппроксимация результата требуется только при скорости перемещения узла модели точечных масс ниже 0,2 м/с.
Ключевые слова: гидродинамическое сопротивление, дискретная математическая модель, компьютерное моделирование, SIMD, ряд Тейлора.
Введение
Имитационное моделирование, в отличие от других видов исследований, таких как натурный эксперимент или физическое моделирование, является наиболее гибким и эффективным с точки зрения временных и материальных затрат [1, 2]. Благодаря постоянному совершенствованию элементной базы современных персональных компьютеров и, одновременно, их удешевлению, последние становятся все более доступными широкому кругу пользователей. В настоящее время персональные компьютеры имеют мощные вычислительные ресурсы, способные эффективно обрабатывать большие объемы данных в реальном времени. Это достигается возможностью одновременной обработки множества данных на одном или нескольких вычислительных модулях, таких как центральный и графический процессоры.
Садки аквакультуры состоят из множества однотипных элементов (рис. 1). Одним из методов моделирования таких объектов является дискретное моделирование. При данном виде моделирования, с учетом строения садка, его математическая модель условно разбивается на множество дискретных элементов, способных обрабатываться по одному алгоритму.
I, н
Рис. 1. Садок аквакультуры при трехмерной постановке задачи
При дискретном моделировании динамики поведения садка под влиянием различных факторов внешней среды (приливы/отливы, неравномерное течение, волнение) важную роль
1 Статья подготовлена в рамках выполнения гранта РФФИ №15-08-464а.
играет эффективность алгоритмов. Одним из трудоемких этапов расчета является вычисление мгновенных значений и направлений сил гидродинамического сопротивления, приложенных к каждому элементу садка.
В связи с вышеизложенным предлагается метод оптимизации вычисления сил гидродинамического сопротивления, основанный на разложении функции компонентов их вектора от скорости (при трехмерной постановке задачи для каждой силы получаем три функции) в ряд Тейлора и применении для вычислений SIMD-инструкций центрального и графического процессоров.
Постановка задачи
Учитывая то, что элементы садка представляют собой сетные участки, проведем расчет сил гидродинамического сопротивления по аналогии с сетью [3]. В модели точечных масс [4] проекции сил гидродинамического сопротивления имеют следующие эвристические зависимости от проекций участка сети:
д=1
X 2
( п—2Л
с -с-Г*
V 90 0)
1 Г Г 1 Г Г , Яу = -С90 v2pSFв, Я2 = 1С90 ^2pSFo, (1)
2 Г 2 г
где Ях, Яу, Я2 - проекции сил гидродинамического сопротивления; с0 - гидродинамический коэффициент для сети, расположенной параллельно потоку воды; с90 - гидродинамический коэффициент для сети, расположенной перпендикулярно потоку воды; г - длина связи между узлами модели точечных масс; гх, гу, г2 - длина проекции связи на оси прямоугольной декартовой системы координат, ось ОХ которой направлена по направлению вектора скорости перемещения узла относительно воды; V - абсолютное значение скорости перемещения узла относительно воды; р - плотность воды; S - площадь участка сети; Fo - сплошность сети.
Гидродинамические коэффициенты находятся по следующим формулам:
Со = ^ Яе-6, С90 = 16 , Ке = ^ , (2)
где Яе - число Рейнольдса для сети; й - диаметр ниток сети; 0 - кинематическая вязкость воды; Ь - коэффициент, учитывающий тип материала сети (Ь = 0,165 для капроновой сети, Ь = 0,16 для полиэтиленовой сети).
Для вычисления гидродинамических коэффициентов в формулах (2) используется математическая операция возведения в степень. Данная операция достаточно трудоемка при вычислениях на персональном компьютере.
Решение
Рассмотрим метод оптимизации вычисления сил гидродинамического сопротивления, основанный на разложении функции компонентов их вектора от скорости (при трехмерной постановке задачи для каждой силы получаем три функции) в ряд Тейлора [5] и применении SIMD-инструкций [6] центрального и графического процессоров.
Подставив (2) в (1), получим
я=V2-6 - V2-6) х, (3)
где Я - одна из компонент сил гидродинамического сопротивления; к0, к - постоянные для участка сетной части садка коэффициенты, зависящие от параметров посадки и материала участка; х - безразмерный параметр, определяющий ориентацию участка относительно вектора относительной скорости V (скорости движения воды относительно участка садка).
Значения коэффициента к90 одинаковы, а значения коэффициента к0 различны для каждого компонента Я. В декартовой прямоугольной системе координат при трехмерной постановке задачи для компонент ЯОг и Ят, исходя из (1), получаем к0 = 0. Данные коэффициенты постоянны в течение всего процесса моделирования и должны вычисляться перед началом работы основного цикла алгоритма.
Параметр х является переменным и вычисляется на каждой итерации алгоритма моделирования для каждой компоненты Я.
В выражении (3) математические операции возведения переменной v в степень запишем
в виде
100 = V, (4)
где p = 2 - Ь для ^ иp = 1,72 для k90.
Разложим выражение (4) в ряд Тейлора четвертого порядка [5]. Общий вид разложения задается выражением
г (v ) = ! ^ (v -«У.
1=0 1 !
где/(г)(а) - значение производной 1-го порядка разлагаемой функции в точке a. При a = 1 для разложения четвертого порядка получим
f (v) = 1 + p(v - 1>fl +1 (p -l)(v - 1>fl + ^(p - 2)(v -1)1 1 + -(p -3)(v -1) I I (5)
v 2 v
Ь _w J, 1
3V 4 4
В выражении (5) используются только операции умножения, сложения и вычитания. Данные операции выполняются на современных процессорах аппаратным способом и быстро.
Реализуем алгоритм вычисления (5) для персонального компьютера с одним или несколькими центральным и графическим процессорами с поддержкой SIMD-инструкций. В современном центральном процессоре архитектуры x86 и большинстве графических процессоров одна SIMD-инструкция способна одновременно обрабатывать четыре вещественных значения одинарной точности. Учитывая это и то, что в выражении присутствуют четыре группы однотипных последовательностей операций, попробуем оптимизировать данный алгоритм. При оптимизации необходимо учесть тот факт, что для вычисления трех компонент вектора силы гидродинамического сопротивления необходимо вычислить (5) дважды - для k0 и k90.
Для анализа описанного алгоритма и определения его основных характеристик (области определения, скорости работы, погрешности вычислений) была разработана компьютерная программа. В программе был реализован предлагаемый алгоритм с использованием вычислительных ресурсов центрального процессора с набором инструкций SSE3. Ниже представлен программный код алгоритма на языках программирования С++ и Ассемблер.
float fastcall Func(float x, float y) { _
static float one = l.f; static float& vals = Values[0]; float v;
asm {
mov eax, vals
movss xmml, x
subss xmml, one
movsldup xmm0, xmml
paddd xmml, [eax+32]
movlhps xmm0, xmm0 //x,x,x,x
movlhps xmml, xmml //x,1,x,1
mulps xmml, xmml //xA2,l,xA2,l
mulps xmml, xmm0 //xA3,x,xA3,x
mulps xmm0, xmml //xA4,xA2,xA4,xA2
mulps xmml, [eax+0]
mulps xmm0, [eax+l6]
addps xmm0, xmml
haddps xmm0, xmm0
addss xmm0, one
movss v, xmm0
}
return v;
}
Эксперимент
При помощи компьютерной программы был проведен численный эксперимент и определены скоростные характеристики оптимизированного алгоритма в сравнении с алгоритмом, использующим операцию возведения в степень. Замер времени выполнения производился для одного миллиарда узлов. В эксперименте использовался персональный компьютер с характеристиками, указанными в табл. 1. Результат замера показан в табл. 2.
Таблица 1
Основные характеристики экспериментального оборудования
Центральный процессор AMD Phenom II
Количество ядер 4
Тактовая частота 3.2 ГГц
Оперативная память 12 Гб
Таблица 2
Параметры численного эксперимента
Параметр Алго1 ритм
Возведение в степень Оптимизированный
Количество узлов 1000000000
Длительность вычислений 5 минут 41,2 секунды 9,7 секунды
Максимальная относительная погрешность при у = [0,08 ... 3,7] м/с - 5 %
По данным табл. 2 можно сделать вывод, что скорость работы оптимизированного алгоритма оказалась в 35 раз выше.
Графики зависимостей Ау) и относительной погрешности вычисления е(у) от относительной скорости V приведены на рис. 2. На графиках сплошной линией изображена зависимость, полученная алгоритмом возведения в степень, пунктирной линией - зависимость, полученная оптимизированным алгоритмом, а штрихпунктирной линией - зависимость относительной погрешности вычислений от V для оптимизированного алгоритма.
0,7 0,6 0,5 0,4 0,3 0,2 0,1 0,0 -0,1 -0,2 -0,3 -0,4
0,00 0,04 0,03 0,12 0,16 0,20 0.24 0,23 0,32 0,36 0,40 0,44
V, и/с
Рис. 2. Графики зависимостей^у) и е(у) при низких значениях скорости
Обсуждение
На рис. 3 видно, что в интервале значений скорости от 0,2 до 1,8 м/с погрешность минимальная и практически нулевая. В интервале значений скорости выше 1,8 до 3,7 м/с погрешность не превышает 5 %. В интервале значений скорости до 0,2 м/с при уменьшении скорости погрешность начинает увеличиваться. Более того, при значениях скорости менее 0,08 м/с (см. рис. 1) у оптимизированного алгоритма результат вычисления Ду) становится отрицательным. Данное обстоятельство следует учитывать в тех случаях, когда необходимо производить исследования при низких или нулевых значениях скорости.
Рис. 3. Графики зависимостейf(v) и e(v) при значениях скорости до 3,7 м/с
Одним из методов устранения указанного недостатка алгоритма является аппроксимация результата при значениях скорости ниже 0,2 м/с. Возможна аппроксимация линейным сплайном как между точками (0, 0) и (0,2, f(0,2)), так и между точками (0, f(0,2)) и (0,2, f(0,2)). Оба вида должны дать приемлемый результат.
Заключение
Сравнение скоростных характеристик оптимизированного алгоритма и результатов, полученных с применением алгоритма, использующего операцию возведения в степень, показало, что скорость работы оптимизированного алгоритма оказалась в 35 раз выше (замер времени выполнения производился для одного миллиарда узлов). Аппроксимация результата требуется только при скорости перемещения узла модели точечных масс ниже 0,2 м/с.
СПИСОК ЛИТЕРА ТУРЫ
1. Недоступ А. А., Наумов В. А., Ражев А. О., Белых А. В. Математическое моделирование орудий и процессов рыболовства: моногр. Калининград: Изд-во КГТУ, 2013. Ч. 1. 253 с.
2. Недоступ А. А., Ражев А. О. Математическое моделирование орудий и процессов рыболовства: моногр. Калининград: Изд-во КГТУ, 2014. Ч. 2. 249 с.
3. Недоступ А. А., Володько Д. А., Ражев А. О. Гидродинамический расчет рыболовной сети // Сб. тр. V Междунар. науч.-практ. конф. «Инженерные системы-2012». М.: РУДН, 2012. С. 27-31.
4. Nedostup A. A., Razhev A. O. A discrete model of gill nets for static and dynamic problems // Contributions on the theory of fishing gears and related marine systems DEMaT-2013. Proseedings of the 11th International workshop (Germany, Rostock, 9-12 October, 2013). Germany, Rostock, 2013. Vol. 8. P. 13-22.
5. Ильин В. А., Садовничий В. А., Сендов Б. Х. Математический анализ. М.: Проспект, ТК Велби, 2004. Ч. 1. 672 с.
6. Patterson D. A., Hennessy J. L. Computer Organization and Design. The Hardware/Software Interface. Morgan Kaufmann, 2011. 914 p.
Статья поступила в редакцию 24.10.2016
ИНФОРМАЦИЯ ОБ АВТОРАХ
Недоступ Александр Алексеевич — Россия, 236022, Калининград; Калининградский государственный технический университет; канд. техн. наук, доцент; зав. кафедрой промышленного рыболовства; nedostup@klgtu.ru.
Ражев Алексей Олегович - Россия, 236022, Калининград; Калининградский государственный технический университет; аспирант кафедры промышленного рыболовства; progacpp@live.ru.
A. A. Nedostup, A. O. Razhev
OPTIMIZATION OF CALCULATION OF FORCES OF HYDRODYNAMIC RESISTANCE IN COMPUTER MODELING OF THE AQUACULTURE CAGES
Abstract. The proposed algorithm for calculating the forces of hydrodynamic resistance of the cage of aquaculture, optimized to run on modern PCs, is based on the decomposition of forces dependencies on the flow rate in a Taylor series of the 4-th order. The calculation was performed using SIMD-instructions simultaneously for two hydrodynamic coefficients: hydrodynamic coefficient for network situated parallel to the water flow rate, and the hydrodynamic network situated perpendicular to the water flow. With the help of the computer program numerical experiment was conducted and performance characteristics of the optimized algorithm, compared with the algorithm, using the operation of exponentiation, were specified. The runtime measurement was made for one billion units. The speed of then optimized algorithm was 35 times higher. Approximation of the result is required only when moving the node of point-mass model at speed below 0.2 m/s.
Key words: hydrodynamic resistance, discrete mathematical model, computer simulation, SIMD, Taylor series.
REFERENCES
1. Nedostup A. A., Naumov V. A., Razhev A. O., Belykh A. V. Matematicheskoe modelirovanie orudii i protsessov rybolovstva [Mathematical modeling of tools and processes of fishing]. Kaliningrad, Izd-vo KGTU, 2013. Part 1. 253 p.
2. Nedostup A. A., Razhev A. O. Matematicheskoe modelirovanie orudii i protsessov rybolovstva [Mathematical modeling of tools and processes of fishing] Kaliningrad, Izd-vo KGTU, 2014. Part 2. 249 p.
3. Nedostup A. A., Volod'ko D. A., Razhev A. O. Gidrodinamicheskii raschet rybolovnoi seti [Hydrodynamic calculation of fishing net]. Sbornik trudov VMezhdunarodnoi nauchno-prakticheskoi konferentsii «Inzhenernye sistemy-2012». Moscow, RUDN, 2012. P. 27-31.
4. Nedostup A. A., Razhev A. O. A discrete model of gill nets for static and dynamic problems. Contributions on the theory of fishing gears and related marine systems DEMaT-2013. Proseedings of the 11th International workshop (Germany, Rostock, 9-12 October, 2013). Germany, Rostock, 2013, vol, pp. 13-22.
5. Il'in V. A., Sadovnichii V. A., Sendov B. Kh. Matematicheskii analiz [Mathematical analysis]. Moscow, Prospekt, TK Velbi Publ., 2004. Part 1. 672 p.
6. Patterson D. A., Hennessy J. L. Computer Organization and Design. The Hardware/Software Interface. Morgan Kaufmann, 2011. 914 p.
The article submitted to the editors 24.10.2016
INFORMATION ABOUT THE AUTHORS
Nedostup Alexander Alekseevich — Russia, 236022, Kaliningrad; Kaliningrad State Technical University; Candidate of Technical Sciences, Assistant Professor; Head of the Department of Commercial Fishery; nedostup@klgtu.ru.
Razhev Alexey Olegovich - Russia, 236022, Kaliningrad; Kaliningrad State Technical University; Postgraduate Student of the Department of Commercial Fishery; progacpp.
