## 3 Experimental Verification
### One-Dimensional Parametric Euler Equation  
The conservation of mass, momentum, and energy in a compressible, inviscid flow can be described by the Euler equation. The conservative form of the Euler equation can be written as follows [32]:  
\[\partial_{t}U+\
abla\cdot f(U)=0,x\in\Omega\subset R^{d},d=1,2,t\in(0,T]. \tag{19}\]  
The one-dimensional Euler equation has the following form:  
\[\text{U}=\left( \begin{array}{c}\rho\\ \rho u\\ \rho E\end{array} \right),f(U)=\left( \begin{array}{c}\rho u\\ \rho u^{2}+p\\ u(\rho E+p)\end{array} \right), \tag{20}\]  
Figure 10: Input sample of \(u(x)\) for Figure 8.  
where \(\rho\) is the density, \(p\) is the pressure, \(u\) is the velocity, and \(E\) is the total energy. To close the equation, a state equation describing the relationship between pressure and energy is introduced. In this paper, we consider the state equation for an ideal gas:  
\[p=(\gamma-1)\bigg{(}\rho E-\frac{1}{2}\rho||u||^{2}\bigg{)}, \tag{21}\]  
where \(\gamma=1.4\) is the adiabatic coefficient. In this class of PDE, even if the initial conditions are smooth, the solutions will be discontinuous at some finite time [33]. In such a case, it is difficult to obtain analytical solutions. Thus, numerical methods are widely used to find approximate solutions to such equations. The goal is to perform non-oscillatory reconstruction around the discontinuities. Therefore, this example aims to verify the ability of PIDeepONet and MulSPIDeepONet to deal with one-dimensional Euler equations with moving contact discontinuities [34]. In this example, the range of the computational domain \(x\) is \([0,1]\), the initial discontinuity is at \(x=0.5\), and the left and right sides of the discontinuity can be described by  
\[(\rho_{L},u_{L},p_{L})=(a,0.1,1.0),(\rho_{R},u_{R},p_{R})=(1.0,0.1,1.0), \tag{22}\]  
where \(a\) is a variable parameter. To reduce the computational overhead, the value range for the training process is \(a\in[1.2,5.4]\). We use Dirichlet boundary conditions and specify the exact solution as follows:  
\[\rho(x,t)=\left\{ \begin{array}{ll}a,\ x<0.5+0.1t,\\ 1.0,\ x>0.5+0.1t,\end{array} \right.u(x,t)=0.1,\ p(x,t)=1.0. \tag{23}\]  
Sampling different points on a hypercube within the range of \(a\), we obtain sufficient training data. The prediction accuracy is represented by the variable parameter, which is not used during the model training. The training goal is to learn the solution operator that maps the variable parameter \(a\) to the solution of Equation (19). Sixty randomly distributed boundary points are used for training. The number of initial training points is 60 and the number of function training points is 1000. The distribution of training points under each set of training parameters is shown in Figure 11.  
Figure 11: Training data points. The blue points are the initial condition training points, orange and green points are the boundary training points, and red points are the equation constraint training points. In this example, the learning rate is set to 0.005, and the Adam optimizer is used to train the model for 60,000 steps. The exponential decay learning strategy is adopted. The learning rate attenuates every 3000 steps, and the attenuation rate is 0.95. The training set includes 50 sets of training data randomly sampled in the range \(a\in[1.2,5.4]\). The specific network size settings are listed in Table 4, which displays the depth of each network and the number of neurons in each layer. The results for the Loss during the training process are shown in Figure 12.  
Figure 13 shows the density and velocity curves of Equation (19) solved by PIDEepONet and MulSPIDeepONet at \(t=2.0\) for \(a=4.4\) without training. The experimental results show that the velocity curves of PIDEepONet and MulSPIDeepONet are almost consistent with the analytical solution, but the density curve given by PIDEepONet at the discontinuity has obvious errors.