Распараллеливание методом декомпозиции области численного решения уравнения диффузии

Дорожкина Алла Алексеевна

Распараллеливание методом декомпозиции области численного решения уравнения диффузии

В данной выпускной бакалаврской работе рассмотрена следующая тема: Распараллеливание области численного решения уравнения диффузии методом декомпозиции. Цель работы: реализовать метод, который позволяет строить решение такого дифференциального уравнения быстро и эффективно. В работе осуществлен метод декомпозиции области, в котором используется быстрый способ решения системы уравнений. Проведено сравнение быстрого метода с методом Гаусса. Код программы распараллелен на несколько потоков. Представлены временные результаты и графики решения дифференциального уравнения Пуассона с разным количество точек для прямой формы и параллельного хода программы. Данная работа содержит 13 рисунков, 1 таблицу. Ссылается на 10 различных источников литературы и содержит 24 страницы. Работа реализована посредством пакета прикладных программ для решения задач технических вычислений и одноимённым языка программирования Matlab R2014a.

Математика

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d364f5f1be77c40d58c33

UUID: 4a7c8fc3-8134-47ca-99b9-33a911406962

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

Опубликовано: почти 8 лет назад

Просмотры: 80

Дорожкина Алла Алексеевна

Источник: Санкт-Петербургский государственный университет

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

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

Скачать - 1514906 bytes

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

Санкт-Петербургский Государственный Университет Математико-механический факультет Кафедра параллельных алгоритмов Дорожкина Алла Алексеевна Распараллеливание методом декомпозиции области численного решения уравнения диффузии Бакалаврская работа Научный руководитель: Профессор Корнеев Вадим Глебович Рецензент: Профессор Демьянович Юрий Казимирович Санкт-Петербург 2016

SAINT-PETERSBURG STATE UNIVERSITY Department of Parallel Algorithms Dorozhkina Alla Parallel solution of diffusion equation by domain decomposition method Bachelor’s Thesis Scientific supervisor: Professor Vadim Glebovich Korneev Reviewer: Professor Yuri Kazimirovich Dem'yanovich Saint-Petersburg 2016

Оглавление Введение......................................................................................................... 4 Постановка задачи ......................................................................................... 7 Решение .......................................................................................................... 8 Программная реализация ............................................................................ 10 Результаты .................................................................................................... 15 Заключение................................................................................................... 23 Список литературы ....................................................................................... 24 3

Введение Уравнение диффузии является частным видом дифференциального уравнения в частных производных. Решая уравнение диффузии, необходимо найти зависимость концентрации вещества (или объектов) от пространственных координат и времени. Если речь идет об общем случае, то данный коэффициент зависит от времени и пространственных координат. В общем виде уравнение записывается следующим образом: 𝜕𝜑(𝑟, 𝑡) = ∇ ∙ [𝐷(𝜑, 𝑟)∇𝜑(𝑟, 𝑡)], 𝜕𝑡 где φ(r, t) — плотность диффундирующего вещества в точке r и во время t и D(φ, r) — обобщённый диффузионный коэффициент для плотности φ в точке r. ∇ - оператор набла. Если коэффициент диффузии зависит от плотности — уравнение нелинейно, в противном случае — линейно. Если D постоянное, то уравнение диффузии сводится к уравнению теплопроводности (линейному дифференциальному уравнению). Уравнение диффузии было открыто в 1855 году известным немецким физиком и филологом Адольфом Фиком (1829 – 1901) по аналогии с уравнением Ж. Фурье для потока тепла [1]. Я. М. Гельфер в книге «История и методология термодинамики и статистической физики» (1981) отмечает: «Явление диффузии было впервые исследовано вюрцбургским ученым А.Фиком на примере соляных растворов. Фик путем тщательных исследований показал, что свободная диффузия соляных растворов происходит по законам, совершенно аналогичным законам распространения тепла в твердых телах». Уравнение диффузии бывает стационарным и нестационарным. Первое представляет собой параболическое 4

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

