Added functions to compute force/torque on objects, added more benchm…

…arks in Readme, added function to translate Mesh, added Stokes drag validation setup

README.md

@@ -138,6 +138,7 @@ Q = &Sigma;<sub><i>&alpha;&beta;</i></sub> <i>&Pi;<sub>&alpha;&beta;</sub></i><s
    - <kbd>V</kbd>: toggle stereoscopic rendering for VR
    - <kbd>B</kbd>: toggle VR-goggles/3D-TV mode for stereoscopic rendering
    - <kbd>N</kbd>/<kbd>M</kbd>: adjust eye distance for stereoscopic rendering
+   - <kbd>H</kbd>: print current camera position/rotation in console as copy/paste command
    - <kbd>Esc</kbd>/<kbd>Alt</kbd>+<kbd>F4</kbd>: quit
 
 
@@ -181,13 +182,19 @@ In consequence, the arithmetic intensity of this implementation is 2.13 (FP32/FP
 | Nvidia Tesla K40m             |               4.29 |          12 |          288 |             1131 (60%) |              1868 (50%) |               912 (24%) |
 | Nvidia Tesla K80  (1 GPU)     |               4.11 |          12 |          240 |              916 (58%) |              1642 (53%) |               943 (30%) |
 | Nvidia Tesla K20c             |               3.52 |           5 |          208 |              861 (63%) |              1507 (56%) |               720 (27%) |
+| AMD Radeon RX 6900 XT         |              20.63 |          16 |          512 |             1968 (59%) |              4227 (64%) |              4207 (63%) |
+| AMD Radeon RX 5700 XT         |               9.75 |           8 |          448 |             1368 (47%) |              3253 (56%) |              3049 (52%) |
 | AMD Radeon RX Vega 64         |              13.35 |           8 |          484 |             1875 (59%) |              2878 (46%) |              3227 (51%) |
+| AMD Radeon RX 580 4GB         |               6.50 |           4 |          256 |              946 (57%) |              1848 (56%) |              1577 (47%) |
 | Nvidia GeForce RTX 3090       |              39.05 |          24 |          936 |             5418 (89%) |             10732 (88%) |             10215 (84%) |
 | Nvidia GeForce RTX 3080 10GB  |              29.77 |          10 |          760 |             4230 (85%) |              8118 (82%) |              7714 (78%) |
+| Nvidia GeForce RTX 3060       |              13.17 |          12 |          360 |             2108 (90%) |              4070 (87%) |              3566 (76%) |
+| Nvidia GeForce RTX 3060M      |              10.94 |           6 |          336 |             2019 (92%) |              4012 (92%) |              3572 (82%) |
 | Nvidia Quadro RTX 8000 Server |              14.93 |          48 |          624 |             2591 (64%) |              5408 (67%) |              5607 (69%) |
 | Nvidia GeForce RTX 2080 Ti    |              13.45 |          11 |          616 |             3194 (79%) |              6700 (84%) |              6853 (86%) |
+| Nvidia GeForce RTX 2080 Super |              11.34 |           8 |          496 |             2434 (75%) |              5284 (82%) |              5087 (79%) |
 | Nvidia GeForce RTX 2060 Super |               7.18 |           8 |          448 |             2503 (85%) |              5035 (87%) |              4463 (77%) |
-| Nvidia Tesla T4               |               8.14 |          16 |          300 |             1356 (69%) |              2869 (74%) |              2887 (74%) |
+| Nvidia Tesla T4               |               8.14 |          15 |          300 |             1356 (69%) |              2869 (74%) |              2887 (74%) |
 | Nvidia GeForce GTX 1660 Ti    |               5.48 |           6 |          288 |             1467 (78%) |              3041 (81%) |              3019 (81%) |
 | Nvidia Titan Xp               |              12.15 |          12 |          548 |             2919 (82%) |              5495 (77%) |              5375 (76%) |
 | Nvidia GeForce GTX 1080       |               9.78 |           8 |          320 |             1623 (78%) |              3100 (75%) |              3182 (77%) |
@@ -200,7 +207,8 @@ In consequence, the arithmetic intensity of this implementation is 2.13 (FP32/FP
 | Nvidia Quadro K2000           |               0.73 |           2 |           64 |              312 (75%) |               444 (53%) |               171 (21%) |
 | Nvidia GeForce GT 630 (OEM)   |               0.46 |           2 |           29 |              151 (81%) |               185 (50%) |                78 (21%) |
 | Apple M1 Pro GPU 16C 16GB     |               4.10 |          11 |          200 |             1204 (92%) |              2329 (90%) |              1855 (71%) |
-| AMD Radeon Vega 8 Graphics    |               1.23 |           7 |           38 |              157 (63%) |               282 (57%) |               288 (58%) |
+| AMD Radeon Vega 8 (4750G)     |               2.15 |          27 |           57 |              263 (71%) |               511 (70%) |               501 (68%) |
+| AMD Radeon Vega 8 (3500U)     |               1.23 |           7 |           38 |              157 (63%) |               282 (57%) |               288 (58%) |
 | Intel UHD Graphics 630        |               0.46 |           7 |           51 |              151 (45%) |               301 (45%) |               187 (28%) |
 | Intel HD Graphics 5500        |               0.35 |           3 |           26 |               75 (45%) |               192 (58%) |               108 (32%) |
 | Intel HD Graphics 4600        |               0.38 |           2 |           26 |              105 (63%) |               115 (35%) |                34 (10%) |

src/kernel.cpp

@@ -1958,7 +1958,7 @@ string opencl_c_container() { return R( // ########################## begin of O
 	const float3 p5 = (float3)(p.x+0.5f, p.y-0.5f, p.z+0.5f); // +-+
 	const float3 p6 = (float3)(p.x+0.5f, p.y-0.5f, p.z-0.5f); // +--
 	const float3 p7 = (float3)(p.x-0.5f, p.y+0.5f, p.z+0.5f); // -++
-	if(!(not_xm||not_ym)) draw_line(p0, p2, c, camera_cache, bitmap, zbuffer); // to draw the entire surface, replace && by ||
+	if(!(not_xm||not_ym)) draw_line(p0, p2, c, camera_cache, bitmap, zbuffer); // to draw the entire surface, replace || by &&
 	if(!(not_xm||not_zm)) draw_line(p0, p4, c, camera_cache, bitmap, zbuffer);
 	if(!(not_ym||not_zm)) draw_line(p0, p6, c, camera_cache, bitmap, zbuffer);
 	if(!(not_xp||not_yp)) draw_line(p1, p3, c, camera_cache, bitmap, zbuffer);

src/lbm.cpp

@@ -81,7 +81,9 @@ void LBM::sanity_checks_constructor() { // sanity checks on grid resolution and
 #ifndef VOLUME_FORCE
 	if(fx!=0.0f||fy!=0.0f||fz!=0.0f) print_error("Volume force is set in LBM constructor in main_setup(), but VOLUME_FORCE is not enabled. Uncomment \"#define VOLUME_FORCE\" in defines.hpp.");
 #else // VOLUME_FORCE
+#ifndef FORCE_FIELD
 	if(fx==0.0f&&fy==0.0f&&fz==0.0f) print_warning("The VOLUME_FORCE extension is enabled but the volume force in LBM constructor is set to zero. You may disable the extension by commenting out \"#define VOLUME_FORCE\" in defines.hpp.");
+#endif // FORCE_FIELD
 #endif // VOLUME_FORCE
 #ifndef SURFACE
 	if(sigma!=0.0f) print_error("Surface tension is set in LBM constructor in main_setup(), but SURFACE is not enabled. Uncomment \"#define SURFACE\" in defines.hpp.");
@@ -243,6 +245,44 @@ void LBM::reset() { // reset simulation (takes effect in following run() call)
 void LBM::calculate_force_on_boundaries() { // calculate forces from fluid on TYPE_S nodes
 	kernel_calculate_force_on_boundaries.set_parameters(2u, t).run();
 }
+float3 LBM::calculate_force_on_object(const uchar flag_marker) { // add up force for all nodes flagged with flag_marker
+	double3 force(0.0, 0.0, 0.0);
+	for(uint n=0u; n<get_N(); n++) {
+		if(flags[n]==flag_marker) {
+			force.x += (double)F.x[n];
+			force.y += (double)F.y[n];
+			force.z += (double)F.z[n];
+		}
+	}
+	return float3(force.x, force.y, force.z);
+}
+float3 LBM::calculate_torque_on_object(const uchar flag_marker) { // add up torque around center of mass for all nodes flagged with flag_marker
+	double3 center_of_mass(0.0, 0.0, 0.0);
+	uint counter = 0u;
+	for(uint n=0u; n<get_N(); n++) {
+		if(flags[n]==flag_marker) {
+			const float3 p = position(n);
+			center_of_mass.x += (double)p.x;
+			center_of_mass.y += (double)p.y;
+			center_of_mass.z += (double)p.z;
+			counter++;
+		}
+	}
+	return calculate_torque_on_object(float3(center_of_mass.x/(double)counter, center_of_mass.y/(double)counter, center_of_mass.z/(double)counter)+center(), flag_marker);
+}
+float3 LBM::calculate_torque_on_object(const float3& rotation_center, const uchar flag_marker) { // add up torque around specified rotation center for all nodes flagged with flag_marker
+	double3 torque(0.0, 0.0, 0.0);
+	const float3 rotation_center_in_box = rotation_center-center();
+	for(uint n=0u; n<get_N(); n++) {
+		if(flags[n]==flag_marker) {
+			const float3 t = cross(position(n)-rotation_center_in_box, float3(F.x[n], F.y[n], F.z[n]));
+			torque.x += (double)t.x;
+			torque.y += (double)t.y;
+			torque.z += (double)t.z;
+		}
+	}
+	return float3(torque.x, torque.y, torque.z);
+}
 #endif // FORCE_FIELD
 
 #ifdef MOVING_BOUNDARIES
@@ -417,7 +457,7 @@ string LBM::device_defines() const { return
 	"\n	#define def_Nx "+to_string(Nx)+"u"
 	"\n	#define def_Ny "+to_string(Ny)+"u"
 	"\n	#define def_Nz "+to_string(Nz)+"u"
-	"\n	#define def_N " +to_string(get_N())+"ul"
+	"\n	#define def_N "+to_string(get_N())+"ul"
 
 	"\n	#define D"+to_string(dimensions)+"Q"+to_string(velocity_set)+"" // D2Q9/D3Q15/D3Q19/D3Q27
 	"\n	#define def_velocity_set "+to_string(velocity_set)+"u" // LBM velocity set (D2Q9/D3Q15/D3Q19/D3Q27)

src/lbm.hpp

@@ -151,6 +151,9 @@ class LBM {
 
 #ifdef FORCE_FIELD
 	void calculate_force_on_boundaries(); // calculate forces from fluid on TYPE_S nodes
+	float3 calculate_force_on_object(const uchar flag_marker=TYPE_S); // add up force for all nodes flagged with flag_marker
+	float3 calculate_torque_on_object(const uchar flag_marker=TYPE_S); // add up torque around center of mass for all nodes flagged with flag_marker
+	float3 calculate_torque_on_object(const float3& rotation_center, const uchar flag_marker=TYPE_S); // add up torque around specified rotation_center for all nodes flagged with flag_marker
 #endif // FORCE_FIELD
 
 #ifdef MOVING_BOUNDARIES

src/setup.cpp

@@ -114,6 +114,52 @@
 
 
 
+/*void main_setup() { // Stokes drag validation
+	const uint T = 100; // check error every T steps
+	const float R = 32.0f; // sphere radius
+	const float Re = 0.01f; // Reynolds number
+	const float nu = 1.0f; // kinematic shear viscosity
+	const float rho = 1.0f; // density
+	const uint L = to_uint(8.0f*R); // simulation box size
+	const float u = units.u_from_Re(Re, 2.0f*R, nu); // velocity
+	// ######################################################### define simulation box size, viscosity and volume force ############################################################################
+	LBM lbm(L, L, L, nu); // flow driven by equilibrium boundaries
+	// #############################################################################################################################################################################################
+	const uint N=lbm.get_N(), Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(uint n=0u, x=0u, y=0u, z=0u; n<N; n++, lbm.coordinates(n, x, y, z)) {
+		// ########################################################################### define geometry #############################################################################################
+		if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_E;
+		if(sphere(x, y, z, lbm.center(), R)) {
+			lbm.flags[n] = TYPE_S|TYPE_X; // flag boundary nodes for force summation additionally with TYPE_X
+		} else {
+			lbm.rho[n] = units.rho_Stokes(lbm.position(x, y, z), float3(-u, 0.0f, 0.0f), R, rho, nu);
+			const float3 un = units.u_Stokes(lbm.position(x, y, z), float3(-u, 0.0f, 0.0f), R);
+			lbm.u.x[n] = un.x;
+			lbm.u.y[n] = un.y;
+			lbm.u.z[n] = un.z;
+		}
+	}	// #########################################################################################################################################################################################
+	double E1=1000.0, E2=1000.0;
+	while(true) {
+		lbm.run(T);
+		lbm.calculate_force_on_boundaries();
+		lbm.F.read_from_device();
+		const float3 force = lbm.calculate_force_on_object(TYPE_S|TYPE_X);
+		const double F_theo = units.F_Stokes(rho, u, nu, R);
+		const double F_sim = (double)length(force);
+		const double E0 = fabs(F_sim-F_theo)/F_theo;
+		print_info(to_string(lbm.get_t())+", expected: "+to_string(F_theo, 6u)+", measured: "+to_string(F_sim, 6u)+", error = "+to_string((float)(100.0*E0), 1u)+"%");
+		if(converged(E2, E1, E0, 1E-4)) { // stop when error has sufficiently converged
+			print_info("Error converged after "+to_string(lbm.get_t())+" steps to "+to_string(100.0*E0, 1u)+"%");
+			wait();
+			break;
+		}
+		E2 = E1;
+		E1 = E0;
+	}
+} /**/
+
+
+
 /*void main_setup() { // cylinder in rectangular duct
 	// ######################################################### define simulation box size, viscosity and volume force ############################################################################
 	const float Re = 25000.0f;

src/utilities.hpp

@@ -4182,6 +4182,16 @@ struct Mesh { // triangle mesh
 		}
 		find_bounds();
 	}
+	inline void translate(const float3& translation) {
+		for(uint i=0u; i<triangle_number; i++) {
+			p0[i] += translation;
+			p1[i] += translation;
+			p2[i] += translation;
+		}
+		center += translation;
+		pmin += translation;
+		pmax += translation;
+	}
 	inline void rotate(const float3x3& rotation) {
 		for(uint i=0u; i<triangle_number; i++) {
 			p0[i] = rotation*(p0[i]-center)+center;
@@ -4190,6 +4200,9 @@ struct Mesh { // triangle mesh
 		}
 		find_bounds();
 	}
+	inline const float3& get_center() const {
+		return center;
+	}
 };
 inline Mesh* read_stl(const string& path, const float3& box_size, const float3& center, const float3x3& rotation, const float size) { // read binary .stl file
 	const string filename = create_file_extension(path, ".stl");