CUDA РЕАЛИЗАЦИЯ НЕЛОКАЛЬНОЙ ЯВНОЙ КОНЕЧНО-РАЗНОСТНОЙ СХЕМЫ ДЛЯ РЕШЕНИЯ ЗАДАЧ ДРОБНОЙ ДИНАМИКИ ПРОЦЕССОВ С НАСЫЩЕНИЕМ

CUDA РЕАЛИЗАЦИЯ НЕЛОКАЛЬНОЙ ЯВНОЙ КОНЕЧНО-РАЗНОСТНОЙ СХЕМЫ ДЛЯ РЕШЕНИЯ ЗАДАЧ ДРОБНОЙ ДИНАМИКИ ПРОЦЕССОВ С НАСЫЩЕНИЕМ

Твёрдый Дмитрий Александрович

канд. физ.-мат. наук, науч. сотр., Камчатский государственный университет им. Витуса Беринга,

РФ, г. Петропавловск-Камчатский

Паровик Роман Иванович

д-р физ.-мат. наук, вед. науч. сотр., Камчатский государственный университет им. Витуса Беринга,

РФ, г. Петропавловск-Камчатский

CUDA IMPLEMENTATION OF A NONLOCAL EXPLICIT FINITE DIFFERENCE SCHEME FOR SOLVING PROBLEMS OF FRACTIONAL DYNAMICS OF PROCESSES WITH SATURATION

Dmitrii Tverdyi

PhD of Phys.-Math. Sci., Researcher, Kamchatka State University,

Russia, Petropavlovsk-Kamchatsky

Roman Parovik

Doctor of Phys.-Math. Sci., Leading Researcher, Kamchatka State University,

Russia, Petropavlovsk-Kamchatsky 

АННОТАЦИЯ

В статье представлены методика и результаты распараллеливания численного алгоритма нелокальной явной конечно-разностной схемы для решения задачи Коши, в основе которой лежит модельное дробное уравнение Риккати с непостоянными коэффициентами. Инструментом для разработки параллельного численного метода была выбрана программно-аппаратная архитектура параллельных вычислений CUDA, на языке программирования C. Реализация и тестирование параллельного численного алгоритма были проведены на серийном ноутбуке HP Pavilion Gaming Laptop 15-dk0xxx. Также рассмотренные в этой статье алгоритмы были запущенны на суперкомпьютере NVIDIA DGX STATION с производительностью в 500 терафлопс, Института математики им. В.И. Романовского. В результате было показано на тестовых примерах, что параллельная версия численного алгоритма может давать существенный рост производительности 2-5 раз в сравнении с последовательной версией алгоритма.

ABSTRACT

The article presents a technique and results of parallelization of the numerical algorithm of a non-local explicit finite-difference scheme for solving the Cauchy problem, which is based on a model fractional Riccati equation with non-constant coefficients. The software and hardware architecture of parallel computing CUDA, in the C programming language, was chosen as a tool for developing a parallel numerical method. Implementation and testing of a parallel numerical algorithm were carried out on a serial HP Pavilion Gaming Laptop 15-dk0xxx. Also, the algorithms discussed in this article were run on the NVIDIA DGX STATION supercomputer with a performance of 500 teraflops, the Institute of Mathematics. V.I. Romanovsky. As a result, it was shown on test examples that the parallel version of the numerical algorithm can provide a significant performance increase of 2-5 times in comparison with the sequential version of the algorithm.

Ключевые слова: динамические процессы, дробная динамика, эффекты памяти, оператор дробного дифференцирования, уравнение Риккати, CUDA

Keywords: dynamic processes, fractional dynamics, memory effects, fractional differentiation operator, Riccati equation, CUDA

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

В настоящей работе был разработан эффективный параллельный численный алгоритм основный на технологии CUDA от компании NVIDIA [3]. Данный алгоритм реализует нелокальную явную конечно-разностную схему (EFDS) на основе метода Эйлера [4], который используется для численного решения задачи Коши для уравнения Риккати.

Уравнение Риккати [5] имеет свои приложения во многих областях науки. Но в рамках наших задач важно то, что уравнение Риккати хорошо описывает процессы, которые подчиняются логистическому закону [6], а также применимо в логистических моделях, целью которых является определение времени насыщения (выхода на плато) [7].  В том числе и дробное уравнение Риккати с оператором типа Герасимова-Капуто [8], как постоянного  порядка [9] так и переменного  порядка во времени [4, 10], в место обычной производной 1-го порядка. Как показано авторами в работах [11, 12] дробное уравнение Риккати, действительно может быть применимо в моделях, где предполагается выход на некий уровень насыщения, а также нужно учитывать эффекты памяти динамической системы. Подробнее об этом можно узнать из статей [4, 12].

