Yes, you should get advantage on getting initial data with grasshopper, but the actual heat diffusion will be done properly/easy only through scripting, it needs iterations, many iterations. (I suggest c# against python )
Also i suspect a sphere with “aliased” surface (due to voxels) will behave wrongly, non uniform.
Total irradiating (naked) area of a voxel sphere is A=6\pi r^2 while a normal sphere is A=4\pi r^2 !
So, unless you calculate the self-irradiation on the “corrugated” surface of your voxel earth, i’d suggest you on using a thin layer/skin of out-of-grid voxels, like the face centers of “Mesh Sphere Ex”, to have a flat and smooth outer surface that receive and irradiate heat.
radiating outwards - count number of voxels in outer layer divide true area of a sphere by count, each voxel made radiate by that amount. It feels like it’s just burying the error, as there would proportionally be more voxels per km walk along the surface in areas where the sides of the voxel nearly line up with the surface compared to locations where it’s at 45 degrees and the voxels are spaced out by their diagonals not their sides.
So that would be another factor to account for. multiply up by the averaged angle between surface and voxel
Maybe something like this:
For each outer voxel, for each of its naked face, compute the angle between that face normal and local sphere normal, and multiply the area of the face by the cosine of that angle; then, sum the areas found that way for every naked voxel. And you have a real irradiating area for each voxel!
Which matrix will you use for the laplacian smoothing?