Сайт Александра Зеленина

Александр Зеленин.

Базовый алгоритм FDTD

Во введении описан (кратко) вывод конечно-разностных уравнений в декартовых координатах для двух компонент сетки Yee из шести.

На рис. 2. приводится полная система для всех векторов сетки Yee.

Рис. 2. Полная система уравнений.

Примечания.1. Полуцелые индексы, которые применены в нотации Yee, часто заменяются на целые путем уменьшения на пол-индекса, например, (j+1/2) заменено на (j), а (j-1/2) заменено на (j-1). Это сделано для удобства программирования, т.к. не бывает полуцелых индексов массивов.

2. Переменные ε, σ и μ повсюду должны быть представлены как ε(i,j,k) σ(i,j,k) μ(i,j,k), но на рисунке индексы (i,j,k) опущены для экономии места.

3. ε и μ здесь - это абсолютные проницаемости, т.е. сразу готовое произведение εε_o и μμ_o, как это представлено ранее (см. введение). Путаница в этих обозначениях - широко распространенное явление в литературе, поэтому надо держать ухо востро.

4. Коэффициенты С1 и С2 из (8) (см. введение), на первый взгляд, не такие, как на этом рисунке. На самом деле это то же самое. Попробуйте разделить числитель и знаменатель коэффициентов С1 и С2 (8) на (ε_(i+1/2,j,k)ε_o).

Шесть уравнений базового алгоритма очень просты в реализации. Для этого необходимо создать шесть массивов Ex, Ey, Ez, Hx, Hy и Hz. А также массивы электрофизических характеристик ε σ μ, которые назовем Eps, Sig и Mu. Этих массивов три штуки, если электрофизические характеристики задаются для ячейки Yee в целом. Здесь возможны разные варианты. Если ε σ и μ планируется задавать для каждого вектора свои, то потребуются массивы Epsx, Epsy, Epsz, Sigx и т.д. Если не планируется применение магнитных материалов, то массив Mu можно вообще не создавать. Если планируется проводить расчеты только для идеальных проводников (PEC) в свободном пространстве, то массивы электрофизических констант не нужны вовсе. В этом случае векторы Е, которые принадлежат PEC-объекту, попросту обнуляются, а все остальные считаются принадлежащими свободному пространству.

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

Раньше значительное ускорение вычислений можно было получить создав еще ряд массивов и поместив в них предварительно вычисленные коэффициенты С, С1 и С2 из уравнений (7) и (8) (см. введение). Но сейчас в гигагерцовых процессорах скорость вычислений арифметики с плавающей точкой выросла больше, чем скорость выборки дополнительных элементов массива из памяти, поэтому смысл в предварительном вычислении этих коэффициентов снизился. В программе ZFDTD я пожертвовал скоростью во имя экономии памяти, хотя в ранних версиях программы было наоборот.

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

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

Ход вычислений рассмотрим на примере, написанном на языке Pascal.

Две процедуры вычисляют поля E и H:

//                     ОПРЕДЕЛЕНИЕ МАГНИТНОГО ПОЛЯ H

Procedure TimeStepForH;

  var I,J,K: Integer;

begin

      For I:=1 to XSize do  for J:=1 to YSize-1 do   for K:=1 to ZSize-1 do

         HX[I,J,K]:=HX[I,J,K]-((EZ[I,J+1,K]-EZ[I,J,K])*TY-(EY[I,J,K+1]-EY[I,J,K])*TZ)*AH[I,J,K];

      For I:=1 to XSize-1 do   for J:=1 to YSize do       for K:=1 to ZSize-1 do

          HY[I,J,K]:=HY[I,J,K]-((EX[I,J,K+1]-EX[I,J,K])*TZ-(EZ[I+1,J,K]-EZ[I,J,K])*TX)*AH[I,J,K];

      For I:=1 to XSize-1 do  for J:=1 to YSize-1 do  for K:=1 to ZSize do

          HZ[I,J,K]:=HZ[I,J,K]-((EY[I+1,J,K]-EY[I,J,K])*TX-(EX[I,J+1,K]-EX[I,J,K])*TY)*AH[I,J,K];

