[LbmBenchmarkKernelsPublic.git] / src / BenchKernelD3Q19ListAaCommon.c

// --------------------------------------------------------------------------
//
// Copyright
//   Markus Wittmann, 2016-2017
//   RRZE, University of Erlangen-Nuremberg, Germany
//   markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
//
//   Viktor Haag, 2016
//   LSS, University of Erlangen-Nuremberg, Germany
//
//  This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
//
//  LbmBenchKernels is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  LbmBenchKernels is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with LbmBenchKernels.  If not, see <http://www.gnu.org/licenses/>.
//
// --------------------------------------------------------------------------
#include "BenchKernelD3Q19ListAaCommon.h"

#include "Memory.h"
#include "Vtk.h"
#include "Padding.h"

#include <math.h>


// Forward definition.
void FNAME(D3Q19ListAaKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);


// -----------------------------------------------------------------------
// Functions which are used as callback by the kernel to read or write
// PDFs and nodes.

static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
{
	Assert(kd != NULL);
	Assert(kd->PdfsActive != NULL);
	Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
	Assert(pdf != NULL);

	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
	Assert(dir >= 0); Assert(dir < N_D3Q19);

	KernelDataList * kdl = (KernelDataList *)kd;

	if (kdl->Iteration % 2 == 0) {
		// Pdfs are stored inverse, local PDFs are located in remote nodes

		uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];

		if (dir != D3Q19_C) {
			uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;

			*pdf = kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + D3Q19_INV[dir]]];
		}
		else {
			*pdf = kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)];
		}

	}
	else {
		*pdf = kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)];
	}


	return;
}

static void FNAME(BCSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf)
{
	Assert(kd != NULL);
	Assert(kd->PdfsActive != NULL);
	Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
	Assert(dir >= 0); Assert(dir < N_D3Q19);

	if (isnan(pdf)) {
		printf("ERROR: setting nan %d %d %d %d %s\n", x, y, z, dir, D3Q19_NAMES[dir]);
		DEBUG_BREAK_POINT();
		exit(1);
	}

	KernelDataList * kdl = (KernelDataList *)kd;

	if (kdl->Iteration % 2 == 0) {
		// Pdfs are stored inverse, local PDFs are located in remote nodes

		uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];

		if (dir != D3Q19_C) {
			uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;

			kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + D3Q19_INV[dir]]] = pdf;
		}
		else {
			kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)] = pdf;
		}

	}
	else {
		kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)] = pdf;
	}

	return;
}


static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
{
	Assert(kd != NULL);
	Assert(kd->PdfsActive != NULL);
	Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
	Assert(pdfs != NULL);
	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);

	KernelDataList * kdl = (KernelDataList *)kd;

	if(kdl->Iteration % 2 == 0){

		uint32_t nodeIndex = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
		uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;

		// Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
		pdfs[D3Q19_C] = kd->PdfsActive[P_INDEX_3(kdl->nCells, nodeIndex, D3Q19_C)];

		#define X(name, idx, idxinv, _x, _y, _z)	pdfs[idx] = kd->PdfsActive[kdl->AdjList[adjListIndex + idxinv]];
		D3Q19_LIST_WO_C
		#undef X

	} else {

		#define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
		#define X(name, idx, idxinv, _x, _y, _z)	pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)];
		D3Q19_LIST
		#undef X
		#undef I

	}

#if 0
	// Detect NaNs
	for (int d = 0; d < 19; ++d) {
		if(isnan(pdfs[d]) || isinf(pdfs[d])) {
			printf("%d %d %d %d nan! get node\n", x, y, z, d);
						for (int d2 = 0; d2 < 19; ++d2) {
							printf("%d: %e\n", d2, pdfs[d2]);
						}
			exit(1);
		}
	}
#endif

	return;
}


static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
{
	Assert(kd != NULL);
	Assert(kd->PdfsActive != NULL);
	Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
	Assert(pdfs != NULL);

	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);

#if 0
	// Detect NaNs
	for (int d = 0; d < 19; ++d) {
		if(isnan(pdfs[d])) {
			printf("%d %d %d %d nan! get node\n", x, y, z, d);
						for (int d2 = 0; d2 < 19; ++d2) {
							printf("%d: %e\n", d2, pdfs[d2]);
						}
			exit(1);
		}
	}
#endif

	KernelDataList * kdl = (KernelDataList *)kd;

	if(kdl->Iteration % 2 == 0){

		uint32_t nodeIndex = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
		uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;

		// Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
		kd->PdfsActive[P_INDEX_3(kdl->nCells, nodeIndex, D3Q19_C)] = pdfs[D3Q19_C];

		#define X(name, idx, idxinv, _x, _y, _z)	kd->PdfsActive[kdl->AdjList[adjListIndex + idxinv]] = pdfs[idx];
		D3Q19_LIST_WO_C
		#undef X

	} else {

		#define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
		#define X(name, idx, idxinv, _x, _y, _z)	kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx];
		D3Q19_LIST
		#undef X
		#undef I

	}

	return;
}

