Machine Learning (ML) methods for solving Partial Differential Equations (PDEs) have recently undergone unprecedented development. Physics-Informed Neural Networks (PINNs) have gained significant attention for their ability to integrate underlying physics into learning frameworks \cite{Raissi2019}. However, classic PINNs rely on Automatic Differentiation (AD) to compute physics-based loss terms. During backpropagation, AD generates large tensors that increase training overhead and slow the optimization process \cite{Li2024a}. Furthermore, PINNs based on Deep Neural Networks (DNNs) often lack geometric adaptability. This limitation makes them suitable for structured regular grids rather than complex geometries. The challenge of reducing the computational overhead due to AD has been tackled using external solvers to compute the discrete residual associated to the target equation and using such a residual in a physics informed loss \cite{Halder2025a, Halder2025b}. Graph Neural Networks (GNNs) have been proved to solve the limitations of DNN based approaches, performing well in contexts where geometric adaptability is the key \cite{Li2024a, Li2024b}. We present a novel approach that couples a Graph Neural Network framework with an external numerical solver. This model pairs the superior geometric adaptability of GNNs with the reduced training overhead of external solvers. By using an external solver, our framework eliminates the complications represented by AD in graph aggregation and message passing contexts. The implementation bridges Python-based GNN architectures with a high-performance C++ solver. This talk covers our model's architecture, mathematical properties, examples, applications, and use cases with a particular focus on Computational Fluid Dynamics (CFD).