Проста реалізація FDTD на Java

FDTD (Finite Difference Time Domain) — метод скінченних різниць у часовій області — самий «чесний» метод рішення завдання електродинаміки від низьких частот до видимого діапазону. Суть — рішення рівнянь Максвелла «в лоб». Тут непогано розписано. Особливо подивіться сітку.

Задача вирішувалася в двовимірному випадку простий явної різницевої схемою. Неявні схеми я не люблю, і вони потребують багато пам'яті. Розрахунок з нормальною точністю вимагає сіток малого кроку, порівняно з більш простими методами потрібно дуже багато часу. Тому максимальний упор зроблений на продуктивність.

Представлена реалізація алгоритму на Java і C++.

image

Передмова
Близько шести років основним моїм мовою для розрахунків і дрібниць був Matlab. Причина — простота написання і візуалізації результату. А так як я перейшов з Borland C++ 3.1, то прогрес можливостей був очевидний. (В Python я ніколи не шарив, а в С++ тоді слабо).

FDTD потрібний був мені для розрахунків як надійний і достовірний метод. Почав вивчати питання в 2010-11 роках. Наявні пакети або було незрозуміло, як використовувати, або не вміли те, що мені потрібно. Вирішив написати свою програму для чіткого контролю над усім. Прочитавши класичну статтю «Чисельного solution of initial boundary value problems involving Maxwell's equations in isotropic media», написав тривимірний випадок, але потім спростив його до двовимірного. Чому 3D — це складно, поясню потім.

Потім максимально оптимізував і спрощував код Matlab. Після всіх поліпшень вийшло, що сітка 2000х2000 проходиться за 107 хвилин. На i5-3.8 ГГц. Тоді такій швидкості вистачало. З корисного було те, що в Матлабе відразу йшов розрахунок комплексних полів, і легко було показувати картинки поширення. Також всі вважалося на даблах, тому що на швидкість Матлаба практично не впливало. Так, стандартний мій розрахунок — один прохід світла по фотонного кристалу.

Проблеми було дві. Знадобилося з високою точністю обчислювати спектри, а для цього використовувати велику ширину сітки. А фізичне збільшення області за обома координатами в 2 рази збільшує час розрахунку в 8 раз (розмір кристала*час проходу).

Matlab я продовжував використовувати, але став Java-програмістом. І, порівнюючи продуктивність різних алгоритмів, почав щось підозрювати. Наприклад, сортування бульбашкою — тільки цикли, масиви і порівняння — в Матлабе працює у 6 разів повільніше, ніж С++ або Java. І цього для нього ще добре. Пятимерный цикл для гіпотези Ейлера в Матлабе в 400 разів довше.

Спочатку почав писати FDTD на C++. Там був вбудований std::complex. Але тоді я закинув цю ідею. (У Матлабе не такі дужки, копіпаст не працював, довелося витратити краплю часу). Зараз перевірив C++ — комплексна математика дає втрати швидкості в 5 разів. Це занадто багато. В результаті написав на Java.

Трохи про питання «Чому на Java?». Деталі продуктивності опишу потім. Якщо коротко — на простому не ООП коді, тільки арифметика і цикли — C++ з оптимізацією О3 або такої ж швидкості, або до +30% швидше. Просто зараз я краще розуміюся на Java, інтерфейсі і роботи з графікою.

FDTD – детальніше
Тепер перейдемо до коду. Спробую показати все, з поясненнями завдання і алгоритму. У двовимірному випадку система рівнянь Максвелла розпадається на дві незалежні підсистеми — ТІ і ТМ хвилі. Код показаний для ТІ. Є 3 компоненти — електричне поле Ez і магнітні Hx, Hy. Для спрощення час має розмірність метрів.

Спочатку я думав, що від float сенсу немає, так як всі обчислення приводяться до double. Тому наведу код для double — в ньому менше зайвого. Всі масиви мають розміри +1, щоб індекси збігалися з матрицями Матлаба (від 1, а не 0).

Де-то в початковому коді:

public static int nx = 4096;// ширина області. Ступінь 2 для Фур'є.
public static int ny = 500;// довжина області

Основний метод:

Ініціалізація початкових змінних
public static double[][] callFDTD(int nx, int ny, String method) {
int i, j;// індекси для циклів
double x; //метрична координата

final double lambd = 1064e-9; // Довжина хвилі в метрах
final double dx = lambd / 15; // розмір кроку сітки. Максимум λ/10. Чим менше, тим вище точність.
final double dy = dx; //реально в коді не використовується. Квадратні клітинки
final double period = 2e-6; // поперечний період решітки фотонного кристала
final double Q = 1.0;// специфічне ставлення
final double n = 1;// мінімальний показник заломлення кристала
final double prodol = 2 * n * period * period / lambd / Q; //поздовжній період

final double omega = 2 * PI / lambd; // циклічна частота світлової хвилі
final double dt = dx / 2; // крок по часу. НЕ більше, ніж по координаті
final double tau = 2e-5 * 999;// час загасання джерела світла. Потрібно для запуску коротких імпульсів. Якщо велике, загасання непомітно.
final double s = dt / dx; //відношення кроку по часу і простору
final double k3 = (1 - s) / (1 + s);// для граничних умов
final double w = 19e-7;// напівширина гаусового пучка
final double alpha = sin(0.0 / 180 * PI);// кут лазерного променя від нормалі. В радіанах.

double[][] Ez = new double[nx + 1][ny + 1]; 
double[][] Hx = new double[nx][ny]; // основні масиви полів
double[][] Hy = new double[nx][ny];


final int begin = 10; // координата, де пустка змінюється кристалом. Без необхідності багато не ставити.
final double mod = 0.008 * 2;//максимальна глибина модуляції діелектричної проникності = 2*Δn;
final double ds = dt * dt / dx / dx;// масив для підгону систем числення

Ініціалізація діелектричної проникності (решітка кристалу). Точніше, зворотна величина з константами. Використовується дискретна функція (0 або 1). Це мені здається ближче до реальних зразків. Звичайно, сюди можна написати що завгодно:

double[][] e = new double[nx + 1][ny + 1];
for (i = 1; i < nx + 1; i++) {
for (j = 1; j < ny + 1; j++) {
e[i][j] = ds / (n + ((j < begin) ? 0 : (mod / 2) * (1 + signum(-0.1 + cos(2 * PI * (i - nx / 2.0 + 0.5) * dx / period) * sin(2 * PI * (j - begin) * dy / prodol)))));
}
}

Решітка кристалу виглядає приблизно так:

image

Масиви, використовувані для граничних умов:

double[][] end = new double[2][nx + 1]; 
double[][] top = new double[2][ny + 1];
double[][] bottom = new double[2][ny + 1];

Межа рахунки за часом. Оскільки у нас крок dt = dx / 2, стандартний коефіцієнт 2. Якщо середовище щільна, чи потрібно йти під кутом більше.

final int tMax = (int) (ny * 2.2);

Починаємо головний цикл:

for (int t = 1; t <= tMax; t++) {
double tt = Math.min(t * s + 10, ny - 1);

Змінна tt тут враховує обмеження швидкості світла, з невеликим запасом. Ми будемо вважати тільки ту сферу, куди світ міг дійти.

Замість комплексних чисел я вважаю окремо дві компоненти — синус і косинус. Я подумав, що для швидкості краще скопіювати шматок, ніж робити вибір всередині циклу. Можливо, заміню на виклик функції або лямбду.

switch (method) {
case "cos":
for (i = 1; i < = nx - 1; i++) {
x = dx * (i - (double) nx / 2 + 0.5); // поперечна координата в метрах
Ez[i][1] = exp(-pow(x, 2) / w / w - (t - 1) * dt / tau) * cos((x * alpha + (t - 1) * dt) * omega);
}
break;

sin
case "sin":
for (i = 1; i < = nx - 1; i++) {
x = dx * (i - (double) nx / 2 + 0.5);
Ez[i][1] = exp(-pow(x, 2) / w / w - (t - 1) * dt / tau) * sin((x * alpha + (t - 1) * dt) * omega);
}
break;
}


Тут у нас на вхід області (в крайню ліву координату) надходить гаусів промінь напівширини w, під кутом alpha, осцилирующий у часі. Саме так і виникає «лазерне» випромінювання потрібної нам частоти/довжини хвилі.

Далі копіюємо тимчасові масиви під поглинаючі граничні умови Мура:

for (i = 1; i < = nx; i++) {
end[0][i] = Ez[i][ny - 1];
end[1][i] = Ez[i][ny];
}
System.arraycopy(Ez[1], 0, top[0], 0, ny + 1);
System.arraycopy(Ez[2], 0, top[1], 0, ny + 1);
System.arraycopy(Ez[nx - 1], 0, bottom[0], 0, ny + 1);
System.arraycopy(Ez[nx], 0, bottom[1], 0, ny + 1);

Тепер переходимо до головного обчисленню — наступного кроку по полю. Особливість рівнянь Максвелла — зміна в часі магнітного поля залежить тільки від електричного, і навпаки. Це дозволяє написати просту різницеву схему. Вихідні формули були такі:

image
Всі зайві константи я перерахував заздалегідь, замінив розмірність Н, вніс і врахував діелектричну проникність. Оскільки оригінальні формули для сітки, зрушеної на 0,5, потрібно не помилитися з індексами масивів Е і Н. вони різної довжини — Е на 1 більше.

Цикл по області для Е:

for (i = 2; i <= nx - 1; i++) {
for (j = 2; j <= tt; j++) {//можна тут оптимізувати порядок доступу до комірок?
Ez[i][j] += e[i][j] * ((Hx[i][j - 1] - Hx[i][j] + Hy[i][j] - Hy[i - 1][j]));
}
}

А тепер нарешті-то застосовуємо граничні умови. Вони потрібні тому, що різницева схема не вважає крайні клітинки — для них немає формул. Якщо нічого не робити, світло буде відбиватися від стінок. Тому застосовуємо метод, який мінімізує відображення при нормальному падінні. Обробляємо три сторони — верх, низ і право. Втрати продуктивності на граничні умови близько 1% (тим менше, чим більше завдання).

for (i = 1; i < = nx; i++) { 
Ez[i][ny] = end[0][i] + k3 * (end[1][i] - Ez[i][ny - 1]);//end
}
for (i = 1; i < = ny; i++) {
Ez[1][i] = top[1][i] + k3 * (top[0][i] - Ez[2][i]);//kray verh
Ez[nx][i] = bottom[0][i] + k3 * (bottom[1][i] - Ez[nx - 1][i]);
}

Зліва кордон особлива — в ній генерується промінь. Також, як в минулий раз. Просто на 1 крок далі по часу.

Лазерний
switch (method) {
case "cos":
for (i = 1; i < = nx - 1; i++) {
x = dx * (i - (double) nx / 2 + 0.5);
Ez[i][1] = exp(-pow(x, 2) / w / w - (t - 1) * dt / tau) * cos((x * alpha + t * dt) * omega);
}
break;
case "sin":
for (i = 1; i < = nx - 1; i++) {
x = dx * (i - (double) nx / 2 + 0.5);
Ez[i][1] = exp(-pow(x, 2) / w / w - (t - 1) * dt / tau) * sin((x * alpha + t * dt) * omega);
}
break;
}


Залишилося тільки порахувати магнітне поле:

for (i = 1; i < = nx - 1; i++) { // main Hx Hy
for (j = 1; j < = tt; j++) {
Hx[i][j] += Ez[i][j] - Ez[i][j + 1];
Hy[i][j] += Ez[i + 1][j] - Ez[i][j];
}
}
} // кінець 3-мірного циклу.

І ще дрібниці: передати кінцевий відрізок для обчислення перетворення Фур'є — знаходження картини в дальній зоні (простір напрямків):

int pos = method.equals("cos") ? 0 : 1; //який компонент записувати
BasicEx.forFurier[pos] = new double[nx]; //там в коді масив для Фур'є
int endF = (int) (ny * 0.95);// довільна координата в кінці
for (i = 1; i < = nx; i++) {
BasicEx.forFurier[pos][i - 1] = Ez[i][endF];
for (j = 1; j < = ny; j++) {
Ez[i][j] = abs(Ez[i][j]);// ABS
} //я ось зараз зрозумів, що брати модуль не треба - я беру квадарат
}
Hx = null; //я вірю, що це очистить пам'ять
Hy = null;
e = null;
return Ez;
}

