авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 2 | 3 ||

«ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ УЧРЕЖДЕНИЕ НАУКИ ИНСТИТУТ ВЫЧИСЛИТЕЛЬНОГО МОДЕЛИРОВАНИЯ СИБИРСКОГО ОТДЕЛЕНИЯ РОССИЙСКОЙ АКАДЕМИИ НАУК ...»

-- [ Страница 4 ] --

внутренняя граница, диэлектрик, проводник. Происходит проверка на непротиворечивость введенных данных.

3. Заполняются массивы координат x, y, z узлов всех граничных элемен тов сетки, точек, сегментов, четырехугольников.

4. Рассчитываются все недостающие сеточные значения массивов x, y, z внутри шестигранников. В случае использования многосеточного мето да, данные вводятся для основной самой мелкой сетки.

5. Проводится проверка невырожденности сетки.

6. Задаются массивы C - коэффициент дифференциального уравнения, например, проводимость, Q - вычисляется правая часть уравнения 1.6, например, распределение плотности источника заряда в пространстве, F - строится начальное приближенное значение искомой функции, на пример, потенциала. Последнее, в частности, нужно для задания гра ничных значений в случае задачи Дирихле. Задать массивы можно, считав их из файла, или вычисляя значение аналитической функции от координат данного узла сетки. В случае использования многосеточ ного метода данные вводятся для основной, самой мелкой, сетки.

7. - вычисляются коэффициенты B вариационно-разностной схемы с ис пользованием массивов x, y, z и C для основной сетки. Для более круп ных сеток, если это необходимо, происходит вычисление по уже насчи танным значениям для более мелких сеток. Проводится проверка на не вырожденность элементов области.

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

9. - полученный результат сравнивается с аналитическим решением, если таковое известно.

10. - визуализация полученного решения.

11. - дополнительная возможность получить более точное решение в выде ленной подобласти. В области решения задачи строится специальный шестигранник. Интерполяцией уже имеющихся значений заполняются соответствующие массивы x, y, z, C, F, Q. Повторяются шаги 6-9 для этого шестигранника, чтобы получить более точное решение в выделен ной подобласти.

Поиск решения системы линейных алгебраических уравнений, соответ ствующих вариационно-разностной схеме, выполняется с помощью одного из двух итерационных методов: метода последовательной верхней релаксации или многосеточного метода, которые описаны ниже в разделах A.2.3. и A.2.2.

Отметим, что при значении релаксационного параметра = 1 метод последо вательной верхней релаксации становится методом Зейделя, который также называется методом Гаусса-Зейделя.

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

1. Самая мелкая сетка ks = 1 задается текущей.

2. Вычисляется невязка H на текущей сетке с помощью одной итерации метода последовательной верхней релаксации с нулевым шагом.

3. Невязка H переносится на более крупную сетку ks + 1 в массив правой части Q.

4. Более крупная сетка задается текущей ks = ks + 1.

5. Операции 2-4 повторяются, пока текущей сеткой не становится самая крупная ks = ksm 6. Начиная с нулевых значений искомой функции F, осуществляются ите рации метода последовательной верхней релаксации на текущей сетке, чтобы получить приближенные значения элементов массива F. Коли чество итераций задается пользователем.

7. Полученные значения F переносятся с текущей сетки на более мелкую ks 1 с добавлением к уже имеющимся значениям.

8. Более мелкая сетка задается текущей ks = ks 9. Операции 6-8 повторяются, пока текущей сеткой не становится основная сетка ks = 10. Осуществляются итерации метода последовательной верхней релакса ции на основной сетке, чтобы уточнить приближенные значения F. Ко личество итераций задается пользователем.

11. Операции 2-10 повторяются заданное количество раз или пока норма невязки не достигнет требуемого значения.

A.2.3. Метод последовательной верхней релаксации Описание и обоснование метода последовательной верхней релаксации, который иногда называется методом верхней релаксации, можно найти в учебнике [40]. Это один из простейших итерационных методов решения си стем линейных алгебраических уравнений. Мы используем его при решении системы (2.38) и аналогичных систем для крупных сеток при реализации многосеточного метода, поскольку сам по себе он дает малую скорость схо димости. Приведем необходимые формулы на примере системы I (A.1) Lij Fj = qi, i = 1,..., I.

j= Каждая итерация метода последовательной верхней релаксации состоит в последовательной замене приближенных значений Fi на новые Fi за счет ре шения отдельно iтого уравнения. При этом используются "старые"значения Fj, j = i + 1,..., I, и "новые"значения Fj, j = 1,..., i 1, которые к моменту вычисления Fi уже получены.

Если просто разрешать iтое уравнение относительно Fi i1 I (A.2) Fi = qi Lij Fj Lij Fj, Lii j=1 j=i+ такие итерации называются методом Зейделя.

Перепишем эту формулу в виде i1 I (A.3) Fi Fi = qi Lij Fj Lij Fi Lij Fj, Lii j=1 j=i+ который можно трактовать как изменение Fi с целью обратить в нуль невязку iтого уравнения системы (A.1). Невязка iтого уравнения при полученных до его разрешения значениях неизвестных и есть выражение в скобках в (A.3).

Обозначим его как Hi.

Зачастую сходимость итерационного процесса можно ускорить, если из менять значение Fi на иную величину Hi (A.4) Fi Fi =.

Lii В учебнике [40] доказано, что такой итерационный процесс сходится при произвольном значении итерационного параметра из интервала 0 для симметричных положительно определенных матриц L, к числу которых принадлежат матрицы L системы (2.38) и аналогичных систем для укруп ненных сеток.

При 1 имеем метод нижней релаксации, при = 1 – метод Зейделя, при 1 – метод верхней релаксации, который мы и используем.

Отметим, что формулы (A.1-A.4) зачастую записывают более компакт но, представляя матрицу L в виде суммы диагональной матрицы с элемен тами Lii, матрицы с элементами L, лежащими ниже ее главной диагонали, и третьей матрицы с элементами L, лежащими выше ее главной диагонали.