static void ParameterUsage()
{
	printf("Kernel parameters:\n");
	printf("  [-blk <n>] [-blk-[xyz] <n>]\n");
#ifdef DATA_LAYOUT_SOA
	printf("  [-pad auto|modulus_1+offset_1(,modulus_n+offset_n)*]\n");
#endif
	return;
}

static void ParseParameters(Parameters * params, int * blk, PadInfo ** padInfo)
{
	Assert(blk != NULL);

	blk[0] = 0; blk[1] = 0; blk[2] = 0;
	*padInfo = NULL;

	#define ARG_IS(param)                   (!strcmp(params->KernelArgs[i], param))
	#define NEXT_ARG_PRESENT() \
		do { \
			if (i + 1 >= params->nKernelArgs) { \
				printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
				exit(1); \
			} \
		} while (0)


	for (int i = 0; i < params->nKernelArgs; ++i) {
		if (ARG_IS("-blk") || ARG_IS("--blk")) {
			NEXT_ARG_PRESENT();

			int tmp = strtol(params->KernelArgs[++i], NULL, 0);

			if (tmp < 0) {
				printf("ERROR: blocking parameter must be >= 0.\n");
				exit(1);
			}

			blk[0] = blk[1] = blk[2] = tmp;
		}
		else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) {
			NEXT_ARG_PRESENT();

			int tmp = strtol(params->KernelArgs[++i], NULL, 0);

			if (tmp < 0) {
				printf("ERROR: blocking parameter must be >= 0.\n");
				exit(1);
			}

			blk[0] = tmp;
		}
		else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) {
			NEXT_ARG_PRESENT();

			int tmp = strtol(params->KernelArgs[++i], NULL, 0);

			if (tmp < 0) {
				printf("ERROR: blocking parameter must be >= 0.\n");
				exit(1);
			}

			blk[1] = tmp;
		}
		else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) {
			NEXT_ARG_PRESENT();

			int tmp = strtol(params->KernelArgs[++i], NULL, 0);

			if (tmp < 0) {
				printf("ERROR: blocking parameter must be >= 0.\n");
				exit(1);
			}

			blk[2] = tmp;
		}
#ifdef DATA_LAYOUT_SOA
		else if (ARG_IS("-pad") || ARG_IS("--pad")) {
			NEXT_ARG_PRESENT();

			*padInfo = PadInfoFromStr(params->KernelArgs[++i]);
		}
#endif
		else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
			ParameterUsage();
			exit(1);
		}
		else {
			printf("ERROR: unknown kernel parameter.\n");
			ParameterUsage();
			exit(1);
		}
	}

	#undef ARG_IS
	#undef NEXT_ARG_PRESENT

	return;
}

void FNAME(D3Q19ListAaInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
{
	KernelData * kd;
	KernelDataList * kdl;
	MemAlloc((void **)&kdl, sizeof(KernelDataList));

	kd = (KernelData *)kdl;
	*kernelData = kd;

#ifdef DEBUG
	kd->Pdfs[0] = NULL;
	kd->Pdfs[1] = NULL;
	kd->PdfsActive = NULL;
	kd->DstPdfs = NULL;
	kd->SrcPdfs = NULL;
	kd->Dims[0] = -1;
	kd->Dims[1] = -1;
	kd->Dims[2] = -1;
	kd->GlobalDims[0] = -1;
	kd->GlobalDims[1] = -1;
	kd->GlobalDims[2] = -1;
	kd->Offsets[0] = -1;
	kd->Offsets[1] = -1;
	kd->Offsets[2] = -1;

	kd->ObstIndices = NULL;
	kd->nObstIndices = -1;
	kd->BounceBackPdfsSrc = NULL;
	kd->BounceBackPdfsDst = NULL;
	kd->nBounceBackPdfs = -1;

	kdl->AdjList = NULL;
	kdl->Coords = NULL;
	kdl->Grid = NULL;
	kdl->nCells = -1;
	kdl->nFluid = -1;
#endif

	int blk[3] = { 0 };
	PadInfo * padInfo = NULL;

	ParseParameters(params, blk, &padInfo);

	// Ajust the dimensions according to padding, if used.
	kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0];
	kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1];
	kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2];

	int * lDims = ld->Dims;

	int lX = lDims[0];
	int lY = lDims[1];
	int lZ = lDims[2];

	int nTotalCells = lX * lY * lZ;
	int nCells = ld->nFluid;
	int nFluid = ld->nFluid;

#ifdef DATA_LAYOUT_SOA
	{
		nCells = PadCellsAndReport(nCells, sizeof(PdfT), &padInfo);
		PadInfoFree(padInfo); padInfo = NULL;
	}
