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

Нефедов Михаил Сергеевич

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

Бакалаврская работа

Математика

Дипломы

Вуз: Мордовский государственный университет имени Н.П. Огарева

ID: 6065e2c4e010f800016c945d

UUID: 8ac70000-752a-0139-ab53-0242ac18000a

Язык: Русский

Опубликовано: больше 3 лет назад

Просмотры: 6

Нефедов Михаил Сергеевич

Источник: федеральное государственное бюджетное образовательное учреждение высшего образования "Национальный исследовательский Мордовский государственный университет им. Н. П. Огарёва"

Комментировать 0

Рецензировать 0

Скачать - 5,0 МБ

Поделиться работой

РЕФЕРАТ Бакалаврская работа содержит 46 страниц, 14 рисунков, 1 таблицу, 11 использованных источников, 2 приложения. ГАЗ, НЕЙРОННАЯ СЕТЬ, ОБУЧЕНИЕ, АППРОКСИМАЦИЯ, МЕТОД HLLC, ЗАДАЧА РИМАНА. Цель работы – разработать нейронную сеть, которая вычисляет дискретные потоки при численном решении задач газовой динамики. При решении задачи использовалась среда разработки Google Colaboratory и языки программирования Python и C++. Степень внедрения – частичная. Область применения – практические занятия по численным методам и машинному обучению. Эффективность – данная работа направлена на увеличение уровня интереса к предмету «Машинное обучение», понимание применения нейронных сетей в решении задач газовой динамики. 4

СОДЕРЖАНИЕ ВВЕДЕНИЕ .............................................................................................................. 6 1. Постановка задачи 7 1.1 Автомодельная картина возникающего течения на плоскости 7 1.2 Алгоритм построения обобщенного решения задачи 9 1.3 Сходимость решения ................................................................................ 12 2. Описание машинного обучения ............................................................... 11118 2.1 Предмет машинного обучения 18 2.2 Описание объектов 18 2.3 Задача регрессии ........................................................................................ 19 2.4 Основные понятия нейронной сети 20 3. Реализация нейронной сети ............................................................................ 23 3.1 Подготовка данных .................................................................................... 23 3.2 Построение и обучение модели ............................................................... 25 3.3 Решение тестовой задачи .......................................................................... 29 3.3.1 Задача Сода........................................................................................... 29 ЗАКЛЮЧЕНИЕ ............................................................................................ 333333 СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ ........................................... 34 ПРИЛОЖЕНИЕ А (обязательное) Листинг программы .................................. 36 ПРИЛОЖЕНИЕ Б (обязательное) Решение методом HLLC ............................ 41 5

ВВЕДЕНИЕ Одной из самых больших проблем при решении задач газодинамики является возможность образования разрывов решения. Данная проблема оказывает огромное влияние на качество решения при конструировании разностных схем. Самым точным и надежным методом численного решения задач газодинамики является схема С. К. Годунова. Метод Годунова базируется на приближении потоков на границах ячеек схемы с помощью точного решения задачи Римана. Зачастую, очень сложно получить точное решение задачи Римана о распаде разрыва. Из-за этого возникло огромное количество более экономичных методов, приближенных к схеме Годунова. Однако, для их применения требуется подробное изучение и доскональное понимание проблемы. Случаются ситуации, что на подробное изучение поставленной задачи нет времени, но посмотреть на приблизительный результат процесса необходимо. В таком случае может помочь машинное обучение. С помощью данного инструмента откроется возможность решения огромного числа задач. Следует понимать, что для проверки адекватности и правильности решения необходимы знания изучаемой сферы. Данная работа стремится продемонстрировать, что такая дисциплина, как машинное обучение и ее методы способны решать некоторые задачи газовой динамики с достаточно высокой точностью. 6

1 Постановка задачи 1.1 Автомодельная картина возникающего течения на плоскости Рассматривается определенная среда (газ). Записывается уравнение состояния (1.1.1) [2, c. 106], в котором определяются следующие переменные: 𝑝 – давление, 𝜌 – плотность, 𝜀 – внутренняя энергия единицы массы, 𝑝0 = 1 𝜘 𝜌0 𝑐02 , 𝜘, 𝜌0 , 𝑐0 – некоторые постоянные величины 𝑝 + 𝜘𝑝0 𝑐02 𝜀= − . (𝜘 − 1)𝜌 𝜘 − 1 (1.1.1) Имеется начальный момент времени (t = 0). Вводятся обозначения: 𝑝1 , 𝜌1 , 𝑢1 , которые характеризуют левое полупространство (x < 0). Параметры 𝑝11 , 𝜌11 , 𝑢11 отвечают за характеристику правого полупространства (x > 0). В данном случае, u – компонента вектора скорости. Для идеального газа 𝑐0 = 0. Известно, что поверхность разрыва образуется, например, при соприкосновении двух масс газа (без перегородки между ними) с различными давлениями. В тоже время, на поверхности разрывов должны выполняться следующие равенства [1]: [𝜌]𝐷 − [𝜌𝑢] = 0, [𝜌𝑢]𝐷 − [𝑝 + 𝜌𝑢2 ] = 0, 2 2 [𝜌(𝜀 + 𝑢 ⁄2)] 𝐷 − [𝜌𝑢 (𝜀 + 𝑢 ⁄2) + 𝑝𝑢] = 0. 7 (1.1.2)

Схематично представляются пять возможных вариантов возникающего течения (рисунок 1). Первые четыре изображения содержат контактный разрыв (КР) [2, c. 107]. Параметром P будут обозначаться значения областей слева, а параметром U справа от КР. Вводятся обозначения для плотности и внутренней энергии, где 𝑅1 , 𝐸1 – значение плотности и внутренней энергии для левой области, а 𝑅11 , 𝐸11 – для правой. УВ КР УВ КР ВР O УВ УВ O КР КР ВР ВР O O ВР Вакуум O Рисунок 1 8 ВР ВР

Исходя из представленных вариантов, вводятся следующие условные обозначения: УВ – ударная волна, ВР – волна разряжения, (𝑝1 , 𝜌1 , 𝑢1 ) – параметры слева, (𝑝11 , 𝜌11 , 𝑢11 ) – параметры справа. 1.2 Алгоритм построения обобщенного решения задачи Для УВ должны выполняться соотношения (1.2.1) [2, c. 108], которые связывают её скорость и величины «за» и «перед» фронтом. Из этого следует равенство 1 1 1 𝐸 − 𝜀 − ( − ) (𝑝 + 𝑃) = 0, 2 𝜌 𝑅 (1.2.1) где (𝑝, 𝜌, 𝜀) – величины перед ударной волной, (P, R, E) – за ударной волной. Исключается 𝜀 и E с помощью (1.1.1) и получается адиабата Гюгонио [2, c. 108] 𝑅=𝜌 (𝜘 + 1)(𝑃 + 𝑝0 ) + (𝜘 − 1)(𝑝 + 𝑝0 ) . (𝜘 − 1)(𝑃 + 𝑝0 ) + (𝜘 + 1)(𝑝 + 𝑝0 ) (1.2.2) Вводится массовая скорость [2, c. 108] 𝑎𝐼 = 𝜌𝐼 (𝑢𝐼 − 𝐷𝐼 ) = 𝑅𝐼 (𝑈 − 𝐷𝐼 ), (1.2.3) затем получается для левой УВ равенство (1.2.4) [2, c. 108] 𝑈 − 𝑢𝐼 + 𝑃 − 𝑝𝐼 = 0, 𝑎𝐼 9 (1.2.4)

