Dev cuda

[Halfmuh]

Cuda realization of spatial mesh (SpatialMeshCu) and field solver (FieldSolver)
map_particle_charge methods commented (conflicting attempts to variables now stored on GPU memory)
todo:inner regions realization

[noooway]

(1 / (1 + is_inside_borders))
я правильно понимаю, что это должно быть 1/2 для точки внутри области и 1 для точки на границе?
У меня есть опасение, что сейчас из-за округления это будет либо 1, либо 0.

[noooway]

is_on_up_border = ((threadIdx.x == 0) && (blockIdx.x == 0));
is_on_low_border = ((threadIdx.x == (blockDim.x - 1)) && (blockIdx.x == (gridDim.x - 1)));

Для e.y эти строчки такие же, как и выше для случая e.x. И дальше для e.z тоже без изменений.
Это точно должно быть так? Не нужно поменять индексацию на что-нибудь типа ((threadIdx.y == 0) && (blockIdx.y == 0)) ?

Можно еще раз, как организована индексация?

Количество узлов в spat_mesh считывается из конфига.

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

dim3 SpatialMeshCu::GetThreads() {
	return dim3(16, 16, n_nodes.z / 16);
}

dim3 SpatialMeshCu::GetBlocks(dim3 nThreads) {
	return dim3(n_nodes.x / nThreads.x, n_nodes.y / nThreads.y, 16);

Дальше мы запускаем ядро с этим числом блоков/процессов на GPU.
Дальше каждый процесс обрабатывает только один "свой" элемент из spat_mesh.

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

__device__ int GetIdx() {
	//int xStepthread = 1;
	int xStepBlock = blockDim.x;
	int yStepThread = d_n_nodes[0].x;
	int yStepBlock = yStepThread * blockDim.y;
	int zStepThread = d_n_nodes[0].x * d_n_nodes[0].y;
	int zStepBlock = zStepThread * blockDim.z;
	return (threadIdx.x + blockIdx.x * xStepBlock) +
		(threadIdx.y * yStepThread + blockIdx.y * yStepBlock) +
		(threadIdx.z * zStepThread + blockIdx.z * zStepBlock);
 }

Т.е. за каждым GPU процессом закреплена своя точка из spat_mesh и пересчет из номера процесса в индекс spat_mesh можно делать по формуле spat_mesh_idx = threadIdx.x + blockIdx.x * blockSize.x (ну или что-то вроде) и аналогично для y и z компонент?

[noooway]

Все эти оптимизации это конечно замечательно. Но сейчас есть какой-то косяк с выставлением граничных условий. И я сходу не могу понять где именно.

[noooway]

Если я правильно понял исходную идею за этими функциями, то они должны вызываться либо с 2 потоками и 1 блоком по Z, либо с 1 потоком и 2 блоками.
Есть подозрение, что когда 2 и в блоках и в потоках - результат будет неверным.

threads = dim3(4, 4, 2);
blocks = dim3(n_nodes.x / 4, n_nodes.z / 4, 2);