The metal casting process, which is one of the key drivers of the manufacturing industry, involves several physical phenomena occurring simultaneously like fluid flow, phase change, and heat transfer which affect the casting yield and quality. Casting process modeling involves numerical modeling of these phenomena on a computer. In recent decades, this has become an inevitable tool for foundry engineers to make defect-free castings. To expedite computational time Graphics Processing Units (GPUs) are being increasingly used in the numerical modeling of heat transfer and fluid flow. Initially, in this work a CPU based implicit solver code is developed for solving the 3D unsteady energy equation including phase change numerically using Finite Volume Method (FVM) which predicts the thermal profile during solidification in the metal casting process in a completely filled mold. To address the computational bottleneck, which is identified as the linear algebraic solver based on the Bi-Conjugate Gradient Stabilized (BiCGSTAB) method, a GPU-based code is developed using CUDA toolkit and was implemented on the GPU. The CPU and GPU based codes are then validated against a commercial casting simulation code FLOW-3D CAST® for a simple casting part and against in-house experimental results for gravity die casting of a simple geometry. Parallel performance is analyzed for grid sizes ranging from 10x10x10 to 210x210x210 and for three time-step sizes. The performance of the GPU code based on occupancy and throughput is also investigated. The GPU code exhibits a maximum speedup of 308x compared to the CPU code for a grid size of 210x210x210 and a time-step size of 2s.