[LbmBenchmarkKernelsPublic.git] / src / BenchKernelD3Q19List.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 "BenchKernelD3Q19ListCommon.h"

#include "Memory.h"
#include "Vtk.h"
#include "LikwidIf.h"

#include <inttypes.h>
#include <math.h>


void FNAME(D3Q19ListKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
{
	Assert(ld != NULL);
	Assert(kernelData != NULL);
	Assert(cd != NULL);

	Assert(cd->Omega > F(0.0));
	Assert(cd->Omega < F(2.0));

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

	PdfT omega = cd->Omega;
	PdfT omegaEven = omega;
	PdfT magicParam = F(1.0) / F(12.0);
	PdfT omegaOdd   = F(1.0) /(F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));

	PdfT evenPart = F(0.0);
	PdfT oddPart  = F(0.0);
	PdfT dir_indep_trm = F(0.0);

	PdfT w_0 = F(1.0) / F( 3.0);
	PdfT w_1 = F(1.0) / F(18.0);
	PdfT w_2 = F(1.0) / F(36.0);

	PdfT w_1_x3 = w_1 * F(3.0);	PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);	PdfT w_1_indep = F(0.0);
	PdfT w_2_x3 = w_2 * F(3.0);	PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);	PdfT w_2_indep = F(0.0);

	PdfT ux, uy, uz, ui;
	PdfT dens;

	// Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
	#define X(name, idx, idxinv, x, y, z)	PdfT JOIN(pdf_,name);
	D3Q19_LIST
	#undef X

	PdfT * src = kd->Pdfs[0];
	PdfT * dst = kd->Pdfs[1];
	PdfT * tmp;

	int maxIterations = cd->MaxIterations;
	int nFluid = kdl->nFluid;
	int nCells = kdl->nCells;

	uint32_t adjListIndex;
	uint32_t * adjList = kdl->AdjList;

	#ifdef VTK_OUTPUT
		if (cd->VtkOutput) {
			kd->PdfsActive = src;
			VtkWrite(ld, kd, cd, 0);
		}
	#endif

	#ifdef STATISTICS
		kd->PdfsActive = src;
		KernelStatistics(kd, ld, cd, 0);
	#endif

	// TODO: outer openmp parallel
	for(int iter = 0; iter < maxIterations; ++iter) {


	X_LIKWID_START("list-os");

	#ifdef _OPENMP
		#pragma omp parallel for default(none) \
				shared(nFluid, nCells, kd, kdl, adjList, src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \
				w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd) \
				private(ux, uy, uz, ui, dens, dir_indep_trm, adjListIndex, \
					pdf_C, \
				  	pdf_N, pdf_E, pdf_S, pdf_W, \
				  	pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
				  	pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
				  	pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \
				  	evenPart, oddPart, w_1_indep, w_2_indep)
	#endif
	#ifdef INTEL_OPT_DIRECTIVES
		#pragma ivdep
	#endif
	for (int index = 0; index < nFluid; ++index) {

			#define I(index, dir)	P_INDEX_3((nCells), (index), (dir))

#ifdef PROP_MODEL_PUSH

			// Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
			#define X(name, idx, idxinv, _x, _y, _z)	JOIN(pdf_,name) = src[I(index, idx)];
			D3Q19_LIST
			#undef X

#elif PROP_MODEL_PULL

			adjListIndex = index * N_D3Q19_IDX;

			pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)];

			// Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_N]]; ...
			#define X(name, idx, idxinv, _x, _y, _z)	JOIN(pdf_,name) = src[adjList[adjListIndex + idx]];
			D3Q19_LIST_WO_C
			#undef X

#else
	#error No implementation for PROP_MODEL_NAME.
#endif

// #define LID_DRIVEN_CAVITY

#ifdef LID_DRIVEN_CAVITY
			int nX = kd->Dims[0];
			int nY = kd->Dims[1];
			int nZ = kd->Dims[2];

			int x = kdl->Coords[C_INDEX_X(index)];
			int y = kdl->Coords[C_INDEX_Y(index)];
			int z = kdl->Coords[C_INDEX_Z(index)];

			if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
				ux = 0.1 * 0.577;
				uy = 0.0;
				uz = 0.0;
			} else {
#endif
				ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
					 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
				uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
					 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
				uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
					 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
#ifdef LID_DRIVEN_CAVITY
			}
