[LbmBenchmarkKernelsPublic.git] / src / Main.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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>  // strcasecmp

#include <inttypes.h>

#ifdef _OPENMP
#include <omp.h>
#endif

#include "Base.h"
#include "Kernel.h"
#include "Memory.h"

#include "Lattice.h"
#include "Geometry.h"
#include "Pinning.h"
#include "LikwidIf.h"

#include "KernelFunctions.h"

#ifdef __x86_64__
	#include <xmmintrin.h>


	#define MXCSR_DAZ					6
	#define MXCSR_FTZ 					15


	int FpIsMxCsrMaskSet(unsigned int mask)
	{
		unsigned int mxcsr;
		unsigned int mxcsrNew;

		mxcsr = _mm_getcsr();

		mxcsrNew = mxcsr & mask;

		return (mxcsrNew == mask);
	}

	int FpGetFtz()
	{
		return FpIsMxCsrMaskSet(1 << MXCSR_FTZ);
	}

	int FpGetDaz()
	{
		return FpIsMxCsrMaskSet(1 << MXCSR_DAZ);
	}
#endif


int ParseDimensions(const char * parameter, int * nX, int * nY, int * nZ)
{
	char * tmp;

	*nX = atoi(parameter);

	if (*nX <= 0) {
		printf("ERROR: parameter for X dimension must be > 0.\n");
		return 0;
	}

	tmp = strchr(parameter, 'x');

	if (tmp == NULL) {
		printf("ERROR: parameter for Y dimension is missing.\n");
		return 0;
	}

	*nY = atoi(tmp + 1);

	if (*nY <= 0) {
		printf("ERROR: parameter for Y dimension must be > 0.\n");
		return 0;
	}

	tmp = strchr(tmp + 1, 'x');

	if (tmp == NULL) {
		printf("ERROR: parameter for Z dimension is missing.\n");
		return 0;
	}

	*nZ = atoi(tmp + 1);

	if (*nZ <= 0) {
		printf("ERROR: parameter for Z dimension must be > 0.\n");
		return 0;
	}

	return 1;
}

int main(int argc, char * argv[])
{
	int dims[3] = { 20, 20, 20 };		// Dimensions in x, y, and z direction
	const char * geometryType = "channel";
	// int latticeDumpAscii = 0;
	int verify = 0; UNUSED(verify);
	char * kernelToUse = NULL;
	int nThreads = 1;
	const char * pinString = NULL;
	int periodic[3] = { 0 };

	CaseData cd;

	cd.MaxIterations 		= 1000;
	cd.RhoIn		 		= 1.0;
	cd.RhoOut 		 		= 1.0;
	cd.Omega		 		= 1.0;
	cd.VtkOutput	 		= 0;
	cd.VtkModulus    		= 100;
	cd.StatisticsModulus 	= 100;
	cd.XForce				= 0.00001;
	kernelToUse				= "push-soa";

	Parameters p;
	p.nArgs        = argc;
	p.Args         = argv;
	p.nKernelArgs  = 0;
	p.KernelArgs   = NULL;

#define LBM_BENCH_KERNELS_VERSION_MAJOR 0
#define LBM_BENCH_KERNELS_VERSION_MINOR 1

    printf("Lattice Boltzmann Benchmark Kernels (LbmBenchKernels) Copyright (C) 2016, 2017 LSS, RRZE\n");
    printf("This program comes with ABSOLUTELY NO WARRANTY; for details see LICENSE.\n");
	printf("This is free software, and you are welcome to redistribute it under certain conditions.\n");
	printf("\n");
	printf("LBM Benchmark Kernels %d.%d, compiled %s %s, type: %s\n",
		LBM_BENCH_KERNELS_VERSION_MAJOR, LBM_BENCH_KERNELS_VERSION_MINOR, __DATE__, __TIME__,
#ifdef VERIFICATION
		"verification"
#else
		"benchmark"
#endif
	);

	// ----------------------------------------------------------------------
	// Parse command line arguments

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

    for (int i = 1; i < argc; ++i) {

        if (ARG_IS("-dims") || ARG_IS("--dims")) {
            NEXT_ARG_PRESENT();


            if (!ParseDimensions(argv[++i], &dims[0], &dims[1], &dims[2])) {
                return 1;
            }
        }
		// else if (ARG_IS("-lattice-dump-ascii") || ARG_IS("--lattice-dump-ascii")) {
		// 	latticeDumpAscii = 1;
		// }
		else if (ARG_IS("-geometry") || ARG_IS("--geometry")) {
			NEXT_ARG_PRESENT();

			geometryType = argv[++i];
		}
		else if (ARG_IS("-iterations") ||ARG_IS("--iterations")) {
			NEXT_ARG_PRESENT();

			cd.MaxIterations = strtol(argv[++i], NULL, 0);

			if (cd.MaxIterations <= 0) {
				printf("ERROR: number of iterations must be > 0.\n");
				return 1;
			}
		}
		else if (ARG_IS("-rho-in") ||ARG_IS("--rho-in")) {
			NEXT_ARG_PRESENT();

			cd.RhoIn = strtod(argv[++i], NULL);
		}
		else if (ARG_IS("-rho-out") ||ARG_IS("--rho-out")) {
			NEXT_ARG_PRESENT();

			cd.RhoOut = strtod(argv[++i], NULL);
		}
		else if (ARG_IS("-omega") ||ARG_IS("--omega")) {
			NEXT_ARG_PRESENT();

			cd.Omega = strtod(argv[++i], NULL);
		}
		else if (ARG_IS("-x-force") ||ARG_IS("--x-force")) {
			NEXT_ARG_PRESENT();

			cd.XForce = strtod(argv[++i], NULL);
		}
		else if (ARG_IS("-verify") || ARG_IS("--verify")) {
#ifdef VERIFICATION

			// Choose this preset for verification. As geometry type "box" is
			// used but x and y direction are made pridoc.
			// Everything else can be altered, but enough iterations should be
			// performed in order to receive a fully developed flow field.
			verify = 1;

			cd.Omega  = 1.0;
			cd.RhoIn  = 1.0;
			cd.RhoOut = 1.0;
			geometryType = "box";
			dims[0] = 16;
			dims[1] = 16;
			dims[2] = 16;
			cd.XForce = 0.00001;
			cd.MaxIterations = 1000;
			periodic[0] = 1;
			periodic[1] = 1;
			periodic[2] = 0;

			printf("#\n");
			printf("# VERIFICATION: verifying flow profile of channel flow.\n");
			printf("#\n");

			// TODO: this is not a good idea as we ignore all other options...

#else
			printf("ERROR: in order to use -verify VERIFICATION must be defined during compilation.\n");
			printf("       Recompile with VERIFICATION=on.\n");
			return 1;
#endif
		}
		else if (ARG_IS("-vtk") || ARG_IS("--vtk")) {
#ifdef VTK_OUTPUT

			cd.VtkOutput = 1;

			// If the next parameter is a number it is used as the itartion count,
			// if not it is probably another parameter.
			if (i + 1 < argc) {

				int vtkModulus = atoi(argv[i+1]);

				if (vtkModulus > 0) {
					cd.VtkModulus = vtkModulus;
					++i;
				}
			}
#else
			printf("ERROR: in order to use -vtk VTK_OUTPUT must be defined during compilation.\n");
			printf("       Recompile with VTK_OUTPUT=on.\n");
			return 1;
#endif
		}
		else if (ARG_IS("-statistics") || ARG_IS("--statistics")) {
#ifdef STATISTICS
			NEXT_ARG_PRESENT();

			cd.StatisticsModulus = atoi(argv[++i]);

			if (cd.StatisticsModulus <= 0) {
				printf("ERROR: the iteration count for -statistics must be > 0.\n");
				return 1;
			}
#else
			printf("ERROR: in order to use -statistics STATISTICS must be defined during compilation.\n");
			printf("       Recompile with STATISTICS=on.\n");
			return 1;
#endif
		}
		else if (ARG_IS("-kernel") || ARG_IS("--kernel")) {
			NEXT_ARG_PRESENT();

			kernelToUse = argv[++i];
		}
		else if (ARG_IS("-list") || ARG_IS("--list")) {
			printf("Available kernels to benchmark:\n");

			for (int j = 0; j < N_ELEMS(g_kernels); ++j) {
				printf("   %s\n", g_kernels[j].Name);
			}

			return 0;
		}
		else if (ARG_IS("-pin") || ARG_IS("--pin")) {
			NEXT_ARG_PRESENT();

			pinString = argv[++i];
		}
		else if (ARG_IS("-t") || ARG_IS("-threads") || ARG_IS("--threads")) {
#ifdef _OPENMP
			NEXT_ARG_PRESENT();

			nThreads = atoi(argv[++i]);

			if (nThreads <= 0) {
				printf("ERROR: number of threads must be > 0.\n");
				return 1;
			}
#else
			printf("ERROR: specifying number of threads is only available when compiled with OpenMP support.\n");
			return 1;
#endif
		}
		else if (ARG_IS("-periodic-x") || ARG_IS("--periodic-x")) {
			periodic[0] = 1;
		}
		else if (ARG_IS("-periodic-y") || ARG_IS("--periodic-y")) {
			periodic[1] = 1;
		}
		else if (ARG_IS("-periodic-z") || ARG_IS("--periodic-z")) {
			periodic[2] = 1;
		}
		else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
			printf("ERROR: unknown argument: %s\n", argv[i]);
			printf("\n");
			printf("Usage:\n");
			printf("./lbmbenchk -list\n");
			printf("./lbmbenchk \n");
			printf("      [-dims XxYyZ] [-geometry box|channel|pipe|porosity[-value]] [-iterations <iterations>] [-lattice-dump-ascii]\n");
			printf("      [-rho-in <density>] [-rho-out <density] [-omega <omega>] [-kernel <kernel>]\n");
			printf("      [-periodic-x]\n");
#ifdef STATISTICS
			printf("      [-statistics <every-n-iteration>]\n");
#endif
#ifdef VTK_OUTPUT
			printf("      [-vtk [<every-n-iteration>]]\n");
#endif
#ifdef _OPENMP
			printf("      [-t <number of threads>]\n");
#endif
			printf("      [-pin core{,core}*]\n");
#ifdef VERIFICATION
			printf("      [-verify]\n");
#endif
			printf("      -- <kernel specific parameters>\n");
			printf("\n");
			printf("-list           List available kernels.\n");
			printf("\n");
			printf("-dims XxYxZ		Specify geometry dimensions.\n");
			printf("\n");
			printf("-geometry porosity-<value>\n");
			printf("                Geometetry with blocks of size <value> regularily layout out.\n");
			printf("\n");
			return 1;
		}
		else if (ARG_IS("--")) {
			// printf("# kernel args start with %s  these are %d args.\n", argv[i + 1], argc - i - 1);
			p.KernelArgs = &argv[++i];
			p.nKernelArgs = argc - i;
			break;
		}
		else {
			printf("ERROR: unknown parameter: %s.\n", argv[i]);
			exit(1);
		}
	}

	#undef ARG_IS
	#undef NEXT_ARG_PRESENT


	// ----------------------------------------------------------------------
	// Check if we exceed our index addressing PDFs.

	{
		uint64_t nPdfs = ((uint64_t)19) * dims[0] * dims[1] * dims[2];

		if (nPdfs > ((2LU << 31) - 1)) {
			printf("ERROR: number of PDFs exceed 2^31.\n");
			exit(1);
		}
	}

	// ----------------------------------------------------------------------

#ifdef _OPENMP
	omp_set_num_threads(nThreads);
#endif

 	LatticeDesc ld;

	GeoCreateByStr(geometryType, dims, periodic, &ld);

	const char * defines[] = {
#ifdef VTK_OUTPUT
	"VTK_OUTPUT",
#endif
#ifdef STATISTICS
	"STATISTICS",
#endif
#ifdef VERIFICATION
	"VERIFICATION",
#endif
#ifdef _OPENMP
	"_OPENMP",
#endif
#ifdef HAVE_LIKWID
	"HAVE_LIKWID",
#endif
	};

	printf("#  defines:           ");
	for (int j = 0; j < N_ELEMS(defines); ++j) {
		printf("%s ", defines[j]);
	}
	printf("\n");

	printf("#  nodes total:       % 10d\n", ld.nObst + ld.nFluid);
	printf("#  nodes fluid:       % 10d (including inlet & outlet)\n", ld.nFluid);
	printf("#  nodes obstacles:   % 10d\n", ld.nObst);
	printf("#  nodes inlet:       % 10d\n", ld.nInlet);
	printf("#  nodes outlet:      % 10d\n", ld.nOutlet);
	printf("#  periodicity:       x: %d y: %d z: %d\n", ld.PeriodicX, ld.PeriodicY, ld.PeriodicZ);

#ifdef VTK_OUTPUT
	printf("#  VTK output:         %d (every %d iteration)\n", cd.VtkOutput, cd.VtkModulus);
#endif
#ifdef STATISTICS
	printf("#  statistics:         every %d iteration\n", cd.StatisticsModulus);
#endif

	printf("#  omega:             %f\n", cd.Omega);
	printf("#  initial density at inlet/outlet:\n");
	printf("#    rho in:          %e\n", cd.RhoIn);
    printf("#    rho out:         %e\n", cd.RhoOut);
	printf("# iterations:         %d\n", cd.MaxIterations);

#ifdef __x86_64__
	printf("# fp status:          DAZ: %d  FTZ: %d\n", FpGetDaz(), FpGetFtz());
#endif

#ifdef _OPENMP
	printf("# OpenMP threads:     %d\n", omp_get_max_threads());

	if (pinString != NULL) {
		#pragma omp parallel
		{
			int threadId = omp_get_thread_num();
			int err;

			err = PinCurrentThreadByCpuList(pinString, 0, 0, threadId);

			if (err) {
				printf("ERROR [thread %d]: pinning failed.\n", threadId);
				exit(1);
			}

			const char * cpuList = PinCpuListAsString();
			Assert(cpuList != NULL);

			// Not so nice hack to print the thread ids ordered.
			#pragma omp for ordered
			for (int i = 0; i < omp_get_num_threads(); ++i) {
				#pragma omp ordered
				printf("# thread %2d  pinned to core(s):  %s\n", threadId, cpuList);
			}

			free((void *)cpuList);
		}
	}
#endif

	KernelData * kd;

	KernelFunctions * kf = NULL;

	if (kernelToUse == NULL) {
		kf = &g_kernels[0];
	}
	else {
		for (int j = 0; j < N_ELEMS(g_kernels); ++j) {

			if (!strcasecmp(kernelToUse, g_kernels[j].Name)) {
				kf = &g_kernels[j];
				break;
			}
		}
	}

	if (kf == NULL) {
		printf("ERROR: requested kernel \"%s\" not found.\n", kernelToUse);
		exit(1);
	}

	printf("#\n");
	printf("# kernel:           %s\n", kf->Name);
	printf("#\n");

	// Initialize kernel by calling its own initialization function
	kf->Init(&ld, &kd, &p);

#ifdef VERIFICATION
	if (verify) {
		KernelSetInitialDensity( &ld, kd, &cd);
		KernelSetInitialVelocity(&ld, kd, &cd);
	}
#endif

	printf("# starting kernel...\n");

	X_LIKWID_INIT();

	double timeStart = Time();

	// Call the LBM kernel
	kd->Kernel(&ld, kd, &cd);

	double duration = Time() - timeStart;

	X_LIKWID_DEINIT();

	// Print some statistics...
	KernelStatisticsAdv(kd, &ld, &cd, cd.MaxIterations, 1 /* force output */);

#ifdef VERIFICATION
	PdfT errorNorm = -1.0;
	KernelVerifiy(&ld, kd, &cd, &errorNorm);
#endif

	// Deinitialize kernel by calling its own deinitialization function
	kf->Deinit(&ld, &kd);


	double perf = (double)ld.nFluid * (double)cd.MaxIterations / duration / 1.e6;

	printf("P:   %f MFLUP/s  t: %d  d: %f s  iter: %d  fnodes: %f x1e6  geo: %s  kernel: %s  %s\n",
		perf, nThreads, duration, cd.MaxIterations, ld.nFluid / 1e6,
		geometryType, kernelToUse,
#ifdef VERIFICATION
		"# VERIFICATION"
#else
		"# benchmark"
#endif
	);

	int exitCode = 0;

#ifdef VERIFICATION

	if (verify) {
		printf("# VERIFICATION: deviation from analytical solution: %e\n", errorNorm);

		if (errorNorm > 0.1) {
			printf("# VERIFICATION FAILED.\n");
			exitCode = 1;
		}
		else {
			printf("# VERIFICATION SUCCEEDED.\n");
		}
	}