Приведенные здесь детальные формулы представляются более удобными для программирования.

В случае использования метода последовательной верхней релаксации как основного метода решения задачи, что мы делаем лишь для демонстрации эффективности многосеточного метода, вычисления проводятся только на основной сетке. Итерационный параметр задается пользователем с учетом ограничений 0 2, которые гарантируют сходимость итераций. При ре шении уравнения Пуассона на сетках, близких к равномерным прямоуголь ным, оптимальное = 1.3. При использовании существенно неравномерных сеток с вытянутыми ячейками и при существенно переменных коэффициен тах происходит ухудшение обусловленности матрицы системы линейных ал гебраических уравнений, соответствующих вариационно-разностной схеме, и тогда оптимальное возрастает почти до = 2.

Итерации метода последовательной верхней релаксации с учетом специ фики разбиения сетки на элементы организованы следующим образом:

1. Установить равным нулю текущее значение нормы невязки hhiter для выполняемой итерации номер iter.

2. Для каждого узла сетки i, кроме граничных i = 1, nn, nn + 1, nn,..., nn nm,..., nn nm lm, внутри каждого шестигранника P ar вычислить величину невязки hi, используя текущее значение в узле Fi и значения Fj в соседних узлах j = i ± 1, i ± (nn ± 1),... с использовани ем коэффициентов Bi,j, соответствующих данной паре узлов, и узловое значение правой части уравнения Qi.

3. Добавить величину значения невязки hi каждого узла к норме невязки hhiter = hhiter + hi. Для всех итераций iter кроме последней iter = itermax вычислить новое значение Fi, используя значение невязки hi, коэффициент Bi,i и параметр.

4. Провести вычисления невязки hi для узлов, лежащих на гранях i = 1, nn, nn + 1, 2 nn,..., nn nm,..., nn nm lm каждого шестигран ника P ar, учитывая только значения Fj соседних узлов j = i ± 1, i ± (nn ± 1),... внутри текущего шестигранника с соответствующим им ко эффициентом массива Bi,j. С этого шага и по 10 учитывается тип ni граничного условия.

5. Перенести рассчитанные граничные значения невязки hi с граней ше стигранников P ar на соответствующие им четырехугольники Rect, сум мируя значения невязки hi = hj + hk для общих узлов j и k на гранях соседствующих шестигранников Rect(i) P ar(j), Rect(i) P ar(k).

6. Для каждого узла сетки i, кроме граничных i = 1, nn, nn+1,..., nnnm, внутри каждого четырехугольника Rect вычислить величину невяз ки, используя текущее значение Fi, значения Fj в соседних узлах j = i ± 1, i ± (nn ± 1), с использованием коэффициентов Bi,j, полу ченный изнутри шестигранников вклад в невязку hi и значение правой части уравнения Qi.

7. Добавить значение невязки hi в норму невязки hhiter = hhiter + hi. Для всех итераций iter кроме последней iter = itermax вычислить новое зна чение Fi, используя значение невязки hi, коэффициент Bi,i и параметр.

8. Провести вычисления вклад в невязки hi изнутри четырехугольника для узлов i = 1, nn, nn + 1, 2 nn,..., nn nm, лежащих на ребрах четырехугольников Rect, учитывая только значения узлов Fj внутри текущего четырехугольника j = i ± 1, i ± (nn ± 1), с соответствующим им коэффициентом Bi,j.

9. Перенести рассчитанные граничные значения невязки hi с ребер че тырехугольников Rect на соответствующие им сегменты Seg, сумми руя значения невязки hi = hj + hk + hl для общих узлов Seg(i) Rect(j), Seg(i) Rect(j), Seg(k) Rect(l).

10. Для каждого узла сетки i внутри каждого сегмента Seg, кроме его концов i = 1, nn, вычислить невязку hi, используя текущее значение Fi, значения Fj в соседних узлах j = i ± 1 и коэффициента Bi,j для каждого соседнего узла, полученный изнутри четырехугольников вклад в невязку hi и значение правой части уравнения Qi.

11. Добавить все значения невязки hi к норме невязки hhiter = hhiter + hi. Для всех итераций iter кроме последней iter = itermax вычислить новое значение Fi, используя значение невязки hi, коэффициент Bi,i и параметр.

12. Провести вычисления вклад изнутри сегмента Seg в невязку hi для граничных узлов i = 1, nn, учитывая только значения Fj в соседних узлах j = i 1, i + 1 внутри текущего сегмента Seg с соответствующим им коэффициентами Bi,j.

13. Перенести рассчитанные граничные значения невязки hi из концов сег ментов Seg в соответствующие им точки P oint, суммируя значения невязки hi = hj + hk + hl для общих узлов P oint(i) Seg(j), P oint(i) Seg(k), P oint(i) Seg(l).

14. Для каждой точки P oint вычислить невязку hP oint, используя текущее значение узла FP oint, полученный изнутри сегментов вклад в невязку hP oint и значение правой части уравнения QP oint.

15. Добавить значение невязки hP oint каждой точки к норме невязки hhiter = hhiter + hP oint. Для всех итераций iter кроме последней iter = itermax вычислить новое значение FP oint, используя значение невязки hP oint, коэффициент BP oint,P oint и параметр.

16. Перенести значение FP oint из каждой точки P oint в соответствующие этой точке граничные узлы сегментов Fi = Fj, P oint(i) Seg(j).

17. Перенести значение Fi с каждого сегмента Seg в соответствующие этим сегментам граничные узлы четырехугольников Rect Fi = Fj, Seg(i) Rect(j).

18. Перенести значение Fi каждого четырехугольника Rect в соответству ющие этим четырехугольникам граничные узлы шестигранников P ar Fi = Fj, Rect(i) P ar(j).

19. Повторить шаги 1-18 нужное количество iter раз или пока значение нормы невязки hhiter для текущей итерации iter не уменьшится до необ ходимой величины hhmin.