end;// КОНЕЦ ОПРЕДЕЛЕНИЯ МАГНИТНОГО ПОЛЯ H

//  ОПРЕДЕЛЕНИЕ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E

Procedure TimeStepForE;

  var I,J,K: Integer;

begin

   for I:=1 to XSize-1 do  for J:=2 to YSize-1 do    for K:=2 to ZSize-1 do

          EX[I,J,K]:=EX[I,J,K]*AE[I,J,K]+((HZ[I,J,K]-HZ[I,J-1,K])*TY-

                                          (HY[I,J,K]-HY[I,J,K-1])*TZ)*BE[I,J,K];

  for I:=2 to XSize-1 do  for J:=1 to YSize-1 do for K:=2 to ZSize-1 do

          EY[I,J,K]:=EY[I,J,K]*AE[I,J,K]+((HX[I,J,K]-HX[I,J,K-1])*TZ-

                                        (HZ[I,J,K]-HZ[I-1,J,K])*TX)*BE[I,J,K]

   for I:=2 to XSize-1 do  for J:=2 to YSize-1 do for K:=1 to ZSize-1 do

        EZ[I,J,K]:=EZ[I,J,K]*AE[I,J,K]+((HY[I,J,K]-HY[I-1,J,K])*TX-

                                        (HX[I,J,K]-HX[I,J-1,K])*TY)*BE[I,J,K];

end;// КОНЕЦ ОПРЕДЕЛЕНИЯ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E

Массивы AE[I,J,K], BE[I,J,K] и AН[I,J,K] - это предварительно вычисленные коэффициенты С1, С2 и С соответственно из (7) и (8). XSize, YSize и ZSize - размеры счетного пространства. А коэффициенты TX, TY и TZ равны 1/∆X, 1/∆Y и 1/∆Z соответственно. Как видно, здесь применены предварительно вычисленные массивы коэффициентов С1, С2 и С.

Приведенный код вычисления полей рассчитан на такое представление объектов, когда считается, что электрофизические характеристики пространства в ячейке (i,j,k) задаются одинаковыми для всех трех Нx,y,z(i,j,k) и трех Ex,y,z(i,j,k) векторов. При таком представлении точность вычислений на границах раздела сред снижается до первого порядка. Обычно это проявляется в снижении на 25..40% вычисленной резонансной частоты, такого же увеличения вычисленных токов в длинных структурах, снижении в разы скорости затухания свободных колебаний и в других ошибках.

Для любителей Фортрана. Важно помнить, что в фортране порядок следования переменных в массивах отличается от нормальных языков (С++, Delphi). Поэтому тройной цикл Паскаля:

For I:=1 to XSize do

      for J:=1 to YSize-1 do

           for K:=1 to ZSize-1 do begin

<операторы цикла>

end;

на Фортране должен выглядеть так:

   do K=1, ZSize-1

   do J=1, YSize-1

   do I=1, XSize

<операторы цикла>

      end do

      end do

    end do

т.е. индексами задом наперед. Иначе скорость счета будет совсем низкой – в 3-4 раза ниже.

Приведенные процедуры вызываются в цикле по времени:

// ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ

T:=0;  //Исходное значение времени

for I:=0 to <заданное количество шагов по времени> do

begin

      TimeStepForH; //Вычисление H по всему объему

      BordersSymmetryH; //Выполнение граничных условий для симметрии по Н

      TimeStepForE; //Вычисление E по всему объему

      BordersABC; // Граничные УСЛОВИЯ ПОГЛОЩЕНИЯ

      BordersSymmetryE;//Выполнение граничных условий для симметрии по Е

      BordersРЕС;//Выполнение граничных условий РЕС

      T:=T+dT; // инкремент времени на один шаг

end; //ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ

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

Граничные условия рассмотрим в следующей главе.