Rays (~1k) intersecting triangle meshes (~100k) with ca. 100 vertices each I am looking for a way to efficiently compute ray-triangle intersections for a large number of rays (~1k) and a large set of triangle-mesh bodies (~100k) with each composed of several hundred faces. My current approach towards optimization consists of 


thread parallelism for each ray
pre-select bodies by testing for ray-bounding box intersections
loop over faces of pre-selected bodies to detect ray-triangle intersections using the Möller-Trumbore algorithm.


Not impressed with performance so far. (Admittedly using python). Are there alternative solutions ?
Your procedure could be optimized in two directions:


If your rays are structured you could possibly take advantage of space coherence: if a ray intersects a given triangle, then a nearby rays has high possibility of intersecting the same triangle or an neighboring one. Therefore simple parallelism over the rays may not be optimal, you could at least process them in batches based on adjacency.
If you have adjacency information on the triangles of your meshes, you could possibly exploit it when testing a ray against the mesh: the Möller-Trumbore algorithm can inform you if the intersection point of the ray and the triangle plane lies inside or outside each of the edges, so you may have a hint at which (neighboring) triangle to check next. That said, your meshes are quite small so the extra overhead may not be worth it.