Грант РФФИ 18-47-340007 р_а

РФФИ 18-47-340007 р_а «Суперкомпьютерное гидродинамическое моделирование русловых процессов, размыва и последствий прорыва заградительных дамб на водных объектах Волгоградской области»

Отчет за 2018-2019 года

  Руководитель проекта
 

Храпов Сергей Сергеевич

  Код и название Конкурса
 

р_поволжье_а Региональный конкурс Поволжье: инициативные

  Аннотация
 

Актуальность проекта по развитию методов численного гидродинамического моделирования и разработке эффективных суперкомпьютерных алгоритмов обусловлена необходимостью решения фундаментальных и практических задач связанных с исследованием русловых процессов и катастрофического затопления территорий в период сезонного и/или ливневого паводка. Разработанные в рамках проекта математические и численные модели самосогласованной динамики поверхностных вод и донных наносов основаны на оригинальном лагранжево-эйлеровом методе CSPH-TVD. Создан программный комплекс, позволяющий проводить быстрое имитационное моделирование русловых процессов и динамики затопления на обширных территориях с высоким пространственным разрешением, основанный на компьютерном алгоритме для гибридных вычислительных систем с графическими процессорами Tesla K20/40/80 и P100. Этот алгоритм разработан с использованием технологий параллельных вычислений OpenMP-CUDA и GPUDirect. Для потенциально опасных и наиболее значимых водных объектов Волгоградской области, построены цифровые модели местности и проведено имитационное гидродинамическое моделирование процессов затопления населенных пунктов и объектов инфраструктуры в период сезонного и/или ливневого паводка на реках и водохранилищах региона (р. Волга, р. Ахтуба, р. Дон, р. Хопер, р. Иловля, р. Медведица, Волгоградское водохранилище и др.).  На основе результатов моделирования определены границы зон затопления для различной интенсивности наводнений.

Развернутый научный отчет

  Полученные в 2018-2019 годах важнейшие результаты
 

1) Построена самосогласованная математическая модель динамики поверхностных вод и влекомых наносов, основанных на уравнениях Сен-Венана и Экснера, с учетом основных физических факторов, влияющих на процесс перемещения воды и поверхностного слоя грунта/почвы: рельефа местности, вращения Земли (сила Кориолиса), испарения и ветрового режима, источников воды (плотины и осадки с заданной параметрами притока), инфильтрации воды в почву (в рамках многослойных моделей впитывания), придонного трения, турбулентной вязкости потока, придонного потока наносов (в рамках моделей Грасса и Шилдса), зависящего от величины придонного трения и уклона дна. В комплексной математической модели учтены законы сохранения массы и импульса взаимодействующих компонент. Построенная математическая модель позволяет корректно описывать динамику рассматриваемых процессов.
2) Разработан устойчивый численный алгоритм решения уравнений самосогласованной динамики поверхностных вод и влекомых наносов, основанной на оригинальном лагранжево-эйлеровой методе CSPH-TVD, и определен критерий устойчивости численного алгоритма в зависимости от параметров математической модели и, в первую очередь, от параметров моделей придонного потока наносов, позволяющий проводить устойчивые гидродинамические расчеты в широком диапазоне параметров течения. 
3) Создан программный комплекс, позволяющий проводить быстрое имитационное моделирование русловых процессов и динамики затопления на обширных территориях (100 - 10000 кв. км) с высоким пространственным разрешением (1-10 м), основанный на компьютерном алгоритме для гибридных вычислительных систем с сопроцессорами Nvidia Tesla K20/40/80, разработанном с использованием технологий параллельных вычислений MPI-OpenMP-CUDA и GPUDirect. Проведено исследование эффективности распараллеливания алгоритма для различных конфигураций вычислительных систем CPU+GPU и определены оптимальные параметры этих конфигураций для одиночных и массовых гидродинамических расчетов.
4) Проведен анализ водных объектов (реки, каналы, озера, водохранилища и т.п.) Волгоградской области и определены потенциально опасные и наиболее значимые для гидрологического режима территорий объекты. Для потенциально опасных и наиболее значимых водных объектов построены актуализированные цифровые модели местности с пространственным разрешением (1-10 м) на основе картографических данных, данных дистанционного зондирования Земли (ДДЗЗ) и топографической съемки с использованием беспилотных летательных аппаратов (БПЛА). Цифровая модель местности (ЦММ) включает: цифровую модель рельефа (ЦМР) и матрицы качества - пространственное распределение коэффициентов шероховатости и инфильтрации, характеристик грунта (пористость, медианная крупность и плотность частиц). 
5) Для потенциально опасных и наиболее значимых водных объектов Волгоградской области, выявленных в результате проведенного анализа, проведено имитационное гидродинамическое моделирование русловых процессов, размыва и последствий прорыва заградительных дамб (плотин) в период сезонного и/или ливневого паводка. На основе результатов моделирования построены карты затопления.

  Сопоставление полученных результатов с мировым уровнем
 

