Integrate the Temperature module with other modeles.

[ArmstrongInCN]

Thank u for your wonderful code! How can the temperature module be integrated with other cases? For instance, I want to simulate the temperature changes during a collision process, but I found that simply enabling the temperature module causes the program to stop running after a few time steps. Why is this happening? Is there something wrong with my settings?

stl.zip

void main_setup() { // Collision test; required extensions in defines.hpp: FP16S, VOLUME_FORCE, FORCE_FIELD, TEMPERATURE, SURFACE, INTERACTIVE_GRAPHICS
	// ################################################################## define simulation box size, viscosity and volume force ###################################################################
	const uint3 lbm_N = resolution(float3(0.2f, 1.0f, 0.5f), 0.2 * 8000u); 
	const float si_nu = 0.148f;
	LBM lbm(lbm_N, units.nu(si_nu));
	// ###################################################################################### define geometry ######################################################################################
	Mesh* target = read_stl(get_exe_path() + "/stl/Submarine_body.stl"); // https://www.thingiverse.com/thing:814319/files
	Mesh* bullet = read_stl(get_exe_path() + "/stl/Submarine_rotor.stl"); // plane and rotor separated with Microsoft 3D Builder
	const float lbm_width = 0.95f * (float)lbm_N.y;
	const float scale = lbm_width / target->get_bounding_box_size().y; // scale plane and rotor to simulation box size
	target->scale(0.2 * scale);
	bullet->scale(0.2 * scale);
	const float3 offset = lbm.center() - target->get_bounding_box_center(); // move plane and rotor to simulation box center
	target->translate(offset);
	bullet->translate(offset);
	target->set_center(target->get_bounding_box_center()); // set center of meshes to their bounding box center
	bullet->set_center(bullet->get_bounding_box_center());
	lbm.voxelize_mesh_on_device(target);
	lbm.voxelize_mesh_on_device(bullet);
	const uint Nx = lbm.get_Nx(), Ny = lbm.get_Ny(), Nz = lbm.get_Nz(); parallel_for(lbm.get_N(), [&](ulong n) { uint x = 0u, y = 0u, z = 0u; lbm.coordinates(n, x, y, z);
	if (lbm.flags[n] == TYPE_S && y >= 0.5 * Ny) {
		lbm.flags[n] = TYPE_F;
		lbm.u.y[n] = -0.25f;
	}
	if (lbm.flags[n] == TYPE_S && y < 0.5 * Ny) {
		lbm.flags[n] = TYPE_F;
		lbm.u.y[n] = 0.025f;
	}
	if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) lbm.flags[n] = TYPE_S; // all non periodic
	});
	 // ####################################################################### run simulation, export images and data ##########################################################################
	lbm.graphics.visualization_modes = lbm.get_D()==1u ? VIS_PHI_RAYTRACE : VIS_PHI_RASTERIZE;
	lbm.run();
} /**/