𝑎𝐼 = √𝜌𝐼 [ 𝜘+1 2 (𝑃 + 𝑝0 ) + 𝜘−1 2 (1.2.5) (𝑝𝐼 + 𝑝0 )] и для правой [2, c. 108] 𝑎𝐼𝐼 = 𝜌𝐼𝐼 (𝐷𝐼𝐼 − 𝑢𝐼𝐼 ) = 𝑅𝐼𝐼 (𝐷𝐼𝐼 − 𝑈), 𝑈 − 𝑢𝐼𝐼 + 𝑎𝐼𝐼 = √𝜌𝐼𝐼 [ Из-за ВР 𝜘+1 2 𝑃−𝑝𝐼𝐼 𝑎𝐼𝐼 (1.2.7) = 0, (𝑃 + 𝑝0 ) + необходимо (1.2.6) 𝜘−1 2 (1.2.8) (𝑝𝐼𝐼 + 𝑝0 )] . использовать условие непрерывности римановских инвариантов [2, c. 108] [𝑢] ± 2 𝜘−1 [с] = 0, [𝜎(𝑆)] = 0, (1.2.9) где c – скорость звука, S – энтропия, а c и 𝜎 вычисляются [2, c. 109] c =√ 𝜘 𝑝+𝑝0 𝜌 , 𝜎(𝑆)= 𝜘(𝑝+𝑝0 ) 𝜌𝜘 . (1.2.10) Следует подстановка их в (1.2.9) [2, c. 109], получается равенство для левой волны 𝑈 − 𝑢𝐼 + 2 𝑐 𝜘−1 𝐼 𝑃+𝑝0 [1 − (𝑝 +𝑝 ) 𝐼 и для правой в виде [2, c. 109] 10 0 𝜘−1 2𝜘 ] = 0, (1.2.11)

𝑈 − 𝑢𝐼𝐼 + 2 𝑐 𝜘−1 𝐼𝐼 𝑃+𝑝0 [1 − (𝑝 𝐼𝐼 +𝑝0 𝜘−1 2𝜘 ) ] = 0. (1.2.12) Вводятся условные массовые скорости, используя, формулы [2, c. 109] 𝜘−1 𝑎𝐼 = 2𝜘 𝜘−1 𝑎𝐼𝐼 = 2𝜘 𝜌𝐼 𝑐𝐼 𝑃+𝑝0 𝑝𝐼 +𝑝0 𝜘−1 𝑃+𝑝0 2𝜘 ( ) 𝑝𝐼 +𝑝0 1− , 𝑃+𝑝0 𝑝𝐼𝐼 +𝑝0 𝜘−1 𝑃+𝑝0 2𝜘 ( ) 𝑝𝐼𝐼 +𝑝0 1− 𝜌𝐼𝐼 𝑐𝐼𝐼 (1.2.13) . (1.2.14) Итерационный процесс будет выглядеть так: пусть 𝑃 (𝑖−1) – полученное приближение для величины 𝑃. (𝑖−1) Вычисляется 𝑎𝐼 = 𝑎𝐼 (𝑃(𝑖−1) ) по формуле (1.2.5) [2, c. 109], если 𝑃(𝑖−1) ≥ 𝑝𝐼 , или (1.2.13) [2, c. 109], если 𝑃(𝑖−1) < 𝑝𝐼 , а параметр (𝑖−1) 𝑎𝐼𝐼 = 𝑎𝐼𝐼 ( 𝑃(𝑖−1) ) – по формуле (1.2.8) [2, c. 109], если 𝑃(𝑖−1) ≥ 𝑝𝐼𝐼 , или (1.2.14) [2, c. 109], если 𝑃 (𝑖−1) < 𝑝𝐼𝐼 . После этого, приближение 𝑃(𝑖) вычисляется [2, c. 110] 𝑃(𝑖) = 𝜑( 𝑃(𝑖−1) ) = (𝑖−1) 𝑎𝐼𝐼 (𝑖−1) 𝑝𝐼 +𝑎𝐼 (𝑖−1) 𝑝𝐼𝐼 +𝑎𝐼 (𝑖−1) 𝑎𝐼𝐼 (𝑖−1) (𝑖−1) 𝑎𝐼 +𝑎𝐼𝐼 (𝑢𝐼 −𝑢𝐼𝐼 ) . (1.2.15) Величина U находится по формуле [2, c. 110] 𝑈= 𝑎𝐼 𝑢𝐼 +𝑎𝐼𝐼 𝑢𝐼𝐼 +𝑝𝐼 −𝑝𝐼𝐼 𝑎𝐼 +𝑎𝐼𝐼 . (1.2.16) Видно, что процесс сходится не всегда. Чтобы сделать его сходящимся во всех случаях, вводится [2, c. 110] 11

