Объяснение для магических чисел в решетке моделирования Больцмана?

Question

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

Iудалось выяснить, что факторы 1/9, 4/9 и 1/36 связаны с длинами векторов, соединяющих центры клеток вдоль разных направлений решетки, и я даже могу найти ресурсы, которые объясняют, какие разные факторыиспользовать для различных типов решеток (с решеткой D2Q9, используемой в этом коде). Но я не смог найти ничего, что объясняло бы, как вы переходите от общего векторного уравнения, которое определяет шаг столкновения решеточного алгоритма Больцмана к конкретным девяти строкам арифметики реализации, показанным ниже. В частности, откуда взялись эти факторы 3, 1,5 и 4,5?

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

function collide(viscosity) { // kinematic viscosity coefficient in natural units
  var omega = 1 / (3*viscosity + 0.5); // reciprocal of relaxation time
  for (var y=1; y<ydim-1; y++) {
    for (var x=1; x<xdim-1; x++) {
      var i = x + y*xdim; // array index for this lattice site
      var thisrho = n0[i]+nN[i]+nS[i]+nE[i]+nW[i]+nNW[i]+nNE[i]+nSW[i]+nSE[i];

      // macroscopic horizontal velocity
      var thisux = (nE[i]+nNE[i]+nSE[i]-nW[i]-nNW[i]-nSW[i])/thisrho;

      // macroscopic vertical velocity
      var thisuy = (nN[i]+nNE[i]+nNW[i]-nS[i]-nSE[i]-nSW[i])/thisrho;

      var one9thrho = thisrho / 9;
      var one36thrho = thisrho / 36;
      var ux3 = 3 * thisux;
      var uy3 = 3 * thisuy;
      var ux2 = thisux * thisux;
      var uy2 = thisuy * thisuy;
      var uxuy2 = 2 * thisux * thisuy;
      var u2 = ux2 + uy2;
      var u215 = 1.5 * u2;

      n0[i]  += omega * (4/9*thisrho * (1                              - u215) - n0[i]);
      nE[i]  += omega * (  one9thrho * (1 + ux3       + 4.5*ux2        - u215) - nE[i]);
      nW[i]  += omega * (  one9thrho * (1 - ux3       + 4.5*ux2        - u215) - nW[i]);
      nN[i]  += omega * (  one9thrho * (1 + uy3       + 4.5*uy2        - u215) - nN[i]);
      nS[i]  += omega * (  one9thrho * (1 - uy3       + 4.5*uy2        - u215) - nS[i]);
      nNE[i] += omega * ( one36thrho * (1 + ux3 + uy3 + 4.5*(u2+uxuy2) - u215) - nNE[i]);
      nSE[i] += omega * ( one36thrho * (1 + ux3 - uy3 + 4.5*(u2-uxuy2) - u215) - nSE[i]);
      nNW[i] += omega * ( one36thrho * (1 - ux3 + uy3 + 4.5*(u2-uxuy2) - u215) - nNW[i]);
      nSW[i] += omega * ( one36thrho * (1 - ux3 - uy3 + 4.5*(u2+uxuy2) - u215) - nSW[i]);
    }
  }
}

Answer 1

