Моделирование электрического поля средствами CUDA

Данная статья написана с целью продемонстрировать как с помощью технологии CUDA можно смоделировать простое взаимодействие заряженых частиц (см. Закон Кулона). Для вывода статической картинки я использовал библиотеку freeglut.
Как пишут частенько на Хабре:

Ещё во времена Turbo Pascal и Borland C 5.02 меня интересовал процесс симуляции простых физических явлений природы — будь то силы гравитационного притяжения или кулоновского взаимодействия. К сожалению, мощности одного процессора с несколькими сотнями МГц не хватало для простых математических операций, коих множество было велико.

Хоть я и не погроммист, мне захотелось однажды вспомнить далёкое детство, да и подучить основы JavaScript, и я как-то смоделировал гравитационное взаимодействие любого разумного числа материальных точек. Так, например, частичка ниже крутится вокруг «солнца» пользуясь Ньютоновским законом:

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

Итак, хватит ностальгировать, давайте погромпрограммировать. Моя платформа — это ноутбук с Ubuntu и видеокартой GeForce GT 540m и установленным CUDA Toolkit 5.0. Если вы начинаете знакомство с CUDA, то могу вам посоветовать интересную книгу CUDA by Example.

Для начала создайте проект из шаблона Nsight и подключите к нему freeglut используя библиотеки GLU, GL, glut. Это позволит нам отобразить полученное векторное поле, которое мы рассчитаем. Приблизительная структура нашего файла с функцией main будет следующая:
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char** argv) {
	// Initialize freeglut
	// ...
	
	// Initialize Host and Device variables
	// ...
	
	// Launch Kernel to render the image
	
	// Display Image
	
	// Free resources
	
	return 0;
}

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

Инициализация окна


Так как я, когда задумывал приложение, не планировал динамическое изменение размера окна, то я решил указать его декларативно, используя глобальные константы:
const int width = 1280;
const int height = 720;

Ах да! Не судите строго за смешивание стилей C++ и C99. Некоторые вещи, типа перегрузки, я порой нахожу полезным, но и от malloc я не отказался, чтобы хоть как-нибудь следовать стилю приложениям с CUDA.

Само отображение окна, обработка клавиши Escape и простой рендер можно осуществить следующим кодом:
#include <GL/freeglut.h>
#include <GL/gl.h>
#include <GL/glext.h>

// ...

void key(unsigned char key, int x, int y) {
	switch (key) {
	case 27:
		printf("Exit application\n");
		glutLeaveMainLoop();
		break;
	}
}

void draw(void) {
	glClearColor(0.0, 0.0, 0.0, 1.0);
	glClear(GL_COLOR_BUFFER_BIT);
	// glDrawPixels(width, height, GL_RGBA, GL_UNSIGNED_BYTE, screen);
	glFlush();
}

// ...
	// Initialize freeglut
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_SINGLE | GLUT_RGBA);
	glutInitWindowSize(width, height);
	glutCreateWindow("Electric field");
	glutDisplayFunc(draw);
	glutKeyboardFunc(key);
	glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE, GLUT_ACTION_CONTINUE_EXECUTION);

	// ...

	// Display Image
	glutMainLoop();

При запуске этого кода, вы ничего кроме чёрного экрана не увидите. Так что я опущу здесь скриншот. Обратите внимание на следующие моменты:
  • glutLeaveMainLoop — когда нажат Escape, то freeglut позволяет прервать главный цикл, чтобы мы смогли освободить ресурсы приложения. Если хочется, то можно использовать и atexit;
  • glDrawPixels — позже я создам буфер из unsigned byte, где я буду хранить все пиксели (их цвета) и потом отображать в окне;
  • glutSetOption — это мне позволить продолжить выполнение кода в main, чтобы очистить ресурсы в виде динамической памяти.

Разбиение вычисление на блоки и потоки


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

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

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

Чтобы узнать количество потоков внутри блока, воспользуемся функцией cudaGetDeviceProperties:
int threadsCount(void) {
	int deviceCount = 0;
	cudaGetDeviceCount(&deviceCount);

	if (deviceCount <= 0) {
		printf("No CUDA devices\n");
		exit(-1);
	}

	cudaDeviceProp properties;
	cudaGetDeviceProperties(&properties, 0);

	return properties.maxThreadsPerBlock;
}