A.3. Алгоритм работы параллельной программы Как было сказано выше, используемое в нашем комплексе программ постро ение сетки позволяет использовать при расчетах параллельные вычисления, так как вычисления внутри каждого шестигранника могут выполняется неза висимо. Но на дополнительных границах, построенных нами внутри области при ее разделении на криволинейные шестигранники, необходим перенос рас считанных значений на границу и вычисление значений на самой границе, а также обратный перенос рассчитанных граничных значений в соответствую щие узлы шестигранников.

Решение этой задачи было выполнено также на языке Fortran с помощью библиотеки MPI.

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

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

A.3.1. Подготовка данных для параллельной программы Подготовка данных для параллельной части выполняется на персональ ном компьютере пользователя. Алгоритм работы показан на Рис. A.2.1. и состоит из следующих операций:

1. Инициализация программы: задаются размеры массивов для хранения данных и под них выделяется память.

2. Задается область: определяется количество всех элементов сетки, зада ются типы прилегания элементов, типы граничных условий для всех граничных элементов: внутренняя граница, условие Неймана (изоля тор), условие Дирихле (проводник). Проверяется непротиворечивость введенных данных.

3. В файл обмена сохраняется общая информация об области и ее элемен тах.

4. Заполняются массивы координат узлов сетки x, y, z, для следующих элементов сетки: точки, сегменты, четырехугольники. Также проводит ся проверка непротиворечивости этих данных. В случае использования многосеточного метода данные вводятся для основной самой мелкой сетки.

5. Значений x, y, z, полученные для четырехугольников, пересылаются на грани всех шестигранников.

6. Задаются массивы C - узловых значений коэффициента дифференци ального уравнения 1.6, например, проводимости, Q - правой части урав нения 1.6, например, плотность источника заряда в задаче электроста тики или плотность источника тока – в задаче электропроводности, F - начального значения неизвестной функции, например, потенциала (необходимо для задания граничных значений при условии Дирихле).

Задать массивы можно, считав их из файла, или используя аналити ческой функцию от координат данного узла сетки. В случае исполь зования многосеточного метода данные вводятся для основной, самой мелкой, сетки.

7. Сохранение значений массивов C, F, Q всех шестигранников.

В случае задания большого числа элементов и внутренних узлов элемен тов, пункты 6-7 выполняются последовательно для каждого шестигранника.

Так как сетка внутри шестигранников в нашем комплексе программ стро ится с помощью интерполяции координат узлов с граней, вместо полных мас сивов x, y, z можно передавать только значения на гранях каждого шести гранника, что заметно уменьшает объем данных, переносимых на суперком пьютер.

A.3.2. Multiproc-программа Многопроцессорная часть программного комплекса выполнена с исполь зованием библиотеки MPI. Использование данной библиотеки связано с вы сокой гибкостью получаемых решений, высокой распространенностью этой библиотеки и простотой использования полученной программы на большин стве суперкомпьютеров с разделенной памятью.

При запуске программы на суперкомпьютере задается число исполь зуемых процессоров. Внутри программы особенно выделяется только 0 процессор, который считается мастер-процессором, все остальные процессоры для программы одинаковы.

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

Общий алгоритм работы показан на Рис. A.2.1., часть II, и детализирован ниже. Уникальные части алгоритма для мастер-процессора помечены буквой М, а для остальных процессоров О.

1. Инициализировать программу, MPI. Определить мастер-процессор.

2M. Считать с файла данные для определения конфигурации области и эле ментов области. Разослать считанные данные на все остальные процес соры. Сохранить считанные данные в выходной файл.

2O. Получить от мастер-процессора данные о конфигурации области и эле ментов области.

3. Рассчитать распределение элементов по процессорам.

4M. Считать с файла значения массивов x, y, z, C, F, Q и разослать дан ные элементов сетки на соответствующие им процессоры, которые бу дут рассчитывать данные внутри элемента сетки. Сохранить считанные значения массивов x, y, z в выходной файл.

4O. Получить от мастер-процессора значения массивов x, y, z, C, F, Q, которые принадлежат элементам сетки, используемых при расчетах на этом процессоре.

5. Рассчитать все недостающие значения массивов x, y, z внутри шести гранников.

6. Вычислить коэффициенты B, используя данные из массивов x, y, z и C для основной сетки. Для более крупных сеток, если это необходи мо, происходит вычисление по уже насчитанным значениям для более мелких сеток.

7. Провести итерацию для уточнения решения задачи в шестиграннике.

8O. Отправить полученное решение внутри шестигранников на мастер процессор.

8M. Получить решение задачи от всех процессоров. Сохранить полученное решение в файл.

Многосеточный алгоритм поиска решения не отличается от однопроцес сорной версии.

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

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

1. В соответствии с распределением элементов сетки по процессорам со брать данные, необходимые для передачи на каждый процессор.

2. Асинхронно отправить собранные данные соответствующим этим дан ными процессорам.

3. Ждать получения данных от других процессоров.

4. Получить данные и поместить их в соответствующие элементы сетки.

Асинхронная передача была выбрана из-за невозможности равномерно распределить задачи на процессорах, а использование барьеров и синхронно го обмена данными будет медленнее и сложнее.

A.3.3. Обработка полученных данных Обработка решения, полученного Multiproc-программой, такое же как и полученного Single-программой:

1. Считать из файла данные для определения конфигурации области и элементов области.

2. Считать из файла значения массивов x, y, z для всех элементов области.

3. Считать из файла значения массива F.

4. Вычислить сеточные значения аналитического решения A.

5. Сравнить полученный результат F с аналитическим решением A.

6. Визуализировать решение F.