Русла рек имеют сложную нерегулярную структуру дна и береговой линии. Их морфометрическими особенностями являются крупномасштабные и мелкомасштабные донные неоднородности, меандрирование [Гидравлические сопротивления и скоростные поля потоков в руслах сложных форм сечения / Н. Б. Барышников, М. С. Дрегваль, Д. И. Исаев, И. С. Гаврилов // Ученые записки Российского государственного гидрометеорологического университета, 2017, 46, 10–20; Pradhan A., Kumar Khatua K. // Sankalp Advances in Civil Engineering, 2018, 1569271; Xu D,. Ji C., Bai Y., Song X., // Journal of Hydroinformatics, 2017, 19 (5), 666-685], вариации сечения канала [Neal J C, Odoni N A, Trigg M A, Freer J E, Garcia-Pintado J, Mason D.C., Bates P.D. // Journal of Hydrology, 2015, 529, 169-183-5; Bahram R., Knight D.W. // Journal of Hydraulic Engineering, 2010, 137 (8), 815-824; Naik B., Khatua K.K. // Aquatic Procedia, 4, 1500-1507] и гидравлического радиуса. Каждый из этих факторов вносит свой вклад в сопротивление потоку канала [Mohammadi S., Kashefipour S.M. // Water resources, 2014, 41 (4), 412-420]. Гидрологическая модель Шези традиционно используется для оценки гидродинамического сопротивления через коэффициент шероховатости Маннинга [Коэффициенты шероховатости пойм / Н. Б. Барышников, Е. С. Субботина, Е. М. Скоморохова, Е. А. Поташко // Ученые записки Российского государственного гидрометеорологического университета, 2012, 23, 13–20; Ghani, Aminuddin A.B., et al. // International Journal of River Basin Management, 2007, 5 (4), 329-346], [2]. Существует несколько различных подходов к определению коэффициента. Поскольку коэффициент шероховатости является феноменологическим параметром, учитывающим большое количество разнородных факторов, возникают определенные трудности при его оценке [2]. Практически при расчете пропускной способности каналов коэффициент шероховатости можно определить из специальных таблиц в соответствии с характеристиками каналов [Hatzigiannakis E., Pantelakis D., Hatzispiroglou I., Arampatzis G., Ilias A., Panagopoulos A. // Environmental Processes, 2016, 3 (1), 263-275]. Для расчета коэффициента Маннинга существует большое количество формул, справедливых для частных случаев, их вывод основан на натурных наблюдениях и лабораторных экспериментах [Habibi, M. An experimental investigation to calculate flow resistance in a steep river / M. Habibi, M. R. Namaee, M. Saneie // KSCE Journal of Civil Engineering, 2014, 18, 4, 1176–1184; Lee, J. S. Composite flow resistance / J. S. Lee, P. Y. Julien // Journal of flood engineering, 2017, 8 (2), 55–75]. В данном проекте при построении моделей самосогласованной динамики поверхностных вод и влекомых наносов расчет коэффициента шероховатости по Маннинга будет осуществляться на основе сравнения результатов численного гидродинамического моделирования в меандрированном и прямом руслах с гладким и содержащим мелкомасштабные неоднородности дном [2, 3].
Выделим ряд задач, связанных с получением теоретической зависимости коэффициентов Манинга и Шези (коэффициенты шероховатости в уравнениях Сен-Венана) от скорости течения, уровня турбулизации потока, глубины, амплитуды и характерных шкал неоднородностей дна, типа дна речного русла или затопляемой поверхности. Метод построения основывается на сравнении расчетов динамики воды с использованием вязких уравнений Сен-Венана и трехмерных моделей с учетом мелкомасштабной неоднородности и пространственной нерегулярности профиля дна. Учет микротурбулентности с характерными пространственными масштабами меньшими размера расчетных ячеек основан на использовании подсеточных моделей турбулентности. Применение подсеточных методов для моделирования различных физических факторов (турбулентности, вязкости, переноса тепла и излучения) основано на определенных правилах усреднения. При использовании классического метода крупных вихрей применяются различные правила подсеточного замыкания [Davidson L. Lecture notes on large eddy simulation. Goeteborg, 2000; Meinke et al., Computer&Fluids, 2002, 31, 695]. Эффективной представляется модель Смагоринского на основе гипотезы Буссинеска с заданием в явном виде турбулентной подсеточной вязкости и модель Германо (SGS или динамическая модель Смагоринского с применениемдвух фильтров) [Germano, Piomelli et al., 1991; Абдибеков и др., 2007; Гарбарук А.В. и др., 2012; Юн А.А., 2009; Darryl D. Holm and Cesare Tronci., 2012].
Уравнения мелкой воды [Agoshkov V.I., Ambrosi D., Pennati V., Quarteroni A., Saleri F., 1993], используемые в проекте для моделирования динамики затопления территорий, широко применяются для решения самых различных задач динамики поверхностных вод. Предложены их различные модификации, связанные с учетом адвективного переноса [Karelsky K.V., Petrosyan A.S., 2006], с многослойными моделями [Чикин А.Л., Бирюков П.А., 2010], обобщениями без приближения гидростатического равновесия, допускающие эффективную численную реализацию [Bristeaua M.-O., Goutalb N., Sainte-Marie J., 2011; Wang K.-H., Li W., Lee H., 2008]. Наибольшее развитие получили гидродинамические модели для описания различных физических процессов как в отдельных областях морей и океанов, так и для всего водоема в целом [Аверкиев А.С., Клеванный К.А., 2007; Чикин А.Л., Бирюков П.А., 2010; Будинова Е.В., Носов В.Н., Терехин А.Т., 1987]. В рамках таких крупномасштабных моделей изучаются различного рода глобальные циркуляции [Аверкиев А.С., Клеванный К.А., 2007], изменчивость уровня, динамика солености [Фомин В.В., 2006] и ледяных заторов [Дебольский В.К., Дебольская Е.И., Масликова О.Я., 2010], волновые движения [Воеводин A.Ф., Никифоровская В.С., Остапенко B.В., 2008], цунами [Заибо Н., Пелиновский Б.Н., Храмушин В.Н., 2001], аварийные ситуации [Еремин М.А., Хоперсков А.В., 2006]. Решение широкого круга задач динамики поверхностных вод для самых различных приложений основывается на модели мелкой воды (уравнениях Сен-Венана) [Ю.И. Шокин, З.И. Федотова, Г.С. Хакимзянов, 2015] и ее модификациях [С. П. Баутин, С. Л. Дерябин, 2012; С. П. Баутин, С. Л. Дерябин, 2013; Е. Н. Пелииотский, 1996; В.А. Прокофьев, 2000; Z. I. Fedotova, G. S. Khakimzyanov, 2008; Z.I. Fedotova, G.S. Khakimzyanov and D. Dutykh, 2014; A.E. Green and P.M. Naghdi, 1976], включая многослойные модели [К.Н. Данилова, В.Ю. Ляпидевский, 2014; K.V. Karelsky, A.S. Petrosyan and A.G. Slavin, 2014], либо с учетом вертикальной структуры на примитивных уравнениях [Gayaz S Khakimzyanov, Oleg I Gusev, Sofya A Beizel, Leonid B Chubarov and Nina Yu Shokina, 2015], с учетом дисперсионных эффектов в различных приближениях [Ю.И. Шокин, З.И. Федотоваи Г.С. Хакимзянов, 2015; Р.Х. Зейтунян, 1995; V. Yu. Liapidevskii, 1998]. Укажем на направления, связанные с моделированием русловых объектов — водохранилищ [О.И. Гусев, Н.Ю. Шокина, В.А. Кутергин and Г.С. Хакимзянов, 2013] и речных систем [M. V. Bolgov, G. F. Krasnozhon and K. Yu. Shatalova, 2014; D. Caviedes-Voullieme, M. Morales-Hernandez, I. Lopez-Marijuan and P. Garca-Navarro, 2014; M.S. Horritt, G. Di Baldassarre, P.D. Bates and A. Brath, 2007; A.V. Pisarev, S.S. Khrapov, E.O. Agafonnikova and A.V. Khoperskov, 2013], крупных озер, дождевых потоков [Pierfranco Costabile, Carmelina Costanzo and Francesco Macchione, 2013], динамики плотных мелкозернистых геофизических потоков [C. Juez, D. Caviedes-Voullieme and P Murillo, J. and. Garca-Navarro, 2014], волн цунами [А.Г. Марчук and П.С. Мошкалев, 2014].
3D модель местности, построенная с высоким разрешением для больших площадей, позволяет решать различные проблемы в области наук о Земле. Эти проблемы связаны с многими приложениями в области охраны окружающей, моделирования и мониторинга. Отметим новые научные направления связанные с археологическими раскопками и реставрацией исторических ландшафтов [Lerma et al., 2010; Mozas-Calvache et al., 2012; Brutto et al., 2014] для крупномасштабных карт с сантиметровым пространственным разрешением [Azhar and Ahmad, 2014]. Рельеф местности с высоким пространственным разрешением имеет большое значение для моделирования динамики поверхностных вод [Wang et al, 2010; Хоперсков и др., 2012]. Кроме того, различные сезонные затопления поймы [Khrapov et al., 2013; Yuan et al., 2013] и катастрофические разливы рек [Bosa and Petti, 2013], наводнения из-за очень сильных дождей [Burguete et al., 2002] и цунами событий [Fernandes, 2009] должны быть приняты во внимание. Недавние трагедии в Крымске [Котляков др., 2013] 6-7 июля 2012 года, очень высокая вода в бассейне Амура [Данилов-Данильян и др., 2014], наводнение на Балканах в мае 2014 года, [Holt, 2014], являются примерами трагических событий и для того, чтобы избежать подобных катастроф в будущем, важно иметь высококачественные цифровые модели рельефа. Рельеф местности является основным фактором, который влияет на динамику воды. К сожалению, в последнее время точность лучших топографических карт является не достаточно высокой. Кроме того, появляются новые проблемы в небольших пространственных масштабах, связанных с изменениями на поверхности местности из-за природных и антропогенных причин [Воронин и др 2012; Ahmed et al., 2014]. Технология мобильного лазерного сканирования(MLS) представляет собой современный эффективный метод для построения трехмерных топографических карт [Kukko др., 2012]. Беспилотные летательные аппараты (БПЛА) используется для гладких поверхностей или арктических территорий, обеспечивающих сантиметровую точность [Stecchi др 2013; Lucieer др 2014]. Есть примеры успешного уточнения цифровой модели рельефа (ЦМР) с использованием MLS для небольших участков рек [Vaaja, 2011]. Однако MLS-подход имеет ограничения для пересеченной местности с большим количеством озер и ручьев на большой площади. Несмотря на успехи в классификации растительного покрова [Saarinen et al., 2013; Wallace et al., 2014] кусты и деревья ограничивают применение мобильного лазерного сканирования для получения высокоточной 3D-топографии. Различные локальные изменения местности, которые трудно обнаружить, могут существенно повлиять на качество гидрологических расчетов. БПЛА является важным современным инструментом для целого ряда проблем, таких как мониторинг поверхности и приповерхностного слоя атмосферы Земли, для решения проблемы стихийных бедствий, для построения цифровой модели рельефа [Jensen et al., 2009; Jensen et al., 2011; Maglione, 2014; Udin and Ahmad, 2014; Prakash et al., 2014], что обусловлено быстрым прогрессом в навигационных системах, обработке данных датчиков, контроля полета в режиме реального времени [Chiang et al., 2012] и мобильных фотограмметрических платформ [Baiocchi, 2014; Udin and Ahmad, 2014]. В настоящем проекте развивается оригинальный подход построения ландшафта ВАП, который основан на наблюдении за динамикой береговой линии для большого числа водоемов во время весеннего паводка. Данный подход лишен недостатков MLS метода при построении ЦМР на пересеченной местности с большим количеством озер и ручьев и позволяет строить ЦМР высокого разрешения для важных в плане гидрологии областей ВАП.
При наличии твердых границ мы имеем возможность корректного описания жидкости, как для Эйлеровых численных схем, так и в случае Лагранжевого подхода, применяя, например, метод модифицированных виртуальных граничных частиц (Modifed Virtual Boundary Particle) [R. Vacondio, B. D. Rogers and P. K. Stansby, 2012]. В случае озер или морей, если водный объект полностью находится в пределах расчетной сетки, проблема граничных условий отсутствует, если применять алгоритмы «движения по сухому дну» [S.P. Bautin, S.L. Deryabin, A.F. Sommer, G.S. Khakimzyanov and N.Yu. Shokina, 2011]. В случае протяженных систем (водохранилищ и рек) практически всегда имеются границы, через которые вода поступает в расчетную область и уходит из нее, что требует аккуратной постановки граничных условий [J. Burguete and P. Garcia-Navarro, 2004; J. Burguete, P. Garcia-Navarroand, J. Murillo, 2006]. Наиболее сложным оказывается режим, при котором граница «жидкость – сухое дно» достигает границы расчетной области. Такая ситуация имеется при численном исследовании динамики поверхностных стоков воды, возникающих при затоплении территории, движении мелкодисперсных потоков вещества [C. Juez, D. Caviedes-Voullieme and P Murillo, J. and. Garca-Navarro, 2014], особенно на сложном горном рельефе [J. Burguete, P. Garcia-Navarro and R. Aliod, 2012], например, в случае прорыва моренно-подпрудных приледниковых и ледниковых озерных систем (of moraine-dammed supraglacial and proglacial lake systems) [M.J. Westoby, N.F. Glasser, J. Brasington, M.J. Hambrey, D.J. Quincey and J.M. Reynolds, 2014], из-за дождевых осадков [J. Singh, M. Altinakar and Y. Ding, 2014]. Сходная проблема возникает при использовании уравнений мелкой воды для моделирования динамики атмосферы [Yu. N. Skiba, 1995]. При прохождении вещества через границу помимо численной устойчивости расчета важным является выполнение законов сохранения и построение неотражающих условий (nonrefecting conditions) [M.A. Ilgamov and A.N. Gilmanov, 2003]. Хорошо известна проблема описания нестационарной линии уреза воды, связанная с необходимостью находить решение в области с подвижной границей (например, [С. П. Баутин, С. Л. Дерябин, А. Ф. Соммер и Г. С. Хакимзянов, 2010; С. П. Баутин и С. Л. Дерябин, 2012] для различных моделей мелкой воды [Jean-Marie Zokagoa and Azzeddine Soulaimani, 2010]). В численных подходах при ее решении используются различные методы [Qiuhua Liang and Alistair GL Borthwick, 2009]. В случае достаточно гладкой функции дна b(x,y) можно успешно использовать регуляризацию [S. Vater, N. Beisiegel and J. Behrens, 2015]. Использование модифицированного закона сохранения полного импульса с учетом образования локальных турбулентно-вихревых структур в поверхностном слое воды позволяет получить удовлетворительное согласие с экспериментальными данными [V.V. Ostapenko, 2007]. Полученные результаты аналитического исследования решений использованы для разработки новых аппроксимаций краевых условий на подвижной линии уреза. Однако, учитывая нерегулярный характер изменений реального дна b(x,y) на различных масштабах возникает необходимость строить более сложные алгоритмы, специально предназначенные для моделирования динамической границы «жидкость – сухое дно» [С.С. Храпов и др., 2011]. В рамках проекта проведено актуальное исследование влияния различных граничных условий в численной модели мелкой воды на структуру руслового потока для ситуаций, требующих аккуратной постановки граничных условий: 1) жидкость свободно проходит через границу в докритическом режиме; 2) вблизи свободной границы имеется участки сильно неоднородного рельефа; 3) возникает поток жидкости внутрь расчетной области; 4) происходят существенно нестационарные процессы вблизи границы.
Развиваемый в проекте оригинальный численный метод CSPH-TVD основан на совместном использовании лагранжевого [Monaghan J.J., 1985; Monaghan J.J., 1992] и эйлерового подходов [van Leer B., 1979; Harten A., 1983], применяемых в современных численных алгоритмах для моделирования динамики сплошных сред. По точности, устойчивости и эффективности CSPH-TVD метод не уступают мировым аналогам [Titarev V.A., Toro E.F., 2007; Castro M.J., Pares C., Pardo A., Toro E.F., 2008; Toro E.F., Hidalgo A., Dumbser M., 2009; Kallemov B., Miller G.H., Trebotich D., 2010; Colella P., Dorr M.R., Hittinger J.A.F., Martin D.F., 2011; Guzik S., McCorquodale P., Colella P., 2012]. CSPH–TVD-метод обладает большей эффективностью и устойчивостью по сравнению с другими численными алгоритмами, применяемыми для моделирования сложных гидродинамических течений в природных и технических системах при наличии сильно неоднородных на мелких масштабах нерегулярных гидродинамических и внешних сил и нестационарных границ «газ/жидкость-вакуум». Подтверждением данного утверждения является успешное и эффективное применение CSPH–TVD метода для описания динамики поверхностных вод на нерегулярном рельефе местности, содержащем изломы и резкие перепады уровней [Храпов С.С. и др. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD-подхода // Вычисл. мет. и программ. 2011. Т. 12. №1. С.282-297]. Построенная для уравнений мелкой воды (Сен-Венана)численная схема является хорошо сбалансированной, консервативной и позволяет проводить сквозной устойчивый расчет при наличии нестационарных границ «вода–сухое дно» на нерегулярном рельефе. Отметим, что возможны различные варианты реализации предложенного алгоритма CSPH-TVD. В численной схеме можно использовать различные сглаживающие ядра, TVD–ограничители (сохраняющие TVD–свойство разностной схемы), методы решения задачи Римана, а также учитывать нестационарные источники/стоки воды, внешние силы и нестационарность рельефа дна. Кроме того, данный алгоритм может быть обобщен на случай неравномерных сеток. Преимуществом рассматриваемого CSPH-TVD метода в случае нерегулярного и/или разрывного рельефа является хорошая точность моделирования, сравнимая с точностью численных схем годуновского типа, основанных на точном решении обобщенной задачи Римана, при меньших вычислительных затратах и более простых расчетных формулах численного алгоритма. По сравнению с различными модификациями SPH-метода численная схема CSPH-TVD имеет более высокую точность, является менее диссипативной и существенно лучше сбалансированной, а также обладает большей вычислительной скоростью при равном количестве частиц. Данные особенности рассматриваемого алгоритма позволяют применять его для моделирования динамики поверхностных вод на значительных территориях с хорошим пространственным разрешением, что необходимо для решения широкого круга задач динамики поверхностных вод, связанных с сезонными затоплениями, проведением экспертиз гидротехнических сооружений, проблемами загрязнения водоемов. Численный алгоритм CSPH-TVD способен описывать сдвиговые неустойчивости, несмотря на то, что на одном из этапов схемы CSPH-TVD используется подход сглаженных частиц SPH [Писарев А.В., Храпов С.С., Хоперсков А.В. Численнаясхема на основе комбинированного подхода SPH-TVD: проблема моделирования сдвиговых течений // Вестник ВолГУ. Сер.1: Математика. Физика. 2011. Т.15. №2. С. 138-141]. Таким образом, метода CSPH-TVD отсутствует недостаток классического SPH алгоритма, не позволяющий корректно моделировать тангенциальные и контактные разрывы.
Современная вычислительная техника позволяет решать самые различные задачи гидрологии с учетом реального рельефа местности и всех значимых физических факторов для больших территорий [Писарев А.В., Храпов С.С., Агафонникова Е.О., Хоперсков А.В. 2013]. Особо выделим проблему управления гидрологическим режимом пойменных ландшафтов во время весеннего паводка на крупных реках [Воронин А.А., Елисеева М.В., Храпов С.С., Писарев А.В., Хоперсков А.В. 2012]. Решение этой задачи требует максимально эффективных численных моделей на основе параллельных технологий [Glinskiy B., Kulikov I., Snytnikov A., Romanenko A., Chernykh I., Vshivkov V. 2014]. Для эколого-экономического управления территорией, например, Волго-Ахтубинской поймы, необходимо решать задачу оптимизации гидрографа для конкретных условий в каждом году, и делать экспертные заключения для различных режимов работы десятков гидросооружений в пойме и по строительству новых [Воронин А.А., Васильченко А.А., Писарева М.В., Писарев А.В., Хоперсков А.В., Храпов С.С., Подщипкова Ю.Е. 2015]. Каждая такая задача требует проведения сотен численных экспериментов по прямому гидродинамическому моделированию весеннего гидрологического режима для территории площадью от 2 тыс. до 20 тыс. кв. км.
Наша практика использования больших суперкомпьютеров (в частности, на вычислительных ресурсах Суперкомпьютерного комплекса МГУ [Sadovnichy V., Tikhonravov A., Voevodin V., Opanasenko V. 2013]) для проведения массовых гидродинамических расчетов выявила ряд проблем, связанных с необходимостью проводить за короткий период большое число расчетов и передачей через сеть больших объемов данных для последующей обработки и анализа [9]. Оперативность проведения расчетов и обработки данных моделирования являются важнейшим фактором применения такого рода моделей в практической деятельности. Дополнительной проблемой является визуализация проведенных расчетов, что представляется общей трудностью для очень высокопроизводительных машин [Moreland K., Larsen M., Childs H. 2015]. Переход на вычислительные ресурсы класса "персональные суперкомпьютеры" на основе GPU позволяет отчасти решить указанные проблемы. В работе рассмотрены результаты создания пакета программ для параллельных расчетов на вычислительных узлах NVIDIA TESLA C2070, K20, K40, K80.
Для эффективного использования в составе проекта систем имитационного моделирования применяются технологии параллельных вычислений, поскольку при проведении гидродинамических расчетов с высоким разрешением на обширных территориях в двумерных и трехмерных моделях необходимы значительные вычислительные ресурсы. Современные широко распространенные технологии параллельных вычислений - OpenMP (системы с общей памятью CPU), MPI (системы с разделенной памятью CPU) и CUDA (графические ускорители NVIDIA GPU) [Антонов А.С., 2002; Букатов А.А. и др., 2003; Антонов А.С., 2009; Gorobets A. at all, 2011; Деги Д.В., Старченко А.В., 2012; Горобец А.В. и др., 2011; Боресков А.В., Харламов А.А., 2011; СандерсДж., Кэндрот Э., 2011; Казѐннов А.М., 2010].

  Методы и подходы, использованные в ходе выполнения Проекта
 

Важным фактором, определяющим динамику поверхностных вод и влекомых наносов, является придонное трение, интенсивность которого зависит от коэффициента шероховатости по Маннингу, который, в свою очередь, зависит от свойств подстилающей поверхности дна, извилистости и структуры русла, мелкомасштабных неоднородностей рельефа дна, турбулентности течения, диссипации энергии потока воды за счет перемещения донных наносов. 
В рамках проекта для расчета эффективного коэффициента Маннинга будут использованы оригинальные подходы, основанные на сравнении результатов численного гидродинамического моделирования в меандрированном и прямом руслах с гладким и содержащим мелкомасштабные неоднородности дном [2, 3]. С усилением степени извилистости русла и амплитуды неоднородностей дна средняя скорость движения жидкости в сечении канала уменьшается. Таким образом, варьируя основное значение коэффициента Маннинга в прямом канале на основе численных гидродинамических моделей, можно оценить вклад в сопротивление потоку от меандрирования русла и мелкомасштабных неоднородностей рельефа дна. 
Развиваемый в рамках проекта численный метод CSPH-TVD (комбинированная схема Smoothed Particle Hydrodynamics – Total Variation Diminishing) основан на совместном использовании лагранжевого и эйлерового подходов, применяемых в современных численных алгоритмах для моделирования динамики сплошных сред. Численные схемы, основанные на данных подходах, обладают как преимуществами, так и недостатками при проведении гидродинамических расчетов.
К достоинствам численных схем, основанных на лагранжевом подходе, можно отнести достаточную простоту/универсальность реализуемого алгоритма и корректное описание динамической границы «вода - сухое дно», что очень важно задач моделирования динамики жидкости со свободной поверхностью. Недостатками схем моделирующих сплошную среду лагранжевыми частицами является невысокая точность вычислений, сильная диссипация (необходимая для повышения устойчивости), слабая сбалансированность получаемого решения, большие вычислительные затраты в случае высокой концентрации частиц, наличие на разрывах осцилляций численного решения, что не позволяет, например, корректно моделировать ударные волны, тангенциальные и контактные разрывы.
Наибольшую распространенность в задачах вычислительной гидродинамики (Computational Fluid Dynamics – CFD) получили методы, основанные на эйлеровом (сеточном) подходе. К настоящему времени разработан целый ряд эйлеровых схем, обеспечивающих высокую точность, устойчивость и отсутствие на разрывах «паразитных» осцилляций при моделировании динамики жидкости, однако при использовании данных методов в ряде случаев возникают следующие проблемы [Храпов С.С., и др. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD-подхода // Вычисл. методы и программ. 2011. Т. 12. №1. С.282-297]: 
1) для корректного описания динамической границы «вода - сухое дно» требуется либо дополнительная регуляризация, что может вносить большую и неконтролируемую погрешность, либо построение точного решения задачи Римана, приводящее к дополнительным вычислительным затратам;
2) при наличии в уравнениях мелкой воды, записанных в консервативной (дивергентной) форме, отличной от нуля правой части (внешние и вязкие силы, источники и стоки вещества и импульса) необходимо строить решение задачи Римана с учетом этих дополнительных слагаемых, что является весьма трудоемкой задачей;
3) для построения хорошо сбалансированных (well balanced) схем требуются дополнительные модификации алгоритма, также приводящие к дополнительным вычислительным затратам;
Одной из целей проекта является создание эффективного численного алгоритма обладающего достоинствами лагранжева и эйлерова подходов и сводящего к минимуму недостатки обоих подходов при гидродинамическом моделировании совместной динамики поверхностных и грунтовых вод, размыва дамб и влекомых наносов на гибридных суперкомпьютерах с графическими ускорителями.
Участниками проекта ранее был предложен лагранжево-эйлеров подход для описания динамики поверхностных вод на нерегулярном рельефе местности [Храпов С.С., и др. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD-подхода // Вычисл. методы и программ. 2011. Т. 12. №1. С.282-297], [14], [1].
В CSPH-TVD-методе расчетная область покрывается неподвижной (эйлеровой) сеткой в центры ячеек которой в начальный момент времени помещаются «жидкие» частицы. Данная схема кардинально отличается от идеологии классического метода частиц в ячейках [Харлоу Ф.Х. М.: Мир, 1967, с.316], а также различных его модификаций [Петренко В.Е., и др. // Числ. методы мех. сплош. среды. 1976. Т.7. № 4. С.130-148; Ким В.В., и др. // Математ. моделир. 2006. Т. 18. № 8. С.5-11].
Основная идея предлагаемого метода заключается в том, чтобы на первом лагранжевом (частицы) этапе учитывать силы (градиент давления, внешние и вязкие) и источники/стоки массы и энергии, а на втором эйлеровом (сетка) этапе рассчитывать потоки массы, импульса и энергии, обусловленные перемещением жидкости через границы ячеек (адвекция).
На первом лагранжевом этапе CSPH-TVD-метода используется SPH подход [Monaghan J.J. // Comp. Phys. reports. 1985. V.3, N2. P.71–124; Monaghan J.J. // Ann. Rev. Astron and Astroph. 1992. V.30. P.543-574], модифицированный нами для обеспечения хорошей сбалансированности и уменьшения численной диссипации алгоритма.
Классический SPH-метод основан на дискретизации сплошной среды конечным набором лагранжевых частиц, движущихся со скоростью потока и допускающих произвольную связность между собой. В данном подходе частицы представляются не в виде точечных объектов, а как «жидкие» частицы, непрерывно распределенные в пространстве и имеющие конечный объем. Пространственное распределение величин, характеризующих жидкие частицы (плотность, энергия), представляется с использованием сглаживающей функции W. В качестве сглаживающего ядра W применяются сплайн-функции различных порядков или Гауссово ядро.
Повышение порядка аппроксимации (по времени и пространству) CSPH-TVD-метода на лагранжевом этапе обеспечивается посредством:
1) применения различных схем интегрирования по времени (Рунге-Кутта, прогноза-коррекции) сохраняющих TVD-свойство метода;
2) увеличения порядка полиномиального представления сглаживающего ядра W.
На втором эйлеровом этапе CSPH-TVD-алгоритма используются сеточные методы удовлетворяющие принципу TVD [Harten A. // J. Comput. Phys. 1983. V.49. N3. P.357–393]. Выполнение условия TVD осуществляется посредством применения функций-ограничителей, не позволяющих антидиффузионным слагаемым разностной схемы становиться слишком большими [Sweby P.K. // SIAM J. Numer. An. 1984. V.21. N5. P.995–1011]. Таким образом, можно строить схемы второго порядка точности и выше, у которых отсутствуют или существенно уменьшены паразитные осцилляции вблизи разрывов.
Хорошо апробированным вариантом реализации CSPH-TVD-метода является использование на эйлеровом этапе подхода MUSCL (Monotonic Upstream Schemes for Conservation Laws) [van Leer B. // J. Comput. Phys. 1979. V.32. N1. P.101–136]. Данный подход является развитием схемы Годунова, основанной на использовании решения задачи Римана (RP – Riemann Problem) при вычислении потоков через границы расчетных ячеек. В отличие от оригинальной схемы Годунова, основанной на кусочно-постоянном распределении вычисляемых величин внутри ячеек, в методе MUSCL решение задачи Римана строится для кусочно-линейного распределения параметров (GRP – Generalized Riemann Problem). Обобщенное решение задачи Римана (GRP) может быть построено, однако оно слишком сложно для вычислений. Поэтому в данном подходе ограничиваются коррекцией значений параметров слева и справа от границы расчетной ячейки, которое дает учет кусочно-линейного распределения параметров [van Leer B. // SIAM J. Sci. Stat. Comput. 1984. V.5. N.1. P.1–20]. Развитием такого подхода стало построение кусочно-параболического метода (Piecewise Parabolic Method — PPM) [Colella P., Woodward P.R. // J. Comput. Phys. 1984. V.54. N.1. P.174–201]. 
Поскольку решение задачи Римана (Riemann solver) является важным моментом при построении такого рода схем, то на данном этапе используются как точные, так и приближенные методы (повышающие производительность вычислений) решения задачи Римана (Куранта-Изаксона-Риса, Роу, Ошера, Хартена-Лакса-ван Лира и др.) [Roe P.L. // J. Comput. Phys. 1981. V.43. N.2. P.357–372; Harten A., Lax P.D., van Leer B. (1983) // SIAM Review. 1983. V.255. N.1. P.32–61; Toro E.F., Chakroborthy A. // Aeronaut. J. 1994. V.98. N.979. P.325–332; Toro E.F. // Phylos. Trans. Royal Soc. London. 1989. Ser A341. N.1662. P.499–530]. Отметим также, что для сохранения свойства монотонности или удовлетворения требованию TVD для наклонов распределения параметров в ячейках используются различные ограничители.
Повышение порядка аппроксимации (по времени и пространству) CSPH-TVD-метода на эйлеровом этапе обеспечивается посредством:
1) использования для расчетов потоков усредненных по времени величин, рассчитанных на первом лагранжевом этапе;
2) увеличения порядка полиномиальной реконструкции распределений вычисляемых величин внутри расчетных ячеек с использованием стандартных процедур или с помощью сглаживающего ядра W.
3) использования либо точного решение задачи Римана (GRP), либо приближенного решения, являющегося нелинейным и полным (в решении учитываются ударные волны, волны разряжения, контактные и тангенциальные разрывы).
Благодаря первому лагранжеву этапу CSPH-TVD-метода при наличии в уравнениях гидродинамики источниковых (правая часть уравнений) слагаемых нет необходимости искать решение задачи Римана с учетом этих слагаемых (в общем случае это может оказаться очень трудоемкой и затратной процедурой), а ограничиться (без потери точности) уже существующими апробированными решениями. Кроме этого, CSPH-TVD-подход позволяет строить новые решения задачи Римана на основе полученных после лагранжева этапа пространственно-временных распределениях вычисляемых величин.
Исследование устойчивости и точности разрабатываемого в рамках проекта численного алгоритма на основе CSPH-TVD-метода будет осуществляться посредством сравнения численного решения с точным на сжимающихся сетках для различных наборов гидродинамических тестов: распад произвольного разрыва глубины и/или скорости течения; распад столба жидкости на сухом дне; распространение малых возмущений поверхностных гравитационных волн при наличии неоднородного дна; исследование процесса длительного распространения солитона, позволяющее оценить дисперсионные свойства численной схемы; расчет гидравлического прыжка; расчет стационарных и нестационарных течений на существенно неоднородном рельефе дна, включая разрывные профили; накат волн типа цунами на пологий берег; колебания уровня жидкости на параболическом дне; течение типа водопад; сдвиговые течения; распад циркуляционной дамбы на нерегулярном рельефе. Верификации численной модели по переносу наносов и размыва берегов будет осуществляться с использованием экспериментальных данных и утвержденных методик расчета размыва дамб.
Важной характеристикой численных схем является эффективность – точность по отношению к вычислительным затратам. Более эффективной будет схема, для которой при заданной точности время вычислений меньше [Titarev V.A., Toro E.F. // J. of Numer. Analysis. 2007. V.27. P.616–630]. Оценка эффективности численного алгоритма на основе CSPH-TVD метода будет проводиться по отношению к порядку аппроксимации и различным методам решения задачи Римана. 
Практическое использование численных схем связано с проведением гидродинамических расчетов с высоким разрешением в двумерных и трехмерных моделях, требующих применения технологий параллельных вычислений. Для оценки эффективности распараллеливания численных алгоритмов CSPH-TVD-метода будут использованы технологии параллельных вычислений OpenMP (системы с общей памятью CPU), MPI (системы с разделенной памятью CPU) и CUDA (графические ускорители NVIDIA GPU) [Антонов А.С. – М.: Изд-во МГУ. 2002; Букатов А.А. и др. – Ростов-на-Дону: Издательство ООО «ЦВВР». 2003; Антонов А.С. – М.: Издательство МГУ. 2009; Gorobets A. at all // Computer & Fluids. 2011. V. 49. №1. P. 101-109; Деги Д.В., Старченко А.В. // Вестник Томского гос. университета. 2012. №2 (18). С.88-98; Горобец А.В. и др. // Вестник Южно-Уральского государственного университета. Серия: Математическое моделирование и программирование. 2011. №25(242). С.76-86; Боресков А.В., Харламов А.А. – М.: ДМК Пресс, 2011. – 232 с.; Сандерс Дж., Кэндрот Э. – М.: ДМК Пресс. 2011. – 232 с.; Казённов А.М. // Комп. исслед. и моделир. 2010. Т.2. №3. С.295-308].
Наиболее перспективным направлением, обеспечивающим высокую производительность при моделировании динамики жидкости и газа, является адаптация вычислительного параллельного кода для гибридных суперкомпьютеров на базе графических процессоров (GPU) NVIDIA, используя технологию CUDA. Для гидродинамических расчетов планируется использовать двухуровневые и трехуровневые гибридные схемы распараллеливания численных алгоритмов OpenMP-CUDA/OpenCL, MPI-OpenMP-CUDA/OpenCL [14], [Богданов П.Б., Горобец А.В., Суков С.А. Журнал вычислительной математики и математической физики. 2013], [9]. Двухуровневая гибридная схема OpenMP-CUDA предназначена для работы на гибридных вычислительных системах с общей памятью, в составе которых имеется несколько GPU. Для проведения расчетов на вычислительном кластере, состоящем из нескольких гибридных систем с общей памятью, целесообразно использовать трехуровневую гибридную схему распараллеливания MPI-OpenMP-CUDA. При использовании нескольких GPU необходимо применять технологию NVIDIA GPUDirectTM (Direct Access, Direct Transfers) для увеличения производительности гидродинамических расчетов.
Для детального описания мелкомасштабных структур будут применяться адаптивные сетки и технология Zoom-in Simulation (расчет внутри расчета). При моделировании в рамках Zoom-in подхода на первом этапе расчет ведется глобально для всей области. На следующем этапе рассчитываются заданные локальные области, с учетом граничных и начальных условий согласованных с предыдущим этапом. В зависимости от конкретной задачи и учета физических механизмов может использоваться несколько таких вложенных расчетов. Данный подход, например, используется в космологических расчетах [например, Semelin, Combes 2005, A&A, 441, 55; Governato et al 2009, MNRAS, 398, 312; Scannapieco et al 2009, MNRAS, 396, 696], имеются примеры использования zoom-in для локальных моделей динамики молекулярной фазы [Haugboelle et al., 2011, LPI Contr. N1639, 9116].

  Вклад каждого члена коллектива в выполнение Проекта в 2018-2019 годах
 

1) Построение самосогласованной математической модели динамики поверхностных вод и влекомых наносов, основанных на уравнениях Сен-Венана и Экснера. (Храпов С.С., Хоперсков А.В.)

2) Разработка численной модели на основе CSPH-TVD метода для решения уравнений самосогласованной динамики поверхностных вод и влекомых наносов. (Храпов С.С., Кузьмин Н.М.)

3) Исследование устойчивости численного алгоритма в зависимости от параметров математической модели. (Храпов С.С., Кузьмин Н.М.)

