Диссертация (1145383), страница 12
Текст из файла (страница 12)
Как легко видеть из оценки (2.32), успешность -уточнениязависит от гладкости решения. Если решение обладает производными высокихпорядков, то скорость сходимости может быть сверхстепенной, и даже экспоненциальной в случае аналитичности решения. Если же гладкость решения ограничена, то увеличение степени не приводит к увеличению точности дискретногорешения. Именно этим соображением объясняется использование полиномовневысоких степеней в технических приложениях.В любом из этих итеративных подходов, важную роль играют методы оценки погрешности построенного решения на каждом КЭ. Основные методы могутбыть разделены на несколько классов [127]:∙ Интерполяционные методы.
После построения численного решения, оцениваются неизвестные константы в априорных оценках погрешностей.∙ Методы оценки невязки. Построение локального, на некоторой подобласти или КЭ, решения с граничными условиями, полученными из предварительно найденного глобального решения. Сравнение этих двух решенийпозволяет оценить погрешность на подобласти.∙ Методы экстраполяции. На основе нескольких полученных решений строится экстраполяция, которая затем сравнивается с самими решениями.∙ Методы восстановления решений. После вычисления решения, можно найти его значения или значения его градиента с более высокой точностью на73некотором подмножестве. Сравнивая исходное решение и решение, восстановленное с этого подмножества, можно получить оценку погрешности.Все перечисленные методы активно и успешно используются в приложениях [129].В данной работе была исследована эффективность использования для оценки локальной погрешности одного из методов четвёртого типа – метода сверхсходимости.
Теоретические и практические аспекты этого метода развиваютсясо середины 80-х годов прошлого века [130], имеются результаты для эллиптических уравнений второго и четвёртого порядков, а также гиперболическихи параболических уравнений. Имеются также работы, посвящённые явлениюсверхсходимости в спектральных задачах [131, 132].Эффект сверхсходимости заключается в том, что для некоторого множества точек скорость сходимости приближённого решения к точному на одинпорядок по ℎ больше, чем это следует из оценки (2.32).
Эти множества известны явно, причём для самого решения и его градиента они различны [133].Важно отметить, что наличие такого эффекта не противоречит неулучшаемости оценки (2.32): она верна для всех точек области, тогда как сверхсходимостьприсутствует только для некоторых точек.Собственные значения исследуемого уравнения определяются матричными элементами, включающими первые производные. Поскольку скорость сходимости первых производных на один порядок по ℎ меньше, именно они даютосновной вклад в погрешность собственных значений, и именно этот вклад необходимо минимизировать.
C этой целью вычисляется градиент приближенногорешения в точках сверхсходимости на рассматриваемом КЭ и на соседних с ним.По этому избыточному множеству точек с помощью метода наименьших квадратов строится новое приближение к градиенту, имеющее на порядок большуюскорость сходимости. Это приближение используется для оценки локальной погрешности в целях ℎ-уточнения.Опишем используемый алгоритм боле подробно.
Каждый шаг адаптивной74процедуры заключается в обработке приближённого градиента ∇ΨFEM полученного решения ΨFEM . В плоскости = 90 на каждом КЭ с помощью метода наименьших квадратов строится уточнённый градиент ∇ΨFEM−LS . Длявычисления уточнённого градиента для каждого КЭ берутся все точки сверхсходимости внутри элемента и ближайшие к границам точки сверхсходимостис соседних элементов. Далее определяется норма ошибки – разности двух нормградиентов на каждом КЭ :Z =⃒ ⃒⃒)︀(︀⃒⃒∇ΨFEM−LS ⃒ − ⃒∇ΨFEM ⃒ 2 (2.37)по которой и происходит оценка КЭ с наибольшими погрешностями.
Пусть –некоторая заранее заданная абсолютная погрешность скалярного квадрата градиента производной. Алгоритм уточнения начинает работу, если выполняетсяусловие < max .Вычисляется средняя по всем КЭ погрешность^ =1 ∑︁ , где – количество КЭ. Разбиение элемента производится, когда величина соответствующей погрешности на этом элементе больше средней:> 1.^Количество новых элементов вычисляется согласно асимптотической формуле оценки погрешности на конечном элементе (2.33) так, чтобы ошибка накаждом из вновь полученных элементов не превышала :[︂(︁ )︁ ]︂ 1/ =,где [] – целая часть числа . Для того, чтобы предотвратить слишком мелкоеделение элементов, дополнительно вводится параметр , определяющий максимально возможное количество подэлементов. Окончательно, количество новых75элементов вычисляется как{︂[︂(︁ )︁ ]︂ }︂ 1/ = min, .(2.38)На рисунке 2.1 на примере вычисления энергии основного состояния тримера неона (см.
раздел 3.3.2) показано функционирование алгоритма. Точкислева направо описывают параметры численной аппроксимации на последовательных шагах итеративного уточнения: размер МКЭ матрицы и погрешностьэнергии, полученной с этой матрицей. График может быть неплохо аппроксимирован прямой линией, что, с учётом его масштаба, означает обратный степенной закон для зависимости погрешности от размера матрицы. Это наблюдениетакже поддерживает предположение (2.35), использованное для вывода экстраполяционной формулы (2.36). Более подробное описание алгоритма и анализего работы представлены в работе [125].Предложенный выше алгоритм не лишён некоторых недостатков.
Основной из них – то, что МКЭ-разбиение оптимизируется для одной волновой функции, в то время как большинство представленных ниже приложений требуетрасчёта нескольких (а иногда – очень многих) волновых функций на одной итой же сетке. В связи с этим, в основном использовалась часть алгоритма, позволяющая оценить локальную погрешность с помощью сверхсходимости, безпоследующего адаптивного ℎ-уточнения.2.5.
Особенности программной реализацииОписанный в данном разделе численный подход был реализован в программе ACESPA, написанной на языке fortran 90. Основная последовательностьработы программы выглядит следующим образом:∙ Ввод начальных данных (хранящихся в отдельном файле), описывающихпараметры исследуемой системы и численной реализации.
Проверка совместимости параметров, генерация МКЭ разбиения и наборов функций.76Рис. 2.1. Погрешность вычисления энергии основного состояния тримера неона как функция размера глобальной МКЭ матрицы. Треугольниками отмечены последовательные шагиадаптивного алгоритма.∙ Вычисление элементов разреженной матрицы T (2.20), отвечающих операторам (2.2,2.4) для каждой компоненты, запись этих элементов на диск.∙ Решение алгебраической задачи для построенной матрицы: обобщённойзадачи на собственные значения для поиска собственных значений илирезонансов, системы линейных алгебраических уравнений для задачи рассеяния.∙ Постобработка полученных решений (если необходимо): восстановлениеволновых функций, вычисление матричных элементов операторов.Программа и файл данных поддерживают простейшие команды, позволяющиеорганизовать многократный расчёт с разными параметрами системы и численной схемы.
Разработанная программа может работать под управлением операционных систем семейств Windows и Linux. Простая переносимость обеспечива77ет как удобную среду разработки и отладки программы, так и возможность еёработы на высокопроизводительных системах.В связи с особенностями современных вычислительных систем, обеспечивающих высокую производительность за счёт параллельного выполнения, были предприняты специальные усилия для выбора наиболее эффективных программных средств для реализации параллельного выполнения программы. Наибольшие вычислительные ресурсы требуются при вычислении матричных элементов операторов и вектора правой части, а также при восстановлении волновой функции. Именно эти части программы были распараллелены. Были исследованы реализации с использованием OpenMP, MPI, а также GPGPU [52, 57, 62,124].
Было показано [57], что наиболее эффективным и масштабируемым средством для выбранного численного метода является MPI. В рамках МКЭ, естественное распараллеливание обеспечивается вычислением всех матричных элементов операторов для одного КЭ на отдельном процессорном элементе (ПЭ).Такой подход позволяет достаточно эффективно использовать любое имеющееся в вычислительной системе количество ПЭ, по крайней мере, если количествоПЭ меньше, чем количество КЭ в МКЭ разбиении.
Если это не так, схема разбиения вычислительной нагрузки на ПЭ должна быть усложнена, однако покарасчёты на таких вычислительных системах не производились.Расчёты рассеяния электрона на атоме водорода и ионе гелия показывают [62], что ускорение расчётов оказывается не менее, чем 0.8 * количество(ПЭ)для вычисления матричных элементов операторов (около 90% общего временивычислений) и не менее, чем 0.5 * количество(ПЭ) для вычисления вектораправой части, при использовании до 24 ПЭ.Для решения разреженных систем линейных алгебраических уравнений использовались два свободно распространяемых программных пакета: PARDISO [134,135] и MUMPS [136, 137]. Оба этих пакета также могут использовать доступные параллельные вычислительные ресурсы, хотя их параллельная эффективность оказывается ниже, чем при вычислении матричных элементов. Для ре78шения обобщённой задачи на собственные значения использовался метод Арнольди [138], а именно его реализация в свободно распространяемом пакетеARPACK [139].
Для поиска собственных значений, ближайших к некоторомузначению 0 , используется спектральное преобразование уравнения (2.1) к виду(︁^ − 0 ^)︁−1^ = Ψ,Ψгде =1. − 0(2.39)^ симметричен, то оператор в левой части этого уравненияЕсли оператор также будет симметричен в скалярном произведении⟨, ⟩ = ⊤ ^ .Такое преобразование чрезвычайно сильно ускоряет сходимость собственныхзначений в методе Арнольди, позволяя быстро вычислять до нескольких десятков собственных значений для выбранного 0 . Недостатком такого подхода явля)︁(︁^^ется необходимость факторизации (впрочем, однократной) матрицы − 0 .2.6. Выводы ко второй главеВ данной главе описан вычислительный подход, разработанный и реализованный для решения систем трёхмерных дифференциальных уравнений вчастных производных, сформулированных в первой главе.
Особенность этихуравнений состоит в том, что входящие в них операторы не являются, вообщеговоря, симметричными, что сужает круг доступных вычислительных методов.Использованный подход основан на записи уравнений в виде вариационногоуравнения (не задачи минимизации!) и применении к этому уравнению подхода Галёркина. В качестве базисных функций для галёркинского разложенияиспользуются функции, построенные с помощью метода конечных элементов.Их основное достоинство – хорошая локальность: они отличны от нуля лишьна малом подмножестве полного конфигурационного пространства. Такое свойство приводит к большой разреженности итоговых матриц систем линейных79уравнений, что существенно понижает вычислительную стоимость их численного решения.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
