We develop an efficient open-source computational framework for the automated search for well placements in hot fracture-controlled reservoirs that sustainably optimize geothermal energy production. This search is performed via 3D simulations of groundwater flow and heat transfer. We model the reservoirs as geologically consistent, randomly generated discrete fracture networks (DFNs) 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 spatial discretization is based on the finite element method stabilized via the algebraic flux correction, treating the wells using the non-matching approach. For the time discretization, we use a semi-implicit approach to enhance the solver efficiency. The optimization is performed via a gradient-free global optimization algorithm. We will present the results of our optimization tests for randomly generated DFNs consisting of thousands of fractures, considering realistic values of physical parameters. Additionally, we will show how our framework can be utilized to analyze the structure of the above DFNs.