Numerical modeling of motion of near-earth objects in a parallel computing environment.pdf Введение Впервые разработанная в Томске «Численная модель движения искусственных спутников Земли (ИСЗ) [1] была представлена научной общественности в 1982 г. на Объединенных научных чтениях по космонавтике. С того времени данная модель неоднократно претерпевала различные модификации, связанные с переходом на новые вычислительные средства, с выходом новых стандартов Международного астрономического союза (МАС), с разработкой новых более эффективных интеграторов, с расширением задач, решаемых численным моделированием. В конце 90-х гг. ХХ в. в связи с переходом на персональные компьютеры численная модель подверглась значительной переработке и в виде, близком к существующему, была представлена в 2007 г. [2], а в 2009 г. программа перенесена [3] на кластер «СКИФ Cyberia» Томского госуниверситета. В результате этого стало возможным одновременно прогнозировать движения большого числа объектов. В 2010 г. была добавлена возможность вычисления параметра MEGNO [4], предназначенного для оценки хаотичности в динамике ИСЗ. В 2017 г. проведено уточнение моделей действующих сил [5] в соответствии с рекомендациями МАС и интегратор Эверхарта [6, 7] заменен интегратором Gauss32 [8]. В 2020 г. программа была дополнена [9] возможностью вычисления в процессе численного интегрирования вековых частот, необходимых для исследования резонансной динамики околоземных объектов. В представляемой здесь версии численной модели проведен ряд важных модификаций, которые позволили улучшить производительность программы. Основными модификациями стали замена интегратора Gauss32 [8] на более эффективный новый коллокационный интегратор Lobbie [10], разработанный в НИИ ПММ ТГУ, и смена способа распараллеливания. В настоящее время численная модель существует в двух версиях: одна для персонального компьютера и другая - для кластера, которую будем называть «Численной моделью движения систем ИСЗ». 1. Структура «Численной модели движения систем ИСЗ» Представляемая версия «Численной модели движения систем ИСЗ» содержит четыре основных программных блока, предназначенных для решения следующих задач: - высокоточного моделирования движения ИСЗ [5]; - исследования хаотичности движения околоземных космических объектов [4, 5] (с расчетом параметра MEGNO [11, 12] по методике [13]); - определения параметров движения ИСЗ по данным измерений [5]; - исследования резонансной динамики околоземных объектов [9]. Поскольку для исследования хаотичности в динамике ИСЗ нет необходимости учитывать тонкие эффекты, которые важны при высокоточном моделировании и при определении орбит, во втором программном блоке список сил уменьшен. Разницу в учете возмущающих сил можно увидеть в табл. 1. Таблица 1 Характеристика составляющих модели сил Возмущающие факторы Составляющие модели сил Блоки для высокоточного прогноза движения ИСЗ и определения орбит Блок для расчета MEGNO Гармоники геопотенциала Модель EGM2008 (рекомендовано IERS Conventions2010), алгоритм Л. Каннингема [14] Дополнительные возмущения геопотенциала (рекомендовано IERS Conventions2010). Вековые изменения первых зональных гармоник; возмущения от приливных деформаций центрального тела; возмущения от океанических и полюсного приливов; возмущения от океанического полюсного прилива [15] Не учитываются Влияние Солнца, Луны и планет Модель притяжения точечной массы [16] с использованием фонда DE/LE421 Гармоники селенопотенциала Модель LP150Q [17] Не учитываются Радиационные силы Конусная модель светового давления с тенью и полутенью; учет искажения диска Солнца и сжатия Земли в атмосфере; учет солнечного затмения с эффектом потемнения диска Солнца к краям; переотраженное и тепловое излучение от Земли [18,19]; эффект Пойнтинга - Робертсона [20, 21] Конусная модель светового давления с тенью и полутенью [20]; эффект Пойнтинга - Робертсона [20, 21] Релятивистские эффекты Шварцшильдовские возмущения; эффекты Лензе-Тирринговой прецессии; релятивистские квадрупольные члены [22, 23] Не учитываются Возмущения от сопротивления атмосферы Модель атмосферы NRLMSISE-00 [24] Статическая модель атмосферы [20] 2. Новый коллокационный интегратор Lobbie для численного решения задач орбитальной динамики ИСЗ Основные формулы интегратора Lobbie Интегратор Lobbie является программной реализацией нового коллокационного метода на разбиениях Лобатто, который предназначен для численного решения смешанных систем дифференциальных уравнений орбитальной динамики второго и первого порядков с возможностью автоматического выбора величины шага интегрирования. Примечательной особенностью интегратора Lobbie является то, что его теоретическая основа, как и программная реализация, универсальна для любого порядка. Практически порядок определяется количеством узловых значений разбиения Лобатто, через которые выражаются все остальные константы пошаговой схемы интегрирования. Интегратор Lobbie обладает геометрическими свойствами, являясь симметричным и орбитально устойчивым. Кроме того, он позволяет на каждом шаге легко конструировать приближенное аналитическое решение (в виде полинома), чем удобно пользоваться для частого вывода результатов на плотной временной сетке. Подробное описание интегратора Lobbie можно найти в работе [10], в данном разделе напомним только основные формулы данного интегратора. Пусть орбитальное движение ИСЗ описывается смешенной системой дифференциальных уравнений второго и первого порядков (1) где - вектор положения ИСЗ; - время; штрих означает полную производную по времени; - вектор вспомогательных переменных (например, для вычисления показателей хаотичности); и - известные вектор-функции времени и динамических переменных , , ; , , - значения динамических переменных на некоторую начальную эпоху . Необходимо определить динамическое состояние ИСЗ и вектор вспомогательных переменных на момент ( - величина шага интегрирования): (2) Схема интегрирования в любом коллокационном методе конструируется посредством прямого интегрирования полиномиальных интерполянтов правых частей дифференциальных уравнений [25-28]. В качестве таких интерполянтов примем полиномы Ньютона (3) на разбиении Лобатто . Здесь - безразмерная временная переменная; - узловые значения ; . Разделенные разности в (3) определяются из узловых значений функций и по рекуррентным формулам (4) Интегрируя интерполянты (3) по на отрезках и и учитывая, что и , получаем приближенные решения (2) системы (1) на момент : (5) где (6) Обозначим (7) -кратный интеграл от -й базисной функции Ньютона. Если величины рассматривать как элементы некой матрицы размера , то константы коллокационного метода - интегралы от базисных функций Ньютона (6) - будут составлять первые ее два столбца: Между тем элементы матрицы для произвольного значения вычисляются построчно через рекуррентные соотношения: (8) Порядок интегратора зависит от количества узловых точек . На разбиении Лобатто . Узловые значения определяются из решения алгебраического уравнения Промежуточные решения , и в (5) задаются неявным образом, поскольку от них же зависят разделенные разности и . Поэтому промежуточные решения вычисляются итерационно. До итераций известны решения , и , а также разделенные разности и на первой точке коллокаций . В начале итерационного процесса на второй точке коллокаций определяется группа решений , и , а из них - и , по которым уточняются разделенные разности и в соответствии с рекуррентными формулами (4). Затем таким же образом уточняются разделенные разности и на третьей точке коллокаций и так далее. После поочередного уточнения всех разделенных разностей на шаге итерация повторяется. Итерационный процесс продолжается до тех пор, пока не выполнится неравенство (9) где - решение на предыдущей итерации, а - малая величина, которая задает точность сходимости. Впрочем, количество итераций на шаге можно ограничить заданным числом, не учитывая выполнение условия (9). В таком случае нужно иметь в виду, что получаемое решение уже не будет соответствовать задаваемому порядку, а метод лишится геометрических свойств. В качестве начальных приближений разделенных разностей и принимаются их оценки, получаемые из рекуррентных формул (4) по экстраполированным значениям функций и : (10) На первом шаге при отсутствии таких оценок итерационный процесс начинается буквально с нуля, т.е. при нулевых значениях разделенных разностей, и продолжается до обязательного выполнения условия (9). При многократном выводе приближенных решений на плотной временной сетке удобно и целесообразно пользоваться подручным коллокационным полиномом, накрывающим на определенном шаге временные моменты сетки, вместо того чтобы выполнять пошаговое интегрирование на каждый из этих моментов с очень мелким шагом. Коллокационные полиномы для системы (1) на произвольный момент времени можно представить в виде (11) где величины и вычисляются по формулам (8). Альтернативный способ получения временных рядов приближенных решений - это полиномиальное интерполирование на промежуточных решениях типа (3): (12) Однако точность интерполяционных полиномов (12) между точками коллокаций оказывается значительно ниже, нежели точность коллокационных полиномов (11), получаемых прямым интегрированием полиномов (3). Длина шага как параметр интегратора задается пользователем. Однако возможен режим автоматического выбора длины шага в процессе пошагового интегрирования. Стратегия выбора переменного шага заимствована из интегратора Эверхарта [6] в редакции В.А. Авдюшева [8]. Тестирование «Численной модели движения систем ИСЗ» с новым интегратором «Численная модель движения систем ИСЗ» после замены интегратора Gauss32 на Lobbie тестировалась на трех модельных объектах с почти экваториальными круговыми орбитами и большими полуосями = 8, 25 и 42 тыс. км, которые условно можно отнести к семействам низколетящих, навигационных и геостационарных космических объектов. В численном эксперименте моделировалось орбитальное движение объектов на годичном интервале времени с применением нового интегратора Lobbie и его предшественника Gauss32. Учитывались следующие возмущающие факторы: притяжение Земли (с точностью до гармоник геопотенциала 5-го порядка), Луны и Солнца, рассматриваемые как гравитирующие точки. Вычисления выполнялись в компьютерной арифметике с двойной точностью. Путем варьирования допустимой ошибки на шаге интегрирования для интеграторов 8-го и 16-го порядков получены характеристики точность - быстродействие, которые приведены на рис. 1. При этом различия между результатами интеграторов Lobbie и Gauss32 составили величины порядка вычислительных ошибок (рис. 2). Точность в векторе положения оценивалась по завершении пошагового интегрирования относительно предвычисленного положения, полученного с предельно высокой (вычислительной) точностью. Ввиду регулярности орбитального движения максимальная ошибка в векторе положения для всех трех объектов достигалась на конце временного интервала. В качестве показателя быстродействия принято количество обращений к процедуре правых частей уравнений ( ). Рис. 1. Характеристики точность - быстродействие интеграторов Lobbie и Gauss32 Как видно из рис. 1, новый интегратор практически на всех уровнях точности интегрирования (вплоть до вычислительной точности, обусловленной ошибками округления) позволяет повысить быстродействие почти в 3 раза. В действительности такая сравнительно высокая эффективность по быстродействию интегратора Lobbie достигается вследствие его лучшей локальной точности. Интеграторы Lobbie и Gauss32 численно решают уравнения разных порядков: Lobbie - второго порядка, тогда как Gauss32 - первого порядка. Фактически уравнения интегратора Gauss32 - это уравнения интегратора Lobbie, приведенные к первому порядку. Вследствие этого точность схемы интегрирования для координат в Lobbie на порядок выше, нежели в Gauss32. Кроме того, скорость сходимости итерационного процесса для определения промежуточных решений в Lobbie лучше, чем в Gauss32. Следовательно, объем вычислений в Lobbie меньше. Каждую характеристику точность - быстродействие условно можно разделить на две части: с доминирующими методической и вычислительной составляющими ошибки. Последняя обусловлена ошибками округления компьютерной арифметики. Таким образом, она задает некий предельный уровень максимально достижимой точности: дециметровая для низколетящего объекта (a = = 8 тыс. км) и приблизительно миллиметровая для двух других. Однако этот уровень точности можно существенно повысить, если вычисления выполнять с четверной точностью, хотя это, разумеется, сопряжено с большими объемами вычислений. На примере геостационарного объекта ( = 42 тыс. км) рис. 3 показывает, что в таком случае точность интегрирования может быть значительно выше нанометрового уровня. Рис. 2. Согласованность между результатами интеграторов Lobbie и Gauss32 Рис. 3. Характеристики точность - быстро¬действие для интегратора Lobbie в компь¬ютерных арифметиках с двойной и четверной точностями на примере геостационарного объекта (a = 42 тыс. км) 3. Особенности реализации «Численной модели движения систем ИСЗ» в среде параллельных вычислений Основной применяемый нами ранее и в настоящее время принцип распределения вычислений по ядрам кластера - это распределение по объектам. В предыдущей версии в программный комплекс «Численная модель движения систем ИСЗ» был заложен механизм, который заключался в том, что на каждое ядро поступал пучок с одинаковым количеством объектов (M/N, где M - число всех рассматриваемых объектов, N - число ядер) и шаг интегрирования выбирался для всего пучка, а не для каждого объекта отдельно (рис. 4, а, кругами обозначены ядра кластера). После считывания начальных условий первое ядро распределяло динамические состояния объектов по ядрам и пересылало информацию о них на все остальные N-1 ядер. Далее каждое ядро интегрировало ту систему уравнения движения совокупности объектов, данные о которых прислало первое ядро. При этом первое ядро также интегрировало свою систему уравнений движения. Рис. 4. Схема распараллеливания: а - старая версия; б - новая версия Такой способ распараллеливания приводил к зависимости точности интегрирования от количества объектов в пучке (рис. 5), что не совсем корректно, поскольку у каждого объекта своя орбитальная эволюция и оптимальный шаг интегрирования индивидуален. Такой подход помимо вытекающих вычислительных ошибок порой приводил и к интенсивному дроблению шага, следствием чего являлось увеличение времени интегрирования. На рис. 5 показана зависимость точности интегрирования от количества ядер, по которым распределялись объекты для старого способа распараллеливания, а также точность при новом методе. Для тестирования рассматривалась 1000 объектов на орбитах, равномерно распределенных по большой полуоси от 8000 до 42000 км, с нулевым наклонением и эксцентриситетом равным 0.001. В табл. 2 приведены оценки временных затрат при различном количестве ядер для старого и нового способов распараллеливания. В новой версии расчеты организованы таким образом, что на каждом ядре прогнозируется движение только одного объекта вместо пучка. На рис. 4, б приведен схематический механизм применяемого способа распределения вычислений. На схеме кругами, как и ранее, обозначены ядра, а квадратами - интегрируемые уравнения движения для объекта с соответствующим номером. Количество запрашиваемых ядер равно N, а объектов - M, причем . Ядра в данной схеме с номерами K1 и K2 . Рис. 5. Распределение ошибки интегрирования при распараллеливании Таблица 2 Время расчетов Число ядер Время, с Новая модель Старая модель 10 24.54 285.69 25 9.29 109.53 50 4.63 40.22 Каждое ядро берет данные только одного объекта из массива. Параметры оставшихся незанятых объектов, если таковые имеются, распределяются следующим образом. Первое ядро K1, которое завершает расчеты для своего объекта, начинает счет для следующего объекта, т.е (N+1)-го, следующее незанятое ядро K2 берет на себя следующий по очереди объект (N+2) и так, пока не закончатся объекты. При данном способе распараллеливания отсутствует зависимость точности от числа выбранных ядер, он дает уверенность, что шаг интегрирования выбирается корректно, исключая ошибки вычисления из-за неправильно подобранного шага для отдельного объекта, а также не допускает увеличение времени интегрирования за счет излишнего дробления шага. Как видно из рис. 5 и табл. 2, при новом способе распараллеливания точность интегрирования более стабильна, а скорость интегрирования возрастает почти на порядок. Заключение Таким образом, в настоящей работе представлена новая версия численной модели движения ИСЗ, состоящая из четырех программных блоков, решающих следующие задачи: - высокоточный прогноз движения ИСЗ; - исследование хаотичности движения околоземных космических объектов, основанное на MEGNO-анализе движения; - определение параметров движения ИСЗ по данным измерений; - исследования резонансной динамики околоземных объектов. Главная особенность новой версии состоит в использование нового более эффективного интегратора, представляющего собой дальнейшее развитие известного интегратора Эверхарта. Приведенные в работе оценки показывают, что при одной и той же точности новый интегратор обладает значительно большим быстродействием по сравнению с предыдущей модификацией того же интегратора. Программный комплекс, предназначенный для использования в среде параллельных вычислений и называемый «Численная модель движения систем ИСЗ», был модифицирован с целью оптимизации процесса распараллеливания вычислений. Как показывают оценки, представленные в работе, при новом способе распараллеливания точность интегрирования более стабильна, а скорость интегрирования возрастает в несколько раз. Замена интегратора в совокупности с новым способ распараллеливания позволяет повысить скорость вычислений в 10 раз при сохранении высокой точности.

Скачать электронную версию публикации