В частности, разработанный и изученный авторами математический аппарат численных решений [4, 12] дробного уравнения Риккати применяется при моделировании динамики изменения объёмной активности радона (RVA). Подробнее об этом можно узнать из [13, 14]. Кратко укажем только то, что RVA как раз имеет тенденции к выходу на уровень насыщения, что может описываться моделью с использованием оператора типа Герасимова-Капуто постоянного порядка [15]. И такая математическая модель модет помочь в понимании процессов, происходящих в подземных водах и подпочвенном воздухе перед землетрясениями, с целью их прогноза. А при использовании оператора типа Герасимова-Капуто переменного порядка, модель уже способна описывать некоторые аномальные режимы RVA, также иногда наблюдаемые перед сильными землетрясениями.

Но дальнейшее развитие этих приложений, требует решения обратной задачи, для установления по известным данным процесса, корректных значений параметров модели их размерностей, а также позволит с большей точностью говорить о процессах, предшествующих выходу радона на поверхность. Такие алгоритмы обратных задач потребуют частого пересчёта численных алгоритмов, реализующих модель. Очевидно, используемые в работах [4, 12] последовательные алгоритмы будут не эффективны, т.к. чрезвычайно сильно увеличат конечное время решения поставленной задачи. А значит возникает необходимость в разработке эффективных параллельных численных алгоритмов, используемых моделью.

Будем кратки т.к. данная постановка задачи Коши и её решение, подробно изучены авторами в работах [4, 12]. А также приводится в статье [16] посвящённой OpenMP реализации [17, 18] этого алгоритма, с которым в том числе будем проводить сравнение. Поэтому представим сразу дискретный вариант нелокальной явной конечно-разностной схемы (EFDS) реализуемый в CUDA алгоритмах:

                   (1)

где C – известная константа.

В работе [4], исследованы вопросы сходимости и устойчивости EFDS для дискретной задачи Коши (1). Показано что EFDS устойчива при условии и сходятся с первым порядком точности .

Для оценок вычислительной эффективности используем пример 5 из работы [4], результат которого представлен на рис. 1, для чего рассмотрим дискретную задачу (1) со следующими параметрами:

- переменный показатель порядка дробной производной, гладкая монотонно убывающая функция;

 - гладкая периодическая функция;

- гладкая монотонно убывающая функция;

- гладкая периодическая функция;

- начальное условие задачи Коши;

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

Рисунок 1. Результат численного решения задачи по схеме (1)

Для решения задачи разработки эффективного параллельного численного алгоритма, одной из технологий была выбрана OpenMP - открытый стандарт программного интерфейса приложений для параллельных систем с общей памятью [19], состоящий из набора директив для компиляторов, библиотек функций и набора переменных окружения, разрабатываемый ARB (ARchitecture Board) примерно с 1997 года. Стандарт OpenMP разрабатывался для эффективного программирования на SMP-системах (многопроцессорные суперкомпьютеры, кластеры). Однако в наше время с разработкой многоядерных процессоров, OpenMP имеет спецификации и для них. В следствии разработка эффективных параллельных программных алгоритмов стала возможная даже на персональных домашних компьютерах.

OpenMP официально поддерживается языками программирования С, С++ и нас интересует именно С, в частности из-за универсальности языка C, разработанного для использования в системном программировании [17]. И как следствие, отсутствие в нем строгих ограничений при работе с памятью. А также из-за возможности совместной компиляции программного кода, содержащего функции и директивы как: OpenMP, так и CUDA С. Части программного кода написанного на CUDA C, компилируются через nvcc 11.8, а части написанные просто на C включая OpenMP через gcc 11.2.0.

Действительно, с программной точки зрения архитектура CUDA это расширение языков программирования C и C++, но с некоторыми дополнительными аннотациями для отличия кода: выполняемого на процессоре и на устройстве. И для реализации программ под GPU компания NVIDIA выпустила свои расширения для языка C, компилятор nvcc для сборки таких программ, ввела в обиход новое расширение .cu для файлов, которые содержат CUDA вызовы.

