Акустический журнал, 2023, T. 69, № 1, стр. 13-21

ВВЕДЕНИЕ

Численный эксперимент на сегодняшний день является неотъемлемой частью развития методов неинвазивной хирургии с использованием нелинейных ультразвуковых волн. Одним из основных объектов, для исследования которого активно используются численные методы, являются нелинейные фокусированные ультразвуковые пучки высокой интенсивности, с помощью которых производится разрушение задаваемых структур внутри тела человека, например, опухолей [1]. Для генерации мощного ультразвука используются преобразователи различной формы и конструкции, при проектировании которых для конкретных клинических приложений необходимо уметь количественно описывать структуру создаваемых ими акустических полей (рис. 1). Такая задача может быть эффективно реализована при помощи методов численного моделирования [2–4]. Одной из наиболее полных волновых моделей для теоретического описания мощных ультразвуковых пучков является уравнение Вестервельта, которое позволяет количественно точно описывать нелинейные ударно-волновые поля, создаваемые фокусированными преобразователями мощного ультразвука в однородных поглощающих средах [5, 6]. При постановке задачи обычно считается, что излучателем генерируется монохроматическая волна, спектр которой в процессе распространения обогащается высшими гармониками за счет эффекта акустической нелинейности. В общем случае для решения данного уравнения, когда эволюционной переменной выступает время, необходимо использование суперкомпьютерных мощностей [7]. Однако ряд упрощений позволяет существенно снизить вычислительные затраты.

Рис. 1.

Схема распространения сфокусированного ультразвукового пучка, создаваемого HIFU излучателем.

Для оптимизации численного моделирования нелинейных фокусированных ультразвуковых пучков создаются специальные алгоритмы. В данной работе рассматривается важный для практики случай направленного трехмерного пучка, что позволяет упростить вычислительную задачу путем перехода в бегущую систему координат, ось которой ориентирована вдоль преимущественного направления распространения волн в пучке. В этом случае эволюционной переменной является координата вдоль оси пучка z (рис. 1). Численное решение эволюционного уравнения такого типа обычно строится с использованием метода расщепления по физическим факторам [8], согласно которому на каждом шаге сетки по оси z каждый физический эффект рассчитывается отдельно с помощью наиболее подходящего численного метода. Нелинейное эволюционное уравнение Вестервельта в такой постановке решалось для излучателей различной формы и размера [6, 9, 10]. При этом типовой размер матриц для хранения параметров поля давления составляет M_x = 10 000 на M_y = 10 000 (рис. 1) для каждой из N = 1000 спектральных составляющих профиля нелинейной волны [9]. В этом случае требуемый объем оперативной памяти, который достаточно велик (несколько сотен Гб), и соответствующая сложность расчетов все еще требуют использования суперкомпьютерных мощностей.

Однако, если учесть особенности пространственной структуры нелинейных полей сильно сфокусированных излучателей, то рассматриваемая задача может быть решена на обычных персональных компьютерах с многоядерными центральными процессорами (CPU). Поскольку нелинейные эффекты в основном проявляются там, где амплитуда волны велика, т.е. вблизи фокуса, то большое количество гармоник необходимо учитывать именно в этой сравнительно небольшой области пространства [9]. Такой способ позволяет на порядок снизить вычислительные затраты и требования к объему оперативной памяти компьютера, однако даже при параллельном исполнении вычислений на нескольких процессорных ядрах (обычно от 2 до 16), позволяющем существенно ускорить расчеты, моделирование поля для одного набора входных параметров может длиться до нескольких суток [9]. С учетом того, что для характеризации поля одного излучателя во всем диапазоне рабочих мощностей требуется выполнить несколько десятков расчетов, такая скорость моделирования не является удовлетворительной.

Потенциальным способом решения проблемы скорости вычислений является быстро развивающаяся в последнее время технология параллельного программирования на графических процессорах (GPU). В отличие от CPU данные процессоры имеют до несколько тысяч узкоспециализированных ядер, способных выполнять широкий спектр математических операций [11]. Такое количество ядер позволяет увеличить скорость расчетов за счет большого количества одновременно запускаемых процессов, которые параллельно обрабатывают однотипные сегменты данных. Получаемое ускорение вычислений необязательно пропорционально числу обрабатывающих потоков, что обусловлено меньшей мощностью, производительностью и скоростью вычисления ядер GPU по сравнению с ядрами CPU, а также особенностями хранения и передачи данных между оперативной памятью CPU и графической памятью. Так как хранение и запись данных при выполнении каждого из шагов алгоритма вдоль эволюционной координаты производится в оперативной памяти CPU, дополнительное время уходит на обмен данными между центральным и графическим процессорами на каждом шаге алгоритма. В силу рассматриваемых особенностей архитектуры GPU важной задачей становится нахождение баланса между максимальными размерами обрабатываемых на данном шаге алгоритма данных, которые помещаются в оперативную память графического процессора во избежание их передачи по частям, и минимизацией обмена данными между процессорами.