#endif

	// TODO: check nCells/nFluid do not exceed 2^31. This actually has to be
	// done during lattice setup.
	kdl->nCells = nCells;
	kdl->nFluid = nFluid;

	PdfT * pdfs[2];

	if (blk[0] == 0) blk[0] = lX;
	if (blk[1] == 0) blk[1] = lY;
	if (blk[2] == 0) blk[2] = lZ;

	printf("# blocking               x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);

	printf("# allocating data for %d fluid LB nodes with padding (%lu bytes = %f MiB for both lattices)\n",
		nCells, 2 * sizeof(PdfT) * nCells * N_D3Q19,
		2 * sizeof(PdfT) * nCells * N_D3Q19 / 1024.0 / 1024.0);

	MemAlloc((void **)&pdfs[0], sizeof(PdfT) * nCells * N_D3Q19);

	kd->Pdfs[0] = pdfs[0];

	// Initialize PDFs with some (arbitrary) data for correct NUMA placement.
	// Here we touch only the fluid nodes as this loop is OpenMP parallel and
	// we want the same scheduling as in the kernel.
	#ifdef _OPENMP
		#pragma omp parallel for
	#endif
	for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
		pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0;
	} }

	// Initialize all PDFs to some standard value.
	for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
		pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0;
	} }

	// ----------------------------------------------------------------------
	// create grid which will hold the index numbers of the fluid nodes

	uint32_t * grid;

	if (MemAlloc((void **)&grid, nTotalCells * sizeof(uint32_t))) {
		printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t));
		exit(1);
	}
	kdl->Grid = grid;

	int latticeIndex;

#ifdef DEBUG
	for(int z = 0; z < lZ; ++z) {
		for(int y = 0; y < lY; ++y) {
			for(int x = 0; x < lX; ++x) {

				latticeIndex = L_INDEX_4(ld->Dims, x, y, z);

				grid[latticeIndex] = ~0;
			}
		}
	}
#endif

	// ----------------------------------------------------------------------
	// generate numbering over grid

	uint32_t * coords;

	if (MemAlloc((void **)&coords, nFluid * sizeof(uint32_t) * 3)) {
		printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3);
		exit(1);
	}

	kdl->Coords = coords;

	// Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes
	// through the ld->Lattice array.
	int counter = 0;

	// Blocking is implemented via setup of the adjacency list. The kernel later will
	// walk through the lattice blocked automatically.
	for (int bX = 0; bX < lX; bX += blk[0]) {
	for (int bY = 0; bY < lY; bY += blk[1]) {
	for (int bZ = 0; bZ < lZ; bZ += blk[2]) {

		int eX = MIN(bX + blk[0], lX);
		int eY = MIN(bY + blk[1], lY);
		int eZ = MIN(bZ + blk[2], lZ);

		for (int x = bX; x < eX; ++x) {
		for (int y = bY; y < eY; ++y) {
		for (int z = bZ; z < eZ; ++z) {

			latticeIndex = L_INDEX_4(lDims, x, y, z);

			if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) {
				grid[latticeIndex] = counter;

				coords[C_INDEX_X(counter)] = x;
				coords[C_INDEX_Y(counter)] = y;
				coords[C_INDEX_Z(counter)] = z;

				++counter;
			}
		} } }
	} } }

	Verify(counter == nFluid);

	uint32_t * adjList;

	// AdjList only requires 18 instead of 19 entries per node, as
	// the center PDF needs no addressing.
	if (MemAlloc((void **)&adjList, nFluid * sizeof(uint32_t) * N_D3Q19_IDX)) {
		printf("ERROR: allocating adjList array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * N_D3Q19_IDX);
		exit(1);
	}

	kdl->AdjList = adjList;

	int x, y, z;

	uint32_t neighborIndex;
	uint32_t dstIndex;

	int nx, ny, nz, px, py, pz;

	// Loop over all fluid nodes and compute the indices to the neighboring
	// PDFs for configure data layout (AoS/SoA).
	#ifdef _OPENMP
		#pragma omp parallel for
	#endif
	for (int index = 0; index < nFluid; ++index) {
		for (int d = 0; d < N_D3Q19_IDX; ++d) {
			adjList[index * N_D3Q19_IDX + d] = -1;
		}
	}

	// #ifdef _OPENMP --> add line continuation
	// 	#pragma omp parallel for default(none)
	// 		shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z,
	// 				stderr,
	// 				lDims, grid, ld, lX, lY, lZ, adjList)
	// 		private(x, y, z, nx, ny, nz, neighborIndex, dstIndex)
	// #endif
	for (int index = 0; index < nFluid; ++index) {
		x = coords[C_INDEX_X(index)];
		y = coords[C_INDEX_Y(index)];
		z = coords[C_INDEX_Z(index)];

		Assert(x >= 0 && x < lX);
		Assert(y >= 0 && y < lY);
		Assert(z >= 0 && z < lZ);

		Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE);

		// Loop over all directions except the center one.
		for(int d = 0; d < N_D3Q19 - 1; ++d) {
			Assert(d != D3Q19_C);

#ifdef PROP_MODEL_PUSH
			nx = x + D3Q19_X[d];
			ny = y + D3Q19_Y[d];
			nz = z + D3Q19_Z[d];

#elif PROP_MODEL_PULL
			nx = x - D3Q19_X[d];
			ny = y - D3Q19_Y[d];
			nz = z - D3Q19_Z[d];
#else
			#error No implementation for this PROP_MODEL_NAME.
#endif
			// If the neighbor is outside the latcie in X direction and we have a
			// periodic boundary then we need to wrap around.
			if ( 	((nx < 0 || nx >= lX) && ld->PeriodicX) ||
					((ny < 0 || ny >= lY) && ld->PeriodicY) ||
					((nz < 0 || nz >= lZ) && ld->PeriodicZ)
																){
				// x periodic

				if (nx < 0) {
					px = lX - 1;
				}
				else if (nx >= lX) {
					px = 0;
				} else {
					px = nx;
				}
				// y periodic
				if (ny < 0) {
					py = lY - 1;
				}
				else if (ny >= lY) {
					py = 0;
				} else {
					py = ny;
				}

				// z periodic
				if (nz < 0) {
					pz = lZ - 1;
				}
				else if (nz >= lZ) {
					pz = 0;
				} else {
					pz = nz;
				}

				if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
					dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
				}
				else {
					neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)];

					AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d);

					dstIndex = P_INDEX_3(nCells, neighborIndex, d);
				}
			}
			else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
				dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
			}
			else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
				dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
			}
			else {
				neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)];

				Assert(neighborIndex != ~0);

				dstIndex = P_INDEX_3(nCells, neighborIndex, d);
			}

			Assert(dstIndex >= 0);
			Assert(dstIndex < nCells * N_D3Q19);

			adjList[index * N_D3Q19_IDX + d] = dstIndex;
		}
	}


	// Fill remaining KernelData structures
	kd->GetNode = GetNode;
	kd->SetNode = SetNode;

	kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf);
	kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf);

	kd->Kernel = FNAME(D3Q19ListAaKernel);

	kd->DstPdfs = NULL;
	kd->PdfsActive = kd->Pdfs[0];

	return;
}