Для более эффективного распараллеливания алгоритмов потребуется вычисление в параллельных областях кода, а главное после, хранение значительного объёма данных в RAM памяти. Это может потребовать хранения данный в многомерных массивах (матрицах), с большими затратами RAM. И при дальнейшем анализе эти затраты RAM мы не будем учитывать.

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

Например в наших рассмотренных случаях, при  численное моделирование по EFDS (1) алгоритмом с использованием OpenMP, потребует расчёта весового коэффициента , где , , представляющего собой в памяти нижне-треугольную матрицу:

                                      (2)

и т.к. для  мы определяем тип данных float, элемент которого занимает в памяти 4 байта, то  из  элементов, займёт в памяти:

Однако при использовании технологии CUDA при работе с многомерными массивами, причем именно напрямую с памятью, массивы представляются в линеаризованном виде т.е., например  при . А значит работать с памятью выделенной в GPU RAM под нижне-треугольную матрицу так просто не получится. Придется определять  как квадратную матрицу  в GPU RAM.

Более того результат вычисленный на GPU и хранимый в его памяти, нужно скопировать уже на CPU RAM т.е. в основную оперативную память компьютера, с помощью специальных функций языка CUDA C. Где CPU сможет продолжить работу алгоритма с этими значениями. А значит нужно и на CPU RAM тоже выделить память под квадратную матрицу . Даже если новая, добавочная половина матрицы, заполнена 0 и её элементы в расчётах никак не участвуют. Тогда затраты памяти существенно возрастают:

А так как среди 2-х используемых нами компьютерных систем и их компонентов, наименьшим объёмом памяти обладает GPU : NVIDIA GeForce GTX 1650 (4102 [Mb]), и учитывая тот факт что эксперименты на всех 2-х системах должны быть с одинаковыми параметрами, получаем что для всех экспериментов, для всех алгоритмов распараллеливания EFDS (1) параметр  будет максимальным. Несмотря на то, что последовательный алгоритм не вычисляет  и способен работать при .

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

Замеры времени, потраченного на выполнение кода на основе CUDA, OpenMP и последовательного алгоритма будем проводить с использованием функции gettimeofday() языка С библиотеки <time.h>, которая возвращает сколько прошло секунд и микросекунд с полуночи 1 января 1970 года. Это один из наиболее точных методов подсчета времени выполнения программ.

Результаты оценки вычислительной эффективности согласно [19] на основе таких понятий как: время работы, ускорение, эффективность и стоимость, и представим на рис. 2-5 далее:

Рисунок 2. Среднее время расчёта в секундах при разном с - числе потоков

Рисунок 3. Ускорение расчёта при разном с - числе потоков

Рисунок 4. Эффективность расчёта при разном с - числе потоков

Рисунок 5. Стоимость алгоритма при разном с - числе потоков

Все графики на рис. 1-6 представленные в статье построены с помощью gnuplot -- свободного программного обеспечения для создания графиков.

В заключении можно сказать, что были разработан эффективный параллельные численный CUDA алгоритм, для решения уравнения Риккати с переменными коэффициентами и с производной дробного переменного порядка типа Герасимова-Капуто. Разработанные алгоритмы были запущены на 2 различных системах : ноутбуке с CPU Intel(R) Core(TM) i5-9300H и GPU NVIDIA GeForce GTX 1650, а также на суперкомпьютере при ИМ им. В.И. Романовского АН РУз с CPU Intel(R) Xeon(R) CPU E5-2698 v4 и GPU 4 х NVIDIA Tesla V100-DGXS.

Был проведён анализ вычислительной эффективности разработанного численного алгоритма явной схемы, где показан существенный прирост производительности в 3-5 раз. Из результатов анализа видно, что при увеличении числа процессоров для параллельной обработки, эффективность использования процессоров существенно снижается. Обе системы выдают искомый результат за приблизительно равное время вычисления, а увеличение числа потоков параллельного алгоритма более чем 15-20 не даёт прироста производительности.

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

Исследование выполнено рамках гранта "Развитие математических моделей дробной динамики с целью исследования колебательных процессов и процессов с насыщением" МД-758.2022.1.1 в КамГУ им. Витуса Беринга.

Список литературы:

