CC BY
52
152
УДК 519.63
В. А. Седайкина
ОПТИМИЗАЦИЯ МЕТОДА FIM ДЛЯ РЕШЕНИЯ ТРЕХМЕРНОГО УРАВНЕНИЯ ЭЙКОНАЛА
Быстрый итерационный метод (FIM), описанный в статье [1 ], является удобным и действенным для решения уравнения эйконала для двумерных и небольших трехмерных задач. Однако для решения больших трехмерных задач требуется довольно продолжительное время и большое количество оперативной памяти. Приведена оптимизация этого метода, сокращающая объем хранимой в памяти информации и время, необходимое для решения задачи. Приведены численные эксперименты.
A Fast Iterative Method (FIM) for eikonal equations [1] are convenient and optimal in two-dimensional and small three-dimensional models. But it takes a lot of time and random access memory for big three-dimensional models. This article describes the time and memory optimization of this method for big 3d models. The results of simulations are presented.
Ключевые слова: уравнение эйконала, численное моделирование, оптимизация, быстрый итерационный метод.
Key words: eikonal equation, numerical modeling, optimization, fast iterative method.
© В. А. Седайкина, 2015
Вестник Балтийского федерального университета им. И. Канта. 2015. Вып. 10. С. 152-159.
Быстрый итерационный метод решения уравнения эйконала
Уравнение эйконала — нелинейное дифференциальное уравнение первого порядка
Кх)|2 --Г^ = 0, Vx eQ (1)
v (x)
где Q — область в R3; т( x) — эйконал; v( x) — функция скорости.
В [1] приведен быстрый итерационный метод (FIM) решения уравнения эйконала.
Краткое описание метода FIM
Дано: скоростная модель и координаты источника (ов).
Найти: решение уравнения (1) во всех точках сетки.
Метод можно разделить на следующие основные пункты.
1. Задается некоторое начальное приближение: точка источника и несколько точек, ближайших к ней, добавляются в «действующий список», представляющий собой массив, элементы которого будут пересчитаны на каждом шаге итерации;
2. Каждый элемент «действующего списка» пересчитывается способом, описанным в [1] и сравнивается со своим прошлым значением (полученным на предыдущей итерации). Если разница между этими значениями по модулю:
а) больше некоторой погрешности (погрешность выбирается пользователем), элемент остается в «действующем списке»;
б) меньше погрешности, то элемент покидает «действующий список». Затем пересчитываются не находящиеся в «действующем списке» соседи этого элемента в основных направлениях (4 соседа в 2D и 6 в 3D случаях). Соседи, чьи значения при пересчете оказались меньше, попадают в «действующий список».
3. Пункт 2 повторяется до тех пор, пока «действующий список» не опустеет. Тогда решение во всей области считается найденным с заданной точностью.
Эффективность данного алгоритма обоснована в статье [1].
Заметим, что узкое место алгоритма в расширяющемся «действующем списке», который постоянно пересчитывается; чем больше становится список, тем больше времени занимает одна итерация. Для двумерной задачи это не является критичным, однако при переходе к трехмерной задаче это становится сильно заметным. Также при решении больших трехмерных задач возникают проблемы, связанные с необходимостью хранения в оперативной памяти больших объемов информации. Рассмотрим этот вопрос подробнее.
Для решения уравнения (1) в области размером NxNyNz (где
Nx, Ny, Nz — количество отсчетов по осям Ox, Oy, Oz) требуются:
1) массив С, содержащий скорость в среде, и массив T, содержащий решения уравнения (1). Эти массивы имеют одинаковый размер 4 NxNyNz (тип данных float) или 8 NxNyNz (тип данных double);
154
2) логический массив x_log (состоящий из 0 и 1), показывающий, находится ли точка области в «действующем списке» или нет. Размер массива NxNyNz (тип данных boolean);
3) массив L, содержащий «действующий список», либо динамический (его размер будет постоянно изменяться при решении задачи), либо статический (с предопределенным размером). В случае изменяющегося размера есть существенная потеря во времени при перезаписи массива, а в случае предопределенного размера — потеря в памяти, так как сразу используется максимально возможный размер, равный 4 NxNyNz (тип данных integer). Дополнительно это значение должно
быть умножено на 3, поскольку хранятся координаты (x, y, z) для каждой точки.
То есть, например, для решения уравнения (1) в области 801 • 801 • 601 в оперативной памяти нужно хранить: 801 • 801 • 601 • 4 • 2 байт (массивы C и Т) + 801 • 801 • 601 • 1 байт (массив x_log) + 801 • 801 • 601 • 3 • 4 байт (массив L максимально возможного размера) 8097646221 байт « 7,5 Гигабайт.
Это затраты только на хранение основных данных, не учитывая самих расчетов и вспомогательных данных. При таких затратах памяти программа очень сильно теряет в скорости или даже может не работать на некоторых компьютерах. Таким образом, требуется оптимизация, позволяющая уменьшить затраты памяти и увеличить скорость выполнения программы.
Оптимизация
Пусть расчетная область имеет форму прямоугольного параллелепипеда и находится в пределах x е [0, Nx ], y е [0, Ny ], z е [0, Nz ].
Первый шаг оптимизации — это ограничение размера «действующего списка»
Сначала рассчитывается начальное приближение — это куб, в центре которого находится источник. Если источников несколько, рекомендуется выбрать «опорный источник», находящийся примерно посередине области. Расстояние от источника до любой из плоскостей куба равно 9N (N — выбранное пользователем число, величина которого определяется опытным путем из следующих соображений: при слишком большом N решение задачи даже в кубе с начальными значениям становится неоправданно долгим, а при слишком маленьком теряется точность. Оптимальным можно считать N е [3, 7]). В этом кубе вычисляются значения уравнения (1) с помощью метода FIM без изменений.
Решение в этой области будет начальным приближением для запуска следующей итерационной схемы:
1. Пусть источник находится в координатах (x0, y0, z0). Расчетная область увеличивается пропорционально по каждой координате в обе стороны на значение N. То есть если
x1 е [x0 - 9N; x0 + 9N], y1 е[y0 - 9N; y0 + 9N], z1 е[z0 - 9N; z0 + 9N]
суть границы области начального приближения, то границы новой увеличенной области будут следующими:
xt _ new e [min(xt) - N; max(xt) + N], y1_ new e [min() - N; max() + N],
z1 _ new e [min( z1) - N; max(z:) + N].
В центре этой новой области находится куб с уже посчитанными значениями, а по бокам — неизвестные значения. 2. Далее расчет происходит по следующей схеме. 2.1. Рассчитываются 8 кубов в углах новой области (рис. 1, а).
Т
s
155
Рис. 1. Расчет начальных кубов
Размер одной грани каждого такого куба БЫ. На рисунке 1, б изображен увеличенный куб, на который указывает стрелка на рисунке 1а. Серым цветом схематично изображено область, данные в которой посчитаны на предыдущем шаге, а белым — данные, которые должны быть рассчитаны на этом. Размер грани серого куба 4Ы.
Задача решается в каждом из этих восьми кубов. Так как они не пересекаются, эти вычисления можно выполнить параллельно.
2.2. Расчет прямоугольных параллелепипедов, расположенных на ребрах области.
Чтобы иметь возможность распараллеливания, предлагается сначала найти решение задачи в областях, указанных на рисунке 2, а (так как они не имеют пересечений), затем в областях, указанных на рисунках 2, б и 2, в
Рис. 2. Расчет прямоугольных параллелепипедов
На рисунке 2, г приведено схематическое изображение одной из этих областей (область со стрелкой на рис. 2, б). Здесь темно-серым цве-
о
а
156
том изображены значения, полученные на предыдущем шаге итерации, а светло-серым — значения, рассчитанные на текущем шаге итерации, в пункте 2.1, белым цветом — неизвестные значения.
Так как с каждым шагом итерации длина прямоугольных параллелепипедов будет увеличиваться вместе с областью, вычисление значений в них тоже должно быть разбито на части.
Область прямоугольного параллелепипеда, начиная с угла (значения в котором посчитаны в пункте 2.1) разбивается на небольшие области по следующему принципу: новая область имеет длину ребра БЫ, причем из них 3Ы — известные значения, а 2Ы — неизвестные. На рисунке 3а светло-серым цветом изображены значения, которые должны быть посчитаны при решении задачи в области, отмеченной на рисунке 3а цифрой 1, а белым — в области 2. То есть для расчета новых значений область 2 использует известные данные, рассчитанные ранее в области 1 и в предыдущей области.
( [
/ /
) /
/ У
// V / V)}
\
/ / ' / /
/ / /
Рис. 3. Расчет отдельного параллелепипеда
Значения начинают рассчитываться от угловых кубов в сторону центра (схема на рис. 3, б). Когда значения во всем параллелепипеде рассчитаны, возле центра (с отступом на БЫ в одну и в другую сторону) формируется общая область и значения в ней пересчитываются для выравнивания (рис. 3, в, область 1).
2.3. Расчет плоскостей.
Чтобы иметь возможность распараллелить вычисление, предлагается рассчитывать попарно независимые участки — две боковых, как на рисунке 4, а. Аналогично верхнюю и нижнюю, и переднюю и заднюю. По контуру этих областей известные значения, посчитанные на шаге 2.2, а с внутренней (относительно всей области) стороны — значения, полученные на предыдущем шаге итерации. Для каждой такой области повторяются шаги 2.1 — 2.2.
// //
// , >У
б
а
в
б
а
в
Рис. 4. Расчет плоскостей
Сначала снова считаются значения в угловых кубах, на рисунке 4б — темно-серые области (на рис. 4б и 4в изображены сечения в плоскости YZ при x = min(x:)), а затем — в прямоугольных параллелепипедах, светло-серые области на рисунке 4а. Далее область уменьшается (рис. 4в). Самый темный цвет — значения, рассчитанные на предыдущем шаге (рис. 4б). И шаги 2.1 — 2.2 повторяются, пока значения во всей области не будут рассчитаны.
3. Решение задачи в увеличенном кубе найдено, теперь этот куб становится начальным приближением:
xt = x1 _ new, y1 = y1 _ new, zt = zt _ new .
И повторяются шаги 1 — 3 до тех пор, пока расчетная область не увеличится до размеров исходной области:
x1 e [0; Nx], y1 e [0; Ny], z1 e [0; Nz].
Второй шаг оптимизации — это разгрузка оперативной памяти
На первом шаге оптимизации размер «действующего списка» был уменьшен. Для сокращения размера остальных массивов предлагается следующая последовательность действий:
1) определить некоторую, не слишком маленькую глубину Nzmin, такую, что Nzmin < Nz. Например, Nzmin = Nz /6 . Выбрать число M, равное примерно Nzmin /3 или Nzmin /4 ^
2) задать размер матриц C, T, x_log равным NxNzNzmin ;
3) считать из файла, содержащего скоростную модель, данные, соответствующие координатам x e [0; Nx ], y e [0; Ny ], z e [0; Nz ];
4) найти решение задачи в этой области;
5) сохранить значения матрицы T в файл;
6) для матриц С и Т переместить значения Cxy,z = Cx, N _M+z, Tx,y,z = Tx,y,Nzmm -M+z для X e [°; N ], y e [0; Ny ], z e [0; M].
Остальные значения матрицы T (то есть для x e[0; Nx ], y e[0; Ny ], z e [ M +1; Nzmin ]) сделать равными значениям по умолчанию, а для матрицы С - считать из файла значения, соответствующие координатам x e [0; Nx ], y e [0; Ny ], z e[#Zmin +1,2* NzKan -M];
7) сделать всю границу z = M границей источников и подготовить область для решения в соответствии с первым шагом оптимизации, где размер расчетной области будет увеличиваться только по оси Z (так как по остальным осям область максимальна);
8) повторять шаги 4) — 8) до тех пор, пока не будет исчерпана вся исходная расчетная область.
Таким образом, массивы С, T, x_log теперь имеют размер NxNyNz /6 (в зависимости от возможностей компьютера, размер может
быть увеличен или уменьшен). Массив L размера 113 • N3 • 3, где N ^ 7.
Повторим произведенные ранее расчеты. Теперь для решения уравнения (1) в области размером 801-801-601 в оперативной памяти нужно хранить: 801-801-601/(6-4-2) байт (массивы C и T) + 801-801-601/(6-1) байт (массив x_log) + 77-77-77-3-4 байт (массив L), что составляет 583881697,5 байт и 0,5 Гбайт. Это значение вполне реально для обычного персонального компьютера.
Если заменить тип float на тип double для повышения точности вычислений, в оперативной памяти будет храниться примерно 1,03 Гб.
158
Численные эксперименты
Эксперимент с моделью размером 801-801-601 не может быть рассчитан для неоптимизированной задачи из-за нехватки памяти. Для оптимизированной задачи модель такого размера вычисляется в течение часа. Далее приведен эксперимент с задачей меньшего размера — 401-401-300 с точечным источником в координатах (201, 201, 0).
На рисунке 5 изображено сечение УЪ при х = 201 (значение в отсчетах). Цветом изображена скорость в м/ с.
50 100 150 200 25) 300 Э50 «О
Рис. 5. Скоростная модель
Рисунок 6, а — это результат неоптимизированного решения, а рисунок 6, б — оптимизированного (изображение в том же сечении, что и на рисунке 5).
Относительная погрешность между этими вычислениями равна примерно 0,34%.
Время вычисления неоптимизированной задачи составляет 3 часа 45 секунд. Время вычисления оптимизированной задачи мало — 10 минут 49 секунд. Прирост во времени выполнения задачи большой — примерно в 16,7 раза.
Рис. 6. Результаты решения
Список литературы
1. Jeong W.-K., Whitaker R. T. A Fast Iterative Method for Eikonal Equations. // SIAM Journal on Scientific Computing. 2008. Vol. 30, № 5. P. 2512-2534.
Об авторе
Валерия Александровна Седайкина — науч. сотр., Балтийский федеральный университет им. И. Канта, Калининград. Email: VSedaikina@kantiana.ru
About the author
Valeria Sedaikina — researcher, I. Kant Baltic Federal University, Kaliningrad. Email: VSedaikina@kantiana.ru
