В данной статье будет рассказано о симуляции распространения на большие расстояния волн цунами, порождённых подводными землетрясениями, при помощи программного пакета Wolfram Mathematica. В качестве математической модели системы используется вращающийся гравитирующий геоид с нетривиальной поверхностью дна и выколотыми областями, которые эмулируют материковые образования, и вязкая несжимаемая жидкость на его поверхности.
Введение
По определению, цунами - это сейсмические поверхностные волны с очень большим отношением длинны волны к глубине водоёма 
(для земных океанов длина волны составляет порядка сотни километров), вызванные резким подъемом или опусканием небольшого по площади столба жидкости одновременно во всей толще. Другими словами, волны цунами классифицируются как вторичные явления после крупных геофизических катаклизмов, таких как подводные землетрясения, оползни, сходы берегов, извержения вулканов и т.п., и образуется из-за резкого изменения профиля донной поверхности.
Как следствие своей физической природы, волны цунами имеют несколько отличительных особенностей.
Во-первых, скорость волн цунами даётся выражением
, где 
- ускорение свободного падения, а 
- глубина. Отсюда сразу видно, что волны распространяются быстрее в открытом океане, нежели чем у побережья.Во-вторых, всё также на основании вышеприведённой формулы, цунами всегда автоматически и вне зависимости от начального угла разворачиваются по направлению к берегу анфас, так как отстающий участок фронта на глубине догоняет опережающий, замедлившийся на мелководье.
В-третьих, энергия волн цунами сильно "размазана" по толще жидкости, поэтому на больших глубинах цунами имеют высоту лишь несколько десятков сантиметров, но при выходе на шельф, в силу закона сохранения энергии, высота гребня может увеличиваться до десятков и даже сотен метров.
Эти отличительные черты делают цунами крайне опасными для человека и существенно усложняют своевременное реагирования на подобного сорта события.
Широко известно, что среди всех природных катаклизмов, одно из лидирующих мест по числу ��еловеческих жертв и урону инфраструктуре занимают именно волны цунами. Одно из самый разрушительных цунами за всю историю наблюдений произошло в Суматре в 2004 году. Общее число единовременных жертв в Индонезийском бедствии превысило 230 тысяч человек. Неприятной особенностью того страшного события была близость источника цунами - подводного землетрясения - к густонаселенному побережью.
К сожалению, геофизики до сих пор не обладают универсальной методикой прогнозирования цунами. Во многих странах мира, которые подвержены риску воздействия волновых явлений, существуют системы экстренного оповещения и эвакуации населения. Однако очень часто подобные системы ошибаются в предсказании параметров и локализации приходящих волн, или же наоборот - совершают переоценку рисков, что приводит к ложным тревогам и, соответственно, экономическим убыткам.
В России опасность цунами имеется в двух регионах: на Камчатке и на Курильских островах (Сахалинская область). Основным источником "наших" цунами является Тихоокеанское вулканическое огненное кольцо - сейсмически активная область по периметру Тихого океана, на которой расположено большое количество действующих вулканов. Из-за особенностей рельефа океанического дна, а также в силу географического положения, чаще всего к нам приходят цунами из Чилийской зоны в Южной Америке.
Таким образом, теоретическая задача исследования и практическая проблема предсказания цунами остро стоят на повестке дня как в общем смысле мирового знания, так и в контексте повышен��я уровня жизни в Дальневосточном регионе России.
Физическая модель цунами
Произвольную динамику классической жидкости всеобъемлюще можно описать при помощи гидродинамических дифференциальных уравнений Навье-Стокса в частных производных. В данной работе в качестве модели рассматривается вязкая несжимаемая жидкость, которая находится на монотонно вращающемся притягивающем массу геоиде с нетривиальной формой дна и выколотыми областями под внешним атмосферным давлением. Для описания распространения волн используется так называемое приближение мелкой воды. Данное упрощение вполне оправдано, так как волны цунами по определению имеют очень большие длины.
Описанная выше модель хорошо применима для описания цунами на больших и средних глубинах. В самой близости от побережья, когда гребень волны начинается опрокидываться, характер движения жидкости меняется с ламинарного на турбулентный, в связи с чем закон эволюция волн во времени качественно изменяется. Однако исходная постановка задачи предполагает лишь определение высоты пришедшей волны вблизи берега в зависимости от координаты и времени, в то время как процесс обрушения волнового гребня является существенно более нетривиальной проблемой и лежит вне рамок текущего обсуждения.
Полноценный вывод окончательных симулируемых уравнений является довольно сложной и кропотливой задачей, лежащей за пределами данной статьи. Тем не менее, для более всестороннего пони��ания сути требуемых преобразований необходимо привести их краткое последовательное описание.
К исходной системе уравнений Навье-Стокса в правую часть важно добавить центробежную силу и силу Кориолиса, действующие на единичный объём жидкости во вращающейся системе отсчёта.
Для учёта кинематической вязкости в правую часть исходной системы дифференциальных уравнений следует дописать слагаемое с лапласианом.
Далее нужно аккуратно перейти из декартовой системы координат в сферическую, попутно переписав все дифференциальные операторы при помощи известной матрицы Якоби.
На последнем этапе необходимо выполнить интегрирование скорости вдоль медленной радиальной компоненты, параллельно отбрасывая слагаемые с порядком малости
, где 
- средний радиус геоида. Это и будет реализацией приближения мелкой воды.
После всех вышеупомянутых действий уравнения Навье-Стокса принимают финальный вид:
![\begin{equation*} \begin{aligned} & r_{0}\partial_{t}\left[h(t, \theta, \varphi) - f(t, \theta, \varphi)\right] + v(t, \theta, \varphi)\partial_{\theta}\left[h(t, \theta, \varphi) - f(t, \theta, \varphi)\right] + \frac{u(t, \theta, \varphi)}{\sin(\theta)}\partial_{\varphi}\left[h(t, \theta, \varphi) - f(t, \theta, \varphi)\right] = \\ & = -\frac{\left[h(t, \theta, \varphi) - f(t, \theta, \varphi) + d\right]}{\sin(\theta)}\left\{\partial_{\theta}\left[v(t, \theta, \varphi)\sin(\theta)\right] + \partial_{\varphi}u(t, \theta, \varphi)\right\}, \\ & r_{0}\partial_{t}v(t, \theta, \varphi) + v(t, \theta, \varphi)\partial_{\theta}v(t, \theta, \varphi) + \frac{u(t, \theta, \varphi)}{\sin(\theta)}\partial_{\varphi}v(t, \theta, \varphi) = u^{2}(t, \theta, \varphi)\frac{\cos(\theta)}{\sin(\theta)} + \\ & + 2\left(\Omega_{0}r_{0}\right)u(t, \theta, \varphi)\cos(\theta) + \left(\Omega_{0}r_{0}\right)^{2}\sin(\theta)\cos(\theta) - g\partial_{\theta}h(t, \theta, \varphi) - \frac{1}{\rho}\partial_{\theta}P(t, \theta, \varphi) + \\ & + \frac{\nu}{r_{0}}\left\{\frac{\partial_{\theta}\left[\sin(\theta)\partial_{\theta}v(t, \theta, \varphi)\right]}{\sin(\theta)} + \frac{\partial_{\varphi}^{2}v(t, \theta, \varphi)}{\sin(\theta)^{2}}\right\}, \\ & r_{0}\partial_{t}u(t, \theta, \varphi) + v(t, \theta, \varphi)\partial_{\theta}u(t, \theta, \varphi) + \frac{u(t, \theta, \varphi)}{\sin(\theta)}\partial_{\varphi}u(t, \theta, \varphi) = -v(t, \theta, \varphi)u(t, \theta, \varphi)\frac{\cos(\theta)}{\sin(\theta)} - \\ & - 2\left(\Omega_{0}r_{0}\right)v(t, \theta, \varphi)\cos(\theta) - g\frac{\partial_{\varphi}h(t, \theta, \varphi)}{\sin(\theta)} - \frac{1}{\rho}\frac{\partial_{\varphi}P(t, \theta, \varphi)}{\sin(\theta)} + \\ & + \frac{\nu}{r_{0}}\left\{\frac{\partial_{\theta}\left[\sin(\theta)\partial_{\theta}u(t, \theta, \varphi)\right]}{\sin(\theta)} + \frac{\partial_{\varphi}^{2}u(t, \theta, \varphi)}{\sin(\theta)^{2}}\right\}. \end{aligned} \end{equation*}](https://habrastorage.org/getpro/habr/formulas/7/77/77d/77d9ed83e23b8d7f25cbb3e33648e3d3.svg)
Обозначения: - средний радиус геоида, 
- угловая скорость вращения, 
- плотность жидкости, 
- кинематическая вязкость, - ускорение свободного падения на поверхности геоида, 
- время, 
- широта, 
- долгота, 
- общий равновесный уровень жидкости от нулевого значения, 
- измение высоты уровня жидкости над равновесным, 
- функция рельефа дна, 
- компонента скорости вдоль широты, 
- компонента скорости вдоль долготы, 
- атмосферное давление.
Заинтересованный читатель может найти вывод указанных уравнений с разной степенью детализации в работах:
Л.Д. Ландау, Е.М. Лифшиц, Теоретическая физика, том 6.
М.А. Носов, Введение в теорию волн цунами. ISBN 978-5-8037-0773-8.
А.А. Черевко, А.П. Чупахин, Уравнения Модели Мелкой Воды на Вращающейся Притягивающей Сфере.
З.И. Федотова, Г.С. Хакимзянов, Иерархия уравнений мелкой воды: вывод, исследование, вычислительные алгоритмы. Прикладная механика и техническая физика. 2009. Т. 50, № 2.
R.O. Popovych, Physics-informed neural networks for the shallow-water equations on the sphere, DOI:10.48550/arXiv.2104.00615.
Программная подготовка модели
Моделирование волн производится на геоиде с близкими к земным параметрами:
r0 = 6.4*10^6; (*Earth radius*) g = 9.81; (*Earth gravitational constant*) \[Omega]0 = 7.3*10^-5; (*Earth angular speed*) d = 3.8*10^3; (*average depth of ocean between surface and floor*) p0 = 1.013*10^5; (*normal atmospheric pressure*) \[Rho] = 10^3; (*water normal density*) \[Nu] = 1.003*10^-6; (*water normal kinematic viscosity*)
Как известно, сферические координаты имеют особенности около полюсов. С точки зрения физики, поведение жидкости при больших широтах ничем качественно не выделяется, но для корректного математического описания при 
необходимо правильно сшивать карты. Во избежании затруднений, связанных с этим фактом, и для минимизации численных ошибок, шапки геоида будут интерпретироваться как некоторые материковые образования. Таким образом, карта долгот и широт имеет следующие границы:
\[Theta]min = 0 + \[Pi]/20; \[Theta]max = \[Pi] - \[Pi]/20; \[Phi]min = 0; \[Phi]max = 2 \[Pi];
В качестве невозмущённого профиля дна будет рассмотрена периодическая функция:
FloorProfile[\[Theta]_, \[Phi]_] := d/4*Sin[2 \[Theta]]*(1 + 1/2 Sin[\[Phi]]);
Из-за наличия вращения средний равновесный уровень жидкости должен иметь форму:
Compensator[\[Theta]_, \[Phi]_] := 1/(2*g)*(\[Omega]0*r0)^2*Sin[\[Theta]]^2;
Исследуемая далее волна цунами будет порождаться землетрясением, модельная суть которого состоит в очень быстром подъёме небольшого участка донной поверхности с координатами:
\[Theta]exmin = 0.45*\[Pi] - \[Pi]/4; \[Theta]exmax = 0.55*\[Pi] - \[Pi]/4; \[Phi]exmin = 0.95*\[Pi] + \[Pi]/2; \[Phi]exmax = 1.05*\[Pi] + \[Pi]/2;
Также для землетрясения нужно задать амплитуду и продолжительность:
aex = d/100; tex = 10;
Сам процесс подводного толчка для примера описывается функцией:
FloorProfileExcitation[t_, \[Theta]_, \[Phi]_] := If[t <= 0, 0, If[t >= tex, aex, aex/tex*t]]* If[\[Theta]exmin < \[Theta] < \[Theta]exmax, Sin[\[Pi]*((\[Theta] - \[Theta]exmin)/(\[Theta]exmax - \ \[Theta]exmin))]^2, 0]* If[\[Phi]exmin < \[Phi] < \[Phi]exmax, Sin[\[Pi]*((\[Phi] - \[Phi]exmin)/(\[Phi]exmax - \ \[Phi]exmin))]^2, 0];
Для большей реалистичности на поверхность геоида можно нанести несколько полигональных фигур, имитирующих материки:
region = Rectangle[{\[Theta]min, \[Phi]min}, {\[Theta]max, \[Phi]max}]; region = RegionDifference[region, Polygon[{{\[Theta]min + \[Pi]/3, (2 \[Pi])/ 3}, {(\[Theta]max + \[Theta]min)/2, (5 \[Pi])/ 6}, {\[Theta]max - \[Pi]/3, (2 \[Pi])/3}}]]; region = RegionDifference[region, RegularPolygon[{(\[Theta]max + \[Theta]min)/ 2, \[Phi]min + \[Pi]/4}, {\[Pi]/10, 0}, 6]]; region = RegionDifference[region, Disk[{(\[Theta]max + \[Theta]min)/2, \[Phi]max - \[Pi]/2}, \[Pi]/ 10]]; region = RegionDifference[ RegionDifference[region, Rectangle[{\[Theta]min, 0.95*\[Pi]}, {\[Theta]min + \[Pi]/4, 1.05*\[Pi]}]], Rectangle[{\[Theta]max - \[Pi]/4, 0.95*\[Pi]}, {\[Theta]max, 1.05*\[Pi]}]];
Для быстрого численного решения уравнений Навье-Стокса при помощи встроенных утилит языка Wolfram Language необходимо произвести триангуляцию домена средствами специального расширения FEM:
<< NDSolve`FEM` meshregion = ToElementMesh[region, "MeshElementType" -> TriangleElement, MeshQualityGoal -> "Maximal", MaxCellMeasure -> {"Length" -> 0.4, "Area" -> 5}, "MaxBoundaryCellMeasure" -> 0.2, "ImproveBoundaryPosition" -> True, "BoundaryMeshGenerator" -> "Continuation"]; meshboundary = ToBoundaryMesh[meshregion, MaxCellMeasure -> {"Length" -> 0.4}, "MaxBoundaryCellMeasure" -> 0.2];
Уравнения и краевая задача
Уравнения Навье-Стокса в терминах синтаксиса языка Wolfram имеют вид:
equations = {r0*\!\( \*SubscriptBox[\(\[PartialD]\), \(t\)]\((h[t, \ \[Theta], \ \[Phi]] - f[t, \ \[Theta], \ \[Phi]])\)\) + v[t, \[Theta], \[Phi]]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((h[ t, \ \[Theta], \ \[Phi]] - f[t, \ \[Theta], \ \[Phi]])\)\) + u[t, \[Theta], \[Phi]]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\((h[ t, \ \[Theta], \ \[Phi]] - f[t, \ \[Theta], \ \[Phi]])\)\)/ Sin[\[Theta]] + (h[t, \[Theta], \[Phi]] - f[t, \[Theta], \[Phi]] + d)*(\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((v[ t, \ \[Theta], \ \[Phi]]*Sin[\[Theta]])\)\) + \!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(u[ t, \ \[Theta], \ \[Phi]]\)\))/Sin[\[Theta]] == 0, r0*\!\( \*SubscriptBox[\(\[PartialD]\), \(t\)]\(v[ t, \ \[Theta], \ \[Phi]]\)\) + v[t, \[Theta], \[Phi]]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(v[ t, \ \[Theta], \ \[Phi]]\)\) + u[t, \[Theta], \[Phi]]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(v[ t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] == u[t, \[Theta], \[Phi]]^2*Cos[\[Theta]]/Sin[\[Theta]] + 2*\[Omega]0*r0*u[t, \[Theta], \[Phi]]*Cos[\[Theta]] - g*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(h[ t, \ \[Theta], \ \[Phi]]\)\) - 1/\[Rho]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(p[ t, \ \[Theta], \ \[Phi]]\)\) + \[Nu]/r0*( \!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((Sin[\[Theta]]*\ \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]v[ t, \ \[Theta], \ \[Phi]])\)\)/Sin[\[Theta]] + \!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]v[ t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]]^2), r0*\!\( \*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[ t, \ \[Theta], \ \[Phi]]\)\) + v[t, \[Theta], \[Phi]]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(u[ t, \ \[Theta], \ \[Phi]]\)\) + u[t, \[Theta], \[Phi]]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(u[ t, \ \[Theta], \ \[Phi]]\)\)/ Sin[\[Theta]] == -v[t, \[Theta], \[Phi]]*u[t, \[Theta], \[Phi]]* Cos[\[Theta]]/Sin[\[Theta]] - 2*\[Omega]0*r0*v[t, \[Theta], \[Phi]]*Cos[\[Theta]] - g*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(h[ t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] - 1/\[Rho]*\!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(p[ t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] + \[Nu]/r0*( \!\( \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((Sin[\[Theta]]*\ \*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]u[ t, \ \[Theta], \ \[Phi]])\)\)/Sin[\[Theta]] + \!\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\( \*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]u[ t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]]^2)};
Начальные значения скоростей и подъема жидкости, а также внешние фиксированные функции поверхности дна и атмосферного давления имеют вид:
v0[\[Theta]_, \[Phi]_] := 0; u0[\[Theta]_, \[Phi]_] := 0; h0[\[Theta]_, \[Phi]_] := 0; f[t_, \[Theta]_, \[Phi]_] := FloorProfile[\[Theta], \[Phi]] + FloorProfileExcitation[t, \[Theta], \[Phi]]; p[t_, \[Theta]_, \[Phi]_] := p0;
Вышеуказанные выражения позволяют сразу же сформулировать начальные условия для системы дифференциальных уравнений:
initialconditions = {h[0, \[Theta], \[Phi]] == h0[\[Theta], \[Phi]], v[0, \[Theta], \[Phi]] == v0[\[Theta], \[Phi]], u[0, \[Theta], \[Phi]] == u0[\[Theta], \[Phi]]};
В силу того, что сфера является периодической структурой, необходимо потребовать граничные условия на циклическую сшивку долготы:
periodicboundaryconditions = {PeriodicBoundaryCondition[ h[t, \[Theta], \[Phi]], \[Phi] == \[Phi]max, TranslationTransform[{0, -\[Phi]max}]], PeriodicBoundaryCondition[ v[t, \[Theta], \[Phi]], \[Phi] == \[Phi]max, TranslationTransform[{0, -\[Phi]max}]], PeriodicBoundaryCondition[ u[t, \[Theta], \[Phi]], \[Phi] == \[Phi]max, TranslationTransform[{0, -\[Phi]max}]], PeriodicBoundaryCondition[ h[t, \[Theta], \[Phi]], \[Phi] == \[Phi]min, TranslationTransform[{0, \[Phi]max}]], PeriodicBoundaryCondition[ v[t, \[Theta], \[Phi]], \[Phi] == \[Phi]min, TranslationTransform[{0, \[Phi]max}]], PeriodicBoundaryCondition[ u[t, \[Theta], \[Phi]], \[Phi] == \[Phi]min, TranslationTransform[{0, \[Phi]max}]]};
И, наконец, для предотвращения расходимостей на границах области симуляции очень важно поставить граничную задачу Дирихле:
dirichletboundaryconditions = DirichletCondition[{h[t, \[Theta], \[Phi]] == h0[\[Theta], \[Phi]], v[t, \[Theta], \[Phi]] == v0[\[Theta], \[Phi]], u[t, \[Theta], \[Phi]] == u0[\[Theta], \[Phi]]}, {\[Theta], \[Phi]} \[Element] meshboundary && \[Phi]min < \[Phi] < \[Phi]max];
Физически такие условия в некотором смысле моделируют высокие стенки бассейна, от которых набегающие волны абсолютно упруго отражаются.
Численное решение СДУ
Уравнения Навье-Стокса на сфере сами по себе являются сложной системой для симуляции. Поэтому, несмотря на всю мощь и автоматизированность языковой среды Wolfram, метод численных вычислений для такого рода непростых задач необходимо настраивать вручную, эмпирически подбирая удачные параметры:
method = {"PDEDiscretization" -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"FiniteElement", "BoundaryTolerance" -> Automatic, "IntegrationOrder" -> Automatic, "InterpolationOrder" -> Automatic, "MeshOptions" -> {"MeshElementType" -> "TriangleElement"}}}};
Wolfram Mathematica из коробки умеет работать с процессорной многопоточностью, однако использование связки CUDA/GPU всё ещё недоступно для целого ряда встроенных методов. Поэтому длительность и параметры симуляции выбирались во многом исходя из разумного времени вычисления. Во время экспериментов использовался довольно взрослый игровой ноутбук Asus N751J (8 Gb RAM, Intel Core i7 4710HQ 2,5 ГГц, 4 physical cores). Для конкретно данной машины примерно за 
минут работы в реальном времени на сетке из ~
элементов удаётся просчитать 
секунд симуляции.
time = 2*10^5;
Этого временного интервала вполне достаточно для практических задач, так как волновой фронт, движущийся со скоростью 
метров в секунду, за указанное время проходит расстояние 
тысяч километров, что соответствует среднему диаметру земного геоида.
Таким образом, численное решение СДУ Навье-Стокса может быть проделано следующим программным модулем:
Module[{timer = 0}, Monitor[solution = NDSolveValue[{equations, initialconditions, dirichletboundaryconditions, periodicboundaryconditions}, {h, v, u}, {t, 0, time}, {\[Theta], \[Phi]} \[Element] meshregion, Method -> method, InitialSeeding -> initialconditions, MaxSteps -> 20000, MaxStepSize -> 100, StartingStepSize -> 5, InterpolationOrder -> All, WorkingPrecision -> MachinePrecision, AccuracyGoal -> Automatic, PrecisionGoal -> Automatic, StepMonitor :> (timer = t)], ProgressIndicator[timer/time]]]; // AbsoluteTiming
Результаты симуляции
На графиках отчётливо видно, что, в целом, симуляция верно передаёт обобщённые характеристики процесса распространения цунами: присутствует и впоследствии разъезжается водяной колокол, четко прослеживается первичный волновой фронт, имеется правильно учтённая периодичность сферы, есть эффект замедления волны на шельфовой части дна, а также корректно реализуется взаимодействие волн с препятствиями. Однако, что особенно хорошо видно на последних этапах симуляции, поверхность жидкости после прохода основного вала обрабатывается весьма посредственно: вместо постепенной релаксации она начинает колебаться с постоянно возрастающей амплитудой. Это паразитическое явление, скорее всего, есть следствие накопления численных ошибок. Другими словами, проведённая симуляция позволяет правильно уловить общие штрихи пропагирования волны, но конкретная детализация динамики поверхности не является релевантной.
Скрытый текст
Nota Bene: Если отказаться от условия периодичности для карты геоида, то паразитные колебания поверхности жидкости заметно уменьшаются, хоть и не пропадает полностью.
Количественные оценки качества симуляции
Для подтверждения вышеупомянутой гипотезы о накоплении численных ошибок, имеет смысл построить график зависимости полной энергии системы от времени:
Monitor[energydispersion = Table[{t, Quiet[NIntegrate[\[Rho]*(((solution[[2]][t, \[Theta], \[Phi]] - u0[\[Theta], \[Phi]])^2 + (solution[[3]][ t, \[Theta], \[Phi]] - v0[\[Theta], \[Phi]])^2) + g/2*Abs[solution[[1]][t, \[Theta], \[Phi]] - h0[\[Theta], \[Phi]]])* Abs[solution[[1]][t, \[Theta], \[Phi]] - h0[\[Theta], \[Phi]]]*r0^2* Sin[\[Theta]], {\[Theta], \[Phi]} \[Element] meshregion, Method -> {Automatic, "SymbolicProcessing" -> False}, Exclusions -> All, WorkingPrecision -> MachinePrecision, AccuracyGoal -> Automatic, PrecisionGoal -> Automatic]]}, {t, 0, time, time/50}], ProgressIndicator[t/time]]; // AbsoluteTiming
Согласно численной оценке, энергия системы за всё время симуляции выросла почти на порядок. Очевидно, что из-за этого эффекта степень достоверности результатов существенно ограничена на больших временах.
Скрытый текст
В случае апериодического домена рост энергии на момент окончания симуляции значительно меньше и составляет около 200%.
Ещё одним способом проверки валидности симуляции динамики жидкости является расчёт интегрального вектор-потока через стенки резервуара. Для непроницаемых стенок этот поток, очевидно, должен зануляться:
Monitor[Total[ Table[time/50* Quiet[NIntegrate[{solution[[2]][t, \[Theta], \[Phi]], solution[[3]][t, \[Theta], \[Phi]]} . {{1, 0}, {0, 1/ Sin[\[Theta]]}} . BoundaryUnitNormal[\[Theta], \[Phi]]* Sin[\[Theta]], {\[Theta], \[Phi]} \[Element] meshboundary, Method -> {Automatic, "SymbolicProcessing" -> False}, Exclusions -> All, WorkingPrecision -> MachinePrecision, AccuracyGoal -> Automatic, PrecisionGoal -> Automatic]], {t, time/50, time, time/50}]], ProgressIndicator[t/time]] // AbsoluteTiming
Оказывается, что для текущих вычислений значение полного потока, действительно, близко к нулю. Отсюда можно сделать вывод, что алгоритм решения СДУ в данной задаче чаще ошибается не при интегрировании пространственных размерностей, а при переходе от одной временной линии к другой.
Комментарий от автора
Численное решение нелинейных уравнений Навье-Стокса во вращающихся криволинейных координатах на домене сложной формы - крайне нетривиальная задача. Очень редко удаётся быстро и безболезненно получить красивый результат. В подавляющем большинстве случаев, наоборот, приходится долго и упорно искать удачный приём или алгоритм среди методов вычислительной физики.
Wolfram Mathematica - очень мощный и "умный" программный комплекс, который по умолчанию "из коробки" имеет широкий функционал для различных символьных и численных проектов. Встроенные функции Wolfram часто автоматически адаптируются под поставленную пользователем задачу, самостоятельно подбирая удачные внутренние параметры и алгоритмы. В этом, с одной стороны, крайне большой плюс и одновременно значительный минус данного языка, особенно в контескте решения многоэтапных и ресурсоёмких задач. Иногда возникают такие ситуации, когда приходится углубляться в недокументированные либо вспомогательные функции среды, чтобы принудительно переключить тот или иной внутренний тумблер. Более того, несмотря на достаточно серьёзную оптимизацию функционала, многие методы не имеют удобного многопоточного или же CUDA-friendly аналога, что препятствует эффективному экстенсивному масштабированию тяжелых проектов на параллельных исполнителях. Основываясь на приведённых фактах, напрашивается логический вывод: Wolfram Language является удобным и комфортным средством для относительно быстрой проверки гипотез и получения "MVP" в короткие сроки. Однако на дальней дистанции имеет смысл отдавать предпочтение GPU-совместимым языкам или библиотекам для Python, таким как Jax, Numba, PyTorch etc.
Отдельным пунктом хочется отметить, что Wolfram Language - это не самый популярный язык в сообществе физиков-программистов, в связи с чем попытки оптимизировать скрипты кода при помощи современных LLM-мок часто приводят к нерабочим решениям как прямое следствие галлюцинаций.
Дальнейшие изыскания по тематике текущей статьи могут быть потенциально направлены в сторону борьбы с патологическим ростом энергии, увеличения числа элементов на триангулированной сетке, а также улучшения детализации границ домена.
Полезные ссылки
Репозиторий с кодом на GitHub.
Дискуссия на StackExchange.