Приложение B Построение сетки в расчетной области Большинство методов решения краевых задач для дифференциальных урав нений в частных производных состоит в замене исходных задач на задачи для сеточных функций, то есть функций, значения которых определены в узлах некоторых сеток. Такие задачи называют дискретным. По существу, основным требованием к сетке является в некотором смысле равномерная по области расстановка узлов. Желательным свойством является регулярная нумерация узлов сетки. Это проще всего достичь, если область является па раллелепипедом. Достаточно выбрать одну из вершин, расставить на одном из выходящих из этой вершины сегментов nm узлов, считая концы, на втором выходящем сегмента mm узлов, на третьем - lm узлов. Сегменты выбираем в таком порядке, чтобы они образовывали правую систему координат, ана логично осям x,y,z декартовых координат. После этого достаточно провести плоскости через построенные узлы, параллельные граням параллелепипеда и взять в качестве сеточных узлов получившиеся nm mm lm точек пересече ния этих плоскостей. Если теперь деформировать область, перемещая узлы сетки, не меняя их нумерации, получатся области, называемые криволиней ными шестигранниками. При этом должны сохраняться такие свойства, как положительность объемов ячеек и т. п.

Построение такой сетки в произвольной области является сложной зада чей.

Полной альтернативой являются нерегулярные сетки, но их недостатком является сложность определения номеров соседних узлов.

6 4 4 1 Рис. B.1. Пример разбиения плоской области на 3 криволинейных четырехугольника B.1. Элементы сетки.

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

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

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

Для наглядности дублирования четырехугольники отодвинуты друг от друга. При этом дубликатом каждого четырехугольника могут быть одна или две грани шестигранников, сегмента – несколько ребер четырехугольников и шестигранников, точек – несколько вершин сегментов, четырехугольников и шестигранников.

2 Рис. B.2. Нумерация вершин в сегменте.

3 3 4 4 Рис. B.3. Нумерация ребер и вершин в четырехугольнике.

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

Все они пронумерованы в произвольном порядке, который фиксируется с по мощью внутренних переменных.

Сегмент, его вершины и направление показаны на рис. B.2.

Нумерация вершин и ребер, а также их направления, в четырехугольни ке изображены на рис. B.3. Для краткости направлением четырехугольника называем направление нормали к нему. Для четырехугольника на рис. B. положительным считается направление на нас. Для граней на рис. B.4 по ложительны направления наружу. Если четырехугольник направлен от нас, это, как обычно, отмечаем кружком с крестом, как это сделано на рис. B.4 и B.5.

Нумерация граней, ребер и вершин в шестиграннике, показана на рис.

B.4.

z, l z, l y, m y, m x, n x, n z, l y, m x, n Рис. B.4. Нумерация граней, ребер и вершин в шестиграннике.

Также на рисунке показаны положительные направления ребер и граней рассматриваемого шестигранника. Различные способы прилегания сегмента к ребру и четырехугольника к грани будут описаны далее.

3 2 2 1 4 1 1 2 1 3 2 3 4 2 1 4 4 3 2 3 3 4 3 1 4 1 2 4 Рис. B.5. Типы прилегания. Показаны 8 возможных положений прилегающего прямоугольника. Его нижняя часть отмечена двойной чертой. В случаях 1-4 прямоугольник "смотрит"на нас, в случаях 5-8 - от нас. Положение грани фиксировано. Они совпадают в случае 1.

B.2. Типы прилегания В общем случае при задании соответствия между гранью шестигранника и дублирующего ее четырехугольника их направления могут не совпадать.

При численной реализации граничных условий подразумевается, что четы рехугольник должен быть направлен наружу области. Например, у области, заданной одним шестигранником, все четырехугольники должны быть на правлены наружу, чтобы соответствовать общепринятому положительному направлению нормали на границе, но в случае когда 2 шестигранника со прикасаются, то четырехугольник, соответствующий 2-м соприкасающимся граням, будет прилегать к этим граням по разному, т.е. к одной грани по типу 1, а к другой по типу 5.

У четырехугольников может быть 8 типов прилегания, см. рис. B.5. Типы прилегания есть и у сегментов, но их всего 2: направления либо совпадают, либо противоположны.

B.3. Переменные Исследуемая область задается следующими основными параметрами:

kparm, krectm, ksegm и kpointm — количество шестигранников, четы рехугольников, сегментов и точек соответственно. Причем количество точек не превышает суммарного количества вершин шестигранников, так как точка может принадлежать не одному шестиграннику. Например у одного шести гранника: 6 граней, 12 сегментов и 8 вершин. А в области, состоящей из 2-х шестигранников с одной общей гранью: 11 четырехугольников, 20 сегментов и 12 вершин.

B.4. Массивы В массивах npar, nrect, nseg и npoint хранится описание всех элементов области, включающее в себя номера под-элементов: например, описание сег мента включает номера 2-х под-элементов – точек, являющихся началом и концом этого сегмента. Часть значений в эти массивы записывает пользова тель, а остальные насчитываются автоматически.

В дополнительных массивах npointseg и nsegrect описывается обрат ные зависимости, например, вершинами скольких и каких сегментов является данная точка.

Внутри программы есть и другие массивы, например, ntoprectrect и nedgepar. Они используются для задания соответствия номеров вершин или сторон четырехугольников номерам вершин и ребер той грани шестигранни ка, которой соответствует этот четырехугольник. Например, вершина 3 четы рехугольника, соответствующего грани 2 шестигранника, соответствует вер шине 7 этого же шестигранника, если грань и четырехугольник прилегают по типу 1.

Все дополнительные массивы насчитываются в программе. Они использу ются в том числе и для определения ошибок при задании элементов области.

B.5. Нумерация узлов сетки Строящаяся в каждом шестиграннике сетка, характеризуется следующими параметрами: nm, mm, lm – означающими количество узлов сетки на трех ребрах, выходящих из одной вершины, имеющей номер 4 на рис. B.4.

Удобно определить аналогичные общие параметры сразу для всех шести гранников, имея в виду, что через них вычисляются параметры отдельного шестигранника.

