CC BY
15УДК 536.2:519.87:622.276
ПАРАЛЛЕЛЬНЫЕ МЕТОДЫ И ТЕХНОЛОГИИ РЕШЕНИЯ ЗАДАЧ НЕФТЕГАЗОВОЙ
ОТРАСЛИ
В. П. Ильин
Институт вычислительной математики и математической геофизики СО РАН, Новосибирский государственный технический университет, ilin@corp.nstu.ru
Рассмотрены фундаментальные и прикладные вопросы высокопроизводительного моделирования процессов и явлений в актуальных проблемах нефтегазовой отрасли. Исследуются стратегические и тактические особенности распараллеливания различных технологических стадий крупномасштабного вычислительного эксперимента на многопроцессорных суперкомпьютерах гетерогенной архитектуры с распределенной и иерархической общей памятью. Особое внимание уделено реализации многомасштабных разрывных методов Галеркина высоких порядков для решения междисциплинарных прямых и обратных задач с реальными данными, включая сложные геометрические конфигурации многосвязных кусочно-гладких границ и контрастных материальных свойств в различных средах, а также алгоритмам решения сверхбольших разреженных СЛАУ. Описаны технологии гибридного программирования с использованием средств передачи межпроцессорных данных, многопотоковых вычислений, векторизации машинных операций и графических ускорителей. Проанализированы характеристики и показатели эффективности распараллеливания алгоритмов и их отображения на архитектуру компьютера, а также вопросы создания математического и программного обеспечения на компонетных принципах формирования больших интегрирующих окружений.
Ключевые слова: масштабируемое распараллеливание, гибридное программирование, многопроцессорные гетерогенные системы, многомерные междисциплинарные прямые и обратные задачи, коммуникационные потери, декомпозиция областей.
PARALLEL METHODS AND TECHNOLOGIES FOR SOLVING OIL&GAS INDUSTRY
PROBLEMS
V. P. Ilyin
Institute for Computational Mathematics and Mathematical Geophysics, Siberian Branch, RAS, Novosibirsk State Technical University, ilin@corp.nstu.ru
The paper deals with basic and applied aspects of high-performance process simulation for oil&gas industry. We study large- and low-scale approaches to managing large parallel simulations with multi-CPU heterogamous computers with distributed hierarchical common RAM. The paper focuses of the implementation of multiscale discontinuous Galerkin high order methods to solve interdisciplinary direct and inverse real-life problems such as complex geometric configurations of multiply connected piecewise-smooth boundaries and contrast material properties in various environments, and algorithms for solving super-large sparse linear equation systems. Hybrid programming technologies a presented that use interCPU data transfer, multi-thread computing, vectorization, and graphics accelerator chips. We also analyze the properties and performance of algorithm paralleling and its mapping to the computer architecture, component-based software development for creating large integrated environments.
Keywords: scalable paralleling, hybrid programming, multiCPU heterogeneous systems, multidimensional interdisciplinary direct and inverse problems, communication losses, area decomposition.
Введение
Математические проблемы нефтегазовой отрасли представляют собой широкий класс междисциплинарных прямых и обратных задач, требующий для проведения предсказательного моделирования с высоким разрешением результатов экстремальных вычислительных ресурсов. Это относится как к фундаментальным исследованиям глубинных процессов и явлений, так и к оптимизации эксплуатационных режимов производственного оборудования. Требования принятия оперативных решений по расчетным данным очевидным образом подразумевают массивное распараллеливание алгоритмов на
современных многопроцессорных вычислительных системах (МВС) с распределенной и иерархической общей памятью. Достижение высокой производительности крупномасштабного вычислительного эксперимента требует квалифицированного отображения структуры алгоритмов на архитектуру суперкомпьютера, а создание эффективного математического программного обеспечения с длительным жизненным циклом подразумевает еще и возможности расширения состава моделей и алгоритмов, а также адаптацию к эволюции машинных платформ. Проблемы моделирования в целом затрагивают широкий круг вопросов [1].
Масштабируемое распараллеливание можно рассматривать в сильном или слабом смысле. Первое означает примерное сохранение времени решения задачи при прямо пропорциональном увеличении числа степеней свободы в алгоритме и количества используемых компьютерных процессоров, а второе — обратно пропорциональное уменьшение времени счета для фиксированной задачи при росте числа процессоров.
Основным принципом ускорения параллельных вычислений является сокращение объемов межпроцессорных обменов данными в сравнении с количеством выполняемых арифметических операций. Вопрос заключается не только в том, что коммуникации намного медленнее непосредственных расчетов, но они еще являются и наиболее энергозатратными, а последний фактор становится критическим при создании МВС экзафлопсной производительности. Данные технические обстоятельства имеют то математическое следствие, что в численных методах приоритет обретают алгоритмы высокой точности. Например, методы конечных элементов (МКЭ) четвертого и большего порядков позволяют на трехмерной сетке с числом узлов порядка 1003 достичь намного меньшей погрешности, чем «стандартные» алгоритмы 1-го или 2-го порядков на сетке 1 0003, а соответствующее сокращение размеров числовых массивов ведет к резкой экономии межпроцессорных обменов, что с лихвой компенсирует увеличение количества арифметических действий в расчете на один узел, которое неизбежно при усложнении аппроксимационной схемы.
Главные современные подходы к решению реальных задач моделирования заключаются в решении систем уравнений в частных производных (как правило, в вариационных или слабых постановках) в многомерных расчетных областях со сложной конфигурацией многосвязных кусочно-гладких граничных поверхностей и контрастными свойствами материальных сред.
Важно отметить, что математические проблемы нефтегазовой отрасли в целом представляют собой междисциплинарные прямые и обратные задачи, в которых, с одной стороны, требуется моделировать процессы различной физической природы, а с другой — решать уравнения не с полностью заданными исходными данными, а зависящими от неизвестных параметров, которые надо оптимизировать по некоторым дополнительным условиям. Характерным примером здесь является проблема идентификации параметров модели на основе сопоставления расчетных результатов и данных натурных измерений, что формально определяется как условная минимизация описываемого заранее целевого функционала. Если он имеет не один минимум, а несколько, то мы приходим к задаче не локальной минимизации, а глобальной (что является гораздо более трудоемкой процедурой).
Стратегии и тактики распараллеливания определяются по-своему на разных технологических стадиях компьютерного эксперимента: геометрическое и функциональное моделирование, построение адаптивных сеток, аппроксимация исходных уравнений, решение алгебраических систем и собственных проблем, методы оптимизации для решения обратных задач, постпроцессинг и визуализация результатов, управление вычислительным процессом и принятие решения по результатам расчетов.
Основным инструментом распараллеливания алгоритмов решения многомерных краевых задач является аддитивный метод декомпозиции областей (МДО), заключающийся в сбалансированном разбиении сеточной расчетной области на сеточные подобласти с параметризованной шириной пересечений и со специальными интерфейсными условиями на внутренних смежных границах. При этом исходная «большая» задача распадается на подзадачи, каждая из которых должна решаться синхронно на «своем» процессоре (с помощью прямых или итерационных алгоритмов). Внешние итерации по подобластям осуществляются блочными методами Шварца — Якоби в подпространствах Крылова, ускоряемыми путем грубосеточной коррекции, или малоранговой аппроксимации исходной матрицы.
Целью данной работы является анализ технологий распараллеливания алгоритмов на наиболее характерных математических постановках, и мы не останавливаемся на многих актуальных проблемах (гидроразрыв пласта, электромагнитные и/или химические воздействия и т. д.), для которых формули-
руемые нами выводы и рекомендации остаются справедливыми.
В статье представлены характерные математические модели решаемых междисциплинарных прямых и обратных задач, описаны принципы построения неявных аппроксимаций начально-краевых задач высокого порядка с помощью многомасштабных разрывных методов Галеркина, а также итерационные методы декомпозиции для решения систем линейных алгебраических уравнений (СЛАУ), которые являются «узким горлышком» моделирования «больших задач», поскольку на данной стадии объем требуемых вычислительных ресурсов растет нелинейно с увеличением общего количества степеней свободы (т. е. размерности СЛАУ). Изложены стратегии и тактики распараллеливания на различных этапах моделирования, а в заключении обсуждаются некоторые открытые вопросы и возможные направления дальнейших исследований по решению рассматриваемых проблем.
Описание математических моделей
Математические модели задач нефтегазовой отрасли достаточно хорошо описаны в литературе [2]. В данной работе для иллюстрации проблем распараллеливания вычислительных методов и технологий мы остановимся на междисциплинарной проблеме, включающей нелинейную систему векторных уравнений Дарси для многофазной фильтрации и скалярное уравнение теплопроводности с фазовыми переходами.
Двухфазная задача теплопроводности. Рассмотрим сначала проблему моделирования температурных полей с фазовыми переходами, которая в классической постановке формулируется как начально-краевая задача Стефана для уравнения теплопроводности, см. [3-5]:
дТ
с— = т) +! (Г), г е п, (1)
где ^ и Г — независимые переменные по времени и пространству, Т — температура, / — функция источника тепла, к — коэффициент теплопроводности, который в анизотропных средах представляется тензором, а с(Т) есть объемная, или удельная теплоемкость.
Здесь расчетная область П предполагается ограниченной, а ее граница Г в общем случае может быть кусочно-гладкой и многосвязной. Для искомой функции — температуры Т(1,Г) — ставятся начальные и краевые условия
ти=0 = Т0(г), гт\т = я(Г, I), (2)
где I есть некоторый оператор граничных условий. Типичная ситуация для краевых условий — это задание температуры То(Г,() на части границы Гд (условия Дирихле), задание теплового потока I (Г,1) = — на другом граничном фрагменте Г N (условия Неймана), а также нелинейное условие
излучения: — к^щ = Т4 — Т4, где п есть внешняя нормаль к П, а Тт — температура внешней среды. Будем считать, что область П состоит из двух подобластей П^ и П5, первая из которых соответствует жидкой фазе среды, а вторая — твердой (ликвидус и солидус соответственно).
На границе раздела фаз, являющейся изотермой Т = Т^, справедливо уравнение баланса энергии (условие Стефана):
£ (дТ/дп)^ к1 (дТ/дпъ = Р5\ъп, (3)
где А — удельная теплота плавления, — плотность среды в твердой фазе, — нормальная к поверхности раздела фаз составляющая скорости движения данной поверхности, а нормаль п выбирается внешней по отношению к твердой фазе. Отметим, что наряду с классической постановкой задачи Стефана (3), может быть сформулирована модель и с переходной температурной зоной Т € [Т^ Д, Т+Д], в пределах которой содержится смесь из жидкой и твердой фазы вещества.
Рассматриваемая задача может быть сформулирована также в энтальпийной постановке:
дН
— = а^^гаа т ) +! (Г), г е п, (4)
о1
где Н есть объемная энтальпия, определяемая формулой
Т
Н (Т )= [ с(Т )йТ + вЬ, (5)
в которой L — это объемная теплота фазового перехода, Tc — произвольно выбранное значение температуры, а величина 0 = 1 в жидкой фазе и 0 = 0 — в твердой. По известной энтальпии температура находится следующим образом:
T (Н) = I
Г ^ < : T = HCHL + Т^ А,
Н < H < Щ : T = [(7> — Д)(H2 — H)+
(6)
+ (Т> + Д)(Н-Н1)]/(Н2- Н1),
Н > Н2 : Т = НСН2 + т + А,
где с = с^ и с = с суть значения теплоемкости в жидкой и твердой фазе соответственно. В данной модели энтальпийный интервал Н € [Н2Н1] фактически определяет переходную двухфазную зону, которой соответствует температурная двухфазная зона Т € [Т> — А, Т> + А]. Рассматриваемая кусочно-линейная зависимость Т(Н) является некоторым обобщением модели [5].
Изучаемые процессы теплопроводности мы будем рассматривать в пористой среде, характеризуемой коэффициентом пористости т € [0, 1], который определяет массовую долю вмещающей породы, или скелета (при этом величина 1 — т характеризует в грунте долю флюида, который может находиться в твердой или жидкой фазе). В этих фазах определяются эффективные, или усредненные коэффициенты теплоемкости и теплопроводности:
СзРз = (1 - т)р5с5 + тр$сС$с, с_т = (1 - т)рс + тр^Сцс,
^ = (1 - т)к + тк,^, ( )
к = (1 — т)к + тк5с,
где р5 и р( — плотности твердой и жидкости фазы соответственно, а символами <«о> обозначаются теплофизические параметры скелета. Отметим, что при учете пористости среды величины с5 и с в (6) надо заменить на с5Ср соответственно.
Относительно всех геометрических и материальных исходных данных задачи предполагаем, что они обеспечивают функциональные свойства искомого решения, необходимые в дальнейшем для обоснования применяемых нами методов аппроксимации.
Перейдем далее к уравнениям Дарси, описывающим течения в пористых средах. Мы приведем в достаточно общей форме трехмерные уравнения многофазной неизотермической (предполагается, что какие-то параметры модели зависят от температуры) модели фильтрации вязкой сжимаемой газожидкостной смеси в анизотропной пористой среде (ниже используются декартовые системы координат х, у, г). Более конкретно, мы описываем постановку для трехфазной среды, актуальной при изучении процессов нефтедобычи, где п/ = 3, а I = 1,2,3 соответствуют воде, нефти и газу. В данной модели «черной (нелетучей) нефти» предполагается, что вода и нефть не смешиваются между собой, а газ растворяется в воде и нефти [6, 7]:
д -к — (^дс) = (Ар/ - 7/АО)) + qc,
с =1,...,Пс; Р2- Р3 = Рс,2,3, Р2 - р1 = Рс,2,1, (8)
51+ 52 + 53 = 1, N1= I1, N2 = В2, N3 = В3+ #3,2 В2,
В1 В2 В3 В2
где п/ и пс — количество фаз и число компонентов в смеси, р/ = р/(¡,х,у,г) — давление в /-й фазе, В/ = В/(р/) — коэффициент объемного расширения фазы /, N = (Д1,...,ДПс) — вектор молярных плотностей, г]с,/ = г]с,/(р1,...,рп/— молярная доля с-го компонента в /-й фазе, £/ = £(р/, N) — молярная плотность фазы /, к = к (р 1,..., п/) — тензор абсолютной проницаемости, /// = р/ (р/) — вязкость фазы /, 5/ = 5/(¡,х,у,г) — насыщенность /-й фазы, qc = qc(р1,...,рп/,N, 1,х,у,г) — источник компоненты с, (р = <р(р1,...,рп,х,у,г) — пористость среды, #3,2 = #3,2(р2) — растворимость газа в нефтяной фазе, кс/ = кс / (51,53) — относительная проницаемость фазы /,7/ = р^ — вертикальный градиент гидростатического давления в фазе /, g — гравитационная постоянная, О = О(х,у,г) — вектор глубины (сверху
вниз), р/ = р/ (р/) — массовая плотность фазы /, Рс,2,з = Рс,2,з(5з) — капиллярное давление в системе нефть-газ, Рс>2,\ = Рс,2,\(5\) — капиллярное давление в системе вода-нефть.
В системе (8) неизвестными являются функции Ыс, 5/ и р/, а в самом распространенном случае количество компонент пс равно трем. Для рассматриваемой проблемы задаются начальные значения искомого решения, а на внешней границе моделируемого резервуара ставятся граничные условия постоянного давления (условия Дирихле) или заданного потока (условия Неймана).
Описанная модель представляет нетривиальную математическую проблему, отягчаемую, как правило, сложной геометрией расчетной области (реальный пример — куст криволинейных скважин в сложной геологической структуре) с сильной разномасштабностью и контрастными свойствами материальных сред, которые на практике известны очень приближенно и по-хорошему требуют еще решения обратной задачи по идентификации параметров модели на основе сравнения расчетных данных с геофизическими измерениями.
Для рассмотренной задачи в реальных геофизических условиях решение уравнений требуется проводить с учетом зависимости всех коэффициентов от температуры среды, являющейся переменной во времени и пространстве. Данный вопрос чрезвычайно актуален, например, при эксплуатации газовых и/или нефтяных скважин в районах Крайнего Севера с многолетней криолитовой зоной, в которой циклически проходят процессы промерзания и растепления слоисто-неоднородных массивов пород.
Оптимизационная постановка обратных задач. Рассмотренные выше модели относятся к постановкам прямых задач, которые при решении обратных задач методами оптимизации формально представляются как параметризованное операторное уравнение состояния, описываемые начально-краевыми задачами вида
Ьи(р) = Т{Ы), £е й, 0 < £ < Т < ос,
(9)
/и = §(х,1), х € Г, и(х,0) = и°(х).
Здесь р = {рк} — т-мерный вектор искомых параметров задачи, Ь является, например, дифференциальным оператором
Ь = А^- + ДВД + СД + О, (10)
о1
с матричными коэффициентами А, В, С, О, а / — некоторый оператор граничных условий, включающий в общем случае краевые условия 1-го и/или 2-го и 3-го рода
и = gо, хе То; Дыи + Ам/\пи = gN, (11)
на участках границы Г о и Г N соответственно.
Оптимизационные методы основываются на минимизации определяемого в описании обратной задачи целевого функционала:
Ф0(и(х,1,рар£)) = ттФ0(и(х4,р)), (12)
при заданных т\ линейных и т2 нелинейных ограничениях на параметры задачи:
ртт<рк<ртах, £ = \,...,т\,
\ )
Ф/(и(х,£,р)) < д/, / = \,...,т2, т\ + т2 = т.
Здесь — некоторые функционалы, определяемые какими-либо дополнительными условиями на исходные данные рассматриваемой задачи условной минимизации. Отметим, что в последние десятилетия методы оптимизации претерпели бурное развитие, приведшее к появлению новых эффективных подходов: модифицированные методы множителей Лагранжа, алгоритмы внутренних точек, последовательного квадратичного программирования, доверительных интервалов и т. д. [8]. Отметим, что данные задачи могут иметь не один, а несколько минимумов, и тогда возникает более трудная проблема глобальной минимизации. Заметим, что практически все обратные задачи, за исключением специальных случаев, реализуются итерационными методами, каждый шаг которых включает решение прямой задачи.
Наукоемкие алгоритмы, технологии и структуры данных
Под понятием «наукоемкие» алгоритмы подразумеваются два ключевых аспекта. Первый круг вопросов относится к многомасштабным (в том числе иммерсионным) разрывным метода Галеркина (РМГ) высоких порядков, которые занимают промежуточное положение между методами конечных объемов и конечных элементов (МКО и МКЭ), являясь в определенном смысле обобщением алгоритмов суперэлементов Р. П. Федоренко [9, 10]. Данный подход заключается в том, что для каждого конечного объема, или элемента, строятся «свои» базисные функции различных порядков, которые определяются аналитически или численно и формируют пространство приближенных решений. Для обеспечения однозначности на внутренних границах смежных элементов ставятся дополнительные условия для потоков субстанции (массы, импульса и т. д.). Данные алгоритмы обладают двумя замечательными свойствами: они в работах различных авторов получили теоретическое обоснование для реальных многомерных начально-краевых задач в расчетных областях с многосвязными кусочно-гладкими граничными поверхностями и контрастными материальными характеристиками сред, а также имеют универсальные хорошо распараллеливаемые технологии реализации на неструктурированных сетках, в том числе иммерсионных. Последнее означает, что в конечном элементе могут иметься детали более мелкого масштаба, а получаемые «подсеточные эффекты» учитываются сеткой второго уровня внутри этого элемента. Данное перспективное направление активно развивается в последние десятилетия, и мы приводим ссылки только на характерные публикации зарубежных и отечественных авторов [11-14].
Второй аспект рассматриваемых вычислительных методов и технологий повышенной сложности касается решения СЛАУ сверхвысоких порядков (до 1010 и выше) с разреженными плохо обусловленными матрицами, возникающими из МКО- или МКЭ-аппроксимаций исходных задач на адаптивных квазиструктурированных сетках. Основным орудием здесь являются блочные предобусловлен-ные итерационные алгоритмы в подпространствах Крылова, среди которых мы отдаем предпочтение мульти-предобусловленным методам полусопряженных или сопряженных невязок (для несимметричных или симметричных матриц соответственно). Главная решаемая здесь проблема заключается в обеспечении высокой скорости итераций по подобластям, когда их количество становится очень большим, что необходимо для массивного параллелизма. При этом приходится привлекать разнообразные приемы ускорения вычислительного процесса, что неизбежно приводит к логическому усложнению алгоритмов.
Разрывные методы Галеркина. Методы аппроксимации высокого порядка для многомерных краевых задач мы для краткости опишем на более простой, в сравнении с (8), квазилинейной модели диффузии с конвекцией:
-V • к(г, и, Vи) + к(г, и,Уи) = >(г), ге О, и = gО, геГО, ()
Уи-п = gN, ге ,
где П = П и Г, Г = Г о и и к (г,£,77) = (г,£,77); I = 1,2,3}. В уравнении (14) мы можем формально предполагать, что в нестационарном случае член с производной по времени перенесен в правую часть, и при этом после пространственной аппроксимации уравнения получаемая система обыкновенных дифференциальных уравнений будет решаться с помощью известных алгоритмов дискретизироваться по времени.
Уравнение (14) перепишем предварительно в смешанной постановке как систему уравнений первого порядка:
—V • 4 + q0 = >, д = к(г, и,о), q0 = к(r, и, о), о = Уи. (15)
Предполагая для простоты изложения область П полигональной, а также полагая go = gN, Ь = 0, обозначим через регулярное семейство ее разбиений на тетраэдры (конечные элементы, или объемы) К с характерными диаметрами Н > 0 и с граничными поверхностями = 1,...,^, где Nt — количество тетраэдров. Обозначим далее через п поле внешних единичных нормалей -го тетраэдра, через Рт(КН) — множество определенных на КН полиномов степени т по совокупности независимых переменных, а через РН — множество «внутренних» граней тетраэдров, т. е. не лежащих на внешней сеточной границе ГН = ГО и . Пространство приближенных решений иН при заданном т > 1
определим как
Ук = ^(О1) : у\кн€ ?т(К\), ^ = 0, П1}. (16)
Аналогично введем пространство для аппроксимации до и векторных компонент д,о:
= Ьж(П11) : Рт{К), К1 е П11}. (17)
Очевидно, что функции из этих пространств не обязаны быть непрерывными (на множестве F'h они являются двузначными) и что Vh С Wh, Ун^н С Wf, где оператор Vh определяется через V поэлементно. Для однозначного определения предельных (односторонних) значений функций, терпящих разрыв на Fh, введем векторное кусочно-постоянное поле p единичных нормалей на F^U Г, которое назовем потоком [13]. Поле p определяется произвольным образом, но с единственным условием: для любого элемента Kh должна существовать выходящая из него нормаль nf, т.е. (p ■ nf ^h = 0. Здесь
и далее через v± = max{±v, 0} обозначается положительная и отрицательная часть функций v, так что v = v+ — v_. Определим также такие понятия, как одностороннее значение и скачок функции w в точке r G F h вдоль нормали p(r):
w±p(r) = lim, [w]p(r = w+p)(r)- w-p(r),
причем на внешней границе Qh скачок полагаем нулевым, т. е. для r € Th [w]p(r)|™ = 0.
h
Для единственного определения значения (следа) функции Uh € Vh на общей грани F двух контактирующих тетраэдров Kh и Kh принимается следующий противопотоковый способ:
Uh\F = Uh-p = (p-Пк)+ Uh|K + (p-ПкГUh|K . (18)
Для векторных функций qh G Wf, наоборот, принимается значение по потоку:
qh|F = qh+p = (p-Пк )+qh\K, + (p-Пк Г qh|K. (19)
Под разрывным смешанным методом Галеркина для приближенного решения задачи (14) на сетке Ùh = Qh U Гн понимается следующая дискретная задача: найти функцию Uh € Vh, а также вектор-функции qh,Oh^ Wf такие, что для выбранного поля p выполняется равенство
ah(Uh,Vh) = fh(vh)V Vh£ Vh, (20)
где билинейная и линейная формы определяются следующим образом:
ah(Uh,Vh) = q0,hVhdr + ^ / UhV ■ Whdr + qh+pp [vh]pds,
i K^hJK i (21)
fh(Vh) = fvhdr,
h
а функции он = Sh(Uh), qh = Qh(un) и q0,h = q0,h(UN) определяются по Uh из соотношений
Oh ■ Whdr + J2 JuhV- Whdr + jUh,- p [Wh] ppds = 0, Vw^ Wjhf,
tth K&lh K rh
[(qh- k(r,Uh,Oh))- Wh + к(Г,щ,о h)w0,h] dr = 0, (22)
h
VW^ Wjf, W0,h£ Wh.
Уравнения (20)-(22) определяют систему алгебраических (в общем случае нелинейных) уравнений для определения искомых приближенных решений Uh, qh и oh Мы не будем вдаваться в детали алгоритмов их решения, а также нетривиального обоснования данного варианта РМГ [13]. Отметим
только, что если выполнены достаточные условия монотонности операторов, а также существования и единственности обобщенного решения задачи (1) в пространстве Соболева Hq(îî), и если к тому же u G Hm+1(Q), q G [Hm(Q)]3, то погрешности построенных приближенных решений в соответствующих нормах суть величины порядка O(hm).
Отметим, что в обширной литературе опубликованы различные варианты разрывных методов Галеркина с многочисленными примерами их успешного применения, подтверждающими перспективность данного направления.
Решение больших разреженных СЛАУ. Пусть требуется решить некоторым итерационным методом невырожденную вещественную несимметричную систему линейных алгебраических уравнений
Au = f, A е ПыN, u,f e TlN, (23)
при начальном приближении u0, определяющем вектор начальной невязки r0 = f — Au0.
Мы рассмотрим возможность применения какого-то из блочных или мульти-предобусловленных, методов полусопряженных направлений (MPSCD, от Multi-Preconditioned Semi-Conjugate Direction [15]), являющихся развитием алгоритма обобщенных сопряженных невязок GCR (Generalized Conjugate Residual). Предлагаемое семейство методов запишем в следующем виде:
um+l = um + Pm âm, rm+1 = rm^APm «m, m =0,1..........(24)
где Pm = (pm...pMm) e ,Mm есть составленная из направляющих векторов pf матрица, а âm = (am,i-..am,Mm)T — вектор итерационных параметров, которые определяются из ортогональных свойств
(Apk A1pm ) =Ж, = WAtf), (25)
7 = 0,1,; m = 0,1,...,m — 1; k,k' = 1,2,...,Mm.
Здесь — символ Кронекера, равный единице при m = m', k = k' и нулю — в противных случаях, а значения 7 = 0,1 определяют в дальнейшем метод полусопряженных градиентов или полусопряженных невязок (MPSCG или MPSCR соответственно). Отметим, что в отличие от обычных методов полусопряженных направлений, в формулах (24) на каждой m-й итерации используется не один, а Mm направляющих векторов, количество которых на разных итерациях может, вообще говоря, меняться.
Векторы коэффициентов âm ={«„/} в (24) для 7 = 0,1 при выполнении (25) из условия экстремума
дФ^/да^ = 0, Ф-V+1)- (rm+1,^ 1rm+1), (26)
определяются формулами
'V/ = (^ "m,/1 ,' )И>m,/
«2 = AB-^rn/m. (27)
Направляющие векторы рт определяются по условиям ортогональности (25) при 7 = 0,1 в следующем блочном виде:
m
P0 = {p0 = B-jr°}, Pm+1 = Pm+1,0 - £ Pk= {РГ+1 =
k=0
m
= B-1 rm+1 M R(l) nk\ (28)
= Bm+1,/r m,k,/p/ ,
k=o/=1
m = 0,1,...; Bm,/£1IN,N, / =1,...,Mm, 7 = 0,1.
Здесь = { pmxi} = )Г e ПМт суть векторы коэффициентов, а KN N -
предобуславливающие матрицы, которые выбираются из соображений невырожденности, легкой обратимости и эффективного ускорения конструируемых итерационных процессов. Отметим, что рассматриваемые «предобуславливатели» Bm,i являются динамическими, или гибкими (flexible, как в методах FGMRES), поскольку они для каждого I зависят от номера итерации m.
Подставляя выражения (28) в условия ортогональности (25), мы получаем формулы для коэффициентов
= (Api, AB-ml+l/m+l)/p{2) = {Api Apmk)lpti, m = 0,1,...-,
к = 0,...,m; l = \,...,Mm,
mi
в которых векторы pt определяются из соотношении
k-1
т,к = nm,k-1 ф) пк-1 = в-1 rm+1 V^ fí^) pm l =1 M
pl = pl ~ Pm,k- 1,lpl = Bm+1,lr - 2^Pm,i,lpl , l = 1,...,Mn
i=0
(30)
к = 0,1,...,m +1; pf,0= Bmlurm+1, pf,m+1= p^+1
на основе модифицированного метода Грама — Шмидта.
Можно показать, что функционалы Ф^ ^ (гт+1) при условиях экстремальности для любых q = 0,1,...,т справедливы соотношения
т Мт
ф(т\гт+1) = Б~1 )2/$. (31)
k=q 1=1
Для 7 = 1 или при симметричности матрицы А в методе МР8СЭ обеспечивается минимальность нормы | |гт+^| = (гт+1 ,гт+1)1/2 в «мультипредобусловленном» подпространстве Крылова
^т+У,А) = врал{ Бо}г0,..,БоМмо г0,
г1,...мтМ,....*тБт\ гт,...,АтБтМтгт},
размерность которого равна ^т+1 = Мо + ... + Мт. В методах полусопряженных градиентов, т. е. при 7 = 0, функционал Фт0^(гт+1) = (А1 гт+1 ,гт+1) при несимметричности А своего минимума в подпространстве Крылова, вообще говоря, не достигает.
Для симметричных матриц методы МР8СЭ существенно упрощаются, так как коэффициенты при к > 0 обращаются в нуль, и формулы (28)-(30) в этой ситуации дают алгоритмы сопряженных направлений (МРСв или МРСЯ). Однако в общем случае рассмотренные блочные методы полусопряженных направлений используют так называемые «длинные» рекурсии, т. е. необходимо хранить все вычисленные на предыдущих итерациях направляющие векторы, или составленные из них матрицы Р1 ,...,Рт. Во избежание слишком больших требований к памяти при этом мы будем использовать известную процедуру рестартов, заключающуюся в том, что через какое-то количество итераций вектор невязки определяется не из обычных рекурсий (24), а из исходного уравнения, и крыловский процесс как бы формируется заново, но с текущего приближения для искомого вектора.
Вид предобуславливающих матриц Бт, мы пока не конкретизируем, и они будут рассмотрены ниже при описании параллельных методов декомпозиции областей.
Стратегии и тактики распараллеливания технологических стадий моделирования
Проблема автоматизации распараллеливания алгоритмов, несмотря на многократные попытки решения, остается до сих пор открытой. Отдельные ее фрагменты решаются на уровне технических средств или опций компилятора, но в целом производительность кода зависит от профессиональных знаний и искусства программиста. Имитационное моделирование гетерогенной МВС практически невозможно, поскольку объем и разнообразие функциональных элементов суперкомпьютера потребовали бы для этого еще больших мощностей. Более того, для уже реализованного алгоритма оценить его быстродействие с достаточной точностью тоже невозможно, ввиду отсутствия адекватной модели машинных вычислений. Практически единственно на какую помощь может рассчитывать прикладной программист-вычислитель, — это динамическая картина интенсивности загрузки различных элементов МВС, формируемая с помощью библиотеки параллельных инструментариев.
Общие характеристики параллелизма. В силу существующих методологий мы можем пользоваться только такими «укрупненными» характеристиками, как ускорение вычислений Бр и эффективность распараллеливания Ер-.
Бр = Т1/Тр, Ер = Бр/Р, (33)
где Tp означает время решения задачи (или реализации какого-то алгоритма) на P процессорах. Понятно, что идеальным вариантом является Sp = P, Ep = 1, хотя на практике приходится мириться с ситуацией, когда эти величины в несколько раз меньше. Отметим, что если отношение объема последовательных вычислений к объему параллельных операций равна , то ускорение от распараллеливания выражается законом Амдала [17]:
Sp = Ti/(0Ti + (1 - 0Ti))/P = P/[1 + 0(P - 1)]. (34)
Отсюда следует, в частности, что при P оо величина Sp остается равной примерной 1/0. Общее календарное время выполнения складывается из времен реализации арифметических действий Ta и межпроцессорных обменов Tc:
T = Ta + Tc, Ta = NaTa, Tc = M(Tq + NCTC), (35)
где a есть среднее время выполнения одной арифметической операции, c - время передачи одного числа, M и tq общее количество коммуникаций и время ожидания (настройки) каждой из них, а Nc — средний объем каждой из их. Поскольку на практике выполняются соотношения Ta ■< tc ■< tq, то отсюда следует, что надо стремиться сокращать не только количество передаваемой информации, но и число самих передач.
Понятно, что такие «средние по больнице» характеристики слишком грубые, но более детальные описания получить трудно, и сама наука распараллеливания фактически является экспериментальным исследованием, базирующимся на многократных расчетах и измерениях, а оптимизация кода заключается в сравнительном анализе производительности различных технологических приемов при разных конфигурациях МВС и опциях компилятора. В целом архитектуру современного суперкомпьютера можно представить кластерной конфигурацией, т. е. сетью из вычислительных узлов, имеющих свою распределенную память и объединяемых по принципу «каждый с каждым» с помощью простой системы передачи данных MPI (Message Passage Interface). Каждый типовой узел содержит несколько центральных многоядерных (до нескольких десятков ядер) процессоров CPU (Central Processor Unit) с общей многоуровневой памятью, а также графические ускорители типа GPGPU (General Purpose Graphic Processor Unit) или Intel Phi с сотнями ядер и достаточно сложными внутренними структурами, управляемыми своими специальными программными системами. На каждом CPU «внутреннее» распараллеливание организуется средствами многопотоковых вычислений (multi-thread, системы типа OpenMP). Кроме того, на сверхбыстрых регистрах, количество которых относительно невелико, с помощью их системы команд типа AVX можно формировать векторизацию вычислений с детальным учетом особенностей доступа к памяти нижнего уровня (например, кэш-1).
Очевидно, что программное обеспечение распараллеливания алгоритмов представляет собой достаточно дорогой овеществленный интеллектуальный потенциал. Поэтому оно должно разрабатываться как кросс-платформенный продукт, адаптируемый к непрерывно происходящей эволюции компьютерных архитектур, а при его создании очень важно придерживаться следующего правила: отсутствие программных ограничений на число степеней свободы задачи и на количество используемых вычислительных процессоров или ядер.
Технологические стадии распараллеливания. В общем случае решение большой содержательной (междисциплинарной, обратной, со сложными физическими и геометрическими данными) задачи включает различные стратегии и тактики распараллеливания на последовательно выполняемых стадиях моделирования.
Геометрическое и функциональное моделирование — это операции с различными математическими объектами, составляющими описание модели изучаемого процесса или явления, а также вычислительной схемы для реализации компьютерного сеанса (что и как считать, содержание и форма представляемых результатов). Сюда включается описание расчетной области, свойств материальных сред, решаемых функциональных уравнений, начально-краевых условий и т. д. [18]. Данный вычислительный этап не является критическим по ресурсоемкости, поскольку количество геометрических и функциональных объектов относительно невелико (например, меньше 104), поэтому здесь распараллеливание целесообразно реализовывать путем копирования данных по всем организуемым MPI-процессам, или соответствующим узлам, которые, в частности, содержат графические ускорители, с
помощью которых можно быстро осуществлять графический интеллектуальный интерфейс с конечным пользователем. Получаемые в результате этой стадии геометрическая и функциональные структуры данных (ГСД и ФСД) также следует идентично распределять по всем MPI-процессам во избежание лишних коммуникаций. Отметим, что если решается обратная задача путем оптимизации параметризованной исходной информации, то данный этап выполняется многократно, но последовательно (при решении каждой прямой задачи).
Дискретизация задачи. Формально говоря, на данном этапе происходит преобразование ГСД + ФСД => ССД (сеточная структура данных). Другими словами, происходит формирование сеточной расчетной области, включающей всю необходимую информацию о решаемой задаче, которая может содержать миллиарды чисел степеней свободы и, соответственно, узлов сетки. Поскольку основной подход к распараллеливанию для таких случаев состоит в декомпозиции всей области на подобласти, каждой из которых сопоставляется свой MPI-процесс, то формирование ССД целесообразно сразу делать распределенным образом. Конечно, напрашивается более простой логических путь — сначала построить ССД на одном процессоре, а потом ее «раскидать» по всем остальным, но в этом случае памяти одного компьютерного узла может не хватить для всей большой сетки.
Что касается алгоритмов построения адаптивных сеток с высокими критериями качества (в простейшем случае это отсутствие почти вырождающихся конечных элементов и сильной неравномерности шагов), то здесь имеется большое разнообразие подходов [19], и мы на этих вопросах останавливаться не будем. Можно сказать, что построение оптимальной сетки в общем случае — это проблема более сложная, чем решить текущую задачу моделирования на просто хорошей, или даже удовлетворительной дискретизации. Поэтому здесь в качестве «золотой середины» используют достаточно простые и быстрые сеточные генераторы, так что в смысле времени счета данная стадия также не является самой критической. Исключение составляют динамические сетки, многократно перестраиваемые в процессе решения задачи, но таких случаев стараются избегать, и для этого существуют различные подходы.
Важно отметить, что современные и наиболее перспективные методы решения многомерных начально-краевых задач используют все усложняющуюся логику вычислительных схем с использованием многосеточных подходов и декомпозиции областей. Естественно, что такие алгоритмические возможности надо предусматривать загодя, т. е. как раз на этапе построения распределенных по процессорам сеточных подобластей.
Аппроксимация уравнений. Основные подходы к построению дискретных сеточных моделей — это классические методы конечных разностей, конечных объемов и конечных элементов (МКР, МКО и МКЭ соответственно), на стыке которых появились и активно развиваются в последние десятилетия разрывные методы Галеркина (РМГ), обладающие уникальной универсальностью: для них существуют логически простые, естественным образом распараллеливаемые и теоретически обоснованные схемы высокого порядка точности, в том числе для многомерных задач с геометрически сложными кусочно-гладкими граничными поверхностями и контрастными материальными свойствами различных сред. Поэлементные технологии МКО, МКЭ и РМГ основываются на единообразном вычислении локальных матриц для различных базисных функций и операторов с последующей сборкой глобальной матрицы. Естественно, что последняя в параллельной версии должна формироваться сразу распределенной по различным MPI-процессам. Предложена концепция построения библиотеки аппроксиматоров CHEBYSHEV, использующая систематизацию параллельных технологий для разных типов сеточных элементов, базисных функций и дифференциальных операторов [20]. Архитектура этой библиотеки организуется как система, открытая для расширения состава реализуемых моделей и алгоритмов, в том числе разными группами разработчиков.
Данная стадия моделирования является, пожалуй, самой наукоемкой, поскольку в ней реализуются результаты фундаментальных исследований высоких аппроксимационных теорий. Математическая техника для базисных функций высоких порядков и нетривиальных конечных элементов действительно сложна, а запись формул для матричных элементов могут насчитывать много страниц. Исторически именно по этой причине, главным образом, в мире получили наибольшее распространение простейшие алгоритмы низких порядков. В данном случае кардинально изменить ситуацию может автоматизация построения алгоритмов на основе машинной реализации аналитических выкладок с помощью, например, программных систем компьютерной алгебры MATHEMATICA и MAPLE. Примером
такой реализации может служить пакет прикладных программ (ППП) HELMHOLTZ 3D [21].
Результатом аппроксимационных процедур является алгебраическая структура данных (АСД), содержащая полную информацию о формируемой дискретной системе уравнений. Поскольку в дальнейшем СЛАУ должна решаться параллельными методами декомпозиции, то АСД необходимо сразу создавать распределенным образом. Типичные матрицы являются сильно разреженными (например, ненулевые элементы составляют доли процента от общего числа), поэтому их хранение осуществляется в компактных форматах (один из наиболее распространенных — это Compressed Sparse Row [22]). Отметим, что в общем случае для квазиструктурированных сеток в различных сеточных подобластях могут использоваться разные матричные форматы.
Решение алгебраических задач. Здесь две главные вычислительные проблемы — это решение СЛАУ и нахождение собственных чисел и соответствующих собственных векторов матриц, из которых мы остановимся только на первой как наиболее распространенной. Рассматривая для общности несимметричные системы, мы остановимся на особенностях параллельной реализации двухуровневых итерационных процессов декомпозиции областей в подпространствах Крылова при использовании внешних итераций с помощью мультипредобусловленных методов полусопряженных невязок, описываемых формулами (24)-(32) и реализованными в составе библиотеки программ KRYLOV [23].
Расчетную область представим в виде совокупности из P непересекающихся подобластей:
P
fi = fiUr=Ufiq; =0 при q = r. (36)
q=i
Общую границу двух контактирующих подобластей обозначим через Yqr, а всю границу Tq подобласти Пq разобьем на две части: внутреннюю rq и внешнюю rq = Г |J Tq (принципиальное отличие состоит в том, что на Гq ставится краевое условие для искомого решения, а на Гq — нет):
Tq,r = ^qfl ^r> rq = U Tq,r = ГqUrq' (37)
r
Граница расчетной области в целом составляется из внешних граничных сегментов подобла-
(38)
стей:
r = U rq = U Tq,o = r ^U1^.
Здесь Г д и Г N означают совокупности граничных сегментов с краевыми условиями Дирихле и с краевыми условиями, содержащими производные по нормали (т. е. 2-го или 3-го рода), соответственно. Кроме того, мы отмечаем, что «внешнее дополнение» замыкания расчетной области П = П У Г ко всему трехмерному пространству может рассматриваться формально как подобласть с нулевым номером (По = К3/й). Отметим, что до сих пор мы рассматривали ограниченные расчетные области, однако методы декомпозиции можно применять и для решения внешних краевых задач, вводя в рассмотрение неограниченную подобласть (с условиями на бесконечности) как внешность некоторой сферы достаточно большого радиуса, где решение ищется с помощью интегральных представлений [24].
Решение дискретизованной (алгебраической) задачи производится в сеточной расчетной области, состоящей из сеточных подобластей, которые могут рассматриваться как полученные в результате независимых генераций сеток соответствующих непрерывных подобластей. В частности, эти сетки могут быть согласованные или несогласованные, т. е. узлы на внутренних границах контактирующих сеточных подобластей могут совпадать между собой или различаться.
При всей алгоритмической сложности и многообразии современных алгоритмов декомпозиции областей основой параллельного решения рассматриваемых СЛАУ (в алгебраических терминах) является блочный итерационный метод Якоби, или другими словами (в геометрическом смысле) — аддитивный метод Шварца. Если обозначить через й,, ¡,, , = 1,...,Р + 1, подвекторы, соответствующие подобласти П,, или соответствующей блочной строке матрицы А, то такой алгоритм записывается в виде
Ад,дйд = 1д ~ Ад,гйП 1 = 1я, Я = 1, . ,Р; п = 1,2,..., (39)
и,
где п есть номер итерации, а каждый диагональный блок Ад,д — фактически матрица автономной вспомогательной СЛАУ для подобласти Пд. Внедиагональные же блоки Ад,г, г = д, отвечают за взаимосвязи между подобластью Од и контактирующими с ней подобластями, имеющими номера г из множества Шд. Отметим, что введенные в (39) блоки Ад,г в совокупности определяют блочную строку матрицы:
Ад = (... Ад,г ..., ге Шд), А = {Ад, д = 1,...,Р} .
Итерационный процесс (39) легко формально обобщить, добавив к левой и правой частям формулы члены с некоторым одинаковым матричным множителем:
Вд,дип = (Ад,д + Сд,д)ип = ^ д + Сд,дип ~ ^ ^ Ад,гип , п = 1,2,... (40)
г£ Шд
Очевидно, что если при п —> оо итерационные процессы (39), (40) сходятся, то предельный вектор и = {и^ у них один и тот же, а именно — решение исходной СЛАУ. Вводимые дополнительные блоки Сд,д не меняют структуру и логику итерационного алгоритма, если они не изменяют портрет матрицы Вд,д, в сравнении с Ад,д. Обычно добавляются члены только к сеточным уравнениям для околограничных узлов из Од, и эта процедура соответствует изменению интерфейсных, или краевых, итерируемых условий между контактирующими подобластями. В частности, если положить Сд,д = 0йд,д, где Оддд — диагональная матрица, каждый элемент которой равен сумме внедиагональных элементов матриц Ад,г, г € Шд, то при в = 0 интерфейс между смежными подобластями соответствует условию Дирихле, а при = 1 — условию Неймана.
Подчеркнем, что мы для простоты рассматриваем до сих пор декомпозицию без пересечения подобластей, чтобы выделить главные алгебро-геометрические особенности алгоритмов. При наличии пересечений структура методов несколько усложняется. Например, блочное представление СЛАУ будет содержать перекрывающиеся блоки матрицы А и подвекторов. Само построение пересечений целесообразно делать последовательно: сначала осуществляется конструирование подобластей Пд без пересечения, а затем каждая из подобластей постепенно расширяется путем включения в нее по одному сеточному слою на очередном этапе. Величину пересечений естественно измерять параметром Д — числом расширяющих (дополнительных) сеточных слоев подобласти Пд, который для простоты будем считать одинаковым для всех д. Последовательность формируемых границ можно описать следующим образом:
Гд= Г0= {1'€Ш1, ¡€Пд , Г /Пд , = ^ = Од У ,
1*-1 (V = о 1 = о
йд , Ъйд Ъйд 1 £, д
гд = , ипд-1, 1'епд-1, од = п д-д- 'и гд- ч.
Характеризуемая формулами (39), (40) блочно-диагональная матрица В = {Вд,д} определяет параллельный метод декомпозиции областей в подпространствах Крылова (24)-(32), которые используются, естественно, вместо слишком медленного стационарного итерационного процесса Якоби -Шварца (39), (40). Обращение матрицы В — это фактически осуществляемое синхронно на разных процессорах решение подсистем в соответствующих подобластях. Отметим, что если эти подсистемы решаются итерационным методом, то фактически на каждой т-й итерации мы имеем разные предобу-славливатели Вт.
Для ускорения сходимости МДО эффективно применяются методы грубосеточной коррекции, имеющие в разных версиях названия алгоритмов дефляции, агрегации и основанные алгебраически на малоранговой аппроксимации исходной матрицы А [15]. Дело заключается в том, что при увеличении числа подобластей Р, необходимого для задействования большего количества процессоров, скорость сходимости итераций типа блочного алгоритма Якоби существенно замедляется в силу локальности связей между подобластями. Суть грубосеточной коррекции заключается в том, что она экономичным образом дополнительно устанавливает на каждой итерации дальние связи между подобластями.
На качественном уровне данный подход заключается в построении грубой сетки с числом узлов Ыс <С N. Формально эта редкая сетка может быть никак не связана с макросеткой, составленной из подобластей, однако наиболее естественным является выбор Ыс = Р, когда каждой подобласти сопоставляется один узел. На построенной редкой сетке определяется некоторая система базисных финитных функций <рд(г), д = 1,...Р, с помощью которой формируется (а затем решается) СЛАУ
порядка N относительно достаточно грубого галеркинского приближения к искомому вектору. Проще всего взять базисные функции нулевого порядка, т. е. кусочно-постоянные: , (г) = 1 при г € О, и <Ря (г) = 0 в противном случае. Тогда с их помощью определяется прямоугольная матрица полного ранга
Г = (ш>1...тМс )е Пы л, тЯ = {ш!Я = <рч (П), ¡е ПЛ},
компоненты вектор-столбцов которой — это значения базисных функций в узлах исходной сетки. Опуская детали (см. подробнее [16]), алгоритм можно представить с помощью «грубосеточной» предобу-славливающей матрицы
Б~1 = , А = КМсл,
простейшее применение которой заключается в построении «хорошего» начального приближения. В любом современном итерационном процессе в подпространствах Крылова вычисления начинаются с построения начального вектора г° и направляющего вектора р0. Пусть й~1 будет произвольный вектор. Тогда до начала итераций определяем величины
й0 = й-Ч Б-хт~ \ т~1= / — Ай г0= ¡^Ай0, р0 = г0^ Б- 1г0,
(41)
что обеспечивает выполнение свойств ортогонализации
^тг(0 = 0, WTAp0 = 0, (42)
по отношению к «пространству дефляции» Брап-^1 ,...,'Мс}. Отметим, что формулы (41), (42) реализуются сравнительно экономично, поскольку требуют только решения вспомогательной «агрегированной» СЛАУ с матрицей А малого порядка, которая в алгебраическом смысле представляет собой малоранговую аппроксимацию исходной матрицы А. Требуемое решение можно получить с помощью прямого или некоторого итерационного метода. Рассмотренные «грубосеточные» предобуславливатели Бс могут использоваться и на каждом шаге крыловского итерационного процесса, составляя в совокупности с блочно-диагональными матрицами {Б,,,} мульти-предобусловленные алгоритмы, описанные в предыдущем пункте.
Общая схема параллельной реализации задачи. В целом вычислительный процесс решения междисциплинарных прямых и обратных задач имеет разветвленную структуру с многократными вложенными циклами, что предоставляет богатые возможности для выбора и оптимизации технологических режимов оптимизации расчетного сеанса. Если рассматривать этапы «снизу — вверх», то процедуры моделирования можно охарактеризовать следующим образом.
Если проблема является междисциплинарной и описывается системой дифференциальных уравнений со многими неизвестными функциями (например, больше пяти), то получаемые после дискретизации СЛАУ имеют крупно-блочный вид со специальной структурой, определяемой конкретным видом взаимодействия моделируемых процессов. Здесь общие рекомендации по алгоритмам дать трудно, кроме того очевидного факта, что для больших задач, так или иначе, применяться должны итерационные аддитивные или мультипликативные методы в подпространствах Крылова.
Для систем нелинейных алгебраических уравнений (СНАУ) итерационный процесс получается двухуровневым и включает какие-либо варианты квазилинеаризации типа метода Ньютона с многократным последовательным решением СЛАУ, при этом зачастую требуется пересчитывать и саму матрицу. В случае решения динамических задач с помощью неявных схем мы получаем еще один внешний цикл — численное интегрирование уравнений по времени, на каждом шаге которого могут решаться СНАУ и выполняться остальные необходимые процедуры (вплоть до пересчета сетки, особенно если при этом сильно меняется структура решения и конфигурация границы). Наконец, если решается обратная задача с помощью методов оптимизации, то формируется очередной внешний цикл с последовательным решением прямых задач и вычислением целевого функционала при разных значениях искомых параметров в исходных данных. Естественно, что на всех вычислительных стадиях подразумевается постобработка и визуализация для контроля пользователем процесса моделирования. Последние операции могут быть очень ресурсоемкими, и здесь свою роль должны сыграть быстрые графические ускорители с большим количеством ядер. А поскольку при этом критическими являются такие качества, как оперативность, наглядность и многоплановость представления информации, то в
данных случаях также эффективно может применяться параллельная обработка данных по подобластям.
Заключение
Представленные нами (далеко не в полном объеме) междисциплинарные прямые и обратные задачи нефтегазовой отрасли являются чрезвычайно актуальными и экстремально ресурсоемкими, а высокопроизводительные вычислительные методы и технологии, обеспечивающие их эффективное и оперативное решение на современных суперкомпьютерах постпетафлопсного уровня, составляют канву крупномасштабного математического эксперимента, проведение которого невозможно без создания интегрированного программного окружения, встраиваемого как неотъемлемый инструмент производственного процесса. Такое математическое и программное обеспечение является очень крупным по объему, с неизбежными трудозатратами во много десятков человеко-лет и должно удовлетворять строгим методическим и техническим нормативам. Такой комплекс должен иметь гибкие средства расширения состава модулей и алгоритмов, чтобы идти в ногу с быстро развивающейся наукой, а также адаптации к эволюции компьютерных архитектур и платформ, для поддержки длительного жизненного цикла продукта, без которой не может быть экономической оправданности всего проекта. Необходимо отметить и такие технические требования к разработке, как возможности открытого переиспользования внешних программных продуктов, в которых овеществлен огромный интеллектуальный потенциал (САПР-овские системы, графические ППП и т. д.), а также согласованного участия различных групп разработчиков. Наконец, при превышении какого-то критического уровня общего объема задействованных программ существование и развитие проекта становится невозможным без применения компонентных технологий и средств управления большими разработками, которые активно формируются в мировом вычислительном сообществе [25].
Важно также понимать, что при всей насыщенности современными математическими достижениями конечный программный код выходит за рамки академических исследований и должен иметь сугубо производственный характер со своими строгими правилами эксплуатации, сопровождения и развития. А для его востребованности и окупаемости необходима большая научно-организационная работа по созданию специальной инфраструктуры и концепции продвижения продукта среди широкого круга потенциальных пользователей, которых надо просвещать, обучать и фактически создавать из них специалистов в профессии нового поколения.
Работа поддержана грантом РФФИ № 16-29-15122, а также грантом РНФ № 15-11-10024.
ЛИТЕРАТУРА
1. Ильин В. П. Фундаментальные вопросы математического моделирования // Вестн. Рос. академии наук. 2016. Т. 86. № 4. С. 26-36.
2. Ильин В. П. О вычислительных методах и технологиях решения задач нефтегазовой отрасли // Вестн. кибернетики. 2016. № 2. 2016. С. 117-127.
3. Вабищевич П. Н., Самарский А. А. Вычислительная теплопередача. М. : Едиториал УРСС. 2003. 784 с.
4. Бычин И. В., Гавриленко Т. В., Галкин В. А., Гореликов А. В., Ряховский А. В. Численное моделирование 3Э-задач теплопроводности с фазовыми переходами на вычислительных системах с распределенной памятью // Вестн. ЮУрГУ. Сер. Вычислительная математика и информатика. 2013. Т. 2. С. 84-93.
5. Вабищевич П. Н., Васильева М. В., Павлова Н. В. Численное моделирование термостабилизации фильтрующих грунтов // Математ. моделирование. 2014. Т. 26. № 9. С. 111-125.
6. Богачев К. Ю., Мельниченко Н. С. О пространственной аппроксимации методом подсеток для задачи фильтрации вязкой сжимаемой жидкости в пористой среде // Вычислительные методы и программирование. 2008. Т. 9. С. 191-199.
7. Бутнев О. И., Горев И. В., Колесников С. С., Кузнецов В. Ю., Пронин В. А., Сидоров М. Л., Яруллин А. Д. Полностью неявная схема решения задач трехфазной фильтрации на неструктурированных сетках в пакете программ НИМФА // Вестн. кибернетики. 2015. № 3 (19). С. 56-72.
8. Ильин В. П. О численном решении прямых и обратных задач электромагнитной электроразведки // Сибир. журн. вычислит. математики. 2003. Т. 6. С. 381-394.
9. Федоренко Р. П. Введение в вычислительную физику. М. : МФТИ, 1994. 530 с.
10. Галанин М. П., Лазарева С. А. Метод конечных суперэлементов и его применение для решения задач науки и техники // Математ. моделирование. 2013. Т. 25. № 6. С. 32-40.
11. Arnold D. N. Unified Analysis of Discontinuous Galerkin methods for Elliptic Problems // SIAM J Numer Fnal. 2002. V. 3. № 5. P. 1749-1779.
12. Bastian P. A fully-coupled discontinuous Galerkin method for two-phase in porous media with discontinuous capillary pressure // Comput Geosci. 2014. V. 18. P. 779-796.
13. Даутов Р. З., Федотов Е. М. Разрывный смешанный метод Галеркина без штрафа для квазилинейных эллиптических уравнений второго порядка // ЖВММФ. 2013. Т. 53. № 11. P. 1791-1803.
14. Жалнин Р. В., Ладонкина М. Е., Масягин В. Ф., Тишкин В. Ф. Решение задач о нестационарной фильтрации вещества с помощью разрывного метода Галеркина на неструктурированных сетках // ЖВММФ. 2016. Т. 56. № 6. P. 989-998.
15. Il'in V. P. Problems of parallel solution of large systems of linear algebraic equations // J of Mathem Sci. 2016. V. 216. № 6. P. 795-804.
16. Гурьева Я. Л., Ильин В. П., Перевозкин Д. В. Алгебро-геометрические и информационные структуры методов декомпозиции областей // Вычислит. методы и программирование. 2016. Т. 17. С. 132-146.
17. Konshin I. N. Parallel computational models for estimation of actual speedup of analyzed algorithm // Russian Supercomputing Days : Proc. of the Int. Conf. MSU Publ, 2016. P. 269-280.
18. Голубева Л. А., Ильин В. П., Козырев А. Н. О программных технологиях в геометрических аспектах математического моделирования // Вестн. НГУ. Сер. Информационные технологии. 2012. Т. 10. С. 25-33.
19. Ильин В. П. DELAUNAY: технологическая среда генерации сеток // СибЖИМ. 2013. Т. 16. С. 83-97.
20. Бутюгин Д. С., Ильин В. П. Chebyshev: принципы автоматизации построения алгоритмов в интегрированной среде для сеточных аппроксимаций начально-краевых задач // Тр. Междунар. конф. ПАВТ-2014. Челябинск : изд-во ЮУрГУ. 2014. С. 42-50.
21. Butyugin D. S., Il'in V. P. Solution of problems of harmonic electromagnetic field simulation in regularized and mixed formulations // RJNAMM. 2014. V. 29. № 1. P. 1-12.
22. Intel Mathematical Kernel Library from Intel. URL: http://software.intel.com/en-us/intel-mkl (дата обращения: 30.01.2017).
23. Бутюгин Д. С., Гурьева Я. Л., Ильин В. П., Перевозкин Д. В., Петухов А. В., Скопин И. Н. Функциональность и технологии алгебраических решателей в библиотеке Krylov // Вестн. ЮУрГУ. Сер. Вычислительная математика и информатика. 2013. Т. 2. С. 92-105.
24. Савченко А. О., Ильин В. П., Бутюгин Д. С. Метод решения внешней трехмерной краевой задачи для уравнения Лапласа // СибЖИМ. 2016. Т. XIX. № 2. С. 88-99.
25. CCA: The Common Component Architecture Forum. URL: www.cca-forum.org/ (дата обращения: 30.01.2017).
