Passive heat sinks cooled by natural convection are reliable, compact, and low-noise. They are widely used in telecommunication devices, LEDs, and so forth. This work builds upon the recent advancements in fluid topology optimization (TO) to present a case study of two-and three-dimensional optimum design and thermal modeling for the natural convection problems using a reaction-diffusion equation (RDE)-based level-set method. To this end, first, a high-fidelity thermal-fluid model is constructed where the full Navier-Stokes equations are strongly coupled with the energy equation through the Boussinesq approximation. We benchmark our simulation solver against the experimental analysis and other numerical analysis methods. Next, we carefully investigate the flow behavior under different Grashof numbers using a fully transient simulation solver. Then, we revisit the RDE-based level-set TO methodology and construct an open-source parallel TO framework. The main findings reveal that using the body-fitted mesh adaptation, the proposed methodology can capture the explicit fluid-solid boundary and we are free of the continuation approach to penalize the design variable to the binary structure. A moderately large-scale TO problem with 3.56 × 10 6 DOFs can be solved in parallel on a standard multi-process platform. Our numerical implementation uses FreeFEM for finite element analysis, PETSc for distributed linear algebra, and Mmg for mesh adaptation. For comparison and for accessing our various techniques, a variety of 2D and 3D benchmarks are presented to support these remarkable features.