Ошибки вычислений в окрестностях машинного нуля

Периодически на хабре возникают замечательные статьи о тонкостях арифметики с плавающей точкой. Собственно, упомянутая публикация стала одним из первых источников, прочитанных при попытках разобраться с проблемой. Яснее от этого сразу не стало, но тем не менее, организация нейронных связей как-то упорядочилась. Ближе к делу.

Проблема образовалась при выполнении расчётов в рамках одного проекта и будущей магистерской диссертации по гидродинамике пористой среды. Не скрою, что корни скрыты отчасти в личной криворукости автора и пренебрежении банальными общеизвестными советами касательно обработки малых чисел, но тем не менее, это привело к достаточно интересным наблюдениям и размышлениям.

Если отбросить физическую суть, то задача состоит в необходимости решения системы семи уравнений в частных производных. Сказано-сделано, сложности в этом особой нет — пишем явную конечно-разностную схему, распараллеливаем на OpenMP и после окончательной синтаксической отладки и оптимизации скорости «машинка начинает шуршать».

Вычислительная конфигурация: видавший виды HP 550 с Core 2 Duo 1.8 ГГц на борту, под управлением Ubuntu 11.04.

Компиляторы: gfortran 4.5.2 и Intel Fortran Compiler 12.1.0.

В начальных условиях предполагается, что внутри расчётной области совсем нет воды в жидком состоянии — она появляется в процессе фазовых превращений. И именно это отсутствие сыграло главную роль в нашем спектакле.

Итак, воды в области нет. Что же хочется написать в начальном условии? Естественно, сразу же было написано 0.0D+00 (программа написана с двойной точностью вещественных чисел). Счёт на первоначальных этапах идёт в непосредственной окрестности машинного нуля. Каковы результаты? Посмотрим на график:

Ноль

Здесь изображено распределение насыщенности пор водой вдоль координаты в наиболее интересной части расчётной области по прошествии чуть более суток физического времени модели.

Нужное отступление о подписях в легенде графика: всего было опробовано 18 различных комбинаций ключей (6 для gfortran и 12 для ifort), однако многие из них давали абсолютно одинаковые результаты, и потому объединены. Квадратные скобки в легенде означают, что могла быть написана любая из заключённых в них опций. Например, «шифровка» -O[1,2,3] [-fp-model [strict] [precise]] говорит о том, что компилятор использовался с оптимизацией всех возможных уровней, и дополнительно могла быть включена одна моделей вычислений с плавающей точкой (а могла быть и не включена). Три варианта (два от -fp-model и один без неё) умножить на три уровня оптимизации — итого девять комбинаций. Все они оказались эквивалентны.

А теперь результат. Нечто реалистичное и физически возможное удалось получить только на gfortran без включения соответствия стандарту IEEE 754 (ключ -ffloat-store). Весь остальной хаос линий не содержит ни капельки физического смысла, потому что даже математически уравнения этого не допускают. Изначально подозревавшаяся неустойчивость разностной схемы была оправдана, поскольку никакие методы борьбы с ней к успеху не привели.

Было замечено, что при наличии воды в начальный момент счёт остаётся устойчивым и различия между опциями и компиляторами скрываются в окрестности шестой значащей цифры. И т.к. характерные порядки величины, судя по полученным графикам, должны быть в районе 1.0e-2, то в начальное условие было вписано некоторое ненулевое значение, но очень маленькое. Подбором удалось установить, что для 8-байтового real оно должно быть не менее 1.0e-21. И тогда:

Почти ноль

Да, здесь, как и написано в легенде, на самом деле пять графиков. Просто мы получаем то самое различие в пределах шестой значащей цифры.

Причины? Вполне очевидны. Они лежат в тонкостях одновременной обработки больших и малых чисел. И тот факт, что явление настолько существенно по своим масштабам, был, в общем-то, ожидаемым. Но, прежде всего, интерес вызывает нестабильность работы действительно отличного инструмента от Intel на фоне относительного успеха gfortran 4.5.2. Стоит отметить, что подобные проблемы были также найдены при расчётах на старой версии gfortran 4.1.2 с включённой оптимизацией и без других опций, установленной на доступном кластере (под управлением Slamd64), однако им тогда не было уделено должного внимания.

Соответствие IEEE 754, как ни странно, сыграло критическую роль для gfortran. Без него счёт достаточно стабильный и точный. Для детища Intel это оказалось не столь существенно, ибо оно и так не работало корректно.