Ранее для данной задачи был реализован алгоритм на GPU, в котором были написаны функции-ядра для вычисления оператора дифракции методом углового спектра, нахождения точного решения для оператора поглощения и решения оператора нелинейности в спектральном представлении методом Рунге–Кутты 4 порядка [12]. Однако при использовании спектрального метода количество вычислительных операций и, соответственно, время вычислений пропорциональны квадрату количества образующихся гармоник. Поэтому данный метод будет эффективен только при слабом проявлении нелинейных эффектов, когда в спектре волны образуется небольшое число высших гармоник. Оператор нелинейности можно также вычислять с помощью различных методов во временном представлении, однако во многих из них для описания области разрыва требуется большое количество узлов сетки по оси времени (50–100), что увеличивает общее количество точек сетки и, таким образом, количество обрабатываемых данных [13, 14]. Для решения данной проблемы эффективнее использовать удароулавливающие схемы, позволяющие моделировать образование ударного фронта с использованием небольшого количества точек на нем (около 3-х), например, консервативную схему типа Годунова [15].

Целью данной работы является реализация алгоритма удароулавливающей схемы типа Годунова на графическом процессоре (GPU) для большего ускорения параллельных вычислений по сравнению с использованием многоядерных центральных процессоров (CPU) и обеспечения возможности моделирования мощных ультразвуковых пучков с ударными фронтами на обычном персональном компьютере.

ТЕОРЕТИЧЕСКАЯ МОДЕЛЬ

Нелинейное уравнение Вестервельта

Уравнение Вестервельта в бегущей системе координат можно записать в эволюционной форме как:

(1)

$\begin{gathered} \frac{{{{\partial }^{2}}p}}{{\partial \tau \partial z}} = \frac{{{{c}_{0}}}}{2}\left( {\frac{{{{\partial }^{2}}p}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}p}}{{\partial {{y}^{2}}}} + \frac{{{{\partial }^{2}}p}}{{\partial {{z}^{2}}}}} \right) + \\ + \,\,\frac{\beta }{{2{{\rho }_{0}}c_{0}^{3}}}\frac{{{{\partial }^{2}}{{p}^{2}}}}{{\partial {{\tau }^{2}}}} + \frac{\delta }{{2c_{0}^{3}}}\frac{{{{\partial }^{3}}p}}{{\partial {{\tau }^{3}}}}, \\ \end{gathered} $

где p(x, y, z, τ) – акустическое давление,${\text{\;}}{{c}_{0}}$ – скорость звука в среде, $\tau = t - {z \mathord{\left/ {\vphantom {z {{{c}_{0}}}}} \right. \kern-0em} {{{c}_{0}}}}$ – время в бегущей системе координат, β – коэффициент нелинейности, $\delta $ – коэффициент термовязкого поглощения [5]. Дифференциальные операторы в правой части уравнения в порядке слева направо описывают эффекты дифракции, нелинейности и термовязкого поглощения.

При численном решении уравнения (1) на каждом шаге вдоль оси z дискретизованное поле давления $p\left( {x,y,z,\tau } \right)$ представляется в памяти ЭВМ в виде трехмерной матрицы, которая содержит набор комплексных амплитуд N гармоник спектра волны в разложении ее профиля в конечный ряд Фурье в каждой пространственной точке плоскости xy на равномерной сетке с числом точек M_x и шагом Δx по оси x и M_y с шагом Δy по оси y:

(2)

$p\left( {x,y,z,\tau } \right) = \mathop \sum \limits_{n = - N}^N {{p}_{n}}\left( {x,y,z} \right){{e}^{{ - i{{\omega }_{n}}\tau }}}.$

Здесь ${{\omega }_{n}} = \omega n$ – круговые частоты гармоник с номером n, ω – круговая частота монохроматического источника, ${{p}_{n}}$ – комплексные амплитуды гармоник. Отметим, что в памяти ЭВМ достаточно хранить только одну половину спектра положительных частот, т.к. амплитуда второй половины отрицательных частот является комплексно-сопряженной к первой.