#endif

			dens = pdf_C +
				   pdf_N  + pdf_E  + pdf_S  + pdf_W  +
				   pdf_NE + pdf_SE + pdf_SW + pdf_NW +
				   pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
				   pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;

			dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz) * F(3.0) / F(2.0);

#ifdef PROP_MODEL_PUSH

			adjListIndex = index * N_D3Q19_IDX;

			// direction: w_0
			dst[I(index, D3Q19_C)             ]  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);

			// direction: w_1
			w_1_indep = w_1*dir_indep_trm;

			ui = uy;
			evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
			dst[adjList[adjListIndex + D3Q19_N]]  = pdf_N - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_S]]  = pdf_S - evenPart + oddPart;

			ui = ux;
			evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
			dst[adjList[adjListIndex + D3Q19_E]]  = pdf_E - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_W]]  = pdf_W - evenPart + oddPart;

			ui = uz;
			evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
			dst[adjList[adjListIndex + D3Q19_T]]  = pdf_T - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_B]]  = pdf_B - evenPart + oddPart;

			// direction: w_2
			w_2_indep = w_2*dir_indep_trm;

			ui = -ux + uy;
			evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
			dst[adjList[adjListIndex + D3Q19_NW]] = pdf_NW - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_SE]] = pdf_SE - evenPart + oddPart;

			ui = ux + uy;
			evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
			dst[adjList[adjListIndex + D3Q19_NE]] = pdf_NE - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_SW]] = pdf_SW - evenPart + oddPart;

			ui = -ux + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
			dst[adjList[adjListIndex + D3Q19_TW]] = pdf_TW - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_BE]] = pdf_BE - evenPart + oddPart;

			ui = ux + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
			dst[adjList[adjListIndex + D3Q19_TE]] = pdf_TE - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_BW]] = pdf_BW - evenPart + oddPart;

			ui = -uy + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
			dst[adjList[adjListIndex + D3Q19_TS]] = pdf_TS - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_BN]] = pdf_BN - evenPart + oddPart;

			ui = uy + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
			dst[adjList[adjListIndex + D3Q19_TN]] = pdf_TN - evenPart - oddPart;
			dst[adjList[adjListIndex + D3Q19_BS]] = pdf_BS - evenPart + oddPart;

