In order to simulate the developed non-linear mathematical model and optimization problem, the continuous Galerkin finite element discretization was employed to solve the governing equations, thereby transforming them into discrete forms suitable for computational analysis. The weak form was derived utilizing a standard first-order finite element scheme for charge transport equations (Eqs. 28 and 29), while a second-order scheme was employed for the Navier-Stokes equations (Eq. 31). The computational mesh necessary for discretization was generated using a collection of first-order tetrahedral elements using Gmsh, ensuring compatibility with our chosen finite element method. The open-source PDE solver Firedrake was utilized for solving the resulting discrete equations, in which Newton's method was employed to linearize the system. For this purpose, part of the TOPE code was adapted and further developed to include the additional physics and requirements. In order to solve the optimization problem (Eq. 45), the method of moving asymptotes (MMA) algorithm implemented in pyMMAopt library was used, and the sensitivities were calculated using the pyadjoint library (dolphin-adjoint) by algorithmic differentiation .