[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>

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

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];
	ld->Name = g_geoTypeStr[type];

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