Используя спектральное представление (2) для поля давления в схеме расщепления по физическим факторам решения уравнения Вестервельта (1), можно вычислить оператор дифракции методом углового спектра и точно рассчитать оператор поглощения для каждой из спектральных компонент волны независимо друг от друга [12]. Особенности параллельных вычислений оператора нелинейности в спектральном и временном представлениях рассмотрены ниже. Рассматривается реализация только нелинейного оператора в уравнении Вестервельта, т.е. фактически решается уравнение простых волн. Так как данное уравнение решается независимо в каждой пространственной координате, для упрощения представления данных в памяти ЭВМ перейдем к рассмотрению строки исходного массива по пространственной координате x с M элементами, что не повлияет на исследование эффективности работы алгоритма. Таким образом, волна распространяется вдоль оси z, нелинейные эффекты рассчитываются на каждом шаге по z независимым образом M раз для всех узлов сетки по координатe x (рис. 2), оптимизация расчетов проводится путем распараллеливания потоков данных по этой координате. В дальнейшем при моделировании полной нелинейно-дифракционной задачи расчеты будут проводиться на трехмерной сетке. Для этого данные с двумерной пространственной сетки в плоскости xy с полным количеством узлов M = M_xM_y можно обработать, разбивая их на описанные выше одномерные массивы размера M.

Рис. 2.

Схема реализации параллельного расчета нелинейного оператора, где M – количество пространственных точек.

Метод Годунова

Численное решение уравнения (4) на каждом шаге сетки по эволюционной координате σ строилось с помощью схемы типа Годунова (рис. 3). Алгоритм расчета представляет собой явную шеститочечную консервативную схему второго порядка точности по времени и первого порядка точности по координате распространения [13, 15]:

(6)

$V_{j}^{{k + 1}} = V_{j}^{k} - \frac{{{{h}_{\sigma }}}}{{{{h}_{\theta }}}}\left( {H_{{j + \frac{1}{2}}}^{k}\left( \sigma \right) - H_{{j - \frac{1}{2}}}^{k}\left( \sigma \right)} \right),$

где потоки $H_{{j + \frac{1}{2}}}^{k}$ через центры ячеек численной сетки на k-том шаге по координате $~\sigma $ задаются следующим образом

(7)

$\begin{gathered} H_{{j + \frac{1}{2}}}^{k}\left( \sigma \right) = - \frac{1}{4}\left( {{{{\left( {V_{{j + \frac{1}{2}}}^{ + }} \right)}}^{2}} + {{{\left( {V_{{j + \frac{1}{2}}}^{ - }} \right)}}^{2}}} \right) - \\ - \,\,\frac{{a_{{j + \frac{1}{2}}}^{k}\left( \sigma \right)}}{2}\left[ {V_{{j + \frac{1}{2}}}^{ + }\left( \sigma \right) - V_{{j + \frac{1}{2}}}^{ - }\left( \sigma \right)} \right]. \\ \end{gathered} $

Рис. 3.

Шаблон для консервативной численной схемы типа Годунова для уравнения простых волн.

Здесь $~{{h}_{\sigma }}$ – шаг схемы вдоль координаты σ, $~{{h}_{\theta }}$ – шаг сетки по оси безразмерного времени ${{\theta }}$, и локальная скорость потока в сеточной ячейке:

(8)

$a_{{j + \frac{1}{2}}}^{k}\left( \theta \right) = \max \left| {V_{{j + \frac{1}{2}}}^{ - },V_{{j + \frac{1}{2}}}^{ + }} \right|.$

Для получения второго порядка точности в данном алгоритме используется кусочно-линейная реконструкция значений давления $V\left( {{{\sigma }_{k}},~{{\theta }_{j}}} \right)$ справа $V_{{j + \frac{1}{2}}}^{ + }\left( \sigma \right)$ и слева $V_{{j + \frac{1}{2}}}^{ - }\left( \sigma \right)$ от узла численной сетки $\left( {k,~j} \right)$:

(9)

$\begin{gathered} V_{{j + \frac{1}{2}}}^{ + } = V_{{j + 1}}^{k}(\sigma ) - \frac{{{{h}_{\theta }}}}{2}\left( {\frac{{\partial V}}{{\partial \theta }}} \right)_{{j + 1}}^{k}, \\ V_{{j + \frac{1}{2}}}^{ - } = V_{j}^{k}(\sigma ) - \frac{{{{h}_{\theta }}}}{2}\left( {\frac{{\partial V}}{{\partial \theta }}} \right)_{j}^{k}. \\ \end{gathered} $