• Код (который, кажется, написан на C #) начинается с вычисления частоты релаксации омега от вязкости , которая устанавливается путем применения определенного числа Рейнольдса Re: = UL / nu. (Как правило, это константа, и поэтому ее можно рассчитать заранее, а не в каждом цикле.)

var omega = 1 / (3*viscosity + 0.5);

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

for (var y=1; y<ydim-1; y++) {
   for (var x=1; x<xdim-1; x++) {`

• Для хранения популяций он полагается на глобальные одномерные массивы для каждого направления решетки \ alpha и, таким образом, превращаются в линейное индексирование , чтобы выбрать правильный узел. Направление в этом случае называется 0 для остального узла и в соответствии с соответствующим направлением (север, восток и т. Д.) Для других групп населения. Обычно коды используют вместо этого соглашение о нумерации.

var i = x + y*xdim;

• Он считывает совокупности с линейным индексированием и суммирует их для получения плотности (импульс нулевого порядка).

var thisrho = n0[i]+nN[i]+nS[i]+nE[i]+nW[i]+nNW[i]+nNE[i]+nSW[i]+nSE[i];

• Затем он определяет скорость , применяя импульс первого порядка \ rho * u_i = \ sum_ \ alpha e_ {i \ alpha} f_ \альфа для текущей решетки. В этом случае x указывает в горизонтальном направлении (восток-запад), а y указывает в вертикальном направлении (север-юг), где восток и север являются положительными направлениями. Это значит \ alpha = {север, северо-восток, северо-запад, юг, юго-восток, юго-запад}.

var thisux = (nE[i]+nNE[i]+nSE[i]-nW[i]-nNW[i]-nSW[i])/thisrho;

var thisuy = (nN[i]+nNE[i]+nNW[i]-nS[i]-nSE[i]-nSW[i])/thisrho;

• Тогда несколько переменных объявлено, что будет повторно использовано в следующем разделе: one9thrho и one36thrho представляют собой комбинацию весов и плотности для горизонтальной / вертикальной и диагональной совокупностей соответственно. Смущает то, что ux3 означает 3 * ux, тогда как ux2 предназначен для представления ux ^ 2.

var one9thrho = thisrho / 9;
var one36thrho = thisrho / 36;
var ux3 = 3 * thisux;
var uy3 = 3 * thisuy;
var ux2 = thisux * thisux;
var uy2 = thisuy * thisuy;
var uxuy2 = 2 * thisux * thisuy;
var u2 = ux2 + uy2;
var u215 = 1.5 * u2;

• Последний шаг обрабатывает столкновение f_ \ alpha ^ t (где t обозначает временный, поскольку за ним следует потоковая передача) = f_ \ alpha + \ omega * (f_ \ alpha ^ {eq} - f_ \ alpha). Каждая отдельная строка отвечает за одно дискретное направление \ alpha, поэтому за одну запись векторного уравнения. Примечание:

  • В этом случае f_ \ alpha ^ t снова сохраняется в f_ \ alpha и, следовательно, достаточно увеличить (+ =) f_ \ alpha (массивы nNWSE) от \ omega * (f_ \ alpha ^ {eq} - f_ \ alpha). Первый член в скобках - это функция равновесия, в то время как последний член соответствует текущему значению f_ \ alpha.

  • Функция несжимаемого равновесия задается как f_ \ alpha ^ {eq} = w_\ alpha \ rho (1 + e_ {i \ alpha} u_i / c_s ^ 2 + (e_i \ alpha u_i) ^ 2 / (2 c_s ^ 4) - u_i ^ 2 / (2 c_s ^ 2)), где мы должны сумма по всем терминам, содержащим индекс i дважды . Скорость решетки звука c_s задается как 1 / \ sqrt (3) * dx / dx и, следовательно, f_ \ alpha ^ {eq} = w_ \ alpha \ rho (1 + 3 e_ {i \ alpha} u_i + 4,5 (e_i\ alpha u_i) ^ 2 - 1,5 u_i ^ 2). Термины one9thr и one36thrho соответствуют термину перед скобками. Сумма ux3 и uy3 соответствует второму члену в скобках 3 e_ {i \ alpha}. Второй последний член 4.5 (e_i \ alpha u_i) ^ 2 соответствует 4.5 u_x ^ 2 или 4.5 u_y ^ 2 для горизонтального или вертикального направления и 4.5 (+ -u_x + - uy) ^ 2 для диагонального направления в качестве обоих компонентовприсутствуют и, следовательно, приводит к 4,5 (u_x ^ 2 + u_y ^ 2 + - u_x u_y). При вычитании u215 = 1,5 * u2 = 1,5 * (ux + uy) соответствует последнему члену. Вам нужно будет принять во внимание, что соответствующая проекция скорости e_ {i \ alpha} скорости решетки \ vec e_ \ alpha в направлении i равна 0 или + -1. Последний отвечает за знаки

n0[i]  += omega * (4/9*thisrho * (1                              - u215) - n0[i]);
nE[i]  += omega * (  one9thrho * (1 + ux3       + 4.5*ux2        - u215) - nE[i]);
nW[i]  += omega * (  one9thrho * (1 - ux3       + 4.5*ux2        - u215) - nW[i]);
nN[i]  += omega * (  one9thrho * (1 + uy3       + 4.5*uy2        - u215) - nN[i]);
nS[i]  += omega * (  one9thrho * (1 - uy3       + 4.5*uy2        - u215) - nS[i]);
nNE[i] += omega * ( one36thrho * (1 + ux3 + uy3 + 4.5*(u2+uxuy2) - u215) - nNE[i]);
nSE[i] += omega * ( one36thrho * (1 + ux3 - uy3 + 4.5*(u2-uxuy2) - u215) - nSE[i]);
nNW[i] += omega * ( one36thrho * (1 - ux3 + uy3 + 4.5*(u2-uxuy2) - u215) - nNW[i]);
nSW[i] += omega * ( one36thrho * (1 - ux3 - uy3 + 4.5*(u2+uxuy2) - u215) - nSW[i]);

Подобный код вряд ли будет очень эффективным. Чтобы ускорить код, вы должны использовать два массива для f_ \ alpha и f_ \ alpha ^ t и выполнять столкновение и потоковую передачу за один шаг вместо двух. Кроме того, функцию равновесия можно переписать, чтобы пересчитать как можно меньше терминов, комбинируя омега и oneXXthrho, и переписать квадратные термины. Обратитесь к коду, который сопровождает «Метод решетчатого Больцмана: принципы и практика» для лучшего написания кода. Это должно уже повысить производительность как минимум в два раза. Чтобы симулировать 3D на вашем компьютере, вам придется применить несколько более сложных оптимизаций . Кроме того, существуют лучшие форумы по этой теме, например, Palabos (Университет Женевы) и OpenLB (Технологический институт Карлсруэ).