При использовании многосеточного метода количество используемых се ток задается параметром ksm. Первая сетка, самая подробная, каждая сле дующая имеет вдвое меньше ячеек по каждому направлению. Чтобы это было возможно, необходимо, чтобы параметры основной самой мелкой сетки: nm, mm, lm – были равны или кратны числу: 2ksm + 1 (например: 3 2ksm + 1).

Процедуры построения сетки эту специфику nm, mm, lm не использу ют. Существенно только совпадение таких параметров на общих граничных объектах.

B.6. Сеточные функции Сеточной функцией называется функция, значения которой определе ны во всех узлах сетки. Она может задаваться заранее, как например, xкоординаты узлов, или быть неизвестной, как например, узловые значе ния потенциала. В обоих случаях имеется в виду совокупность значений в узлах сетки.

Каждая сеточная функция хранится в отдельном одномерном массиве со специальной структурой, показанной на рис. B.6.

1-ая сетка (основная, все par все rect все seg все point самая мелкая) далее 2-ая сетка все par все rect все seg все point.

.

.

ksm сетка все par все rect все seg все point (самая крупная) Рис. B.6. Устройство массива значений сеточной функции.

.

Сначала хранится сеточная функция, определенная на первой, самой мел кой сетке. Значения на каждом элементе области хранятся цельно и последо вательно, а именно: сначала записаны значения во всех узлах первого шести гранника, затем второго, и т.д., после kparm-шестигранника находятся дан ные первого четырехугольника, затем второго и т.д. Некоторая информация дублируется. Например, значения в узлах четырехугольника совпадают со значениями на грани соответствующего шестигранника, или 2-х шестигран ников, если четырехугольник лежит внутри области.

Поскольку каждый элемент области содержит определенное число узлов сетки, то создаются массивы индекса (быстрого поиска): numpar, numrect, numseg и numpoint. Эти массивы содержат информацию о смещении каж дого элемента области относительно начала массива и информацию о коли честве узлов сетки в элементе: например, nm, mm, lm и istart для шестигран ника.

B.7. Алгоритм построения сетки.

1. Разделить область на шестигранники таким образом, чтобы каждое ребро и грань были ненулевой длинны и площади соответственно. При разбиении области также необходимо учитывать, что тип границы зада ется для целой грани. Вычислить количество полученных шестигранни ков, четырехугольников, сегментов и точек, сохранить их в переменных kparm, krectm, ksegm и kpointm соответственно.

2. Для каждого шестигранника заполнить в массиве npar соответствие его граней четырехугольникам в общей нумерации и тип их прилегания.

3. Для каждого четырехугольника заполнить в массиве nrect поля о типе границы для всего этого четырехугольника и поля нумерации ребер:

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

4. Для каждого сегмента заполнить в массиве nseg тип границы этого сегмента и соответствие его вершин точкам в общей нумерации.

5. Для каждой точки заполнить тип границы в массиве npoint.

6. Вызвать процедуру controln для проверки введенных данных и за полнения всех остальных полей массивов npar, nrect, nseg, npoint и данных о соответствии элементов сетки их подэлементам.

7. Задать количество ячеек сетки в каждом направлении nm, mm и lm у каждого шестигранника. Так как шестигранники соприкасаются друг с другом, то необходимо задавать одинаковое количество и расположение ячеек на соприкасающихся гранях. Также нужно вычислить количество узлов в каждом шестиграннике. Заполнить полученными данными мас сив numpar.

8. Вызвать процедуру controlm для проверки непротиворечивости вве денных данных и для заполнении массивов numrect, numseg и numpoint.

9. Задать координаты каждой точки в массивах x, y и z.

10. Перенести координаты точек на соответствующие им вершины сегмен тов процедурой shiftseg.

11. Задать координаты узлов на каждом сегменте. Если сегмент лежит на прямой, то для расположения узлов можно использовать специальную процедуру line. Если сегмент сетки лежит на дуге - можно использовать процедуру circle.

12. Перенести координаты узлов сегментов на соответствующие им ребра четырехугольников процедурой shiftrect.

13. Задать координаты узлов внутри каждого четырехугольника. Для рас положения узлов на плоскости можно использовать процедуру quadr.

14. Перенести координаты узлов четырехугольников на соответствующие им грани шестигранников процедурой shiftpar.

15. Задать координаты всех узлов внутри шестигранников. Для простых сеток можно воспользоваться процедурой par.

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

Тогда каждая грань естественным образом делится на четыре треугольника, и ячейка - на 24 тетраэдра. Значения кусочно-линейной функции в этих цен трах получаем некоторой интерполяцией узловых значений, обычно, как их среднее арифметическое.

Используя данный алгоритм для задания координат узлов сетки, можно описать большое количество реальных 3х-мерных областей. Главным услови ем является невырожденность элементов, то есть положительность объемов тетраэдров, площадей треугольников и длин ребер.

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

Приложение C Фактическая сходимость численного метода В настоящем приложении на примере трех задач Дирихле для уравнения электропроводности (1.17) демонстрируется фактическая сходимость числен ных решений к точным и скорость сходимости многосеточных итераций. За дачи отличаются видом области и распределением проводимости. Численные решения строятся на различных сетках. Используются различные нормы, чтобы характеризовать погрешности.

Для проверки созданного нами программного продукта используем мо дель атмосферы, представляющую собой область, окруженную идеальным проводником. Рассматриваем параллелепипед или часть шарового слоя. На нижней и боковых границах потенциал задаем равным нулю, а на верхней границе - произведению синусов. Проводимость внутри области задаем со гласно (4.1). Для таких задач несложно найти аналитические решения, что дает возможность наиболее полной проверки созданного численного метода.

В качестве меры разности между численным и точным решением исполь зуем значение максимума разности их значений в узлах сутки, которое обо значаем V. Используем также энергетическую норму этой разности. Для простоты ее вычисления предварительно выполняем линейную интерполя цию узловых значений точного решения внутри ячеек сетки, то есть строим его кусочно-линейное восполнение. Обозначаем ее (C.1) W = [Vh V, Vh V ], где обозначаемое квадратными скобками энергетическое скалярное произве дение определено формулой (2.14).


