Создание задачи

Постановка задачи

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

u_{tt} = a^2 \Delta u = a^2 (u_{xx} + u_{yy}).

Для моделирования используем явную трёхслойную схему:

\frac{u^{n + 1}_i - 2 u^{n}_i + u^{n - 1}_i}{\Delta t^2} = \frac{a^2}{V_i} \sum_{j} \frac{u^{n}_j - u^{n}_i}{|r_j - r_i|} S_{ij},

которая записана методом конечных объемов. Здесь $\Delta t$ – временной шаг, $V_i$ – объем ячейки, $S_{ij}$ – площадь грани ячейки, $r_i$ – центр ячейки. Для устойчивости схемы требуется выполнение условия Куранта (приближенная формула):

\Delta t \le \frac{\min(\Delta x, \Delta y)}{\sqrt{2} a}.

Добавление цели

Добавим исходный файл wave.cpp в папку problems. В файле problems/CMakeLists.txt добавим строку

list(APPEND PROBLEMS wave)

Для сборки цели wave необходимо убедиться, что включена опция cmake ZEPHYR_ENABLE_PROBLEMS. В этом случае цель wave появится в выдаче cmake в списке компилируемых задач. Перейдем к написанию файла wave.cpp.

Полный листинг задачи

Подключение основных модулей

Подключим основные заголовочные файлы:

#include <zephyr/geom/generator/rectangle.h>
#include <zephyr/mesh/euler/eu_mesh.h>
#include <zephyr/io/pvd_file.h>

Файл generator/rectangle.h включает сеточный генератор, euler/eu_mesh.h содержит класс для работы с эйлеровыми сетками, io/pvd_file.h содержит функции для сохранения расчетных данных.

Разрешим использовать все пространства имен:

using namespace zephyr::geom;
using namespace zephyr::mesh;
using namespace zephyr::io;
using generator::Rectangle;

Создание сетки

Cоздадим сеточный генератор Rectangle:

int main() {
    Rectangle gen(-2.0, 1.0, -2.0, 1.0);  // [x_min, x_max, y_min, y_max]
    gen.set_nx(400);
    gen.set_boundaries({.left=Boundary::WALL, .right=Boundary::WALL,
                        .bottom=Boundary::WALL, .top=Boundary::WALL});

Генератор Rectangle в качестве аргументов принимает границы области по осям, в нашем случае это границы квадрата $[-2, 1] \times [-2, 1]$ . Далее выставляется число ячеек по оси $x$ , разбиение по оси $y$ выбирается автоматически. В последней строке задаются граничные условия. Флаги граничных условий:

Boundary::UNDEFINED – пропускать грань (по умолчанию);
Boundary::ORDINARY – внутренняя грань между двумя ячейками;`
Boundary::WALL – отражающие граничные условия (стенка);
Boundary::ZOE – повторяет внутреннюю ячейку (zero order extrapolation);
Boundary::PERIODIC – периодические граничные условия.

В нашем случае (для закрепленной мембраны) целесообразно указать Boundary::WALL. Флаги граничных условий потребуется обрабатывать в решателе.

Создаем расчетную сетку EuMesh и добавляем расчетные поля:

    EuMesh mesh(gen);

    Storable<double> u_prev = mesh.add<double>("u_prev");
    Storable<double> u_curr = mesh.add<double>("u_curr");
    Storable<double> u_next = mesh.add<double>("u_next");

Схема является трёхслойной, поэтому необходимо хранить на сетке три параметра u_prev – предыдущий временной слой, u_curr – текущий, u_next – следующий. На сетку можно добавлять не только тип double, но и пользовательские типы.

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

Storable<T> key = mesh.add<T>(name);

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

Запись в файл

Создадим экземпляр PvdFile для сохранения данных расчета:

    PvdFile pvd("wave", "output");
    pvd.variables.append("u", u_curr);

В данном случае "output" – название директории, "wave" – название файла без расширения. Итоговый файл будет иметь название "wave.pvd". В переменные для сохраненияи добавлен ключ u_curr под именем "u". Результаты выполнения программы можно визуализировать в ParaView.

Задание начальных условий

Выставим начальные условия:

    for (auto cell: mesh) {
        double r = cell.center().norm();
        cell[u_prev] = 0.5 + 0.5 * std::tanh(50 * (0.5 - r));
        cell[u_curr] = cell[u_prev];
    }

Для обхода ячеек сетки доступен range based цикл. Чтобы получить центр ячейки используется функция cell.center(), последующий вызов .norm() возвращает радиус центра ячейки. Доступ к полям ячейки выполняется через ключи u_prev, u_curr и оператор []. Другие функции ячеек.

Определим шаг интегрирования по времени с учетом условия Куранта:

    double speed = 1.0;
    double h = std::min(mesh[0].hx(), mesh[0].hy());
    double dt = 0.95 * h / (std::sqrt(2.0) * speed);

Вызов mesh[0] позволяет получить первую ячейку сетки, функции hx() и hy() – размеры ячейки.

Основной цикл программы

Каждые 10 шагов будем делать вывод в консоль и сохранение в файл:

    size_t n_step = 0;
    double curr_time = dt;
    while (curr_time < 5.0 && n_step < 5000) {
        if (n_step % 10 == 0) {
            std::cout << "Step: " << n_step << ";\tTime: " << curr_time << "\n";
            pvd.save(mesh, curr_time);
        }

В основном цикле по ячейкам необходимо посчитать значения u_next в соответствии с численной схемой:

        for (auto cell: mesh) {
            double uc = cell[u_curr];

            double sum = 0.0;
            for (auto face: cell.faces()) {
                double un = face.neib(u_curr);
                if (face.is_boundary()) {
                    un = -uc;
                }

                double dr = (face.center() - cell.center()).norm();
                sum += (un - uc) * face.area() / (2.0 * dr);
            }

            cell[u_next] = 2.0 * cell[u_curr] - cell[u_prev] +
                std::pow(speed * dt, 2) * sum / cell.volume();
        }

Для обхода по граням ячейки также используется range based цикл. Другие функции граней:

face.normal() – внешняя нормаль грани;
face.area() – площадь (длина) грани;
face.is_boundary() – граничное условие?
face.flag() – флаг граничных условий;
face.neib() – возвращает соседнюю ячейку через грань;
face.neib(Storable<T> key) – данные соседней ячейки, эквивалентно вызову face.neib()[key].

Получение соседней ячейки face.neib() через грань всегда корректно. Если соседней ячейки не существует (для ячеек на границе), тогда функция вернет исходную ячейку. Аналогично для функции face.neib(Storable<T>): если соседа не существует, то мы получим данные самой ячейки.

Записываем значение текущей ячейки в переменную uc, значение соседней ячейки в un. Если грань расположена на границе (проверка face.is_boundary()), то в качестве un используем отраженное значение uc.

После вычисления новых значений u_next для всех ячеек, смещаем временные слои: u_curr $\to$ u_prev, u_next $\to$ u_curr.

        for (auto cell: mesh) {
            cell[u_prev] = cell[u_curr];
            cell[u_curr] = cell[u_next];
        }

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

        curr_time += dt;
        ++n_step;
    } // конец основного цикла
    pvd.save(mesh, curr_time);
    return 0;
} // конец main

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

Постановка задачи​

Добавление цели​

Полный листинг задачи​

Подключение основных модулей​

Создание сетки​

Запись в файл​

Задание начальных условий​

Основной цикл программы​