[LbmBenchmarkKernelsPublic.git] / src / Geometry.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 "Geometry.h"
#include "Memory.h"

#include <strings.h>
#include <math.h>

#include <errno.h>

void GeoCreateByStr(const char * geometryType, int dims[3], int periodic[3], LatticeDesc * ld)
{
	int type = -1;
	void * typeDetails = NULL;
	int tmp;

	if (strncasecmp("channel", geometryType, 7) == 0) {
		type = GEO_TYPE_CHANNEL;
	}
	else if (strncasecmp("box", geometryType, 3) == 0) {
		type = GEO_TYPE_BOX;
	}
	else if (strncasecmp("pipe", geometryType, 4) == 0) {
		type = GEO_TYPE_PIPE;
	}
	else if (strncasecmp("blocks", geometryType, 6) == 0) {
		type = GEO_TYPE_BLOCKS;

		// Default block size
		tmp = 8;

		if (strlen(geometryType) > 7) {
			int blockSize = atoi(&geometryType[7]);

			int dimMin = dims[0];

			if (dims[1] < dimMin) dimMin = dims[1];
			if (dims[2] < dimMin) dimMin = dims[2];

			if (blockSize < 0 || blockSize > dimMin / 2) {
				printf("ERROR: block size for geometry must be > 0 and smaller than half of the smalest dimension.\n");
				// TODO: find a better solution for handling errors in here.
				Verify(0);
			}

			tmp = blockSize;
		}

		typeDetails = &tmp;
	}
	else if (strncasecmp("fluid", geometryType, 5) == 0) {
		type = GEO_TYPE_FLUID;
	}
	else {
		printf("ERROR: unknown geometry specified.\n");
		Verify(0);
	}

	GeoCreateByType(type, typeDetails, dims, periodic, ld);

	return;
}

