The demand for surrogate models in Computer-Aided Engineering (CAE) tools has grown rapidly in recent years. Among many methods for surrogate modeling, DeepONet and its derivatives are one of the most widely used frameworks. DeepONet can be physics-informed (PI-DeepONet) by introducing physics-informed loss, which has been innovated in Physics-informed neural networks (PINNs), and it has been shown to improve the generalization of the trained model. During the training of PI-DeepONet or data-hybrid PINN, a decrease in the total loss function often stops at an unacceptable level. One of the reasons is due to the truncation error involved in the CAE-generated data. The partial differential equation (PDE), which is supposed to be solved in CAEs, is discretized into a linear or nonlinear algebraic equation and solved with limited spatial and/or time resolution. Thus, the equation actually solved is not the same as originally supposed to be solved, and the obtained solution is also different from the solution of the original PDE. In the surrogate modeling using physics losses, the data loss is evaluated using truncated data, whereas the physics loss is evaluated based on the original PDE. Therefore, the data and physics losses cannot be completely satisfied. To improve this problem, we propose a structure-preserving methodology to compensate for the fidelity of the data in physics-informed surrogate modeling. We assumed that a discretized PDE can be written as F*(u* + u') = F(u*) + G(u*), where F and F* are the original and discretized PDEs, respectively. u* and u' are the truncated solution and compensation. G is an unknown term corresponding to the difference between the original and the discretized equations, and is designed to be discovered from the data. In the machine learning framework, two neural networks are introduced to predict u* and u', and trained to minimize data loss and two physics losses defined by evaluating F(u* + u') = 0 and F(u*) + G(u*) = 0.