нестационарном уравнении. Общий вид стационарного уравнения таким образом: При , не зависящем от , стационарное уравнение диффузии становится уравнением Пуассона (данное уравнение неоднородное, описывает электростатическое поле, стационарное поле температуры, поле давления, поле потенциала скорости в гидродинамике) или уравнение Лапласа (однородное, т е при ) Численное решение уравнения Пуассона является важным элементом многих задач вычислительной физики. Так, например, при решении уравнений гидродинамики несжимаемой жидкости решение уравнения Пуассона требует значительных вычислительных усилий, особенно для расчетов с большим числом ячеек (от сотен тысяч и выше). [2] Метод декомпозиции является одни из методов быстрого решения данного уравнения. [3] 6

Постановка задачи Необходимо решить уравнение Пуассона – для области Ω = Ω1 ∪ Ω2, являющейся объединением двух подобластей с непустым пересечением (Рис.1). На всей границе δΩ устанавливается однородное краевое условие Дирихле u|δΩ= 0. Провести распараллеливание полученного решения. Сравнить время выполнения прямой формы и параллельной. Рис. 1. Расчетная область 7

Решение Проведем построение триангуляций области Ω = Ω1 ∪ Ω2 [4], обозначаемые Th посредством таких четырехугольных сеток, что сетка на каждом участке области Ωk топологически эквивалентна квадратной сетке с шагом h для прямоугольника, который обозначим через Ω*k. За сеточный параметр h принимаем среднее расстояние между двумя соседними узлами расчетной сетки МКЭ [5] с четырехугольными или треугольными конечными элементами первого порядка, имеющими прямолинейные стороны. Пусть KhUh = Fh – система алгебраических уравнений МКЭ для рассматриваемой задачи и Kh,k – блок на диагонали матрицы Kh , отвечающий задаче Дирихле на Ωk. Предобуславливатель [6] для матрицы Kh обозначаемый здесь Bh определяется через его обратный по формуле где Bh,k – прпедобуславливатель для блока Kh,k. Псевдообратные матрицы для матриц Bh,k имеют вид: В качестве предобуславливателей для области декомпозиции выбирается матрица Bh,k = Δh,k , где Δh,k пятиточечная матрица жесткости МКЭ для уравнения Лапласа в прямоугольнике на регулярной квадратной сетке шага h. Данные матрицы спектрально эквивалентны матрицам Kh,k . Для решения задачи применяется предобусловленный метод сопряженных градиентов (PCG) [7] c предобуславливателем Bh. 8

Уравнение для итераций: Вектор решения уравнения Пуассона состоит из двух частей, соответствующих I и II области. Индексы значений для первой и второй области перекрываются, так как перекрываются сами области. Получаем систему уравнений для нахождения решений для первой и второй подобласти итерационным методом. Где VI = - KI UI + FI и VII = - KII UII + FII - невязки для I и II подобластей Введем обозначения: 𝑊𝐼 = 𝐵−1 (−𝐾𝐼 𝑈 + 𝐹𝐼 ) 𝑊𝐼𝐼 = 𝐵−1 (−𝐾𝐼𝐼 𝑈 + 𝐹𝐼𝐼 ) Невязки будем находить из следующих уравнений: 𝐵𝐼 𝑊𝐼 = 𝑉𝐼 𝐵𝐼𝐼 𝑊𝐼𝐼 = 𝑉𝐼𝐼 Таким образом система уравнений принимает вид: 9

Программная реализация Программная реализация решения выполнена с использованием пакета MATLAB [8]. Структура реализации выбрана следующая: имеется главная функция main, которая отвечает за построение области и вывод решений. Все вычисления инкапсулированы во второй функции – solver, которая собственно и производит решение уравнения Пуассона. Главная функция: main, передает в функцию solver область и получает решения, которые отображает в виде графиков. В качестве исходных данных в функции main задается область, на которой решается задача и строится регулярная сетка для данной области. После этого задается область, соответствующая исходной, состоящая из секторов, необходимая для вычисления матриц предобуславливателей. После задания области функция main вызывает функцию solver передавая области в качестве параметров. Функция solver действует по описанному далее алгоритму. Происходит запуск цикла, который на каждом шаге увеличивает число разбиений области. Внутри цикла имеется: 1) Счетчик времени обработки одной итерации цикла, т.е. подсчета времени потраченного на построение решения. 2) Блок, который позволяет программе выполнять расчет как в параллельной форме, так и прямой в зависимости от значения параметра is_parallel if is_parallel parfor i = 1:3 if i == 1 dim = [1 2 3]; elseif i == 2 dim = [3 4]; else dim = [1 2 3 4]; 10