void GeoCreateByType(GEO_TYPES type, void * typeDetails, int dims[3], int periodic[3], LatticeDesc * ld)
{
	Assert(dims != NULL);
	Assert(dims[0] > 0);
	Assert(dims[1] > 0);
	Assert(dims[2] > 0);

	Assert(periodic != NULL);
	Assert(periodic[0] >= 0);
	Assert(periodic[1] >= 0);
	Assert(periodic[2] >= 0);

	Assert(ld != NULL);

	Assert(type >= GEO_TYPE_MIN);
	Assert(type <= GEO_TYPE_MAX);

	const char * geoTypeStr[] = { "box", "channel", "pipe", "blocks", "fluid" };

	printf("# geometry: %d x %d x %d nodes, type %d %s\n", dims[0], dims[1], dims[2], type, geoTypeStr[type]);

	ld->Dims[0] = dims[0];
	ld->Dims[1] = dims[1];
	ld->Dims[2] = dims[2];
	ld->nCells = dims[0] * dims[1] * dims[2];
	ld->PeriodicX = periodic[0];
	ld->PeriodicY = periodic[1];
	ld->PeriodicZ = periodic[2];

	LatticeT * lattice;
	MemAlloc((void **)&lattice, sizeof(LatticeT) * dims[0] * dims[1] * dims[2]);

	ld->Lattice = lattice;

	for (int z = 0; z < dims[2]; ++z) {
		for (int y = 0; y < dims[1]; ++y) {
			for (int x = 0; x < dims[0]; ++x) {
				lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_FLUID;
			}
		}
	}

	if (type == GEO_TYPE_CHANNEL || type == GEO_TYPE_BLOCKS || type == GEO_TYPE_PIPE) {
		periodic[0] = 1;
	}
	else if (type == GEO_TYPE_FLUID) {
		// Actually this geometry is non-periodic, but we need fluid at the boundary.
		periodic[0] = 1; periodic[1] = 1; periodic[2] = 1;
	}

	// Walls or periodic on first and last x plane.
	for (int z = 0; z < dims[2]; ++z) {
		for (int y = 0; y < dims[1]; ++y) {
			if(periodic[0]){
				lattice[L_INDEX_4(dims, 0, y, z)] 				= LAT_CELL_FLUID;
				lattice[L_INDEX_4(dims, dims[0] - 1, y, z)] 	= LAT_CELL_FLUID;
			} else {
				lattice[L_INDEX_4(dims, 0, y, z)] 				= LAT_CELL_OBSTACLE;
				lattice[L_INDEX_4(dims, dims[0] - 1, y, z)] 	= LAT_CELL_OBSTACLE;
			}
		}
	}

	// Walls or periodic on first and last y plane.
	for (int z = 0; z < dims[2]; ++z) {
		for (int x = 0; x < dims[0]; ++x) {
			if(periodic[1]){
				lattice[L_INDEX_4(dims, x, 0, z)] 				= LAT_CELL_FLUID;
				lattice[L_INDEX_4(dims, x, dims[1] - 1, z)] 	= LAT_CELL_FLUID;
			} else {
				lattice[L_INDEX_4(dims, x, 0, z)] 				= LAT_CELL_OBSTACLE;
				lattice[L_INDEX_4(dims, x, dims[1] - 1, z)] 	= LAT_CELL_OBSTACLE;
			}
		}
	}

	// Walls or periodic on first and last z plane.
	for (int y = 0; y < dims[1]; ++y) {
		for (int x = 0; x < dims[0]; ++x) {
			if(periodic[2]){
				lattice[L_INDEX_4(dims, x, y, 0)] 				= LAT_CELL_FLUID;
				lattice[L_INDEX_4(dims, x, y, dims[2] - 1)] 	= LAT_CELL_FLUID;
			} else {
				lattice[L_INDEX_4(dims, x, y, 0)] 				= LAT_CELL_OBSTACLE;
				lattice[L_INDEX_4(dims, x, y, dims[2] - 1)] 	= LAT_CELL_OBSTACLE;
			}
		}
	}

	if (type == GEO_TYPE_CHANNEL) {
		// Nothing todo here.
	}
	else if (type == GEO_TYPE_PIPE) {
		#define SQR(a) ((a)*(a))
		double centerZ = dims[2] / 2.0 - 0.5;
		double centerY = dims[1] / 2.0 - 0.5;
		double minDiameter = MIN(dims[1], dims[2]);
		double minRadiusSquared = SQR(minDiameter / 2 - 1);

		for (int z = 0; z < dims[2]; ++z) {
			for (int y = 0; y < dims[1]; ++y) {
				if((SQR(z - centerZ) + SQR(y - centerY)) >= minRadiusSquared) {
					for (int x = 0; x < dims[0]; ++x) {
						lattice[L_INDEX_4(dims, x, y, z)] 	= LAT_CELL_OBSTACLE;
					}
				}
			}
		}
		#undef SQR
	}
	else if (type == GEO_TYPE_BLOCKS) {

		int blockSize = *((int *)typeDetails);

#if 0
		if (blockSize == 0) {
			blockSize = 8;
		}
#endif
		if (blockSize > 0) {

			int dimMin = dims[0];

			if (dims[1] < dimMin) dimMin = dims[1];
			if (dims[2] < dimMin) dimMin = dims[2];

			if (blockSize < 0 || blockSize > dimMin / 2) {
				printf("ERROR: block size for geometry must be > 0 and smaller than half of the smalest dimension.\n");
				// TODO: find a better solution for handling errors in here.
				Verify(0);
			}

			// Number of blocks in x, y, and z direction.
			int nbx = blockSize, nby = blockSize, nbz = blockSize;

			for (int z = 0; z < dims[2]; ++z) {
					if ((z % (2 * nbz)) < nbz) continue;

				for (int y = 0; y < dims[1]; ++y) {
					if ((y % (2 * nby)) < nby) continue;

					for (int x = 0; x < dims[0]; ++x) {

						if ((x % (2 * nbx)) >= nbx) {
							lattice[L_INDEX_4(dims, x, y, z)] 	= LAT_CELL_OBSTACLE;
						}
					}
				}
			}
		}
	}

// 	if (latticeDumpAscii) {
// 		const char strLatCellType[] = "X.IxO"; // X = Obstacle, . = Fluid, I = inlet, O = outlet
// 		for (int z = dims[2] - 1; z >= 0; --z) {
// 			printf("plane % 2d\n", z);
//
// 			for (int y = dims[1] - 1; y >= 0; --y) {
// 				printf(" %2d  ", y);
// 				for (int x = 0; x < dims[0]; ++x) {
// 					printf("%c", strLatCellType[lattice[L_INDEX_4(dims, x, y, z)]]);
// 				}
// 				printf("\n");
// 			}
// 		}
//	}

// Lattice Helper Function

	ld->nObst = 0;
	ld->nFluid = 0;
	ld->nInlet = 0;
	ld->nOutlet = 0;

	for (int z = 0; z < dims[2]; ++z) {
		for (int y = 0; y < dims[1]; ++y) {
			for (int x = 0; x < dims[0]; ++x) {
				switch (lattice[L_INDEX_4(dims, x, y, z)]) {
					case LAT_CELL_OBSTACLE:					ld->nObst++; break;
					case LAT_CELL_FLUID:					ld->nFluid++; break;
					case LAT_CELL_INLET:					ld->nInlet++;  ld->nFluid++; break;
					case LAT_CELL_OUTLET:					ld->nOutlet++; ld->nFluid++; break;
					default:
						Verify(0);
				}
			}
		}
	}

	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 "Geometry.h"
	28	#include "Memory.h"
	29
	30	#include <strings.h>
	31	#include <math.h>
	32
	33	#include <errno.h>
	34
	35	void GeoCreateByStr(const char * geometryType, int dims[3], int periodic[3], LatticeDesc * ld)
	36	{
	37	int type = -1;
	38	void * typeDetails = NULL;
	39	int tmp;
	40
	41	if (strncasecmp("channel", geometryType, 7) == 0) {
	42	type = GEO_TYPE_CHANNEL;
	43	}
	44	else if (strncasecmp("box", geometryType, 3) == 0) {
	45	type = GEO_TYPE_BOX;
	46	}
	47	else if (strncasecmp("pipe", geometryType, 4) == 0) {
	48	type = GEO_TYPE_PIPE;
	49	}
	50	else if (strncasecmp("blocks", geometryType, 6) == 0) {
	51	type = GEO_TYPE_BLOCKS;
	52
	53	// Default block size
	54	tmp = 8;
	55
	56	if (strlen(geometryType) > 7) {
	57	int blockSize = atoi(&geometryType[7]);
	58
	59	int dimMin = dims[0];
	60
	61	if (dims[1] < dimMin) dimMin = dims[1];
	62	if (dims[2] < dimMin) dimMin = dims[2];
	63
	64	if (blockSize < 0 \|\| blockSize > dimMin / 2) {
65	printf("ERROR: block size for geometry must be > 0 and smaller than half of the smalest dimension.\n");
66	// TODO: find a better solution for handling errors in here.
67	Verify(0);
68	}
69
70	tmp = blockSize;
71	}
72
73	typeDetails = &tmp;
74	}
e3f82424 MW	75	else if (strncasecmp("fluid", geometryType, 5) == 0) {
	76	type = GEO_TYPE_FLUID;
	77	}
10988083 MW	78	else {
	79	printf("ERROR: unknown geometry specified.\n");
	80	Verify(0);
	81	}
	82
	83	GeoCreateByType(type, typeDetails, dims, periodic, ld);
	84
	85	return;
	86	}
	87
	88	void GeoCreateByType(GEO_TYPES type, void * typeDetails, int dims[3], int periodic[3], LatticeDesc * ld)
	89	{
	90	Assert(dims != NULL);
	91	Assert(dims[0] > 0);
	92	Assert(dims[1] > 0);
	93	Assert(dims[2] > 0);
	94
	95	Assert(periodic != NULL);
	96	Assert(periodic[0] >= 0);
	97	Assert(periodic[1] >= 0);
	98	Assert(periodic[2] >= 0);
	99
	100	Assert(ld != NULL);
	101
	102	Assert(type >= GEO_TYPE_MIN);
	103	Assert(type <= GEO_TYPE_MAX);
	104
e3f82424	105	const char * geoTypeStr[] = { "box", "channel", "pipe", "blocks", "fluid" };
10988083 MW	106
	107	printf("# geometry: %d x %d x %d nodes, type %d %s\n", dims[0], dims[1], dims[2], type, geoTypeStr[type]);
	108
	109	ld->Dims[0] = dims[0];
	110	ld->Dims[1] = dims[1];
	111	ld->Dims[2] = dims[2];
	112	ld->nCells = dims[0] * dims[1] * dims[2];
	113	ld->PeriodicX = periodic[0];
	114	ld->PeriodicY = periodic[1];
	115	ld->PeriodicZ = periodic[2];
	116
	117	LatticeT * lattice;
	118	MemAlloc((void *)&lattice, sizeof(LatticeT) dims[0] * dims[1] * dims[2]);
	119
	120	ld->Lattice = lattice;
	121
	122	for (int z = 0; z < dims[2]; ++z) {
	123	for (int y = 0; y < dims[1]; ++y) {
	124	for (int x = 0; x < dims[0]; ++x) {
	125	lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_FLUID;
	126	}
	127	}
	128	}
	129
	130	if (type == GEO_TYPE_CHANNEL \|\| type == GEO_TYPE_BLOCKS \|\| type == GEO_TYPE_PIPE) {
	131	periodic[0] = 1;
	132	}
e3f82424 MW	133	else if (type == GEO_TYPE_FLUID) {
	134	// Actually this geometry is non-periodic, but we need fluid at the boundary.
	135	periodic[0] = 1; periodic[1] = 1; periodic[2] = 1;
	136	}
10988083 MW	137
	138	// Walls or periodic on first and last x plane.
	139	for (int z = 0; z < dims[2]; ++z) {
	140	for (int y = 0; y < dims[1]; ++y) {
	141	if(periodic[0]){
	142	lattice[L_INDEX_4(dims, 0, y, z)] = LAT_CELL_FLUID;
	143	lattice[L_INDEX_4(dims, dims[0] - 1, y, z)] = LAT_CELL_FLUID;
	144	} else {
	145	lattice[L_INDEX_4(dims, 0, y, z)] = LAT_CELL_OBSTACLE;
	146	lattice[L_INDEX_4(dims, dims[0] - 1, y, z)] = LAT_CELL_OBSTACLE;
	147	}
	148	}
	149	}
	150
	151	// Walls or periodic on first and last y plane.
	152	for (int z = 0; z < dims[2]; ++z) {
	153	for (int x = 0; x < dims[0]; ++x) {
	154	if(periodic[1]){
	155	lattice[L_INDEX_4(dims, x, 0, z)] = LAT_CELL_FLUID;
	156	lattice[L_INDEX_4(dims, x, dims[1] - 1, z)] = LAT_CELL_FLUID;
	157	} else {
	158	lattice[L_INDEX_4(dims, x, 0, z)] = LAT_CELL_OBSTACLE;
	159	lattice[L_INDEX_4(dims, x, dims[1] - 1, z)] = LAT_CELL_OBSTACLE;
	160	}
	161	}
	162	}
	163
	164	// Walls or periodic on first and last z plane.
	165	for (int y = 0; y < dims[1]; ++y) {
	166	for (int x = 0; x < dims[0]; ++x) {
	167	if(periodic[2]){
	168	lattice[L_INDEX_4(dims, x, y, 0)] = LAT_CELL_FLUID;
	169	lattice[L_INDEX_4(dims, x, y, dims[2] - 1)] = LAT_CELL_FLUID;
	170	} else {
	171	lattice[L_INDEX_4(dims, x, y, 0)] = LAT_CELL_OBSTACLE;
	172	lattice[L_INDEX_4(dims, x, y, dims[2] - 1)] = LAT_CELL_OBSTACLE;
	173	}
	174	}
	175	}
	176
	177	if (type == GEO_TYPE_CHANNEL) {
e3f82424	178	// Nothing todo here.
10988083 MW	179	}
	180	else if (type == GEO_TYPE_PIPE) {
	181	#define SQR(a) ((a)*(a))
	182	double centerZ = dims[2] / 2.0 - 0.5;
	183	double centerY = dims[1] / 2.0 - 0.5;
	184	double minDiameter = MIN(dims[1], dims[2]);
	185	double minRadiusSquared = SQR(minDiameter / 2 - 1);
	186
	187	for (int z = 0; z < dims[2]; ++z) {
	188	for (int y = 0; y < dims[1]; ++y) {
	189	if((SQR(z - centerZ) + SQR(y - centerY)) >= minRadiusSquared) {
	190	for (int x = 0; x < dims[0]; ++x) {
	191	lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_OBSTACLE;
	192	}
	193	}
	194	}
	195	}
	196	#undef SQR
	197	}
	198	else if (type == GEO_TYPE_BLOCKS) {
	199
	200	int blockSize = ((int )typeDetails);
	201
e3f82424	202	#if 0
10988083 MW	203	if (blockSize == 0) {
	204	blockSize = 8;
	205	}
e3f82424 MW	206	#endif
e3f82424 MW	207	if (blockSize > 0) {
10988083	208
e3f82424	209	int dimMin = dims[0];
10988083	210
e3f82424 MW	211	if (dims[1] < dimMin) dimMin = dims[1];
e3f82424 MW	212	if (dims[2] < dimMin) dimMin = dims[2];
10988083	213
e3f82424 MW	214	if (blockSize < 0 \|\| blockSize > dimMin / 2) {
	215	printf("ERROR: block size for geometry must be > 0 and smaller than half of the smalest dimension.\n");
	216	// TODO: find a better solution for handling errors in here.
	217	Verify(0);
	218	}
10988083	219
e3f82424 MW	220	// Number of blocks in x, y, and z direction.
e3f82424 MW	221	int nbx = blockSize, nby = blockSize, nbz = blockSize;
10988083	222
e3f82424 MW	223	for (int z = 0; z < dims[2]; ++z) {
e3f82424 MW	224	if ((z % (2 * nbz)) < nbz) continue;
10988083	225
e3f82424 MW	226	for (int y = 0; y < dims[1]; ++y) {
e3f82424 MW	227	if ((y % (2 * nby)) < nby) continue;
10988083	228
e3f82424	229	for (int x = 0; x < dims[0]; ++x) {
10988083	230
e3f82424 MW	231	if ((x % (2 * nbx)) >= nbx) {
	232	lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_OBSTACLE;
	233	}
10988083 MW	234	}
	235	}
	236	}
	237	}
	238	}
	239
	240	// if (latticeDumpAscii) {
	241	// const char strLatCellType[] = "X.IxO"; // X = Obstacle, . = Fluid, I = inlet, O = outlet
	242	// for (int z = dims[2] - 1; z >= 0; --z) {
	243	// printf("plane % 2d\n", z);
	244	//
	245	// for (int y = dims[1] - 1; y >= 0; --y) {
	246	// printf(" %2d ", y);
	247	// for (int x = 0; x < dims[0]; ++x) {
	248	// printf("%c", strLatCellType[lattice[L_INDEX_4(dims, x, y, z)]]);
	249	// }
	250	// printf("\n");
	251	// }
	252	// }
	253	// }
	254
	255	// Lattice Helper Function
	256
	257	ld->nObst = 0;
	258	ld->nFluid = 0;
	259	ld->nInlet = 0;
	260	ld->nOutlet = 0;
	261
	262	for (int z = 0; z < dims[2]; ++z) {
	263	for (int y = 0; y < dims[1]; ++y) {
	264	for (int x = 0; x < dims[0]; ++x) {
	265	switch (lattice[L_INDEX_4(dims, x, y, z)]) {
	266	case LAT_CELL_OBSTACLE: ld->nObst++; break;
	267	case LAT_CELL_FLUID: ld->nFluid++; break;
	268	case LAT_CELL_INLET: ld->nInlet++; ld->nFluid++; break;
	269	case LAT_CELL_OUTLET: ld->nOutlet++; ld->nFluid++; break;
	270	default:
	271	Verify(0);
	272	}
	273	}
	274	}
	275	}
	276
	277	return;
	278	}