High-fidelity 3D electromagnetic (EM) simulations play a key role in subsurface exploration and Earth system studies, yet they remain computationally demanding due to the use of large meshes, heterogeneous resistivity models, and dense source-receiver configurations. In this work, we present recent developments in an upgraded version of PETGEM (Parallel Exascale Toolkit for Geophysical Electromagnetic Modeling), building on previous contributions and focusing on practical challenges encountered when targeting large-scale HPC systems. A major effort has been the transition from a Python-based implementation to a fully C-based code using PETSc, enabling improved control over data structures, solvers, and parallel execution. We report on our experience with large-scale forward modeling for controlled-source EM applications, including performance analysis using Extrae to identify computational bottlenecks and guide targeted optimizations. These efforts have resulted in improved scalability and robustness on modern supercomputing platforms. Beyond performance, we emphasize the importance of software sustainability for long-term scientific development. We discuss the integration of best programming practices such as modular code design, automated testing, and continuous integration workflows, which are essential as PETGEM evolves toward inverse modeling capabilities. Overall, this contribution shares lessons learned from advancing an established EM modeling code toward extreme-scale computing, highlighting the interplay between numerical methods, performance engineering, and sustainable software development in computational geosciences.