𝑃 (𝑖) 𝑎𝑖−1 𝑃(𝑖−1) + 𝜑( 𝑃(𝑖−1) ) = , 1 + 𝑎𝑖−1 𝜘−1 3𝜘 𝑎𝑖−1 = 1 − 𝑧𝑖−1 𝜘+1 2𝜘 𝑧𝑖−1 (1 − 𝜘−1 − 2𝜘 𝑧𝑖−1 ) 1, если это выражение больше нуля, 0 в противном случае, { 𝑧𝑖−1 𝑃 (𝑖−1) + 𝑝0 = . (𝑝𝐼 + 𝑝0 ) + (𝑝𝐼𝐼 + 𝑝0 ) Таким образом, был рассмотрен алгоритм для получения обобщенного решения. 1.3 Сходимость решения Далее описывается алгоритм, который использует метод касательных Ньютона и обеспечивает быструю сходимость процесса. Необходимо получить уравнение для давления P. Уравнение будет выглядеть следующим образом [2, c. 110]: 𝐹(𝑃) ≡ 𝑓(𝑃, 𝑝𝐼 , 𝜌𝐼 ) + 𝑓(𝑃, 𝑝𝐼𝐼 , 𝜌𝐼𝐼 ) = 𝑢𝐼 − 𝑢𝐼𝐼 , где для значений 𝑘 = 𝐼, 𝐼𝐼 [2, c. 110] 12 (1.3.1)

𝑃−𝑝𝑘 𝜌𝑘 𝑐𝑘 √ 𝑓(𝑃, 𝑝𝑘 , 𝜌𝑘 ) = 2 𝑐 𝜘−1 𝑘 { где 𝜋𝑘 = 𝑃+𝑝0 𝑝𝑘 +𝑝0 , 𝑐𝑘 = √𝜘 𝑝𝑘 +𝑝0 𝜌𝑘 𝜘+1 𝜘−1 𝜋 + 2𝜘 𝑘 2𝜘 𝜘−1 2𝜘 (𝜋𝑘 , если 𝑃 ≥ 𝑝𝑘 , (1.3.2) − 1) , если 𝑃 ≤ 𝑝𝑘 , . Вычисляется первая и вторая производная функции 𝑓(𝑃, 𝑝𝑘 , 𝜌𝑘 ) по переменной P, при всех 𝑃 > − 𝑝0 получится [2, c. 111] 𝑓′(𝑃, 𝑝𝑘 , 𝜌𝑘 ) > 0, 𝑓′′(𝑃, 𝑝𝑘 , 𝜌𝑘 ) < 0, (1.3.3) В общем случае график 𝐹(𝑃) схематично изображен ниже (рисунок 2) [2, c. 111]. 𝐹(𝑃) 𝑈уд −𝑝0 𝑝𝐼 𝑝𝐼𝐼 𝑈раз 𝑝 𝑈вак Рисунок 2 Значения функции 𝐹(𝑃) в точках 𝑃 = −𝑝0 , 𝑝𝐼 , 𝑝𝐼𝐼 позволяют определить, какой из вариантов на (рисунок 1) возникает при распаде разрыва. Получается [2, c. 112]: 13

𝐹(𝑝𝐼𝐼 ) = 𝑈уд = 𝑝𝐼𝐼 − 𝑝𝐼 √𝜌𝐼 [𝜘 + 1 (𝑝𝐼𝐼 + 𝑝0 ) + 𝜘 − 1 (𝑝𝐼 + 𝑝0 )] 2 2 , 𝜘−1 𝐹(𝑝𝐼 ) = 𝑈раз 2𝑐𝐼𝐼 𝑝𝐼 + 𝑝0 2𝜘 =− [1 − ( ) ], 𝜘−1 𝑝𝐼𝐼 + 𝑝0 −𝐹(𝑝0 ) = 𝑈вак = − (1.3.4) 2𝑐𝐼 2𝑐𝐼𝐼 − . 𝜘−1 𝜘−1 В зависимости от 𝑢𝐼 − 𝑢𝐼𝐼 возможны ситуации: - 𝑢𝐼 − 𝑢𝐼𝐼 > 𝑈уд . Тогда 𝑃 > 𝑝𝐼𝐼 и 𝑃 > 𝑝𝐼 . УВ распространяется влево и вправо; - 𝑈раз < 𝑢𝐼 − 𝑢𝐼𝐼 < 𝑈уд . Тогда 𝑝𝐼 < 𝑃 < 𝑝𝐼𝐼 . УВ распространяется влево, ВР вправо; - 𝑈вак < 𝑢𝐼 − 𝑢𝐼𝐼 < 𝑈раз . Тогда 𝑝0 < 𝑃 < 𝑝𝐼 . ВР распространяется влево и вправо; - 𝑢𝐼 − 𝑢𝐼𝐼 < 𝑈вак . Получается вакуум, уравнение (1.3.1) [7] не имеет вещественного корня. Уравнение (1.3.1) имеет единственный корень, если выполняется [2, c. 112] 𝑢𝐼 − 𝑢𝐼𝐼 ≥ 𝑈вак = − 2 (𝑐 + 𝑐𝐼𝐼 ). 𝜘−1 𝐼 (1.3.5) В данном случае наиболее целесообразным является метод касательных Ньютона. Если 𝑃 (𝑖−1) – приближенное значение корня, то [2, c. 112] 𝑃 (𝑖) =𝑃 (𝑖−1) 𝑓(𝑃(𝑖−1) , 𝑝𝐼 , 𝜌𝐼 ) + 𝑓(𝑃(𝑖−1) , 𝑝𝐼𝐼 , 𝜌𝐼𝐼 ) − (𝑢𝐼 − 𝑢𝐼𝐼 ) − . 𝑓′(𝑃(𝑖−1) , 𝑝𝐼 , 𝜌𝐼 ) + 𝑓(𝑃(𝑖−1) , 𝑝𝐼𝐼 , 𝜌𝐼𝐼 ) 14 (1.3.6)

По скорости сходимости метод Ньютона является итерационным процессом второго порядка, убывание погрешности корня 𝑃∗ получается следующим образом [2, c. 113]: 𝑃 (𝑖) 1 𝐹 ′′ (𝑃̃) (𝑖−1) ∗ 2 −𝑃 = − 𝑃 [𝑃 ] , 2 𝐹 ′ (𝑃(𝑖−1) ) ∗ (1.3.7) где 𝑃̃ – некоторая точка на отрезке [𝑃(𝑖−1) , 𝑃(𝑖) ]. Вместо (1.3.1) рассматривается результат его линеаризации [2, c. 113]: 𝑃 − 𝑝𝐼 𝑃 − 𝑝𝐼𝐼 + = 𝑢𝐼 − 𝑢𝐼𝐼 , 𝜌𝐼 𝑐𝐼 𝜌𝐼𝐼 𝑐𝐼𝐼 из него определяется [2, c. 113] 𝑃(0) = 𝑝𝐼 𝜌𝐼𝐼 𝑐𝐼𝐼 + 𝑝𝐼𝐼 𝜌𝐼 𝑐𝐼 + (𝑢𝐼 − 𝑢𝐼𝐼 )𝜌𝐼 𝑐𝐼 𝜌𝐼𝐼 𝑐𝐼𝐼 . 𝜌𝐼 𝑐𝐼 + 𝜌𝐼𝐼 𝑐𝐼𝐼 (1.3.8) Равенство (1.3.8) [2, c. 113] удовлетворяет формуле (1.3.7) [2, c. 113] и дает приближение для искомого корня 𝑃∗ снизу. Благоприятные результаты дает замена (1.3.1) на уравнение [2, c. 113] 𝑃 − 𝑝𝐼 √𝜃𝐼 𝜌𝐼 (𝑃 + 𝑝0 ) + 𝑃 − 𝑝𝐼𝐼 √𝜃𝐼𝐼 𝜌𝐼𝐼 (𝑃 + 𝑝0 ) = 𝑢𝐼 − 𝑢𝐼𝐼 , (1.3.9) где 𝜃𝐼 , 𝜃𝐼𝐼 – const. Тогда равенство (1.3.1) [2, c. 114] принимает вид − 2 𝜘−1 𝑐𝐼 [1 − ( 𝑃+𝑝0 𝑝𝐼 +𝑝0 𝜘−1 2𝜘 ) 2 𝑃+𝑝0 ] − 𝜘−1 𝑐𝐼𝐼 [1 − (𝑝 𝐼𝐼 +𝑝0 15 ) 𝜘−1 2𝜘 ] = 𝑢𝐼 − 𝑢𝐼𝐼 , (1.3.10)

и его решение 𝑃∗ вычисляется по формуле [2, c. 114] 2𝜘 𝑢𝐼 − 𝑢𝐼𝐼 − 𝑈вак 𝜘−1 𝑃∗ = (𝑝𝐼 + 𝑝0 ) [ − 𝑝0 . ] 𝑈раз − 𝑈вак (1.3.11) Чтобы завершить описание алгоритма, нужно привести остальные формулы для других величин, описывающие схемы на (рис. 1). Их можно получить, после нахождения давления P на контактном разрыве. Понадобится уравнение (1.2.16) [2, c. 114], приведенное выше. В нем 𝑎𝐼 , 𝑎𝐼𝐼 – массовые скорости. Если УВ слева, то ее скорость вычисляется как [2, c. 114] 𝐷𝐼 = 𝑢𝐼 − 𝑎𝐼 , 𝜌𝐼 (1.3.12) а плотность (𝑅𝐼 ) вычисляется по адиабате Гюгонио (1.2.2) [2, c. 115]: 𝑅𝐼 = 𝜌𝐼 (𝜘 + 1)(𝑃 + 𝑝0 ) + (𝜘 − 1)(𝑝𝐼 + 𝑝0 ) 𝜌𝐼 𝑎𝐼 = . (𝜘 − 1)(𝑃 + 𝑝0 ) + (𝜘 + 1)(𝑝𝐼 + 𝑝0 ) 𝑎𝐼 − 𝜌𝐼 (𝑎𝐼 − 𝑈) (1.3.13) Если УВ слева, то ее скорости вычисляются [2, c. 115]: 𝐷𝐼 = 𝑢𝐼 − 𝑐𝐼 , где 𝑐𝐼∗ = 𝑐𝐼 + 𝜘−1 2 𝐷𝐼∗ = 𝑈 − 𝑐𝐼∗ , (1.3.14) (𝑢𝐼 − 𝑈), а плотность 𝑅𝐼 слева – по формуле [2, c. 115] 𝑅𝐼 = 𝜘 𝑃 + 𝑝0 . (𝑐𝐼∗ )2 16 (1.3.15)

Если УВ справа, то ее характеристики вычисляются [2, c. 115]: 𝐷𝐼𝐼 = 𝑢𝐼𝐼 − 𝑎𝐼𝐼 𝜌𝐼𝐼 , 𝑅𝐼𝐼 = 𝜌𝐼𝐼 𝑎𝐼𝐼 𝑎𝐼𝐼 −𝜌𝐼𝐼 (𝑎𝐼𝐼 −𝑈) . (1.3.16) Если волна является ВР, то [2, c. 115] 𝐷𝐼𝐼 = 𝑢𝐼𝐼 + 𝑐𝐼𝐼 , 𝑐𝐼𝐼∗ = 𝑐𝐼𝐼 + 𝜘−1 (𝑢𝐼𝐼 − 𝑈), 2 (1.3.17) 𝐷𝐼𝐼∗ = 𝑈 + 𝑐𝐼𝐼∗ , 𝑅𝐼𝐼 = 𝜘 𝑃 + 𝑝0 . (𝑐𝐼𝐼∗ )2 Таким образом, был описан алгоритм построения обобщенного решения задачи о распаде начального разрыва. 17

2 Описание машинного обучения 2.1 Предмет машинного обучения Машинное обучение (Machine Learning) – это класс методов искусственного интеллекта, чертой которого является поиск методов решения задач путем обучения в процессе решения сходных задач [3]. Для построение таких методов используются различные разделы математики, например: численные методы, методы оптимизации, теория вероятностей и др. Как и любая другая дисциплина, машинное обучение имеет определенную специфику. Одна из самых глобальных – переобучение. Переобучение в МО – это явление, при котором построенная модель показывает адекватную работу и правильные ответы, используя данные, полученные в процессе обучения, но работает некорректно на конкретных примерах. Многие методы МО разрабатывались, как альтернатива уже известным подходам решения. 2.2 Описание объектов Одно из направлений машинного обучения связано с различными задачами, например, даны: - набор 𝑋; - набор 𝑌. Набор X представляет собой изначальные объекты, Y – ответы. 18

Имеется некая зависимость (функциональная) 𝑓: 𝑋 → 𝑌 между набором X и Y, но она неизвестна. Дана только совокупность 𝑆 пар вида (X, Y), которая называется обучающая выборка (training sample): 𝑆 = {(𝑥𝑖 , 𝑦𝑥𝑖 = 𝑓(𝑥𝑖 )) ∈ X × Y | i = 1, … , l}. Нужно найти приближенный вид 𝑓 путем построения аппроксимирующей функции 𝑎𝑆 ∶ 𝑋 → 𝑌, такой, что ∀ 𝑥 ∈ 𝑋 𝑎𝑆(𝑥) ≈ 𝑓(𝑥). Объекты в машинном обучении могут иметь различное происхождение, например: стулья, люди, часы и т.д. Эти признаки определяются решаемыми задачами. Важно правильно подобрать адекватное множество объектов, что является нетривиальной задачей. 2.3 Задача регрессии В задачах регрессии выборка имеет вид 𝑆 = {(𝑥1 , 𝑦1 ), … , (𝑥𝑙 , 𝑦𝑙 )} ∈ 𝑅𝑛 × 𝑅, т.е. ответами являются действительные числа либо вектора действительных чисел. Задача регрессии заключается в том, чтобы по известным значениям 𝑦𝑖 определенной неизвестной функции в заданных точках 𝑥𝑖 (𝑖 = 1, . . . , 𝑙) предсказать значения функции в других точках. Основу регрессионного анализа составляет метод наименьших квадратов. Для примера рассматривается задача линейной регрессии: Пусть дана выборка 𝑆 = {(𝑥1 , 𝑦1 ), … , (𝑥𝑙 , 𝑦𝑙 )} ∈ 𝑅𝑛 × 𝑅. Нужно найти аппроксимирующую функцию 𝑎𝑆 (𝑥) в виде 〈𝑥, 𝑤〉 + 𝑤0 , где 𝑤 ∈ 𝑅𝑛 , 𝑤0 ∈ 𝑅, и должно выполнять условие: 19

𝑙 𝑙 𝑖=1 𝑖=1 1 1 𝑄(𝑎𝑆 ) = ∑ 𝐿(𝑎𝑆 , 𝑥𝑖 ) = ∑(𝑎𝑆 (𝑥𝑖 ) − 𝑦𝑖 )2 → 𝑚𝑖𝑛. 𝑙 𝑙 2.4 Основные понятия нейронной сети Нейронная сеть — это последовательность нейронов, соединенных между собой синапсами [5]. Данный термин был позаимствован из такой науки, как биология. Работа нейронной сети схожа с работой человеческого мозга. Соответственно, компьютер, как и мозг способен анализировать, получать и запоминать информацию, подающуюся на входе. Нейронные сети помогают решать сложные математические задачи, которые зачастую сложно или невозможно решить аналитически в данный момент. Существуют различные виды задач, которые решаются при помощи нейронных сетей. Необходимо рассмотреть самые распространенные и важные из них: 1. Задача классификации – распределение данных по разным группам (параметрам, структуре, применимости и т.д.); 2. Задача предсказания – возможность предугадывать следующий шаг, из, каких-либо, наборов данных; 3. Задача распознавания – возможность находить или определять определенный объект, группу из нескольких. Чаще всего это используется в распознавании человеческих лиц, номеров автомобилей и т. д.; 4. Задача регрессии – это задача наилучшего приближения функции отображения f от входных переменных X к непрерывной выходной переменной y. 20

Рассмотрение состава нейронной сети: Нейрон – это некоторая единица, получающая информацию, производящая вычисления и, передающая данные дальше. Они бывают трех базовых видов: 1. входной (синий); 2. скрытый (черный); 3. выходной (красный). Более наглядное взаимодействие нейронов представлено на рисунке 3. Рисунок 3 Работает это по принципу «слоеного пирога». Информация поступает на входной слой, затем поступает в скрытый слой, где она обрабатывается. Количество скрытых слоев зависит от задачи, их количество может быть бесконечно большим. Выходной слой выдает и выводит результат работы сети. Числовое значение нейрона находится в интервале от 0 до 1, либо от -1 до 1. Для обработки числового значения, выходящего за данный интервал, требуется единицу разделить на это число. Сами нейроны связываются между собой синапсами. Синапсы обладают таким параметром как вес, из-за него информация, подающаяся на вход, изменяется, когда передается от одного нейрона к другому. Приоритетность передачи информации зависит от веса нейрона, у кого он больше, та информация будет доминирующей в следующем нейроне. 21

На рисунке 4 изображена основная идея работы нейронной сети, где I – входные нейроны, O – скрытые нейроны, w – веса нейронов. Например, путь на входе будет 0 и 1. 𝑤1 = 0.3, 𝑤2 = 0.7, тогда входные данные нейрона 𝑂1 будут выглядеть: 0.3 ∗ 1 + 0 ∗ 0.7 = 0.3. Функция активации представляет из себя способ нормализации входных данных. Самые распространенные функции активации: - линейная; - сигмоид; - гиперболический тангенс. Отличие функций активаций – диапазон значений [4]. 𝐼1 𝑤1 𝑂1 𝐼2 𝑤2 𝑂1𝑖𝑛𝑝𝑢𝑡 = (𝐼1 ∗ 𝑤1 ) + (𝐼2 ∗ 𝑤2 ) 𝑂1𝑜𝑢𝑡𝑝𝑢𝑡 = 𝑓𝑎𝑐𝑡𝑖𝑣𝑎𝑡𝑖𝑜𝑛 (𝑂1𝑖𝑛𝑝𝑢𝑡 ) Рисунок 4 22

3 Реализация программы 3.1 Подготовка данных При подготовки данных для обучения нейронной сети, использовалась подпрограмма для решения задачи о распаде разрыва, разработанная в ИПМ им. М. В. Келдыша РАН [7]. На вход программы подается 10 параметров: RB, PB, UB, VB, WB, RE, PE, UE, VE, WE, где R – плотность, P – давление, E – энергия, UB, VB, WB и UE, VE, WE – компоненты вектора скорости. Индексы B и E – соответственно параметры «слева» и «справа». На выходе получаются параметры: RI, EI, PI, UI, VI, WI, где E – энергия. Переводя на более простой язык, возникает задача: по 10 параметрам, полученным на входе, необходимо предсказать выходные данные. Как было описано выше в теоретическом материале в пункте 2.3, данную проблему поможет решить задача регрессии. В задаче регрессии, необходимо дать прогноз какого-либо непрерывного значения, в данном случае – выходных параметров. В данной подпрограмме сгенерируем случайным образом входные параметры, например, 2000 значений, чтобы получить выходные параметры. Рассматривается полученный набор данных (рисунок 5). Рисунок 5, лист 1 – Начальный набор данных 23

Рисунок 5, лист 2 – Начальный набор данных После того как отделили от сгенерированных данных прогнозные значения, необходимо разделить полученный data set на обучающую и тестовую выборки. Методом проб и ошибок замечено, что самая благоприятная размерность тестовой выборки составляет 20 %. Для более точного обучения модели необходимо нормализовать признаки (рисунок 6). Для того, чтобы сделать нормализацию, необходимо стандартизировать данные. Стандартизация означает преобразование данных таким образом, чтобы их среднее значение являлось нулем, а дисперсия единицей. Выражаясь языком программы, требуется: из каждого признака вычесть его среднее значение и разделить на стандартное отклонение. Нужно учесть, что модель может сходится и без нормализации признаков, но тогда итоговая модель может зависеть от единиц измерения и обучение заметно усложнится. Рисунок 6 Выражение axis=0 на рисунке 6 означает, что рассчитывается среднее значение и стандартное отклонение не по всему набору данных, а отдельно для 24

каждого признака. Данные операции проделываются на данных для обучения и тестирования. Таким образом, в данном подразделе были подготовлены данные для последующего обучения модели. 3.2 Построение и обучение модели Вся реализация нейронной сети была выполнена в Google Colaboratory [9]. Основные библиотеки, упрощающие построение нужной нейронной сети на языке Python: keras, tensorflow. Tensorflow – это библиотека программного обеспечения с открытым исходным кодом для сложных вычислений. Главная особенность в данной работе – поддержка машинного обучения и глубокого обучения [10]. Keras – один из самых распространенных и известных фреймворков глубокого обучения [8], который может работать поверх Tensorflow [11]. Для построения модели использовалась Sequential (последовательная) модель с двумя полносвязными скрытыми слоями (рисунок 7). Рисунок 7 25

После построение модели необходимо ее скомпилировать. При компиляции мы должны указать функцию ошибки и метрики качества работы сети. В качестве функции ошибки, которая будет использоваться оптимизатором в алгоритме обратного распространения ошибки выбирается: mse – Mean squared error (среднеквадратичная ошибка). Данная функция ошибки часто используется для задачи регрессии. В качестве метрики выберем: mae – Mean absolute error (средняя абсолютная ошибка). С помощью mae можно узнать, насколько сильно отличается число, которое выдала нам нейронная сеть от правильного ответа. Обучение полученной модели произведем на 1000 эпохах (рисунок 8). Рисунок 8 После обучения модели можно заметить, что средняя абсолютная ошибка равняется 0.0056 (рисунок 9). Рисунок 9 Для наглядности можно визуализировать обучение модели. Получаются графики (рисунок 10). 26

Рисунок 10 Данный график показывает достаточно серьезную деградацию ошибки валидации на отрезке (0; 20) эпох обучения. После двадцатой эпохи результат стабилизируется в лучшую сторону. Для устранения данной проблемы, можно обновить метод model.fit так, чтобы автоматически прекращалось обучение, когда ошибка валидации Val loss прекращала улучшаться. Используется функция EarlyStopping callback, которая проверяет показатели улучшения после каждой эпохи (рисунок 11). 27

Рисунок 11 Получаются следующие графики (рисунок 12). Рисунок 12 28

Последним этапом рассмотрения полученной модели является – прогнозирование полученных значений, используя данные из тестовой выборки (рисунок 13). Рисунок 13 Данная модель дает достаточно хорошие предсказания. Для проверки ее на практике, необходимо применить полученную нейронную сеть в решении некоторой конкретной задачи. 3.3 Решение тестовой задачи 3.3.1 Задача Сода В качестве примера, применим полученную нейронную сеть для решения задачи Сода. Данная задача позволяет проверить отсутствие нефизических осцилляций решения в расчетах движущихся ударных волн и контактного разрыва, а также корректность воспроизведения профиля плотности [5]. 29

Обученная нейронная сеть была внедрена в программу, написанная на C++, которая решает задачу Сода с помощью метода HLLC [6]. Требуется заменить вызов функции calc_flux_hllc (приложение Б) на полученную нейронную сеть. (а) (б) Рисунок 14, лист 1 – Решение задачи Сода, с помощью нейронной сети: (а) – давление, (б) – плотность, (в) – энергия, (г) – скорость 30

(в) (г) Рисунок 14, лист 2 – Решение задачи Сода, с помощью нейронной сети: (а) – давление, (б) – плотность, (в) – энергия, (г) – скорость Красным цветом обозначается точное решение задачи, синим – решение с помощью нейронной сети на сетке с 800 ячейками. 31

Для вычисления порядка аппроксимации необходимо произвести серию расчетов на сетках различной разрешающей способности. В Таблице 1, для сеток с размерностью 100, 200, 400, 800 ячеек, приводится порядок сходимости решения. Можно заметить, что при очень мелком разбиении сетки, полученное решение имеет достаточно высокие порядки точности: около двух и более. Т а б л и ц а 1 – Порядок сходимости решения Сетка (100; 200) (200; 400) (400; 800) 𝜌 (плотность) 1.812 2.243 2.963 u (скорость) 1.174 2.123 3.323 p (давление) 1.823 2.351 2.991 e (энергия) 1.799 1.997 2.985 Полученные результаты показывают, что нейронная сеть вполне неплохо справилась с данной задачей. Получилось достаточно точно решение, близкое к решению, полученному методом HLLC. 32

ЗАКЛЮЧЕНИЕ В ходе выполнения выпускной квалификационной работы была разработана и реализована нейронная сеть, способная вычислять дискретные потоки при численном решении задач газовой динамики. Вычисление таким способом не требует особого познания в конкретной области поставленной задачи. Для этого необходимо понимать принцип и структуру работы нейронной сети, правильное применение её методов. Вычисление нейронной сетью не имеет проблем, присущих другим схемам на основе приближенного решения задачи Римана. Применение используемого метода на тестовой задаче – задаче Сода, показал, что решение имеет достаточно высокие порядки точности: около двух и более. 33

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ 1. Годунов С. К. Разностные схемы. Введение в теорию / С. К. Годунов, В. С. Рябенький. – 2-е изд. – М. : 1977. – 440 с. 2. Годунов С. К. Численное решение многомерных задач газовой динамики / С. К. Годунов, А. В. Забродин, М. Я. Иванов, А. Н. Крайко, Г. П. Прокопов. – М. : Наука, 1976. – С. 105-115 3. Машинное обучение [Электронный ресурс]. – Режим доступа: https://ru.wikipedia.org/wiki/%D0%9C%D0%B0%D1%88%D0%B8%D0%BD% D0%BD%D0%BE%D0%B5_%D0%BE%D0%B1%D1%83%D1%87%D0%B5% D0%BD%D0%B8%D0%B5 4. Нейронные сети для начинающих [Электронный ресурс]. – Режим доступа: https://habr.com/ru/post/312450 (03.06.2020) 5. Одномерные задачи газовой динамики и их решение при помощи разностных схем высокой разрешающей способности [Электронный ресурс]. – Режим доступа: https://ntv.ifmo.ru/file/article/13702.pdf 6. Сравнение схем с расщеплением потока для численного решения уравнений Эйлера сжимаемого газа [Электронный ресурс]. – Режим доступа: https://mipt.ru/upload/medialibrary/388/13_stranitsy-iz-trudy-mfti_1_37_2018_c122_141.pdf 7. Тишкин В. Ф. Разностные схемы трехмерной газовой динамики задачи о развитии неустойчивости Рихтмайера–Мешкова / В. Ф. Тишкин, В. В. Никишин, И. В. Попов, А. П. Фаворский // Математическое моделирование. – Т.7, №5, 1995. – С. 15–25 8. Deep learning vs machine learning: a simple way to understand the difference [Электронный ресурс]. – https://www.zendesk.com/blog/machinelearning-and-deep-learning/ 34

9. Google [Электронный Colaboratory ресурс]. – https://colab.research.google.com/drive/1lz8P9v7KVdbawtXw0ik6eOLZguUx55F 0 10. Install TensorFlow with pip [Электронный ресурс]. – Режим доступа: https://www.tensorflow.org/install/pip?hl=ru 11. Keras: the Python deep learning API [Электронный ресурс]. – https://keras.io/ 35

ПРИЛОЖЕНИЕ А (обязательное) Листинг программы # Установим библиотеку seaborn для построения парных графиков !pip install -q seaborn # To determine which version you're using: !pip show tensorflow # For the current version: !pip install --upgrade tensorflow # For a specific version: !pip install tensorflow==1.2 # For the latest nightly build: !pip install tf-nightly # -*- coding: utf-8 -*from __future__ import absolute_import, division, print_function, unic ode_literals import import import import import pathlib numpy as np matplotlib.pyplot as plt pandas as pd seaborn as sns import tensorflow as tf from tensorflow import keras from tensorflow.keras import layers print(tf.__version__) from from from from from from from from from pandas import read_csv, DataFrame sklearn.neighbors import KNeighborsRegressor sklearn.linear_model import LinearRegression, LogisticRegression sklearn.svm import SVR sklearn.ensemble import RandomForestRegressor sklearn.metrics import r2_score sklearn.model_selection import train_test_split sklearn.model_selection import cross_val_score tensorflow.keras.models import load_model 36

Продолжение ПРИЛОЖЕНИЯ А dataset_path = keras.utils.get_file("outpuuti800.csv", "https://s155my t.storage.yandex.net/rdisk/7900f3bc6860b94fc8229d9ec336f30393be8f a75525ad1050b2ca922b9121fe/5ed3b5d1/hX1zYKNUYTSLzQn0DUnmH1839wsBXlCJEKHaxuFrTteZxHi1cjj_SZhqbgBTbELC3XbOxuC8_8zYqhzs TuuJw==?uid=344471608&filename=outpuuti800.csv&disposition=attach ment&hash=&limit=0&content_type=text%2Fplain&owner_uid=344471608& fsize=178281&hid=163c3ebb90ec585740f606d525f4e2ce&media_type=spre adsheet&tknv=v2&etag=7df62481e2d328068097950bf5d79149&rtoken=NSj4 hZ9BbneG&force_default=yes&ycrid=naf952cfa489ecdb653b60486f48a27f6fdownloader1e&ts=5a6f1f10c5640&s=8b87d0abf266110859b5d35db64efc5e3 1b9ffad691aafd72704d19bb6b6603e&pb=U2FsdGVkX1_yyQduhORnYuXt6zhBXI Ms4o9E4MHv_gizzNqSvtZkTitKV3iuBXPuJFg4iCTfAJN52AaNTKHSkLLFmcHWGPI8b6YAeLM3Xg") dataset_path # Нумеруем столбцы column_names = ['RB', 'UB', 'PB', 'RE', 'UE', 'PE', 'RI', 'UI', 'PI', 'EI'] raw_dataset = pd.read_csv(dataset_path, names=column_names) dataset = raw_dataset.copy() dataset.head() trg = dataset[['RI', 'UI', 'PI', 'EI']] trn = dataset.drop(['RI', 'UI', 'PI', 'EI'], axis=1) dataset.tail() Xtrn, Xtest, Ytrn, Ytest = train_test_split(trn, trg, test_size=0.2) def norm(x): return (x - Xtrn.mean(axis=0)) / Xtrn.std(axis=0) normed_Xtrn = norm(Xtrn) normed_Xtest = norm(Xtest) normed_Xtrn[:5] def build_model(): model = keras.Sequential([ layers.Dense(64, activation='relu', input_shape=[len(normed_Xtrn.k eys())]), layers.Dense(64, activation='relu'), layers.Dense(4) ]) optimizer = tf.keras.optimizers.RMSprop(0.001) model.compile(loss='mse', optimizer=optimizer, 37

Продолжение ПРИЛОЖЕНИЯ А metrics=['mae', 'mse']) return model import numpy as np k = 4 # кол-во блоков num_val_samples = len(normed_Xtrn) // k # колвл строк в валидационной выборке num_epochs = 500 # кол-во эпох обучения all_scores = [] # сохраняем получившиеся оценки for i in range(k): print('processing fold #', i) # выделение валидационной выборки val_data = normed_Xtrn[i * num_val_samples: (i + 1) * num_val_samp les] val_targets = Ytrn[i * num_val_samples: (i + 1) * num_val_samples] # остальные данные для обучения partial_normed_Xtrn = np.concatenate( [normed_Xtrn[:i * num_val_samples], normed_Xtrn[(i + 1) * num_val_samples:]], axis=0) partial_Ytrn = np.concatenate( [Ytrn[:i * num_val_samples], Ytrn[(i + 1) * num_val_samples:]], axis=0) model = build_model() model.fit(partial_normed_Xtrn, partial_Ytrn, epochs=num_epochs, batch_size=1, verbose=0) # Оцениваем модель model.evaluate(normed_Xtest, Ytest) from keras import backend as K # Some memory clean-up K.clear_session() from keras import models num_epochs = 1000 all_mae_histories = [] for i in range(k): print('processing fold #', i) val_data = normed_Xtrn[i * num_val_samples: (i + 1) * num_val_samp les] val_targets = Ytrn[i * num_val_samples: (i + 1) * num_val_samples] 38

Продолжение ПРИЛОЖЕНИЯ А partial_normed_Xtrn = np.concatenate( [normed_Xtrn[:i * num_val_samples], normed_Xtrn[(i + 1) * num_val_samples:]], axis=0) partial_Ytrn = np.concatenate( [Ytrn[:i * num_val_samples], Ytrn[(i + 1) * num_val_samples:]], axis=0) model = build_model() history = model.fit(partial_normed_Xtrn, partial_Ytrn, validation_data=(val_data, val_targets), epochs=num_epochs, batch_size=1, verbose=0) class PrintDot(keras.callbacks.Callback): def on_epoch_end(self, epoch, logs): if epoch % 100 == 0: print('') print('.', end='') history = model.fit( partial_normed_Xtrn, partial_Ytrn, epochs=num_epochs, validation_split = 0.2, verbose=0, callbacks=[PrintDot()]) hist = pd.DataFrame(history.history) hist['epoch'] = history.epoch hist.tail() model.summary() model.evaluate(normed_Xtest, Ytest) print() def plot_history(history): hist = pd.DataFrame(history.history) hist['epoch'] = history.epoch plt.figure() plt.xlabel('Epoch') plt.ylabel('Mean Abs Error [RI UI PI EI]') plt.plot(hist['epoch'], hist['mae'], label='Train Error') plt.plot(hist['epoch'], hist['val_mae'], label = 'Val Error') plt.ylim([0,0.075]) plt.legend() plt.figure() 39

Окончание ПРИЛОЖЕНИЯ А plt.xlabel('Epoch') plt.ylabel('Mean Square Error [$[RI UI PI EI]^2$]') plt.plot(hist['epoch'], hist['mse'], label='Train Error') plt.plot(hist['epoch'], hist['val_mse'], label = 'Val Error') plt.ylim([0,0.075]) plt.legend() plt.show() plot_history(history) model = build_model() # Параметр patience определяет количество эпох, проверяемых на улучшен ие early_stop = keras.callbacks.EarlyStopping(monitor='val_loss', patienc e=10) history = model.fit(normed_Xtrn, Ytrn, epochs=EPOCHS, validation_split = 0.2, verbose=0, callbacks=[earl y_stop, PrintDot()]) plot_history(history) test_predictions = model.predict(normed_Xtest) plt.scatter(Ytest, test_predictions) plt.xlabel('True Values [RI UI PI EI]') plt.ylabel('Predictions [RI UI PI EI]') plt.axis('equal') plt.axis('square') plt.xlim([0,plt.xlim()[1]]) plt.ylim([0,plt.ylim()[1]]) _ = plt.plot([-10000000000, 10000000000], [-10000000000, 10000000000]) err = np.mean(abs(test_predictions - Ytest)) print(err) f = open("prediction800.txt", "w") with np.printoptions(threshold=np.inf): for arr in test_predictions: for v in arr: f.write(f"{v:.3f}\t") f.write("\n") f.close() 40

ПРИЛОЖЕНИЕ Б (обязательное) Решение методом HLLC # -*- coding: utf-8 -*import sys import math import numpy as np from tensorflow.keras.models import load_model N = 100 K_WENO = 3 lo = K_WENO hi = K_WENO + N - 1 XMIN = 0.0 XMAX = 1.0 H = (XMAX - XMIN) / N ro = [0.] * (N+2*K_WENO) ru = [0.] * (N+2*K_WENO) re = [0.] * (N+2*K_WENO) ro_old = [0.] * (N+2*K_WENO) ru_old = [0.] * (N+2*K_WENO) re_old = [0.] * (N+2*K_WENO) int_ro = [0.] * (N+2*K_WENO) int_ru = [0.] * (N+2*K_WENO) int_re = [0.] * (N+2*K_WENO) r_ = [0.] * (2*K_WENO) p_ = [0.] * (2*K_WENO) u_ = [0.] * (2*K_WENO) TMAX = 0.2 TAU = 1.e-05 * 0.5 TH = TAU / H * 0.5 GAM = 1.4 model = load_model('regression.h5') def main(): 41

Продолжение ПРИЛОЖЕНИЯ Б # /*** Init ***/ for i in range(lo, hi + 1): x = XMIN + (i - lo + 0.5) * H if x < 0.5: ro[i] = 1.0 ru[i] = 0.0 re[i] = 1.0 / (GAM - 1.0) else: ro[i] = 0.125 ru[i] = 0.0 re[i] = 0.1 / (GAM - 1.0) t = 0.0 step = 0 while t < TMAX: step += 1 t += TAU for i in range(lo, hi + 1): ro_old[i] = ro[i] ru_old[i] = ru[i] re_old[i] = re[i] int_ro[i] = 0.0 int_ru[i] = 0.0 int_re[i] = 0.0 # /*** Boundary conditions ***/ for i in range(K_WENO): ro[lo-i-1] = ro[lo+i] ru[lo-i-1] = ru[lo+i] re[lo-i-1] = re[lo+i] ro[hi+i+1] = ro[hi-i] ru[hi+i+1] = ru[hi-i] re[hi+i+1] = re[hi-i] # /*** Fluxes ***/ #print("steps from", (lo-1), "to", hi+1) for i in range(lo-1, hi + 1): rl, pl, ul = 0., 0., 0. rr, pr, ur = 0., 0., 0. fr, fu, fe = 0., 0., 0. for k in range(-K_WENO+1, K_WENO + 1): r_[k+K_WENO-1] = ro[i+k] u_[k+K_WENO-1] = ru[i+k]/ro[i+k] 42

Продолжение ПРИЛОЖЕНИЯ Б p_[k+K_WENO-1] = (re[i+k]0.5*ru[i+k]*ru[i+k]/ro[i+k])*(GAM-1.0) rl, rr = WENO(r_, rl, rr) pl, pr = WENO(p_, pl, pr) ul, ur = WENO(u_, ul, ur) #if i % 100 == 0: #print("step=", i) #fr, fu, fe = calc_flux_hllc(rl, ul, pl, GAM, rr, ur, pr, GAM, fr, fu, fe) answer = model.predict([[rl, ul, pl, rr, ur, pr]]) fr, fu, fp, fe = answer[0] int_ro[i] -= fr int_ru[i] -= fu int_re[i] -= fe int_ro[i+1] += fr int_ru[i+1] += fu int_re[i+1] += fe # /*** Update fields ***/ for i in range(lo, hi + 1): ro[i] += TH*int_ro[i] ru[i] += TH*int_ru[i] re[i] += TH*int_re[i] if step % 100 == 0: print("STEP = " + str(step) + '\n') f = open('result8.txt', 'w') f.write("x,r,u,p,e\n") for i in range(lo, hi + 1): x = XMIN+(i-lo+0.5)*H r = ro[i] u = ru[i] / ro[i] e = (re[i]/ro[i]-0.5*u*u) p = r*e*(GAM-1.) f.write(str(x) + ", " + str(r) + ", " + str(u) + ", " + str(p) + ", " + str(e) + "\n") def WENO(UU, Um, Up): BETA = [0.] * 3 ALPHA = [0.] * 3 eps = 1.0e-6 if (UU[2]-UU[1])*(UU[3]-UU[2]) < 0.0: Um = UU[2] else: 43

Продолжение ПРИЛОЖЕНИЯ Б BETA[0] = (13./12.)*(UU[2]-2*UU[3]+UU[4])*(UU[2]2*UU[3]+UU[4])+0.25*(3*UU[2]-4*UU[3]+UU[4])*(3*UU[2]-4*UU[3]+UU[4]) BETA[1] = (13./12.)*(UU[1]-2*UU[2]+UU[3])*(UU[1]2*UU[2]+UU[3])+0.25*(UU[1]-UU[3])*(UU[1]-UU[3]) BETA[2] = (13./12.)*(UU[0]-2*UU[1]+UU[2])*(UU[0]2*UU[1]+UU[2])+0.25*(UU[0]-4*UU[1]+3*UU[2])*(UU[0]-4*UU[1]+3*UU[2]) ALPHA[0] = 0.3/((eps+BETA[0])*(eps+BETA[0])) ALPHA[1] = 0.6/((eps+BETA[1])*(eps+BETA[1])) ALPHA[2] = 0.1/((eps+BETA[2])*(eps+BETA[2])) Um = (ALPHA[0]*(2*UU[2]+5*UU[3]-UU[4])+ALPHA[1]*(UU[1]+5*UU[2]+2*UU[3]) + ALPHA[2]*(2*UU[0]7*UU[1]+11*UU[2]))/((ALPHA[0]+ALPHA[1]+ALPHA[2])*6) if (UU[3]-UU[2])*(UU[4]-UU[3]) < 0.0: Up = UU[3] else: BETA[0] = (13./12.)*(UU[3]-2*UU[4]+UU[5])*(UU[3]2*UU[4]+UU[5])+0.25*(3*UU[3]-4*UU[4]+UU[5])*(3*UU[3]-4*UU[4]+UU[5]) BETA[1] = (13./12.)*(UU[2]-2*UU[3]+UU[4])*(UU[2]2*UU[3]+UU[4])+0.25*(UU[2]-UU[4])*(UU[2]-UU[4]) BETA[2] = (13./12.)*(UU[1]-2*UU[2]+UU[3])*(UU[1]2*UU[2]+UU[3])+0.25*(UU[1]-4*UU[2]+3*UU[3])*(UU[1]-4*UU[2]+3*UU[3]) ALPHA[0] = 0.1/((eps+BETA[0])*(eps+BETA[0])) ALPHA[1] = 0.6/((eps+BETA[1])*(eps+BETA[1])) ALPHA[2] = 0.3/((eps+BETA[2])*(eps+BETA[2])) Up = (ALPHA[0]*(11*UU[3]7*UU[4]+2*UU[5])+ALPHA[1]*(2*UU[2]+5*UU[3]-UU[4]) + ALPHA[2]*(UU[1]+5*UU[2]+2*UU[3]) )/((ALPHA[0]+ALPHA[1]+ALPHA[2])*6) return Um, Up def F_HLLC_U(UK, FK, SK, SS, PK, RK, VK): return (((SS)*((SK)*(UK)-(FK)) + (SK)*( (PK)+(RK)*((SK)(VK))*((SS)-(VK)) )) / ((SK)-(SS))) def F_HLLC_R(UK, FK, SK, SS, PK, RK, VK): return (((SS)*((SK)*(UK)-(FK))) / ((SK)-(SS))) def F_HLLC_E(UK, FK, SK, SS, PK, RK, VK): return (((SS)*((SK)*(UK)-(FK)) + (SK)*( (PK)+(RK)*((SK)(VK))*((SS)-(VK)) )*(SS)) / ((SK)-(SS))) def calc_flux_hllc(rl, ul, pl, gaml, rr, ur, pr, gamr, qr, qu, qe): sl, sr, p_star, s_star, p_pvrs, ql_, qr_, tmp = 0., 0., 0., 0., 0. , 0., 0., 0. czl = math.sqrt(pl*gaml/rl) czr = math.sqrt(pr*gamr/rr) el = pl/(rl*(gaml-1.)) 44

Продолжение ПРИЛОЖЕНИЯ Б er = pr/(rr*(gamr-1.)) e_tot_l = el+0.5*(ul*ul) e_tot_r = er+0.5*(ur*ur) p_pvrs = 0.5*(pl+pr)-0.5*(ur-ul)*0.25*(rl+rr)*(czl+czr) p_star = p_pvrs if p_pvrs > 0. else 0. ql_ = 1 if p_star <= pl else math.sqrt(1.+(gaml+1.)*(p_star/pl1.)/(2.*gaml)) qr_ = 1 if p_star <= pr else math.sqrt(1.+(gamr+1.)*(p_star/pr1.)/(2.*gamr)) sl = ul-czl*ql_ sr = ur+czr*qr_ if sl > sr: tmp = sl sl = sr sr = tmp s_star = pr-pl s_star += rl*ul*(sl-ul) s_star -= rr*ur*(sr-ur) s_star /= (rl*(sl-ul)-rr*(sr-ur)) if s_star < sl: s_star = sl if s_star > sr: s_star = sr if not ((sl <= s_star) and (s_star <= sr)): print("HLLC: inequaluty SL <= S* <= SR is FALSE.\n") sys.exit() if sl >= 0.: qu = rl*ul*ul+pl qe = (rl*e_tot_l+pl)*ul qr = rl*ul elif sr <= 0.: qu = rr*ur*ur+pr qe = (rr*e_tot_r+pr)*ur qr = rr*ur else: if s_star >= 0: # UK, FK, SK, SS, PK, RK, VK 45

Окончание ПРИЛОЖЕНИЯ Б qu = F_HLLC_U(rl * ul, rl * ul * ul + pl, sl, s_star, pl, rl, ul) # UK, FK, SK, SS, PK, RK, VK qe = F_HLLC_E(rl * e_tot_l, (rl * e_tot_l + pl)*ul, sl, s_ star, pl, rl, ul) # UK, FK, SK, SS, PK, RK, VK qr = F_HLLC_R(rl, rl * ul, sl, s_star, pl, rl, ul) else: # UK, FK, SK, SS, PK, RK, VK qu = F_HLLC_U(rr * ur, rr * ur * ur + pr, sr, s_star, pr, rr, ur) # UK, FK, SK, SS, PK, RK, VK qe = F_HLLC_E(rr * e_tot_r, (rr * e_tot_r + pr)*ur, sr, s_ star, pr, rr, ur) # UK, FK, SK, SS, PK, RK, VK qr = F_HLLC_R(rr, rr * ur, sr, s_star, pr, rr, ur) return qr, qu, qe main() 46

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв