We develop an efficient open-source numerical framework for the automated search for the placements of injection and production wells in hot fracture-controlled reservoirs that sustainably optimize geothermal energy production. This search is carried out via 3D simulations of groundwater flow and heat transfer. We model the reservoirs as geologically consistent randomly generated discrete fracture networks (DFN) in which the fractures are 2D manifolds with polygonal boundaries embedded in a 3D porous medium. The wells are modeled via the immersed boundary method. The flow and heat transport in the DFN-matrix system are modeled by solving the balance equations for mass, momentum, and energy. The fully developed numerical framework combines the finite element method with semi-implicit time-stepping, algebraic flux correction, and approximation of the wells via the non-matching approach. To perform the optimization, we use several gradient-free algorithms. We will present our latest results, considering geologically and physically realistic scenarios.