Чтобы характеризовать скорость сходимости итераций, вычисляем сумму модулей невязок всех решаемых линейных алгебраических уравнений, кото рую обозначаем H, а также корень из суммы их квадратов H2.

Чтобы полнее характеризовать скорость сходимости итераций, хорошо было бы знать точное решение системы линейных алгебраических уравне ний, которая при этом решается. Поскольку мы его не знаем, вместо него используем численное решение, получающееся после очень большого числа итераций, которое условно называем точным численным решением. В каче стве меры разности между численным решением после данной итерации и точным численным решением используем величину максимума разности их значений в узлах сутки, которую обозначаем Vh.

C.1. Задача в параллелепипеде.

Рассматриваемую область зададим в виде параллелепипеда x0 x0 zup с размерами горизонтальных сторон x0 = 10 000км и высотой zup = 80км.

Внутри области проводимость зададим функцией (4.1) с характерным для атмосферы параметром z0 = 6 км.

На нижней и боковых границах области зададим нулевой потенциал, а на верхней границе в виде функции:

V (x, y, z)|z=zup = sin(nx/x0) sin(my/x0 ), где n и m - целые числа, определяющие масштаб гармоник вдоль осей x и y соответственно.

Аналитическое решение уравнения электропроводности (1.17) имеет вид, лишь обозначениями отличающийся от (4.12):

exp (1 z) exp (2z) V (x, y, z) = sin(nx/x0) sin(my/x0 ), exp (1 zup ) exp (2zup ) 1/(2z0)2 + (n2 + m2 ) 2 /x2. (C.2) 1,2 = 1/(2z0) ± Для такой конфигурации области решим несколько задач с различными наборами параметров:

z0 =, n = 1, m = 1, z0 = 6 км, n = 1, m = 1, z0 = 6 км, n = 5, m = 5, (C.3) первый из которых соответствует постоянной проводимости, когда уравнение электропроводности становится уравнением Лапласа.

В качестве основных, самых мелких в многосеточном методе, используем равномерные прямоугольные сетки, имеющие следующее число шагов вдоль каждой декартовой оси:

(а) 16 16 16, (б) 32 32 32, (C.4) (в) 64 64 64.

На Рис. C.1 представлена работа многосеточного метода при решении сформулированной задачи Дирихле в области 400 км 400 км 80 км для уравнения электропроводности для каждого из трех наборов параметров (C.3) на основных сетках вида (C.4). В левой колонке показано убывание нор мы невязки H. В правой колонке показано убывание разности между реше нием после данной итерации и точным численным решением Vh (сплошная линия), а также разности между численным и аналитическими решениями V (пунктир). По горизонтали i - номер итерации.

Как видно по Рис. C.1, для всех задач итерации сходятся. При этом невяз ка достигает своего минимального уровня раньше, чем прекращает меняться решение. Точность чисел на современных компьютеров ограничена согласно стандарту IEEE Standard 754, что не позволяет достигать точности большей чем 1016. Этим объясняется отчетливо видный на графиках минимальный уровень погрешности около 1015. Он повышается с ростом числа узлов сетки 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа, (z) = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z б) (z) = e, n = m = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z в) (z) = e, n = m = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.1. Cходимость многосеточных итераций. Область 10 000 км 10 000 км 80 км. В левой колонке показано убывание нормы невязки H. В правой колонке показано убывание разности между решением после данной итерации и точным численным решением Vh (сплошная линия), а также разности между численным и аналитическими решениями V (пунктир). По горизонтали i - номер итерации.

просто из-за суммирования большего числа невязок, соответствующих каж дому узлу сетки.

При сравнении численного решения с аналитическим видно, что точность решения достигает некоторого предела, а все дальнейшие итерации изменя ют решение в некоторой области около точного решения, не приближаясь и не удаляясь от него. Поэтому большую часть операций можно было не про водить, но определить этот момент на практике затруднительно.

Также из сравнения с аналитическим решением видно, что при уменьше нии размера ячейки сетки (увеличении количества ячеек) в 8 раз (вдвое по каждому направлению), точность численного решения возрастает примерно вчетверо. Это соответствует теоретически известной квадратичной сходимо сти для метода конечных элементов, когда кусочно-линейными функциями аппроксимируется гладкое решение. Таким образом, выбранный метод поз воляет получить решение с требуемой точностью, если только она не превос ходит возможности вычислительной техники.

В верхней части Рис. C.1 представлена работа многосеточного метода при решении уравнения Лапласа. Для задач такого типа обычно итерации сходятся в несколько раз быстрее, чем у более сложных задач, однако в на шем случае имеем замедление сходимости. Такое несоответствие объясняет ся большой вытянутостью ячеек по горизонтали. Так, в показанном решении отношении вертикального размера ячейки сетки к горизонтальному равно 1 : 125. В качестве обоснования этого утверждения на Рис. C.2 представле ны решения этих же задач, но в области со сравнимыми горизонтальным и вертикальным размером 400 400 80 км. В такой области решение урав нения Лапласа находится примерно за столь же малое число итераций, как и решение остальных задач. Экспоненциальный рост проводимости внутри области уменьшает скорость сходимости итераций и уменьшает точность по лученного решения, что видно на Рис. C.1 и Рис. C.2, но на решение таких задач меньше влияет вытянутость ячеек области.

0 log H log V -5 - -10 - -15 - -20 - 0 4 8 12 16 20 0 4 8 12 16 а) уравнение Лапласа (z) = 0 log H log V -5 - -10 - -15 - -20 - 0 4 8 12 16 20 0 4 8 12 16 z/z б) (z) = e, n = m = 0 log H log V -5 - -10 - -15 - -20 - 0 4 8 12 16 20 0 4 8 12 16 z/z в) (z) = e, n = m = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.2. Cходимость многосеточных итераций. Область 400 км 400 км 80 км. В левой колонке показано убывание нормы невязки H. В правой колонке показано убывание разности между решением после данной итерации и точным численным решением Vh (сплошная линия), а также разности между численным и аналитическими решениями V (пунктир). По горизонтали i - номер итерации.