void FNAME(D3Q19ListAaDeinit)(LatticeDesc * ld, KernelData ** kernelData)
{
	KernelDataList ** kdl = (KernelDataList **)kernelData;

	MemFree((void **)&((*kernelData)->Pdfs[0]));

	MemFree((void **)&((*kdl)->AdjList));
	MemFree((void **)&((*kdl)->Coords));
	MemFree((void **)&((*kdl)->Grid));

	MemFree((void **)kernelData);

	return;
}
Commit	Line	Data
10988083 MW	1	// --------------------------------------------------------------------------
	2	//
	3	// Copyright
	4	// Markus Wittmann, 2016-2017
	5	// RRZE, University of Erlangen-Nuremberg, Germany
	6	// markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
	7	//
	8	// Viktor Haag, 2016
	9	// LSS, University of Erlangen-Nuremberg, Germany
	10	//
	11	// This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
	12	//
	13	// LbmBenchKernels is free software: you can redistribute it and/or modify
	14	// it under the terms of the GNU General Public License as published by
	15	// the Free Software Foundation, either version 3 of the License, or
	16	// (at your option) any later version.
	17	//
	18	// LbmBenchKernels is distributed in the hope that it will be useful,
	19	// but WITHOUT ANY WARRANTY; without even the implied warranty of
	20	// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	21	// GNU General Public License for more details.
	22	//
	23	// You should have received a copy of the GNU General Public License
	24	// along with LbmBenchKernels. If not, see <http://www.gnu.org/licenses/>.
	25	//
	26	// --------------------------------------------------------------------------
	27	#include "BenchKernelD3Q19ListAaCommon.h"
	28
	29	#include "Memory.h"
	30	#include "Vtk.h"