end [ks{i}, fs{i}] = assempde(b, p, e, t, 1, 0, 1, [], dim); end else for i = 1:3 if i == 1 dim = [1 2 3]; elseif i == 2 dim = [3 4]; else dim = [1 2 3 4]; end [ks{i}, fs{i}] = assempde(b, p, e, t, 1, 0, 1, [], dim); end end В данных циклах рассматривается пересечение участков области 1, 2, 3 и участков 3 и 4 (см. рис. 2). Выбор участков можно изменить, задав в качестве подобласти Ω1, например, участки 1 и 2, а в качестве Ω2 участки под номерами 2, 3 и 4. Рис. 2. 3) Вычисление матрицы жесткости и вектора FI для уравнения KI UI = FI, которое решается для первой подобласти Ω1, выделенной красным цветом на рис. 3 11

Рис. 3 4) Вычисление матрицы жесткости и вектора FII для уравнения KII UII = FII, которое решается для второй подобласти Ω2, выделенной синим цветом на рис. 4. Рис. 4 5) Вычисление матрицы жесткости для предобуславливателя BI для области, представленной на рис. 3 12

6) Вычисление матрицу жесткости для предобуславливателя BII для области, представленной на рис. 4 7) На каждом разбиении сетки строится итерационный процесс, который осуществляет решение уравнений. Задается первое приближение решения U для каждой подобласти. На каждой итерации происходит расчет нового значения U. Цикл работает до тех пор, пока норма разности решений на текущем и предыдущем шаге не станет меньше заданной погрешности. По достижении нормой этого числа – осуществляется выход из цикла и запоминается решение уравнения. Использованы среднеквадратичная норма и норма, равная max|u(i)|. В данном случае используется погрешность 1e-4. Содержание цикла приведено ниже: while d > 1e-4 u_old = u; v1 = -KI * u1 + FI; w1 = BI \ v1; u1 = u1 + 1 * w1; u(i_up) = u1; u2 = u(i_down); v2 = -KII * u2 + FII; w2 = BII \ v2; u2 = u2 + 1 * w2; u(i_down) = u2; u1 = u(i_up); d = norm(u - u_old); iters(iter) = iters(iter) + 1; end После завершения работы вышеуказанного цикла в solver функция main строит следующие графики: 1) Набор решений исходного уравнения для различного количества точек и итераций 2) Зависимости числа итераций от количества вершин 13

3) Зависимость времени, потраченного на решение с разным количеством вершин для параллельной формы и прямой Следующие важные моменты которые уменьшают время работы программы. Нахождение w из системы Bw = v, где B - предобуславливатель, а v – невязка лучше осуществлять с помощью быстрых методов. В настоящей реализации для нахождения используется двухслойная схема [9]. Для решения можно также применить многосеточный метод [10], который является более быстрым, но он сложнее в реализации. Результаты использование в данном случае, например, метода Гаусса, будут значительно хуже вышеприведенных методов. Так же важным моментом, который существенно увеличивает производительность методов, является тот факт, что надо избегать заполнения и сохранения нулевых элементов матриц, а сохранять и умножать только ненулевые элементы. Например, матрица B – матрица, получаемая из дискретизации оператора Лапласа и в ней только 5 отличных от нуля диагоналей. Поэтому двухслойная схема для такой матрицы будет работать быстро, если сохранить только ненулевые элементы (или эти 5 диагоналей). 14

Результаты Ниже представлены графики решения уравнения для разного количества точек (рис. 5-12) Рис. 5 Рис. 6 15

Рис. 7 Рис. 8 16

Рис. 9 Рис. 10 17

Рис. 11 Рис 12. На графиках, представленных ниже (рис. 13,14), показана зависимость количества вершин от времени в секундах. Пустыми кружочками обозначены точки для параллельного метода. Первоначальная инициализация занимает длительное время, в течение которого процессор создает потоки, что вносит дополнительные накладные расходы. В случае небольшого числа точек, 18