Моя карта возвращает число 1024. К сожалению, я не могу использовать все потоки, потому что в режиме отладки, при делении float чисел, я получаю ошибку cudaErrorLaunchOutOfResources. Я так и не смог отыскать, почему мне не хватило ресурсов (или регистров, как пишет интернет), поэтому единственные три решения, что я смог найти — это снижение числа потоков на блок, использование release режима при компиляции, использование опции maxrregcount с любым числом больше нуля.
Если вы дочитали мою статью до этого момента, то я рад, что не зря сижу ночью перед компом. А если вы знаете ответ, как узнать в чем именно была моя ошибка, и сообщите мне его, то я буду рад вдвойне.

Теперь разобьём наше поле на блоки и потоки:
void setupGrid(dim3* blocks, dim3* threads, int maxThreads) {
	threads->x = 32;
	threads->y = maxThreads / threads->x - 2; // to avoid cudaErrorLaunchOutOfResources error

	blocks->x = (width + threads->x - 1) / threads->x;
	blocks->y = (height + threads->y - 1) / threads->y;
}

// ...
	// Initialize Host
	dim3 blocks, threads;
	setupGrid(&blocks, &threads, threadsCount());
	printf("Debug: blocks(%d, %d), threads(%d, %d)\nCalculated Resolution: %d x %d\n",
		blocks.x, blocks.y, threads.x, threads.y, blocks.x * threads.x,
		blocks.y * threads.y);

Результат будет следующим (я заведомо увеличил количество блоков на единицу в тех случаях, когда размер экрана не кратен числу потоков):
Debug: blocks(40, 24), threads(32, 30)
Calculated Resolution: 1280 x 720

Теория и практика


Заряды, расстояния, кулоны и коэффициенты — это звучит хорошо. Но когда дело доходит до моделирования, то тут мы не можем уже пользоваться супер-маленькими значениями, потому что они могут пропасть в процессе вычисления. Поэтому мы упростим задачу:

Заряд — это структура с координатами (как на экране).
struct charge {
	int x, y, q;
};

Напряжённость электрического поля в точке — это вектор, который мы выразим проекциями и функции, вычисляющей их (вот именно здесь я и получил свою ошибку):
const float k = 50.0f;
const float minDistance = 0.9f; // not to divide by zero
const float maxForce = 1e7f;

struct force {
	float fx, fy;

	__device__ force() :
			fx(0.0f), fy(0.0f) {
	}

	__device__ force(int fx, int fy) :
			fx(fx), fy(fy) {
	}

	__device__ float length2() const {
		return (fx * fx + fy * fy);
	}

	__device__ float length() const {
		return sqrtf(length2());
	}

	__device__ void calculate(const charge& q, int probe_x, int probe_y) {
		// F(1->2) = k * q1 * q2 / r(1->2)^2 * vec_r(1->2) / abs(vec_r(1->2))
		// e = vec_F / q2
		fx = probe_x - q.x;
		fy = probe_y - q.y;

		float l = length();
		if (l <= minDistance) {
			return;
		}

		float e = k * q.q / (l * l * l);
		if (e > maxForce) {
			fx = fy = maxForce;
		} else {
			fx *= e;
			fy *= e;
		}
	}

	__device__ force operator +(const force& f) const {
		return force(fx + f.fx, fy + f.fy);
	}

	__device__ force operator -(const force& f) const {
		return force(fx - f.fx, fy - f.fy);
	}

	__device__ force& operator +=(const force& f) {
		fx += f.fx;
		fy += f.fy;
		return *this;
	}

	__device__ force& operator -=(const force& f) {
		fx -= f.fx;
		fy -= f.fy;
		return *this;
	}
};

Задумка моя простая — я создам двухмерный массив из векторов напряжённостей. В каждой точке экрана, напряжённость будет равна сумме векторов кулоновских сил. Звучит заумно, но я иначе не смог сформулировать. Проще сказать, код будет вот таким (вы его прочтите, но не пишите ещё в свой проект):
	force temp_f;
	for (int i = 0; i < chargeCount; i++) {
		temp_f.calculate(charges[i], x, y);
		*f += temp_f;
	}

Если мы попадём в заряд, то я пропускаю вычисление напряженности в этой точке (ибо делить на ноль нельзя). Так что на экране вы должны увидеть чёрную точку — заряд.

Теперь, когда мы создали инфраструктуру — давайте кодить под CUDA. Для начала создадим необходимые переменные хоста и устройства и освободим их!

Создаем и освобождаем ресурсы


Нам понадобятся следующие переменные:
  • экран — хост — цветные пиксели;
  • экран — устройство — они же;
  • заряды — устройство — координаты и величины зарядов в постоянной памяти устройства.


