Зависимость расчетных высот волн цунами от разрешения сетки

Полный текст

1. Введение Важность численного моделирования волн цунами, начиная с момента их возникновения и заканчивая оценкой зоны затопления на берегу, чрезвычайно высока из-за их возможного воздействия на прибрежные территории. Например, волна цунами, вызванная Великим землетрясением Тохоку 11 марта 2011 года, привела к человеческим жертвам (по оценкам, более 20 тысяч человек) и огромному экономическому ущербу (около 1/4 ВВП Японии [18]). Если эпицентр землетрясения расположен в районе одной из зон субдукции у берегов Японии, то волна цунами, вызванная таким землетрясением, достигает ближайшего к очагу прибрежного района примерно за 20 минут. Поэтому предлагаемые инструменты анализа и оценки параметров волн должны давать результаты очень быстро, в течение 1-2 минут. В то же время оценка ожидаемых параметров волны у побережья (главным из которых является максимальная высота волны) должна быть достаточно точной. Очевидно, что значительно заниженное значение ожидаемой амплитуды может привести к человеческим жертвам. В случае значительного завышения максимальных значений стоимость мер по эвакуации, остановке предприятий и дорожного движения может привести к большим экономическим потерям. Существует несколько доступных пакетов программного обеспечения для численного моделирования цунами. Среди наиболее известных мы упоминаем следующие: метод расщепления цунами MOST (Method Of Splitting Tsunamis, Тихоокеанская лаборатория морской среды NOAA, Сиэтл, Вашингтон, США) [1, 12]; COMCOT (Корнеллский университет, Итака, Нью-Йорк, США; NS Science, Новая Зеландия) [13]; TUNAMI-N1/TUNAMI-N2 (Университет Тохоку, Сендай, Япония) [10, 14]; «ВОЛНА» и «НАМИ-ДАНС» (программное обеспечение для моделирования цунами, разработанное Специальным бюро Сахалина, Россия, и METU, Турция) [15]. Для получения результатов в кратчайшие сроки некоторые алгоритмы реализованы на графических процессорах (GPU) и высокопроизводительных кластерах. Однако ни один из этих пакетов не дает результата за требуемый промежуток времени в 1-2 минуты. Как показано в [3-5], использование аппаратного ускорения выполнения компьютерного кода (калькулятор на основе FPGA для системы уравнений мелкой воды) позволяет выполнять корректное (т. е. практически точное) численное моделирование двух часов распространение волны на сетке 2500 × 3000 узлов менее чем за 1 минуту. Время, необходимое для численного расчета распространения волн на компьютере, зависит от следующих факторов: (1) характеристики вычислительной платформы, (2) количество узлов сетки, (3) общее время распространения цунами по акватории и (4) шаг по времени. Таким образом, время вычислений можно сократить либо за счет увеличения производительности компьютера, либо за счет уменьшения количества вычислительных узлов благодаря увеличению шага пространственной сетки. Произвольное увеличение значения шага по времени невозможно из-за возникновения вычислительной нестабильности. Существует максимально допустимое значение этого шага, зависящее от максимальной глубины расчетной области и длины пространственных шагов [6]. Численное моделирование распространения волн цунами обычно (в частности, как это сделано в упомянутых выше программных комплексах) осуществляется с использованием разностных схем на сетках, чаще всего с привязкой к географическим координатам. Разрешение расчетных сеток обычно выбирается исходя из требуемого уровня детализации результатов (высот цунами) вдоль береговой линии. Если высокий уровень детализации не требуется, длину пространственного шага можно выбирать произвольно. Небольшой шаг сетки существенно увеличивает количество вычислительных узлов и, следовательно, время расчета одного временного шага. Кроме того, по условию устойчивости Куранта для явных разностных схем [6] короткий пространственный шаг требует меньшего шага по времени, что также увеличивает время вычислений. Для точной оценки ожидаемой максимальной высоты волн вдоль побережья обычно требуется шаг расчетной сетки от 5 до 15 метров. Поэтому для того, чтобы иметь возможность выполнить расчеты в разумные сроки, в ряде исследований [7, 9] используется метод вложенных (адаптивных) сеток. Возникает вопрос о влиянии длины пространственного шага на результаты моделирования. Это важно как для анализа результатов численного моделирования, так и для правильного выбора регионов в методе вложенных сеток. В настоящей работе расчеты распространения цунами от одного и того же источника выполнены на трех сетках различного разрешения. В расчетах используются реальные батиметрические данные на шельфе Японии. Остальная часть статьи устроена следующим образом. В разделе 2 сначала представлена система уравнений мелкой воды в рассматриваемой форме вместе с конечно-разностной схемой Мак-Кормака, используемой для численных расчетов. Затем кратко описывается архитектура предлагаемого калькулятора на базе FPGA. Результаты численных испытаний приведены в разделе 3. Они обсуждаются в заключительном разделе. 2. Материалы и методы 2.1. Системы уравнений мелкой воды. Следуя большинству исследователей явлений цунами (см., напр., [2, 10], мы моделируем распространение волн на основе версии системы уравнений мелкой воды [11] в следующем виде: ∂H ∂(uH) ∂(vH) + + = 0, ∂t ∂x ∂y ∂u ∂u ∂u ∂H ∂D + u + v + g = g , (2.1) ∂t ∂x ∂y ∂x ∂x ∂v ∂v ∂v ∂H ∂D + u + v + g = g , ∂t ∂x ∂y ∂y ∂y где H(x,y,t) = η(x,y,t)+D(x,y) - расстояние от поверхности моря до дна; η - смещение морской поверхности относительно среднего уровня моря (высота волны); D(x,y) - глубина (значение которой предполагается известной во всех точках сетки); u(x,y,t) и v(x,y,t) -компоненты вектора скорости водного потока в направлениях x и y соответственно; g - ускорение силы тяжести. Система уравнений мелкой воды (2.1) не учитывает дисперсионные эффекты, поскольку при локальном распространении достаточно длинных волн цунами эти эффекты не успевают существенно изменить параметры волн вблизи побережья. Использование аппроксимации уравнений мелкой воды для изучения явления типично среди исследователей цунами. Существует большой список литературы по численному изучению нелинейных систем дифференциальных уравнений в частных производных гиперболического типа, к которым также принадлежит система (2.1). В вышеупомянутом программном обеспечении MOST реализовано разделение шагов расчета по пространственным направлениям. 2.2. Схема Мак-Кормака. Для численной аппроксимации системы уравнений мелкой воды (2.1) использовалась двухполушаговая схема Мак-Кормака (см. [8]). Этот алгоритм представляет собой явную разностную схему, использующую трехточечный шаблон типа «крест». Расчет на каждой итерации проводится в два этапа. Таким образом, на каждой итерации значения высоты волны и компоненты вектора скорости вычисляются с использованием необходимых данных в той же точке и двух соседних точках сетки при предыдущем значении времени; описание шаблонов расчета см. в [3]. Ниже представлена версия схемы Мак-Кормака, используемая для численных экспериментов. Шаг 1: Hˆijn+1 - Hijn Hijn unij - Hin-1j uni-1j Hijn unij - Hijn -1 unij-1 τ Δx Δy uˆnij+1 - unij + un unij - uin-1j + vijn uijn - uijn-1 + g ηijn - ηin-1j = 0, τ ij Δx Δy Δx vˆijn+1 - vijn + un vijn - vin-1j + vijn vijn - vijn -1 + g ηijn - ηijn -1 = 0. τ ij Δx Δy Δy + + = 0, Шаг 2: Hijn+1 - Hˆijn+12+Hijn Hˆin+1+1j uˆni+1+1j - Hˆijn+1 uˆijn+1 Hˆijn+1+1 vˆijn+1+1 - Hˆijn+1 vˆijn+1 + + = 0, τ/2 Δx Δy - uˆnij+1+unij n+1 - uˆnij+1 n uˆnij+1+1 - uˆnij+1 + g ηˆin+1+1j - ηˆijn+1 = 0, unij+1 2 n uˆi+1j + vij + uij τ/2 Δx Δy Δx vijn+1 - vˆijn+12+vijn n vˆin+1+1j - vˆijn+1 + vijn vˆijn+1+1 - vˆijn+1 + g ηˆin+1+1j - ηˆijn+1 = 0. + uij τ/2 Δx Δy Δy Здесь Hijn , unij, vijn и ηijn - сеточные переменные, соответствующие функциям H, u и v в дифференциальной системе (2.1). Параметры τ, Δx и Δy - это временной и пространственные шаги расчетной сетки. Для учета сферической формы Земного шара использовался уменьшающийся шаг сетки по долготе для больших значений широты. Обозначение описывает переменные на временном слое n, описывает промежуточные значения, а соответствует переменным на временном слое n + 1. 2.3. Калькулятор на базе FPGA. Конечно-разностная схема Мак-Кормака, описанная в предыдущем подразделе, была реализована как конвейер, позволяющий параллельно вычислять искомые параметры волны на нескольких временных слоях. Для реализации алгоритма Мак-Кормака (конечно-разностной схемы) на платформе FPGA [4, 5] в качестве базового элемента калькулятора (ускорителя численного решения системы уравнений мелкой воды) были спроектированы так называемые процессорные элементы (ПЭ). Один ПЭ отражает реализацию одного временного шага алгоритма расчета. ПЭ были реализованы с использованием технологии HLS (High Level Synthesize) [19]. На вход ПЭ поступает поток значений H, u, v и D на n-м шаге по времени, что представляет собой последовательный обход вычислительного стека. На выходе получается тот же поток со значениями на временном шаге n + 1. ПЭ работает в конвейерном режиме и позволяет обрабатывать одну точку сетки за один такт компьютера. Соединив выход ПЭ со входом другого ПЭ, можно организовать цепочки (конвейеры) разной длины, позволяющие рассчитывать несколько шагов алгоритма одновременно по мере поступления результатов от предыдущего элемента. Проведем краткое сравнение ключевых параметров FPGA и процессора x86, который де-факто является отраслевым стандартом для вычислительных платформ. Количество рабочих блоков FPGA достигает 12888 при тактовой частоте до 1 ГГц. Процессор x86 имеет менее 100 рабочих блоков при тактовой частоте до 5 ГГц. Пятикратное преимущество перед FPGA в тактовой частоте не компенсируется меньшим более чем в 100 раз количеством рабочих блоков. Еще более серьезное преимущество платформы FPGA заключается в объеме памяти и допустимом количестве обращений к памяти за один такт компьютера. На современных кристаллах FPGA объем памяти достигает 100 Мб при количестве операций чтения/записи за один такт компьютера, равном удвоенному числу рабочих блоков, т. е. достигает более 25000. Процессор x86 с аналогичным объёмом памяти в 100 Мб допускает только одну операцию чтения/записи за такт компьютера. Загрузить все рабочие блоки в компьютер можно только в том случае, если в качестве источника/приемника обрабатываемых данных используется регистровый файл. Он встроен в архитектуру x86. Практика показывает, что при реализации сложных алгоритмов потоковой обработки данных допустимого размера регистрового файла недостаточно. В таких случаях узким местом является путь доступа к памяти (в том числе к кэшу), вплоть до того, что часть рабочих блоков находится в состоянии простоя из-за ожидания получения данных из памяти или ожидания записи обработанных данных в память. В случае использования FPGA эта проблема решается тем, что для каждого рабочего блока создается независимый блок памяти, способный выполнять две операции чтения/записи за один такт компьютера. Это позволяет практически всем операционным блокам загружать новые операции в каждом такте. Общий объем внутренней памяти FPGA сравним с размером кэша процессора x86 и позволяет реализовать достаточно сложные алгоритмы обработки без обращения к внешней памяти. Ограничивающим фактором использования FPGA вместо центрального процессора (ЦП) является сравнительная ограниченность таких инструментов, как среда разработки и реализация алгоритмов обработки данных. Разработка приложений для FPGA ведется на языках описания аппаратуры, таких как Verilog или VHDL, что требует, с одной стороны, ценных знаний в области проектирования цифровых схем, а с другой - дополнительного времени на разработку, поскольку языки описания аппаратуры являются достаточно низкоуровневыми. В настоящее время широкое распространение получил подход HLS (High Level Synthesize): синтез цифровых схем из компьютерных языков высокого уровня. Использование компиляторов из высокоуровневых языков позволяет значительно ускорить процесс разработки. Таким образом, описанный выше калькулятор был реализован одним разработчиком на языке программирования C++ за 1 в месяц. Рис. 1. Архитектура калькулятора для ускорения численного решения системы уравнений мелкой воды. Fig. 1. Calculator architecture for accelerating the numerical solution of a system of shallow water equations. 3. Численные результаты Чтобы исследовать влияние разрешения расчетной сетки на численные результаты, были проведены три вычислительных эксперимента по моделированию распространения цунами от одного и того же модельного источника, расположенного в Японском море у северо-западного побережья острова Хонсю. Для этой цели были построены 3 цифровые батиметрии на регулярных сетках, привязанных к географическим координатам и покрывающих одну и ту же территорию от 137◦ до 141◦ восточной долготы и от 37◦ до 42◦ северной широты. Цифровая батиметрия с координатной сеткой была построена на основе данных JODC [16] и ETOPO-2022 [17]. Длина пространственного шага была выбрана одинаковой в обоих направлениях и равна 0,001 географического градуса (111,3 м в направлении юг-север) для наиболее детальной батиметрии. Для батиметрии со средней и крупной сеткой длина пространственного шага в обоих направлениях составила 0,002◦ (222,6 м) и 0,004◦ (445,2 м). География и топография расчетной области показаны на рис. 2. Рис. 2. География и батиметрия расчетного района. Fig. 2. Geography and bathymetry of the calculated area. Начальное вертикальное смещение водной поверхности имело форму эллипса с центром в точке с координатами x0 = 138◦ восточной долготы, y0 = 40◦ северной широты и определялось следующей формулой: , (3.1) . Здесь r1 и r2 задают длины полуосей эллипса, а η0 -высоту в центре источника. В расчетах r1 и r2 были равны 35 км и 90 км соответственно, а максимальная высота смещения η0 в центре источника принималась равной 1 м. Длинная ось эллипсообразного центра источника была ориентирована под углом 30◦ к направлению юг-север (рис. 3). На рис. 3 показано распределение максимумов высоты цунами, рассчитанное по разностной схеме Мак-Кормака на персональном компьютере с калькулятором на базе FPGA. На рис. 3 видно, что рассчитанные максимальные амплитуды цунами у берега наблюдаются в районах, приблизительно расположенных в направлении короткой оси эллипса. Результаты детальной, промежуточной и грубой сетки визуально не отличаются на изображениях, показывающих всю область, поэтому показан только один рисунок, соответствующий промежуточной сетке (а именно, рис. 3). Разницу высот цунами, рассчитанных на разных сетках по всей области, можно увидеть на рис. 4-7. На рис. 4 и 5 показано распределение разностей максимумов высоты волн по всей области, полученное при сравнении результатов на детальной и промежуточной сетках. Рис. 3. Распределение максимальных высот волн по всей расчетной области в результате численного моделирования распространения цунами от модельного источника эллипсовидной формы на сетке промежуточного разрешения. Fig. 3. Distribution of maximum wave heights over the entire computational domain as a result of numerical simulation of tsunami propagation from a model ellipsoidal source on an intermediate resolution grid. На рис. 4 черным цветом показаны изолинии пространственного распределения разностей между детальной и промежуточной сетками, соответствующие уровню 0,001 м. На рис. 5, демонстрирующем прибрежную часть расчетной области, черной линией обозначены места, где разница максимальных высот цунами превышает 1 см. На рисунке видно, что разница максимальных высот в узлах сеток с одинаковым расположением почти везде не превышает 1 мм, за исключением некоторых прибрежных районов. На рис. 6 и 7 также показано распределение максимальных перепадов высот по результатам моделирования, но в виде сравнения высот на промежуточной и грубой сетках. Сравнительно с рис. 4, 5, на рис. 6, 7 наблюдаются несколько б´ольшие площади акваторий с различиями в результатах более 1 мм (рис. 6) и с различиями более 1 см (рис. 7). Но площадь этих зон существенно меньше площади той части региона, где разница незначительна (менее 1 мм на максимальных высотах у берега порядка 1-2,5 м). Для более детального сравнения результатов на рис. 8 сравниваются максимальные высоты волн в узлах расчетной сетки вблизи участка береговой линии от точки A до точки B, указанного на рис. 4. Анализ распределения максимальных высот на рис. 8 показывает достаточную близость результатов, полученных на детальной и промежуточной сетках. В целом максимальные высоты на детальной сетке незначительно превышают (на 1%-3%) высоты на промежуточной сетке. Аналогичное сравнение было проведено между результатами на промежуточной и грубой сетках. Это отражено на рис. 9. Чтобы количество отсчетов максимальной высоты волн на участке береговой линии от точки A до точки B было одинаковым для промежуточной и грубой сеток, к результатам на грубой сетке добавлены повторы. В результате оба набора данных содержат 1301 значений максимальной высоты. Рис. 4. Распределение разностей между максимальными высотами, рассчитанными на детальной и промежуточной сетках. Изолинии черного цвета соответствуют уровню 1 мм. Fig. 4. Distribution of differences between maximum heights calculated on the detailed and intermediate grids. Black isolines correspond to the 1 mm level. Наконец, на рис. 10 показано сравнение максимальных высот цунами, полученных по детальной сетке (сетка 1) и грубой сетке (сетка 3) вдоль участка береговой линии от точки A до точки B. Распределение максимумов высот цунами (в см) показано на рис. 10. Чтобы уравнять количество отсчетов максимального возвышения над уровнем моря между точками A и B в узлах детальной и грубой сеток, подробные данные были прорежены удалением каждого второго значения, а количество отсчетов в «грубых» данных было удвоено путем повторения каждого значения дважды. 4. Обсуждение Вычислительные эксперименты на выбранной акватории показали, что практически везде увеличение пространственного шага примерно с 110 м до почти 500 м не привело к сколько-нибудь заметному изменению расчетных высот волн цунами. Даже в шельфовой зоне разница в результатах составляет 1-2% от высоты волны. Более заметная разница в амплитудах цунами наблюдается в некоторых точках вблизи береговой линии (см. рис. 8-10). При сравнении результатов численных расчетов на сетках разного разрешения, представленных на рис. 4-10, видно, что в большей части прибрежных вод расчетные максимальные высоты выше для более детальных сеток (рис. 5 и 7). Ту же тенденцию можно наблюдать и по распределениям максимальных высот вдоль береговой линии, представленным на рис. 8-10. В то же время вдоль береговой линии имеется несколько точек, где высоты цунами, рассчитанные по Рис. 5. Распределение разностей между максимальными высотами, рассчитанными на детальной и промежуточной сетках. На этом рисунке более подробно показан сегмент расчетной области, обозначенный белым прямоугольником на рис. 4. Изолинии черного цвета соответствуют уровню 1 см. Fig. 5. Distribution of differences between maximum heights calculated on the detailed and intermediate grids. This figure shows in more detail the segment of the computational domain, indicated by the white rectangle in Fig. 4. Black isolines correspond to the level 1 cm. грубой сетке, оказались выше, чем высоты волн в тех же узлах сетки при использовании более детальной сетки. Это можно объяснить используемыми при численном моделировании граничными условиями у берега, реализующими полное отражение волны в ближайшем к берегу узле расчетной сетки, где глубина составляет не менее 8 м. В связи с тем, что сравнение высот производилось в узлах более грубой из двух сравниваемых сеток, расчетный узел более грубой сетки не всегда является точкой отражения волны при расчете на более детальной сетке. На более детальной сетке может быть узел, находящийся еще ближе к берегу, в котором глубина также не менее 8 м, где моделируется полное отражение, увеличивающее высоту волны вдвое. Из сравнения рис. 8 и 9 видно, что разница высот волн в одних и тех же точках побережья при использовании детальной и промежуточной сеток существенно меньше, чем такая же разница для промежуточной и грубой сеток. Рис. 6. Распределение разностей между максимальными высотами, рассчитанными на промежуточной и грубой сетках. На этом рисунке изолинии черного цвета соответствуют уровню 1 мм. Fig. 6. Distribution of differences between maximum heights calculated on intermediate and coarse grids. In this figure, the black isolines correspond to the 1 mm level. При дальнейшем увеличении шага пространственной сетки до 1-2 км, помимо ухудшения детальности результатов моделирования, будет возрастать численная диссипация разностной схемы, что становится заметным в результатах, если количество узлов сетки на длину волны составляет менее 15-20. В проведенных вычислительных экспериментах волна, генерируемая рассматриваемым модельным источником, изначально имеет длину порядка 60 км на глубине около 3000 м. Затем ее длина уменьшается до 10-15 км вблизи побережья, что приводит к существенному искажению результатов при использовании сеток с пространственным шагом более 1 км. 5. Заключение Численные эксперименты показали, что при 15-20 узлах сетки на длину волны расстояние между точками сетки оказывает незначительное влияние на расчетную максимальную высоту волны. Следовательно, расчетные сетки с более чем 25 узлами на длину волны могут использоваться для повышения точности результирующих значений амплитуды цунами. Этот вывод следует использовать и при применении метода вложенных сеток для расчета параметров волн цунами. В частности, нет необходимости использовать более детальную сетку на достаточно больших глубинах, поскольку это увеличивает время расчета, но не способствует большей точности получаемых значений максимальной высоты волн. Рис. 7. Распределение разностей между максимальными высотами, рассчитанными на промежуточной и грубой сетках. На этом рисунке более подробно показан сегмент расчетной области, обозначенный белым прямоугольником на рис. 6. Черные изолинии соответствуют уровню 1 мм. Fig. 7. Distribution of differences between maximum heights calculated on intermediate and coarse grids. This figure shows in more detail the segment of the computational domain, indicated by the white rectangle in Fig. 6. Black isolines correspond to the 1 mm level. Рис. 8. Распределение максимальных высот в см вдоль побережья острова Хонсю между точками A и B (см. рис. 4), рассчитанное на детальной сетке (серый цвет) и на промежуточной сетке (черный цвет). По горизонтальной оси указано расстояние по широте (с шагом 0,002◦) от точки A. Fig. 8. Distribution of maximum heights in cm along the coast of Honshu Island between points A and B (see Fig. 4), calculated on the detailed grid (gray) and on the intermediate grid (black). The horizontal axis indicates the latitude distance (with 0.002◦ steps) from the point A. Рис. 9. Распределение максимальных высот вдоль участка побережья острова Хонсю от точки A до точки B (см. рис. 4), рассчитанное на промежуточной сетке (черный цвет) и на грубой сетке (серый цвет). Fig. 9. Distribution of maximum heights along the coastline of Honshu Island from point A to point B (see Fig. 4), calculated on an intermediate grid (black) and on a coarse grid (gray). Рис. 10. Распределение максимальных высот вдоль участка побережья острова Хонсю от точки A до точки B, рассчитанное на детальной сетке (черный цвет) и на грубой сетке (серый цвет). Fig. 10. Distributions of maximum heights along the coastline of Honshu Island from point A to point B, calculated on a fine grid (black) and on a coarse grid (gray).

Об авторах

М. М. Лаврентьев

Институт автоматики и электрометрии СО РАН

Автор, ответственный за переписку.
Email: mmlavrentiev@gmail.com
Новосибирск, Россия

К. Ф. Лысаков

Институт автоматики и электрометрии СО РАН

Email: lysakov@sl.iae.nsk.su
Новосибирск, Россия

Ан. Г. Марчук

Институт вычислительной математики и математической геофизики СО РАН

Email: mag@omzg.sscc.ru
Новосибирск, Россия

К. К. Облаухов

Институт автоматики и электрометрии СО РАН

Email: oblaukhov.konstantin@gmail.com
Новосибирск, Россия

М. Ю. Шадрин

Институт автоматики и электрометрии СО РАН

Email: mikesha@sl.iae.nsk.su
Новосибирск, Россия

Список литературы