4) Верификация численной модели с использованием известных точных решений уравнений, экспериментальных данных и утвержденных методик расчета размыва дамб. (Храпов С.С., Хоперсков А.В.)

5) Разработка компьютерного алгоритма численной модели для гибридных вычислительных систем с сопроцессорами Nvidia Tesla K20/40/80 на основе технологий параллельных вычислений MPI-OpenMP-CUDA и GPUDirect. (Дьяконова Т.А., Храпов С.С.)

6) Создание программного комплекса (модуль подготовки и запуска расчетов, расчетный модуль CPU+GPU, модуль визуализации результатов моделирования), позволяющего проводить быстрое имитационное моделирование русловых процессов и динамики затопления на обширных территориях с высоким пространственным разрешением. (Сиволобов С.В., Дьяконова Т.А., Храпов С.С., Астахов А.С.)

7) Исследование эффективности распараллеливания алгоритма для различных конфигураций вычислительных систем CPU+GPU. (Дьяконова Т.А., Храпов С.С.)

8) Анализ водных объектов (реки, каналы, озера, водохранилища и т.п.) Волгоградской области и определение потенциально опасных и наиболее значимых для гидрологического режима территорий объектов. (Агафонникова Е.О., Радченко В.П.)

9) Построение актуализированной цифровой модели исследуемого участка местности (матрицы пространственного распределения абсолютных высот рельефа, коэффициентов шероховатости, инфильтрации, пористости, медианной крупности и плотности частиц грунта), включающего наиболее значимые и потенциально опасные территории с водными объектами Волгоградской области. (Агафонникова Е.О., Радченко В.П., Астахов А.С.)