Для большей устойчивости численного алгоритма, задействованные в (8) производные решения по времени на шаге k по координате σ выбираются таким образом, чтобы их значения были минимальными по модулю из возможных значений производных – правой, левой и центральной с весовым коэффициентом 1 ≤ $~b$ ≤ 2:

(10)

$\begin{gathered} \left( {\frac{{\partial V}}{{\partial \theta }}} \right)_{j}^{k} = \\ = \,\min \bmod ~\left( {\frac{{b\left( {V_{j}^{k}\, - \,V_{{j - 1}}^{k}} \right)}}{{{{h}_{\theta }}}},\frac{{V_{{j + 1}}^{k}\, - \,V_{{j - 1}}^{k}}}{{2{{h}_{\theta }}}},\frac{{b\left( {V_{{j + 1}}^{k}\, - \,V_{j}^{k}} \right)}}{{{{h}_{\theta }}}}} \right). \\ \end{gathered} $

Весовой коэффициент $b = 2$ соответствует наиболее точному решению с минимальным сеточным поглощением, а $b = 1$ – более устойчивому численному алгоритму, что достигается за счет увеличения сеточного поглощения. В данной работе величина b полагалась равной единице. Представленная схема позволяет рассчитывать распространение узких ударных фронтов с большой точностью при использовании всего лишь трех узлов сетки на ударном фронте [15]. Отметим, что при решении трехмерных задач обычно дополнительно вводится искусственное поглощение, которое позволяет размыть ширину ударного фронта до нужных значений (обычно 7–8 точек на фронт) с целью уменьшения больших пространственных градиентов поля давления по поперечным координатам [6].

Расчеты проводились на двух различных компьютерах: на ноутбуке с CPU Intel Core i5-8250U с видеокартой GPU Nvidia GeForce MX150 (380 ядер, 1.5 ГГц), и на персональном компьютере (ПК) с CPU Intel Core i7 4790 с видеокартой Nvidia GTX1070 (1980 ядра, 1.7 ГГц). Для реализации вычислений на графических процессорах была написана программа на языке C, содержащая функцию-ядро CUDA, которая реализовала алгоритмы (6)–(10). Для определения эффективности работы данной программы на языке С были реализованы однопоточный аналогичный алгоритм для CPU, а также алгоритм для GPU на основе спектрального метода, со скоростью вычислений которых и производилось сравнение. При численном решении уравнений (6)–(10) на каждом шаге по σ дискретизованное поле давления представлялось в памяти ЭВМ в виде двумерной матрицы, которая содержала набор M профилей волн, заданных в M пространственных точках. Такое представление данных необходимо для реализации распараллеливания по данным, при котором каждый профиль волны из числа M обрабатывается параллельно с другими профилями различными ядрами графического процессора. Вычисления проводились для различного количества пространственных координат M (100–10 000 точек), а также для различного числа точек на периоде волны (100–3000 точек) с шагами временной сетки $\Delta t = 0.0063,~$ $\Delta \sigma = 0.2378$, удовлетворяющими устойчивости схемы типа Годунова при максимальном рассматриваемом количестве точек на период волны. Для определения правильности реализованного алгоритма во временном представлении использовалось аналитическое решение уравнения простых волн с бесконечно узким ударным фронтом (4) с начальным профилем (5), представленное в неявном виде [17]:

(11)

${{V}_{{{\text{th}}}}} = \sin \left( {\theta + \sigma {{V}_{{{\text{th}}}}}} \right).$

Положение и величина ударного фронта определялись по правилу равенства площадей.

РЕЗУЛЬТАТЫ

На рис. 4 представлены профили волны на различных расстояниях, полученные численно с использованием метода Годунова, реализованного на GPU ноутбука, и теоретическое решение, описываемое уравнением (11). Наблюдается хорошее согласие численного и аналитического решений с ошибкой не более 0.03% в амплитуде ударного фронта. Расхождение наблюдается только в области разрыва, так как, в отличие от аналитического решения, численная схема приводит к формированию ударного фронта конечной ширины с тремя точками на разрыв. Таким образом, ширина фронта при использовании 1000 точек на период составляет около 0.2% от периода волны. Данные результаты подтверждают правильность реализации алгоритма для GPU. Идентичные результаты были получены для алгоритмов, реализованных на ПК с использованием одного из ядер CPU и на ноутбуке с использованием и без использования видеокарты.

Рис. 4.