Итак, выводы и мысли о причинах увиденного.

  • Наиболее вероятным кандидатом на роль причины наблюдаемого поведения расчётов представляются тонкости округления чисел. Т.к. в начальном распределении величины задано нулевое значение, то на первых шагах счёт производится практически на границе машинной точности. Соответственно, это приводит к накоплению заметных погрешностей, которые и проявляются в конечном результате.
  • Потери точности в алгоритме вызываются, очевидно, и тем фактом, что решается размерная система семи дифференциальных уравнений в частных производных, переменные в каждом из которых имеют собственные характерные значения, существенно отличающиеся от остальных. При правильном выборе масштабных множителей, в безразмерной системе уравнений получить близкие хотя бы по порядку величины значения всех переменных, хотя при первоначальных попытках провести обезразмеривание в системе возникали малые коэффициенты перед производными, что и стало поводом для отказа от данной процедуры.
  • Вопрос же о том, почему gfortran, не соответствующий стандарту вычислений с плавающей точкой, способен выдавать приемлемый результат, остаётся открытым. Разумно предполагать наличие каких-то собственных, отличных от стандарта, правил округления, которые и обеспечивают сохранение стабильного счёта, а также их корректировку и уточнение в процессе развития компилятора. «Хрупкий баланс ошибок» либо продуманное исправление подхода? Увы, на моём уровне знаний об инструментарии, нацеленном в первую очередь именно на применение компиляторов, а не их тестирование и изучение свойств, это неизвестно. Но заставляет задуматься и вспомнить предостережения о возможных потерях точности в тех или иных местах программ, данные ещё на первых этапах обучения численным методам в вузе.


Тестирование компиляторов на соответствие IEEE 754 проводилось с помощью FORTRAN-версии «Floating point paranoia» Уильяма Кэхэна.
Ads
AdBlock has stolen the banner, but banners are not teeth — they will be back

More

