CC BY
140
Вычислительные технологии
Том 7, № 2, 2002
НЕКОТОРЫЕ ТЕНДЕНЦИИ РАЗВИТИЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
В. М. КОВЕНЯ
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: [email protected]
The state of the art in mathematical modeling and prospects of its development are analyzed using the problems of continuum mechanics and physics by way of example.
Введение
Математическое моделирование как новый способ исследования и получения новых знаний сформировалось в 70-х годах прошлого столетия на основе широкого применения математических методов для решения теоретических и практических проблем естествознания. Стимулом к его созданию и развитию послужило появление электронно-вычислительных машин, способных производить арифметические и логические вычисления со скоростью, недоступной для отдельных человеческих индивидуумов и их сообществ. Необходимость решения все более сложных задач, возникающих в науке, технике и народном хозяйстве, потребовала разработки и обоснования математических моделей, отражающих основные закономерности исследуемых явлений, и создания эффективных численных алгоритмов их решения. В свою очередь реализация этих алгоритмов на ЭВМ привела к созданию новых языков программирования, операционных систем и систем поддержки программного обеспечения, разработке новых подходов в программировании и информационных технологиях. Как следствие, появились новые типы архитектур ЭВМ. Такой взаимосвязанный процесс разработки математических моделей, численных алгоритмов, программирования и создания комплексов и пакетов программ для решения на ЭВМ этих задач, их анализа, хранения и вывода результатов расчетов и послужил основой нового научного направления в исследованиях — математического моделирования (см., например, монографии и основополагающие работы [1-10]).
Появившись как новое направление исследований, математическое моделирование породило и свои проблемы, от решения которых и зависит его развитие. Математические модели и алгоритмы, программы и комплексы программ, ЭВМ и системы поддержки для решения задач являются элементами моделирования. Их роль и место могут быть правильно оценены лишь в рамках всей цепочки моделирования, которую назовем технологической [11]. Под технологической цепочкой моделирования будем понимать совокупность ее элементов, выполняемых в определенной последовательности и составляющих полный цикл. Разумеется, для различных областей исследования эти элементы могут различаться,
© В. М. Ковеня, 2002.
поэтому за основу технологической цепочки должны быть выбраны общие для всех областей моделирования элементы. Тогда процесс моделирования может быть представлен в виде такой последовательности: исследуемое явление — математические модели — численные алгоритмы — программирование — ЭВМ — расчеты — 'результаты и их анализ, дополняющий известную триаду математического моделирования — модель — алгоритм — программа [2, 5].
Очевидно, все элементы технологической цепочки взаимосвязаны, и эта связь нелинейна, а изменение в одном из ее элементов может привести к изменению не только последующих, но и предыдущих элементов. Вначале моделирования исследователь явно или неявно проводит анализ всей цепочки, исходя из знаний об исследуемом явлении, наличия ресурсов, возможностей ЭВМ и т.д. Конечно, для различных классов исследуемых задач отдельные элементы цепочки могут быть опущены. В качестве примера приведем известный взгляд Н. Н. Яненко на разностную схему (численный алгоритм) как математическую модель для описания физического явления. Стремление получить ЭВМ с максимальной производительностью привело разработчиков ЭВМ к созданию многопроцессорных и параллельных ЭВМ. Численные алгоритмы, эффективные для однопроцессорных ЭВМ, могут не оказаться таковыми для ЭВМ новых архитектур. Разработка новых параллельных алгоритмов в свою очередь может повлиять и на выбор модели. Разумеется, примеры взаимосвязи и взаимного влияния элементов моделирования могут быть продолжены.
Сегодня с полным правом можно утверждать, что математическое моделирование наряду с физическим и натурным экспериментами является основным способом исследования и получения новых знаний в различных областях естествознания. Можно ожидать, что его роль в дальнейшем возрастет, но оно не заменит физический или натурный эксперименты, так как опыт всегда останется основой исследования. При изучении явления следует ожидать сближения различных форм исследования, взаимно дополняющих друг друга. Широкое внедрение результатов математического моделирования в различных областях естествознания обусловлено многими факторами, основными из которых являются:
— усложнение класса исследуемых задач, для изучения которых необходимо создание новых дорогостоящих экспериментальных установок или модельных объектов;
— большая стоимость обслуживания экспериментальных установок и объектов и высокие энергетические затраты на их работу и обслуживание;
— необходимость решения экологических, социальных и других проблем при их эксплуатации;
— невозможность проведения физического (химического, экономического и т.д.) или натурного моделирования в ряде областей исследования;
— сокращение сроков исследования и получения результатов, возможность его многократного и быстрого повторения или уточнения, хранения и т. д.
Развитие математического моделирования приводит к созданию автоматизированных систем для управления производством, что позволит резко увеличить производительность труда и избежать субъективного влияния "человеческого" фактора при принятии решений. Таким образом, сегодня математическое моделирование становится основным способом исследования и получения новых знаний. С другой стороны, результаты математического моделирования могут стать (и становятся) основой для широкого применения их в производстве и других областях человеческой деятельности (например, при создании САПР, экспертных систем и т.д. [12]).
В настоящей работе рассматриваются некоторые тенденции развития математического моделирования. Понимая, что охватить все области его применения невозможно, ниже
будет сделан акцент на задачах механики сплошной среды. В этой области исследования математическое моделирование получило наибольшее распространение как в силу невозможности получения решений на основе других подходов, так и в силу важности этого класса задач для развития производства (см., например, [1-12]). Как отмечалось выше, эффективность математического моделирования может быть правильно оценена в рамках технологической цепочки. Поэтому при рассмотрении ее проблем и тенденций остановимся на анализе состояния и развития ее отдельных элементов и их взаимодействии. Главное внимание будет уделено этапам выбора моделей и численных алгоритмов.
1. Физико-математические модели
Для задач механики сплошных сред в наиболее полной постановке физико-математические модели могут быть описаны интегральными законами сохранения
адУ + У Ж ¿8 = I РдУ, (1)
V Я V
выражающими связь между изменениями во времени в замкнутом объеме V некоторых величин (потоков), их изменениями через границы Б и взаимодействием с внешними источниками (стоками). Интегральные законы сохранения (например, массы, импульсов и энергии для моделей сплошной среды) являются наиболее общей формой описания движения сред, и они справедливы как для непрерывных, так и разрывных решений. Наряду с интегральной формой используется их дифференциальное представление
+ ^уЖ = ^ (2)
полученное из (1), но справедливое лишь для непрерывных решений.
Многообразие и многопараметричность исследуемых задач и их приближенных математических постановок, разномасштабность процессов приводят нас к цепочке физико-математических моделей, каждая из которых получена при определенных предположениях о характере изучаемого явления и описывает основные его закономерности. Характерной особенностью такого подхода является и многообразие уравнений, описывающих эти модели. Полученные уравнения могут принадлежать к различным типам (гиперболическим, параболическим, эллиптическим или уравнениям переменного типа), что приводит к различной постановке начальных и краевых задач. Более того, при исследовании одного класса задач тип уравнений может меняться в зависимости от характера решения. Так, например, стационарные уравнения газовой динамики являются уравнениями эллиптического типа для дозвуковых и гиперболического — для сверхзвуковых скоростей, т. е. в отдельных подобластях расчетной области необходимо решать уравнения различных типов. Естественно, что это накладывает дополнительные требования на применяемые численные методы.
Большинство процессов в механике жидкости, газа и плазмы являются нелинейными и эволюционными, и, как следствие, эти же свойства присущи описываемым ими системам уравнений. Теория таких уравнений изучена недостаточно, для большинства задач не доказаны теоремы существования и единственности, и, более того, их решения могут быть неединственными и разрывными даже при гладких начальных данных. Переход к многомерным задачам и усложнение расчетных областей (рассмотрение реальных
геометрий) вносят дополнительные трудности их постановок. При отсутствии строгих доказательств существования и единственности решений всегда остается вопрос о соответствии физико-математической модели исследуемому явлению. При недостаточной информации об исследуемом процессе возникает необходимость рассмотрения спектра моделей, учитывающих основные закономерности изучаемого явления при различных диапазонах изменения основных параметров. Таким образом, выбор и формулировка физико-математических моделей становится многопараметрической задачей, требующей для своего разрешения анализа всей цепочки моделей. С другой стороны, выбор и обоснование физико-математических моделей зависят от ряда факторов, выходящих за данный этап моделирования. До начала моделирования исследователь должен ответить на ряд вопросов. Например, какова конечная цель исследования и какие данные мы надеемся получить, решив эту задачу? Насколько важна данная задача и каково ее место в общей проблеме? Какова точность модели, с какой точностью требуется получить ответ? Какие ресурсы необходимо привлечь для решения задачи и какие математические и технологические возможности имеются у исследователя? Фактически исследователь проводит анализ всей технологической цепочки моделирования, на основании которого и делает вывод о возможности решения задачи или проблемы на основе существующих моделей, численных алгоритмов и технических средств (ЭВМ) или формулирует условия или требования, необходимые для ее разрешения. Это могут быть требования к созданию дополнительных моделей и численных методов решения или создания новых ЭВМ и т. д.
Выбор моделей, адекватно описывающих исследуемое явление или изучаемый процесс, включает их математическое обоснование и корректную постановку начально-краевых задач. По современным представлениям для задач механики и физики все классы моделей могут быть условно разделены на четыре группы (уровня):
1) аналитические приближения и линеаризованные уравнения;
2) нелинейные уравнения без учета диссипативных процессов;
3) нелинейные уравнения с учетом диссипативных процессов;
4) полные нестационарные модели, описываемые уравнениями с учетом реальных эффектов, как, например, уравнениями Навье — Стокса сжимаемого и теплопроводного газа, уравнениями многокомпонентных и многофазных сред, магнитогидродинамических моделей различного уровня и т. д.
Сегодня в зависимости от целей исследования и классов решаемых задач, важности задачи в общей проблеме, необходимой точности решения, имеющихся математических и технических возможностей исследователей и других факторов используются все группы моделей, начиная от простейших (1-й уровень) до самых полных (3-4-й уровни). Обратимся к некоторым примерам из области вычислительной аэродинамики. В приближении потенциального течения (1-й уровень) панельным методом удается получать решение задач обтекания самолетов реальных конфигураций (например, самолета Г4Г с подвесками для оружия) и распределение параметров течения на его поверхности [13]. Первые такие расчеты были выполнены в 70-х годах прошлого столетия на ЭВМ малой производительности. Однако эта модель не учитывает реальные эффекты в газе, например сжимаемость, вязкость и теплопроводность, и не позволяет получить распределения газодинамических потоков около тела. Переход к большим скоростям потребует использования нелинейных моделей, описываемых уравнениями газовой динамики (2-й уровень). Их решение может включать разрывы газодинамических величин, требует специальных методов расчета и ЭВМ большой производительности (порядка 109 флоп).
С появлением ЭВМ высокой производительности в 80-е годы расчет обтекания сверх-
звукового самолета (например, типа Г 16)в приближении уравнений Эйлера мог быть проведен за время, сравнимое со временем его полета. Однако и сегодня решить задачу нестационарного обтекания с учетом реальных свойств газа, таких как турбулентная вязкость и теплопроводность (3-4-й уровень), удается лишь для модельных форм летательных аппаратов даже с использованием современных суперЭВМ. Эти сложности связаны не только с большими затратами ЭВМ в связи с многопараметричностью задачи (различные числа Маха, Рейнольдса, геометрии обтекаемых тел и т. д.), но и в первую очередь с недостаточным обоснованием математических моделей (например, моделей турбулентности и областей их применимости) и их замыкающих соотношений. Подобные трудности возникают и в области гиперзвуковой аэродинамики, где наряду с решением отмеченных выше проблем необходимо правильно учесть и оценить влияние химических реакций, протекающих в газе около обтекаемого тела и на его поверхности при больших температурах, оценить разрушение поверхности и учесть прочностные и другие характеристики.
Конечно, при решении других задач в рамках различных приближений возникают свои проблемы, успешное решение которых определяет прогресс в развитии математического моделирования. Многопараметричность и разномасштабность исследуемых задач, их нелинейность и многомерность не позволяют сформулировать общие подходы для их постановки и получения решений. Можно говорить лишь о направлениях исследований и основных тенденциях развития, которые включают:
— использование моделей всех уровней в зависимости от целей исследования;
— применение все более сложных моделей с целью учета большего числа реальных физических эффектов исследуемых явлений;
— анализ моделей, их систематизацию и выявление некоторых классов общих моделей, пригодных для описания широкого комплекса проблем;
— дальнейшее математическое обоснование физико-математических моделей и корректных постановок начально-краевых задач.
Отметим, что упрощенные модели получают, как правило, из моделей более высокого уровня при различных предположениях о характере исследуемого явления. Таким образом, взяв за основу более полную модель, можно получить из нее цепочку упрощенных моделей. Такие полные модели, из которых могут быть получены их упрощенные приближения, будем называть накрывающими моделями. Ярким примером такого подхода служит модель, описываемая уравнениями Навье — Стокса сжимаемого теплопроводного газа (модель 4-й группы). При пренебрежении эффектами вязкости и теплопроводности получаем модель газовой динамики, справедливую для описания многих физических задач. Для сильно вязких течений может быть использовано приближение пограничного слоя, полученное из уравнения Навье — Стокса при сохранении членов порядка 0(1/^) и пренебрежении членами более высокого порядка малости. В рамках этого же подхода получают модели вязкого ударного слоя, модели "параболизованных" уравнений Навье — Стокса и т.д. (см., например, [15]). Можно ожидать, что и для других классов физико-математических моделей могут быть построены цепочки упрощенных моделей на основе базовой накрывающей модели. Такой подход позволит сократить число рассматриваемых моделей и сосредоточиться на изучении базовых моделей, справедливых для описания целых классов задач.
В заключение данного раздела заметим, что достигнутый сегодня уровень математического моделирования (по крайней мере, в области механики сплошной среды) основывается на гигантском заделе работ математиков-теоретиков, механиков и физиков предыдущих поколений. Новые достижения в математическом моделировании будут основаны на новых
теоретических результатах в этих областях и их освоении математиками-вычислителями.
2. Численные алгоритмы
Нелинейность большинства исследуемых задач и соответствующих им систем дифференциальных уравнений, как отмечалось выше, не позволяет получить их точные решения, кроме некоторых частных случаев. Более того, такие решения не всегда существуют, поэтому для их нахождения применяются приближенные и численные методы. Первые основаны на некотором представлении решения, исходя из известных предположений о характере решения. Так, например, асимптотические методы в аэродинамике используют представление решения в виде его разложения по малому параметру (по обратным числам Рейнольдса или числам Маха М для гиперзвуковых скоростей и т.д.). Несмотря на простоту получения таких решений, эти подходы применимы лишь для сравнительно простых задач с преобладанием одного типа течений. Всегда остается вопрос об области применимости этих подходов и их обосновании. Асимптотические методы широко использовались на начальном этапе моделирования. Переход к более сложным моделям потребовал разработки и применения численных алгоритмов, справедливых для решения многомерных задач в рамках различных физико-математических моделей.
На современном этапе развития математического моделирования широкое применение получили различные численные методы, основными из которых являются:
— конечно-разностные методы (МКР) и метод конечных объемов (МКО);
— метод конечных элементов, метод граничных элементов;
— специальные методы (метод частиц в ячейках, метод статистических испытаний и т. д.).
Остановимся на анализе развития численных алгоритмов на примере метода конечных разностей и конечных объемов, получивших в силу своей универсальности широкое распространение при решении задач механики сплошной среды. Основы теории конечно-разностных и конечно-объемных методов изложены во многих работах (см., например, [14-22]). Напомним основные понятия. Пусть в области П = П(х) с границей 7 требуется найти решение краевой задачи
Ьи = {, /и7 = ф, (3)
где Ь, I — дифференциальные, интегральные или алгебраические операторы х=(£, Х\..., хN); и = и(х) — вектор искомых функций, а {, ф — векторы правых частей. Предположим, что задача (3) поставлена корректно. Переход от математической постановки задачи к ее численному решению включает следующие этапы:
— замену области П непрерывного аргумента х его дискретным аналогом Пн;
— замену (аппроксимацию) функций и, {, ф непрерывного аргумента х дискретными функциями ин,/н,<Рн;
— аппроксимацию исходных операторов Ь,1 их дискретными аналогами Ьн,1н■
После их выполнения приходим к системе дискретных уравнений
Ьнин = ¡н, lh.Uk |7й= <рн, (4)
которые и называют разностной схемой. Как следствие, возникает вопрос о соотношении (близости) решений исходной дифференциальной (3) и дискретной (4) задач. Существуют и другие требования к постановке задачи (4). Так как каждый этап постановки дискретной задачи — аппроксимация области и функций непрерывного аргумента и исходных
операторов их конечномерными аналогами — может проводиться различными способами, то в итоге для решения задачи (3) может быть получено целое семейство постановок (4), обладающих, вообще говоря, различными свойствами. Требования к разностным (конечно-объемным) схемам также могут быть различными. Сформулируем основные из них:
— сходимость решения дискретной задачи (4) к решению исходной задачи (3);
— достаточная точность расчета;
— экономичность алгоритма;
— универсальность алгоритма, т. е. возможность его адаптации к различным физико-математическим моделям;
— адаптация алгоритма к различным архитектурам ЭВМ.
Очевидно, эти требования могут быть дополнены, например, адекватностью или близостью свойств разностной схемы свойствам исходной задачи [21] или свойствам консервативности, однородности алгоритма и т. д. Более того, сами требования могут быть интерпретированы по-разному. И, конечно, эффективность моделирования на этапе "численные алгоритмы" может быть полно оценена лишь в рамках всей технологической цепочки. Заметим, что требования к численным алгоритмам противоречивы и в зависимости от целей исследования часть из них может быть опущена. Например, при необходимости получения решения с высокой точностью экономичность алгоритма или его универсальность могут оказаться не столь важными. Таким образом, удовлетворение всех требований приводит нас к задаче оптимизации. Как и всякая задача на оптимизацию, в зависимости от ограничений (критериев) она может иметь одно или несколько решений или не иметь их совсем. Отсюда следует очевидный вывод о невозможности построения единого универсального алгоритма для решения классов задач и необходимости создания различных алгоритмов в зависимости от целей исследования.
Проведем анализ основных требований к численным алгоритмам. Для обеспечения сходимости численного решения к решению исходной задачи необходимо удовлетворять, как следует из теории разностных схем [16], условиям аппроксимации и устойчивости (корректности) разностного решения. Доказательство этих утверждений является достаточно сложной задачей, особенно для нелинейных уравнений, и является одной из актуальных в теории разностных схем. Отметим лишь, что строгие доказательства сходимости, например, методом априорных оценок или на основе дифференциального приближения на сегодня получены лишь для узкого класса разностных схем для отдельных систем уравнений. Поэтому сходимость к решению получают, как правило, на основе численных расчетов, что, конечно, не может служить доказательством.
Требования к точности расчета для различных физико-математических задач могут быть разными в зависимости от цели моделирования. И, конечно, точность расчета должна быть согласована с точностью выбранной физико-математической модели. Повышение точности расчета как одной из важнейших характеристик численного расчета также может быть достигнуто различными способами:
— измельчением шагов сетки в разностных схемах или расчетных ячеек в МКО;
— применением неравномерных и подвижных сеток;
— построением схем повышенного порядка.
Отметим также и другие способы повышения точности расчетов, которые развиваются в последнее время: повышение точности за счет экстраполяции численных решений, полученных на последовательности сеток [23]; использование информации о гладкости решения (алгоримы без насыщения [24]); использование точных решений в ячейке с
кусочно-постоянными или кусочно-линейными начальными данными при конструировании разностных схем [25]; выделение главных особенностей решения, например, выделение головной ударной волны в задачах сверхзвукового обтекания и т. д.
Измельчение шагов сетки или расчетных ячеек решения многомерных задач не является эффективным способом повышения точности расчета из-за степенного возрастания числа узлов в расчетной области и соответствующего роста числа арифметических операций. Хотя такой подход не требует переделки алгоритма и программ, на практике он используется сравнительно редко. При решении задач чаще всего применяются неравномерные сетки, в том числе сетки, сгущающиеся в областях больших градиентов. При известной информации о поведении решения вводится преобразование координат, сгущающее ячейки в областях, содержащих особенности решения (пограничные слои, ударные волны и т.д.). Такой подход оказывается очень эффективным, так как позволяет при сокращении числа расчетных ячеек существенно повысить точность расчета, хотя число арифметических операций возрастает в силу усложнения аппроксимации уравнений в новых преобразованных координатах. Однако для большинства задач области больших градиентов и других особенностей решения, как правило, априори неизвестны и могут быть получены лишь в процессе решения. Еще более сложная ситуация возникает в нестационарных задачах, где решение меняется во времени и расчетные сетки должны зависеть от времени. Для решения таких задач исходные уравнения в новых преобразованных координатах должны дополняться нестационарными уравнениями для определения закона движения сетки [25, 27-30]. На начальном этапе моделирования сетки строились, как правило, вручную. Усложнение расчетных областей и переход к решению многомерных задач потребовали разработки специальных методов построения или генерации сеток, удовлетворяющих определенным требованиям. Сегодня при решении многомерных задач в областях со сложной геометрией задача построения эффективных сеток становится одной из центральных и, по оценкам специалистов, основные затраты ресурсов приходятся как раз на построение преобразований с заданными свойствами. Некоторые подходы решения этих проблем приведены в монографиях и обзорах [25-32].
Отметим еще один способ повышения точности расчета — на вложенных или адаптивно-встраивающихся сетках. Получив предварительное решение, например, на равномерной сетке и установив основные особенности, проводят перерасчет решения в этих подобластях на сгущающихся сетках. В качестве граничных условий на границах этих подобластей сохраняются значения величин, полученных на начальной сетке путем их интерполяции на вновь введенные узлы. Уточненное решение вновь проектируется на старую сетку. Затем процесс вычислений повторяется. Хотя данный подход недостаточно обоснован, он может послужить основой для получения предварительного решения, а затем процесс вычислений может повторяться во всей расчетной области на таких неструктурированных нерегулярных сетках до получения решений с требуемой точностью. Некоторые методы решения уравнений на неструктурированных сетках изложены в монографии [33].
Наряду с неравномерными сетками в последние десятилетия широко используются схемы повышенного порядка, главным образом двух типов — схемы на расширенном шаблоне и так называемые компактные схемы высокого порядка на трехточечном шаблоне. В первом случае при построении схем сквозного счета в задачах с разрывными данными для уменьшения осцилляций численного решения применяют различные способы монотонизации: добавляют в исходные уравнения или разностную схему диссипативные члены, которые наряду с погашением осцилляцией сглаживают и решение. В схемах типа ТУБ монотонность решения достигается путем введения нелинейной численной диссипации,
обеспечивающей обратную связь схемы с решением [34-36]. Эти схемы, как правило, согласованы с энтропийным неравенством и имеют второй или более высокие порядки аппроксимации. Сглаживание решения происходит на небольшом числе ячеек. Имеется множество модификаций схем (ЕМО-схемы, схемы с корреляцией, схемы ограниченной антидиффузии и т.д. [35]), позволяющих сохранить их монотонность, повышенный порядок на разрывах, консервативность и т. д. Отметим, что использование расширенного шаблона в этих схемах приводит к необходимости задания фиктивных слоев и дополнительных граничных условий, отсутствующих в исходной постановке задач или изменения аппроксимации в приграничных узлах, нарушающей однородность схем.
В компактных схемах аппроксимация повышенного порядка (от 3-го и выше) достигается на трехточечном шаблоне за счет специальной аппроксимации уравнения и введения операторов осреднения, исключающих погрешности низших порядков (подробнее см. [31, 32]). Полученные на этой основе разностные схемы слабо немонотонны и при неявной аппроксимации реализуются трехточечными прогонками для одномерных задач. В качестве недостатка схем повышенного порядка отметим значительное усложнение аппроксимации, особенно в многомерном случае, что приводит к увеличению затрат ресурсов ЭВМ при их реализации. Однако для большинства многомерных задач применение схем повышенного порядка является единственной возможностью их решения с достаточной точностью. Отметим также многослойные схемы, имеющие суммарный повышенный порядок [37], и схемы, использующие следствия из исходных уравнений, т.е. гладкость решения [38].
Разновидностью метода конечных разностей является метод конечных объемов, основанный на аппроксимации исходных уравнений в интегральной форме (см., например, [1620]). Возможность выбора различных форм расчетных ячеек при аппроксимации расчетных областей обеспечила его широкое применение при решении задач в областях сложной геометрии, в том числе в многосвязных областях. Как правило, в МКО не производятся преобразования координат, приводящие к усложнению уравнений. Исходные уравнения аппроксимируются для каждой исходной ячейки, т. е. получаемые схемы являются консервативными. Порядок аппроксимации зависит от точности аппроксимации объемных и поверхностных интегралов, что позволяет проще строить схемы повышенного порядка. Получаемая в результате система алгебраических уравнений является нелинейной и может быть реализована одним из известных итерационных методов или после линеаризации — методами, разработанными для решения конечно-разностных схем, в том числе схем типа ТУБ, ЕМО и т. д.
Требование экономичности алгоритма всегда являлось одним из главных при построении разностных схем и схем МКО и понималось как минимизация числа арифметических операций на решение задачи. Приведем оценки затрат ресурсов ЭВМ при решении многомерных нестационарных задач [11]. Пусть М — размерность задачи по пространству; т — количество уравнений (неизвестных функций); I] — число узлов в направлении X]; д — усредненное число арифметических операций на узел сетки; N — число шагов по времени. Тогда общее число операций, необходимых для решения задачи, равно
Я = mqNI, (5)
где I = II • 12 •... • 1м. Конечно, эта формула не включает ряд факторов, влияющих на затраты ресурсов, таких как число внутренних итераций, заданная точность расчета, влияние начального приближения и т. д. Будем полагать, что они априори учтены в коэффициенте д. Выберем для примера средние значения параметров при решении задач существующими методами: ] = I и 102 ^ 103, т и 10 ^ 102, N и 103 ^ 104, д и 102 ^ 104. Для
M = 3 (трехмерный случай) число операций равно Q = 1012 ^ 1019. Обозначим через R быстродействие ЭВМ, т. е. количество операций, выполняемых за 1 с. Тогда
T = Q/R (6)
есть время решения задачи. Соответственно затраты оперативной памяти могут быть определены формулой
L = mil, (7)
где I — число временных слоев, необходимых для хранения в оперативной памяти ЭВМ (обычно l œ 2). По оценкам специалистов, использование результатов математического моделирования в проектировании может быть оптимальным, если расчет одного варианта не превышает 15 мин (~ 103 с). Для выполнения этого условия быстродействие ЭВМ и оперативная память для решения трехмерных задач должны составлять соответственно R = Q/T œ 109 ^ 1016 флоп, L = 2 х (107 ^ 1011) байт.
Заметим, что в проекте NAS, разрабатываемом в 80-е годы и ориентированном на решение задач аэродинамики, параметры ЭВМ должны были составлять R œ 109 ^ 1010 флоп, L œ 2.5 х 108 байт, т. е. планировалось создать вычислительную систему решения пространственных задач аэродинамики в приближении уравнений газовой динамики и упрощенных уравнений Навье — Стокса (2-3-я группы моделей).
При определении быстродействия R = Q(T) нигде не оговаривалось, используются ли однопроцессорная ЭВМ или многопроцессорный вычислительный комплекс. Полагалось, что быстродействие R не зависит ни от типа алгоритмов, размерности задачи и потребляемых ресурсов памяти, ни от других факторов. На практике, конечно, быстродействие ЭВМ зависит от многих факторов. Для приближенного их учета формула быстродействия может быть подправлена и представлена в виде
R' = sR, (8)
где R — пиковая производительность системы, а s — средний коэффициент загрузки процессоров (s < 1), значение которого может значительно изменяться для различных многопроцессорных вычислительных систем.
Сокращение затрат на решение задачи согласно (6) может быть достигнуто за счет уменьшения Q или увеличения быстродействия R. Число арифметических операций Q и определяет экономичность алгоритма. Дадим одно из определений экономичности, следуя соотношению (5). Пусть для исследуемого класса задач, решаемых каким-либо численным алгоритмом, известны число уравнений m и среднее число арифметических операций на узел сетки q. Введем I0 = NI — общее число шагов по времени и пространству. Тогда Q = Q(I0). Численный алгоритм назовем экономичным, если Q = Q(I0) — линейная функция. Как известно, классы разностных схем делятся на явные и неявные, причем существуют несколько критериев для их определения. Воспользуемся одним из них (см. [11]). Пусть разностная схема (4) есть двухслойная схема, т. е. представима в виде
CX+1 = Cml (9)
где C0 и C1 — разностные операторы, связывающие значения функций Uh в узлах сетки на временных или итерационных слоях — n-м и (n + 1)-м соответственно. Назовем схему (9) явной, если C0 = I, и неявной в противном случае. В соответствии с этим определением в явных схемах значения функций u^+1 на новом (n + 1)-м слое определяются лишь через
известные значения на слое п. Подобное определение может быть дано и для многослойных схем и т. д.
Явные схемы, как правило, условно устойчивы, т. е. разностная задача корректна при выполнении определенных соотношений между временным (итерационным) т и пространственными к] шагами сетки. Для неявных схем эти ограничения отсутствуют или они более слабые, чем для явных схем. В соответствии с введенным выше определением экономичности явные схемы не являются экономичными. Проиллюстрируем вышесказанное на простейшем примере. Пусть (4) — явная разностная схема, аппроксимирующая нестационарное уравнение теплопроводности в области П = {0 < £ < 1, 0 < X < 1}. Введем шаги разностной сетки т = 1/^ к = 1/1, где N,1 — число узлов сетки в направлениях £ и X. Явная схема устойчива при выполнении условия т < к2/2 = К1-2. Тогда число арифметических операций на решение задачи равно Я = qNI и Кд13. При уменьшении шага пространственной сетки к (увеличении числа узлов I) в р раз число операций увеличивается в р3. Подобная ситуация возникает и при использовании явных схем для решения гиперболических уравнений, для которых число опрераций увеличивается по квадратичному закону при линейном увеличении числа шагов сетки.
Остановимся на вопросах универсальности алгоритма и возможности его обобщения на многомерный случай или более сложные физико-математические модели. Для явных схем обобщение алгоритмов на многомерный случай не представляет значительных сложностей, однако проблема экономичности выходит на первый план, так как увеличение размерности задачи приводит для явных схем к более жестким ограничениям на устойчивость. Для неявных схем обобщение алгоритмов на многомерный случай также может быть проведено подобно тому, как это сделано для явных схем. Однако для нелинейных задач получающаяся система разностных уравнений является нелинейной, что требует разработки специальных алгоритмов, например введения итераций по нелинейности или линеаризации схем. Реализация неявных схем даже для многомерных линейных задач значительно усложняется. Для иллюстрации этого сошлемся на разностную схему с весами для численного решения линейного уравнения теплопроводности (см., например, [14]). Эта схема безусловно устойчива при а > 0.5 и аппроксимирует уравнение теплопроводности с порядком 0(т2 + к2) при а = 0.5. В одномерном случае она реализуется скалярными трехточечными прогонками и требует независимо от узлов сетки 8 арифметических операций на узел сетки. В двумерном случае ее решение также может быть получено методом прогонки — методом матричной прогонки. Реализация методом прогонки требует уже обращения матриц размером I х I в каждом узле сетки, т.е. — д!3 операций на решение задачи на одном временном слое. И, очевидно, для большого числа узлов эта схема становится неэкономичной. Подобная ситуация возникает и для гиперболических уравнений.
Таким образом, прямое обобщение схем, эффективных для решения одномерных задач, на многомерный случай чаще всего приводит к усложнению их реализации и потере экономичности. В действительности переход от решения одномерных уравнений к многомерным потребовал разработки новых численных алгоритмов, пригодных для эффективного решения задач любой размерности. И такие алгоритмы активно разрабатывались в 70 - 80-е годы (см., например, [1-3, 11-21, 32]). При их построении использовались два основных метода — метод факторизации и метод расщепления (см., например, [12-22]). Их суть состоит в сведении исходной многомерной задачи к последовательности их одномерных аналогов или более простых задач. В методе факторизации это проводится после аппроксимации исходных уравнений и получения базовых схем, которые затем представляются в виде приближенного произведения простых одномерных разностных аналогов. В методе
расщепления исходная задача представляется в виде совокупности (слабой аппроксимации) более простых одномерных задач, которые затем аппроксимируются явными и неявными разностными схемами [3, 14, 15]. Отметим, что неявные схемы безусловно устойчивы для уравнений любой размерности независимо от типа уравнений, а схемы приближенной факторизации теряют это свойство в трехмерном случае для гиперболических уравнений. Заметим, что введение расщепления или факторизации в разностных схемах приводит к появлению дополнительных членов порядка 0(т2), отсутствующих в исходных уравнениях, что может приводить к понижению точности решения или увеличению числа итераций для получения стационарного решения методом установления.
Наконец, обратимся к последнему, но не по важности, требованию к численным алгоритмам — адаптации алгоритмов к различным архитектурам вычислительных систем.
Основные численные алгоритмы, широко использующиеся для решения различных классов задач, были созданы в эпоху однопроцессорных ЭВМ и предполагали последовательное выполнение арифметических операций. Создание векторных ЭВМ, использующих однородные вычисления для векторов большой размерности (на ряде арифметических устройств), потребовало незначительной модификации алгоритмов, реализуемой на программном уровне. Появление многопроцессорных вычислительных систем различных типов потребовало переоценки разработанных численных алгоритмов, оценки их эффективности и возможности их адаптации к новым архитектурам. Сложность решения этой проблемы состоит еще в том, что новые архитектуры ЭВМ создаются на основе разных подходов в их проектировании. Поэтому на сегодня возможности адаптации большинства существующих методов к различным типам вычислительных систем изучены недостаточно. Более простая ситуация с явными схемами, которые достаточно легко могут быть распараллелены. Например, в силу явного способа вычислений на новом шаге исходная задача может быть разбита на отдельные ее сегменты или модули, каждый из которых может быть решен независимо друг от друга (с точностью до вычислений на границах сегментов и передачи данных от одного к другим). Таким образом, появление новых вычислительных средств привело к возникновению нового направления в вычислительной математике — параллельных алгоритмов решения многомерных задач.
Сформулируем основные проблемы, возникающие на этапе разработки численных алгоритмов:
— построение расчетных сеток с заданными свойствами (например, сгущающихся в областях больших градиентов ) для решения многомерных задач в сложных многосвязных геометрических областях;
— адаптация существующих и разработка новых экономичных численных алгоритмов решения задач в рамках различных физико-математических моделей, в том числе и для ЭВМ различных архитектур;
— развитие математического аппарата и его применение для обоснования численных алгоритмов.
Касаясь элементов технологической цепочки моделирования ЭВМ — программирование, отметим, что и в этих областях произошли, быть может, более значительные изменения. Современное состояние языков программирования и перспективы их развития изложены в [39]. Быстродействие ЭВМ и вычислительных систем постоянно возрастает с экспоненциальной скоростью как за счет совершенствования элементной базы, так и путем создания принципиально новых архитектур вычислительных структур и систем моделирования. Наряду с созданием суперЭВМ с пиковой производительностью до 1012 — 1013 флоп разрабатываются вычислительные средства кластерного типа, объединяющие в сети
ЭВМ различной производительности: персональные ЭВМ, рабочие станции и т.д. Следует ожидать создания новых вычислительных средств и усиления их влияния на численные алгоритмы и модели, создания экспертных систем и систем моделирования в различных областях человеческой деятельности.
В заключение хотелось бы отметить, что вышеизложенное является лишь взглядом автора на проблему и, конечно, не является полным. Автор будет благодарен всем читателям за замечания или дополнения, которые могут быть направлены по электронному адресу: [email protected].
Список литературы
[1] ЯнЕнко Н. Н. Математика. Механика: Избр. тр. М.: Наука, 1991. 416 с.
[2] Самарский А. А. Математическое моделирование и вычислительный эксперимент // Вестник АН СССР. 1979. №5. C. 38-49.
[3] Белоцерковский О. М. Численное моделирование в механике сплошных сред. М.: Наука. 1984. 520 с.
[4] ЯнЕнко Н. Н., Коновалов А. Н. Технологические аспекты численных методов математической физики // Acta Univ. Carot. Math. Phys. 1974. Vol. 2. P. 47-53.
[5] ЯнЕнко Н. Н., Коновалов А.Н. Модульный принцип построения программ как основа создания пакета прикладных программ решения задач механики сплошной среды // Комплексы программ мат. физики. ВЦ СО АН СССР. Новосибирск, 1972. C. 46-54.
[6] Комплексы программ математической физики / Под редакцией Н. Н. Яненко. Новосибирск: ВЦ СО АН СССР, 1972, 1973, 1976.
[7] Адрианов А. Л., Бабенко К. И., Забродин А. В. и др. О структуре вычислителя для решения задач обтекания. Комплексный подход к проектированию // Вычисл. процессы и системы. М.: Наука, 1985. Вып. 2. C. 13-62.
[8] Ершов А. П., Ильин В. П. Пакеты программ, технология решения прикладных задач. Новосибирск, 1978 (Препр. / АН СССР. Сиб. от-ние. ВЦ; №21).
[9] Карпов В. Я., Корягин Д. А., Самарский А. А. Принципы разработки пакетов прикладных программ для решения задач математической физики // Журн. вычисл. математики и мат. физики. 1978. Т. 18, №2. С. 458-467.
[10] КовЕня В.М., ЯнЕнко Н. Н. Некоторые проблемы развития пакетов программм для решения задач аэродинамики // Числ. методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. от-ние. ВЦ; ИТПМ. 1979. Т. 10, №3. С. 89-98.
[11] КовЕня В. М. Методы вычислений (дополнительные главы): Учеб. пособие. Новосибирск: НГУ. 1995. 92 с.
[12] Белоцерковский О. М. Математическое моделирование на суперкомпьютерах (опыт и тенденции) // Журн. вычисл. математики и мат. физики. 2000. Т. 40, №8. С. 1221-1236.
[13] Современные проблемы газовой динамики / Под ред. У. Х. Т. Лоха. М.: Мир, 1971. 403 с.
[14] ЯнЕнко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967. 196 с.
[15] КовЕня В. М., ЯнЕнко Н. Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981. 304 с.
[16] Самарский А. А. Теория разностных схем. М.: Наука, 1989. 614 с.
[17] Годунов С. К. РяБЕнький В. С. Разностные схемы. (Введение в теорию). М.: Наука, 1973. 400 с.
[18] Марчук Г. И. Методы вычислительной математики. М.: Наука, 1980. 536 с.
[19] Флетчер К. Вычислительные методы в динамике жидкостей. М.: Мир, 1991. Т. 1,2.
[20] Самарский А. А., Николаев е. С. Методы решения сеточных уравнений. М.: Наука, 1978. 501 с.
[21] Шокин Ю.И., ЯнЕнко Н. Н. Метод дифференциального приближения. Новосибирск: Наука, 1985. 364 с.
[22] Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Наука, 1987. 598 с.
[23] Марчук Г. И., ШАйдуров В. В. Повышение точности решений разностных схем. М.: Наука, 1979. 318 с.
[24] БАБЕнко К. И. Основы численного анализа. М.: Наука, 1986. 744 с.
[25] Годунов С. К., Забродин А. В., Иванов М. Я. и др. Численные методы решения многомерных задач газовой динамики. М.: Наука, 1989. 614 с.
[26] Головачев Ю. П. Численное моделирование течений вязкого газа в ударном слое. М.: Наука, 1996. 374 с.
[27] LisEIKIN V. D. Grid Generation Method. Berlin: Springer-Verl., 1999. 362 p.
[28] ЛисЕйкин В. Д. Обзор методов построения структурных адаптивных сеток // Журн. вычисл. математики и мат. физики. 1996. Т. 36, №1. С. 3-41.
[29] КругляковА Л. В., Неледова А. В., Тишкин В.Ф., Филатов А. Ю. Неструктурированные адаптивные сетки для задач математической физики (обзор) // Мат. моделирование: Сб. науч. тр. / РАН Сиб. от-ние. ВЦ; ИТПМ. 1998. Т. 10, №3. С. 93-115.
[30] ГильмАнов А. Н. Методы адаптивных сеток в задачах газовой динамики. М.: Наука, 2000. 248 с.
[31] Толстых А. И. Компактные разностные схемы и их применение в задачах гидродинамики. М.: Наука, 1996. 230 с.
[32] Паасонен В. И. Схема третьего порядка аппроксимации на неравномерной сетке для уравнений Навье — Стокса // Вычисл. технологии. 2000. Т. 5, №5. С. 78-85.
[33] Самарский А. А. и др. Разностные схемы на нерегулярных сетках. Минск: Критерий, 1996. 273 с.
[34] HARTEN A. A high resolution scheme for the computation of werk solution of hyperbolic conservatin laws //J. Comput. Phys. 1983. Vol. 49, No. 3. P. 357-393.
[35] Карамышев В. Б. Монотонные схемы и их приложения в газовой динамике. Новосибирск: НГУ, 1994. 109 с.
[36] Пинчуков В. И., Шу Ч.-В. Численные методы высших порядков для задач аэрогидродинамики. Новосибирск: Изд-во СО РАН, 2000. 231 с.
[37] Русанов В. В. Разностная схема третьего порядка точности для сквозного расчета разрывных решений // Докл. АН СССР. 1968. Т. 180, №2. С. 1303-1305.
[38] Тушева Л. А., Шокин Ю. И., ЯнЕнко Н. Н. О построении разностных схем повышенного порядка аппроксимации на основе дифференциальных следствий // Некоторые проблемы вычислительной и прикладной математики. Новосибирск: Наука, 1975. С. 184-195.
[39] Кауфман В.Ш. Языки программирования. Концепции и принципы. М.: Радио и связь, 1993. 432 с.
Поступила в редакцию 4 декабря 2001 г.