Потім складаємо квадрати двох компонент, і виводимо зображення інтенсивності.
for (int i = 0; i < nx + 1; i++) {
for (int j = 0; j < ny + 1; j++) {
co.E[i][j] = co.E[i][j] * co.E[i][j] + si.E[i][j] * si.E[i][j];
}
}


Окремо беремо перетворення Фур'є:

fft.fft(forFurier[0], forFurier[1]);

Оскільки в швидкому Фур'є не розбираюся, взяв перший-ліпший. Мінус — ширина тільки ступінь двійки.

Про продуктивності
Переходимо до найцікавішого для мене — що доброго ми отримали, перейшовши з Матлаба на Java. У Матлабе я в свій час оптимізував все, що зміг. У Java — в основному внутрішній цикл (складність n^3). У Матлаба вже враховано те, що він вважає відразу дві компоненти. Швидкість на першому етапі (більше — краще):
Matlab 1
Java double 50
C++ gcc double 48
C++ MSVS double 55
C++ gcc float 73
C++ MSVS float 79
Опишемо основних учасників змагання.

  • Matlab — 2011b і 2014. Перехід на 32-бітні числа дає дуже малий приріст швидкості.
  • Java — спочатку 7u79, але в основному 8u102. Мені здалося, що 8 трохи краще, але детально не порівняв.
  • Microsoft VisualStudio 2015 конфігурація release
  • MinGW gcc — фіг знає який, грудень 2014, завжди оптимізація O3. Оптимізація під Core i, без AVX. AVX у мене на ноут немає, там де є приросту не давав.