10) Имитационное гидродинамическое моделирование русловых процессов, размыва и последствий прорыва заградительных дамб (плотин) в период сезонного и/или ливневого паводка для исследуемого участка местности включающего наиболее значимые и потенциально опасные территории с водными объектами Волгоградской области и анализ результатов моделирования. (Агафонникова Е.О., Радченко В.П., Дьяконова Т.А., Хоперсков А.В., Храпов С.С.)

  Количество научных работ по Проекту за 2018-2019 года
 

Всего: 3

Из них в изданиях, включенных в перечень ВАК: 1

Из них в изданиях, включенных в библиографическую базу данных РИНЦ: 2

Из них в изданиях, включенных в международные системы цитирования: 2

  Участие в 2018-2019 годах в научных мероприятиях по тематике Проекта
 

Доклад (постерный) на V МЕЖДУНАРОДНОЙ КОНФЕРЕНЦИИ И МОЛОДЕЖНОЙ ШКОЛЕ "ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ И НАНОТЕХНОЛОГИИ", Самара, 21-24 мая 2019 г.

  Библиографический список всех публикаций по Проекту
 

1) Klikunova A.Yu., Khoperskov A.V., Agafonnikova E.O., Kuz’mich A.S., Dyakonova T.A., Khrapov S.S., Gusev I.M. Creation of cadastral maps of flooding based on numerical modeling // Journal of computational and engineering mathematics, 2019, vol. 6, no. 2, pp. 3-17. DOI: 10.14529/jcem190201. https://jcem.susu.ru/jcem/article/view/195/152

2) Klikunova A.Yu., Khoperskov A.V. Creation of digital elevation models for river floodplains // CEUR Workshop Proceedings, 2019, vol.2391, pp. 275-284. http://ceur-ws.org/Vol-2391/paper38.pdf

3) Кликунова А.Ю., Хоперсков А.В. Построение цифровых моделей рельефа местности для речных пойм / Сборник трудов ИТНТ-2019 под редакцией В.В. Мясникова. 2019. Издательство: Новая техника. С. 110-115. https://elibrary.ru/item.asp?id=37344175