1. Rauber T., Runger G. Parallel Programming for Multicore and Cluster Systems. Second edition // New York: Springer. 2013. P.513.
2. Miller R., Boxer L. Algorithms Sequential and Parallel: A Unified Approach. Third edition // Boston: Cengage Learning. 2013. P. 417.
3. Sanders J., Kandrot E. CUDA by Example: An Introduction to General-Purpose GPU Programming // London: Addison-Wesley Professional. 2010. Vol. 1. P. 311. ISBN: 0131387685
4. Tverdyi D.A., Parovik R.I. Investigation of Finite-Difference Schemes for the Numerical Solution of a Fractional Nonlinear Equation // Fractal and Fractional. 2022. Vol. 6(1). No. 23. P. 1-27. DOI:10.3390/fractalfract6010023
5. Taogetusang, Sirendaoerji, Li S. New application to Riccati equation // Chinese Physics B. 2010. Vol. 19. P. 080303. DOI: 10.1088/1674-1056/19/8/080303
6. Drozdyuk A.V. Logistic curve // Toronto: Choven. 2019. P. 270. ISBN: 978-0-9866300-2-6
7. Kurkin A.A., Kurkina O.E., Pelenovskij E.N. Logistic models of the spread of epidemics // Proceedings of NSTU named R.E. Alekseev. 2020. Vol. 129. P. 9-18.
8. Kilbas A.A., Srivastava H.M., Trujillo J.J.  Theory and Applications of Fractional Differential Equations // Amsterdam: Elsevier Science Limited. 2006. P. 523. ISBN: 9780444518323
9. Tverdyi D.A., Parovik R.I. Mathematical modeling of some logistic laws using the hereditary dynamical Riccati system // In: Proceedings of the 11th All-Russian Scientific Conference with International Participation (May 27–30, 2019). Mathematical modeling and boundary value problems. Samara: Samara State Technical Univ. 2019. P. 348-352.
10. Tvyordyj D.A. Hereditary Riccati equation with fractional derivative of variable order // Journal of Mathematical Sciences. 2021. Vol. 253. No. 4. P. 564-572.  DOI: 10.1007/s10958-021-05254-0
11. Tverdyi D.A., Parovik R.I. Fractional Riccati equation to model the dynamics of COVID-19 coronavirus infection // Journal of Physics: Conference Series. 2021. Vol. 2094. P. 032042. DOI: 10.1088/1742-6596/2094/3/032042
12. Parovik R.I. Tverdyi D.A. Application of the Fractional Riccati Equation for Mathematical Modeling of Dynamic Processes with Saturation and Memory Effect // Fractal and Fractional. 2022. Vol. 6(3). No. 163. P. 1-35. DOI: 10.3390/fractalfract6030163
13. Фирстов П.П., Макаров Е.О., Глухова И.П. и др. Поиск предвестниковых аномалий сильных землетрясений по данным мониторинга подпочвенных газов на Петропавловск-Камчатском геодинамическом полигоне // Геосистемы переходных зон. 2018. Т. 2. № 1. С. 16-32.
14. Hauksson E. Radon content of groundwater as an earthquake precursor: evaluation of worldwide data and physical basis // Journ. Geophys. Res. 1981. Vol. 86. Р. 9397–9410
15. Tverdyi D.A., Parovik R.I., Makarov E.O., Firstov P.P. Research of the process of radon accumulation in the accumulating chamber taking into account the nonlinearity of its entrance // E3S Web Conference. 2020. Vol. 196. No. 02027. P. 1-6. DOI: 10.1051/e3sconf/2020196020278
16. Твёрдый Д.А., Паровик Р.И. Эффективный параллельный численный алгоритм для решения задач дробной динамики: дробная модель процессов насыщения // Проблемы вычислительной и прикладной математики. 2022. Т. 42. № 4. С.19-35.
17. Калиткин Н.Н. Численные методы. 2-е изд. // Санкт-Петербург: БХВ. 2011. C. 592.
18. Гергель В.И. Высокопроизводительные вычисления для многоядерных многопроцессорных систем // Москва: Издательство МГУ. 2010. C.544.
19. Борзунов С.В., Кургалин С.Д., Флегель А.В. Практикум по параллельному программированию: учеб. Пособие // Санкт-Петербург: БХВ. 2017. C. 236. ISBN: 978-5-9909805-0-1
20. Brent R.P. The parallel evaluation of general arithmetic expressions // Journal of the Association for Computing Machinery. 1974. Vol. 21. No. 2. DOI: 10.1145/321812.321815