e3f82424	31	#include "Padding.h"
10988083 MW	32
	33	#include <math.h>
	34
	35
	36	// Forward definition.
	37	void FNAME(D3Q19ListAaKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
	38
	39
	40
	41
	42	// -----------------------------------------------------------------------
	43	// Functions which are used as callback by the kernel to read or write
	44	// PDFs and nodes.
	45
	46	static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
	47	{
	48	Assert(kd != NULL);
	49	Assert(kd->PdfsActive != NULL);
	50	Assert(kd->PdfsActive == kd->Pdfs[0] \|\| kd->PdfsActive == kd->Pdfs[1]);
	51	Assert(pdf != NULL);
	52
	53	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
	54	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
	55	Assert(dir >= 0); Assert(dir < N_D3Q19);
	56
	57	KernelDataList * kdl = (KernelDataList *)kd;
	58
	59	if (kdl->Iteration % 2 == 0) {
	60	// Pdfs are stored inverse, local PDFs are located in remote nodes
	61
	62	uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];
	63
	64	if (dir != D3Q19_C) {
	65	uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
	66
	67	*pdf = kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + D3Q19_INV[dir]]];
	68	}
	69	else {
	70	*pdf = kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)];
	71	}
	72
	73	}
	74	else {
	75	*pdf = kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)];
	76	}
	77
	78
	79	return;
	80	}
	81
	82	static void FNAME(BCSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf)
	83	{
	84	Assert(kd != NULL);
	85	Assert(kd->PdfsActive != NULL);
	86	Assert(kd->PdfsActive == kd->Pdfs[0] \|\| kd->PdfsActive == kd->Pdfs[1]);
	87	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
	88	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
	89	Assert(dir >= 0); Assert(dir < N_D3Q19);
	90
	91	if (isnan(pdf)) {
	92	printf("ERROR: setting nan %d %d %d %d %s\n", x, y, z, dir, D3Q19_NAMES[dir]);
	93	DEBUG_BREAK_POINT();
	94	exit(1);
	95	}
96
97	KernelDataList * kdl = (KernelDataList *)kd;
98
99	if (kdl->Iteration % 2 == 0) {
100	// Pdfs are stored inverse, local PDFs are located in remote nodes
101
102	uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];
103
104	if (dir != D3Q19_C) {
105	uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
106
107	kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + D3Q19_INV[dir]]] = pdf;
108	}
109	else {
110	kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)] = pdf;
111	}
112
113	}
114	else {
115	kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)] = pdf;
116	}
117
118	return;
119	}
120
121
122	static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
123	{
124	Assert(kd != NULL);
125	Assert(kd->PdfsActive != NULL);
126	Assert(kd->PdfsActive == kd->Pdfs[0] \|\| kd->PdfsActive == kd->Pdfs[1]);
127	Assert(pdfs != NULL);
128	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
129	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
130
131	KernelDataList * kdl = (KernelDataList *)kd;
132
133	if(kdl->Iteration % 2 == 0){
134
135	uint32_t nodeIndex = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
136	uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
137
138	// Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
139	pdfs[D3Q19_C] = kd->PdfsActive[P_INDEX_3(kdl->nCells, nodeIndex, D3Q19_C)];
140
141	#define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[kdl->AdjList[adjListIndex + idxinv]];
142	D3Q19_LIST_WO_C
143	#undef X
144
145	} else {
146
147	#define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
148	#define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)];
149	D3Q19_LIST
150	#undef X
151	#undef I
152
153	}
154
155	#if 0
156	// Detect NaNs
157	for (int d = 0; d < 19; ++d) {
158	if(isnan(pdfs[d]) \|\| isinf(pdfs[d])) {
159	printf("%d %d %d %d nan! get node\n", x, y, z, d);
160	for (int d2 = 0; d2 < 19; ++d2) {
161	printf("%d: %e\n", d2, pdfs[d2]);
162	}
163	exit(1);
164	}
165	}
166	#endif
167
168	return;
169	}
170
171
172	static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
173	{
174	Assert(kd != NULL);
175	Assert(kd->PdfsActive != NULL);
176	Assert(kd->PdfsActive == kd->Pdfs[0] \|\| kd->PdfsActive == kd->Pdfs[1]);
177	Assert(pdfs != NULL);
178
179	Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
180	Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
181
182	#if 0
183	// Detect NaNs
184	for (int d = 0; d < 19; ++d) {
185	if(isnan(pdfs[d])) {
186	printf("%d %d %d %d nan! get node\n", x, y, z, d);
187	for (int d2 = 0; d2 < 19; ++d2) {
188	printf("%d: %e\n", d2, pdfs[d2]);
189	}
190	exit(1);
191	}
192	}
193	#endif
194
195	KernelDataList * kdl = (KernelDataList *)kd;
196
197	if(kdl->Iteration % 2 == 0){
198
199	uint32_t nodeIndex = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
200	uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
201
202	// Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
203	kd->PdfsActive[P_INDEX_3(kdl->nCells, nodeIndex, D3Q19_C)] = pdfs[D3Q19_C];
204
205	#define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[kdl->AdjList[adjListIndex + idxinv]] = pdfs[idx];
206	D3Q19_LIST_WO_C
207	#undef X
208
209	} else {
210
211	#define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
212	#define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx];
213	D3Q19_LIST
214	#undef X
215	#undef I
216
217	}
218
219	return;
220	}
221
222	static void ParameterUsage()
223	{
224	printf("Kernel parameters:\n");
225	printf(" [-blk <n>] [-blk-[xyz] <n>]\n");
e3f82424 MW	226	#ifdef DATA_LAYOUT_SOA
	227	printf(" [-pad auto\|modulus_1+offset_1(,modulus_n+offset_n)*]\n");
	228	#endif
10988083 MW	229	return;
	230	}
	231