#elif PROP_MODEL_PULL

			// direction: w_0
			dst[I(index, D3Q19_C )]  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);

			// direction: w_1
			w_1_indep = w_1*dir_indep_trm;

			ui = uy;
			evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
			dst[I(index, D3Q19_N )]  = pdf_N - evenPart - oddPart;
			dst[I(index, D3Q19_S )]  = pdf_S - evenPart + oddPart;

			ui = ux;
			evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
			dst[I(index, D3Q19_E )]  = pdf_E - evenPart - oddPart;
			dst[I(index, D3Q19_W )]  = pdf_W - evenPart + oddPart;

			ui = uz;
			evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
			dst[I(index, D3Q19_T )]  = pdf_T - evenPart - oddPart;
			dst[I(index, D3Q19_B )]  = pdf_B - evenPart + oddPart;

			// direction: w_2
			w_2_indep = w_2*dir_indep_trm;

			ui = -ux + uy;
			evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
			dst[I(index, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
			dst[I(index, D3Q19_SE)] = pdf_SE - evenPart + oddPart;

			ui = ux + uy;
			evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
			dst[I(index, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
			dst[I(index, D3Q19_SW)] = pdf_SW - evenPart + oddPart;

			ui = -ux + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
			dst[I(index, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
			dst[I(index, D3Q19_BE)] = pdf_BE - evenPart + oddPart;

			ui = ux + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
			dst[I(index, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
			dst[I(index, D3Q19_BW)] = pdf_BW - evenPart + oddPart;

			ui = -uy + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
			dst[I(index, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
			dst[I(index, D3Q19_BN)] = pdf_BN - evenPart + oddPart;

			ui = uy + uz;
			evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
			oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
			dst[I(index, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
			dst[I(index, D3Q19_BS)] = pdf_BS - evenPart + oddPart;

#endif
			#undef I
		} // loop over fluid nodes

		X_LIKWID_STOP("list-os");

		#ifdef VERIFICATION
			kd->PdfsActive = dst;
			KernelAddBodyForce(kd, ld, cd);
		#endif

		#ifdef VTK_OUTPUT
			if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
				kd->PdfsActive = dst;
				VtkWrite(ld, kd, cd, iter);
			}
		#endif

		#ifdef STATISTICS
			kd->PdfsActive = dst;
			KernelStatistics(kd, ld, cd, iter);
		#endif

		// swap grids
		tmp = src;
		src = dst;
		dst = tmp;

	} // for (int iter = 0; ...

#ifdef VTK_OUTPUT
	if (cd->VtkOutput) {
		kd->PdfsActive = src;
		VtkWrite(ld, kd, cd, maxIterations);
	}
#endif

#ifdef STATISTICS
	kd->PdfsActive = src;
	KernelStatistics(kd, ld, cd, maxIterations);
#endif

	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 "BenchKernelD3Q19ListCommon.h"
	28
	29	#include "Memory.h"
	30	#include "Vtk.h"
e3f82424	31	#include "LikwidIf.h"
10988083 MW	32
	33	#include <inttypes.h>
	34	#include <math.h>
	35
	36
	37	void FNAME(D3Q19ListKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
	38	{
	39	Assert(ld != NULL);
	40	Assert(kernelData != NULL);
	41	Assert(cd != NULL);
	42
0fde6e45 MW	43	Assert(cd->Omega > F(0.0));
0fde6e45 MW	44	Assert(cd->Omega < F(2.0));
10988083 MW	45
	46	KernelData * kd = (KernelData *)kernelData;
	47	KernelDataList * kdl = (KernelDataList *)kernelData;
	48
	49	PdfT omega = cd->Omega;
	50	PdfT omegaEven = omega;
0fde6e45 MW	51	PdfT magicParam = F(1.0) / F(12.0);
0fde6e45 MW	52	PdfT omegaOdd = F(1.0) /(F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
10988083	53
0fde6e45 MW	54	PdfT evenPart = F(0.0);
	55	PdfT oddPart = F(0.0);
	56	PdfT dir_indep_trm = F(0.0);
10988083	57
0fde6e45 MW	58	PdfT w_0 = F(1.0) / F( 3.0);
	59	PdfT w_1 = F(1.0) / F(18.0);
	60	PdfT w_2 = F(1.0) / F(36.0);
10988083	61
0fde6e45 MW	62	PdfT w_1_x3 = w_1 * F(3.0); PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0); PdfT w_1_indep = F(0.0);
0fde6e45 MW	63	PdfT w_2_x3 = w_2 * F(3.0); PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0); PdfT w_2_indep = F(0.0);
10988083 MW	64
	65	PdfT ux, uy, uz, ui;
	66	PdfT dens;
	67
	68	// Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
	69	#define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
	70	D3Q19_LIST
	71	#undef X
	72
	73	PdfT * src = kd->Pdfs[0];
	74	PdfT * dst = kd->Pdfs[1];
	75	PdfT * tmp;
	76
	77	int maxIterations = cd->MaxIterations;
	78	int nFluid = kdl->nFluid;
	79	int nCells = kdl->nCells;
	80
	81	uint32_t adjListIndex;
	82	uint32_t * adjList = kdl->AdjList;
	83
	84	#ifdef VTK_OUTPUT
	85	if (cd->VtkOutput) {
	86	kd->PdfsActive = src;
	87	VtkWrite(ld, kd, cd, 0);
	88	}
	89	#endif
	90
	91	#ifdef STATISTICS
	92	kd->PdfsActive = src;
	93	KernelStatistics(kd, ld, cd, 0);
	94	#endif
	95
	96	// TODO: outer openmp parallel
	97	for(int iter = 0; iter < maxIterations; ++iter) {
	98
	99
e3f82424	100	X_LIKWID_START("list-os");
10988083 MW	101
	102	#ifdef _OPENMP
	103	#pragma omp parallel for default(none) \
	104	shared(nFluid, nCells, kd, kdl, adjList, src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \
	105	w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd) \
	106	private(ux, uy, uz, ui, dens, dir_indep_trm, adjListIndex, \
	107	pdf_C, \
	108	pdf_N, pdf_E, pdf_S, pdf_W, \
	109	pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
	110	pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
	111	pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \
	112	evenPart, oddPart, w_1_indep, w_2_indep)
	113	#endif
0fde6e45 MW	114	#ifdef INTEL_OPT_DIRECTIVES
	115	#pragma ivdep
	116	#endif
10988083 MW	117	for (int index = 0; index < nFluid; ++index) {
	118
	119	#define I(index, dir) P_INDEX_3((nCells), (index), (dir))
	120
	121	#ifdef PROP_MODEL_PUSH
	122
	123	// Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
	124	#define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(index, idx)];
	125	D3Q19_LIST
	126	#undef X
	127
	128	#elif PROP_MODEL_PULL
	129
	130	adjListIndex = index * N_D3Q19_IDX;
	131
	132	pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)];
	133
	134	// Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_N]]; ...
	135	#define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[adjList[adjListIndex + idx]];
	136	D3Q19_LIST_WO_C
	137	#undef X
	138
	139	#else
	140	#error No implementation for PROP_MODEL_NAME.
	141	#endif
	142
	143	// #define LID_DRIVEN_CAVITY
	144
	145	#ifdef LID_DRIVEN_CAVITY
	146	int nX = kd->Dims[0];
	147	int nY = kd->Dims[1];
	148	int nZ = kd->Dims[2];
	149
	150	int x = kdl->Coords[C_INDEX_X(index)];
	151	int y = kdl->Coords[C_INDEX_Y(index)];
	152	int z = kdl->Coords[C_INDEX_Z(index)];
	153
	154	if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
	155	ux = 0.1 * 0.577;
	156	uy = 0.0;
	157	uz = 0.0;
	158	} else {
	159	#endif
	160	ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
	161	pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
	162	uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
	163	pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
	164	uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
	165	pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
	166	#ifdef LID_DRIVEN_CAVITY
	167	}
	168	#endif
	169
	170	dens = pdf_C +
	171	pdf_N + pdf_E + pdf_S + pdf_W +
	172	pdf_NE + pdf_SE + pdf_SW + pdf_NW +
	173	pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
	174	pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
	175
0fde6e45	176	dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz) * F(3.0) / F(2.0);
10988083 MW	177
	178	#ifdef PROP_MODEL_PUSH
	179
	180	adjListIndex = index * N_D3Q19_IDX;
	181
	182	// direction: w_0
	183	dst[I(index, D3Q19_C) ] = pdf_C - omegaEven(pdf_C - w_0dir_indep_trm);
	184
	185	// direction: w_1
	186	w_1_indep = w_1*dir_indep_trm;
	187
	188	ui = uy;
0fde6e45 MW	189	evenPart = omegaEven( F(0.5)(pdf_N + pdf_S) - uiuiw_1_nine_half - w_1_indep );
0fde6e45 MW	190	oddPart = omegaOdd(F(0.5)(pdf_N - pdf_S) - ui*w_1_x3 );
10988083 MW	191	dst[adjList[adjListIndex + D3Q19_N]] = pdf_N - evenPart - oddPart;
	192	dst[adjList[adjListIndex + D3Q19_S]] = pdf_S - evenPart + oddPart;
	193
	194	ui = ux;
0fde6e45 MW	195	evenPart = omegaEven( F(0.5)(pdf_E + pdf_W) - uiuiw_1_nine_half - w_1_indep );
0fde6e45 MW	196	oddPart = omegaOdd(F(0.5)(pdf_E - pdf_W) - ui*w_1_x3 );
10988083 MW	197	dst[adjList[adjListIndex + D3Q19_E]] = pdf_E - evenPart - oddPart;
	198	dst[adjList[adjListIndex + D3Q19_W]] = pdf_W - evenPart + oddPart;
	199
	200	ui = uz;
0fde6e45 MW	201	evenPart = omegaEven( F(0.5)(pdf_T + pdf_B) - uiuiw_1_nine_half - w_1_indep );
0fde6e45 MW	202	oddPart = omegaOdd(F(0.5)(pdf_T - pdf_B) - ui*w_1_x3 );
10988083 MW	203	dst[adjList[adjListIndex + D3Q19_T]] = pdf_T - evenPart - oddPart;
	204	dst[adjList[adjListIndex + D3Q19_B]] = pdf_B - evenPart + oddPart;
	205
	206	// direction: w_2
	207	w_2_indep = w_2*dir_indep_trm;
	208
	209	ui = -ux + uy;
0fde6e45 MW	210	evenPart = omegaEven( F(0.5)(pdf_NW + pdf_SE) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	211	oddPart = omegaOdd(F(0.5)(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083 MW	212	dst[adjList[adjListIndex + D3Q19_NW]] = pdf_NW - evenPart - oddPart;
	213	dst[adjList[adjListIndex + D3Q19_SE]] = pdf_SE - evenPart + oddPart;
	214
	215	ui = ux + uy;
0fde6e45 MW	216	evenPart = omegaEven( F(0.5)(pdf_NE + pdf_SW) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	217	oddPart = omegaOdd(F(0.5)(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083 MW	218	dst[adjList[adjListIndex + D3Q19_NE]] = pdf_NE - evenPart - oddPart;
	219	dst[adjList[adjListIndex + D3Q19_SW]] = pdf_SW - evenPart + oddPart;
	220
	221	ui = -ux + uz;
0fde6e45 MW	222	evenPart = omegaEven( F(0.5)(pdf_TW + pdf_BE) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	223	oddPart = omegaOdd(F(0.5)(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083 MW	224	dst[adjList[adjListIndex + D3Q19_TW]] = pdf_TW - evenPart - oddPart;
	225	dst[adjList[adjListIndex + D3Q19_BE]] = pdf_BE - evenPart + oddPart;
	226
	227	ui = ux + uz;
0fde6e45 MW	228	evenPart = omegaEven( F(0.5)(pdf_TE + pdf_BW) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	229	oddPart = omegaOdd(F(0.5)(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083 MW	230	dst[adjList[adjListIndex + D3Q19_TE]] = pdf_TE - evenPart - oddPart;
	231	dst[adjList[adjListIndex + D3Q19_BW]] = pdf_BW - evenPart + oddPart;
	232
	233	ui = -uy + uz;
0fde6e45 MW	234	evenPart = omegaEven( F(0.5)(pdf_TS + pdf_BN) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	235	oddPart = omegaOdd(F(0.5)(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083 MW	236	dst[adjList[adjListIndex + D3Q19_TS]] = pdf_TS - evenPart - oddPart;
	237	dst[adjList[adjListIndex + D3Q19_BN]] = pdf_BN - evenPart + oddPart;
	238
	239	ui = uy + uz;
0fde6e45 MW	240	evenPart = omegaEven( F(0.5)(pdf_TN + pdf_BS) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	241	oddPart = omegaOdd(F(0.5)(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083 MW	242	dst[adjList[adjListIndex + D3Q19_TN]] = pdf_TN - evenPart - oddPart;
	243	dst[adjList[adjListIndex + D3Q19_BS]] = pdf_BS - evenPart + oddPart;
	244
	245	#elif PROP_MODEL_PULL
	246
	247	// direction: w_0
	248	dst[I(index, D3Q19_C )] = pdf_C - omegaEven(pdf_C - w_0dir_indep_trm);
	249
	250	// direction: w_1
	251	w_1_indep = w_1*dir_indep_trm;
	252
	253	ui = uy;
0fde6e45 MW	254	evenPart = omegaEven( F(0.5)(pdf_N + pdf_S) - uiuiw_1_nine_half - w_1_indep );
0fde6e45 MW	255	oddPart = omegaOdd(F(0.5)(pdf_N - pdf_S) - ui*w_1_x3 );
10988083 MW	256	dst[I(index, D3Q19_N )] = pdf_N - evenPart - oddPart;
	257	dst[I(index, D3Q19_S )] = pdf_S - evenPart + oddPart;
	258
	259	ui = ux;
0fde6e45 MW	260	evenPart = omegaEven( F(0.5)(pdf_E + pdf_W) - uiuiw_1_nine_half - w_1_indep );
0fde6e45 MW	261	oddPart = omegaOdd(F(0.5)(pdf_E - pdf_W) - ui*w_1_x3 );
10988083 MW	262	dst[I(index, D3Q19_E )] = pdf_E - evenPart - oddPart;
	263	dst[I(index, D3Q19_W )] = pdf_W - evenPart + oddPart;
	264
	265	ui = uz;
0fde6e45 MW	266	evenPart = omegaEven( F(0.5)(pdf_T + pdf_B) - uiuiw_1_nine_half - w_1_indep );
0fde6e45 MW	267	oddPart = omegaOdd(F(0.5)(pdf_T - pdf_B) - ui*w_1_x3 );
10988083 MW	268	dst[I(index, D3Q19_T )] = pdf_T - evenPart - oddPart;
	269	dst[I(index, D3Q19_B )] = pdf_B - evenPart + oddPart;
	270
	271	// direction: w_2
	272	w_2_indep = w_2*dir_indep_trm;
	273
	274	ui = -ux + uy;
0fde6e45 MW	275	evenPart = omegaEven( F(0.5)(pdf_NW + pdf_SE) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	276	oddPart = omegaOdd(F(0.5)(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083 MW	277	dst[I(index, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
	278	dst[I(index, D3Q19_SE)] = pdf_SE - evenPart + oddPart;
	279
	280	ui = ux + uy;
0fde6e45 MW	281	evenPart = omegaEven( F(0.5)(pdf_NE + pdf_SW) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	282	oddPart = omegaOdd(F(0.5)(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083 MW	283	dst[I(index, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
	284	dst[I(index, D3Q19_SW)] = pdf_SW - evenPart + oddPart;
	285
	286	ui = -ux + uz;
0fde6e45 MW	287	evenPart = omegaEven( F(0.5)(pdf_TW + pdf_BE) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	288	oddPart = omegaOdd(F(0.5)(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083 MW	289	dst[I(index, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
	290	dst[I(index, D3Q19_BE)] = pdf_BE - evenPart + oddPart;
	291
	292	ui = ux + uz;
0fde6e45 MW	293	evenPart = omegaEven( F(0.5)(pdf_TE + pdf_BW) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	294	oddPart = omegaOdd(F(0.5)(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083 MW	295	dst[I(index, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
	296	dst[I(index, D3Q19_BW)] = pdf_BW - evenPart + oddPart;
	297
	298	ui = -uy + uz;
0fde6e45 MW	299	evenPart = omegaEven( F(0.5)(pdf_TS + pdf_BN) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	300	oddPart = omegaOdd(F(0.5)(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083 MW	301	dst[I(index, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
	302	dst[I(index, D3Q19_BN)] = pdf_BN - evenPart + oddPart;
	303
	304	ui = uy + uz;
0fde6e45 MW	305	evenPart = omegaEven( F(0.5)(pdf_TN + pdf_BS) - uiuiw_2_nine_half - w_2_indep );
0fde6e45 MW	306	oddPart = omegaOdd(F(0.5)(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083 MW	307	dst[I(index, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
	308	dst[I(index, D3Q19_BS)] = pdf_BS - evenPart + oddPart;
	309
	310	#endif
	311	#undef I
	312	} // loop over fluid nodes
	313
e3f82424 MW	314	X_LIKWID_STOP("list-os");
e3f82424 MW	315
10988083 MW	316	#ifdef VERIFICATION
	317	kd->PdfsActive = dst;
	318	KernelAddBodyForce(kd, ld, cd);
	319	#endif
	320
	321	#ifdef VTK_OUTPUT
	322	if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
	323	kd->PdfsActive = dst;
	324	VtkWrite(ld, kd, cd, iter);
	325	}
	326	#endif
	327
	328	#ifdef STATISTICS
	329	kd->PdfsActive = dst;
	330	KernelStatistics(kd, ld, cd, iter);
	331	#endif
	332
	333	// swap grids
	334	tmp = src;
	335	src = dst;
	336	dst = tmp;
	337
	338	} // for (int iter = 0; ...
	339
	340	#ifdef VTK_OUTPUT
	341	if (cd->VtkOutput) {
	342	kd->PdfsActive = src;
	343	VtkWrite(ld, kd, cd, maxIterations);
	344	}
	345	#endif
	346
	347	#ifdef STATISTICS
	348	kd->PdfsActive = src;
	349	KernelStatistics(kd, ld, cd, maxIterations);
	350	#endif
	351
	352	return;
	353	}