Для такой простой области как параллелепипед, несложно проводить ис следования разных способов увеличения скорости сходимости итераций. Пе ред первой итерацией при решении указанных выше задач в граничных узлах значения потенциала V задавались в соответствии с условиями Дирихле, а значения во всех внутренних узлах задавались нулевыми. Если во внутренних узлах перед первой итерацией задавать некоторое начальное приближение, то при его удачном выборе можно получить результат за меньшее количества итераций.


На Рис. C.3 показана сходимость итераций при использовании в качестве начального приближения решения одномерных в направлении z задач, в ко торых в качестве верхнего и нижнего граничного значений потенциала ис пользованы значения из соответствующих узлов трехмерной сетки, лежащих на верхней и нижней границах:

V (x, y, z) = sin(nx/x0) sin(my/x0 )f (z), (1 exp (z/z0 ))(1 exp (zup/z0)), если (z) ez/z (C.5) f (z) = z/zup, если (z) = const Как видно на Рис. C.3, одномерное начальное приближение позволяет сократить время вычислений за счет уменьшения количества необходимых итераций.

В качестве еще одной тестовой задачи в параллелепипеде рассмотрим воз можности использования неравномерной сетки, когда в одной части области линии сетки расположены более часто, чем в другой. Такую неравномерность можно использовать для более детальной передачи неоднородностей в неко торой подобласти. Например, перед землетрясением на поверхности Земли возможно сильное изменение электрического поля возле разлома, и представ ляет интерес возмущение поля в окружающей части атмосферы. При этом область расчета не ограничена по горизонтали, хотя очевидно, что решение быстро убывает с удалением от источника, и поэтому можно рассматривать 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа (z) = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z б) (z) = e, n = m = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z в) (z) = e, n = m = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.3. Cходимость многосеточных итераций при ненулевом начальном приближении. Область 10 000 км 10 000 км 80 км. В левой колонке показано убывание нормы невязки H. В правой колонке показано убывание разности между решением после данной итерации и точным численным решением Vh (сплошная линия), а также разности между численным и аналитическими решениями V (пунктир). По горизонтали i - номер итерации.

ограниченную область. В ее удаленных частях нет необходимости строить решение точно, и поэтому там может быть использована более грубая сетка.

Рассмотрим задачу такого рода.

Исходную область — параллелепипед |x| b/2, |y| b/2, 0 z zup, разобьем на три части по обоим горизонтальным направлениям. Начало ко ординат находится в центре нижней грани области, оси x и y горизонтальные, z - вертикальная. На нижней и боковой границе области потенциал положим нулевым. В получившемся центральном параллелепипеде на верхней границе |x| a, |y| a, z = zup зададим (C.6) V (x, y, zup) = (1 + cos(2x/a))(1 + cos(2y/a))/4, где a - размер горизонтального ребра центрального параллелепипеда, а вне этого квадрата положим V (x, y, zup) = 0.