e3f82424	232	static void ParseParameters(Parameters * params, int * blk, PadInfo ** padInfo)
10988083 MW	233	{
	234	Assert(blk != NULL);
	235
	236	blk[0] = 0; blk[1] = 0; blk[2] = 0;
e3f82424	237	*padInfo = NULL;
10988083 MW	238
	239	#define ARG_IS(param) (!strcmp(params->KernelArgs[i], param))
	240	#define NEXT_ARG_PRESENT() \
	241	do { \
	242	if (i + 1 >= params->nKernelArgs) { \
	243	printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
	244	exit(1); \
	245	} \
	246	} while (0)
	247
	248
	249	for (int i = 0; i < params->nKernelArgs; ++i) {
	250	if (ARG_IS("-blk") \|\| ARG_IS("--blk")) {
	251	NEXT_ARG_PRESENT();
	252
	253	int tmp = strtol(params->KernelArgs[++i], NULL, 0);
	254
e3f82424 MW	255	if (tmp < 0) {
e3f82424 MW	256	printf("ERROR: blocking parameter must be >= 0.\n");
10988083 MW	257	exit(1);
	258	}
	259
	260	blk[0] = blk[1] = blk[2] = tmp;
	261	}
	262	else if (ARG_IS("-blk-x") \|\| ARG_IS("--blk-x")) {
	263	NEXT_ARG_PRESENT();
	264
	265	int tmp = strtol(params->KernelArgs[++i], NULL, 0);
	266
e3f82424 MW	267	if (tmp < 0) {
e3f82424 MW	268	printf("ERROR: blocking parameter must be >= 0.\n");
10988083 MW	269	exit(1);
	270	}
	271
	272	blk[0] = tmp;
	273	}
	274	else if (ARG_IS("-blk-y") \|\| ARG_IS("--blk-y")) {
	275	NEXT_ARG_PRESENT();
	276
	277	int tmp = strtol(params->KernelArgs[++i], NULL, 0);
	278
e3f82424 MW	279	if (tmp < 0) {
e3f82424 MW	280	printf("ERROR: blocking parameter must be >= 0.\n");
10988083 MW	281	exit(1);
	282	}
	283
	284	blk[1] = tmp;
	285	}
	286	else if (ARG_IS("-blk-z") \|\| ARG_IS("--blk-z")) {
	287	NEXT_ARG_PRESENT();
	288
	289	int tmp = strtol(params->KernelArgs[++i], NULL, 0);
	290
e3f82424 MW	291	if (tmp < 0) {
e3f82424 MW	292	printf("ERROR: blocking parameter must be >= 0.\n");
10988083 MW	293	exit(1);
	294	}
	295
	296	blk[2] = tmp;
	297	}
e3f82424 MW	298	#ifdef DATA_LAYOUT_SOA
	299	else if (ARG_IS("-pad") \|\| ARG_IS("--pad")) {
	300	NEXT_ARG_PRESENT();
	301
	302	*padInfo = PadInfoFromStr(params->KernelArgs[++i]);
	303	}
	304	#endif
10988083 MW	305	else if (ARG_IS("-h") \|\| ARG_IS("-help") \|\| ARG_IS("--help")) {
	306	ParameterUsage();
	307	exit(1);
	308	}
	309	else {
	310	printf("ERROR: unknown kernel parameter.\n");
	311	ParameterUsage();
	312	exit(1);
	313	}
	314	}
	315
	316	#undef ARG_IS
	317	#undef NEXT_ARG_PRESENT
	318
	319	return;
	320	}
	321
	322	void FNAME(D3Q19ListAaInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
	323	{
	324	KernelData * kd;
	325	KernelDataList * kdl;
	326	MemAlloc((void **)&kdl, sizeof(KernelDataList));
	327
	328	kd = (KernelData *)kdl;
	329	*kernelData = kd;
	330
	331	#ifdef DEBUG
	332	kd->Pdfs[0] = NULL;
	333	kd->Pdfs[1] = NULL;
	334	kd->PdfsActive = NULL;
	335	kd->DstPdfs = NULL;
	336	kd->SrcPdfs = NULL;
	337	kd->Dims[0] = -1;
	338	kd->Dims[1] = -1;
	339	kd->Dims[2] = -1;
	340	kd->GlobalDims[0] = -1;
	341	kd->GlobalDims[1] = -1;
	342	kd->GlobalDims[2] = -1;
	343	kd->Offsets[0] = -1;
	344	kd->Offsets[1] = -1;
	345	kd->Offsets[2] = -1;
	346
	347	kd->ObstIndices = NULL;
	348	kd->nObstIndices = -1;
	349	kd->BounceBackPdfsSrc = NULL;
	350	kd->BounceBackPdfsDst = NULL;
	351	kd->nBounceBackPdfs = -1;
	352
	353	kdl->AdjList = NULL;
	354	kdl->Coords = NULL;
	355	kdl->Grid = NULL;
	356	kdl->nCells = -1;
	357	kdl->nFluid = -1;
	358	#endif
	359
e3f82424 MW	360	int blk[3] = { 0 };
	361	PadInfo * padInfo = NULL;
	362
	363	ParseParameters(params, blk, &padInfo);
	364
10988083 MW	365	// Ajust the dimensions according to padding, if used.
	366	kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0];
	367	kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1];
	368	kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2];
	369
	370	int * lDims = ld->Dims;
	371
	372	int lX = lDims[0];
	373	int lY = lDims[1];
	374	int lZ = lDims[2];
	375
	376	int nTotalCells = lX * lY * lZ;