Comments 39

    +11
    Погрешность в числах с плавающей запятой есть всегда.

    Откровенно говоря, тут даже и обсудить нечего, ибо исходников нет.
    От того, как написан код, используются ли длинные суммы, используется ли «притягивание» к нулю, выполняется ли суммирование чисел разного порядка, зависит погрешность.
      +1
      Основная причина погрешностей скрыта именно в разных порядках переменных. Особенно если начальным условием по одной из переменной задан ноль.

      Исходники вряд ли несут большой интерес — там простое решение системы сеточных уравнений явным методом. Переменные не масштабировались, порядки их сохранились, отсюда и эффект.
        0
        Их можно попытаться оптимизировать или хотя бы поискать узкие места.
          0
          Так может сменить решатель сетки на что-нибудь более стабильное?
            0
            Что-нибудь конечно-элементное? Был бы рад освоить тот же ANSYS, на самом деле. Хотя прежде всего, в силу отсутствия надобности писать весь код самому от и до, а не столько в силу его массовости и работоспособности.
              0
              Ну в нашей шарашке на МКЭ как раз этому учили, в целом, мне понравилось, когда одна и та же СЛАУ разными решателями давали весьма такой неплохой разброс решения, как от начальных условий, так и просто от тяготения некоторых к локальным экстремумам.
              А в результатае, почти все в дипломах пользовались чьим-нибудь одним проверенным качественным решателем, чтобы не реизобретать велосипеды.
            0
            Явням методом???
          +11
          На первый взгляд, кажется, что это связано с тем, что в некоторых режимах работы x87 FPU хранит данные и проводит вычисления не в double precision, а в extended precision. При выгрузке из регистра происходит сужение числа до double. Если же число не покидает регистры FPU, то оно так и остается в extended precision.
          Более подробно это описано в en.wikipedia.org/wiki/SSE2#Differences_between_x87_FPU_and_SSE2
            +3
            При сложении двух чисел с плавающей запятой происходит выравнивание порядков в сторону старшего числа, поэтому мантисса смещается влево, к младшим разрядом и теряет значимость.

            Решением проблемы при сложении больших наборов чисел может быть:
            — сложение сначала самых маленьких чисел, а затем больших (в идеале по возрастанию).
            — блочное сложение (по строкам, подматрицам, ещё как угодно), главное, чтобы числа не отличались сильно.

            При различии на 2^23 для float и 2^54 для double меньшее число прибавки не даст.

            Ключевые слова: потеря значимости, accuracy problem

            Если где-то ошибся, то скажите, пожалуйста.
              0
              Добавлю, что при «удалении» числа от нуля увеличивается шаг между числами, что автоматически сказывается на потери точности. Поэтому при вычислениях, в которых участвуют числа одного порядка, иногда полезно их переносить в окрестность нуля.
              +1
              Вот что еще может быть:
              1. вы запускаете алгоритм вне допустимых значений параметров (например, для некоторых алгоритмов жидкость должна быть несжимаемой, то есть числа Маха очень маленькими, меньше 0,001, то есть предельная скорость жидкости маленькой; некоторые по тем же или другим причинам требуют малых чисел Рейнольдса/Рэлея/Нуссельта); поэтому решения расходятся.
              2. если вы используете генератор случайных чисел (для фазовых превращений, например), то попробуйте использовать более медленный, но лучшего качества
              3. решение может расходиться, если разностная схема низкого порядка (например, если считать якобианы/производные не по квадратичной схеме, а по линейной)

              Вот в этом комментарии есть множество полезных способов дебага ошибок моделирования гидродинамики.
                0
                Вот еще несколько идей:
                4. некоторые алгоритмы требуют, чтобы результаты расчетов брались не на каждой итерации, а, например, через одну; или как полусумма значений на соседних итерациях; от этого могут быть ошибки
                5. могут быть проблемы в параллелизме. Если программа хорошо спроектирована и вы можете заменить параллельную версию на непараллельную быстро (то есть заменить вызовы OpenMP вызовами методов-заглушек (через директиву прекомпиляции), которые ничего не делают), то попробуйте сделать так. Если нет, попробуйте дампить состояние системы на каждой итерации, и смотреть как оно меняется. Может быть, вы заметите что-то странное именно там, где присутствует распараллеливание. Либо менять число процессорных узлов.
                  0
                  Спасибо за комментарий со способами, пригодится.

                  С недопустимыми значениями параметров вопрос актуален, т.к. скорости фильтрации редко превышают мкм/с. Со случайными числами проблем нет, т.к. фазовый переход описан кинетическим уравнением вида dx/dt ~ Cx (вполне типичное уравнение Аррениуса).

                  Порядок точности схемы, как ни странно, на результат радикально не влиял. Изменение происходило на уровне 3й-4й значащей цифры. Тестировал и первый, и второй порядки. Хотя была нереализованная идея сделать трёхслойную схему по времени, дабы получить второй же порядок в d/dt.

                  Суммирование результата с разных итераций чаще всего требует определённой неявности, а писать разные громоздкие методы решения таких сеточных систем и отлаживать их довольно противно, простите лентяя.

                  Параллелизм считать источником ошибки не склонен, поскольку результаты двухпоточного совпадают с вычислениями на одном потоке. А больше, чем на два потока, пока распараллеливать не на чем. Доступа на новый вузовский кластер по статусу не дают (ибо ещё магистр), а старый кластер кафедры, упомянутый в тексте, содержит gcc/gfortran 4.1.2, где OpenMP ещё не было.
                    0
                    ну на старых фортранах в нашей шарашке считали на MPI.
                    openMP разве научился на несколько узлов разноситься?
                    что-то на моей памяти только несколько узлов с общей памятью
                      0
                      Кластер из четырёхядерников, так что OpenMP мог бы работать. Только как-то странно жизнеспособность системы поддерживается, народ либо на CUDA переходит, либо на маплы и математики. И с MPI там тоже проблемы. :)

                      А новый недавно установленный кластер и то, и другое, и третье позволяет, есть возможности, в будущем надеюсь обеспечить аккаунт.
                      0
                      Как пишут в комментарии ниже, скорее всего, вы что-то неправильно закодировали, потому что вы далеко не первый, кто решает подобные задачи. И ошибки должны появляться намного позже шестого знака после запятой (в бенчмарках типа потока Пуазейля, по крайней мере).

                      И как пишет небезызвестный Стив Макконнел в книжке «Совершенный код», если вы видите следы копыт, не думайте что это зебра (в смысле «ошибка компилятора/базы данных/реализации стандартов/и т.п.»).

                      Кстати, вот что вы для начала можете рассчитать: бенчмарки а ля поток Пуазейля, течение Куэтта, конвекция Рэлея-Бенара и т.п.
                        0
                        Не стоит забывать, что течения Пуазейля и Куэтта, Рэлей-Бенар и прочее — это совершенно другие системы. В них: Навье-Стокс + несжимаемость + уравнение теплопроводности (если конвекция).

                        Здесь: закон Дарси + уравнения неразрывности + теплопроводность. Система совершенно иная. Свойства течений в пористой среде во многом нетривиальны, особенно если есть намёки на конвекцию и теплоперенос. Потому упомянутые тесты не сгодятся, хотя они и хороши, прежде всего возможностью хотя бы какой-то аналитики.
                          0
                          Я для примера привел) Как ни странно, как раз полгода назад я начал делать PhD по близкой теме, для продолжения работы своих коллег. Но они пользовались методом решеток Больцмана для моделирования гидродинамики в пористой среде, без всяких приближений типа закона Дарси.
                  0
                  устойчивость явной схемы не зависит ни от начальных условий, ни от правой части.
                    0
                    Только если эта явная схема — линейная. А бывают еще и нелинейные схемы, которые возникают при решении нелинейных уравнений.
                      0
                      Если Вас не затруднит, укажите, пожалуйста, пример явной схемы, устойчивость которой зависит от начальных условий.
                        0
                        Смотря что понимать под устойчивостью. Если устойчивость в том смысле, что схема не расходилась или расходилась бы по прошествии некоторого времени из-за экспоненциального роста ошибок, то начальные условия здесь ни при чём.

                        Если в начальном условии гиперболической задачи есть разрыв (например, расчёт ударной волны), и в схеме слабая диффузия, то эффект вполне известен — вместо красивого ровного фронта может получиться нечто колеблющееся и на ударную волну не похожее (немонотонность схемы).

                        Стоит сказать, что борьба с немонотонностью в том числе и в описываемой здесь задаче тоже заняла немалое место. В схеме болтается парочка регуляризаторов, снижающих точность и «размазывающих» тем самым колебания.
                          0
                          Вроде бы существует определение устойчивости разностной схемы :)

                          По поводу немонотонности. Вот тут уже есть разные определения монотонности. Например, в последнее время в задачах с фронтами популярны монотонные по Хартену схемы (TVD-схемы).

                          В расчете же ударной волны критично также условие Куранта, которое является там более жестким, чем условие устойчивости явной схемы.
                          +1
                          Схема решения уравнения u_t=-u*u_x.
                          u_{t+1,i}=u_{t,i}-\tau * u_{t,i}*(u_{t,i}-u{t,i-1})/h
                          Сходимость этой схемы зависит от соотношения между h/\tau и u. Если h/\tau меньше чем max(u), то возмущения в схеме распространяются медленнее чем в уравнении, а значит схема не может быть устойчивой. Если u отрицательно, то возмущения в схеме вообще распространяются в другую сторону и схема тем более не может быть устойчивой.
                            +1
                            Вы только что сформулировали условие Куранта.
                              0
                              Знаю. Просто читателям может быть лень лезть в википедию и читать что это такое — а так глядишь кто-нибудь его запомнит.
                        0
                        автору статьи. скажите, а что у Вас на графиках отложено? И что считали, вода-воздух в пористой среде? с диффузией? Явная схема по насыщенности, неявная по давлению?
                          0
                          на графиках — насыщенность водой в зависимости от координаты. Задача одномерная, ибо это стадия тестирования модели.

                          Верно, вода и идеальный газ в пористой среде, фильтрация по закону Дарси, но без учёта диффузии как отдельных компонент в поры, так и друг в друга. В ближайшее время, возможно, эти эффекты будут учтены — после того, как свои задачи потестирует другой участник проекта, примем окончательное решение о том, как вести описание.

                          Схема полностью явная по всем переменным. И, т.к. уравнения содержат произведения величин вида d(u1 u2) / dt, то они расписывались в форме u2 d(u1)/dt + u1 d(u2)/dt, образуя небольшую системку для производных уже от одной переменной. Решение её даёт сразу полный набор приращений всех переменных на текущем шаге, но здесь опять же происходит столкновение с совместной обработкой величин разных порядков. Вполне возможно, это и есть слабое место алгоритма.
                            0
                            В чем причина трудно сказать не имея перед глазами расчетных формул и их программной реализации :)

                            P.S. вы случайно не пальцеобразование исследуете?
                              0
                              Причина и её устранение это дело наживное. Чаще действую по принципу «Работает — не трогай».

                              Пальцеобразование это там где что-то подобное неустойчивости Рэлея-Тейлора, но в пористой среде? Условно говоря, капли, тянущие за собой струйку жидкости?

                              Нет, на самом деле, работа идёт по гидратам метана. Достаточно продолжительный уже проект. Или проектик. Скорее второе, масштабы деятельности невелики. :)
                          0
                          Насколько я понимаю, зависит еще и от того, как эта схема получена. Уравнения гидродинамики (Навье-Стокса) нелинейны сами по себе, и перевод их в явную конечно-разностную схему--нетривиальная задача. То есть неустойчивость появляется не на этапе решения схемы, а на этапе перехода к этой схеме.

                          Например, вычислительная схема метода решеток Больцмана полностью явная, но устойчивость сильно зависит от начальных условий.

                          Даже если моделировать простой поток Пуазейля (по трубе) со слишком большим градиентом давления--решение расходится.
                            0
                            Как раз перевод в явную схему обычно сложностей не представляет. Всё, что вне производной по времени — с текущего шага, а в ней (при простейшей аппроксимации) — в числителе возникнет всего одна лишь скорость с нового шага. Вот если схема полу-неявная, тогда есть с чем повозиться. Какое слагаемое отправить на новый шаг, какое на текущий, а может и на предыдущий (в трёхслойной схеме), это уже вопрос действительно нетривиальный.

                            Чтобы этого избежать, конечно, можно решать полностью неявную схему. Но для такого требуется писать сеточный метод Ньютона или ему подобный для нелинейных систем. Более жуткой деятельности при работе с разностными задачами представить трудно.

                            А течение Пуазейля с большим градиентом давления — не получаются ли там числа Рейнольдса такими, что начинается турбулентность и неустойчивость самого течения?
                              +1
                              «Как раз перевод в явную схему обычно сложностей не представляет.» Уверяю вас, вы очень сильно заблуждаетесь. Для простейших уравнений--не представляет. Для Навье-Стокса, Фоккера-Планка, Больцмана--представляет. Даже простое уравнение диффузии вы нормально не решите при сложных граничных условиях; поэтому приходится пользоваться вот этим методом: en.wikipedia.org/wiki/Single_particle_tracking.

                              «начинается турбулентность и неустойчивость самого течения» В том-то и фокус, что до турбулентности еще далеко, а схема начинает расходиться.
                                +1
                                А метод частиц всегда интересовал и интриговал красивыми картинками в литературе, вот что вспомнилось-то. В частности, «Вычислительные методы в физике» Поттера, главы о гидродинамике.

                                Признаю, отсутствие опыта симуляции именно течений и динамики, сказывается, и во многом связь здесь со спецификой кафедры — основным направлением является конвективная устойчивость, где расчёт течений как таковых не особенно важен. По крайней мере, могу уверенно сказать, что значительная часть исследований ведётся на простых задачах, которые порой допускают и аналитику для равновесных и стационарных решений. Работаем в слоях (горизонтальных-вертикальных), в полостях несложных форм (куб-параллелепипед), а то и вовсе в ячейках Хеле-Шоу. Про сферическую полость ходит легенда, что её проклял тот, кто исследовал, и с тех пор больше никто за сферы не брался. :)
                          +2
                          Скорее всего вы что-то неправильно закодировали. Подобные задачи считают тысячи приматов каждый год, но мало кому удается нарваться на проблемы с машинной точностью.
                            0
                            Да уж, в точку! Я как-то долго ломал голову, почему у меня черт знает что отображается в окне openGL после кое-каких расчетов. Оказалось, что из-за использования float, а не double. Но т.к. переходить на double не представлялось возможным, пришлось перегруппировать вычисления так, чтобы погрешность в результате была бы минимальной.
                              +3
                              Для операций с плавающей точкой!

                              1. Погрешность преобразования из одной системы в другую (например из десятичной в двоичную).

                              2. Нарушение ассоциативности операций сложения (или вычитания) за счет округления. Часто
                              (a + b) + c != a + (b + c)
                              поэтому может помочь «правильное» переписывание формул по которым считаете. С умножением все тип-топ. Я понимаю что у всех в башке матанализ, а должна быть алгебра. Но это исторический дефект нашего математического образования.

                              Читайте Кнута «Искусство программирования» т.2 «Получисленные алгоритмы».

                              3. Обычно меняют само представление. Их обычно стандартно: обычной точности, удвоенной и long double и смотрят на результаты, затем делают выводы.

                              4. Советую уйти с фортана на scipy, он его все равно вызывает и еще C. Легче будет писать программы и экспериментировать с точностью.
                                0
                                А может что-то не понял, но почему нельзя тот же -ftz использовать? И брать -fp-model fast для численных расчетов ИМХО как-то неправильно.
                                  0
                                  -fp-model fast был чисто для проверки взят, насколько плохо всё будет в таком режиме. Результат виден. Кстати говоря, gfortran с -ffast-math отработал лишь немногим хуже, чем без этого ключа только на одной оптимизации — также получился достаточно гладкий график с небольшими колебаниями.

                                  С -ftz попытки были, тоже малоуспешные.

                                Only users with full accounts can post comments. Log in, please.