данные расходы занимают значительную часть от всего времени выполнения алгоритма. Поэтому при небольшом количестве точек расчетное время выполнения при использовании прямого метода несколько меньше чем при использовании параллельной формы. Для большого же числа точек выгоднее считать параллельной формой. Рис 13. 19

Рис. 14 Рассмотрев графики, можно понять, что время параллельной формы меньше, в данном случае незначительно, так как результаты представлены для работы на обыкновенном компьютере, имеющем только 2 ядра, следовательно, процессы распараллелены только на 2 потока. В таблице 1 показана зависимость времени выполнения от количества итераций для параллельного и прямого метода. Время решения уравнения методом декомпозиции в обоих случаях достаточно мало. Для 8-ого построения, т. е. в том случае, когда количество точек равно 263425, время работы прямой формы равно 3.12 секундам, для параллельной 2.7 20

Итер. Прямая форма (с.) Параллельная форма (с.) 1 0 0.05 2 0.01 0.05 3 0.01 0.06 4 0 0.07 5 0.01 0.05 6 0.11 0.08 7 0.67 0.53 8 3.12 2.74 Таблица 1 Данные результаты можно сравнить с результатами, полученными с помощью метода Гаусса, который находит w, решая систему Bw = v, в которой B – предобуславливатель, а v – невязка. Уже на 5-ой итерации построение решения занимает почти 200 секунд. Зависимость времени от количества вершин для метода Гаусса показана на графике на рис. 15. 21

Рис. 15 22

Заключение В ходе настоящей работы выполнена задача нахождения численного решения уравнения диффузии. Выявлено, что распараллеливание области методом декомпозиции является эффективным, что подтверждается полученными графиками решений и полученными результатами зависимости времени построения решения от количества точек. В соответствии с результатами, время, затраченное на решение уравнения Пуассона для достаточно большого количества точек (например, 263 тысячи, что является достаточно точным решением) является малым (для данного примера 2,7 секунды при параллельной форме и 3.1 секунды для прямой формы). Проведено сравнение применения метода декомпозиции с использованием быстрого метода решения системы уравнений и с использованием метода Гаусса. Второй способ менее эффективный. Для 4257 точек метод, разработанный в данной работе выполняет построение решения за 0.1 секунду при прямой форме работы программы и 0.5 при параллельной форме (разница с прямой формой такая, т к изначально происходит долгая инициализация всех данных, нужных для работы программы). Для метода декомпозиции в совокупности с методом Гаусса же время работы равняется 200 секундам. Таким образом, применение метода декомпозиции с быстрым решением системы, отсутствием хранения неиспользуемых ячеек матриц и параллельным выполнением функций является эффективным в задачи нахождения области численного решения уравнения диффузии. 23

Список литературы 1. A. Fick, Ueber Diffusion, Pogg. Ann. Phys. Chem. 170 (4. Reihe 94), 59-86 (1855). 2. A Tutorial on Elliptic PDE Solvers and their Parallelization Craig C. Douglas, Gundolf Haase, Ulrich Langer 3. V. Korneev and U. Langer. Domain Decomposition Methods and Preconditioning. In: Encyclopedia of Computational Mechanics,V.1. E. Stein, R. de Borst and Th.J.R. Hudges eds. 2004 John Wiley \& Sons,Ltd. 4. А. А. Клячин, Равномерная триангуляция плоских областей, 1-6 (2011) 5. В.Г. Корнеев, С. Е. Енсен Эффективное предобуславливание методом декомпозиции области для p-версии с иерархическим базvисом. I // Изв. вузов. Математика. - 1999. - №5. - C. 37-56. 6. В. Г. Корнеев, “Схемы метода конечных элементов высоких порядков точности для эллиптических уравнений второго порядка в трехмерных областях. 1979 год 7. Preconditioned conjugate gradients method. http://www.mathworks.com/help/matlab/ref/pcg.html. 2016 8. MathWorks – Makers of MATLAB and Simulink http://www.mathworks.com/ 9. А.А. Самарский, Е.С. Николаев Методы решения сеточных уравнений Глава 6, Москва «НАУКА» Главная редакция физико-математической литературы 1978 год 10. В. В. Мареев, Е. Н. Станкова, Многосеточные методы. Введение в стандартные методы.: Издательство СПбГУ, 2012 24

Рецензии:

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

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

Отзывы:

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