e3f82424	377	int nCells = ld->nFluid;
10988083 MW	378	int nFluid = ld->nFluid;
10988083 MW	379
e3f82424 MW	380	#ifdef DATA_LAYOUT_SOA
	381	{
	382	nCells = PadCellsAndReport(nCells, sizeof(PdfT), &padInfo);
	383	PadInfoFree(padInfo); padInfo = NULL;
	384	}
	385	#endif
	386
8cafd9ea MW	387	// TODO: check nCells/nFluid do not exceed 2^31. This actually has to be
8cafd9ea MW	388	// done during lattice setup.
10988083 MW	389	kdl->nCells = nCells;
	390	kdl->nFluid = nFluid;
	391
	392	PdfT * pdfs[2];
	393
10988083 MW	394	if (blk[0] == 0) blk[0] = lX;
	395	if (blk[1] == 0) blk[1] = lY;
	396	if (blk[2] == 0) blk[2] = lZ;
	397
	398	printf("# blocking x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);
	399
	400	printf("# allocating data for %d fluid LB nodes with padding (%lu bytes = %f MiB for both lattices)\n",
	401	nCells, 2 * sizeof(PdfT) * nCells * N_D3Q19,
	402	2 * sizeof(PdfT) * nCells * N_D3Q19 / 1024.0 / 1024.0);
	403
	404	MemAlloc((void *)&pdfs[0], sizeof(PdfT) nCells * N_D3Q19);
	405
	406	kd->Pdfs[0] = pdfs[0];
	407
	408	// Initialize PDFs with some (arbitrary) data for correct NUMA placement.
	409	// Here we touch only the fluid nodes as this loop is OpenMP parallel and
	410	// we want the same scheduling as in the kernel.
	411	#ifdef _OPENMP
	412	#pragma omp parallel for
	413	#endif
	414	for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
	415	pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0;
	416	} }
	417
	418	// Initialize all PDFs to some standard value.
	419	for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
	420	pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0;
	421	} }
	422
	423	// ----------------------------------------------------------------------
	424	// create grid which will hold the index numbers of the fluid nodes
	425
	426	uint32_t * grid;
	427
	428	if (MemAlloc((void *)&grid, nTotalCells sizeof(uint32_t))) {
	429	printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t));
	430	exit(1);
	431	}
	432	kdl->Grid = grid;
	433
	434	int latticeIndex;
	435
	436	#ifdef DEBUG
	437	for(int z = 0; z < lZ; ++z) {
	438	for(int y = 0; y < lY; ++y) {
	439	for(int x = 0; x < lX; ++x) {
	440
	441	latticeIndex = L_INDEX_4(ld->Dims, x, y, z);
	442
	443	grid[latticeIndex] = ~0;
	444	}
	445	}
	446	}
	447	#endif
	448
	449	// ----------------------------------------------------------------------
	450	// generate numbering over grid
	451
	452	uint32_t * coords;
	453
	454	if (MemAlloc((void *)&coords, nFluid sizeof(uint32_t) * 3)) {
	455	printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3);
	456	exit(1);
	457	}