Создаём:
uchar4* screen = NULL;

// ...

	screen = (uchar4*) malloc(width * height * sizeof(uchar4));
	memset(screen, 0, width * height * sizeof(uchar4));

	uchar4 *dev_screen = NULL;
	cudaMalloc((void**) &dev_screen, width * height * sizeof(uchar4));
	cudaMemset(dev_screen, 0, width * height * sizeof(uchar4));

// ...

	// Free resources
	free(screen);
	screen = NULL;

	cudaFree(dev_screen);
	dev_screen = NULL;

Заряды проинициализируем просто. Я не сбрасываю генератор случайных чисел, чтобы получить одну и ту же картину. Если вы хотите случайное поле каждый раз, то можете добавить вызов srand:
const int chargeCount = 10;
__constant__ charge dev_charges[chargeCount];
const int maxCharge = 1000;
const int minCharge = -1000;

// ...

void prepareCharges(void) {
	charge* charges = (charge*) malloc(chargeCount * sizeof(charge));
	for (int i = 0; i < chargeCount; i++) {
		charges[i].x = rand() % width;
		charges[i].y = rand() % height;
		charges[i].q = rand() % (maxCharge - minCharge) + minCharge;

		printf("Debug: Charge #%d (%d, %d, %d)\n", i, charges[i].x,
				charges[i].y, charges[i].q);
	}
	cudaMemcpyToSymbol(dev_charges, charges, chargeCount * sizeof(charge));
}

// ...
	prepareCharges();

Помните закомментированную строку, рендерящую экран? Теперь вы можете её раскомментировать:
	glDrawPixels(width, height, GL_RGBA, GL_UNSIGNED_BYTE, screen);

Создадим ядро


Теперь создадим функцию, которую мы будем вызывать в каждом потоке, чтобы срендерить каждый пиксель на экране. Эта функция будет выполнять следующее: она вычислит значение напряжённости поля, и вызовет другую функцию, чтобы сопоставить какой-нибудь цвет точке.

__device__ uchar4 getColor(const force& f) {
	uchar4 color;
	color.x = color.y = color.z = color.w = 0;

	float l = f.length();
	color.x = (l > maxLengthForColor ? 255 : l * 256 / maxLengthForColor);

	return color;
}

__global__ void renderFrame(uchar4* screen)
{
	// build the field and render the frame
}

Чтобы вызвать ядро… а впрочем и так всё ясно:
	// Launch Kernel to render the image
	renderFrame<<<blocks, threads>>>(dev_screen);

	cudaMemcpy(screen, dev_screen, width * height * sizeof(uchar4),
			cudaMemcpyDeviceToHost);

Теперь добавим немного кода в renderFrame и получим то, что нам нужно:
__global__ void renderFrame(uchar4* screen)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;

	force f, temp_f;
	for (int i = 0; i < chargeCount; i++) {
		temp_f.calculate(dev_charges[i], x, y);
		f += temp_f;
	}

	screen[x + y * width] = getColor(f);
}

Вывод


Часто народ пропускает статью, потому что важен вывод. Мне здесь нечего сказать, кроме как опубликовать картинку, что у меня получилась. Ссылку на почти пустой репозиторий найдёте здесь.

