Convergent Operator-Splitting Scheme for Viscosity Solutions: A Foundation for Learning Domain-to-Solution Maps

This work introduces and rigorously analyzes a novel operator-splitting finite element scheme for approximating viscosity solutions of a broad class of constrained second-order partial differential equations. By decoupling the primary PDE evolution from the enforcement of constraints, the proposed method combines a stabilized finite element method for spatial discretization with an efficient semi-implicit time-stepping strategy. The cornerstone of our analysis is a proof that the scheme satisfies a discrete comparison principle. We demonstrate that under a mild time-step restriction and with appropriate stabilization, the discrete operator yields an M-matrix, which is sufficient to guarantee the scheme's monotonicity and consequent $L^\infty$-stability. These properties -- consistency, stability, and monotonicity -- are shown to be sufficient to prove convergence of the numerical approximation to the unique viscosity solution within the celebrated Barles--Souganidis framework. For solutions with enhanced regularity, we further establish an optimal-order error estimate of $O(\Delta t + h^2)$. The rigorously established stability of the scheme provides a blueprint for a novel Physics-Constrained Neural Operator (PCNO) architecture. We prove that by emulating the scheme's structure, the PCNO can provably break the curse of dimensionality for the challenging class of domain-to-solution mapping problems with complex topological variations, a problem for which standard learning approaches often fail. Numerical experiments for both a Hamilton-Jacobi equation with state constraints and a controlled reaction-diffusion system are presented to validate the theoretical findings and demonstrate the scheme's effectiveness.