458
459	kdl->Coords = coords;
460
461	// Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes
462	// through the ld->Lattice array.
463	int counter = 0;
464
465	// Blocking is implemented via setup of the adjacency list. The kernel later will
466	// walk through the lattice blocked automatically.
10988083	467	for (int bX = 0; bX < lX; bX += blk[0]) {
e3f82424 MW	468	for (int bY = 0; bY < lY; bY += blk[1]) {
e3f82424 MW	469	for (int bZ = 0; bZ < lZ; bZ += blk[2]) {
10988083 MW	470
	471	int eX = MIN(bX + blk[0], lX);
	472	int eY = MIN(bY + blk[1], lY);
	473	int eZ = MIN(bZ + blk[2], lZ);
	474
10988083	475	for (int x = bX; x < eX; ++x) {
e3f82424 MW	476	for (int y = bY; y < eY; ++y) {
e3f82424 MW	477	for (int z = bZ; z < eZ; ++z) {
10988083 MW	478
	479	latticeIndex = L_INDEX_4(lDims, x, y, z);
	480
	481	if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) {
	482	grid[latticeIndex] = counter;
	483
	484	coords[C_INDEX_X(counter)] = x;
	485	coords[C_INDEX_Y(counter)] = y;
	486	coords[C_INDEX_Z(counter)] = z;
	487
	488	++counter;
	489	}
	490	} } }
	491	} } }
	492
	493	Verify(counter == nFluid);
	494
	495	uint32_t * adjList;
	496
	497	// AdjList only requires 18 instead of 19 entries per node, as
	498	// the center PDF needs no addressing.
	499	if (MemAlloc((void *)&adjList, nFluid sizeof(uint32_t) * N_D3Q19_IDX)) {
	500	printf("ERROR: allocating adjList array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * N_D3Q19_IDX);
	501	exit(1);
	502	}
	503
	504	kdl->AdjList = adjList;
	505
	506	int x, y, z;
	507
	508	uint32_t neighborIndex;
	509	uint32_t dstIndex;
	510
	511	int nx, ny, nz, px, py, pz;
	512
	513	// Loop over all fluid nodes and compute the indices to the neighboring
	514	// PDFs for configure data layout (AoS/SoA).
e3f82424 MW	515	#ifdef _OPENMP
	516	#pragma omp parallel for
	517	#endif
	518	for (int index = 0; index < nFluid; ++index) {
	519	for (int d = 0; d < N_D3Q19_IDX; ++d) {
	520	adjList[index * N_D3Q19_IDX + d] = -1;
	521	}
	522	}
	523
10988083 MW	524	// #ifdef _OPENMP --> add line continuation
	525	// #pragma omp parallel for default(none)
	526	// shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z,
	527	// stderr,
	528	// lDims, grid, ld, lX, lY, lZ, adjList)
	529	// private(x, y, z, nx, ny, nz, neighborIndex, dstIndex)
	530	// #endif
	531	for (int index = 0; index < nFluid; ++index) {
	532	x = coords[C_INDEX_X(index)];
	533	y = coords[C_INDEX_Y(index)];
	534	z = coords[C_INDEX_Z(index)];
	535
	536	Assert(x >= 0 && x < lX);
	537	Assert(y >= 0 && y < lY);
	538	Assert(z >= 0 && z < lZ);
	539
	540	Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE);
	541
	542	// Loop over all directions except the center one.
	543	for(int d = 0; d < N_D3Q19 - 1; ++d) {
	544	Assert(d != D3Q19_C);
	545
	546	#ifdef PROP_MODEL_PUSH
	547	nx = x + D3Q19_X[d];
	548	ny = y + D3Q19_Y[d];
	549	nz = z + D3Q19_Z[d];
	550
	551	#elif PROP_MODEL_PULL
	552	nx = x - D3Q19_X[d];
	553	ny = y - D3Q19_Y[d];
	554	nz = z - D3Q19_Z[d];
	555	#else
	556	#error No implementation for this PROP_MODEL_NAME.
	557	#endif
	558	// If the neighbor is outside the latcie in X direction and we have a
	559	// periodic boundary then we need to wrap around.
	560	if ( ((nx < 0 \|\| nx >= lX) && ld->PeriodicX) \|\|
	561	((ny < 0 \|\| ny >= lY) && ld->PeriodicY) \|\|
	562	((nz < 0 \|\| nz >= lZ) && ld->PeriodicZ)
	563	){
	564	// x periodic
	565
	566	if (nx < 0) {
	567	px = lX - 1;
	568	}
	569	else if (nx >= lX) {
	570	px = 0;
	571	} else {
	572	px = nx;
	573	}
	574	// y periodic
	575	if (ny < 0) {
	576	py = lY - 1;
	577	}
	578	else if (ny >= lY) {
	579	py = 0;
	580	} else {
	581	py = ny;
	582	}
	583
	584	// z periodic
	585	if (nz < 0) {
	586	pz = lZ - 1;
	587	}
588	else if (nz >= lZ) {
589	pz = 0;
590	} else {
591	pz = nz;
592	}
593
594	if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
595	dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
596	}
597	else {
598	neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)];
599
600	AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d);
601
602	dstIndex = P_INDEX_3(nCells, neighborIndex, d);
603	}
604	}
605	else if (nx < 0 \|\| ny < 0 \|\| nz < 0 \|\| nx >= lX \|\| ny >= lY \|\| nz >= lZ) {
606	dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
607	}
608	else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
609	dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
610	}
611	else {
612	neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)];
613
614	Assert(neighborIndex != ~0);
615
616	dstIndex = P_INDEX_3(nCells, neighborIndex, d);
617	}
618
619	Assert(dstIndex >= 0);
620	Assert(dstIndex < nCells * N_D3Q19);
621
622	adjList[index * N_D3Q19_IDX + d] = dstIndex;
623	}
624	}
625
626
627	// Fill remaining KernelData structures
628	kd->GetNode = GetNode;
629	kd->SetNode = SetNode;
630
631	kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf);
632	kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf);
633
634	kd->Kernel = FNAME(D3Q19ListAaKernel);
635
636	kd->DstPdfs = NULL;
637	kd->PdfsActive = kd->Pdfs[0];
638
639	return;
640	}
641
642	void FNAME(D3Q19ListAaDeinit)(LatticeDesc * ld, KernelData ** kernelData)
643	{
644	KernelDataList kdl = (KernelDataList )kernelData;
645
646	MemFree((void *)&((kernelData)->Pdfs[0]));
647
648	MemFree((void *)&((kdl)->AdjList));
649	MemFree((void *)&((kdl)->Coords));
650	MemFree((void *)&((kdl)->Grid));
651
652	MemFree((void **)kernelData);
653
654	return;
655	}
656