C ++ | Моделирование частиц в ячейке для физики плазмы — решатель уравнения Пуассона над стандартной сеткой
Я создаю свой собственный набор процессов для задачи, которая стоит передо мной, однако в настоящее время я застрял в трехмерном уравнении Пуассона для решения потенциалов и магнитных полей, создаваемых движущимися заряженными частицами.
Короче говоря, Particle In Cell состоит из 4 основных процессов:
- Интегрирование движущихся частиц (из уравнения Лоренца),
- распределение плотностей заряда в узлах сетки (и плотности тока для магнитных полей),
- решение уравнения Пуассона для вычисления потенциалов и электрических / магнитных полей в узлах сетки,
- распределение электрических / магнитных полей на частицы.
И цикл повторяется. Я могу выполнить шаги (1), (2) и (4), но не часть Пуассона. Я интегрирую частицы по схеме Бориса, сетка просто построена как cell[i][j][k]с размерами K x L x M, и процесс распределения немного сложен, но в конце концов не очень сложен.
Однако, к моему большому удивлению, я не могу найти подходящую библиотеку, которая бы решала трехмерное уравнение Пуассона с помощью метода конечных разностей (FDM). Я начинаю думать, что мне придется написать свой собственный решатель, но решение дифференциальных уравнений в частных производных само по себе является наукой, и мне сказали не пробовать это решение, поскольку результаты могут быть нестабильными, если я не сделаю это должным образом. БПФ также может быть решением, но мне нужно что-то, что будет работать, даже если граничные условия не являются периодическими.
Мой вопрос:
Знаете ли вы какую-либо подходящую библиотеку для C ++, которая возьмет правую часть уравнения Пуассона и вычислит результирующие потенциалы? Вычисление должно быть сделано через FDM. Кроме того, правая часть уравнения Пуассона фактически означает всю сетку из точек K x L x M, которая интегрируется по конечным разностям — ребрам сетки.
Сложно ли написать свой собственный решатель для уравнения Пуассона? Какой алгоритм вы рекомендуете?
Знаете ли вы какие-либо другие решения? Свободное программное обеспечение, хорошо написанная статья или документация, дипломная работа, что угодно.
Заранее спасибо, в настоящее время это очень серьезная проблема для меня, так как я не могу продолжать свою работу.
Решение
Вы можете попробовать некоторые библиотеки, такие как libmesh или же сделка II. Они сделаны из такого рода проблем …
Другие решения
Вы могли бы пойти на самодельный подход.
Сначала напишите систему уравнений, которую вам придется решить.
Определите, сколько у вас будет переменных, например, 3 переменные на ячейку x N ячеек.
Определите, сколько уравнений свяжет эти переменные.
Вы должны иметь одинаковое количество переменных и количество уравнений, скажем, N.
Затем посмотрите, являются ли ваши уравнения линейными. Если это так, вы сможете написать свою систему как AX = B, где A — матрица NxN, B — вектор, а X содержит все переменные (неизвестные), которые вы определили ранее.
Тогда решить линейную систему A * X = B довольно просто, например, с помощью MATLAB или с помощью LU-решателя (код C / C ++ для этого легко найти и занимает, может быть, 200 строк).
Если ваша система нелинейна, вы можете использовать метод Ньютона-Рафсона.
Вы должны написать все свои уравнения в виде … = 0.
Затем вы должны вычислить все производные всех ваших уравнений и поместить их в так называемую матрицу Якоби.
Производная уравнения n ° i по переменной n ° j перейдет в строку n ° i и столбец n ° j якобиана. То есть J [i] [j] = df [i] _dX [j].
Вы должны угадать некоторые начальные значения для всех ваших переменных в X, которые должны быть не слишком далеко от решения.
Это предположение не является решением, поэтому оно не приведет к тому, что все ваши уравнения будут = 0.
Вектор, который содержит все эти значения «должен быть 0», представляет собой F. F [i] — это «остаток» уравнения n ° i;
Когда вы вычислили J и F_X, вы можете решить линейную систему J * D = F_X. Тогда D содержит разницу между вашими текущими переменными в векторе X и решением, которое вы ищете.
Затем вы вычисляете новый X, равный X — D.
Если система не линейна, новый X будет ближе к решению, но вам придется повторить процесс несколько раз (вычислить J, F_X, решить J * D = F_X и установить X = X — D), чтобы приблизиться достаточно для решения.
Это может показаться слишком сложным для установки, но на самом деле это выполнимо.
Что ж, в 1D это почти легко, но в 3D граничные условия могут быть утомительными.
редактировать
Кроме того, вам придется выбирать, как выражать пространственные производные ваших переменных, чтобы вы могли преобразовать свои дифференциальные уравнения в дискретные уравнения, дискретизацию АКА. Вы можете начать со схем первого порядка:
- df / dx = (f [i + 1] — f [i]) / dx
- d2f / dx2 = (2 * f [i] — f [i + 1] — f [i-1]) / dx²