#else
//	printf("# VERIFICATION: deviation from analytical solution: %e\n", errorNorm);
//	printf("# VERIFICATION: this is only valid for pipe geometry with enough iterations performed.\n");
#endif

	MemFree((void **)&ld.Lattice);

	return exitCode;
}
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 <stdio.h>
	28	#include <stdlib.h>
	29	#include <string.h>
	30	#include <strings.h> // strcasecmp
	31
	32	#include <inttypes.h>
	33
	34	#ifdef _OPENMP
	35	#include <omp.h>
	36	#endif
	37
	38	#include "Base.h"
	39	#include "Kernel.h"
	40	#include "Memory.h"
	41
	42	#include "Lattice.h"
	43	#include "Geometry.h"
	44	#include "Pinning.h"
	45	#include "LikwidIf.h"
	46
	47	#include "KernelFunctions.h"
	48
	49	#ifdef __x86_64__
	50	#include <xmmintrin.h>
	51
	52
	53	#define MXCSR_DAZ 6
	54	#define MXCSR_FTZ 15
	55
	56
	57	int FpIsMxCsrMaskSet(unsigned int mask)
	58	{
	59	unsigned int mxcsr;
	60	unsigned int mxcsrNew;
	61
	62	mxcsr = _mm_getcsr();
	63
	64	mxcsrNew = mxcsr & mask;
65
66	return (mxcsrNew == mask);
67	}
68
69	int FpGetFtz()
70	{
71	return FpIsMxCsrMaskSet(1 << MXCSR_FTZ);
72	}
73
74	int FpGetDaz()
75	{
76	return FpIsMxCsrMaskSet(1 << MXCSR_DAZ);
77	}
78	#endif
79
80
81	int ParseDimensions(const char * parameter, int * nX, int * nY, int * nZ)
82	{
83	char * tmp;
84
85	*nX = atoi(parameter);
86
87	if (*nX <= 0) {
88	printf("ERROR: parameter for X dimension must be > 0.\n");
89	return 0;
90	}
91
92	tmp = strchr(parameter, 'x');
93
94	if (tmp == NULL) {
95	printf("ERROR: parameter for Y dimension is missing.\n");
96	return 0;
97	}
98
99	*nY = atoi(tmp + 1);
100
101	if (*nY <= 0) {
102	printf("ERROR: parameter for Y dimension must be > 0.\n");
103	return 0;
104	}
105
106	tmp = strchr(tmp + 1, 'x');
107
108	if (tmp == NULL) {
109	printf("ERROR: parameter for Z dimension is missing.\n");
110	return 0;
111	}
112
113	*nZ = atoi(tmp + 1);
114
115	if (*nZ <= 0) {
116	printf("ERROR: parameter for Z dimension must be > 0.\n");
117	return 0;
118	}
119
120	return 1;
121	}
122
123	int main(int argc, char * argv[])
124	{
125	int dims[3] = { 20, 20, 20 }; // Dimensions in x, y, and z direction
126	const char * geometryType = "channel";
127	// int latticeDumpAscii = 0;
128	int verify = 0; UNUSED(verify);
129	char * kernelToUse = NULL;
130	int nThreads = 1;
131	const char * pinString = NULL;
132	int periodic[3] = { 0 };
133
134	CaseData cd;
135
136	cd.MaxIterations = 1000;
137	cd.RhoIn = 1.0;
138	cd.RhoOut = 1.0;
139	cd.Omega = 1.0;
140	cd.VtkOutput = 0;
141	cd.VtkModulus = 100;
142	cd.StatisticsModulus = 100;
143	cd.XForce = 0.00001;
144	kernelToUse = "push-soa";
145
146	Parameters p;
147	p.nArgs = argc;
148	p.Args = argv;
149	p.nKernelArgs = 0;
150	p.KernelArgs = NULL;
151
152	#define LBM_BENCH_KERNELS_VERSION_MAJOR 0
153	#define LBM_BENCH_KERNELS_VERSION_MINOR 1
154
155	printf("Lattice Boltzmann Benchmark Kernels (LbmBenchKernels) Copyright (C) 2016, 2017 LSS, RRZE\n");
156	printf("This program comes with ABSOLUTELY NO WARRANTY; for details see LICENSE.\n");
157	printf("This is free software, and you are welcome to redistribute it under certain conditions.\n");
158	printf("\n");
159	printf("LBM Benchmark Kernels %d.%d, compiled %s %s, type: %s\n",
160	LBM_BENCH_KERNELS_VERSION_MAJOR, LBM_BENCH_KERNELS_VERSION_MINOR, __DATE__, __TIME__,
161	#ifdef VERIFICATION
162	"verification"
163	#else
164	"benchmark"
165	#endif
166	);
167
168	// ----------------------------------------------------------------------
169	// Parse command line arguments
170
171	#define ARG_IS(param) (!strcmp(argv[i], param))
172	#define NEXT_ARG_PRESENT() \
173	do { \
174	if (i + 1 >= argc) { \
175	printf("ERROR: argument %s requires a parameter.\n", argv[i]); \
176	return 1; \
177	} \
178	} while (0)
179
180	for (int i = 1; i < argc; ++i) {
181
182	if (ARG_IS("-dims") \|\| ARG_IS("--dims")) {
183	NEXT_ARG_PRESENT();
184
185
186	if (!ParseDimensions(argv[++i], &dims[0], &dims[1], &dims[2])) {
187	return 1;
188	}
189	}
190	// else if (ARG_IS("-lattice-dump-ascii") \|\| ARG_IS("--lattice-dump-ascii")) {
191	// latticeDumpAscii = 1;
192	// }
193	else if (ARG_IS("-geometry") \|\| ARG_IS("--geometry")) {
194	NEXT_ARG_PRESENT();
195
196	geometryType = argv[++i];
197	}
198	else if (ARG_IS("-iterations") \|\|ARG_IS("--iterations")) {
199	NEXT_ARG_PRESENT();
200
201	cd.MaxIterations = strtol(argv[++i], NULL, 0);
202
203	if (cd.MaxIterations <= 0) {
204	printf("ERROR: number of iterations must be > 0.\n");
205	return 1;
206	}
207	}
208	else if (ARG_IS("-rho-in") \|\|ARG_IS("--rho-in")) {
209	NEXT_ARG_PRESENT();
210
211	cd.RhoIn = strtod(argv[++i], NULL);
212	}
213	else if (ARG_IS("-rho-out") \|\|ARG_IS("--rho-out")) {
214	NEXT_ARG_PRESENT();
215
216	cd.RhoOut = strtod(argv[++i], NULL);
217	}
218	else if (ARG_IS("-omega") \|\|ARG_IS("--omega")) {
219	NEXT_ARG_PRESENT();
220
221	cd.Omega = strtod(argv[++i], NULL);
222	}
223	else if (ARG_IS("-x-force") \|\|ARG_IS("--x-force")) {
224	NEXT_ARG_PRESENT();
225
226	cd.XForce = strtod(argv[++i], NULL);
227	}
228	else if (ARG_IS("-verify") \|\| ARG_IS("--verify")) {
229	#ifdef VERIFICATION
230
231	// Choose this preset for verification. As geometry type "box" is
232	// used but x and y direction are made pridoc.
233	// Everything else can be altered, but enough iterations should be
234	// performed in order to receive a fully developed flow field.
235	verify = 1;
236
237	cd.Omega = 1.0;
238	cd.RhoIn = 1.0;
239	cd.RhoOut = 1.0;
240	geometryType = "box";
241	dims[0] = 16;
242	dims[1] = 16;
243	dims[2] = 16;
244	cd.XForce = 0.00001;
245	cd.MaxIterations = 1000;
246	periodic[0] = 1;
247	periodic[1] = 1;
248	periodic[2] = 0;
249
250	printf("#\n");
251	printf("# VERIFICATION: verifying flow profile of channel flow.\n");
252	printf("#\n");
253
254	// TODO: this is not a good idea as we ignore all other options...
255
256	#else
257	printf("ERROR: in order to use -verify VERIFICATION must be defined during compilation.\n");
258	printf(" Recompile with VERIFICATION=on.\n");
259	return 1;
260	#endif
261	}
262	else if (ARG_IS("-vtk") \|\| ARG_IS("--vtk")) {
263	#ifdef VTK_OUTPUT
264
265	cd.VtkOutput = 1;
266
267	// If the next parameter is a number it is used as the itartion count,
268	// if not it is probably another parameter.
269	if (i + 1 < argc) {
270
271	int vtkModulus = atoi(argv[i+1]);
272
273	if (vtkModulus > 0) {
274	cd.VtkModulus = vtkModulus;
275	++i;
276	}
277	}
278	#else
279	printf("ERROR: in order to use -vtk VTK_OUTPUT must be defined during compilation.\n");
280	printf(" Recompile with VTK_OUTPUT=on.\n");
281	return 1;
282	#endif
283	}
284	else if (ARG_IS("-statistics") \|\| ARG_IS("--statistics")) {
285	#ifdef STATISTICS
286	NEXT_ARG_PRESENT();
287
288	cd.StatisticsModulus = atoi(argv[++i]);
289
290	if (cd.StatisticsModulus <= 0) {
291	printf("ERROR: the iteration count for -statistics must be > 0.\n");
292	return 1;
293	}
294	#else
295	printf("ERROR: in order to use -statistics STATISTICS must be defined during compilation.\n");
296	printf(" Recompile with STATISTICS=on.\n");
297	return 1;
298	#endif
299	}
300	else if (ARG_IS("-kernel") \|\| ARG_IS("--kernel")) {
301	NEXT_ARG_PRESENT();
302
303	kernelToUse = argv[++i];
304	}
305	else if (ARG_IS("-list") \|\| ARG_IS("--list")) {
306	printf("Available kernels to benchmark:\n");
307
308	for (int j = 0; j < N_ELEMS(g_kernels); ++j) {
309	printf(" %s\n", g_kernels[j].Name);
310	}
311
312	return 0;
313	}
314	else if (ARG_IS("-pin") \|\| ARG_IS("--pin")) {
315	NEXT_ARG_PRESENT();
316
317	pinString = argv[++i];
318	}
319	else if (ARG_IS("-t") \|\| ARG_IS("-threads") \|\| ARG_IS("--threads")) {
320	#ifdef _OPENMP
321	NEXT_ARG_PRESENT();
322
323	nThreads = atoi(argv[++i]);
324
325	if (nThreads <= 0) {
326	printf("ERROR: number of threads must be > 0.\n");
327	return 1;
328	}
329	#else
330	printf("ERROR: specifying number of threads is only available when compiled with OpenMP support.\n");
331	return 1;
332	#endif
333	}
334	else if (ARG_IS("-periodic-x") \|\| ARG_IS("--periodic-x")) {
335	periodic[0] = 1;
336	}
337	else if (ARG_IS("-periodic-y") \|\| ARG_IS("--periodic-y")) {
338	periodic[1] = 1;
339	}
340	else if (ARG_IS("-periodic-z") \|\| ARG_IS("--periodic-z")) {
341	periodic[2] = 1;
342	}
343	else if (ARG_IS("-h") \|\| ARG_IS("-help") \|\| ARG_IS("--help")) {
344	printf("ERROR: unknown argument: %s\n", argv[i]);
345	printf("\n");
346	printf("Usage:\n");
347	printf("./lbmbenchk -list\n");
348	printf("./lbmbenchk \n");
349	printf(" [-dims XxYyZ] [-geometry box\|channel\|pipe\|porosity[-value]] [-iterations <iterations>] [-lattice-dump-ascii]\n");
350	printf(" [-rho-in <density>] [-rho-out <density] [-omega <omega>] [-kernel <kernel>]\n");
351	printf(" [-periodic-x]\n");
352	#ifdef STATISTICS
353	printf(" [-statistics <every-n-iteration>]\n");
354	#endif
355	#ifdef VTK_OUTPUT
356	printf(" [-vtk [<every-n-iteration>]]\n");
357	#endif
358	#ifdef _OPENMP
359	printf(" [-t <number of threads>]\n");
360	#endif
361	printf(" [-pin core{,core}*]\n");
362	#ifdef VERIFICATION
363	printf(" [-verify]\n");
364	#endif
365	printf(" -- <kernel specific parameters>\n");
366	printf("\n");
367	printf("-list List available kernels.\n");
368	printf("\n");
369	printf("-dims XxYxZ Specify geometry dimensions.\n");
370	printf("\n");
371	printf("-geometry porosity-<value>\n");
372	printf(" Geometetry with blocks of size <value> regularily layout out.\n");
373	printf("\n");
374	return 1;
375	}
376	else if (ARG_IS("--")) {
377	// printf("# kernel args start with %s these are %d args.\n", argv[i + 1], argc - i - 1);
378	p.KernelArgs = &argv[++i];
379	p.nKernelArgs = argc - i;
380	break;
381	}
382	else {
383	printf("ERROR: unknown parameter: %s.\n", argv[i]);
384	exit(1);
385	}
386	}
387
388	#undef ARG_IS
389	#undef NEXT_ARG_PRESENT
390
391
392	// ----------------------------------------------------------------------
393	// Check if we exceed our index addressing PDFs.
394
395	{
396	uint64_t nPdfs = ((uint64_t)19) * dims[0] * dims[1] * dims[2];
397
398	if (nPdfs > ((2LU << 31) - 1)) {
399	printf("ERROR: number of PDFs exceed 2^31.\n");
400	exit(1);
401	}
402	}
403
404	// ----------------------------------------------------------------------
405
406	#ifdef _OPENMP
407	omp_set_num_threads(nThreads);
408	#endif
409
410	LatticeDesc ld;
411
412	GeoCreateByStr(geometryType, dims, periodic, &ld);
413
414	const char * defines[] = {
415	#ifdef VTK_OUTPUT
416	"VTK_OUTPUT",
417	#endif
418	#ifdef STATISTICS
419	"STATISTICS",
420	#endif
421	#ifdef VERIFICATION
422	"VERIFICATION",
423	#endif
424	#ifdef _OPENMP
425	"_OPENMP",
426	#endif
427	#ifdef HAVE_LIKWID
428	"HAVE_LIKWID",
429	#endif
430	};
431
432	printf("# defines: ");
433	for (int j = 0; j < N_ELEMS(defines); ++j) {
434	printf("%s ", defines[j]);
435	}
436	printf("\n");
437
438	printf("# nodes total: % 10d\n", ld.nObst + ld.nFluid);
439	printf("# nodes fluid: % 10d (including inlet & outlet)\n", ld.nFluid);
440	printf("# nodes obstacles: % 10d\n", ld.nObst);
441	printf("# nodes inlet: % 10d\n", ld.nInlet);
442	printf("# nodes outlet: % 10d\n", ld.nOutlet);
443	printf("# periodicity: x: %d y: %d z: %d\n", ld.PeriodicX, ld.PeriodicY, ld.PeriodicZ);
444
445	#ifdef VTK_OUTPUT
446	printf("# VTK output: %d (every %d iteration)\n", cd.VtkOutput, cd.VtkModulus);
447	#endif
448	#ifdef STATISTICS
449	printf("# statistics: every %d iteration\n", cd.StatisticsModulus);
450	#endif
451
452	printf("# omega: %f\n", cd.Omega);
453	printf("# initial density at inlet/outlet:\n");
454	printf("# rho in: %e\n", cd.RhoIn);
455	printf("# rho out: %e\n", cd.RhoOut);
456	printf("# iterations: %d\n", cd.MaxIterations);
457
458	#ifdef __x86_64__
459	printf("# fp status: DAZ: %d FTZ: %d\n", FpGetDaz(), FpGetFtz());
460	#endif
461
462	#ifdef _OPENMP
463	printf("# OpenMP threads: %d\n", omp_get_max_threads());
464
465	if (pinString != NULL) {
466	#pragma omp parallel
467	{
468	int threadId = omp_get_thread_num();
469	int err;
470
471	err = PinCurrentThreadByCpuList(pinString, 0, 0, threadId);
472
473	if (err) {
474	printf("ERROR [thread %d]: pinning failed.\n", threadId);
475	exit(1);
476	}
477
478	const char * cpuList = PinCpuListAsString();
479	Assert(cpuList != NULL);
480
481	// Not so nice hack to print the thread ids ordered.
482	#pragma omp for ordered
483	for (int i = 0; i < omp_get_num_threads(); ++i) {
484	#pragma omp ordered
485	printf("# thread %2d pinned to core(s): %s\n", threadId, cpuList);
486	}
487
488	free((void *)cpuList);
489	}
490	}
491	#endif
492
493	KernelData * kd;
494
495	KernelFunctions * kf = NULL;
496
497	if (kernelToUse == NULL) {
498	kf = &g_kernels[0];
499	}
500	else {
501	for (int j = 0; j < N_ELEMS(g_kernels); ++j) {
502
503	if (!strcasecmp(kernelToUse, g_kernels[j].Name)) {
504	kf = &g_kernels[j];
505	break;
506	}
507	}
508	}
509
510	if (kf == NULL) {
511	printf("ERROR: requested kernel \"%s\" not found.\n", kernelToUse);
512	exit(1);
513	}
514
515	printf("#\n");
516	printf("# kernel: %s\n", kf->Name);
517	printf("#\n");
518
519	// Initialize kernel by calling its own initialization function
520	kf->Init(&ld, &kd, &p);
521
522	#ifdef VERIFICATION
523	if (verify) {
524	KernelSetInitialDensity( &ld, kd, &cd);
525	KernelSetInitialVelocity(&ld, kd, &cd);
526	}
527	#endif
528
529	printf("# starting kernel...\n");
530
531	X_LIKWID_INIT();
532
533	double timeStart = Time();
534
535	// Call the LBM kernel
536	kd->Kernel(&ld, kd, &cd);
537
538	double duration = Time() - timeStart;
539
540	X_LIKWID_DEINIT();
541
542	// Print some statistics...
543	KernelStatisticsAdv(kd, &ld, &cd, cd.MaxIterations, 1 /* force output */);
544
545	#ifdef VERIFICATION
546	PdfT errorNorm = -1.0;
547	KernelVerifiy(&ld, kd, &cd, &errorNorm);
548	#endif
549
550	// Deinitialize kernel by calling its own deinitialization function
551	kf->Deinit(&ld, &kd);
552
553
554	double perf = (double)ld.nFluid * (double)cd.MaxIterations / duration / 1.e6;
555
556	printf("P: %f MFLUP/s t: %d d: %f s iter: %d fnodes: %f x1e6 geo: %s kernel: %s %s\n",
557	perf, nThreads, duration, cd.MaxIterations, ld.nFluid / 1e6,
558	geometryType, kernelToUse,
559	#ifdef VERIFICATION
560	"# VERIFICATION"
561	#else
562	"# benchmark"
563	#endif
564	);
565
566	int exitCode = 0;
567
568	#ifdef VERIFICATION
569
570	if (verify) {
571	printf("# VERIFICATION: deviation from analytical solution: %e\n", errorNorm);
572
573	if (errorNorm > 0.1) {
574	printf("# VERIFICATION FAILED.\n");
575	exitCode = 1;
576	}
577	else {
578	printf("# VERIFICATION SUCCEEDED.\n");
579	}
580	}
581	#else
582	// printf("# VERIFICATION: deviation from analytical solution: %e\n", errorNorm);
583	// printf("# VERIFICATION: this is only valid for pipe geometry with enough iterations performed.\n");
584	#endif
585
586	MemFree((void **)&ld.Lattice);
587
588	return exitCode;
589	}