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

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

Первые шаги с вариационными формами

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

Код начинается с создания сетки и функционального пространства:

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

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

Далее следует определение граничных условий. Пусть на границах x = 0 и x = 1 решение должно обращаться в ноль:

def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

Функция boundary возвращает булево значение: она истинна на левой и правой границах квадрата. DOLFIN_EPS обеспечивает корректное сравнение вещественных чисел с учетом машинной точности.

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

u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)

a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

Билинейная форма a записывается как скалярное произведение градиентов, проинтегрированное по области dx. Линейная форма L включает источниковый член f и граничный поток g. Обратите внимание: grad(), inner(), dx, ds - все это операторы UFL, языка вариационных форм. Математическая нотация напрямую транслируется в код.

Решение задачи теперь тривиально:

u = Function(V)
solve(a == L, u, bc)

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

Погружение в нелинейность

Линейные задачи хороши для обучения, но реальный мир нелинеен. Рассмотрим модифицированное уравнение Пуассона, где коэффициент диффузии зависит от самого решения. Математически это записывается как дивергенция произведения (1 + u²) на градиент u равна источнику.

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

u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])", degree=2)

F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx

Здесь u уже не TrialFunction, а Function. Это сигнализирует FEniCS о нелинейности задачи. Форма F должна обращаться в ноль при истинном решении. Метод Ньютона требует якобиан - производную F по u. FEniCS вычисляет его автоматически через символьное дифференцирование:

solve(F == 0, u, bc)

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

solver_parameters = {"newton_solver": {"relative_tolerance": 1e-6}}
solve(F == 0, u, bc, solver_parameters=solver_parameters)

Адаптация сетки через индикаторы ошибки

Равномерная сетка расточительна. Зачем детально разрешать области, где решение почти постоянно? Адаптивное разбиение концентрирует элементы там, где происходит действие. Классический подход в FEniCS использовал встроенные адаптивные решатели, но современная практика опирается на явное управление циклом адаптации.

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

mesh = UnitSquareMesh(8, 8)

for iteration in range(5):
    V = FunctionSpace(mesh, "Lagrange", 1)
    # Решение задачи на текущей сетке
    u = solve_problem(V)
    
    # Вычисление индикатора ошибки
    DG = FunctionSpace(mesh, "DG", 0)
    error = project(sqrt(inner(grad(u), grad(u))), DG)
    
    # Маркировка элементов
    cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
    threshold = 0.5 * error.vector().max()
    
    for cell in cells(mesh):
        cell_markers[cell] = error.vector()[cell.index()] > threshold
    
    # Измельчение сетки
    mesh = refine(mesh, cell_markers)

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

Энергетический подход использует плотность энергии деформации. Для задач механики это естественная мера локальной точности:

energy_density = 0.5 * inner(grad(u), grad(u))
error = project(energy_density, DG)

Критерий Дёрфлера обеспечивает, что разбиваются элементы, вносящие заданную долю суммарной ошибки. Сортируем элементы по убыванию индикатора, накапливаем сумму, останавливаемся при достижении порога:

error_values = [(cell.index(), error.vector()[cell.index()]) for cell in cells(mesh)]
error_values.sort(key=lambda x: x[1], reverse=True)

total_error = sum(e[1] for e in error_values)
threshold_error = 0.5 * total_error
accumulated = 0.0

for idx, err in error_values:
    cell_markers[idx] = True
    accumulated += err
    if accumulated >= threshold_error:
        break

Этот алгоритм гарантирует контролируемый рост числа элементов при максимальном снижении ошибки.

Работа со смешанными формулировками

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

Элементы Тейлора-Худа решают проблему: квадратичные функции для скорости, линейные для давления. В FEniCS смешанное пространство создается через умножение:

P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

Пробные и тестовые функции извлекаются через split:

(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

Вариационная форма уравнений Стокса записывается естественно:

a = (inner(grad(u), grad(v)) - div(v)*p - div(u)*q)*dx
L = inner(f, v)*dx

Граничные условия применяются к компонентам смешанного пространства:

bc_velocity = DirichletBC(W.sub(0), Constant((0, 0)), boundary)

W.sub(0) выделяет подпространство скорости. Решение возвращает объединенную функцию, которую можно разделить:

w = Function(W)
solve(a == L, w, bc_velocity)
(u, p) = w.split()

Временная зависимость и адаптация

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

dt = 0.01
u_n = interpolate(u_initial, V)
u = TrialFunction(V)
v = TestFunction(V)

a = u*v*dx + dt*inner(grad(u), grad(v))*dx
L = u_n*v*dx

u = Function(V)
t = 0
while t < T:
    solve(a == L, u, bc)
    u_n.assign(u)
    t += dt

На каждом шаге решается линейная система. Билинейная форма включает массовую матрицу (произведение функций) и матрицу жесткости (произведение градиентов). Коэффициент dt балансирует их вклад.

Адаптация во времени требует перестроения пространства на каждом шаге. После измельчения сетки старое решение проецируется на новое пространство:

if should_adapt:
    mesh = refine(mesh, cell_markers)
    V_new = FunctionSpace(mesh, "Lagrange", 1)
    u_n_new = Function(V_new)
    u_n_new.assign(interpolate(u_n, V_new))
    u_n = u_n_new
    V = V_new

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

Современный FEniCSx и его возможности

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

import dolfinx
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

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

from dolfinx.fem.petsc import LinearProblem

problem = LinearProblem(a, L, bcs=[bc])
uh = problem.solve()

Интеграция с PETSc обеспечивает доступ к продвинутым решателям и предобуславливателям. Параметры передаются через словарь:

petsc_options = {"ksp_type": "cg", "pc_type": "hypre"}
problem = LinearProblem(a, L, bcs=[bc], petsc_options=petsc_options)

Адаптивное разбиение в FEniCSx поддерживается через внешние инструменты. Генерация сеток с помощью Gmsh, адаптация через библиотеку p4est. Интерфейс стал более модульным, позволяя комбинировать компоненты по необходимости.

Практические наблюдения

Квадратурные правила влияют на точность и скорость. FEniCS автоматически выбирает степень квадратуры, достаточную для точного интегрирования полиномов заданной степени. Для элементов первого порядка используются точки Гаусса второго порядка. Явное управление квадратурой возможно через метаданные:

dx_custom = dx(metadata={"quadrature_degree": 4})
a = inner(grad(u), grad(v))*dx_custom

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

Сохранение результатов для визуализации использует формат XDMF, понятный ParaView и VisIt:

from dolfinx import io

with io.XDMFFile(mesh.comm, "solution.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(uh)

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

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