[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 {
		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" };

	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;
	}

	// 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) {
		periodic[0] = 1;
	}
	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 (blockSize == 0) {
			blockSize = 8;
		}

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