Введение

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

Исходная задача

Как и в прошлый раз, в основу статьи лягут несколько задач на моделирование из программы курса «Цифровизация физических процессов» на факультете ВШПИ МФТИ.

М1 Численно вычислить интеграл Френеля для симметричных отверстий различной формы (эллипс, квадрат, прямоугольник, ромб – по выбору) для точки, расположенной на оси симметрии отверстия на различных расстояниях от отверстия. Выбрать оптимальный способ расчета исходя из формы отверстия.

М2. Вычислить значение интеграла Френеля для точек, смещённых от оси круглого отверстия.

М3. Численно вычислить дифракционную картину по схеме Юнга исходя из спектра падающего света, углового размера источника, геометрических размеров щелей. Предусмотреть возможность добавить входную щель.

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

Теория

Иллюстрация к формуле
Иллюстрация к формуле

Для понимания того, что мы будем дальше делать, стоит немного поговорить о теории. Согласно теории дифракции Рэлея-Зоммерфельда, дифракционная картина электрического поля в точке (x, y, z) определяется следующим решением уравнения Гельмгольца:

E(x,y,z) = \frac{1}{i\lambda} \int{\int_{-\infty}^{+\infty}{E(x',y',0)\frac{e^{ikr}z}{r^{2}}(1+\frac{i}{kr})dx'}dy'}   (1),

где

E(x', y', 0) - значения поля в точке плоскости "с которой светим",
λ‎ - длина волны,

k = \frac{2\pi}{\lambda}, r = \sqrt{(x-x')^{2} + (y-y')^{2} + z^{2}}

Пусть

\rho^{2} = (x-x')^{2}+(y-y')^2

Применим аппроксимацию Френеля:

r = \sqrt{(x-x')^{2} + (y-y')^{2} + z^{2}} =\sqrt{\rho^{2}+z^{2}} = z\sqrt{1+\frac{\rho^{2}}{z^{2}}} = z(1 + \frac{\rho^{2}}{2z^{2}} + ...) \approx\approx z + \frac{\rho^{2}}{2z}

Положим также λ << ρ << z , тогда из формулы (1) получим:

E(x,y,z) = \frac{e^{ikz}}{i\lambda z} \int{\int_{-\infty}^{+\infty}{E(x',y',0)e^{\frac{ik}{2z}\rho^{2}}dx'}dy'} (2)

Возможно, Вы заметите сходство полученного результата с довольно важным объектом гармонического анализа — преобразованием Фурье, и действительно, этот интеграл может быть довольно легко переписан в таком виде, можете попробовать, это несложно.

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

Моделирование

Выбор технологий

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

Pipeline программы

Заметим, что из теории следует, что каждая следующая перегородка использует только выходной слой предыдущей для того, чтобы получить картинку за собой, таким образом, появляется естественная идея представить программу в виде некоторого конвейера, который на вход получает input-слой и cover-слой, а на выход выдаёт output-слой — то, что получается на экране.

Ускорение с помощью torch

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

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

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

За параллельные вычисления в решении отвечала библиотека Torch, соответственно, все вычисления происходили в матричном виде. Даже распараллелив только преобразования над точками, удалось достичь ускорения в сотню раз, уменьшив время на один слой примерно до 10 минут (при удовлетворительной дискретизации) вместо ~10–15 часов при полностью последовательном вычислении.

Например, используя torch, функция, вычисляющая значение точки входного слоя после преобразований, выглядит так:

def func(self, x, y, z, px, py, c):
  k = 2 * pi / CONFIG["lmbd"]
  kf = exp(complex(0, 1) * k * z) / (complex(0, 1) * CONFIG["lmbd"] * z)
  mult = torch.exp(torch.mul(torch.square(torch_l2(px, py, x, y)), (complex(0, 1) * k / 2 / z)))
  ind = c
  return torch.mul(torch.mul(kf, mult), ind)

Информирование

Для отслеживания прогресса выполнения расчётов использовалась библиотека tqdm.

Примеры использования

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

Везде далее λ = 500 нм.

Круглая щель

На примере этой задачи будет продемонстрирован полный pipeline, далее будут приводиться только результаты на последнем экране.

r_{\text{щели}} = 0.001\text{м}, z = 2\text{м}
pipeline приложения на примере круглой щели
pipeline приложения на примере круглой щели

Двухщелевой эксперимент

d_{между} = 0.005\text{м}, r_{\text{щели}} = 0.001\text{м}, z = 5\text{м}


Определите расстояние до первого пика. (Ответ d = 0.0005 совпадает с моделью)

Изображения на экранах для различного приближения
Изображения на экранах для различного приближения
Графики значений на срезе
Графики значений на срезе

Пятна Пуассона

Пусть

\frac{d^{2}}{l\lambda} \gtrsim 1,

где d — диаметр объекта с круглой тенью, l — расстояние до экрана, λ — длина волны, тогда в центре тени, отбрасываемой объектом, должна наблюдаться светлое пятно — пятно Пуассона (или Араго, или Френеля).

В этой статье мы не будем разбирать причину его появления, скажем только, что его существование следовало из теории Френеля и приводилось Пуассоном в качестве опровержения до момента, пока не было получено экспериментально. Получение пятна Пуассона стало веским доказательством в пользу волновой теории Френеля.

r_{\text{диска}} = 0.001\text{м},0.002\text{м},0.004\text{м}, z = 1\text{м}
Изображения на экранах для различных радиусов
Изображения на экранах для различных радиусов
Графики значений на срезе
Графики значений на срезе

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

Заключение

В результате работы была написана программа, позволяющая моделировать дифракцию через численное нахождение интеграла Френеля методом Монте-Карло. Решение было опробовано на заслонках различных форм, были получены пятна Пуассона.

Ссылки

GitHub проекта
Статья про маятники
Статья про броски объектов под углом к горизонту

Комментарии (7)


  1. mobi
    23.05.2024 05:05
    +1

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

    А если заменить E(x,y,0) на E_{пусто}(x,y,0)-E_{объект}(x,y,0) и посчитать первый интеграл аналитически для E_{пусто}(x,y,0)=A (будет, очевидно, A\exp(i\pi k))?


    1. yawaedafreven Автор
      23.05.2024 05:05

      Если я правильно вас понял, вы считаете, что E(x', y', 0)- это индикатор (поправьте, если это не так), думаю, в такой постановке задачи предложенное решение довольно удачно и действительно решает проблему. Однако относительно общей задачи, где возможно наличие нескольких перегородок, входной слой не всегда представляет собой индикатор, из-за чего аналитическое решение кажется мне затруднительным. Наверное, можно было бы добавить особенное поведение, включающее ваше замечание, для случая, когда pipeline замечает лишь одну перегородку — это помогло бы повысить качество получаемых изображений.


      1. mobi
        23.05.2024 05:05
        +1

        Я просто про то, что в случае отверстия логично брать интеграл по отверстию, а в случае тени от объекта - интеграл по тени и вычитать его из незатененного значения, которое обычно известно.


  1. khe404
    23.05.2024 05:05
    +1

    Интересная задачка

    Вы поставили перед собой цель смоделировать дифракцию через численное нахождение интеграла Френеля методом Монте-Карло.

    И визуально сравнили полученный результат с наблюдаемым в результате эксперимента.

    Полученный результат интегрирования ожидаемо дал значения заложенные в теоретическую модель.

    В статье смущает вот что: длительность вычисления зависит от количества точек для которых проводится расчет. Очевидно -- взяв меньше точек вы сможете оптимизировать трудозатраты при сохранении визуализации. (Как вы выбирали дискретность вашего расчета?)

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

    В отношении пятна Пуассона, вы получаете интенсивность исходного поля, которая похоже задается как 1, просто исходя из того что exp(0) == 1;


    1. yawaedafreven Автор
      23.05.2024 05:05

      Вы правы, я забыл поменять подпись к оси, и да, интенсивность исходного поля равна 1, впрочем, этот параметр настраиваемый.

      Касательно выбора дискретизации, на данный момент она задаётся вручную для каждого слоя. Для представленных результатов она подбиралась из расчёта, что на каждую из интересующих зон Френеля попадает хотя бы 5–7 точек, если смотреть по диаметру.


      1. khe404
        23.05.2024 05:05

        Тогда не очень понятно как ваши вычисления могут занимать минуты, если исходить из того что в слое 200 на 200 точек то он будет обрабатываться за доли секунды.


        1. yawaedafreven Автор
          23.05.2024 05:05

          Есть несколько возможных причин, например, float32, в целом Python или не самое производительное железо (R3600X). Думаю, основная проблема, помимо железа — это, конечно, Python, про что я написал в статье.

          Возможно, вы не в полной мере осознаёте, насколько много вычислений требуется. Возьмём за пример ваши 200 на 200 точек — это один слой, но для каждой точки считается интеграл по предыдущему (в реализации учитывались все точки предыдущего слоя), нетрудными вычислениями получаем 1,6 * 10^9 в однопоток, пусть у меня 12 потоков, это ~ 10^8, и это не учитывая сложности вычислений. Проведя сравнение на i++ 10^8 раз, я получаю различие в ~78 раз (6с против 0.077c) по сравнению с C++. Сравнительные скорости работы именно torch-операций на разных языках я не сравнивал, но предполагаю, что в этом аспекте языки скорее покажут более схожие результаты в силу особенностей Python-библиотек.

          Таким образом, даже в идеальной вселенной, где Python показывает себя в 75 раз хуже, время работы с 10 минут сократится до 8 секунд, но никак не до долей секунды. В реальном же мире, где Torch не наследуют проблем языка по части производительности и в силу того, что фактически сами Python операции занимают колоссально меньшую часть от времени работы программы, по сравнению с остальными вычислениями, сомневаюсь, что прирост будет больше, чем в 10–15 раз (даже такие числа я нахожу крайне преувеличенными), что даёт около минуты в C++ реализации. Это, конечно, куда лучше, но с вашим тезисом про «доли секунды» я не могу согласиться.

          В прочем, я готов его принять, если вы приведёте примеры реализации, где подобная производительность достигается на сравнимом железе.