[Из песочницы] Простая реализация FDTD на Java10.01.2017 16:48

Задача решалась в двумерном случае простой явной разностной схемой. Неявные схемы я не люблю, и они требуют много памяти. Расчет с нормальной точностью требует сеток малого шага, по сравнению с более простыми методами требуется очень много времени. Поэтому максимальный упор был сделан на производительность.

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

Предисловие

FDTD нужен был мне для расчетов как надежный и достоверный метод. Начал изучать вопрос в 2010–11 годах. Имеющиеся пакеты или было непонятно, как использовать, или не умели то, что мне нужно. Решил написать свою программу для четкого контроля над всем. Прочитав классическую статью »Numerical 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 — детальнее

Изначально я думал, что от 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;// константа для подгона систем исчисления

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)))));
    }
}

Массивы, используемые для граничных условий:

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

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

for (int t = 1; t <= tMax; t++) {
    double tt = Math.min(t * s + 10, ny - 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 - 1) * 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 - 1) * dt) * omega);
            }
            break;
    }

Дальше копируем временные массивы под поглощающие граничные условия Мура:

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);

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

    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]));
        }
    }

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]);//verh kray
        Ez[nx][i] = bottom[0][i] + k3 * (bottom[1][i] - Ez[nx - 1][i]);
    }

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]);

О производительности

• Matlab — 2011b и 2014. Переход на 32-битные числа дает очень малый прирост скорости.
• Java — вначале 7u79, но в основном 8u102. Мне показалось, что 8 немного лучше, но детально не сравнил.
• Microsoft VisualStudio 2015 конфигурация release
• MinGW gcc — фиг знает какой, декабрь 2014, всегда оптимизация O3. Оптимизация под Core i, без AVX. AVX у меня на ноуте нет, там где есть прироста не давал.

2-ядерный этап

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

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

Текущие скорости (2 threads):

4-ядерный этап

Вначале написал его на 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) {
            }

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). Чем меньше ширина, тем лучше — судя по скорости, мы влезаем в кеш-память. Наибольшая польза будет, если считать длинные узкие задачи.

Как мы видим, в самом лучшем случае прирост 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. Текущие затраты времени:

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

#include 
#include 
#include 
#include 
#include 
using namespace std;

void thE (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, else sin
    int i,j;
    float x;
    const float dx=0.5*1e-7/1.4;
    const float dy=dx;
    const float period=1.2e-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