Сравнение аналитического (серая сплошная линия) и численного (черная пунктирная линия) решений для плоских волн с гармоническим начальным профилем на разных расстояниях от излучателя: (а) – σ = 0.5, (б) – σ = 1.0, (в) – σ = 4.0. (г) – Сравнение аналитического и численного решений вблизи ударного фронта при σ = 4.0.

В табл. 1 и 2 сравниваются результаты вычислений на ноутбуке по времени работы однопоточного алгоритма, основанного на временном методе и реализованного на CPU, с аналогичным алгоритмом для GPU, а также с алгоритмом, реализующим спектральный метод на GPU, в зависимости от количества пространственных точек в буфере M и количества гармоник N. Как и ожидалось, время вычислений с использованием спектрального метода на GPU увеличивается квадратично с увеличением числа гармоник, что демонстрирует его неэффективность в случае сильных нелинейных эффектов и большого числа гармоник. При этом время вычислений, проведенных на основе метода Годунова, пропорционально количеству гармоник. Для размеров матрицы поля давления с 10 000 пространственных точек и при 1000 гармоник расчеты, реализованные на GPU с использованием временного подхода, оказываются в 26 раз быстрее по сравнению с CPU и в 2.6 раз быстрее по сравнению со спектральным алгоритмом, реализованном на GPU.

В табл. 3 и 4 представлено сравнение скорости работы тех же программ на более мощном компьютере (ПК). Наблюдаются аналогичные зависимости времени, затрачиваемого на вычисления, от числа используемых гармоник для обоих методов. Для размеров матрицы с 10 000 пространственных точек и 1000 гармоник расчеты, реализованные на GPU с использованием временного подхода, в 33 раза быстрее по сравнению с CPU и в 2.3 раза быстрее по сравнению со спектральным алгоритмом, реализованном на GPU.

На рис. 5, где данные по сравнению скорости работы алгоритмов сведены в виде графиков, видно, что ускорение для ноутбука при увеличении как количества пространственных точек, так и числа гармоник (рис. 6) достигает определенного значения и далее не изменяется. Это обусловлено тем, что графический процессор на ноутбуке имеет небольшое количество ядер, способных выполнять параллельные процессы, а также небольшой объем памяти, ограничивающий вместимость обрабатываемых данных. Поэтому при больших параметрах входного массива возникает необходимость передавать его из памяти CPU в память GPU по частям, что занимает дополнительное время. Также небольшое количество ядер дает возможность обрабатывать лишь соответствующее число потоков одновременно. Ускорение вычислений, реализованных на более мощном компьютере (ПК), растет для всех видов построенных зависимостей, что обусловлено большей мощностью и объемом памяти его GPU.

Рис. 5.

Зависимость ускорения работы алгоритма на основе метода Годунова, реализованного на GPU: (а) – по сравнению с версией для CPU, (б) – по сравнению со спектральным алгоритмом на GPU в зависимости от размера буфера M.

Рис. 6.

Зависимость ускорения работы алгоритма на основе метода Годунова, реализованного на GPU: (а) – по сравнению с версией для CPU, (б) – по сравнению со спектральным алгоритмом на GPU в зависимости от числа гармоник N.

В случае сравнения скорости работы алгоритмов, реализованных с использованием метода Годунова на CPU и GPU, достигаемое ускорение на ПК меньше, чем на ноутбуке. Это обусловлено большей мощностью центрального процессора ПК. Важным результатом является значительное ускорение, достигаемое в расчетах на ПК при увеличении числа гармоник в случае использования временного метода по сравнению со спектральным методом.

ЗАКЛЮЧЕНИЕ

В работе реализован численный алгоритм расчета нелинейных акустических эффектов на графических процессорах (GPU) с использованием удароулавливающей схемы типа Годунова. Расчеты на GPU позволили ускорить вычисления по сравнению с алгоритмом для CPU на порядок и более как для ноутбука, так и для более мощного персонального компьютера для всех рассматриваемых размеров массива данных. Кроме того, для характерного числа гармоник (N = 1000), необходимого для описания высокоамплитудных ударных фронтов, реализация схемы Годунова на GPU дала выигрыш по времени на порядок по сравнению со спектральным алгоритмом в расчетах на персональном компьютере и в 5 раз на ноутбуке. Таким образом, было показано, что проведение вычислений на графических процессорах может эффективно применяться для решения нелинейных волновых задач с использованием удароулавливающих численных схем, работающих во временном представлении.

Работа выполнена при поддержке гранта РНФ 20-12-00145.