Comments 39
Погрешность в числах с плавающей запятой есть всегда.
Откровенно говоря, тут даже и обсудить нечего, ибо исходников нет.
От того, как написан код, используются ли длинные суммы, используется ли «притягивание» к нулю, выполняется ли суммирование чисел разного порядка, зависит погрешность.
Откровенно говоря, тут даже и обсудить нечего, ибо исходников нет.
От того, как написан код, используются ли длинные суммы, используется ли «притягивание» к нулю, выполняется ли суммирование чисел разного порядка, зависит погрешность.
Основная причина погрешностей скрыта именно в разных порядках переменных. Особенно если начальным условием по одной из переменной задан ноль.
Исходники вряд ли несут большой интерес — там простое решение системы сеточных уравнений явным методом. Переменные не масштабировались, порядки их сохранились, отсюда и эффект.
Исходники вряд ли несут большой интерес — там простое решение системы сеточных уравнений явным методом. Переменные не масштабировались, порядки их сохранились, отсюда и эффект.
Их можно попытаться оптимизировать или хотя бы поискать узкие места.
Так может сменить решатель сетки на что-нибудь более стабильное?
Что-нибудь конечно-элементное? Был бы рад освоить тот же ANSYS, на самом деле. Хотя прежде всего, в силу отсутствия надобности писать весь код самому от и до, а не столько в силу его массовости и работоспособности.
Ну в нашей шарашке на МКЭ как раз этому учили, в целом, мне понравилось, когда одна и та же СЛАУ разными решателями давали весьма такой неплохой разброс решения, как от начальных условий, так и просто от тяготения некоторых к локальным экстремумам.
А в результатае, почти все в дипломах пользовались чьим-нибудь одним проверенным качественным решателем, чтобы не реизобретать велосипеды.
А в результатае, почти все в дипломах пользовались чьим-нибудь одним проверенным качественным решателем, чтобы не реизобретать велосипеды.
Явням методом???
На первый взгляд, кажется, что это связано с тем, что в некоторых режимах работы x87 FPU хранит данные и проводит вычисления не в double precision, а в extended precision. При выгрузке из регистра происходит сужение числа до double. Если же число не покидает регистры FPU, то оно так и остается в extended precision.
Более подробно это описано в en.wikipedia.org/wiki/SSE2#Differences_between_x87_FPU_and_SSE2
Более подробно это описано в en.wikipedia.org/wiki/SSE2#Differences_between_x87_FPU_and_SSE2
При сложении двух чисел с плавающей запятой происходит выравнивание порядков в сторону старшего числа, поэтому мантисса смещается влево, к младшим разрядом и теряет значимость.
Решением проблемы при сложении больших наборов чисел может быть:
— сложение сначала самых маленьких чисел, а затем больших (в идеале по возрастанию).
— блочное сложение (по строкам, подматрицам, ещё как угодно), главное, чтобы числа не отличались сильно.
При различии на 2^23 для float и 2^54 для double меньшее число прибавки не даст.
Ключевые слова: потеря значимости, accuracy problem
Если где-то ошибся, то скажите, пожалуйста.
Решением проблемы при сложении больших наборов чисел может быть:
— сложение сначала самых маленьких чисел, а затем больших (в идеале по возрастанию).
— блочное сложение (по строкам, подматрицам, ещё как угодно), главное, чтобы числа не отличались сильно.
При различии на 2^23 для float и 2^54 для double меньшее число прибавки не даст.
Ключевые слова: потеря значимости, accuracy problem
Если где-то ошибся, то скажите, пожалуйста.
Вот что еще может быть:
1. вы запускаете алгоритм вне допустимых значений параметров (например, для некоторых алгоритмов жидкость должна быть несжимаемой, то есть числа Маха очень маленькими, меньше 0,001, то есть предельная скорость жидкости маленькой; некоторые по тем же или другим причинам требуют малых чисел Рейнольдса/Рэлея/Нуссельта); поэтому решения расходятся.
2. если вы используете генератор случайных чисел (для фазовых превращений, например), то попробуйте использовать более медленный, но лучшего качества
3. решение может расходиться, если разностная схема низкого порядка (например, если считать якобианы/производные не по квадратичной схеме, а по линейной)
Вот в этом комментарии есть множество полезных способов дебага ошибок моделирования гидродинамики.
1. вы запускаете алгоритм вне допустимых значений параметров (например, для некоторых алгоритмов жидкость должна быть несжимаемой, то есть числа Маха очень маленькими, меньше 0,001, то есть предельная скорость жидкости маленькой; некоторые по тем же или другим причинам требуют малых чисел Рейнольдса/Рэлея/Нуссельта); поэтому решения расходятся.
2. если вы используете генератор случайных чисел (для фазовых превращений, например), то попробуйте использовать более медленный, но лучшего качества
3. решение может расходиться, если разностная схема низкого порядка (например, если считать якобианы/производные не по квадратичной схеме, а по линейной)
Вот в этом комментарии есть множество полезных способов дебага ошибок моделирования гидродинамики.
Вот еще несколько идей:
4. некоторые алгоритмы требуют, чтобы результаты расчетов брались не на каждой итерации, а, например, через одну; или как полусумма значений на соседних итерациях; от этого могут быть ошибки
5. могут быть проблемы в параллелизме. Если программа хорошо спроектирована и вы можете заменить параллельную версию на непараллельную быстро (то есть заменить вызовы OpenMP вызовами методов-заглушек (через директиву прекомпиляции), которые ничего не делают), то попробуйте сделать так. Если нет, попробуйте дампить состояние системы на каждой итерации, и смотреть как оно меняется. Может быть, вы заметите что-то странное именно там, где присутствует распараллеливание. Либо менять число процессорных узлов.
4. некоторые алгоритмы требуют, чтобы результаты расчетов брались не на каждой итерации, а, например, через одну; или как полусумма значений на соседних итерациях; от этого могут быть ошибки
5. могут быть проблемы в параллелизме. Если программа хорошо спроектирована и вы можете заменить параллельную версию на непараллельную быстро (то есть заменить вызовы OpenMP вызовами методов-заглушек (через директиву прекомпиляции), которые ничего не делают), то попробуйте сделать так. Если нет, попробуйте дампить состояние системы на каждой итерации, и смотреть как оно меняется. Может быть, вы заметите что-то странное именно там, где присутствует распараллеливание. Либо менять число процессорных узлов.
Спасибо за комментарий со способами, пригодится.
С недопустимыми значениями параметров вопрос актуален, т.к. скорости фильтрации редко превышают мкм/с. Со случайными числами проблем нет, т.к. фазовый переход описан кинетическим уравнением вида dx/dt ~ Cx (вполне типичное уравнение Аррениуса).
Порядок точности схемы, как ни странно, на результат радикально не влиял. Изменение происходило на уровне 3й-4й значащей цифры. Тестировал и первый, и второй порядки. Хотя была нереализованная идея сделать трёхслойную схему по времени, дабы получить второй же порядок в d/dt.
Суммирование результата с разных итераций чаще всего требует определённой неявности, а писать разные громоздкие методы решения таких сеточных систем и отлаживать их довольно противно, простите лентяя.
Параллелизм считать источником ошибки не склонен, поскольку результаты двухпоточного совпадают с вычислениями на одном потоке. А больше, чем на два потока, пока распараллеливать не на чем. Доступа на новый вузовский кластер по статусу не дают (ибо ещё магистр), а старый кластер кафедры, упомянутый в тексте, содержит gcc/gfortran 4.1.2, где OpenMP ещё не было.
С недопустимыми значениями параметров вопрос актуален, т.к. скорости фильтрации редко превышают мкм/с. Со случайными числами проблем нет, т.к. фазовый переход описан кинетическим уравнением вида dx/dt ~ Cx (вполне типичное уравнение Аррениуса).
Порядок точности схемы, как ни странно, на результат радикально не влиял. Изменение происходило на уровне 3й-4й значащей цифры. Тестировал и первый, и второй порядки. Хотя была нереализованная идея сделать трёхслойную схему по времени, дабы получить второй же порядок в d/dt.
Суммирование результата с разных итераций чаще всего требует определённой неявности, а писать разные громоздкие методы решения таких сеточных систем и отлаживать их довольно противно, простите лентяя.
Параллелизм считать источником ошибки не склонен, поскольку результаты двухпоточного совпадают с вычислениями на одном потоке. А больше, чем на два потока, пока распараллеливать не на чем. Доступа на новый вузовский кластер по статусу не дают (ибо ещё магистр), а старый кластер кафедры, упомянутый в тексте, содержит gcc/gfortran 4.1.2, где OpenMP ещё не было.
ну на старых фортранах в нашей шарашке считали на MPI.
openMP разве научился на несколько узлов разноситься?
что-то на моей памяти только несколько узлов с общей памятью
openMP разве научился на несколько узлов разноситься?
что-то на моей памяти только несколько узлов с общей памятью
Кластер из четырёхядерников, так что OpenMP мог бы работать. Только как-то странно жизнеспособность системы поддерживается, народ либо на CUDA переходит, либо на маплы и математики. И с MPI там тоже проблемы. :)
А новый недавно установленный кластер и то, и другое, и третье позволяет, есть возможности, в будущем надеюсь обеспечить аккаунт.
А новый недавно установленный кластер и то, и другое, и третье позволяет, есть возможности, в будущем надеюсь обеспечить аккаунт.
Как пишут в комментарии ниже, скорее всего, вы что-то неправильно закодировали, потому что вы далеко не первый, кто решает подобные задачи. И ошибки должны появляться намного позже шестого знака после запятой (в бенчмарках типа потока Пуазейля, по крайней мере).
И как пишет небезызвестный Стив Макконнел в книжке «Совершенный код», если вы видите следы копыт, не думайте что это зебра (в смысле «ошибка компилятора/базы данных/реализации стандартов/и т.п.»).
Кстати, вот что вы для начала можете рассчитать: бенчмарки а ля поток Пуазейля, течение Куэтта, конвекция Рэлея-Бенара и т.п.
И как пишет небезызвестный Стив Макконнел в книжке «Совершенный код», если вы видите следы копыт, не думайте что это зебра (в смысле «ошибка компилятора/базы данных/реализации стандартов/и т.п.»).
Кстати, вот что вы для начала можете рассчитать: бенчмарки а ля поток Пуазейля, течение Куэтта, конвекция Рэлея-Бенара и т.п.
Не стоит забывать, что течения Пуазейля и Куэтта, Рэлей-Бенар и прочее — это совершенно другие системы. В них: Навье-Стокс + несжимаемость + уравнение теплопроводности (если конвекция).
Здесь: закон Дарси + уравнения неразрывности + теплопроводность. Система совершенно иная. Свойства течений в пористой среде во многом нетривиальны, особенно если есть намёки на конвекцию и теплоперенос. Потому упомянутые тесты не сгодятся, хотя они и хороши, прежде всего возможностью хотя бы какой-то аналитики.
Здесь: закон Дарси + уравнения неразрывности + теплопроводность. Система совершенно иная. Свойства течений в пористой среде во многом нетривиальны, особенно если есть намёки на конвекцию и теплоперенос. Потому упомянутые тесты не сгодятся, хотя они и хороши, прежде всего возможностью хотя бы какой-то аналитики.
Я для примера привел) Как ни странно, как раз полгода назад я начал делать PhD по близкой теме, для продолжения работы своих коллег. Но они пользовались методом решеток Больцмана для моделирования гидродинамики в пористой среде, без всяких приближений типа закона Дарси.
устойчивость явной схемы не зависит ни от начальных условий, ни от правой части.
Только если эта явная схема — линейная. А бывают еще и нелинейные схемы, которые возникают при решении нелинейных уравнений.
Если Вас не затруднит, укажите, пожалуйста, пример явной схемы, устойчивость которой зависит от начальных условий.
Смотря что понимать под устойчивостью. Если устойчивость в том смысле, что схема не расходилась или расходилась бы по прошествии некоторого времени из-за экспоненциального роста ошибок, то начальные условия здесь ни при чём.
Если в начальном условии гиперболической задачи есть разрыв (например, расчёт ударной волны), и в схеме слабая диффузия, то эффект вполне известен — вместо красивого ровного фронта может получиться нечто колеблющееся и на ударную волну не похожее (немонотонность схемы).
Стоит сказать, что борьба с немонотонностью в том числе и в описываемой здесь задаче тоже заняла немалое место. В схеме болтается парочка регуляризаторов, снижающих точность и «размазывающих» тем самым колебания.
Если в начальном условии гиперболической задачи есть разрыв (например, расчёт ударной волны), и в схеме слабая диффузия, то эффект вполне известен — вместо красивого ровного фронта может получиться нечто колеблющееся и на ударную волну не похожее (немонотонность схемы).
Стоит сказать, что борьба с немонотонностью в том числе и в описываемой здесь задаче тоже заняла немалое место. В схеме болтается парочка регуляризаторов, снижающих точность и «размазывающих» тем самым колебания.
Вроде бы существует определение устойчивости разностной схемы :)
По поводу немонотонности. Вот тут уже есть разные определения монотонности. Например, в последнее время в задачах с фронтами популярны монотонные по Хартену схемы (TVD-схемы).
В расчете же ударной волны критично также условие Куранта, которое является там более жестким, чем условие устойчивости явной схемы.
По поводу немонотонности. Вот тут уже есть разные определения монотонности. Например, в последнее время в задачах с фронтами популярны монотонные по Хартену схемы (TVD-схемы).
В расчете же ударной волны критично также условие Куранта, которое является там более жестким, чем условие устойчивости явной схемы.
Схема решения уравнения 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 отрицательно, то возмущения в схеме вообще распространяются в другую сторону и схема тем более не может быть устойчивой.
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 отрицательно, то возмущения в схеме вообще распространяются в другую сторону и схема тем более не может быть устойчивой.
автору статьи. скажите, а что у Вас на графиках отложено? И что считали, вода-воздух в пористой среде? с диффузией? Явная схема по насыщенности, неявная по давлению?
на графиках — насыщенность водой в зависимости от координаты. Задача одномерная, ибо это стадия тестирования модели.
Верно, вода и идеальный газ в пористой среде, фильтрация по закону Дарси, но без учёта диффузии как отдельных компонент в поры, так и друг в друга. В ближайшее время, возможно, эти эффекты будут учтены — после того, как свои задачи потестирует другой участник проекта, примем окончательное решение о том, как вести описание.
Схема полностью явная по всем переменным. И, т.к. уравнения содержат произведения величин вида d(u1 u2) / dt, то они расписывались в форме u2 d(u1)/dt + u1 d(u2)/dt, образуя небольшую системку для производных уже от одной переменной. Решение её даёт сразу полный набор приращений всех переменных на текущем шаге, но здесь опять же происходит столкновение с совместной обработкой величин разных порядков. Вполне возможно, это и есть слабое место алгоритма.
Верно, вода и идеальный газ в пористой среде, фильтрация по закону Дарси, но без учёта диффузии как отдельных компонент в поры, так и друг в друга. В ближайшее время, возможно, эти эффекты будут учтены — после того, как свои задачи потестирует другой участник проекта, примем окончательное решение о том, как вести описание.
Схема полностью явная по всем переменным. И, т.к. уравнения содержат произведения величин вида d(u1 u2) / dt, то они расписывались в форме u2 d(u1)/dt + u1 d(u2)/dt, образуя небольшую системку для производных уже от одной переменной. Решение её даёт сразу полный набор приращений всех переменных на текущем шаге, но здесь опять же происходит столкновение с совместной обработкой величин разных порядков. Вполне возможно, это и есть слабое место алгоритма.
В чем причина трудно сказать не имея перед глазами расчетных формул и их программной реализации :)
P.S. вы случайно не пальцеобразование исследуете?
P.S. вы случайно не пальцеобразование исследуете?
Причина и её устранение это дело наживное. Чаще действую по принципу «Работает — не трогай».
Пальцеобразование это там где что-то подобное неустойчивости Рэлея-Тейлора, но в пористой среде? Условно говоря, капли, тянущие за собой струйку жидкости?
Нет, на самом деле, работа идёт по гидратам метана. Достаточно продолжительный уже проект. Или проектик. Скорее второе, масштабы деятельности невелики. :)
Пальцеобразование это там где что-то подобное неустойчивости Рэлея-Тейлора, но в пористой среде? Условно говоря, капли, тянущие за собой струйку жидкости?
Нет, на самом деле, работа идёт по гидратам метана. Достаточно продолжительный уже проект. Или проектик. Скорее второе, масштабы деятельности невелики. :)
Насколько я понимаю, зависит еще и от того, как эта схема получена. Уравнения гидродинамики (Навье-Стокса) нелинейны сами по себе, и перевод их в явную конечно-разностную схему--нетривиальная задача. То есть неустойчивость появляется не на этапе решения схемы, а на этапе перехода к этой схеме.
Например, вычислительная схема метода решеток Больцмана полностью явная, но устойчивость сильно зависит от начальных условий.
Даже если моделировать простой поток Пуазейля (по трубе) со слишком большим градиентом давления--решение расходится.
Например, вычислительная схема метода решеток Больцмана полностью явная, но устойчивость сильно зависит от начальных условий.
Даже если моделировать простой поток Пуазейля (по трубе) со слишком большим градиентом давления--решение расходится.
Как раз перевод в явную схему обычно сложностей не представляет. Всё, что вне производной по времени — с текущего шага, а в ней (при простейшей аппроксимации) — в числителе возникнет всего одна лишь скорость с нового шага. Вот если схема полу-неявная, тогда есть с чем повозиться. Какое слагаемое отправить на новый шаг, какое на текущий, а может и на предыдущий (в трёхслойной схеме), это уже вопрос действительно нетривиальный.
Чтобы этого избежать, конечно, можно решать полностью неявную схему. Но для такого требуется писать сеточный метод Ньютона или ему подобный для нелинейных систем. Более жуткой деятельности при работе с разностными задачами представить трудно.
А течение Пуазейля с большим градиентом давления — не получаются ли там числа Рейнольдса такими, что начинается турбулентность и неустойчивость самого течения?
Чтобы этого избежать, конечно, можно решать полностью неявную схему. Но для такого требуется писать сеточный метод Ньютона или ему подобный для нелинейных систем. Более жуткой деятельности при работе с разностными задачами представить трудно.
А течение Пуазейля с большим градиентом давления — не получаются ли там числа Рейнольдса такими, что начинается турбулентность и неустойчивость самого течения?
«Как раз перевод в явную схему обычно сложностей не представляет.» Уверяю вас, вы очень сильно заблуждаетесь. Для простейших уравнений--не представляет. Для Навье-Стокса, Фоккера-Планка, Больцмана--представляет. Даже простое уравнение диффузии вы нормально не решите при сложных граничных условиях; поэтому приходится пользоваться вот этим методом: en.wikipedia.org/wiki/Single_particle_tracking.
«начинается турбулентность и неустойчивость самого течения» В том-то и фокус, что до турбулентности еще далеко, а схема начинает расходиться.
«начинается турбулентность и неустойчивость самого течения» В том-то и фокус, что до турбулентности еще далеко, а схема начинает расходиться.
А метод частиц всегда интересовал и интриговал красивыми картинками в литературе, вот что вспомнилось-то. В частности, «Вычислительные методы в физике» Поттера, главы о гидродинамике.
Признаю, отсутствие опыта симуляции именно течений и динамики, сказывается, и во многом связь здесь со спецификой кафедры — основным направлением является конвективная устойчивость, где расчёт течений как таковых не особенно важен. По крайней мере, могу уверенно сказать, что значительная часть исследований ведётся на простых задачах, которые порой допускают и аналитику для равновесных и стационарных решений. Работаем в слоях (горизонтальных-вертикальных), в полостях несложных форм (куб-параллелепипед), а то и вовсе в ячейках Хеле-Шоу. Про сферическую полость ходит легенда, что её проклял тот, кто исследовал, и с тех пор больше никто за сферы не брался. :)
Признаю, отсутствие опыта симуляции именно течений и динамики, сказывается, и во многом связь здесь со спецификой кафедры — основным направлением является конвективная устойчивость, где расчёт течений как таковых не особенно важен. По крайней мере, могу уверенно сказать, что значительная часть исследований ведётся на простых задачах, которые порой допускают и аналитику для равновесных и стационарных решений. Работаем в слоях (горизонтальных-вертикальных), в полостях несложных форм (куб-параллелепипед), а то и вовсе в ячейках Хеле-Шоу. Про сферическую полость ходит легенда, что её проклял тот, кто исследовал, и с тех пор больше никто за сферы не брался. :)
Скорее всего вы что-то неправильно закодировали. Подобные задачи считают тысячи приматов каждый год, но мало кому удается нарваться на проблемы с машинной точностью.
Да уж, в точку! Я как-то долго ломал голову, почему у меня черт знает что отображается в окне openGL после кое-каких расчетов. Оказалось, что из-за использования float, а не double. Но т.к. переходить на double не представлялось возможным, пришлось перегруппировать вычисления так, чтобы погрешность в результате была бы минимальной.
Для операций с плавающей точкой!
1. Погрешность преобразования из одной системы в другую (например из десятичной в двоичную).
2. Нарушение ассоциативности операций сложения (или вычитания) за счет округления. Часто
(a + b) + c != a + (b + c)
поэтому может помочь «правильное» переписывание формул по которым считаете. С умножением все тип-топ. Я понимаю что у всех в башке матанализ, а должна быть алгебра. Но это исторический дефект нашего математического образования.
Читайте Кнута «Искусство программирования» т.2 «Получисленные алгоритмы».
3. Обычно меняют само представление. Их обычно стандартно: обычной точности, удвоенной и long double и смотрят на результаты, затем делают выводы.
4. Советую уйти с фортана на scipy, он его все равно вызывает и еще C. Легче будет писать программы и экспериментировать с точностью.
1. Погрешность преобразования из одной системы в другую (например из десятичной в двоичную).
2. Нарушение ассоциативности операций сложения (или вычитания) за счет округления. Часто
(a + b) + c != a + (b + c)
поэтому может помочь «правильное» переписывание формул по которым считаете. С умножением все тип-топ. Я понимаю что у всех в башке матанализ, а должна быть алгебра. Но это исторический дефект нашего математического образования.
Читайте Кнута «Искусство программирования» т.2 «Получисленные алгоритмы».
3. Обычно меняют само представление. Их обычно стандартно: обычной точности, удвоенной и long double и смотрят на результаты, затем делают выводы.
4. Советую уйти с фортана на scipy, он его все равно вызывает и еще C. Легче будет писать программы и экспериментировать с точностью.
А может что-то не понял, но почему нельзя тот же -ftz использовать? И брать -fp-model fast для численных расчетов ИМХО как-то неправильно.
-fp-model fast был чисто для проверки взят, насколько плохо всё будет в таком режиме. Результат виден. Кстати говоря, gfortran с -ffast-math отработал лишь немногим хуже, чем без этого ключа только на одной оптимизации — также получился достаточно гладкий график с небольшими колебаниями.
С -ftz попытки были, тоже малоуспешные.
С -ftz попытки были, тоже малоуспешные.
Sign up to leave a comment.
Ошибки вычислений в окрестностях машинного нуля