Тестові машини:
  • Pentium 2020m 2.4 GHz, ddr3-1600, 1 канал
  • Core i5-4670 3.6–3.8 GHz, ddr3-1600 2 каналу

  • Core i7-4771 3.7–3.9 GHz, ddr3-1333 2 каналу
  • Athlon x3 3.1 GHz, ddr3-1333, дуже повільний контролер пам'яті.

2-ядерний етап
Спочатку я вважав ТІ і ТМ компоненти послідовно. До речі, це єдиний варіант, при нестачі пам'яті. Потім написав два потоки — прості Runnable. Ось тільки прогресу була мало. Всього на 20-22% швидше, ніж 1-потокове. Почав шукати причини. Потоки працювали нормально. 2 ядра стабільно вантажилися на 100%, не даючи жити ноута.

Потім я порахував продуктивність. Вийшло, що вона впирається в швидкість оперативної пам'яті. Код вже працював на межі читання. «Довелося переходити на float. Перевірка точності показала, що ніяких смертельних похибок немає. Візуальної різниці однозначно ні. Сумарна енергія відрізнялася 8 знаку. Після інтегрального перетворення максимум спектра — на 0,7 e-6.

Але головне — отримав ривок продуктивності. Загальний ефект від 2 ядер і переходу на float +87-102%. (Чим швидше пам'ять і більше ядер, тим краще приріст). Athlon x3 дав мало приросту.
Реалізація на С++ аналогічна — через std::thread (дивись в кінці).

Поточні швидкості (2 threads):
Java double 61
Java float 93-101
C++ gcc float 99-110
Всі оцінки проводилися для одноразового запуску розрахунку. Тому що, якщо повторно запускати з GUI, не закриваючи програму, далі швидкість зростає. Імхо, jit-оптимізація рулить.

4-ядерний етап
Мені здавалося, що ще є куди прискорюватися. Я взявся за 4-поточний варіант. Суть — область ділиться навпіл і вважається в 2 потоках, спочатку по Е, потім по Н.

Спочатку написав його на Runnable. Вийшло жахливо — прискорення було тільки для дуже великої ширини області. Занадто великі витрати на породження нових потоків. Потім освоїв java.util.concurrent. Тепер у мене був фіксований пул потоків, яким давалися завдання:

public static ExecutorService service = Executors.newFixedThreadPool(4);
//......
cdl = new CountDownLatch(2);
NewThreadE first = new NewThreadE(Ez, Hx, Hy, e, 2, nx / 2, tt, cdl);
NewThreadE second = new NewThreadE(Ez, Hx, Hy, e, nx / 2 + 1, nx - 1, tt, cdl);
service.execute(first);
service.execute(second);
try {
cdl.await();
} catch (InterruptedException ex) {
}

І для Н:
cdl = new CountDownLatch(2);
NewThreadH firstH = new NewThreadH(Ez, Hx, Hy, 1, nx / 2, tt, cdl);
NewThreadH secondH = new NewThreadH(Ez, Hx, Hy, nx/2+1, nx-1, tt, cdl);
service.execute(firstH);
service.execute(secondH);
try {
cdl.await();
} catch (InterruptedException ex) {
}


Всередині потоку виконуються підлогу циклу.

class NewThreadE
public class NewThreadE implements Runnable {

float[][] Ez;
float[][] Hx;
float[][] Hy;
float[][] e;
int iBegin;
int iEnd;
float tt;

CountDownLatch cdl;

public NewThreadE(float[][] E, float[][] H, float[][] H2, float[][] eps,
int iBegi, int iEn, float ttt, CountDownLatch cdl) {
this.cdl = cdl;

Ez = E;
Hx = H;
Hy = H2;
e = eps;
iBegin = iBegi;
iEnd = iEn;
tt = ttt;
}

@Override
public void run() {
for (int i = iBegin; i <= iEnd; i++) { // main Ez 
for (int j = 2; j <= tt; j++) {
Ez[i][j] += e[i][j] * ((Hx[i][j - 1] - Hx[i][j] + Hy[i][j] - Hy[i - 1][j]));
}
}
cdl.countDown();
}
}


Аналогічний клас для Н — просто інші кордони і свій цикл.

Тепер приріст від 4 потоків є завжди — від 23 до 49% (.jar). Чим менше ширина, тим краще — судячи по швидкості, ми залазимо в кеш-пам'ять. Найбільша користь буде, якщо вважати довгі вузькі завдання.
Java float 4потока 124-151
C++ gcc float 4потока 115-162
Реалізація на С++ поки що містить тільки прості std::thread. Тому для неї чим ширше, тим краще. Прискорення С++ від 5% при ширині 1024 до 47% на 16384.

Як ми бачимо, в найкращому випадку приріст 49% — майже 3 ядра. Тому маємо кумедний факт — при малих завдання краще всього ставити newFixedThreadPool(3);

Для великих — за кількістю потоків в процесорі — 4 або 6-8. Вкажу ще цікавий факт. На Атлоні х3 був дуже слабкий прогрес від переходу на float і потоків 2 — 32% від обох оптимізацій. Приріст від «4»-ядерного коду на С++ також невеликий — 67% між 1 і 4 ядрами (обидва float). Списувати можна на повільний контролер пам'яті і 32-бітну Вінду.

Але 4-ядерний Java-код відпрацював на відмінно. При 3 потоках Екзекутора +50,2% від 2-ядерної версії для великих завдань. Чомусь гірша 2-ядерна реалізація посилилася максимально можливої багатоядерної.

Останнє зауваження по 4-ядерному кодом Java. Поточні витрати часу:
2 Основні-мірні цикли Е і Н 83%
Інше, вважаючи початкову ініціалізацію 16%
Породження (4) потоків Близько 1% (0,86 по профилировщику)
Я постарався максимально оптимізувати основні цикли, думаючи що інше витрачає мало часу.

Викладаю також повний код 4-ядерного випадку на С++:

.cpp
#include < iostream>
#include <complex>
#include < stdio.h>
#include < sys/time.h>
#include <thread>
using namespace std;

thE void (float** &Ez, float** &Hx, float** &Hy, float** &e,
int iBegin, int iEnd, float tt)
{
for (int i = iBegin; i <= iEnd; i++) // main Ez
{
for (int j = 2; j <= tt; j++)
{
Ez[i][j] += e[i][j] * ((Hx[i][j - 1] - Hx[i][j] + Hy[i][j] - Hy[i - 1][j]));
}
}
}
void thH (float** &Ez, float** &Hx, float** &Hy,
int iBegin, int iEnd, float tt)
{
for (int i = iBegin; i <= iEnd; i++)
{
for (int j = 1; j < = tt; j++)
{
Hx[i][j] += Ez[i][j] - Ez[i][j + 1];
Hy[i][j] += Ez[i + 1][j] - Ez[i][j];
}
}
}

void FDTDmainFunction (char m)
{
//m=c: cos, sin else
int i,j;
float x;
const float dx=0.5*1e-7/1.4;
const float dy=dx;
const float period=1.2 e-6;
const float Q=1.5;
const float n=1;//ne використовується в Епсилон матриці поки =1
const float lambd=1064e-9;
const float prodol=2*n*period*period/lambd/Q;
const int nx=1024;
const int ny=700;

float **Ez = new float *[nx+1];
for (i = 0; i < nx+1; i++)
{
Ez[i] = new float [ny+1];
for (j=0; j<ny+1; j++)
{
Ez[i][j]=0;
}
}
float **Hx = new float *[nx];
for (i = 0; i < nx; i++)
{
Hx[i] = new float [ny];
for (j=0; j<ny; j++)
{
Hx[i][j]=0;
}
}
float **Hy = new float *[nx];
for (i = 0; i < nx; i++)
{
Hy[i] = new float [ny];
for (j=0; j<ny; j++)
{
Hy[i][j]=0;
}
}

const float omega=2*3.14159265359/lambd;
const float dt=dx/2;
const float s=dt/dx;//for MUR

const float w=40e-7;
const float alpha =tan(15.0/180*3.1416);

float** e = new float *[nx+1];
for (i = 0; i < nx+1; i++)
{
e[i] = new float [ny+1];
for (j=0; j<ny+1; j++)
{
e[i][j]=dt*dt / dx/dx/1;
}
}
const int tmax= (int) ny*1.9;
for (int t=1; t<=tmax; t++)
{
float tt=min(t*s+10,(float)(ny-1));
if (m == 'c')
{
for (i=1; i < =nx-1; i++)
{
x = dx*(i-(float)nx/2+0.5);
Ez[i][1]=exp(-pow(x,2)/w/w)*cos((x*alpha+(t-1)*dt)*omega);
}
}
else
{
for (i=1; i < =nx-1; i++)
{
x = dx*(i-(float)nx/2+0.5);
Ez[i][1]=exp(-pow(x,2)/w/w)*sin((x*alpha+(t-1)*dt)*omega);
}
}
std::thread thr01(thE, std::ref(Ez), std::ref(Hx), std::ref(Hy), std::ref(e), 2, nx / 2, tt);
std::thread thr02(thE, std::ref(Ez), std::ref(Hx), std::ref(Hy), std::ref(e), nx / 2 + 1, nx - 1, tt);
thr01.join();
thr02.join();
// H
for (i=1; i < =nx-1; i++)
{
x = dx*(i-(float)nx/2+0.5);
Ez[i][1]=exp(-pow(x,2)/w/w)*cos((x*alpha+t*dt)*omega);
}
std::thread thr03(thH, std::ref(Ez), std::ref(Hx), std::ref(Hy), 1, nx / 2, tt);
std::thread thr04(thH, std::ref(Ez), std::ref(Hx), std::ref(Hy), nx / 2 + 1, nx - 1, tt);
thr03.join();
thr04.join();
}
}
int main()
{
struct timeval tp;
gettimeofday(&tp, NULL);
double ms = tp.tv_sec * 1000 + (double)tp.tv_usec / 1000;

std::thread thr1(FDTDmainFunction, 'c');
std::thread thr2(FDTDmainFunction, 's');
thr1.join();
thr2.join();

gettimeofday(&tp, NULL);
double ms2 = tp.tv_sec * 1000 + (double)tp.tv_usec / 1000;
cout << ms2-ms << endl;
return 0;
}

Це найпростіший код, без нормальної решітки і граничних умов. Цікаво, чи можна в З++ ефективніше оголошувати 2-мірні масиви і викликати потоки?

8 ядер
А більше ядер у мене немає. Якщо б мої завдання не упиралися в пам'ять, ділив би далі. Начебто, ThreadPool дає малі витрати. Або переходити на Fork-Join Framework.

Я міг би перевірити — відмовитися від перших двох потоків, і ділити одну задачу на 4-8 шматків на i7. Але сенс буде, тільки якщо хтось протестує на машині з швидкою пам'яттю — DDR-4 або 4 канали.

Кращий спосіб позбутися від браку швидкості пам'яті відеокарти. Переходити на Cuda мені заважає брат, який забороняє оновлювати відеодрайвер (З незнання і Cuda).

Підсумок і післямови
Я можу швидко вирішити будь-яку потрібну мені 2-вимірну задачу з хорошою точністю. Сітка 4096х2000 проходиться на 4-ядернике за 106 секунд. Це буде 300 мікрон х 40 шарів — максимальні зразки у нас.

У 2D при 32-бітної точності потрібно трохи пам'яті — 4байта*4массива*2 комплексних компоненти = 32 байта / піксель в гіршому випадку.

У 3D все набагато гірше. Компонент вже шість. Можна відмовитися від двох потоків — вважати компоненти послідовно і писати на гвинт. Можна не зберігати масив діелектричної проникності, а вважати в циклі або обійтися дуже маленьким періодичним ділянкою. Тоді в 16 ГБ оперативы (максимум у мене на роботі) влізе область 895х895х895. Це буде нормально «щоб подивитися». Зате буде вважатися «всього» 6-7 годин один прохід. Якщо завдання добре поділиться на 4 паралельних потоку. І якщо знехтувати обчисленням ε

Тільки для спектру мені не вистачить. При ширині 1024 я не бачу необхідні деталі. Потрібно 2048. А це 200+ ГБ пам'яті. Тому тривимірний випадок — це складно. Якщо не розробляти код з кеширующими ССД.

P.S. Оцінки швидкості були досить приблизні. Я перевіряв Матлаб тільки на малих завданнях. Зараз перевірив Java — завдання 2048*1976 (аналог 2000*2000) на 4-ядернике. Час розрахунку 45,5 секунд. Прискорення від Матлаба 141 разів (точно).

Можливі плани на майбутнє:

*) Перевірити швидкість З чистого (не ++). За benchmarksgame він завжди швидше.
1) Перевірити комплексні класи в С і Java. Може, вони реалізовані досить швидко. Правда, боюся, всі вони будуть більше 8 байт.
2) Закинути всі 2 — і 4-ядерні версії в MSVS, знайти параметри оптимізації.
3) Перевірити, чи можуть лябмды/стріми прискорити основний цикл або додаткові.
4) Зробити нормальний GUI для вибору всього і візуалізації результатів.
99) Написати Cuda-версію.
Якщо комусь цікаво, опишу FDTD та інші методи розрахунків, фотонні кристали.

На Github виклав 2 версії:

1) 2-потокова з зачатками інтерфейсу вибору параметрів
2) 4-потокова

Обидві малюють картинку і спектр. Просто поки піксель в піксель — не беріть ширину вище 2048. Ще вміють приймати розміри області з консолі.
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.