В принципе, есть куда развиваться — можно например добавить анимацию, показывающую как меняется поле при перемещении зарядов. Но опять же, если это интересно.
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 29

  • НЛО прилетело и опубликовало эту надпись здесь
      +1
      На Udacity есть курс дифференциальных прикладных уравнений, движение под действием гравитации они там тоже рассматривают.
      http://www.udacity.com/overview/Course/cs222/CourseRev/1
        +1
        Когда я ковырял простое выражение: F =ma = G*m*M / r^2, то картина получается гораздо проще, чем кажется — планета не падает на солнце, при некоторой минимальной скорости и ее начальном векторе. Согласно законам Кепплера, она вращается по обычным эллипсам.
        Если интересно, то могу написать об этом. Правда, физик из меня никакой :-) и любой человек из физмата меня сразу на место поставит.
        У меня на первой картинке точка движется по окружности с расчитанной скорость. из a = v^2 / r
        • НЛО прилетело и опубликовало эту надпись здесь
            0
            Даже если честно вычислять новую скорость, когда очередное смещение вдруг меняет радиус? Или лишь если считать скорость константой?
            • НЛО прилетело и опубликовало эту надпись здесь
                0
                При непрянолинейном движении ускорение есть всегда, меняется направление, а не длина вектора
                • НЛО прилетело и опубликовало эту надпись здесь
                    +1
                    Если модуль постоянен, то, как правильно заметил x0wl, может спокойно изменяться направление скорости, ибо она вектор . Второй закон Ньютона формулируется для векторов: (приобретаемое ускорение равно силе, делёной на массу). Проще говоря, куда сила, туда и ускорение. Бывает, что есть ускорение без силы, которое необходимо вводить, если мы из лабораторной системы переходим в ускоренно движущуюся систему отсчета (она будет неинерциальна). Чтобы соблюсти правило сила->ускорение, то вводят нереальную силу инерции. Именно она Вас толкает в метро, когда вагон разгоняется или тормозит, именно она придавливает в пол ускоренно поднимающегося лифта и т.п. Существует принцип эквивалентности Эйнштейна, гласящий, что локально сила инерции неотличима от гравитационной. Именно поэтому есть идея создавать искусственное тяготение на космическом корабле путем его вращения вокруг своей оси. Ведь, при вращении в с/о, связанной с вращающимся космическим кораблем, возникают центробежные ускорение и сила (эффективные), и центростремительное ускорение которое их уравновешивает и создается реальными силами (к примеру, трения или реакции опоры).
                    Если Вы говорите, что «скорость равна константе», то это может понимается, что сам вектор постоянен или только его длина (модуль скорости). Неплохо бы уточнять :)

                    Что касается описания движения взаимодействующих тел, самым простым является случай движения в заданном поле (к примеру, тяготения) — это задача Кеплера. Там три степени свободы (одно тело), одно векторное уравнение движения (три скалярных — на каждую координату). Решением уравнения движения является уравнение кривой второго порядка. То есть, траекторией небесного тела в поле тяготения одного лишь солнца может быть эллипс (планета), гипербола (комета) и парабола (камень брошен под углом к горизонту). Какая же будет траектория — зависит от величины эксцентриситета кривой в сравнении с 1. Сам эксцентриситет выражается через физические параметры системы. Подробнее смотрите в Ландау-Лившиц, том 1 или в Интернетах.
                    Задача двух тел — когда есть два свободных от внешних сил, но взаимодействующих друг с другом тела сравнимых масс — уже сложнее, но имеет аналитическое решение. Задача же трех и более тел аналитически пока не решается. Во всяком случае, для Кулоновского (гравитационного) потенциала.
                    Также я бы был очень аккуратен в плане применения метода Рунге-Кутта (в явной схеме) — он, емнип, с физической точки зрения не сохраняет энергию системы, и решать им уравнение, к примеру, колебаний очень плохо. Для кеплеровой задачи, весьма вероятно, тоже могут быть проблемы. Однако, когда-то еще в школе я пробовал рисовать на Pascal траекторию планеты самым очевидным методом решения уравнения движения: просто расписывая производные как отношение «дельт». Получался отличный эллипс которые ооочень медленно корежился. Как узнал уже в ВУЗе, я тогда применял метод Эйлера, а он еще проще метода Рунге-Кутта. Этим я хочу сказать, что физически задача движения тел во внешнем поле устойчива и по начальным значениям, и вообще :) Однако численный метод может внести сильные погрешности и неустойчивости.
                    • НЛО прилетело и опубликовало эту надпись здесь
                        +2
                        Уточнили. И я уточнил — вдруг кому будет полезно посмотреть на формулы :)

                        Все численные методы вносят ошибку. Это очевидно. Во-вторых, это ошибка машинного округления, но для типа double ею обычно можно уже пренебречь. А во-первых, это изменение исходной задачи в процессе применения численной схемы. Ведь, кода Вы применяете тот самый Рунге-Кутта метод, Вы а) вместо исходной аналитической функции рассматриваете сеточную функцию такую, что: . А это уже совсем другой объект с другими математическими свойствами; б) как следствие, Вы аппросимируете производные разностями (строите разностные схемы). Значение этих «разностей» уже отличаются от настоящей производной. Примером может служить уже упомянутое уравнение колебаний . Решением должен быть синус-косинус, но на практике явный метод Рунге-Кутта любого порядка точности даст нарастающую гармоническую функцию (к примеру, некоторое время назад этот эффект наблюдали тут). Дело в том, что при разностной аппроксимации Вы вносите в исходное уравнение член, пропорциональный . Это трение с отрицательным коэффициентом. То есть, описывает приток энергии. Вещь, конечно же, хорошо известная. Например, тут мужик хорошо это рассказывает. Какой бы численный метод решения чего бы то ни было Вы ни брали, всегда будет раздел «устойчивость метода». Она зависит от шага сетки, порядков аппроксимации и т.д. Если Вы выбрали параметры метода, не удовлетворяющие критериям устойчивости, то численное решение не будет сходится, пусть даже само решение исходной задачи есть и устойчиво во всех смыслах.
                        • НЛО прилетело и опубликовало эту надпись здесь
                  0
                  Конечно же, я говорил про модуль скорости.

                  А система явно должна иметь устойчивость, иначе бы солнечной не было б. Вы что-то упускаете в расчетах.
                  • НЛО прилетело и опубликовало эту надпись здесь
                      0
                      Отнюдь, таких систем миллиарды, это практически самая устойчивая из возможных конфигураций, самостабилизирующаяся к тому же. Это потенциальная яма, из которой очень затратно выбираться.
                        0
                        вот, комметарием ниже, я запостил видео, где без диффуров расчитывал вектора скоростей и ускорений. Вполне можно запустить планетку вращаться без всяких формул.
                0
                Если вам интересно, то я выложил флеш-видео моего прототипа, там как раз демонстрируется то, что вы говорили. Как раз для этого я и изучаю CUDA, чтобы попробовать придать моей математике на яваскипте чуть большую производительность.
                Видео 1 — уравновешенные четыре планеты
                Видео 2 — случайная, но удачная демонстрация предполагаемого апокалипсиса :)
                В случае с первым видео, там палочками отмечены вектора скоростей и ускорений. Кстати, вы увидите, что как по науке — ускорение направлено к центру окружности.

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

                p.s. если вам нужны сами SWF файлы с видюшками, могу выложить на хабрасторейдж (просто не хочу захламлять его лишний раз)
              0
              Вроде всего два диффура: \dot x и \ddot x. Это уже не говоря про то, что планеты сваливаться не будут.
              0
              Или в CUDA что-то изменилось, или люди по-прежнему забывают делать синхронизацию между запуском ядра и копированием результатов с видеокарты.
                0
                Так они же вроде в стандартный нулевой stream оба попадают и тем самым сериализуются?
                  0
                  странно, Джейсон Сандрерс не использовал cudaDeviceSynchronize функцию для синхронизации устройства. Единственное, что он применял — это синхронизацию тредов внутри одного блока.

                  Спасибо вам за заметку, погуглю на эту тему более детально, а то в инете пишут, что не для всех случаев (даже когда используется один поток) эту функцию следует использовать.
                  0
                  Мне кажетя в вашем случае нет нужды гнаться за большим размером блока, т.к. все равно не используется оптимизации вроде разбивки на tiles массива данных (да и вряд-ли тут она возможна), главное что-бы блок имел размер (т.е. общее число тредов) кратный 32, — размер warp-а, что-бы утилизировать ресурсы на полную.
                    0
                    скажите пожалуйста, из-за чего, я получал ошибку о нехватки ресурсов? из-за предельного размера блока, или нехватки регистров? и она возникала именно при выполнении операции деления, которая, как пишут в инете, тяжелая сама по себе.

                    Когда проект собирается в release режиме, я вижу, что количество используемых регистров ниже, чем максимальное на блок, но в этом случае такой ошибки нет. Забавный тот факт, что и в режиме отладки я использую максимум десять процентов от доступного числа.
                      0
                      Признаться — не знаю. Так-бы сразу написал ;) Я сам еще только изучаю CUDA по курсам на coursera.org.
                        0
                        Кроме "формального" ограничения на количество нитей в блоке, есть еще ограничение на количество памяти, отводимой на блок. Если одна нить использует много памяти, то запустить заявленное количество блоков (1024 в вашем случае) не получится.
                          0
                          Спасибо, что набрели на мою статью. Правильно ли я понимаю, что в режиме отладки количество используемых регистров было больше, чем в "релизной", чисто из-за отладочных целей?
                            0
                            В таких тонкостях я уже не силен, но вы скорее всего правы. Во время отладки видимо видеокарте требуется больше ресурсов. Но знаю точно, что реальное количество нитей в блоке надо определять, исходя из того, сколько памяти занимает одна нить. Это можно либо рассчитать, либо посмотреть в каком-нибудь отладчике типа Parallel Nsight или CUDA Profiler.
                      +2
                      Когда прозвучало слово «CUDA», я сразу подумал о решении уравнения Пуассона. А тут оказалось все в рамках школьной физики.
                        0
                        На самом деле, Nbody и SPH вполне популярны на CUDA.

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

                      Самое читаемое