Для построения аналитического решения продолжим задачу на весь плос кий слой 0 z zup. Добавим к функции на верхней границе в направлении осей x и y такие же неоднородности, что и в центре (C.6), на расстоянии b, (1 + cos( 2(xb) ))(1 + cos( 2y ), 4 |x b| a, |y| a a a |x| a, |y b| a (C.7) V (x, y, zup) = 1 (1 + cos( 2x ))(1 + cos( 2(yb) ), 4 a a (1 + cos( 2(xb) ))(1 + cos( 2(yb) ), |x b| a, |y b| a 4 a a Продолжим получившуюся функцию вдоль каждой из осей x и y с пери одом 2b:

V (x, y, zup) = V (x + 2b, y, zup) = V (x, y + 2b, zup) = V (x + 2b, y + 2b, zup).

Такой вид функции позволяет найти аналитическое решение исходной за дачи в виде ряда Фурье по x, y с зависящими от z коэффициентами, которые находим аналогично (C.2), решив обыкновенные дифференциальные уравне ния, получающиеся при подстановке этого ряда в уравнение электропровод ности:

(2n 1)x (2m 1)y V (x, y, z) = fnm cos ( ) cos ( ) b b n=1 m= (e1 z e2 z )/(e1 zup e2 zup ) где n и m - целые числа, fnm - коэффициенты Фурье периодической функции (C.6, C.7) (2n 1)x (2m 1)y V (x, y, zup) = fnm cos ( ) cos ( ).

b b n=1 m= Гармоник с четными индексами и членов с синусами в этом разложении нет из-за симметрии функции относительно прямых x = 0, y = 0 и антисим метрии относительно прямых x = b/2,..., y = b/2,... Последняя обеспечивает равенство потенциала нулю на боковых сторонах исходной области.

Рис. C.4. Горизонтальная проекция области с неравномерно заданной сеткой Внутри центрального параллелепипеда построим более мелкую по гори зонтали сетку, как показано на Рис. C.4.

На Рис. C.5 показана сходимость итераций при решении этой задачи.

Из рисунка видно, что решение на неравномерной сетке с меньшей скоро стью, но сходится.

0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа (z) = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z б) (z) = e, n = m = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.5. Cходимость многосеточных итераций на неравномерных сетках.

Область 30 000 км 30 000 км 80 км. В левой колонке показано убывание нормы невязки H. В правой колонке показано убывание разности между решением после данной итерации и точным численным решением Vh (сплошная линия), а также разности между численным и аналитическими решениями V (пунктир). По горизонтали i - номер итерации.

Рис. C.6. Горизонтальная проекция сетки в призме с шестью боковыми гранями C.2. Задача в призме с шестью боковыми гранями Теперь рассмотрим области, получающуюся при соединении трех шестигран ников таким образом, что в горизонтальной проекции получается правиль ный шестиугольник, показанный на Рис. C.6. Зададим те же граничные усло вия, что и в предыдущей задаче. Данная модель была создана для отработки методов расчета в сложных областях.

На Рис. C.6, что над центральной точкой три шестигранника стыкуют ся углами 2/3. Как показывает анализ решения, представленный на Рис.

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

Разность между численным и аналитическим решениями максимальна на центральной линии соприкосновения шестигранников. Для проверки сходи мости решения в области используем две интегральных нормы. Первая норма – энергетическая W (C.1), вторая – норма в пространстве L2():

(C.8) (Vh V )2d, L = где - область, в которой ищется решение.

Рис. C.8 демонстрирует, что решение с уменьшением размера ячейки схо дится по всей области в этих нормах.

Если же интересует решение именно в выделенной подобласти, то ее мож но вырезать, как это показано на Рис. C.9, и дополнительно пересчитать ре шение. Для вырезанной подобласти значения на верхней и нижней границе берем такими же как и в исходной области, а на боковые грани интерполи руется решение, полученное в исходной области.

В центральной части вырезанной подобласти, не включающей несколько слоев сетки вблизи боковой границы, решение сходится и в более сильной норме – максимум модуля разности, как это видно на Рис. C.11. Выделение центральной части в уже вырезанной подобласти необходимо, так как на ее боковых границах были заданы неточные значения, полученные в исходной области, а они влияют на значений в соседних узлах. Влияние неточных бо ковых граничных значений на центральную подобласть мало из-за больших горизонтальных размеров.

C.3. Задача в шаровом слое При использовании комплекса программ для расчетов глобального распре деления электрического потенциала в атмосфере Земли областью явля ется слой между земной поверхностью и ионосферой. Считаем их сферами с радиусами RE = 6400 и Rion = 6480 км соответственно. При симметрии полушарий можно ограничиться северным полушарием с дополнительным граничным условием на экваториальной плоскости.

Всю область разрезаем на пять криволинейных шестигранников: один в центре и четыре вокруг него. На Рис. C.10 показана проекция построенной сетки на экваториальную плоскость. При таком разбиении сетка возле стыков шестигранников имеет четыре особенности, аналогичных описанной в преды дущем пункте, а в центральной подобласти получается более равномерной.

На внутренней полусфере, представляющей собой поверхность Земли, и на боковой границе, лежащей в экваториальной плоскости, полагаем потен циал равным нулю. На верхней полусфере, представляющей собой границу между атмосферой и ионосферой, зададим в виде полинома Лежандра:

V (Rion,, ) = Pn (cos()).

Проводимость внутри области зададим функцией (r) = exp ((r RE )/z0) с характерным для атмосферы Земли значением z0 = 6 км.

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

V (r,, ) = f (r)Pn(), где функция f (r) должна быть решением краевой задачи для обыкновенного дифференциального уравнения d2 df (r) (r (r) ) n(n + 1)(r)f (r) = dr dr (C.9) f (RE ) = 0, f (Rion) = 1.

Эту задачу решаем с помощью численного метода, описанного в разделе 4.4., а при постоянной проводимости, когда уравнение электропроводности становится уравнением Лапласа, двумя решениями уравнения (C.9) являются степенные функции rn и 1/rn+1, следующая линейная комбинация которых удовлетворяет краевым условиям f (r) = (rn rE /rn+1)/(rion rE /rion ).

2n+1 n 2n+1 n+ На Рис. C.12 и C.13 показана сходимость итераций в различных нормах.

Она мало отличается от сходимости для рассмотренной выше задачи в призме с шестью боковыми гранями.

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

0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа (z) = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 б) (z) = ez/z0, n = m = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 в) (z) = ez/z0, n = m = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.7. Cходимость многосеточных итераций на неравномерных сетках в призме с шестью боковыми гранями. Представление аналогично Рис. C. 0 log L log W -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа (z) = 0 log L log W -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z б) (z) = e, n = m = 0 log L log W -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z в) (z) = e, n = m = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.8. Cходимость многосеточных итераций на неравномерных сетках в призме с шестью боковыми гранями. В левой колонке показано убывание разности между численным и точным решениями в норме пространстве L2, в правой колонке – в энергетической норме. По горизонтали i - номер итерации.

Рис. C.9. Горизонтальная проекция сетки в призме с шестью боковыми гранями с выделенной подобластью Рис. C.10. Горизонтальная проекция сетки в шаровом слое 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа (z) = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z б) (z) = e, n = m = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z в) (z) = e, n = m = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.11. Cходимость многосеточных итераций на неравномерных сетках в области, вырезанная в центре призмы с шестью боковыми гранями. В левой колонке показано убывание нормы невязки H. В правой колонке показано убывание разности между решением после данной итерации и точным численным решением Vh (сплошная линия), а также разности между численным и аналитическими решениями V (пунктир). По горизонтали i номер итерации.

0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа (z) = 1, n = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z б) (z) = e, n = 0 log H log V -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z в) (z) = e, n = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.12. Cходимость многосеточных итераций на неравномерных сетках в шаровом слое. В левой колонке показано убывание нормы невязки H. В правой колонке показано убывание разности между решением после данной итерации и точным численным решением Vh (сплошная линия), а также разности между численным и аналитическими решениями V (пунктир).

По горизонтали i - номер итерации.

0 log(L ) log W -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 а) уравнение Лапласа (z) = 1, n = 0 log L log W -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z б) (z) = e, n = 0 log L log W -5 - -10 - -15 - -20 - 0 50 100 150 200 250 0 50 100 150 200 z/z в) (z) = e, n = сетка 16 16 16 сетка 32 32 32 сетка 64 64 Рис. C.13. Cходимость многосеточных итераций на неравномерных сетках в шаровом слое. В левой колонке показано убывание разности между численным и точным решениями в норме пространстве L2, в правой колонке – в энергетической норме. По горизонтали i - номер итерации.



Pages:     | 1 |   ...   | 2 | 3 ||
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.